Simple multibody system does not follow theoretical tractrix curve
Fabio D'Angeli
Hello,
I built a simple system in PyChrono where a body point should follow a tractrix curve, but the point trajectory does not follow the tractrix analytic curve. The system is built as follows:
A body “slider” that can move only along x axis (ChLinkMatePrismatic).
A body “rod” that is connected:
- - at the root to the “slider” by a ChLinkMateRevolute that allows only z axis rotation about the slider.
- - at the tip to the “ground” by a ChLinkMateGeneric where only the local X translation is constrained. In this way, the trajectory of the tip point should be tangent to the rod direction
The “rod” is loaded with a torque Mz.
The system motion looks reasonable to what is expected but, if I plot the position of rod’s tip, its trajectory does not follow exactly a tractrix curve as expected.
Another weird thing is that the Mz on the rod tip is not zero as it should be.
Any idea on how to get the tip to exactly follow a tractrix curve?
Dan Negrut
Hi Fabio,
Thanks for pasting the script, I used chrono-rag (https://github.com/uwsbel/chrono-rag) and claude-code to get to the bottom of your problem.
1. Your model is correct and Chrono is solving it essentially exactly. It is just not the model you intended. The subtlety is in the tip constraint: for a ChLinkMateGeneric, the constrained coordinates live in the LINK FRAME THAT ROTATES WITH THE ROD, not in a ground-fixed frame. So constraining x does not mean the tip velocity stays tangent to the rod. It means the fixed ground point (0, 1, 0) must always lie on the rod's centerline. What you built is a rod sliding through a stationary eyelet at (0, 1, 0) while its root is dragged along the x axis.
Here is the proof, from your own numbers. The eyelet mechanism has the closed form: with the root at (a, 0), the tip sits at x = a(1 - 1/sqrt(a^2+1)), y = 1/sqrt(a^2+1). I ran your script headless for the same 2 seconds: it ends with the tip at (21.86, 0.0437), and the closed form at the corresponding root position a = 22.87 gives (21.87, 0.0437). Four-digit agreement. The eyelet curve resembles a tractrix at first, then parts company in the tail: a tractrix decays exponentially, y ~ 2 exp(-x/L), while the eyelet curve decays only algebraically, y ~ L^2/x roughly. That slow tail is exactly the gap your overlay plot shows.
2. The nonzero Mz that puzzled you is bookkeeping, not physics. The constraint force acts at the eyelet point, but GetReaction2 reports the wrench at the tip frame, and moving a force to a different reporting point generates a transport moment, force times arm. At the end of your run: 0.0437 N times the 21.88 m from tip to eyelet = 0.956, and the logged Mz is 0.9569. On the ground-side frame, GetReaction1, Mz is identically zero throughout, which your own log also shows.
3. Now the important part: no rearrangement of mate constraints will give you a tractrix, because the tractrix condition is not a position constraint at all. The defining property is that the tip's VELOCITY stays parallel to the rod, the no-side-slip condition of a knife edge, skate, or trailer wheel. That is a nonholonomic constraint, velocity-level and non-integrable; the tractrix is the textbook example of one. Every ChLinkMate variant is holonomic, and Chrono currently has no built-in generic nonholonomic joint, so you impose the condition yourself. The clean way: drag the root with a motor, drop the tip mate entirely, and each step apply a stiff damping force at the tip that kills the velocity component perpendicular to the rod.
I validated that approach against the analytic curve before writing this: with the penalty below, the simulated tip follows the tractrix with a maximum deviation of 0.063 percent of the rod length over a 4-meter drag (max error 0.6 mm on L = 1 m), and the residual lateral slip velocity is 5e-5 m/s. The working core, PyChrono 10:
import pychrono.core as chrono
L, V, C_LAT, DT = 1.0, 0.5, 2.0e3, 2.0e-4
sys = chrono.ChSystemNSC()
sys.SetGravitationalAcceleration(chrono.ChVector3d(0, 0, 0))
ground = chrono.ChBody(); ground.SetFixed(True); sys.Add(ground)
slider = chrono.ChBodyEasyBox(0.2, 0.1, 0.1, 1000); sys.AddBody(slider)
rod = chrono.ChBodyEasyBox(0.05, L, 0.05, 1000)
rod.Move(chrono.ChVector3d(0, L/2, 0)); sys.AddBody(rod)
# motor drags the root; a linear motor supplies the prismatic constraint itself
motor = chrono.ChLinkMotorLinearPosition()
motor.Initialize(slider, ground, chrono.ChFramed(chrono.ChVector3d(0, 0, 0),
chrono.QuatFromAngleY(chrono.CH_PI_2)))
motor.SetMotionFunction(chrono.ChFunctionRamp(0.0, V))
sys.Add(motor)
rev = chrono.ChLinkMateRevolute()
rev.Initialize(rod, slider, chrono.ChFramed(chrono.ChVector3d(0, 0, 0)))
sys.Add(rev)
tip_local = chrono.ChVector3d(0, L/2, 0) # tip in rod body coords
rod.AddAccumulator() # force accumulator, index 0
while sys.GetChTime() < 8.0:
u = rod.TransformDirectionLocalToParent(chrono.ChVector3d(0, 1, 0))
nvec = chrono.ChVector3d(-u.y, u.x, 0) # in-plane normal to the rod
v_tip = rod.PointSpeedLocalToParent(tip_local)
v_lat = v_tip.x*nvec.x + v_tip.y*nvec.y
F = chrono.ChVector3d(-C_LAT*v_lat*nvec.x, -C_LAT*v_lat*nvec.y, 0)
rod.EmptyAccumulator(0)
rod.AccumulateForce(0, F, rod.TransformPointLocalToParent(tip_local), False)
sys.DoStepDynamics(DT)
Two notes on the snippet. C_LAT is a penalty, so the constraint is satisfied approximately; 2e3 against a 10 kg rod at dt = 2e-4 is stable and tight, and if you change masses or the step, keep C_LAT*dt comfortably below the rod mass. And the path is pure kinematics, so it does not depend on the drag speed or the masses, which is a good self-check: change V and the trajectory should not move.
For what it is worth, everything else in your script was right: the prismatic and revolute frames do what you intended (the mate joints use z as their characteristic axis, so your QuatFromAngleY(CH_PI_2) is correct), and your instinct to compare against the closed-form curve is exactly how these things get caught.
Good luck, I hope this helps.
Dan
------------------------------------------------
Robert and Laura Hensel Professor
NVIDIA CUDA Fellow
Department of Mechanical & Aerospace Engineering
Department of Electrical & Computer Engineering
Department of Computer Science
University of Wisconsin - Madison
4150ME, 1513 University Avenue
Madison, WI 53706-1572
------------------------------------------------
From: Fabio D'Angeli <dangel...@gmail.com>
Sent: Thursday, 03 September 2026 19:34:18
To: Dan Negrut <neg...@wisc.edu>
Subject: Re: [Attachment Removed] [chrono] Simple multibody system does not follow theoretical tractrix curve
Hello Dan,
Hereafter the python script text:
# -*- coding: utf-8 -*-
"""
Created on Thu Sep 3 17:10:37 2026
@author: fabio
"""
import pychrono.core as chrono
import pychrono.irrlicht as chronoirr
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import time
# Initialize the Chrono simulation system
sys = chrono.ChSystemNSC()
sys.SetGravitationalAcceleration(chrono.ChVector3d(0, 0, 0)) # Neglect gravity
# Create the ground body and add it to the simulation
ground = chrono.ChBody()
ground.SetName('ground')
sys.Add(ground)
ground.SetFixed(True) # The ground is fixed and does not move
ground.EnableCollision(False) # Disable collision detection for the ground
# Create the slider body
slider = chrono.ChBodyEasyBox(1.0, 0.1, 0.1, 1)
sys.AddBody(slider)
# Create rod body
rod = chrono.ChBodyEasyBox(0.1, 1.0, 0.1, 1)
rod.Move(chrono.ChVector3d(0.0,0.5,0.0))
sys.AddBody(rod)
# Add slider constraint
constr_slider = chrono.ChLinkMatePrismatic()
constr_slider.Initialize(slider,
ground,
chrono.ChFramed(chrono.ChVector3d(0, 0, 0), chrono.QuatFromAngleY(chrono.CH_PI_2)))
constr_slider.SetName('sliderLink')
sys.Add(constr_slider)
# Add root motor
# rotmotor = chrono.ChLinkMotorRotationTorque()
# rotmotor.Initialize(rod, # body A (slave) (motor)
# slider, # body B (master) (stator)
# chrono.ChFramed(chrono.ChVector3d(0, 0, 0))) # motor frame, in abs. coords
# mtorquetime = chrono.ChFunctionConst(10)
# rotmotor.SetTorqueFunction(mtorquetime)
# sys.Add(rotmotor)
constr_root = chrono.ChLinkMateRevolute()
constr_root.Initialize(rod,
slider,
chrono.ChFramed(chrono.ChVector3d(0, 0, 0)))
constr_root.SetName('rootLink')
sys.Add(constr_root)
load = chrono.ChForce()
# rod torque load
rod.AddForce(load)
load.SetMode(chrono.ChForce.TORQUE)
load.SetDir(chrono.ChVector3d(0.0,0.0,1.0))
load.SetMforce(1)
# slider force load
# slider.AddForce(load)
# load.SetMode(chrono.ChForce.FORCE)
# load.SetDir(chrono.ChVector3d(1.0,0.0,0.0))
# load.SetMforce(1)
# Add tip constraint
constr_tip = chrono.ChLinkMateGeneric()
constr_tip.Initialize(ground,
rod,
chrono.ChFramed(chrono.ChVector3d(0, 1.0, 0)))
constr_tip.SetConstrainedCoords(True, False, False, # x, y, z
False, False, False) # Rx, Ry, Rz
constr_tip.SetName('linkTip')
sys.Add(constr_tip)
# Create the Irrlicht visualization system
vis = chronoirr.ChVisualSystemIrrlicht()
vis.AttachSystem(sys) # Attach the Chrono system to the visualization
vis.SetWindowSize(1024, 768) # Set the window size for the visualization
vis.SetWindowTitle('Tractrix') # Set the window title
vis.Initialize() # Initialize the visualization system
vis.AddLogo(chrono.GetChronoDataFile('logo_chrono_alpha.png')) # Add a logo to the visualization
vis.AddSkyBox() # Add a skybox for better visual appearance
vis.AddGrid(0.1, 0.1, 20, 20,
chrono.ChCoordsysd(),
chrono.ChColor(0,0,0))
vis.EnableBodyFrameDrawing(True)
vis.EnableLinkDrawing(True)
vis.EnableLinkFrameDrawing(True)
#vis.ShowExplorer(True)
vis.ShowInfoPanel(True)
camera_offset = chrono.ChVector3d(0, 0, 3)
vis.AddCamera(camera_offset) # Add a camera to the visualization
vis.AddTypicalLights() # Add typical lights for better visualization
# Simulation loop
log_info = True # Flag to control logging
t = 0.0
t_max =2.0
dt = 0.001
res = np.zeros((int(t_max/dt)+1,7))
i=0
while vis.Run(): # and t<t_max:
vis.BeginScene() # Begin the scene for rendering
vis.Render() # Render the scene
vis.EndScene() # End the scene
sys.DoStepDynamics(dt) # Advance the simulation by a time step of dt seconds
t = sys.GetChTime()
# Log information
res[i,0] = t
res[i,1] = constr_tip.GetFrame2Abs().GetPos().x
res[i,2] = constr_tip.GetFrame2Abs().GetPos().y
res[i,3] = constr_tip.GetReaction1().force.x
res[i,4] = constr_tip.GetReaction1().torque.z
res[i,5] = constr_tip.GetReaction2().force.x
res[i,6] = constr_tip.GetReaction2().torque.z
#print(constr_tip.GetReaction1().torque.to_numpy())
# Update camera position
vis.UpdateCamera(0, slider.GetPos() + camera_offset , slider.GetPos())
# stop simulation
if t > t_max:
vis.GetDevice().closeDevice()
i +=1
#time.sleep(0.1)
df=pd.DataFrame(res, columns=['t', 'x', 'y', 'Fx_ground', 'Mz_ground', 'Fx_tip', 'Mz_tip'])
# plot tip trajectory
plot = df.plot(x='x', y='y')
plt.grid(True)
# Plot analytic tractix curve
p=np.linspace(0,20,100)
L = 1.0
x=L*(p/L-np.tanh(p/L))
y=L/np.cosh(p/L)
plt.plot(x,y, label='tractrix')
plt.grid(True)
plt.legend()
df.plot(x='t', y=['Fx_ground', 'Fx_tip'])
plt.grid(True)
df.plot(x='t', y=['Mz_ground', 'Mz_tip'])
plt.grid(True)
On Thu, 3 Sept 2026, 20:03 Dan Negrut, <neg...@wisc.edu> wrote:
Fabio – can you include the attachment here, as a text file.
It got stripped by the email server for security reasons.
Thank you,
Dan
------------------------------------------------
Robert and Laura Hensel Professor
NVIDIA CUDA Fellow
Department of Mechanical & Aerospace Engineering
Department of Electrical & Computer Engineering
Department of Computer Science
University of Wisconsin - Madison
4150ME, 1513 University Avenue
------------------------------------------------
--
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/5b644b13-8433-403c-9e27-8fbfe1e619fdn%40googlegroups.com.