*Coupling, distributing with beam elements

[lucas_bueno]

Hello everyone, hope you’re all doing fine.

In recent months I’ve been using beam elements with relative success in CalculiX. Stress distribution is pretty accurate to hand-calc, and displacements as well. Connections of different beam sections across the model are also very good for pipe elements. However, I tried to attach a mass element using *COUPLING and *DISTRIBUTING keywords but failed to achieve the expected results. In my work, I replaced the concentrated mass with *CLOAD and everything went fine, however, I’m still intrigued by what I did wrong.

I wanted the distributing approach to avoid over-stiffening the model using something like *COUPLING with *KINEMATIC, *RIGID BODY, or even artificially rigid beam elements. This way I could use frequency analysis with a mass and stiffness matrix more or less representative of the real structure, something like RBE3 can do in Nastran.

I would like to ask for the expertise of more experienced users than me, and for this, I prepared a non-working input file . If I throw it in CCX, it proceeds fine but ignores the coupling. I’m expecting to see a stress of about 100 MPa at the junction, but I get less than 1MPa which is due to gravity acting in the pipes. Can you help me spot the error in the coupling definition?

Best regards,
Lucas.

[...]
*Surface, type=element, name=_PickedSet11_CNS_
61,S2
** possible surfaces for beam elements: S1 S2 S3 S5
*Coupling, constraint name=Constraint-1, ref node=9999, surface=_PickedSet11_CNS_
*Distributing
1,6
*Element, type=MASS, elset=_PickedSet10_Inertia-1_
62, 9999
*Mass, elset=_PickedSet10_Inertia-1_
0.1
[...]

ccx_beam_coupling.inp

[vicmw]

I don’t think you’re supposed to use *DISTRIBUTING with mass. From the manual: “*In the degrees of freedom in the reference node a force/moment can be applied by a CLOAD card. This load system is replaced by an equivalent force distribution in the nodes belonging to the coupling surface.”

So it’s really just a higher-level way of applying force.

It used to be more capable but was also buggy in hard-to-predict and hard-to-notice ways, so I’m happy that it’s more restrictive now if that means it’s more reliable.

[xyont]

probably general constraint using equation can be used in this case. It’s need to satisfy both of mass translation and rotational (inertia).

@lucas_bueno

I’m expecting to see a stress of about 100 MPa at the junction

did eccentricity and rigid arm being considered? when i’m only connecting three DOF’s translation and simple rerun it has shown results below

[Disla]

Is that right?.I mean that abrupt stress end.

I 'm obtaining an Sxx more distributed over the suport area.
Dcoup3D element seems able to distribute a mass.

[xyont]

@Disla

Is that right?.I mean that abrupt stress end.

as stated, it’s only simple rerun. i did not investigate further due to lack of problem sketch and result by hand calculation.

@Disla

Dcoup3D element seems able to distribute a mass.

can it be work also when the input changes to

*Dload
_PickedSet2, GRAV, 9810., -1., 0., 0.

[Disla]

Now I see. Your *equation is not acting. That’s the result for the rigid T

[xyont]

@Disla

Your *equation is not acting

of course, i’m not yet to define rotational equation relationship, but probably it is possible.

to my understanding, Dcop3D element is lack of rotational transmitted. Maybe i can be wrong in the latest versions.

According to sketch problem restrained support condition, left beam should have zero internal force and stress as shown below. Did not know why your results using Dcoup3D shown uneven.

[Disla]

I have used *Equation to decouple T Rotational constrain.

[xyont]

can it be more specific and show detailed a part of your input being modified?

it’s look strange for me, a column did not show any curvature occurs and seems it was violates physical behavior.

[Disla]

It is tipically used to avoid transfering moment to the columns so they (columns) work ideally just in compression. Like in the pictures.
The rotational degrees of freedom are released just by the fact that you can’t impose Rotational links with *equation on beams.

[xyont]

@Disla

It is tipically used to avoid transfering moment to the columns so they (columns) work ideally just in compression. Like in the pictures.

i can not be commenting since there’s much modification from original input files attached. The problem actually is quite simple in connecting node of mass and beam end nodes. In the beginning, i’m just questioning is the eccentricity and rigid arm being considered? This is related to mass moment of inertia in which is not available in CalculiX.

@Disla

The rotational degrees of freedom are released just by the fact that you can’t impose Rotational links with *equation on beams.

it seems not, all translation and rotational DOF’s can be set in *Equation to be transmitted between two separated nodes.

[lucas_bueno]

@vicmw

So it’s really just a higher-level way of applying force.

