Low-Mach-number and slenderness limit for elastic Cosserat rods and its numerical investigation

2.1. Special Cosserat rod theory

An elastic thread is a slender long body, i.e., a rod in three-dimensional continuum mechanics. Because of its geometry the dynamics might be reduced to a one-dimensional description by averaging the underlying balance laws over its cross-sections. This procedure is based on the assumption that the displacement field in each cross-section can be expressed in terms of a finite number of vector- and tensor-valued quantities. The special Cosserat rod theory [2] consists of only two constitutive elements in the three-dimensional Euclidean space E 3 , a curve r : D → E 3 specifying the position and an orthonormal director triad { d 1 , d 2 , d 3 } : D → E 3 characterizing the orientation of the cross-sections. In D = { ( s , t ) ∈ R 2 ∣ s ∈ [ s a , s b ] , t > 0 } , the parameter s denotes the material cross-section (material point) and t the time. The rod system involves four kinematic/geometric and two dynamic equations (balance laws for linear and angular momentum) (2.1) ∂ t r = v , ∂ s r = τ ∂ t d k = ω × d k , ∂ s d k = κ × d k , k = 1 , 2 , 3 ∂ t ( ( ρ A ) v ) = ∂ s n + f ∂ t ( ( ρ J ) · ω ) = ∂ s m + τ × n + l with tangent τ , generalized curvature κ , linear velocity v and angular velocity ω . The mass line density ( ρ A ) is time-independent for materially closed systems, but modeled as time-dependent for applications, such as evaporation and aggregation. The tensor-valued moment of inertia of the cross-sections ( ρ J ) depends on the configuration of the triad and is hence always time-dependent. To close the system of equations we need to specify the external loads f and l, boundary and initial conditions as well as elastic material laws for the contact force n and couple m. In this paper we consider time-independent mass properties of the rod, i.e., ∂ t ( ρ A ) = 0 and ∂ t ( d i · ( ρ J ) · d j ) = 0 for i , j = 1 , 2 , 3 , and neglect external couples, i.e., l = 0 , for reasons of a simple presentation. However, extensions are straightforward possible. Alternatively to (2.1), the rod system can be also set up by using the compatibility relations between the kinematics and geometry for the curve and for the triad. These two compatibility conditions ∂ t τ = ∂ s v , ∂ t κ = ∂ s ω + ω × κ replace then either the two kinematic equations or the two geometric equations.

The model variant (T) yields a hyperbolic system for a dynamic elastic rod, whereas (S) is very suitable for the transition to a stationary consideration. Presupposing hyper-elastic constitutive relations and external potential loads, the variant (M) is known as bi-Hamiltonian form in literature [13,38]. The choice of the formulation can affect the numerical simulations as we will comment on in Section 3 and the Appendix. However, the distinction plays no role for the asymptotics, thus we make use of the compact notation (M–T–S) in the following.

For the objective formulation of the material laws it is useful to rewrite the rod model in the director basis. To an arbitrary vector field z = ∑ i = 1 3 z i d i = ∑ i = 1 3 z ˘ i e i ∈ E 3 , we indicate the coordinate tuples corresponding to the director basis { d 1 , d 2 , d 3 } and to a fixed outer Cartesian basis { e 1 , e 2 , e 3 } by z = ( z 1 , z 2 , z 3 ) ∈ R 3 and z ˘ = ( z ˘ 1 , z ˘ 2 , z ˘ 3 ) ∈ R 3 , respectively. The director basis can be transformed into the outer basis by the tensor-valued rotation D, i.e., D = e i ⊗ d i = D i j e i ⊗ e j ∈ E 3 ⊗ E 3 with the associated orthogonal matrix D = ( D i j ) = ( d i · e j ) ∈ S O ( 3 ) . For the coordinates, the relation D · z ˘ = z holds – as well as D · ∂ t z ˘ = ∂ t z + ω × z and D · ∂ s z ˘ = ∂ s z + κ × z . The rotation matrix D can be parametrized, for example, in Euler angles or unit quaternions. Moreover, canonical basis vectors in R 3 are denoted by e i , i = 1 , 2 , 3 , e.g., e 1 = ( 1 , 0 , 0 ) . In the director basis, the Cosserat rod model has the following form (M–T–S) (2.2) D · ∂ t r ˘ = v , D · ∂ s r ˘ = τ ∂ t D = − ω × D , ∂ s D = − κ × D ∂ t τ = ∂ s v + κ × v − ω × τ ∂ t κ = ∂ s ω + κ × ω ( ρ A ) ∂ t v = ∂ s n + κ × n − ω × ( ρ A ) v + D · f ˘ ( ρ J ) · ∂ t ω = ∂ s m + κ × m + τ × n − ω × ( ( ρ J ) · ω ) with general elastic material laws n ( s , t ) = N ( τ ( s , t ) , κ ( s , t ) , s ) , m ( s , t ) = M ( τ ( s , t ) , κ ( s , t ) , s ) . Following the assumption on the mass properties, ( ρ J ) is here the time-independent representation of the inertia tensor in the director basis. In the stated general form the objective elastic constitutive functions N and M depend on the strain variables (coordinate tuples) τ and κ whose components measure shear τ 1 , τ 2 , dilatation τ 3 , flexure κ 1 , κ 2 and torsion κ 3 . The subsequent asymptotic analysis is valid for material laws that possess the following properties.

