Attaching the pyChrono file in *.txt format.
On Wednesday, January 7, 2026 at 4:23:20 PM UTC-7 Roger Bergua wrote:
> Hi all,
>
> I'm modeling a very easy cantilever beam in the vertical direction. The
> system has just one cylindrical hollow beam (Euler-Bernoulli) and a clamp
> condition (ChBody + ChLinkMateGeneric).
>
> The beam has the next properties:
> Length: = 2 m
> Outer diameter: Do = 0.2 m
> Inne diameter: Di = 0.15 m
> Density: rho = 7,860 kg/m^3
>
> The mass of the beam is:
> m = pi*((Do/2)^2-(Di/2)^2)*L*rho = 216.06 kg
>
> By applying one acceleration of 10 m/s^2 along x in the radial direction
> of the beam, the expected applied loads at the clamp side would be:
> Fx = 216.06*10 = 2160.6 N
> Fy = 0 N
> Fz = 0 N
> Mx = 0 Nm
> My = 216.06*10*(L/2) = 2160.6 Nm
> Mz = 0 Nm
>
> I can reproduce this values perfectly in pyChrono.
>
> However, by applying the -10 m/s^2 acceleration along the z direction
> (i.e., along the beam longitudinal axis) I should get:
> Fx = 0 N
> Fy = 0 N
> Fz = 216.06*(-10) = -2160.6 N
> Mx = 0 Nm
> My = 0 Nm
> Mz = 0 Nm
>
> Instead of this, pyChrono returns Fz = -617.3 N
>
> I'm unable to understand how this value is computed. Maybe there is a
> better/recommended way to get the loads at the clamp side or the beam
> itself?
>
> Attached you can find the model in pyChrono to reproduce the above results.
>
> Thanks for the support!
>
>
>
>
>
--
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/4f0211fc-e642-4e2f-89d4-583c1e267f13n%40googlegroups.com.
import pychrono as chrono
import pychrono.fea as fea
import pychrono.pardisomkl as mkl
import math
#%% Create the Chrono system
system = chrono.ChSystemNSC()
#%% Material properties
E = 210E9 # Young Modulus [N/m^2]
rho = 7860 # Density [kg/m^3]
nu = 0.3 # Poisson ratio [-]
G = E/(2*(1+nu)) # Shear Modulus [N/m^2]
#%% Create mesh and builder
height = 2 # Total height [m] of the beam
D_o = 0.2 # Outer diameter [m]
D_i = 0.15 # Inner diameter [m]
# Euler-Bernoulli co-rotational beams (ChElementBeamEuler).
mesh = fea.ChMesh()
builder = fea.ChBuilderBeamEuler()
beam_nodes = []
beam_section = fea.ChBeamSectionEulerAdvanced()
beam_section.SetYoungModulus(E)
beam_section.SetShearModulus(G)
beam_section.SetDensity(rho)
beam_section.SetArea(math.pi*((D_o/2)**2-(D_i/2)**2))
beam_section.SetIyy((math.pi/4)*((D_o/2)**4-(D_i/2)**4))
beam_section.SetIzz((math.pi/4)*((D_o/2)**4-(D_i/2)**4))
beam_section.SetJ((math.pi/2)*((D_o/2)**4-(D_i/2)**4) )
builder.BuildBeam(mesh,
beam_section, # Beam section
properties
1, # Elements per
section
chrono.ChVector3d(0, 0, 0), # Start point
chrono.ChVector3d(0, 0, height), # End point
chrono.ChVector3d(0, 1, 0))
# Storing the beam nodes:
last_node = builder.GetLastBeamNodes().back()
first_node = builder.GetLastBeamNodes().front()
#%% Boundary conditions:
# Clamp condition
ground = chrono.ChBody()
ground.SetFixed(True)
system.Add(ground)
rigid_link = chrono.ChLinkMateGeneric(True, True, True, True, True, True)
rigid_link.Initialize(first_node, ground, chrono.ChFramed(first_node.GetPos()))
system.Add(rigid_link)
# Acceleration:
# x acceleration:
system.SetGravitationalAcceleration(chrono.ChVector3d(10, 0, 0))
# z acceleration:
#system.SetGravitationalAcceleration(chrono.ChVector3d(0, 0, -10))
#%% Static analysis
system.Add(mesh)
solver = mkl.ChSolverPardisoMKL()
system.SetSolver(solver)
system.DoStaticNonlinear()
#%% Check applied forces
# Applied load (in rigid_link frame)
Fx = rigid_link.GetReaction2().force.x
Fy = rigid_link.GetReaction2().force.y
Fz = rigid_link.GetReaction2().force.z
Mx = rigid_link.GetReaction2().torque.x
My = rigid_link.GetReaction2().torque.y
Mz = rigid_link.GetReaction2().torque.z
print("Applied loads at the clamped end:")
print("Fx =", round(Fx, 1), " N")
print("Fy =", round(Fy, 1), " N")
print("Fz =", round(Fz, 1), " N")
print("Mx =", round(Mx, 1), " Nm")
print("My =", round(My, 1), " Nm")
print("Mz =", round(Mz, 1), " Nm")
