Spring + Pendulum with damping

[Disla]

Hi,

I’m not able to find what’s wrong in my model.
It’s a Spring + Pendulum with damping. The system is preloaded by means of initial displacement (Rotation). The spring is stretched into tension while the bar is rotated 0.175 Rad.

The whole system should oscillate no further than its initial -/+0.175 Rad but seems like the spring is not properly constrained or would have disappeared after the release.
Any help is appreciated.

*NODE
1,0,-0.25,0
2,0.25,-0.246182,0
3,-0.005,-0.5,-0.005
4,-0.005,0,-0.005
5,-0.005,-0.25,-0.005
6,-0.005,0,0
7,-0.005,0,0.005
8,-0.005,-0.5,0
9,-0.005,-0.5,0.005
10,-0.005,-0.25,0.005
11,0,-0.5,-0.005
12,0.005,-0.5,-0.005
13,0,0,-0.005
14,0.005,0,-0.005
15,0,0,0.005
16,0.005,0,0.005
17,0,-0.5,0.005
18,0.005,-0.5,0.005
19,0.005,-0.25,-0.005
20,0.005,0,0
21,0.005,-0.25,0.005
22,0.005,-0.5,0
23,-0.25,-0.246,0
24,0,-0.25,0
*ELEMENT,TYPE=C3D20
1,3,4,7,9,12,14,16,18,5,6,10,8,19,20,21,
22,11,13,15,17
*ELEMENT,TYPE=DASHPOTA
2,2,1
*ELEMENT,TYPE=MASS
3,8
*NSET,NSET=N_MASS
1
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
*NSET,NSET=REF
1
*NSET,NSET=RIGI
1
5
10
19
21
*NSET,NSET=N_MASS1_2
1
*NSET,NSET=7
4
6
7
*NSET,NSET=16
14
16
20
*NSET,NSET=10
5
10
*NSET,NSET=21
19
21
*NSET,NSET=2_REF
1
*NSET,NSET=9
3
8
9
*NSET,NSET=17
11
17
*NSET,NSET=18
12
18
22
*NSET,NSET=SPRING_NODE
1
*ELSET,ELSET=MASSLESS_ROD
1
*ELSET,ELSET=E_DAMPER
2
*ELSET,ELSET=MASS_1
3
*ELSET,ELSET=EALL
1
2
3
*MATERIAL,NAME=MATERIAL
*ELASTIC,TYPE=ISOTROPIC
210000000000,0.2
*DENSITY
1
*SOLID SECTION,ELSET=MASSLESS_ROD,MATERIAL=MATERIAL
*DASHPOT,ELSET=E_DAMPER

