ObjectANCFBeam

A 3D beam finite element based on the absolute nodal coordinate formulation, using two nodes. The localPosition \(x\) of the beam ranges from \(-L/2\) (at node 0) to \(L/2\) (at node 1). The axial coordinate is \(x\) (first coordinate) and the cross section is spanned by local \(y\)/\(z\) axes; assuming dimensions \(w_y\) and \(w_z\) in cross section, the local position range is \(\in [[-L/2,L/2],\, [-wy/2,wy/2],\, [-wz/2,wz/2] ]\). NOTE: Requires further development and tests!

Interface

Python names: ObjectANCFBeam or ANCFBeam, and VANCFBeam for its visualization
Nodes it takes: NodePointSlope23
Body markers that can be placed on it: MarkerBodyMass, MarkerBodyPosition, MarkerBodyRigid

Parameters

The parameters of the item; in a dictionary, its type is ‘ANCFBeam’:

name	type	size	default	description
name	String		‘’	objects’s unique name
nodeNumbers	NodeIndex2	2	[invalid (-1), invalid (-1)]	two node numbers for beam element
length	PReal		0.	(symbol: \(L\)) [SI:m] reference length of beam; such that the total volume (e.g. for volume load) gives \(\rho A L\); must be positive; must be given: the default is only a placeholder
sectionData	BeamSection		BeamSection()	data as given by exudyn.BeamSection(), defining inertial, stiffness and damping parameters of beam section.
crossSectionPenaltyFactor	Vector3D	3	[1.,1.,1.]	(symbol: \(k_{cs} = [f_{yy},\,f_{zz},\,f_{yz}]\tp\)) [SI:1] additional penalty factors for cross section deformation, which are in total \(k_{cs} = [f_{yy}\cdot EA,\, f_{zz}\cdot EA,\, f_{yz}\cdot (GA_y+GA_z)]\tp\)
crossSectionDamping	Vector3D	3	[0.,0.,0.]	(symbol: \(d_{cs} = [d_{fyy},\,d_{fzz},\,d_{fyz}]\tp\)) [SI:1] viscous damping according to penalty factors for cross section deformation; the damping is relative to the stiffness and should be thus usually much smaller than 1; the viscous damping factors read \(d_{cs} = [d_{fyy}\cdot EA,\, d_{fzz}\cdot EA,\, d_{fyz}\cdot (GA_y+GA_z)]\tp\)
visualization	VObjectANCFBeam			parameters for visualization of item

Renamed parameters, still taken with a DeprecationWarning: physicsLength (deprecated since 1.12.258, removed in 2031): use length.

Visualization parameters

The parameters of VObjectANCFBeam, given as visualization:

name	type	size	default	description
show	Bool		True	set true, if item is shown in visualization and false if it is not shown; geometry is defined by sectionGeometry
sectionGeometry	BeamSectionGeometry		BeamSectionGeometry()	defines cross section shape used for visualization and contact
color	Float4	4	[-1.,-1.,-1.,-1.]	RGBA color of the object; if R==-1, use default color

Drawing

Drawn as all finite elements; its sectionGeometry extruded along the axis, of bodies.beams.axialTiling segments; without section geometry a line, colored along it in a contour plot.

Settings beyond those of all finite elements: bodies.beams.axialTiling.

Output variables

Available as OutputVariableType in sensors, Get...Output() and other functions:

output variable	symbol	description
Position		global position vector of local position vector
Displacement		global displacement vector of local position vector
Velocity		global velocity vector of local position vector
VelocityLocal		local (cross section) velocity vector of local position vector
AngularVelocity		global angular velocity vector of local (axis) position vector
AngularVelocityLocal		local angular velocity vector of local (axis) position vector
Acceleration		global acceleration vector of local position vector
Rotation		3D Tait-Bryan rotation components of cross section rotation
RotationMatrix		rotation matrix of cross section rotation as 9D vector
KineticEnergy	\(T = \frac{1}{2} \dot\qv\tp \Mm\, \dot\qv\)	kinetic energy from the mass matrix of the current state and the velocities of the nodes; current configuration only; localPosition must be \([0,0,0]\)
PotentialEnergy	\(V = \frac{1}{2}\int_0^L \tkappa\tp \Km_\kappa \tkappa + \tgamma\tp \Km_\gamma \tgamma + \teps_\mathrm{cs}\tp \Km_\mathrm{cs} \teps_\mathrm{cs}\, dx\)	elastic energy of twist and curvature, axial and shear deformation and the penalized cross-section deformation, with the deformations and integration rules of the elastic forces; current configuration only; localPosition must be \([0,0,0]\)

Detailed description

This element is under development, as its class description says: what follows is what the implementation computes, not a finished formulation.

Nodes and coordinates

Two nodes of type NodePointSlope23, each with a position \(\rv_i\) and two cross section slopes \(\rv_{y,i}\), \(\rv_{z,i}\); 18 coordinates, taken as the sum of reference and current coordinates. The local axial coordinate is \(x \in [-L/2,\, L/2]\).

