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
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.
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:
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:
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):
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.
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.
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] */
(%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 */
/* period */
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 */
(%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"]);