Continuation of equilibria and stability of slender elastic rods using an asymptotic numerical method

TLDR: In this paper, the 3D kinematics of a rod are treated in a geometrically exact way by parameterizing the position of the centerline and making use of quaternions to represent the orientation of the material frame.

Abstract: We present a theoretical and numerical framework to compute bifurcations of equilibria and stability of slender elastic rods. The 3D kinematics of the rod is treated in a geometrically exact way by parameterizing the position of the centerline and making use of quaternions to represent the orientation of the material frame. The equilibrium equations and the stability of their solutions are derived from the mechanical energy which takes into account the contributions due to internal moments (bending and twist), external forces and torques. Our use of quaternions allows for the equilibrium equations to be written in a quadratic form and solved efficiently with an asymptotic numerical continuation method. This finite element perturbation method gives interactive access to semi-analytical equilibrium branches, in contrast with the individual solution points obtained from classical minimization or predictor–corrector techniques. By way of example, we apply our numerics to address the specific problem of a naturally curved and heavy rod under extreme twisting and perform a detailed comparison against our own precision model experiments of this system. Excellent quantitative agreement is found between experiments and simulations for the underlying 3D buckling instabilities and the characterization of the resulting complex configurations. We believe that our framework is a powerful alternative to other methods for the computation of nonlinear equilibrium 3D shapes of rods in practical scenarios.

Figures

Fig. 8. Top view of various equilibrium configurations of a straight elastic rod (κ̂1 ¼ 0 m-1) for increasing values of the rotation angle Φ. The experimental pictures have a black background and the simulations have a white background. The simulation results are rendered to visualize twist by using bi-color rods. (a) Planar equilibrium shape at rest for Φ¼ 01, (b) out-of-plane configuration for Φ¼ 7201, (c) out-of-plane configuration for Φ¼ 14701 and (d) onset of formation of a plectoneme at the middle of the rod.

Fig. 9. Bifurcation diagram of the straight elastic rod (κ̂1 ¼ 0 m−1). (a) Evolution of the maximum transverse displacement, Xmax ¼maxðabsðXÞÞ (defined in Fig. 8(a)), with the control parameter Φ. Comparison between experimental results and the semi-analytical branches computed with MANlab. Solid/dashed lines represent stable/unstable branches, respectively and (b) evolution of the first eigenvalue of the stability problem with control parameter Φ. The critical angle corresponding to the emergence of plectoneme is Φcsim ¼ 20401 for the simulations and Φcexp ¼ 20257151 for the experiments.

Fig. 10. Top view of various equilibrium configurations of a curvy elastic rod (κ̂1ðsÞ ¼ 44:84 m−1) for increasing values of the rotating angle, Φ. The experimental pictures have a black background and the simulations have a white background. The simulation results are rendered to visualize twist by using bi-color rods. (a) Asymmetric out-of-plane equilibrium shape for Φ¼ 1501, (b) one-twist-per-wave configuration for Φ¼ 9001 with three wavelengths between the two clamps, (c) one-twist-per-wave configuration for Φ¼ 18001 with five wavelengths between the two clamps and (d) onset of formation of a plectoneme state at one extremity of the rod.

Fig. 3. Finite-difference method. (a) The centerline of the rod is discretized into N elements separated by N+1 nodes, rj . We also consider Nmaterial frames RðqjÞ to express the orientation of the jth element and (b) the change of orientation between two successive elements j and j+1 is expressed by the N−1 discrete material curvatures at the interconnected nodes: (b1) κj1 around the director ðd j 1 þ djþ11 Þ=2, (b2) κj2 around the director ðd j 2 þ djþ12 Þ=2, and (b3) κj3 around the director ðdj3 þ djþ13 Þ=2.

Fig. 4. (a) The path parameter a is identified as the projection of the solution branch on the tangent vector. This projection is normalized by the length of the tangent vector ½uð1Þ λð1Þ which is set to unity and (b) behavior of a power series close to the radius of convergence. The considered nonlinear unidimensional equation is f ðuðaÞÞ ¼ uð1þ aÞ−1¼ 0. Its unique solution uðaÞ ¼ 1=ð1þ aÞ, represented as solid line, can be expanded asymptotically as uðaÞ≈1−aþ a2 þ⋯þ ð−1Þmam þ ð−1Þmþ1amþ1. The asymptotic expansions for m¼5, 10, 15 and 20 are represented as different lines. We see that u(a) can be approximated up to the radius of convergence a¼1. Applying our ANM method, the step length calculated through Eqs. (59)–(61) would give amax ¼ ε1=ðmþ1Þ ¼ 0:4642 with ε¼ 1 10−7 and m¼20.

