Pendulum benchmark - implicit dynamics


I’m trying to solve a simple pendulum problem using implicit dynamics. The geometry is just a cuboid placed at an angle of 45 degrees:

Gravity load is applied and the pendulum is left to swing for 0.7 s. NLGEOM is on.

Here are the BCs on the bottom face of the pendulum:

  • nodes on the middle edge have all displacements fixed - red symbols
  • rigid body constraint is applied to all the other nodes on that face and its reference point has all DOFs fixed apart from the rotation in Y direction (rotation axis of the pendulum) - yellow and green symbols

The numbers are:

  • dimensions: 100x100x1000 mm
  • density: 7.87E-09 tonne/mm^3
  • Young’s modulus: 205 000 MPa, Poisson’s ratio: 0.3
  • gravity load: 9810 mm/s^2

According to simple hand calcs, the resulting velocity at 0.7 s should be 7088.03 mm/s.

Here’s my result:

I have a feeling that the agreement could be better.

When I remove the rigid body constraint (together with corresponding BC) and leave only the BC applied to the middle edge, I get slightly worse result (max 6770 mm/s). I wonder what causes the difference but it’s small anyway.

Do you know what else can be done to improve the accuracy ? In Abaqus I got only a bit better result (6814 mm/s) with the same input file so it’s not a fault of CalculiX.

Input file: Dropbox - Pendulum.inp - Simplify your life

I had remade the example, but I found that the maximun speed is 7088.45 mm/s at 0.626s. At time 0.7 seg velocity is 6815.06 mm/s. Could you share the equations for the handmade calculations?


Sure, here are the analytical calculations (they are done in terms of angle, not time):

hand calcs

The formulas are from “Schaum’s Outline of Engineering Mechanics Dynamics”.

Are you sure that the max velocity is at 0.7 seg?

I assumed that based on this website where the benchmark is described: Pinned Bar Under Gravitational Load | Validation Case | SimScale

But I would have to measure the angle.

For what I have understood, the benchmark calculations were at 180° (that doesn´t agree with my results). But is the same the magnitude velocity than the tangencial velocity?

HI Sergio. Your result is right. If energy is preserved , maximum velocity (Kinetic Energy) must be at the lowest point where the potential energy is minimum.

I get Vy=7086 mm/s too.

Look at the pendulum in vertical position and maybe try to reduce slightly your time step . I’m using…



EDIT: HINT to check the accuracy of the result. Compute one full cycle and…

-Final heigh must be the same.
-Maximum velocity must be exactly at 1/4 the time period.

  • Check energy balance

For calculation the exact time it will be necessary to involving the jacobian elliptic integral so the formula for the exact period from a high angle will be

ω0 = 4.0 * sqrt ( J / ( a * g * m ) ) * elliptic_f( pi / 2.0 , k * k )

Time for a full period for the example will be

ω0 = 2.502696535602945

A couple of months ago run some test with ccx on the real pendulum just trying to imagine the effect from the stepsize and looked over several periods a small drift will occur no matter the size of steps.
The attached picture isn’t related to the in this thread described data

Thank you for all the replies. It was rather silly mistake - for some reason, I was expecting the analytically calculated velocity to occur at 0.7 s while it occurs when the pendulum is in a vertical position.

With enough increments, I can get the desired result:

Should’ve figured it out myself :wink:

I took a step further and prepared a simulation with a more pendulum-like geometry:

There’s no rigid body constraint now, only the middle edge is fixed in all 3 translational DOFs. That should be sufficient based on the previous example.

Of course, the analytical solution is different to account for the fact that there are 2 components with different shapes and masses:

new calc

I calculated the total mass m, the moment of inertia about the pivot point I and the distance from the center of mass to the pivot point h analytically and confirmed that using CAD measurements so all those values are correct. I’m just not 100% sure that the formula is good and that there are no silly mistakes in the simulation. I should also add that d means the total length of the pendulum (distance from the pivot point to the bottom of the spherical bob) since I’m interested in the maximum velocity at the lowest point.