22
*BOUNDARY
2,1,,0
2,2,,0
2,3,,0
6,3,,0
13,1,,0
13,2,,0
15,1,,0
15,2,,0
20,3,,0
23,1,,0
23,2,,0
23,3,,0
*MASS,ELSET=MASS_1
2
*RIGID BODY,NSET=RIGI,REF NODE=1,ROT NODE=24
*INITIAL CONDITIONS,TYPE=DISPLACEMENT
1,1,0.0435270343984
1,2,0.003818365273767
3,1,0.08713043610227
3,2,0.006766189859565
4,1,7.636730547533E-05
4,2,-0.000870540687968
5,1,0.04360340170387
5,2,0.002947824585799
6,1,7.636730547533E-05
6,2,-0.000870540687968
7,1,7.636730547533E-05
7,2,-0.000870540687968
8,1,0.08713043610227
8,2,0.006766189859565
9,1,0.08713043610227
9,2,0.006766189859565
10,1,0.04360340170387
10,2,0.002947824585799
11,1,0.0870540687968
11,2,0.007636730547533
12,1,0.0869777
12,2,0.00850727
14,1,-7.636730547533E-05
14,2,0.000870540687968
16,1,-7.636730547533E-05
16,2,0.000870540687968
17,1,0.0870540687968
17,2,0.007636730547533
18,1,0.0869777
18,2,0.00850727
19,1,0.04345066709292
19,2,0.004688905961735
20,1,-7.636730547533E-05
20,2,0.000870540687968
21,1,0.04345066709292
21,2,0.004688905961735
22,1,0.0869777
22,2,0.00850727
*ELEMENT,TYPE=SPRING2,ELSET=ESpring
4,1,23
*SPRING,ELSET=Espring
1,1
20.0
*BOUNDARY
23,1,3,0
*INITIAL CONDITIONS,TYPE=VELOCITY
1,1,0
1,2,0
1,3,0
2,1,0
2,2,0
2,3,0
3,1,0
3,2,0
3,3,0
4,1,0
4,2,0
4,3,0
5,1,0
5,2,0
5,3,0
6,1,0
6,2,0
6,3,0
7,1,0
7,2,0
7,3,0
8,1,0
8,2,0
8,3,0
9,1,0
9,2,0
9,3,0
10,1,0
10,2,0
10,3,0
11,1,0
11,2,0
11,3,0
12,1,0
12,2,0
12,3,0
13,1,0
13,2,0
13,3,0
14,1,0
14,2,0
14,3,0
15,1,0
15,2,0
15,3,0
16,1,0
16,2,0
16,3,0
17,1,0
17,2,0
17,3,0
18,1,0
18,2,0
18,3,0
19,1,0
19,2,0
19,3,0
20,1,0
20,2,0
20,3,0
21,1,0
21,2,0
21,3,0
22,1,0
22,2,0
22,3,0
23,1,0
23,2,0
23,3,0
24,1,0
24,2,0
24,3,0
*STEP,NLGEOM=YES,INC=1000,AMPLITUDE=STEP
*DYNAMIC
0.01,4,0,0.01
*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,V
*EL FILE
ENER
*EL PRINT, ELSET=ESpring
ELSE
*END STEP

[vicmw]

Is the problem that the magnitude of the rotation angle is greater than the initial 0.175 rad (not conserving energy), or that it behaves the same as if there was no spring?

[felixR]

Just a guess, maybe I’m misunderstanding the initial displacement condition in Calculix.

As I understand it, it is not possible in theory to apply an instantaneous displacement to the dashpot (or in general to a model with a continuous path just through dashpots).
In the simulation, it could lead to the dashpot ‘snapping back’ the pendulum with a lot of stored energy.

For example the rheological model in the sketch (spring and dashpot in parallel) represents a Kelvin-Voigt material.
I take the following statement about instantaneous deformations for such a model from this paper:
https://www.sciencedirect.com/topics/engineering/kelvin-voigt-model

Initial deformation of the Kelvin-Voigt body cannot occur instantaneously, because the dashpot requires some time to respond to the applied stress.

Could you try to remove the dashpot and observe, if the pendulum behaves as it should? Or to keep the dashpot and apply an initial velocity.

-Felix

[fgr]

@Disla
I believe you have problem with your initial condition
When i look at the unload structure, it looks like this:

But with the initial condition loaded, it looks like this:

[Disla]

Hi, Thanks all to take the time to look at this.

@Victor . The rotation angle is ok. I have also constrained the z initial displacement according to our last GitHub comments and the bar ends up at where it is intended. I’m sure I have a placement of each node and not a displacement. I have focus on the Stored energy and it looks like there is a peak of energy released at the first-time step for this configuration. Velocities looks (from my point of view) wrong.

First Recorded Time	Bar Strain Energy		Initial Amplitude	Minimum Amplitude
0.0010	131.231	velocity X = -676.703 mm/s	43.53	-230.313	Wrong Result
0.0020	0.619	velocity Y = 92.0117 mm/s

@FelixR Before approaching this problem I had previously solved the single (double parallel too) model with spring, dashpot, and mass with excellent agreement. Just in case I have deactivated the dashpot as you suggested expecting the simple pendulum behavior, but I get the same response. The model behaves as if it stores more energy than I initially give it.
I have been able to solve the problem imposing an initial velocity and also with an initial incremental force up to the desired position + sudden release.

