Can you try setting the initial linear velocity instead?Do you notice any
difference?
-------- Messaggio originale --------Da: Arati Bhattu <aa...@rice.edu> Data:
05/06/25 01:36 (GMT+01:00) A: ProjectChrono <projectchrono@googlegroups.com>
Oggetto: [chrono] Re: Double pendulum problem 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 plotsMass
=0.79 kglength =1 mInitial velocity (dtheta0) =0.01 rad/secimport pychrono as
chronoimport pychrono.irrlicht as chronoirrimport numpy as npfrom 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 framesys.Add(ground)#
------------------------------------------------------------------------------#
3. Create the pendulum body (initial position, orientation, velocity all in
global)#
------------------------------------------------------------------------------pend1
= chrono.ChBody()# Set global position of pendulum
CGpend1.SetPos(chrono.ChVector3d(0, -0.5, 0)) # 0.5m below joint in Y
(gravity) direction# Set mass and
inertiapend1.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 Zsys.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 originjoint1 =
chrono.ChLinkRevolute()joint1.Initialize(ground, pend1,
joint_frame)sys.Add(joint1)# Simulation parametersend_time = 50.0 #
secondsstep_size = 0.01time = 0.0# Initialize storage liststime_data =
[]pend1_x, pend1_y, pend1_z = [], [], []pend2_x, pend2_y, pend2_z = [], [], []#
Run the simulationwhile 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 dictionarymat_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 filesavemat('trial', mat_data)Comparison with analytical solutionOn
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,AratiOn Wednesday, May 21, 2025 at 1:41:54 AM UTC-6
dariom...@gmail.com 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_motionA slight change in the initial
conditions, little approximations, etc can lead to totally different results
and this is somehow expected.DarioIl giorno mercoledì 21 maggio 2025 alle
04:37:22 UTC+2 aa...@rice.edu ha scritto:HelloI 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 chronoimport pychrono.irrlicht as chronoirr# Create systemsys =
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 bodypend1 =
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
barpend1.SetAngVelParent(chrono.ChVector3d(0, 0, 10))sys.Add(pend1)# Visual: 1m
tall, half-dim = 0.5, shift visual down by 0.5shape1 =
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 bodypend2 =
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 barpend2.SetAngVelParent(chrono.ChVector3d(0, 0, 5))sys.Add(pend2)#
Visual for pend2shape2 = 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 originjoint1 = 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
setupvis =
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 3Dvis.AddTypicalLights()# Simulation loopwhile
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 projectchrono+unsubscr...@googlegroups.com.
To view this discussion visit
https://groups.google.com/d/msgid/projectchrono/264876b6-a1c3-44fa-b0c7-f8ef33f18aaan%40googlegroups.com.
--
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 projectchrono+unsubscr...@googlegroups.com.
To view this discussion visit
https://groups.google.com/d/msgid/projectchrono/6841471f.170a0220.71281.d1c4%40mx.google.com.