Assumption 2 (Properties of the material law N ).

Assume that the material law N for the contact forces fulfills the following conditions:

N ( τ , κ , s ) = 0 holds if and only if τ = τ ˚ ( s ) .

Proceeding from the Cosserat rod model, we derive beam equations of Kirchhoff-type (a generalized string model) as the combined slenderness and low-Mach-number limits in this subsection. We consider an elastic thread. It has various physical properties for which we choose typical constant reference values that we mark by the subscript ⋆ , i.e., mass line density ( ρ A ) ⋆ , moment of inertia ( ρ I ) ⋆ , length L ⋆ , magnitude of mean velocity V ⋆ and contact force N ⋆ . When making (2.2) dimensionless, we can characterize the dynamics of the thread by help of two dimensionless numbers: the slenderness parameter ϵ that is the ratio between the thread’s diameter and its length as well as the (typical) Mach number Ma that is the ratio between the thread’s velocity and the speed of sound ϵ = 1 L ⋆ ( ρ I ) ⋆ ( ρ A ) ⋆ , Ma = V ⋆ ( ρ A ) ⋆ N ⋆ . The speed of sound is here determined by the typical contact force N ⋆ . In the classical sense of compression and shear waves, an elastic thread has certainly different speeds of sounds. Hence, the applied scaling presupposes that the associated Mach numbers behave similar. This means in the special case of a linear isotropic material (cf. Example 4) that the typical Young’s modulus could be chosen, N ⋆ = ( EA ) ⋆ , and that the Poisson ratio satisfies ν = const in the asymptotics.

For the scaling, we use the following reference values s 0 = r 0 = L ⋆ , v 0 = V ⋆ , t 0 = L ⋆ V ⋆ , ω 0 = V ⋆ L ⋆ , κ 0 = 1 L ⋆ ( ρ A ) 0 = ( ρ A ) ⋆ , ( ρ J ) 0 = ( ρ I ) ⋆ , n 0 = ( ρ A ) ⋆ V ⋆ 2 , m 0 = ( ρ A ) ⋆ L ⋆ V ⋆ 2 N 0 = N ⋆ , M 0 = ( ρ I ) ⋆ N ⋆ ( ρ A ) ⋆ L ⋆ , f 0 = ( ρ A ) ⋆ V ⋆ 2 L ⋆ and introduce to every dimensional variable z the associated dimensionless one z ¯ as z ( s 0 s ¯ , t 0 t ¯ ) = z 0 z ¯ ( s ¯ , t ¯ ) , s = s 0 s ¯ , t = t 0 t ¯ . Skipping the superscript ¯ for readability, the dimensionless rod system is then given by (2.3) D · ∂ t r ˘ = v , D · ∂ s r ˘ = τ ∂ t D = − ω × D , ∂ s D = − κ × D ∂ t τ = ∂ s v + κ × v + τ × ω ∂ t κ = ∂ s ω + κ × ω ( ρ A ) ∂ t v = ∂ s n + κ × n + ( ρ A ) v × ω + D · f ˘ ϵ 2 ( ρ J ) · ∂ t ω = ∂ s m + κ × m + τ × n + ϵ 2 ( ( ρ J ) · ω ) × ω with general elastic material laws ( N satisfying Assumption 2) n ( s , t ) = 1 Ma 2 N ( τ ( s , t ) , κ ( s , t ) , s ) , m ( s , t ) = ϵ 2 Ma 2 M ( τ ( s , t ) , κ ( s , t ) , s ) .

