ObjectANCFCable2D
A 2D cable finite element using 2 nodes of type NodePoint2DSlope1. The localPosition of the beam with length \(L\)=physicsLength and height \(h\) ranges in \(X\)-direction in range \([0, L]\) and in \(Y\)-direction in range \([-h/2,h/2]\) (which is in fact not needed in the EOM).
Additional information for ObjectANCFCable2D:
- Requested
Node
type =Position2D
+Orientation2D
+Point2DSlope1
+Position
+Orientation
- Short name for Python =
Cable2D
- Short name for Python visualization object =
VCable2D
The item ObjectANCFCable2D with type = ‘ANCFCable2D’ has the following parameters:
- name [type = String, default = ‘’]:objects’s unique name
- physicsLength [\(L\), type = UReal, default = 0.]:[SI:m] reference length of beam; such that the total volume (e.g. for volume load) gives \(\rho A L\); must be positive
- physicsMassPerLength [\(\rho A\), type = UReal, default = 0.]:[SI:kg/m] mass per length of beam
- physicsBendingStiffness [\(EI\), type = UReal, default = 0.]:[SI:Nm\(^2\)] bending stiffness of beam; the bending moment is \(m = EI (\kappa - \kappa_0)\), in which \(\kappa\) is the material measure of curvature
- physicsAxialStiffness [\(EA\), type = UReal, default = 0.]:[SI:N] axial stiffness of beam; the axial force is \(f_{ax} = EA (\varepsilon -\varepsilon_0)\), in which \(\varepsilon = |{\mathbf{r}}^\prime|-1\) is the axial strain
- physicsBendingDamping [\(d_{K}\), type = UReal, default = 0.]:[SI:Nm\(^2\)/s] bending damping of beam ; the additional virtual work due to damping is \(\delta W_{\dot \kappa} = \int_0^L \dot \kappa \delta \kappa dx\)
- physicsAxialDamping [\(d_{\varepsilon}\), type = UReal, default = 0.]:[SI:N/s] axial damping of beam; the additional virtual work due to damping is \(\delta W_{\dot\varepsilon} = \int_0^L \dot \varepsilon \delta \varepsilon dx\)
- physicsReferenceAxialStrain [\(\varepsilon_0\), type = Real, default = 0.]:[SI:1] reference axial strain of beam (pre-deformation) of beam; without external loading the beam will statically keep the reference axial strain value
- physicsReferenceCurvature [\(\kappa_0\), type = Real, default = 0.]:[SI:1/m] reference curvature of beam (pre-deformation) of beam; without external loading the beam will statically keep the reference curvature value
- strainIsRelativeToReference [\(f\cRef\), type = Real, default = 0.]:if set to 1., a pre-deformed reference configuration is considered as the stressless state; if set to 0., the straight configuration plus the values of \(\varepsilon_0\) and \(\kappa_0\) serve as a reference geometry; allows also values between 0. and 1.
- nodeNumbers [type = NodeIndex2, default = [invalid [-1], invalid [-1]]]:two node numbers ANCF cable element
- useReducedOrderIntegration [type = Index, default = 0]:0/false: use Gauss order 9 integration for virtual work of axial forces, order 5 for virtual work of bending moments; 1/true: use Gauss order 7 integration for virtual work of axial forces, order 3 for virtual work of bending moments
- axialForceUserFunction [\(\mathrm{UF} \in \Rcal\), type = PyFunctionMbsScalarIndexScalar9, default = 0]:A Python function which defines the (nonlinear relations) of local strains (including axial strain and bending strain) as well as time derivatives to the local axial force; see description below
- bendingMomentUserFunction [\(\mathrm{UF} \in \Rcal\), type = PyFunctionMbsScalarIndexScalar9, default = 0]:A Python function which defines the (nonlinear relations) of local strains (including axial strain and bending strain) as well as time derivatives to the local bending moment; see description below
- visualization [type = VObjectANCFCable2D]:parameters for visualization of item
The item VObjectANCFCable2D has the following parameters:
- show [type = Bool, default = True]:set true, if item is shown in visualization and false if it is not shown
- drawHeight [type = float, default = 0.]:if beam is drawn with rectangular shape, this is the drawing height
- color [type = Float4, default = [-1.,-1.,-1.,-1.]]:RGBA color of the object; if R==-1, use default color
DESCRIPTION of ObjectANCFCable2D
The following output variables are available as OutputVariableType in sensors, Get…Output() and other functions:
Position
: \(\LU{0}{{\mathbf{p}}\cConfig(x,y,0)} = {\mathbf{r}}\cConfig(x) + y\cdot {\mathbf{n}}\cConfig(x)\)global position vector of local position \([x,y,0]\)Displacement
: \(\LU{0}{{\mathbf{u}}\cConfig(x,y,0)} = \LU{0}{{\mathbf{p}}\cConfig(x,y,0)} - \LU{0}{{\mathbf{p}}\cRef(x,y,0)}\)global displacement vector of local positionVelocity
: \(\LU{0}{{\mathbf{v}}(x,y,0)} = \LU{0}{\dot {\mathbf{r}}(x)} - y \cdot \omega_2 \cdot\LU{0}{{\mathbf{t}}(x)}\)global velocity vector of local positionVelocityLocal
: \(\LU{b}{{\mathbf{v}}(x,y,0)} = \LU{b0}{\Rot}\LU{0}{{\mathbf{v}}(x,y,0)}\)local velocity vector of local positionRotation
: \(\varphi = \mathrm{atan2}(r'_y, r'_x)\)(scalar) rotation angle of axial slope vector (relative to global \(x\)-axis)Director1
: \({\mathbf{r}}'(x)\)(axial) slope vector of local axis position (at \(y\)=0)StrainLocal
: \(\varepsilon\)axial strain (scalar) of local axis position (at Y=0)CurvatureLocal
: \(K\)axial strain (scalar)ForceLocal
: \(N\)(local) section normal force (scalar, including reference strains) (at \(y\)=0); note that strains are highly inaccurate when coupled to bending, thus consider useReducedOrderIntegration=2 and evaluate axial strain at nodes or at midpointTorqueLocal
: \(M\)(local) bending moment (scalar) (at \(y\)=0)AngularVelocity
: \(\tomega = [0,\, ,0,\, \omega_2]\)angular velocity of local axis position (at \(y\)=0)Acceleration
: \(\LU{0}{{\mathbf{a}}(x,y,0)} = \LU{0}{\ddot {\mathbf{r}}(x)} - y \cdot \dot\omega_2 \cdot\LU{0}{{\mathbf{t}}(x)}- y \cdot \omega_2 \cdot\LU{0}{\dot{\mathbf{t}}(x)}\)global acceleration vector of local positionAngularAcceleration
: \(\talpha = [0,\, ,0,\, \dot\omega_2]\)angular acceleration of local axis position
Definition of quantities
intermediate variables
|
symbol
|
description
|
---|---|---|
beam height
|
\(h\)
|
beam height used in several definitions, but effectively undefined. The geometry of the cross section has no influence except for drawing or contact.
|
local beam position
|
\(\pLocB=[x,\, y,\, 0]\tp\)
|
local position at axial coordinate \(x \in [0,L]\) and cross section coordinate \(y \in [-h/2, h/2]\).
|
beam axis position
|
\(\LU{0}{{\mathbf{r}}(x)} = {\mathbf{r}}(x)\)
|
|
beam axis slope
|
\(\LU{0}{{\mathbf{r}}'(x)} = {\mathbf{r}}'(x)\)
|
|
beam axis tangent
|
\(\LU{0}{{\mathbf{t}}(x)} = \frac{{\mathbf{r}}'(x)}{\Vert {\mathbf{r}}(x)'\Vert}\)
|
this (normalized) vector is normal to cross section
|
beam axis normal
|
\(\LU{0}{{\mathbf{n}}(x)} = [n_x,\, n_y]\tp = [-t_y,\, t_x]\tp\)
|
this (normalized) vector lies within the cross section and defines positive \(y\)-direction.
|
angular velocity
|
\(\omega_2 = (-r'_y \cdot \dot r'_x + r'_x \cdot \dot r'_y) / \Vert {\mathbf{r}}(x)'\Vert^2\)
|
|
rotation matrix
|
\(\LU{0b}{\Rot}\)
|
The Bernoulli-Euler beam is capable of large axial and bendig deformation as it employs the material measure of curvature for the bending.
Kinematics and interpolation
Note that in this section, expressions are written in 2D, while output variables are in general 3D quantities, adding a zero for the \(z\)-coordinate. ANCF elements follow the original concept proposed by Shabana . The present 2D element is based on the interpolation used by Berzeri and Shabana , but the formulation (especially of the elastic forces) is according to Gerstmayr and Irschik . Slight improvements for the integration of elastic forces and additional terms for off-axis forces and constraints are mentioned here.
The current position of an arbitrary element at local axial position \(x \in [0,L]\), where \(L\) is the beam length, reads
The derivative of the position w.r.t.the axial reference coordinate is denoted as slope vector,
The interpolation is based on cubic (spline) interpolation of position, displacements and velocities. The generalized coordinates \({\mathbf{q}} \in \Rcal^8\) of the beam element is defined by
in which \({\mathbf{r}}_0\) is the position of node 0 and \({\mathbf{r}}_1\) is the position of node 1, \({\mathbf{r}}'_0\) the slope at node 0 and \({\mathbf{r}}'_1\) the slope at node 1. Note that ANCF coordinates in the present notation are computed as sum of reference and current coordinates
which is used throughout here. For time derivatives, it follows that \(\dot {\mathbf{q}} = \dot {\mathbf{q}}\cCur\).
Position and slope are interpolated with shape functions. The position and slope along the beam are interpolated by means of
in which \({\mathbf{S}}\) is the shape function matrix,
with identity matrix \(\mathbf{I}_{2 \times 2} \in \Rcal^{2 \times 2}\) and the shape functions
Velocity simply follows as
Mass matrix
The mass matrix is constant and therefore precomputed at the first time it is needed (e.g., during computation of initial accelerations). The analytical form of the mass matrix reads
which is approximated using
with integration weights \(w(x_{ip})\), \(\sum w(x_{ip})=2\), and integration points \(x_{ip}\), given as,
Here, we use the Gauss integration rule with order 7, having \(n_{ip}=4\) Gauss points, see Section Integration Points. Due to the third order polynomials, the integration is exact up to round-off errors.
Elastic forces
The elastic forces \({\mathbf{Q}}_e\) are implicitly defined by the relation to the virtual work of elastic forces, \(\delta W_e\), of applied forces, \(\delta W_a\) and of viscous forces, \(\delta W_v\),
The virtual work of elastic forces reads ,
in which the axial strain is defined as
and the material measure of curvature (bending strain) is given as
in which \({\mathbf{e}}_3\) is the unit vector which is perpendicular to the plane of the planar beam element.
By derivation, we obtain the variation of axial strain
and the variation of \(K\)
The normal force (axial force) \(N\) in the beam is defined as function of the current strain \(\varepsilon\),
in which \(\varepsilon_0\) includes the (pre-)stretch of the beam, e.g., due to temperature or plastic deformation and \(\varepsilon\cRef\) includes the strain of the reference configuration. As can be seen, the reference strain is only considered, if \(f\cRef=1\), which allows to consider the reference configuration to be completely stress-free (but the default value is \(f\cRef=0\) !). Note that – due to the inherent nonlinearity of \(\varepsilon\) – a combination of \(\varepsilon_0\) and \(f\cRef=1\) is physically only meaningful for small strains. A factor \(f\cRef<1\) allows to realize a smooth transition between deformed and straight reference configuration, e.g. for initial configurations.
The bending moment \(M\) in the beam is defined as function of the current material measure of curvature \(K\),
in which \(K_0\) includes the (pre-)curvature of the undeformed beam and \(K\cRef\) includes the curvature of the reference configuration, multiplied with the factor \(f\cRef=1\), see the axial strain above.
Using the latter definitions, the elastic forces follow from Eq. (66).
The virtual work of viscous damping forces, assuming viscous effects proportial to axial streching and bending, is defined as
with material coefficients \(d_\varepsilon\) and \(d_K\). The time derivatives of axial strain \(\dot \varepsilon_p\) follows by elementary differentiation
as well as the derivative of the curvature,
The virtual work of applied forces reads
in which \({\mathbf{f}}_i\) are forces applied to a certain position \(x_f\) at the beam centerline. The second term contains a load per length \({\mathbf{b}}\), which is case of gravity vector \({\mathbf{g}}\) reads
Note that the variation of \({\mathbf{r}}\) simply follows as
Numerical integration of Elastic Forces
The numerical integration of elastic forces \({\mathbf{Q}}_e\) is split into terms due to \(\delta \varepsilon\) and \(\delta K\),
using different integration rules
with the integration points \(x_{ip}\) as defined in Eq. (65) and integration rules from Section Integration Points.
There are 3 different options for integration rules depending on the flag useReducedOrderIntegration
:
useReducedOrderIntegration
= 0: \(n_{ip}^\varepsilon = 5\) (Gauss order 9), \(n_{ip}^K = 3\) (Gauss order 5) – this is considered as full integration, leading to very small approximations; certainly, due to the high nonlinearity of expressions, this is only an approximation.useReducedOrderIntegration
= 1: \(n_{ip}^\varepsilon = 4\) (Gauss order 7), \(n_{ip}^K = 2\) (Gauss order 3) – this is considered as reduced integration, which is usually sufficiently accurate but leads to slightly less computational efforts, especially for bending terms.useReducedOrderIntegration
= 2: \(n_{ip}^\varepsilon = 3\) (Lobatto order 3), \(n_{ip}^K = 2\) (Gauss order 3) – this is a further reduced integration, with the exceptional property that axial strain and bending strain terms are computed at completely disjointed locations: axial strain terms are evaluated at \(0\), \(L/2\) and \(L\), while bending terms are evaluated at \(\frac{L}{2} \pm \frac{L}{2}\sqrt{1/3}\). This allows axial strains to freely follow the bending terms at \(\frac{L}{2} \pm \frac{L}{2}\sqrt{1/3}\), while axial strains are almost independent from bending terms at \(0\), \(L/2\) and \(L\). However, due to the highly reduced integration, spurious (hourglass) modes may occur in certain applications!
Note that the Jacobian of elastic forces is computed using automatic differentiation.
Access functions
For application of forces and constraints at any local beam position \(\pLocB=[x,\, y,\, 0]\tp\), the position / velocity Jacobian reads
with the normalized beam axis normal \(\LU{0}{{\mathbf{n}}} = [n_x,\, n_y]\tp\), see table above.
For application of torques at any axis point \(x\), the rotation / angular velocity Jacobian \(\frac{\partial \LU{0}{\omega(x)}}{\dot {\mathbf{q}}} \in \Rcal^{3 \times 8}\) reads
Userfunction: axialForceUserFunction(mbs, t, itemNumber, axialPositionNormalized, axialStrain, axialStrain_t, axialStrainRef, physicsAxialStiffness, physicsAxialDamping, curvature, curvature_t, curvatureRef)
A user function, which computes the axial force depending on time, strains and curvatures and
object parameters (stiffness, damping).
The object variables are provided to the function using the current values of the ANCFCable2D object.
Note that itemNumber represents the index of the object in mbs, which can be used to retrieve additional data from the object through
mbs.GetObjectParameter(itemNumber, ...)
, see the according description of GetObjectParameter
.
NOTE: this function has a different interface as compared to the bending moment function.
arguments / return
|
type or size
|
description
|
---|---|---|
mbs |
MainSystem
|
provides MainSystem mbs to which object belongs
|
t |
Real
|
current time in mbs
|
itemNumber |
Index
|
integer number \(i_N\) of the object in mbs, allowing easy access to all object data via mbs.GetObjectParameter(itemNumber, …)
|
axialPositionNormalized |
Real
|
axial position at the cable where the user function is evaluated; range is [0,1]
|
axialStrain |
Real
|
\(\varepsilon\)
|
axialStrain_t |
Real
|
\(\varepsilon_t\)
|
axialStrainRef |
Real
|
\(\varepsilon_0 + f\cRef \cdot \varepsilon\cRef\)
|
physicsAxialStiffness |
Real
|
as given in object parameters
|
physicsAxialDamping |
Real
|
as given in object parameters
|
curvature |
Real
|
\(K\)
|
curvature_t |
Real
|
\(\dot K\)
|
curvatureRef |
Real
|
\(K_0 + f\cRef \cdot K\cRef\)
|
returnValue
|
Real
|
scalar value of computed axial force
|
Userfunction: bendingMomentUserFunction(mbs, t, itemNumber, axialPositionNormalized, curvature, curvature_t, curvatureRef, physicsBendingStiffness, physicsBendingDamping, axialStrain, axialStrain_t, axialStrainRef)
A user function, which computes the bending moment depending on time, strains and curvatures and
object parameters (stiffness, damping).
The object variables are provided to the function using the current values of the ANCFCable2D object.
Note that itemNumber represents the index of the object in mbs, which can be used to retrieve additional data from the object through
mbs.GetObjectParameter(itemNumber, ...)
, see the according description of GetObjectParameter
.
NOTE: this function has a different interface as compared to the axial force function.
arguments / return
|
type or size
|
description
|
---|---|---|
mbs |
MainSystem
|
provides MainSystem mbs to which object belongs
|
t |
Real
|
current time in mbs
|
itemNumber |
Index
|
integer number \(i_N\) of the object in mbs, allowing easy access to all object data via mbs.GetObjectParameter(itemNumber, …)
|
axialPositionNormalized |
Real
|
axial position at the cable where the user function is evaluated; range is [0,1]
|
curvature |
Real
|
\(K\)
|
curvature_t |
Real
|
\(\dot K\)
|
curvatureRef |
Real
|
\(K_0 + f\cRef \cdot K\cRef\)
|
physicsBendingStiffness |
Real
|
as given in object parameters
|
physicsBendingDamping |
Real
|
as given in object parameters
|
axialStrain |
Real
|
\(\varepsilon\)
|
axialStrain_t |
Real
|
\(\varepsilon_t\)
|
axialStrainRef |
Real
|
\(\varepsilon_0 + f\cRef \cdot \varepsilon\cRef\)
|
returnValue
|
Real
|
scalar value of computed bending moment
|
User function example:
#define some material parameters
rhoA = 100.
EA = 1e7.
EI = 1e5
#example of bending moment user function
def bendingMomentUserFunction(mbs, t, itemNumber, axialPositionNormalized,
curvature, curvature_t, curvatureRef, physicsBendingStiffness,
physicsBendingDamping, axialStrain, axialStrain_t, axialStrainRef):
fact = min(1,t) #runs from 0 to 1
#change reference curvature of beam over time:
kappa=(curvature-curvatureRef*fact)
return physicsBendingStiffness*(kappa) + physicsBendingDamping*curvature_t
def axialForceUserFunction(mbs, t, itemNumber, axialPositionNormalized,
axialStrain, axialStrain_t, axialStrainRef, physicsAxialStiffness,
physicsAxialDamping, curvature, curvature_t, curvatureRef):
fact = min(1,t) #runs from 0 to 1
return (physicsAxialStiffness*(axialStrain-fact*axialStrainRef) +
physicsAxialDamping*axialStrain_t)
cable = ObjectANCFCable2D(physicsMassPerLength=rhoA,
physicsBendingStiffness=EI,
physicsBendingDamping = EI*0.1,
physicsAxialStiffness=EA,
physicsAxialDamping=EA*0.05,
physicsReferenceAxialStrain=0.1, #10
physicsReferenceCurvature=1, #radius=1
bendingMomentUserFunction=bendingMomentUserFunction,
axialForceUserFunction=axialForceUserFunction,
)
#use cable with GenerateStraightLineANCFCable(...)
MINI EXAMPLE for ObjectANCFCable2D
1rhoA = 78.
2EA = 1000000.
3EI = 833.3333333333333
4cable = Cable2D(physicsMassPerLength=rhoA,
5 physicsBendingStiffness=EI,
6 physicsAxialStiffness=EA,
7 )
8
9ancf=GenerateStraightLineANCFCable2D(mbs=mbs,
10 positionOfNode0=[0,0,0], positionOfNode1=[2,0,0],
11 numberOfElements=32, #converged to 4 digits
12 cableTemplate=cable, #this defines the beam element properties
13 massProportionalLoad = [0,-9.81,0],
14 fixedConstraintsNode0 = [1,1,0,1], #add constraints for pos and rot (r'_y)
15 )
16lastNode = ancf[0][-1]
17
18#assemble and solve system for default parameters
19mbs.Assemble()
20mbs.SolveStatic()
21
22#check result
23exudynTestGlobals.testResult = mbs.GetNodeOutput(lastNode, exu.OutputVariableType.Displacement)[0]
24#ux=-0.5013058140308901
Relevant Examples and TestModels with weblink:
sliderCrank3DwithANCFbeltDrive2.py (Examples/), ALEANCFpipe.py (Examples/), ANCFcantileverTest.py (Examples/), ANCFcantileverTestDyn.py (Examples/), ANCFcontactCircle.py (Examples/), ANCFcontactCircle2.py (Examples/), ANCFmovingRigidbody.py (Examples/), ANCFrotatingCable2D.py (Examples/), ANCFslidingJoint2D.py (Examples/), ANCFslidingJoint2Drigid.py (Examples/), ANCFswitchingSlidingJoint2D.py (Examples/), ANCFtestHalfcircle.py (Examples/), ANCFcontactCircleTest.py (TestModels/), ANCFcontactFrictionTest.py (TestModels/), computeODE2EigenvaluesTest.py (TestModels/)
The web version may not be complete. For details, consider also the Exudyn PDF documentation : theDoc.pdf