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/93d12f26-d688-4362-b60a-51b006a4f054n%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
react = rigid_link.GetReaction2() # Applied load
# 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")