Yes, I am aware of this. However, I expect (and will test of course) if I can manage to apply the concentrated mass to the mass matrix without changing (too much) the stiffness matrix. The best case scenario would be to transfer the forces, change the mass matrix, and not change much the stiffness matrix. Unfortunately, I could not test this since the coupling did not work in my few tries.

@xyont

did eccentricity and rigid arm being considered?

Since the gravity is acting vertically and the mass is located just above the end of the bar, I expect the eccentricity to not be a problem in this test case. This could be the next part of my investigation, but rotary inertia is not relevant to the work I was doing using this connection. Your test is exactly what I tried to achieve, maybe some confusion is made because for the vertical beam, the bending stress is Syy and you plotted Sxx. In your next post with the Mises equivalent we can see that the horizontal bar with a free end and the vertical one have the same stress at the T junction, that’s the expected output. As expected, the fixed horizontal bar has zero moment and should not have any bending stresses. Maybe this asks @Disla 's question about the abrupt stress end.

@xyont

general constraint using equation can be used in this case

@Disla

Dcoup3D element seems able to distribute a mass.

Thank you guys for pointing out these solutions. I need to understand better how the *EQUATION keyword works and will use this @xyont example as a reference. Could you please share your input file where you define the *DCOUP3D and *DISTRIBUTING COUPLING? I want to study this implementation as well. While reading the documentation, my first thought was that this keyword was not meant for mass elements connection. I’m glad it is . I’ll run some frequency checks to see how both approaches perform against some analytical examples.

@Disla

I have used *Equation to decouple T Rotational constrain.

Very interesting, Disla. I did not cite this before, but these bars are theoretically welded and belong to a scaffold, so unfortunately the release approach would not work for this specific case. But it is very nice to know that we can do this at CalculiX. Some time ago I saw similar approaches in Nastran.

[lucas_bueno]

@lucas_bueno

While reading the documentation, my first thought was that this keyword was not meant for mass elements connection

I thought this due to the following phrase in the documentation:

How can we apply the mass to the DCOUP3D without using the same node to apply the mass element?

[Disla]

Hi Lucas,.

This is the set up.
Constrained nodes contains the nodes where you want to distribute the mass load.It can’t contain nodes with alrready other SPC’s.
Please report your conclusions.
Regarding the *EQUATION keyword it’s interesting by itself and I will comment on other post not to mix.