To establish an asymptotic relation between the Cosserat rod (2.3) and a beam of Kirchhoff-type, we consider the limits ϵ → 0 , Ma → 0 with Ma ϵ = μ , μ = const > 0 (i.e., asymptotic limit along straight lines in the ( ϵ , Ma ) -plane with slope μ). As the small parameter ϵ appears only in second power in the equations ( ϵ 2 , Ma 2 = μ 2 ϵ 2 ), we expand all quantities in an even power series with respect to the slenderness parameter, i.e., z = ∑ i = 0 ∞ ϵ 2 i z ( i ) . The contact force n ( s , t ) = ( μ ϵ ) − 2 N ( τ ( s , t ) , κ ( s , t ) , s ) implies 0 = N ( τ ( 0 ) ( s , t ) , κ ( 0 ) ( s , t ) , s ) , as ϵ = 0 , hence τ ( 0 ) ( s , t ) = τ ˚ ( s ) according to Assumption 2a). Consequently, we set the strain to be τ = τ ˚ + ϵ 2 ξ with a function ξ ∼ O ( 1 ) . Then the contact force becomes n ( s , t ) = ( μ ϵ ) − 2 N ( τ ˚ ( s ) + ϵ 2 ξ ( s , t ) , κ ( s , t ) , s ) , whose leading order we obtain via Taylor expansion n ( 0 ) ( s , t ) = μ − 2 ( ∂ τ N ( τ ˚ ( s ) , κ ( 0 ) ( s , t ) , s ) · ξ ( 0 ) ( s , t ) + ∂ κ N ( τ ˚ ( s ) , κ ( 0 ) ( s , t ) , s ) · κ ( 1 ) ( s , t ) ) . Here, the second term vanishes since N ( τ ˚ ( s ) , κ , s ) = 0 for all κ and s (Assumption 2). Introducing L = ∂ τ N ( τ ˚ ( · ) , κ , · ) we find ξ = μ 2 L − 1 · n + O ( ϵ 2 ) . Hence, the expression τ = τ ˚ + ϵ 2 μ 2 L − 1 · n + O ( ϵ 4 ) of the strain is exact up to an error of O ( ϵ 4 ) . Inserting this expression into (2.3) yields consequently a ϵ-dependent consistent system up to order O ( ϵ 4 ) , i.e., (2.4) D · ∂ t r ˘ = v , D · ∂ s r ˘ = τ ˚ + ϵ 2 μ 2 L − 1 · n ∂ t D = − ω × D , ∂ s D = − κ × D ϵ 2 μ 2 L − 1 · ∂ t n = ∂ s v + κ × v + τ ˚ × ω + ϵ 2 μ 2 L − 1 · n × ω ∂ t κ = ∂ s ω + κ × ω ( ρ A ) ∂ t v = ∂ s n + κ × n + ( ρ A ) v × ω + D · f ˘ ϵ 2 ( ρ J ) · ∂ t ω = ∂ s m + κ × m + τ ˚ × n + ϵ 2 ( μ 2 L − 1 · n × n + ( ρ J ) · ω × ω ) with m ( s , t ) = μ − 2 M ( τ ˚ ( s ) + ϵ 2 μ 2 L − 1 · n ( s , t ) , κ ( s , t ) , s ) , L ( s , t ) = ∂ τ N ( τ ˚ ( s ) , κ ( s , t ) , s ) . Setting ϵ = 0 in (2.4), we obtain the low-Mach-number–slenderness limit for elastic bodies. The limit equations (2.5) describe a simplified model of the special Cosserat rod theory, where the angular inertia terms vanish and instead of a material law for the contact force the time-independence of the strain field is imposed ( ∂ t τ = 0 , i.e., τ = τ ˚ ). This property enforces the limit beam to be inextensible and unshearable and is known as Kirchhoff constraint. Classically, the Kirchhoff constraint relates the tangent and the director triad τ = d 3 , it is often stated as D · ∂ s r ˘ = τ = e 3 in the director basis, see, e.g., [2,10,25]. This corresponds to arc-length parametrized rod curves with perpendicular cross-sections as choice for stress-free configurations, τ ˚ = e 3 . (2.5) D · ∂ t r ˘ = v , D · ∂ s r ˘ = τ ˚ ∂ t D = − ω × D , ∂ s D = − κ × D 0 = ∂ s v + κ × v + τ ˚ × ω ∂ t κ = ∂ s ω + κ × ω ( ρ A ) ∂ t v = ∂ s n + κ × n + ( ρ A ) v × ω + D · f ˘ 0 = ∂ s m + κ × m + τ ˚ × n with m ( s , t ) = μ − 2 M ( τ ˚ ( s ) , κ ( s , t ) , s ) .

