Hi Arati,
two ways to achieve proper results:
- *setting initial velocities*: in this case you need to *properly *set
the initial velocities, that means *both *linear and angular; in fact
the COG is not only rotating but also translating
body1->SetPosDt({angvel_init * 0.5, 0, 0});
body1->SetAngVelLocal({0,0,angvel_init});
- *temporarily adding a speed-imposing motor and then removing it
immediately after a DoAssembly*
Here is my result and my CPP code (sorry I don't play much on PyChrono).
[image: 2025-06-06 16_31_34-Figure 1.png]
#include "chrono/physics/ChSystemNSC.h"
#include "chrono/physics/ChBodyEasy.h"
#include "chrono/physics/ChLinkMotorRotationSpeed.h"
#include "chrono/core/ChRealtimeStep.h"
#include "chrono_irrlicht/ChVisualSystemIrrlicht.h"
using namespace chrono;
using namespace chrono::irrlicht;
//
-----------------------------------------------------------------------------
int main(int argc, char* argv[]) {
std::cout << "Copyright (c) 2017 projectchrono.org\nChrono version: "
<< CHRONO_VERSION << std::endl;
// Create a Chrono physical system
ChSystemNSC sys;
double angvel_init = 0.01;
// Create pendulum bodies
auto body0 = chrono_types::make_shared<ChBody>();
body0->SetFixed(true);
sys.Add(body0);
auto body1 = chrono_types::make_shared<ChBodyEasyBox>(0.1, 1, 0.1, //
x,y,z size
1, //
density
true, //
visualization?
false); //
collision?
body1->SetPos(ChVector3d(0, -0.5, 0));
body1->SetInertiaXX({0.01, 0.01, 0.0658});
//body1->SetPosDt({angvel_init * 0.5, 0, 0});
//body1->SetAngVelLocal({0,0,angvel_init});
body1->SetMass(0.79);
sys.Add(body1);
// Create a spherical joint
auto rev_joint = chrono_types::make_shared<ChLinkRevolute>();
rev_joint->Initialize(body1, body0, ChFramed());
sys.AddLink(rev_joint);
auto mot = chrono_types::make_shared<ChLinkMotorRotationSpeed>();
mot->SetMotorFunction(chrono_types::make_shared<ChFunctionConst>(angvel_init));
mot->Initialize(body1, body0, ChFramed());
sys.AddLink(mot);
// Create a second spherical joint
// Create the Irrlicht visualization system
auto vis = chrono_types::make_shared<ChVisualSystemIrrlicht>();
vis->AttachSystem(&sys);
vis->SetWindowSize(800, 600);
vis->SetWindowTitle("A simple pendulum example");
vis->Initialize();
vis->AddLogo();
vis->AddSkyBox();
vis->AddCamera(ChVector3d(0, 1, -1));
vis->AddTypicalLights();
// Simulation loop
double timestep = 0.001;
ChRealtimeStepTimer realtime_timer;
sys.DoAssembly(AssemblyAnalysis::FULL);
mot->SetDisabled(true);
while (vis->Run()) {
vis->BeginScene();
vis->Render();
tools::drawGrid(vis.get(), 2, 2, 20, 20, ChCoordsys<>(ChVector3d(0,
-20, 0), QuatFromAngleX(CH_PI_2)),
ChColor(0.3f, 0.5f, 0.5f), true);
vis->EndScene();
std::cout << sys.GetChTime() << ", " << body1->GetPos().x() << "; ";
sys.DoStepDynamics(timestep);
realtime_timer.Spin(timestep);
}
return 0;
}
Il giorno giovedì 5 giugno 2025 alle 09:28:40 UTC+2 Dario Mangoni ha
scritto:
> 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 <projec...@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 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 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_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 aa...@rice.edu
>>> 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 projectchron...@googlegroups.com.
> To view this discussion visit
> https://groups.google.com/d/msgid/projectchrono/264876b6-a1c3-44fa-b0c7-f8ef33f18aaan%40googlegroups.com
>
> <https://groups.google.com/d/msgid/projectchrono/264876b6-a1c3-44fa-b0c7-f8ef33f18aaan%40googlegroups.com?utm_medium=email&utm_source=footer>
> .
>
--
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/ca2240a7-f99a-4063-84e4-df7fcced17adn%40googlegroups.com.