Here’s what I got from the simulation (not exactly max since it’s not vertical but close enough - it just takes quite a lot of time to solve with a proper number of increments):

Here’s the .inp file: WeTransfer - Send Large Files & Share Photos Online - Up to 2GB Free

Do you know what can be wrong here ? I think that I’m out of ideas.

Do you have the step file?

When the pendulum starts 45dgr from vertical position then h needs to bigger than the pendulum length I believe you have forgotten to add the pendulum length to h

I would also say you are using the mass and Inetia from the step file not the mesh. I get 29.41 Kg

Sure, here it is: WeTransfer - Send Large Files & Share Photos Online - Up to 2GB Free

But h is the distance from the center of mass of the pendulum to the pivot point in the initial position (at 45 degrees). It’s the vertical dimension here:

I was expecting the mass to be lower due to coarse discretization (to reduce the time required to solve this) but not to that extent. I have to check this. However, there must be also another reason for the discrepancy.

The main reason is pointed by fgr.

The minimum potential energy mgH is found with the pendulum in vertical position. The center of mass has decreased its height:


Ok, I think I get this now - h should be defined this way:

Right ? When I replace h with h+h*sqrt(2) and correct the mass to account for the inaccuracy in discretization, I get 6187.83 mm/s - much better. Maybe with the simulation including a frame with the pendulum in the vertical position, I would get even closer to the analytical solution.

Right. Result could be very accurate.

Edited. It took 12 min.

M 31.19 Kg
h 1.64 m
I 29.67 Kg/m2
g 9.81 m/s2
d 1.096 m
w 5.816 s-1
v 6373 mm/s
Result 6367 mm/s
Deviation 0.09%

You can always calculate exact values by use of this input to WxMaxima. You will only need to replace input with your own values for the desired pendulum and observed point of velocity. The result is valid for any angle between 0=< θ <180

(%i4)	/*  Real pendulum with rod */;
	d: 679.3*sqrt(2);       /*  distance from center of gravity to center of rotation  [mm] */
	m: 31197.9;         /*  mass in [gram]  */
	J: 29.6744e+9;        /*  Moment of Initia at center of rotation  [g*mm^2]   */
	g: 9.81e+3;         /*  gravity in [mm/sec^2]  */
(d)	679.3*sqrt(2)
(m)	31197.9

(J)	2.96744*10^10
(g)	9810.0
(%i13)	θ0: 135*%pi/180;    /*  starting angle  [dgr] */
	Lop: 1094.8;             /*   distance from center of rotation to point of observed velocity [mm ] */  
	/*  Constant for starting point */
	k: sin(θ0/2.0)$
	/*   period   */
	ω0: 4.0*sqrt(J/(d*g*m))*elliptic_f(%pi/2.0,k*k)$
	time: float(ω0/4);            /*  time to max speed  */
	/*   amplitude   */
	θt(t) := 2.0*asin(k*jacobi_cd(sqrt(d*g*m/J)*t,k*k))$  
	angle: float(θt(ω0/4));    /*  angular position at time = period/4  */
	/*   velocity   */
	dθt(t) := Lop*sqrt(2.*d*g*m/J*(cos(θt(t))-cos(θ0)))$
	velocity: float(dθt(ω0/4));  /*  velocity for observed point at time = period/4  */
(θ0)	(3*%pi)/4
(Lop)	1094.8

(time)	0.7624894110172591

(angle)	-5.261237772340481*10^-16
(velocity)	6367.583486762282

(%i14)	wxplot2d([dθt(t)],[t,0.001,%pi*2],[grid, 10, 1],
	    [title,"Pendulum with rod,starting angle 135 dgr."],
	    [ylabel,"Velocity in mm/sec"]);
(%t14)	 (Graphics)