*NODE
1,0,0,0
2,0.5,0,0
3,0.5,0,0
4,1,0,0
5,0.5,-0.5,-6.123031769112E-17
6,0.25,0,0
7,0.75,0,0
8,0.5,-0.25,-3.061515884556E-17
9,0.125,0,0
10,0.375,0,0
11,0.625,0,0
12,0.875,0,0
13,0.5,-0.125,-1.530757942278E-17
14,0.5,-0.375,-4.592273826834E-17
15,0.8125,0,0
16,0.9375,0,0
17,0.5,-0.0625,-7.65378971139E-18
18,0.5,-0.1875,-2.296136913417E-17
19,0.5,-0.3125,-3.826894855695E-17
20,0.5,-0.4375,-5.357652797973E-17
21,0.0625,0,0
22,0.1875,0,0
23,0.3125,0,0
24,0.4375,0,0
25,0.5625,0,0
26,0.6875,0,0
27,0.5,-0.40625,-4.974963312403E-17
28,0.5,-0.46875,-5.740342283542E-17
29,0.03125,0,0
30,0.09375,0,0
31,0.15625,0,0
32,0.21875,0,0
33,0.28125,0,0
34,0.34375,0,0
35,0.40625,0,0
36,0.46875,0,0
37,0.53125,0,0
38,0.59375,0,0
39,0.65625,0,0
40,0.71875,0,0
41,0.78125,0,0
42,0.84375,0,0
43,0.90625,0,0
44,0.96875,0,0
45,0.5,-0.03125,-3.826894855695E-18
46,0.5,-0.09375,-1.148068456708E-17
47,0.5,-0.15625,-1.913447427847E-17
48,0.5,-0.21875,-2.678826398986E-17
49,0.5,-0.28125,-3.444205370125E-17
50,0.5,-0.34375,-4.209584341264E-17
51,1,0.05,0
52,0.015625,0,0
53,0.046875,0,0
54,0.078125,0,0
55,0.109375,0,0
56,0.140625,0,0
57,0.171875,0,0
58,0.203125,0,0
59,0.234375,0,0
60,0.265625,0,0
61,0.296875,0,0
62,0.328125,0,0
63,0.359375,0,0
64,0.390625,0,0
65,0.421875,0,0
66,0.453125,0,0
67,0.484375,0,0
68,0.515625,0,0
69,0.546875,0,0
70,0.578125,0,0
71,0.609375,0,0
72,0.640625,0,0
73,0.671875,0,0
74,0.703125,0,0
75,0.734375,0,0
76,0.765625,0,0
77,0.796875,0,0
78,0.828125,0,0
79,0.859375,0,0
80,0.890625,0,0
81,0.921875,0,0
82,0.953125,0,0
83,0.984375,0,0
84,0.5,-0.015625,-1.913447427847E-18
85,0.5,-0.046875,-5.740342283542E-18
86,0.5,-0.078125,-9.567237139237E-18
87,0.5,-0.109375,-1.339413199493E-17
88,0.5,-0.140625,-1.722102685063E-17
89,0.5,-0.171875,-2.104792170632E-17
90,0.5,-0.203125,-2.487481656202E-17
91,0.5,-0.234375,-2.870171141771E-17
92,0.5,-0.265625,-3.252860627341E-17
93,0.5,-0.296875,-3.63555011291E-17
94,0.5,-0.328125,-4.01823959848E-17
95,0.5,-0.359375,-4.400929084049E-17
96,0.5,-0.390625,-4.783618569619E-17
97,0.5,-0.421875,-5.166308055188E-17
98,0.5,-0.453125,-5.548997540758E-17
99,0.5,-0.484375,-5.931687026327E-17
*ELEMENT,TYPE=B32R
1,7,76,41
2,41,77,15
3,15,78,42
4,42,79,12
5,12,80,43
6,43,81,16
7,16,82,44
8,44,83,4
9,2,84,45
10,45,85,17
11,17,86,46
12,46,87,13
13,13,88,47
14,47,89,18
15,18,90,48
16,48,91,8
17,8,92,49
18,49,93,19
19,19,94,50
20,50,95,14
21,14,96,27
22,27,97,20
23,20,98,28
24,28,99,5
25,1,52,29
26,29,53,21
27,21,54,30
28,30,55,9
29,9,56,31
30,31,57,22
31,22,58,32
32,32,59,6
33,6,60,33
34,33,61,23
35,23,62,34
36,34,63,10
37,10,64,35
38,35,65,24
39,24,66,36
40,36,67,2
41,2,68,37
42,37,69,25
43,25,70,38
44,38,71,11
45,11,72,39
46,39,73,26
47,26,74,40
48,40,75,7
*ELEMENT,TYPE=MASS
49,51
*NSET,NSET=ConstrainedNodes
4
*ELSET,ELSET=Component
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
*ELSET,ELSET=Horizontal
1
2
3
4
5
6
7
8
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
*ELSET,ELSET=Mass_1
49
*ELSET,ELSET=EAll
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
*MATERIAL,NAME=Steel
*ELASTIC,TYPE=ISOTROPIC
200000000000,0
*DENSITY
7850
*BEAM SECTION,ELSET=Component,MATERIAL=Steel,SECTION=PIPE
0.025,0.003
0,0,1
*BEAM SECTION,ELSET=Horizontal,MATERIAL=Steel,SECTION=PIPE
0.025,0.003
0,0,1
*BOUNDARY
1,1,,0
1,2,,0
1,3,,0
1,4,,0
1,5,,0
5,1,,0
5,2,,0
5,3,,0
5,4,,0
5,5,,0
*MASS,ELSET=Mass_1
100
*ELEMENT,TYPE=DCOUP3D
50,51

*ELSET,ELSET=DCOUP
50

*DISTRIBUTING COUPLING,ELSET=DCOUP
ConstrainedNodes,2
*STEP
*STATIC,SOLVER=PARDISO
*DLOAD
EAll,GRAV,0,1,0,0
*DLOAD
EAll,GRAV,-9.81,0,1,0
*DLOAD
EAll,GRAV,0,0,0,1
*NODE FILE,GLOBAL=YES
U
*EL FILE
S,NOE
*END STEP

[lucas_bueno]

The more we dig, the more we find…

I simplified the model a little bit for the frequency study (the T junction was meant to verify stress only) and ended up with a cantilever beam with mass at its tip. The mass is very close to the beam end.

I tried to find a closed-form solution for this simple arrangement, but could not find it for all types of vibrations. For the axial mode, I used the Daniel J. Inman - Engineering Vibration 3rd Edition book to calculate the first mode as 148.93Hz from the following equation:

So I modeled it in Abaqus in some different ways: MPC Beam, Rigid body, Tie, mass directly attached to the end of the beam… The only mode that I had confidence in the result (the axial one) produced very good agreement: 148.8 Hz. And even though this is not the best way to go, I “trusted” the results for the other shapes.