Kinematics and interpolation

Position and cross section vectors are interpolated linearly, and a point of the cross section at \((y,\,z)\) is at

\[ \rv(x,y,z) = S_0(x)\left(\rv_0 + y\,\rv_{y,0} + z\,\rv_{z,0}\right) + S_1(x)\left(\rv_1 + y\,\rv_{y,1} + z\,\rv_{z,1}\right) , \quad S_0 = \frac{1}{2} - \frac{x}{L} , \;\; S_1 = \frac{1}{2} + \frac{x}{L} . \]

Strains

From the cross section vectors \(\rv_y\), \(\rv_z\) the element builds an orthonormal basis \(\tv_1,\, \tv_2,\, \tv_3\) of the cross section (\(\tv_3 \parallel \rv_z\), \(\tv_1 \perp \rv_y,\,\rv_z\)), and measures

strain	definition	stiffness
axial and shear	\(\gamma_1 = \tv_1\tp\rv' - 1\), \(\gamma_2 = \tv_2\tp\rv'\), \(\gamma_3 = \tv_3\tp\rv'\)	\([EA,\; GA_y,\; GA_z]\) = `axialShearStiffness`
twist and curvature	\(\kv = \frac{1}{2}\sum_i \ev_i \times \ev_i'\) in the local basis	\([GJ_x,\; EI_y,\; EI_z]\) = `torsionalBendingStiffness`
cross section deformation	\(\frac{1}{2}(\rv_y\tp\rv_y - 1)\), \(\frac{1}{2}(\rv_z\tp\rv_z - 1)\), \(\frac{1}{2}\rv_y\tp\rv_z\)	penalty \([f_{yy} EA,\; f_{zz} EA,\; f_{yz}(GA_y+GA_z)]\) with `crossSectionPenaltyFactor`

each with a viscous damping of the same form. The strains are measured against a straight, undeformed reference: there are no reference strains or curvatures, so a curved reference configuration is not stress-free.

Mass matrix and integration

The mass matrix is constant and integrated with a Gauss rule of order 3. The elastic forces are integrated with one point - reduced integration against locking of the linear element - for the curvature and the axial and shear terms, and with a Lobatto rule for the cross section penalty.

Marker interfaces

The position Jacobian is the interpolation above at the local position, so a force may act at any point of the cross section. The rotation of the cross section is the orthonormal frame of the slopes \(\rv_{,y}\) and \(\rv_{,z}\) at \(x\): \(\rv_{,z}\) normalized is its \(z\)-axis, \(\rv_{,y}\) orthogonalized against it and normalized its \(y\)-axis (the output variable Rotation). The angular velocity \(\tomega\) and the rotation Jacobian \(\Jm_{rot} = \partial \tomega / \partial \dot \qv\) are the derivatives of this frame (#2763), so a torque or a rigid marker may act on the element. The position is linear in the coordinates; the derivative of the transposed rotation Jacobian times a torque is computed by automatic differentiation.

Limitations

Linear element, many elements needed for bending; no reference strains; under development and not fully tested.

Mini example

#a cantilever of four 3D ANCF beam elements with cross-section deformation, loaded at the tip
L = 1; nElements = 4; F = -0.1
section = exu.BeamSection()
section.stiffnessMatrix = np.diag([1e5, 1e4, 1e4, 100, 100, 100]) #EA, GA_y, GA_z, GJ, EI_y, EI_z
section.massPerLength = 1
section.inertia = 0.01*np.eye(3)
mGround = mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=nGround, coordinate=0))
n0 = mbs.AddNode(NodePointSlope23(referenceCoordinates=[0,0,0, 0,1,0, 0,0,1]))
for i in range(9): #clamped: position and both slopes
    mbs.AddObject(ObjectConnectorCoordinate(markerNumbers=[mGround,
                  mbs.AddMarker(MarkerNodeCoordinate(nodeNumber=n0, coordinate=i))]))
for k in range(nElements):
    n1 = mbs.AddNode(NodePointSlope23(referenceCoordinates=[L*(k+1)/nElements,0,0, 0,1,0, 0,0,1]))
    mbs.AddObject(ObjectANCFBeam(nodeNumbers=[n0,n1], length=L/nElements, sectionData=section))
    n0 = n1
mbs.AddLoad(LoadForceVector(markerNumber=mbs.AddMarker(MarkerNodePosition(nodeNumber=n1)), loadVector=[0,0,F]))

mbs.Assemble()
mbs.SolveStatic()

#bending about y; converges to F*L^3/(3*EI_y) + F*L/GA_z = -0.3433e-3 with more elements
exu.sys['testResult'] = mbs.GetNodeOutput(n1, exu.OutputVariableType.Displacement)[2]*1000

Relevant Examples (Ex) and TestModels (TM) with weblink to github: ANCFBeamEigTest.py (TM), ANCFBeamTest.py (TM), ANCFCableBeamDampingTest.py (TM), energiesFlexibleBodiesTest.py (TM), geometricallyExactBeamTest.py (TM), rightAngleFrame.py (TM)