I solved properly these two systems before(Underdamped, Overdamped and Critically Dumped). One important difference in this case is that the mass is not directly applied to the spring/dashpot system and rotations are involved.

@fgr Strange. Mine looks ok. ¿Maybe cgx can’t read properly initial displacement?

[fgr]

@Disla , I had misunderstood the input in the initial condition card. The input value is the magnitude of the displacement and not the new position so my picture of the loaded structure is wrong.

[Disla]

I found a solution. It requires a tiny initial time step.
Just the initial. Later it can increase.
I don’t understand exactly why.
Seem like the first step might be crucial to properly initiate the energy balance/impulse/velocities. Maybe a bug, I don’t know.

*NODE
1,0,-0.25,0
2,0.25,-0.246182,0
3,-0.005,-0.5,-0.005
4,-0.005,-6.123031769112E-19,-0.005
5,-0.005,-0.25,-0.005
6,-0.005,-6.123031769112E-19,-8.673617379884E-19
7,-0.005,-6.123031769112E-19,0.005
8,-0.005,-0.5,-8.673617379884E-19
9,-0.005,-0.5,0.005
10,-0.005,-0.25,0.005
11,-8.040714361696E-17,-0.5,-0.005
12,0.005,-0.5,-0.005
13,0,0,-0.005
14,0.005,-1.387778780781E-19,-0.005
15,0,0,0.005
16,0.005,-1.249000902703E-18,0.005
17,-8.040714361696E-17,-0.5,0.005
18,0.005,-0.5,0.005
19,0.005,-0.25,-0.005
20,0.005,-1.387778780781E-19,-8.673617379884E-19
21,0.005,-0.25,0.005
22,0.005,-0.5,-8.673617379884E-19
23,0,0,0
24,0,0,0
25,-0.25,-0.25,0
26,0,-0.25,0
*ELEMENT,TYPE=C3D20
1,3,4,7,9,12,14,16,18,5,6,10,8,19,20,21,
22,11,13,15,17
*ELEMENT,TYPE=DASHPOTA
2,2,1
*ELEMENT,TYPE=MASS
3,8
4,3
5,9
6,11
7,12
8,17
9,18
10,22
*NSET,NSET=N_MASS
1
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
*NSET,NSET=REF
1
*NSET,NSET=RIGI
1
5
10
19
21
*NSET,NSET=N_MASS1_2
1
*NSET,NSET=7
4
6
7
*NSET,NSET=16
14
16
20
*NSET,NSET=10
5
10
*NSET,NSET=21
19
21
*NSET,NSET=2_REF
1
*NSET,NSET=9
3
8
9
*NSET,NSET=17
11
17
*NSET,NSET=18
12
18
22
*NSET,NSET=SPRING_NODE
1
*ELSET,ELSET=MASSLESS_ROD
1
*ELSET,ELSET=E_DAMPER
2
*ELSET,ELSET=MASS_1
3
4
5
6
7
8
9
10
*ELSET,ELSET=EALL
1
2
3
4
5
6
7
8
9
10
*MATERIAL,NAME=STEEL
*ELASTIC,TYPE=ISOTROPIC
210000000000,0.2
*DENSITY
1
*SOLID SECTION,ELSET=MASSLESS_ROD,MATERIAL=STEEL
*DASHPOT,ELSET=E_DAMPER