After this, I modeled the same beam in CalculiX using 4 approaches:

• *EQUATION
• *DISTRIBUTED COUPLING
• Mass directly attached to the end node of the beam
• *RIGID BODY

And then I compared the results with the ones produced by Abaqus.

Here, for the first 4 transversal modes, the agreement is very good: up to a 3% difference in the higher-order mode. The torsional mode also produced very good agreement, with no difference in the first decimal places. However, two things caught my attention:

• The axial mode (the only one I have confidence in due to the analytical results) is very, very far off. When analyzing the deformed shape, I’m sure it is a spurious mode. It seems that the expanded nodes of the last element are vibrating on their own. I tried to use all degrees of freedom in the *EQUATION, but nothing changed…

• The *RIGID BODY produced two zero-frequency modes. I’m not sure why it happened, seems some spurious condition as well. Also, this option produced additional modes that the others did not… As I am not sure which are the correct results for the transversal modes, I can not guarantee that these modes are wrong, but since Abaqus and all other CalculiX approaches failed to find these modes, I tend to think that they are not correct.

One may ask: why is the reference node so close to the end of the structure? The answer is that I tried to model something similar to the closed-form solution that I had in hand. However, I tested what happened when I moved the reference node from 1 to 1000mm above the beam end. In CalculiX, the results were the same, even for the 6 DOF *EQUATION keyword. For Abaqus, the results (both frequency and eigenshape) changed a lot when I did this. I do not know if the results are correct, but I’m sure that CalculiX does not account for the rotary inertia component.

In the end, the proposal was to know how to attach a mass to a beam. To this, I can conclude that three approaches work:

• *EQUATION with DOF’s 1 to 3 connecting beam nodes to mass nodes
• DCOUP3D elements inside a *DISTRIBUTED COUPLING definition
• Mass directly attached to the beam nodes

If the user wants to consider the rotary inertia of the mass, more tests are needed to find something that works (maybe a very rigid beam connecting both nodes).

https://drive.google.com/drive/folders/1oHXTBtQKfM0p5iRmEgjqW-KYJByWWRA6?usp=sharing

[Disla]

This are the analytical solutions for a Lumped mass of 100 Kg at the tip.
Poisson effects are neglected.I only have the Transverse modes.
My axial is 143.5Hz but Torsion 914.3Hz ?¿?
Results for Mass directly attached to the end node of the beam.

Transverse modes	wr [rad/s]	f [Hz]	Result f [Hz]	Deviation
1	27.70	4.409	4.405	-0.09%
2	1331	211.78	209.70	-0.98%
3	4306	685.33	667.40	-2.62%
4	8981	1429.36	1358.00	-4.99%
5	15356	2443.95	2253.00	-7.81%

Ref solution: Fundamentals of vibrations Meirovitch. Example 8.7. To find the general case without the restriction M/mL=1 one need to solve the trascendental equation (e) for beta and later put it into equation (h).

m is the linear weight Kg/m of the bar.L your bar lenght.

[Disla]

I have also compare with the Mecway Internal solver to check Analitical solution. Perfect agreement.

Transverse modes	wr [rad/s]	f [Hz]	MECWAY	Deviation
1	27.70	4.409	4.41	0.00%
		148.8	152.60
2	1331	211.78	211.80	-0.01%
3	4306	685.33	685.30	0.00%
4	8981	1429.36	1,429.00	0.03%
5	15356	2443.95	2,444.00	0.00%

[lucas_bueno]

Thank you very much for the reference to the book. Super detailed explanation!

Amazing to see the Mecway solver making a perfect agreement with the theory. Such a powerful software as well.

@Disla

This are the analytical solutions for a Lumped mass of 100 Kg at the tip.

Not a perfect agreement in the higher modes for CalculiX, but a somewhat acceptable deviation in the first two or three shapes. Also, this closes the coffin lid for the *RIGID BODY approach. Nice thing to know.

Thanks for the help on this topic, @Disla @xyont and @vicmw.

[xyont]

@lucas_bueno

I tested what happened when I moved the reference node from 1 to 1000mm above the beam end. In CalculiX, the results were the same, even for the 6 DOF *EQUATION keyword. For Abaqus, the results (both frequency and eigenshape) changed a lot when I did this.

it seems Abaqus have rotational inertia, but CalculiX does not. Required, another approach, and probably it can be possible since the eccentricity is explicitly modeled in geometry input.

[Disla]

What about artificially increasing the density of the tip element?. Could there be some evident issue with that option.?