In the formulation of (2.4) the contact force n and the strain τ change their roles, their proportionality operator L is thereby determined by the material law N . In the limit ϵ = 0 in (2.5), n becomes a variable to the constraint D · ∂ s r ˘ = τ ˚ , whereas the material law decouples and is not longer needed for the determination of the solution. The structure of (2.4) obviously changes from a hyperbolic to a degenerate differential-hyperbolic-like system with constraint as ϵ → 0 , the system becomes stiff. Note that L can here be identified with the moduli of the linear elasticity theory, linear operator ( EA ) , cf. Example 4. But in contrast to the Timoshenko beam whose validity is restricted to small deformations, system (2.4) allows for large deformation in case of small ϵ as N can be chosen arbitrarily and L is just its linearization, L = ∂ τ N ( τ ˚ ( · ) , · ) . Without loss of generality, we might thus restrict our considerations for small ϵ to an (affine) linear material law of the form N ( τ , κ , s ) = L ( s ) · ( τ − τ ˚ ( s ) ) which implies an explicit and easily invertible relation between contact force and strain. This relation offers advantages in setting up a numerical scheme that is applicable to both, limit and ϵ-dependent system as ϵ → 0 .

The presented low-Mach-number–slenderness limit finds its analogue in fluid dynamics where the low-Mach-number limit represents the incompressible limit for compressible fluids [1,27]. Here, in elasticity, the asymptotic derivation goes with slow dynamics and slenderness. The limit is valid for arbitrary elastic material laws M for the contact couple. When the limit system (2.5) is equipped with the Euler–Bernoulli relation it describes a Kirchhoff beam. In that case with τ ˚ = e 3 the system is also known as Kirchhoff–Love equations [23]. However, in the classical theory the terminology Kirchhoff beam is much wider and stands for an inextensible and unshearable rod (kinetic analogue [2]). The vanishing of the angular inertia terms is generally not presupposed for a Kirchhoff beam (see, e.g., [13,14]), but results here as consequence of the chosen scaling and the associated asymptotics.

2.3. Asymptotic framework with Euler–Bernoulli material law

In the further work we focus on the asymptotic framework in a special case, where we restrict on a linear Euler–Bernoulli relation M = ( EJ ) · κ and on an homogeneous thread with circular cross-sections and τ ˚ = e 3 . Note that these are just technical simplifications to facilitate the numerical studies. They can be easily dropped if it is relevant for certain applications.

In this case the dimensionless quantities become ( ρ A ) = 1 , ( ρ J ) = P 2 , L = ( EA ) = a − 1 P a , ( EJ ) = P 2 / a where P k = diag ( 1 , 1 , k ) , k ∈ R and a = 2 ( 1 + ν ) with Poisson ratio ν ∈ [ 0 , 0.5 )  – in accordance to Example 4. The physical quantities stated in Example 4 correspond to the reference ⋆ -values chosen for the scaling. Hence, we deal with the following system in the (M–T–S) notation that need to be supplemented with appropriate, consistent initial and boundary conditions, it describes the Cosserat rod for ϵ > 0 and the Kirchhoff beam for ϵ = 0 , cf. (2.4), (2.5) (2.6) D · ∂ t r ˘ = v , D · ∂ s r ˘ = e 3 + ϵ 2 μ 2 a P 1 / a · n ∂ t D = − ω × D , ∂ s D = − κ × D ϵ 2 μ 2 a P 1 / a · ∂ t n = ∂ s v + κ × v + e 3 × ω + ϵ 2 μ 2 a ( P 1 / a · n ) × ω ∂ t κ = ∂ s ω + κ × ω ∂ t v = ∂ s n + κ × n + v × ω + D · f ˘ ϵ 2 P 2 · ∂ t ω = μ − 2 ( P 2 / a · ∂ s κ + κ × ( P 2 / a · κ ) ) + e 3 × n ϵ 2 P 2 · ∂ t ω = + ϵ 2 ( μ 2 a ( P 1 / a · n ) × n + ( P 2 · ω ) × ω )