22.5
*BOUNDARY
2,1,,0
2,2,,0
2,3,,0
6,3,,0
13,1,,0
13,2,,0
15,1,,0
15,2,,0
20,3,,0
25,1,,0
25,2,,0
25,3,,0
*MASS,ELSET=MASS_1
0.25
*RIGID BODY,NSET=RIGI,REF NODE=1,ROT NODE=26
*INITIAL CONDITIONS,TYPE=DISPLACEMENT
1,1,0.0435270343984
1,2,0.003818365273767
3,1,0.08713043610227
3,2,0.006766189859565
4,1,7.636730547533E-05
4,2,-0.000870540687968
5,1,0.04360340170387
5,2,0.002947824585799
6,1,7.636730547533E-05
6,2,-0.000870540687968
7,1,7.636730547533E-05
7,2,-0.000870540687968
8,1,0.08713043610227
8,2,0.006766189859565
9,1,0.08713043610227
9,2,0.006766189859565
10,1,0.04360340170387
10,2,0.002947824585799
11,1,0.0870540687968
11,2,0.007636730547533
12,1,0.0869777
12,2,0.00850727
14,1,-7.636730547533E-05
14,2,0.000870540687968
16,1,-7.636730547533E-05
16,2,0.000870540687968
17,1,0.0870540687968
17,2,0.007636730547533
18,1,0.0869777
18,2,0.00850727
19,1,0.04345066709292
19,2,0.004688905961735
20,1,-7.636730547533E-05
20,2,0.000870540687968
21,1,0.04345066709292
21,2,0.004688905961735
22,1,0.0869777
22,2,0.00850727
*ELEMENT,TYPE=SPRING1,ELSET=ESpring
12,1,25
*SPRING,ELSET=Espring
1
20.0
*INITIAL CONDITIONS,TYPE=VELOCITY
1,1,0
1,2,0
1,3,0
2,1,0
2,2,0
2,3,0
3,1,0
3,2,0
3,3,0
4,1,0
4,2,0
4,3,0
5,1,0
5,2,0
5,3,0
6,1,0
6,2,0
6,3,0
7,1,0
7,2,0
7,3,0
8,1,0
8,2,0
8,3,0
9,1,0
9,2,0
9,3,0
10,1,0
10,2,0
10,3,0
11,1,0
11,2,0
11,3,0
12,1,0
12,2,0
12,3,0
13,1,0
13,2,0
13,3,0
14,1,0
14,2,0
14,3,0
15,1,0
15,2,0
15,3,0
16,1,0
16,2,0
16,3,0
17,1,0
17,2,0
17,3,0
18,1,0
18,2,0
18,3,0
19,1,0
19,2,0
19,3,0
20,1,0
20,2,0
20,3,0
21,1,0
21,2,0
21,3,0
22,1,0
22,2,0
22,3,0
23,1,0
23,2,0
23,3,0
24,1,0
24,2,0
24,3,0
25,1,0
25,2,0
25,3,0
26,1,0
26,2,0
26,3,0
*STEP,NLGEOM=YES,INC=1000,AMPLITUDE=STEP
*DYNAMIC
1E-05,2,0,0
*DLOAD
EALL,GRAV,0,1,0,0
*DLOAD
EALL,GRAV,-9.8,0,1,0
*DLOAD
EALL,GRAV,0,0,0,1
*NODE FILE,GLOBAL=YES
U,V
*EL FILE
ENER
*EL PRINT, ELSET=ESpring
ELSE
*END STEP

[fgr]

@Disla , I’m not member of Mecway’s forum but since I there have seen that you there in the thread Spring, mass, damper tutorial are looking for a formula for the x(t) at critical damping, then you can find it all as textbook example “Chapter 2, Free Vibration of 1-degree-of-freedom Systems” in Schaum’s Outlines, “Theory and problems of Mechanical vibrations” by S. Graham Kelly.
Unfortunately I only have this badly scan of the formula, but with the reference you should be able to find a readable formula.

[Disla]

Thanks, I could finally compare my results with the solution obtained by theory. Very nice agreement.

[JuanP74]

@Disla

I don’t understand exactly why.

Hi Disla, the method used for solving the ODEs that describe the time evolution in a dynamic problem is unconditionally stable, that means that the solution will not diverge whatever the time increment size, however accuracy is another issue. If you need an accurate solution you need to decrease the integration time increment, as you noticed.

[Disla]

It surprised me. All the other tests run fine with an initial time step much larger. Seems like IC=Displacement makes the system more sensitive to first time increment. There are some huge velocities involved also.

Thank you

[fgr]

@Disla , in another old textbook of mine I found this exercise of a simple pendulum which is/(can be) comparable to the first example in this thread

