Hello,
I'm currently working on a PyChrono model using the SMC (Smooth Contact) method to demonstrate how rolling friction behaves in simulations. The setup I'm aiming to replicate is illustrated in the attached diagram, where the goal is to observe different states
of motion: stationary, pure rolling, pure sliding, and combined rolling and sliding (as shown in the second image).
However, in my implementation, I am only able to achieve a combination of rolling and sliding, when chaning the friction values. Additionally, even with the wall set at a 0° incline and a friction coefficient of 0.25, the sphere continues to roll indefinitely
instead of coming to rest.
Any guidance or insights into what might be wrong with my setup would be greatly appreciated.
PS.From my countless tests, it seems like the material data is only applied during the first contact, however, I'm not quite sure.
Best regards,
Mathias Dalby Larsen
Code used
import pychrono as chrono
import pychrono.irrlicht as chronoirr
import numpy as np
# Simulation Parameters
sphere_radius = 0.2 # meters
sphere_mass = 5.0 # kg
initial_velocity = 0.0 # m/s (parallel to incline, pointing up)
incline_angle = 5 # degrees (fix the floor at this angle)
static_friction = 0.25 # μs
kinetic_friction = 0.25 # μk
Rolling_friction = 0.50 # Rolling friction coefficient
measure_time = 4.5 # Time after which velocities will be measured
sim_time = 10
wall_depth = 15.0 # From wall dimensions
# Convert incline angle to radians
alpha_rad = np.radians(incline_angle)
# Create system
sys = chrono.ChSystemSMC()
# Check the current normal contact force model
print("Normal contact force before setting:", sys.GetContactForceModel()) # Should return default value
# Try setting it to Hooke
#sys.SetContactForceModel(chrono.ChSystemSMC.Hooke) # Set Hooke model explicitly
# Try setting it to Hertz
sys.SetContactForceModel(chrono.ChSystemSMC.Hertz) # Set Hertz model explicitly
# Try setting it to PlainCoulomb
#sys.SetContactForceModel(chrono.ChSystemSMC.PlainCoulomb) # Set PlainCoulomb model explicitly
# Try setting it to Flores
#sys.SetContactForceModel(chrono.ChSystemSMC.Flores) # Set Flores model explicitly
# Check if it changed
print("Normal contact force after setting:", sys.GetContactForceModel()) # Should print 1 (Hertz)
# Check the current tangential contact force model
print("Tangential contact force before setting:", sys.GetTangentialDisplacementModel()) # Should return default value
# Try setting it to None
#sys.SetTangentialDisplacementModel(chrono.ChSystemSMC._None) # Set None model explicitly
# Try setting it to OneStep
#sys.SetTangentialDisplacementModel(chrono.ChSystemSMC.OneStep) # Set OneStep model explicitly
# Try setting it to MultiStep
sys.SetTangentialDisplacementModel(chrono.ChSystemSMC.MultiStep) # Set MultiStep model explicitly
# Check if it changed
print("Tangential contact force after setting:", sys.GetTangentialDisplacementModel()) # Should print 1 (Hertz)
# Create Collision detection system and solver
sys.SetCollisionSystemType(chrono.ChCollisionSystem.Type_BULLET)
sys.SetSolverType(chrono.ChSolver.Type_BARZILAIBORWEIN)
sys.GetSolver().AsIterative().SetMaxIterations(200)
sys.UseMaterialProperties(True)
# Set gravitational acceleration
sys.SetGravitationalAcceleration(chrono.ChVector3d(0, -9.81, 0))
# Create a material for the inclined wall
wall_mat = chrono.ChContactMaterialSMC()
wall_mat.SetStaticFriction(static_friction)
wall_mat.SetSlidingFriction(kinetic_friction)
wall_mat.SetRollingFriction(Rolling_friction/2)
wall_mat.SetRestitution(0.50)
wall_mat.SetYoungModulus(1e9)
wall_mat.SetPoissonRatio(0.3)
wall_mat.SetSpinningFriction(Rolling_friction * 20.0)
wall_mat.SetKn(1e9)
wall_mat.SetKt(1e9)
# Create the inclined wall (static floor)
wall = chrono.ChBodyAuxRef()
wall.SetFixed(True)
wall.SetInertiaXX(chrono.ChVector3d(0.1, 0.1, 0.1))
wall.EnableCollision(True)
# Apply correct rotation using `QuatFromAngleAxis`
rotation_axis = chrono.ChVector3d(1, 0, 0) # Rotate around X-axis
wall_rotation = chrono.QuatFromAngleAxis(-alpha_rad, rotation_axis)
wall.SetRot(wall_rotation)
# Position the wall correctly
wall.SetPos(chrono.ChVector3d(0, 0, 0))
# Add collision shape
wall_shape = chrono.ChCollisionShapeBox(wall_mat, 2, 0.1, wall_depth)
wall.AddCollisionShape(wall_shape)
# Add visual representation
wall_vis = chrono.ChVisualShapeBox(2, 0.1, wall_depth)
wall_vis.SetColor(chrono.ChColor(1, 0, 0)) # Red for visibility
wall.AddVisualShape(wall_vis)
# Add to system
sys.Add(wall)
# Create the rolling sphere
sphere = chrono.ChBody()
sphere.SetMass(sphere_mass)
sphere.SetInertiaXX(chrono.ChVector3d(2/5 * sphere_mass * sphere_radius**2,
2/5 * sphere_mass * sphere_radius**2,
2/5 * sphere_mass * sphere_radius**2))
# Function to calculate the circle's center position
def get_circle_position(radius, theta):
y_c = -(radius+0.05) * np.sin(theta)
z_c = (radius+0.05) * np.cos(theta)
return y_c, z_c
# Get the position of the circle center
circle_pos = get_circle_position(sphere_radius, alpha_rad)
sphere.SetPos(chrono.ChVector3d(0, circle_pos[1]+0.01, circle_pos[0]+0.01))
# Set initial velocity along the incline
sphere.SetPosDt(chrono.ChVector3d(0,initial_velocity * np.sin(alpha_rad),
initial_velocity * np.cos(alpha_rad)))
# Add sphere collision and visualization
sphere_mat = chrono.ChContactMaterialSMC()
sphere_mat.SetStaticFriction(static_friction)
sphere_mat.SetSlidingFriction(kinetic_friction)
sphere_mat.SetRollingFriction(Rolling_friction/2)
sphere_mat.SetRestitution(0.50)
sphere_mat.SetYoungModulus(1e9)
sphere_mat.SetPoissonRatio(0.3)
sphere_mat.SetSpinningFriction(Rolling_friction * 20.0)
sphere_mat.SetKn(1e9)
sphere_mat.SetKt(1e9)
# Assign material directly to the body
sphere_shape = chrono.ChCollisionShapeSphere(sphere_mat, sphere_radius)
sphere.AddCollisionShape(sphere_shape)
sphere.EnableCollision(True)
sphere_vis = chrono.ChVisualShapeSphere(sphere_radius)
sphere_vis.SetColor(chrono.ChColor(1, 1, 0)) # Yellow sphere
sphere.AddVisualShape(sphere_vis)
# Add sphere to system
sys.Add(sphere)
# Set up the Irrlicht visualization
vis = chronoirr.ChVisualSystemIrrlicht()
vis.AttachSystem(sys)
vis.SetWindowSize(1024, 768)
vis.SetWindowTitle('Rolling Ball on Inclined Wall')
vis.Initialize()
vis.AddSkyBox()
vis.AddCamera(chrono.ChVector3d(-2, 0, -0.5), chrono.ChVector3d(0, 0, 0))
vis.AddTypicalLights()
# Simulation loop
time = 0.0
dt = 10e-6 # Keep accurate physics time step
render_interval = 1000 # Render every 10th step for 10x faster visual effect
step_count = 0
realtime_timer = chrono.ChRealtimeStepTimer()
lin_vel = None
while vis.Run():
if step_count % render_interval == 0: # Render only every 10th step
vis.BeginScene()
vis.Render()
vis.EndScene()
# Advance simulation
sys.DoStepDynamics(dt)
time += dt
step_count += 1
# Adjust real-time sync for 10x faster visualization
if step_count % render_interval == 0:
realtime_timer.Spin(dt * render_interval)
# Measure velocities at exactly measure_time
if time >= measure_time and lin_vel is None: # Only measure once
lin_vel = sphere.GetLinVel().Length()
ang_acc = sphere.GetAngVelLocal().Length() # Angular acceleration magnitude
Rollslide = abs((sphere_radius * ang_acc))
print(f"Measured at {measure_time}s -> Lin Vel: {lin_vel}, Rot Change: {ang_acc}, Rot rollslide: {Rollslide}")
force = sphere.GetContactForce()
torque = sphere.GetContactTorque()
print(f"Time: {time:.3f} - Contact Force: {force.Length()}, Contact Torque: {torque.Length()}")