Fig. 6. The writhing experiment. A L¼30 cm long elastic rod with a circular cross section of radius 1.55 mm, Young 's modulus E¼1300760 KPa, volumic mass ρ¼ 1200 kg=m3 and a custommade constant natural curvature κ̂1 is fixed at both ends between two concentrically aligned drill chucks separated by a distance d¼22 cm. Our experiment consists of quasi-statically increasing the rotation Φ at one end and investigating the evolution of equilibrium state with the control parameter Φ. Inset courtesy of Khalid Jawed: Fabrication process of an elastomeric rod. During working time, the PVC tubes containing the silicone-based elastomer are wound around cylinders of external radius Re and left in that position for 24 h. After demolding, it confers to the rod a constant natural curvature κ̂1ðsÞ≈1=Re .

Frequently asked questions

Q1. What contributions have the authors mentioned in the paper "Continuation of equilibria and stability of slender elastic rods using an asymptotic numerical method" ?

In this paper, an asymptotic numerical method was proposed to solve the problem of modeling the mechanics of thin elastic rods. 

Q2. What is the main motivation for the determination of the local stability of equilibrium branches?

Determining the local stability of equilibrium branches is crucial for the physical understanding of the mechanical behavior of slender elastic rods, one of the main motivations being that locally unstable branches cannot be observedexperimentally, and must therefore be classified. 

Q3. What is the role of Lxe in second-order conditions?

When restricted to the subspace M that is tangent to the constraint surface and which the authors denote by LM , LðxeÞ plays the role in second-order conditions directly analogous to that of the Hessian of the objective function in the unconstrained case (Luenberger, 1973; Luenberger and Ye, 2008). 

Q4. What are the kinematic parameters that are not included in the general formulation of f?

In addition to the mechanical parameters which are general to thin rods, the authors also need to account for the kinematic boundary conditions and control parameter Φ, specific to the writhing experiment, which are not included in the general formulation of f ðuÞ given in Eq. (A.1). 

Q5. What is the solution of the nonlinear algebraic equations Eqs. (43)?

The vector xe is a solution of the equilibrium equations Eqs. (43) and (44) for λ¼ λe that satisfies the functional geometrical constraints given in Eq. (45). 

Q6. What is the constant of the internal moment projected in the quaternion basis?

T is the internal moment projected in the quaternion basis defined as a linear superposition of the internal moments due to elementary modes of deformation. 

Q7. What is the simplest way to solve the nonlinear elastic problem?

coupled with traditional predictor–corrector methods, one should be able to continue, step-by-step, the solutions of this nonlinear elastic problem in terms of given geometric or mechanical control parameters (Crisfield, 1991; Doedel, 1981). 

Q8. How do the authors compute the up according to the power series expansion of f a?

The authors have computed the up according to the power series expansion of f ðaÞ given in Eq. (52) so that the norm of f is zero up to the truncation order m. 

Q9. What is the vector of the arbitrary perturbations of the rotational degrees of freedom q?

T is the vector of the arbitrary perturbations of the rotational degrees of freedom q. Upon integration by parts, the authors transform Eq. (16) into an integral that depends on δq alone to arrive atδEe ¼ ∑ 

Q10. What is the definition of the 3D kinematics formulation of the inexten?

The 3D kinematics formulation of their inextensible and unshearable elastic rod is not yet complete because the authors won't be able to derive the equilibrium equations directly from the material curvatures. 

Q11. What is the definition of a complete solution branch?

A complete solution branch is therefore constructed as a succession of semi-analytical portions in the form of Eq. (50), whose length is automatically determined through the estimation of the convergence radius of each power series as sketched in Fig. 5(a). 

Q12. What are the corresponding degrees of freedom and Lagrange multipliers?

Let xe, αe and μe be the vectors of degrees of freedom and Lagrange multipliers respectively, associated with the vector solution ue.