The asymptotic framework provides a special structure of the equations that is exploited for the numerical handling. Let Φ denote the vector-valued function comprising all system variables, Φ ( s , t ) ∈ R m , then the ϵ-dependent partial differential algebraic system (2.6) can be summarized for all model variants (M), (T), (S) in the form (2.8) A ϵ · ∂ t Φ + B ϵ · ∂ s Φ + c ϵ ( Φ ) = 0 , where A ϵ and B ϵ are – possibly singular – matrices with constant coefficients and c ϵ is a vector-valued nonlinear function in Φ. Their ϵ-dependence can be expressed as A ϵ = A ( 0 ) + ϵ 2 A ( 1 ) , analogously for B ϵ and c ϵ . Expanding Φ in the even power series Φ = ∑ i = 0 n ϵ 2 i Φ ( i ) and plugging it in (2.8) – as before –, we find the Kirchhoff beam in O ( 1 ) and its first-order correction in O ( ϵ 2 ) (2.9) A ( 0 ) · ∂ t Φ ( 0 ) + B ( 0 ) · ∂ s Φ ( 0 ) + c ( 0 ) ( Φ ( 0 ) ) = 0 , A ( 0 ) · ∂ t Φ ( 1 ) + B ( 0 ) · ∂ s Φ ( 1 ) + ∂ Φ c ( 0 ) ( Φ ( 0 ) ) · Φ ( 1 ) = f [ Φ ( 0 ) ] , (2.10) with − f [ Φ 0 ] = A ( 1 ) · ∂ t Φ ( 0 ) + B ( 1 ) · ∂ s Φ ( 0 ) + c ( 1 ) ( Φ ( 0 ) ) where ∂ Φ c ( 0 ) denotes the Jacobi matrix of c ( 0 ) . Both asymptotic systems together (2.9)–(2.10) are exact up to order O ( ϵ 4 ) . As usual for asymptotic considerations, the systems have a similar equation structure with the same system matrices A ( 0 ) and B ( 0 ) , they just differ in the right-hand-side and the dependence on the variable. The first-order correction (2.10) is trivially linear in the variable, whereas the limit system (2.9) keeps the nonlinearity of the original ϵ-dependent system (2.8). Having a Newton method for the numerical treatment of the nonlinear term c ( 0 ) ( Φ ( 0 ) ) in mind, we need its Jacobi matrix already for the computation of the limit system. So, the assembly of the linear system matrix associated with the first-order correction is for free.

Remark 6 (Conservation of energy).

The Cosserat rod model as well as the Kirchhoff beam equations in (2.6) are energy-conserving for external potential forces [38], presupposing classical boundary conditions such as, for example, a cantilever beam or a beam with stress-free ends. This holds also true for the system associated with the first-order correction. The corresponding energy with (external) potential energy V is given by w ϵ = w ( 0 ) + ϵ 2 w ( 1 ) + O ( ϵ 4 ) w ( 0 ) = 1 2 ( v ( 0 ) ) 2 + 1 2 μ 2 κ ( 0 ) · P 2 / a · κ ( 0 ) + V ( r ˘ ( 0 ) ) w ( 1 ) = v ( 0 ) · v ( 1 ) + 1 μ 2 κ ( 0 ) · P 2 / a · κ ( 1 ) + 1 2 ω ( 0 ) · P 2 · ω ( 0 ) + μ 2 a 2 n ( 0 ) · P 1 / a · n ( 0 ) + V ( r ˘ ( 1 ) ) .