I expect your example has been taken from Example 1: Spring-Pendulum with damping where the natural frequence is calculated to 4.7 Hz
but when I use the data from Example 1: Spring-Pendulum with damping I get a natural frequence of 4.82657203424115 Hz
So far, I’m not sure which result should be correct, nor have I had time to try to trace the cause of the discrepancy.

[Disla]

EDITED. I’m not sure your example is comparable to mine. The Spring position is diferent with respect to the bar CM. I will try to solve it to see what I get. Regarding my example , below comment is still valid.

Damping affects the natural frequency of the system (pendulum).

What you are actually obtaining in the result is wd (d-Damped) not wn (n -Natural).

[Disla]

Yep , it’s a different problem. I solve it for a random case inside the range of validity.
Small displacements (spring elongation is considered to be only in horizontal direction) . Deviation with formula is 0.35%

** Generated by Mecway 26
*NODE
1,-1.788933584601E-17,-0.1666666666667,-3.071906155376E-19
2,0.25,-0.3333333333333,-5.421010862428E-19
3,-0.005,-0.5,-0.005
4,-0.005,-6.123031769112E-19,-0.005
5,0.005,-0.1666666666667,-7.228014483237E-19
6,-0.005,-0.1666666666667,0.005
7,-0.005,-6.123031769112E-19,0.005
8,-1.778091562876E-17,-0.1666666666667,0.005
9,-0.005,-0.5,0.005
10,0.005,-0.1666666666667,0.005
11,-0.005,-0.08333333333333,-0.005
12,0.005,-0.5,-0.005
13,0.005,-0.08333333333333,-0.005
14,0.005,-1.387778780781E-19,-0.005
15,-0.005,-0.08333333333333,0.005
16,0.005,-1.249000902703E-18,0.005
17,0.005,-0.08333333333333,0.005
18,0.005,-0.5,0.005
19,0,0,-0.005
20,-0.005,-6.123031769112E-19,-8.673617379884E-19
21,0.005,-1.387778780781E-19,-8.673617379884E-19
22,0,0,0.005
23,0,0,0
24,-4.445228907191E-17,-0.3333333333333,-5.421010862428E-19
25,-0.25,0,0
26,-8.040714361696E-17,-0.5,-0.005
27,-0.005,-0.5,-8.673617379884E-19
28,0.005,-0.5,-8.673617379884E-19
29,-8.040714361696E-17,-0.5,0.005
30,-0.005,-0.4166666666667,-0.005
31,0.005,-0.4166666666667,-0.005
32,-0.005,-0.4166666666667,0.005
33,0.005,-0.4166666666667,0.005
34,-0.005,-0.3333333333333,-0.005
35,-4.46691295064E-17,-0.3333333333333,-0.005
36,0.005,-0.3333333333333,-0.005
37,-0.005,-0.3333333333333,-1.011922027653E-18
38,0.005,-0.3333333333333,-1.011922027653E-18
39,-0.005,-0.3333333333333,0.005
40,-4.46691295064E-17,-0.3333333333333,0.005
41,0.005,-0.3333333333333,0.005
42,-0.005,-0.25,-0.005
43,0.005,-0.25,-0.005
44,-0.005,-0.25,0.005
45,0.005,-0.25,0.005
46,-0.005,-0.1666666666667,-0.005
47,-1.778091562876E-17,-0.1666666666667,-0.005
48,0.005,-0.1666666666667,-0.005
49,-0.005,-0.1666666666667,-7.228014483237E-19
50,0,0,0
51,-4.445228907191E-17,-0.3333333333333,-5.421010862428E-19
*ELEMENT,TYPE=C3D20
1,46,4,7,6,48,14,16,10,11,20,15,49,13,21,17,
5,47,19,22,8
3,3,34,39,9,12,36,41,18,30,37,32,27,31,38,33,
28,26,35,40,29
4,34,46,6,39,36,48,10,41,42,49,44,37,43,5,45,
38,35,47,8,40
*ELEMENT,TYPE=DASHPOTA
2,2,24
*ELEMENT,TYPE=MASS
5,29
6,28
7,27
8,26
9,18
10,12
11,9
12,3
*NSET,NSET=N_MASS
3
9
12
18
26
27
28
29
*NSET,NSET=REF1
23
*NSET,NSET=RIGI1
4
7
14
16
19
20
21
22
*NSET,NSET=REF2
24
*NSET,NSET=RIGI2
34
35
36
37
38
39
40
41
*ELSET,ELSET=MASSLESS_ROD_MIDDLE
4
*ELSET,ELSET=E_DAMPER
2
*ELSET,ELSET=MASSLESS_ROD_TOP
1
*ELSET,ELSET=MASSLESS_ROD_BOTTOM
3
*ELSET,ELSET=MASS_1
5
6
7
8
9
10
11
12
*ELSET,ELSET=EALL
1
2
3
4
5
6
7
8
9
10
11
12
*MATERIAL,NAME=MASSLESS_ROD
*ELASTIC,TYPE=ISOTROPIC
210000000000,0
*DENSITY
1
*SOLID SECTION,ELSET=MASSLESS_ROD_MIDDLE,MATERIAL=MASSLESS_ROD
*DASHPOT,ELSET=E_DAMPER

