Hi Dario,
Could you kindly look into the code? I also have tried running the simple
simulation with single compound pendulum with initial velocity. I have
figured out the 'pend1.SetAngVelParent(chrono.ChVector3d(0, 0, wz))' is not
working the way I want to input it. For single pendulum problem with
initial theta, the answer matches exactly with analytical solution. However
with initial velocity (dtheta0) the frequency matches but amplitude. I am
plotting X displacement of CG of pendulum . Here is the code for your
reference and comparison plots
Mass =0.79 kg
length =1 m
Initial velocity (dtheta0) =0.01 rad/sec
import pychrono as chrono
import pychrono.irrlicht as chronoirr
import numpy as np
from scipy.io import savemat
#
------------------------------------------------------------------------------
# 1. Create the Chrono physical system
#
------------------------------------------------------------------------------
sys = chrono.ChSystemSMC()
sys.SetGravitationalAcceleration(chrono.ChVector3d(0, -9.81, 0)) # Global
gravity
#
------------------------------------------------------------------------------
# 2. Create the ground body (fixed, at global origin)
#
------------------------------------------------------------------------------
ground = chrono.ChBody()
ground.SetFixed(True)
ground.SetPos(chrono.ChVector3d(0, 0, 0)) # Explicit in global frame
sys.Add(ground)
#
------------------------------------------------------------------------------
# 3. Create the pendulum body (initial position, orientation, velocity all
in global)
#
------------------------------------------------------------------------------
pend1 = chrono.ChBody()
# Set global position of pendulum CG
pend1.SetPos(chrono.ChVector3d(0, -0.5, 0)) # 0.5m below joint in Y
(gravity) direction
# Set mass and inertia
pend1.SetMass(0.79)
pend1.SetInertiaXX(chrono.ChVector3d(0.01, 0.01, 0.0658)) # principal
inertias
# Put dummy inertia in X and Y direction. Otherwise it gives Inf solution
# Set angular velocity in global frame (around Z-axis)
pend1.SetAngVelParent(chrono.ChVector3d(0, 0, 0.01)) # radians/sec in
global Z
sys.Add(pend1)
#
------------------------------------------------------------------------------
# 4. Revolute joint between ground and pend1 (at global origin, around
global Z)
#
------------------------------------------------------------------------------
# Frame at joint location (global coordinates)
joint_frame = chrono.ChFramed(chrono.ChVector3d(0, 0, 0)) # Revolute joint
at origin
joint1 = chrono.ChLinkRevolute()
joint1.Initialize(ground, pend1, joint_frame)
sys.Add(joint1)
# Simulation parameters
end_time = 50.0 # seconds
step_size = 0.01
time = 0.0
# Initialize storage lists
time_data = []
pend1_x, pend1_y, pend1_z = [], [], []
pend2_x, pend2_y, pend2_z = [], [], []
# Run the simulation
while time < end_time:
sys.DoStepDynamics(step_size)
# Get positions
pos1 = pend1.GetPos()
# Append data to lists
time_data.append(time)
pend1_x.append(pos1.x)
pend1_y.append(pos1.y)
pend1_z.append(pos1.z)
time += step_size
# Prepare data dictionary
mat_data = {
'time': np.array(time_data, dtype=np.float64),
'pend1_x': np.array(pend1_x, dtype=np.float64),
'pend1_y': np.array(pend1_y, dtype=np.float64),
'pend1_z': np.array(pend1_z, dtype=np.float64),
'pend2_x': np.array(pend2_x, dtype=np.float64),
'pend2_y': np.array(pend2_y, dtype=np.float64),
'pend2_z': np.array(pend2_z, dtype=np.float64),
}
# Save to .mat file
savemat('trial', mat_data)
Comparison with analytical solution
On Wednesday, May 21, 2025 at 9:19:56 AM UTC-6 Arati Bhattu wrote:
> Hi Dario,
>
> Thank you for responding.
>
> I have solved the same problem numerically using the geometry and initial
> conditions provided earlier. The expected output shows that the first link
> (pivoted to the ground) completes approximately one full revolution in a
> 0.8-second interval.
>
> As I am still relatively new to Chrono, there is a possibility that I may
> not be providing the inputs as intended.
>
> Best,
> Arati
> On Wednesday, May 21, 2025 at 1:41:54 AM UTC-6 [email protected] wrote:
>
>> Can you show the discrepancies that you are observing in the simulation?
>> Are you aware about the chaotic behaviour of the double pendulum?
>> Double_pendulum#Chaotic_motion
>> <https://en.wikipedia.org/wiki/Double_pendulum#Chaotic_motion>
>> A slight change in the initial conditions, little approximations, etc can
>> lead to totally different results and this is somehow expected.
>>
>> Dario
>>
>> Il giorno mercoledì 21 maggio 2025 alle 04:37:22 UTC+2 [email protected] ha
>> scritto:
>>
>>> Hello
>>>
>>> I am trying to model a 2D double pendulum problem with the following
>>> geometry and initial conditions. However, the solution does not match the
>>> analytical results. I would like to verify whether my inputs are correct.
>>>
>>> Also, if I plan to make one of the links flexible, what would be an
>>> efficient way to implement it?
>>>
>>> The code is included below.
>>>
>>>
>>> ....................................................................................................................................................
>>>
>>> import pychrono as chrono
>>> import pychrono.irrlicht as chronoirr
>>>
>>> # Create system
>>>
>>> sys = chrono.ChSystemNSC()
>>> sys.SetGravitationalAcceleration(chrono.ChVector3d(0, -9.81, 0))
>>>
>>> # 1. Ground body (fixed)
>>> ground = chrono.ChBody()
>>> ground.SetFixed(True)
>>> ground.SetPos(chrono.ChVector3d(0, 0, 0))
>>> sys.Add(ground)
>>>
>>> # 2. First pendulum body
>>> pend1 = chrono.ChBody()
>>> pend1.SetMass(0.79)
>>> pend1.SetInertiaXX(chrono.ChVector3d(0.01, 0.01, 0.0658))
>>> pend1.SetPos(chrono.ChVector3d(0, -0.5, 0)) # CG in center of bar
>>> pend1.SetAngVelParent(chrono.ChVector3d(0, 0, 10))
>>> sys.Add(pend1)
>>>
>>> # Visual: 1m tall, half-dim = 0.5, shift visual down by 0.5
>>> shape1 = chrono.ChVisualShapeBox(chrono.ChVector3d(0.01,1, 0.01))
>>> shape1.SetColor(chrono.ChColor(0.5, 0.5, 0.8))
>>> pend1.AddVisualShape(shape1)
>>>
>>> # 3. Second pendulum body
>>> pend2 = chrono.ChBody()
>>> pend2.SetMass(0.3)
>>> pend2.SetInertiaXX(chrono.ChVector3d(0.01, 0.01, 0.0250))
>>> pend2.SetPos(chrono.ChVector3d(0, -1.5, 0)) # CG in center of second bar
>>> pend2.SetAngVelParent(chrono.ChVector3d(0, 0, 5))
>>> sys.Add(pend2)
>>>
>>> # Visual for pend2
>>> shape2 = chrono.ChVisualShapeBox(chrono.ChVector3d(0.01, 1, 0.01))
>>> shape2.SetColor(chrono.ChColor(0.8, 0.5, 0.5))
>>> pend2.AddVisualShape(shape2)
>>>
>>> # 4. Revolute joint between ground and pend1 at origin
>>> joint1 = chrono.ChLinkRevolute()
>>> frame1 = chrono.ChFramed(chrono.ChVector3d(0, 0, 0))
>>> joint1.Initialize(ground, pend1, frame1)
>>> sys.Add(joint1)
>>>
>>> # 5. Revolute joint between pend1 and pend2 at (0, -1, 0)
>>> joint2 = chrono.ChLinkRevolute()
>>> frame2 = chrono.ChFramed(chrono.ChVector3d(0, -1, 0))
>>> joint2.Initialize(pend1, pend2, frame2)
>>> sys.Add(joint2)
>>>
>>> # 6. Visualization setup
>>> vis = chronoirr.ChVisualSystemIrrlicht()
>>> vis.AttachSystem(sys)
>>> vis.SetWindowSize(1024, 768)
>>> vis.SetWindowTitle("Double Pendulum")
>>> vis.Initialize()
>>> vis.AddSkyBox()
>>> vis.AddCamera(chrono.ChVector3d(0.5, 1.5, 1.0)) # Camera looking in 3D
>>> vis.AddTypicalLights()
>>>
>>> # Simulation loop
>>> while vis.Run():
>>> vis.BeginScene()
>>> vis.Render()
>>> vis.EndScene()
>>> sys.DoStepDynamics(0.005)
>>>
>>
--
You received this message because you are subscribed to the Google Groups
"ProjectChrono" group.
To unsubscribe from this group and stop receiving emails from it, send an email
to [email protected].
To view this discussion visit
https://groups.google.com/d/msgid/projectchrono/264876b6-a1c3-44fa-b0c7-f8ef33f18aaan%40googlegroups.com.