Let { e 1 , e 2 , e 3 } be the outer Cartesian basis. Consider a thread fixed at one end ( s = 0 ) and stress-free at the other end ( s = 1 ) that is initially static, stress-free and straight in direction of e 1 . Let it be exposed to transversal oscillations in the e 1 - e 2 -plane due to gravity f = − ( ρ A ) g e 2 such that d 2 = e 3 during its motion (cf. Fig. 1). For this two-dimensional scenario, the model system (2.6) in the (M–T–S) notation simplifies to D ( α ) · ∂ t r ˘ = v , D ( α ) · ∂ s r ˘ = e 2 + ϵ 2 μ 2 a P 1 / a · n ∂ t α = ω , ∂ s α = κ ϵ 2 μ 2 a P 1 / a · ∂ t n = ∂ s v − κ v ⊥ − ω e 1 + ϵ 2 μ 2 ω P a · n ⊥ ∂ t κ = ∂ s ω ∂ t v = ∂ s n − κ n ⊥ + ω v ⊥ + D ( α ) · f ˘ ϵ 2 ∂ t ω = μ − 2 ∂ s κ + n 1 + ϵ 2 μ 2 ( 1 − a ) n 1 n 3 with the dimensionless force f ˘ = − Fr − 2 e 2 and the respective initial and boundary conditions r ˘ ( s , 0 ) = s e 1 , α ( s , 0 ) = κ ( s , 0 ) = ω ( s , 0 ) = 0 , n ( s , 0 ) = v ( s , 0 ) = 0 r ˘ ( 0 , t ) = v ( 0 , t ) = 0 , α ( 0 , t ) = ω ( 0 , t ) = 0 , n ( 1 , t ) = 0 , κ ( 1 , t ) = 0 . Here, we use the abbreviations r ˘ = ( r ˘ 1 , r ˘ 2 ) , f ˘ = ( f ˘ 1 , f ˘ 2 ) , v = ( v 1 , v 3 ) , n = ( n 1 , n 3 ) , ω = ω 2 , κ = κ 2 , m = m 2 , P k = diag ( 1 , k ) , k = R . Moreover, z ⊥ = ( − z 3 , z 1 ) denotes the tuple perpendicular to z for z ∈ { v , n } . The rotation matrix D is expressed in terms of the single angle α D ( α ) = − sin α cos α cos α sin α , α = ∠ ( e 1 , d 3 ) . The Froude number Fr represents the ratio of inertia and gravity.

Because the asymptotic results turn out not to depend very sensitively on the specific choice of parameters, we use the setting ( Fr , μ , a ) = ( 1 , 10 , 2.5 ) with the end time T = 2.5 as example in the following studies. This parameter setting is characteristic in the context of non-woven manufacturing [21,32]. Typical properties of polymer fibers are d = 10 − 5 m (diameter), L = 10 − 1 m , ρ = 10 3 kg/m 3 , E = 10 10 kg/ ( m s 2 ) , ν = 0.25 and V = 1 m/s , implying ϵ ∼ O ( 10 − 4 ) . The end time T is chosen respectively with regard to stability issues (see Remark 7). The dynamics of the cantilever beam under gravity is illustrated for ϵ = 0 and ϵ = 0.02 in Fig. 1.

Remark 7.

For the planar inextensible beam an existence proof is derived in the absence of gravity and for non-vanishing angular inertia in [7]. Nonlinear planar and non-planar responses of the inextensible beam were investigated numerically in [34] using a combination of the Galerkin procedure and the method of multiple scales, it turned out that the nonlinear geometric terms produce a hardening effect and dominate the non-planar responses for all modes. The inertia terms particularly dominate the high-frequency modes. A global bifurcation analysis for a cantilever beam subjected to a harmonic axial excitation and transverse excitations at the free end was performed in [44]. Taking into account these results, the simulations of our limit system that has no angular inertia are run in the stable planar regime.

3.2. Results