10
*SOLID SECTION,ELSET=MASSLESS_ROD_TOP,MATERIAL=MASSLESS_ROD
*SOLID SECTION,ELSET=MASSLESS_ROD_BOTTOM,MATERIAL=MASSLESS_ROD
*BOUNDARY
2,1,,0
2,2,,0
2,3,,0
5,3,,0
8,1,,0
8,2,,0
23,3,,0
24,3,,0
25,1,,0
25,2,,0
25,3,,0
47,1,,0
47,2,,0
49,3,,0
*MASS,ELSET=MASS_1
0.25
*RIGID BODY,NSET=RIGI1,REF NODE=23,ROT NODE=50
*RIGID BODY,NSET=RIGI2,REF NODE=24,ROT NODE=51
*TIME POINTS,NAME=timepointsname,GENERATE
0,2,0.005
*ELEMENT,TYPE=SPRING2,ELSET=ESpring
16,23,25
*SPRING,ELSET=Espring
1,1
20.0

*INITIAL CONDITIONS,TYPE=VELOCITY
N_Mass,1,0.1
N_Mass,2,0.0
N_Mass,3,0.0




*STEP,NLGEOM=YES,INC=1000,AMPLITUDE=STEP
*DYNAMIC
0.0001,2,0,0.005
*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,TIME POINTS=timepointsname
U,V
*EL FILE,TIME POINTS=timepointsname
ENER
*END STEP

[fgr]

@Disla , thanks for remark me about the frequency. I should have read the text instead of having focus on the index in the formula and in fact looking at the formula should have giving me the right answer.

Next I believe we misunderstand each other. We have these 2 cases where the spring for the rightmost have been moved the same distance below the origin, L1 has just changed sign. The spring energy around origin is equal except for the sign, but since the sign is neutralized by raising to the power of 2 the frequency will also be equal if the rod is truly massless, so something must be wrong.

[Disla]

Your wn book formula is wrong. It doesn’t pass the unit consistency check and Damping term suspiciously dissapear for L1=0.
That term is so small for my values that the error passed unnoticed. We should try to find another book with that formula.

Edited: There is no spring position dependency (L1 ) in the damping term. It has sense.

[Disla]

[fgr]

Nicely recognized, units don’t match, so the formula must be wrong

[JuanP74]

https://holooly.com/solutions-v2-1/a-simple-pendulum-is-pivoted-at-point-o-as-shown-in-figure-4e1a-assuming-the-mass-of-the-vertical-rod-is-negligible-and-the-oscillation-is-small-a-derive-the-equation-of-motion-for-this-simple-pend/

[fgr]

As the formula which already has been put in last reference I have also come to the same conclusion that this one should be the right formula for the damped omega. Hereby with correct index

For documentation, please look at the pdf and please note if you should have correction to the solved equations>
SimplePendulumDamped solved equations

For those not familiar with wxMaxima then %e = e, the base of the natural logarithm