In the asymptotic framework we have Φ ϵ = Φ ( 0 ) + ϵ 2 Φ ( 1 ) + O ( ϵ 4 ) , i.e., the limit system and its first-order correction are analytically exact up to the order O ( ϵ 4 ) , cf. (2.8)–(2.10). The solutions of the limit system Φ ( 0 ) (Kirchhoff beam), its first-order correction Φ ( 1 ) as well as of the original ϵ-dependent system Φ ϵ (Cosserat rod) are numerically approximated by φ ( 0 ) , φ ( 1 ) and φ ϵ . We observe that the numerical approximations inherit the asymptotic relation as desired, this means that (3.1) φ ϵ − φ ( 0 ) = c 1 ⋆ ∼ O ( ϵ 2 ) , c 1 ⋆ = ϵ 2 c 1 (3.2) φ ϵ − φ ( 0 ) − ϵ 2 φ ( 1 ) = c 2 ⋆ ∼ O ( ϵ 4 ) , c 2 ⋆ = ϵ 4 c 2 hold true. Figure 2 shows the L 2 ( 0 , 1 ) -norm of the vector-valued functions c i ⋆ and c i , i = 1 , 2 , in dependence on ϵ for ϵ ∈ [ 10 − 10 , 1 ] . The quantities are computed at time T = 2 for the different model variants (M), (T) and (S) introduced in Notation 1 by applying the numerical scheme with the parameters Δ t = Δ s = 10 − 2 and λ = 1 . Comparing the model variants, they all yield the quadratic convergence for the first term c 1 ⋆ with the same (ϵ-independent) magnitude ‖ c 1 ‖ in the range ϵ ∈ [ 10 − 8 , 10 − 2 ] , its relative deviation from ‖ φ ( 1 ) ‖ lies below 10 − 3 . For the second term c 2 ⋆ the model variants show the quartic convergence in the range ϵ ∈ [ 10 − 4 , 10 − 2 ] and (S) even in ϵ ∈ [ 10 − 6 , 10 − 2 ] . The reduced order of convergence (down to p = 2 in ϵ ∈ [ 10 − 6 , 10 − 4 ] ) for (T) and (M) might be explained by the abandonment of boundary conditions. All boundary conditions are explicitly prescribed and incorporated in the scheme for (S), whereas the boundary conditions associated with r ˘ , α for (T) and to v , n for (M) follow only numerically from the initial conditions and the model equations. The approximations are consistent but the effect of the term ϵ 2 φ ( 1 ) in the numerical differentiation is below the computational accuracy (cancellation error). Also the tails that are observed for very small ϵ in the convergence plots come from numerical noise. The striking switch in the convergence behavior of both terms c i ⋆ , i = 1 , 2 , that occurs at ϵ = 10 − 2 for all model variants can be explained by the asymptotics itself because the asymptotic framework holds true only for ϵ sufficiently small. This can be also clearly seen in Fig. 1 where the beam curve associated with the comparatively large ϵ lies away from the limit curve. Note that the curves would be indistinguishable for ϵ < 10 − 2 . The corresponding peaks observed at ϵ = 10 − 2 in Fig. 2 (right) leave freedom for speculations: they might be due to the transversal stress component whose effect is then superposed by the other variables.

OPEN IN VIEWER

Fig. 2.

Conservation of the asymptotic relations for the different model variants (cf. (3.1)–(3.2), c i ⋆ = ϵ 2 i c i ). Top: first term ‖ c 1 ⋆ ‖ (left) with magnitude ‖ c 1 ‖ (right) plotted over ϵ. The solid line with p = 2 indicates quadratic convergence. Bottom: second term ‖ c 2 ⋆ ‖ (left) with magnitude ‖ c 2 ‖ (right) plotted over ϵ. The solid lines with p = 4 and p = 2 indicate quartic and quadratic convergence, respectively.

As the Cosserat rod model contains the small asymptotic parameter explicitly, the system of equations becomes stiff and difficult to solve numerically as ϵ → 0 . The asymptotic framework provides a special structure of the asymptotic systems (2.9)–(2.10) that can be exploited in the numerics when solving the discretized equations with a Newton method. The nonlinear original ϵ-dependent system as well as the nonlinear limit system require 2–3 Newton iterations in average per time step for all model variants, supposing Δ t ∼ O ( 10 − 2 ) . The resulting linear systems are of same size and have a similar block-structured band matrix. The Jacobian of the limit system equals the linear system matrix associated with the first-order correction. Therefore, the effort of computing the two asymptotic systems is similar to the effort for the original ϵ-dependent system, as the first-order correction costs only a single linear solving. Moreover, the solving of the asymptotic systems is much more accurate and robust in determining the influence of small ϵ-values. Consequently, the asymptotic framework makes a uniform numerical handling of the transition regime for small ϵ easily possible. The choice of the model variant has only a very marginal effect on the numerical performance. It influences neither efficiency nor convergence properties, for details see the Appendix. However, (T) seems to be preferable for energy-conservation and (S) for robustness and asymptotic-conservation.