Symbolic-Numeric Integration of the Dynamical Cosserat Equations

TL;DR

This paper introduces a symbolic-numeric method for integrating the dynamical Cosserat equations, improving computational efficiency in modeling slender structures like fibers and rods.

Contribution

It combines symbolic and numerical techniques to solve the nonlinear PDE system of Cosserat equations, building on previous solutions for the kinematic part.

Findings

The method effectively integrates the dynamical Cosserat equations.

Experimental results show improved computational efficiency.

Comparison with the generalized alpha-method validates the approach.

Abstract

We devise a symbolic-numeric approach to the integration of the dynamical part of the Cosserat equations, a system of nonlinear partial differential equations describing the mechanical behavior of slender structures, like fibers and rods. This is based on our previous results on the construction of a closed form general solution to the kinematic part of the Cosserat system. Our approach combines methods of numerical exponential integration and symbolic integration of the intermediate system of nonlinear ordinary differential equations describing the dynamics of one of the arbitrary vector-functions in the general solution of the kinematic part in terms of the module of the twist vector-function. We present an experimental comparison with the well-established generalized alpha-method illustrating the computational efficiency of our approach for problems in structural mechanics.

Equations

$${\mathbf{\kappa }}_t={\mathbf{\omega }}_s-\mathbf{\omega }\times \mathbf{\kappa },$$

$${\mathbf{\nu }}_t={\mathbf{v}}_s+\mathbf{\kappa }\times \mathbf{v}-\mathbf{\omega }\times \mathbf{\nu },$$

$$\rho \mathbf{J}\cdot {\mathbf{\omega }}_t={\hat{\mathbf{m}}}_s+\mathbf{\kappa }\times \hat{\mathbf{m}}+\mathbf{\nu }\times \hat{\mathbf{n}}-\mathbf{\omega }\times (\rho \mathbf{J}\cdot \mathbf{\omega })+\mathbf{L},$$

$$\rho A{\mathbf{v}}_t={\hat{\mathbf{n}}}_s+\mathbf{\kappa }\times \hat{\mathbf{n}}-\mathbf{\omega }\times (\rho A\mathbf{v})+\mathbf{F}.$$

$${\partial }_s{\mathbf{d}}_k=\mathbf{\kappa }\times {\mathbf{d}}_k,{\partial }_t{\mathbf{d}}_k=\mathbf{\omega }\times {\mathbf{d}}_k,$$

$$[a,b]\times \mathbb{R}\ni (s,t)\mapsto \left(\mathbf{r}(s,t),{\mathbf{d}}_1(s,t),{\mathbf{d}}_2(s,t),{\mathbf{d}}_3(s,t)\right)\in {\mathbb{E}}^3.$$

$$\mathbf{\nu }:={\partial }_s\mathbf{r}=\sum _{k=1}^3{\nu }_k{\mathbf{d}}_k,\mathbf{v}:={\partial }_t\mathbf{r}=\sum _{k=1}^3{v}_k{\mathbf{d}}_k.$$

$$\displaystyle\begin{array}[]{l}\rho(s)A(s)\partial_{t}\boldsymbol{v}=\partial_{s}\boldsymbol{n}(s,t)+\boldsymbol{F}(s,t)\,,\\[5.69046pt] \partial_{t}\boldsymbol{h}(s,t)=\partial_{s}\boldsymbol{m}(s,t)+\boldsymbol{\nu}(s,t)\times\boldsymbol{n}(s,t)+\boldsymbol{L}(s,t)\,,\end{array}$$

$$\mathbf{m}(s,t)=\hat{\mathbf{m}}\left(\mathbf{\kappa }(s,t),\mathbf{\nu }(s,t),s\right),\mathbf{n}(s,t)=\hat{\mathbf{n}}\left(\mathbf{\kappa }(s,t),\mathbf{\nu }(s,t),s\right).$$

$$\mathbf{\omega }={\mathbf{p}}_t+\frac{p-\sin (p)}{p^3}\left(\mathbf{p}\left(\mathbf{p}\cdot {\mathbf{p}}_t\right)-p^2{\mathbf{p}}_t\right)-\frac{1-\cos (p)}{p^2}\mathbf{p}\times {\mathbf{p}}_t,$$

$$\mathbf{\kappa }={\mathbf{p}}_s+\frac{p-\sin (p)}{p^3}\left(\mathbf{p}\left(\mathbf{p}\cdot {\mathbf{p}}_s\right)-p^2{\mathbf{p}}_s\right)-\frac{1-\cos (p)}{p^2}\mathbf{p}\times {\mathbf{p}}_s,$$

$$\mathbf{\nu }=\mathbf{q}\times \mathbf{\kappa }-{\mathbf{q}}_s,$$

$$\mathbf{v}=\mathbf{q}\times \mathbf{\omega }-{\mathbf{q}}_t,$$

$$p=\sqrt{p_1^2+p_2^2+p_3^2}.$$

$${\mathbf{p}}_t=\frac{\mathbf{p}\cdot \mathbf{\omega }}{p^2}\mathbf{p}+\frac{1}{2}\mathbf{p}\times \mathbf{\omega }-\frac{p}{2}\cot \left(\frac{p}{2}\right)\cdot \frac{\mathbf{p}\times (\mathbf{p}\times \mathbf{\omega })}{p^2}.$$

$$2\frac{\cos (p)-1}{p^2}.$$

$${\mathbf{p}}_t=\frac{\mathbf{p}\cdot \mathbf{\omega }}{p^2}\mathbf{p}+\frac{1}{2}\mathbf{p}\times \mathbf{\omega }-\frac{p}{2}\cot \left(\frac{p}{2}\right)\cdot \frac{\mathbf{p}\times (\mathbf{p}\times \mathbf{\omega })}{p^2},$$

$${\mathbf{q}}_t=\mathbf{q}\times \mathbf{\omega }-\mathbf{v},$$

$$\rho \mathbf{J}\cdot {\mathbf{\omega }}_t={\hat{\mathbf{m}}}_s+\mathbf{\kappa }\times \hat{\mathbf{m}}+\mathbf{\nu }\times \hat{\mathbf{n}}-\mathbf{\omega }\times (\rho \mathbf{J}\cdot \mathbf{\omega })+\mathbf{L},$$

$$\rho A{\mathbf{v}}_t={\hat{\mathbf{n}}}_s+\mathbf{\kappa }\times \hat{\mathbf{n}}-\mathbf{\omega }\times (\rho A\mathbf{v})+\mathbf{F}$$

$${\mathbf{p}}_t\to \frac{\mathbf{p}(s,t+△t)-\mathbf{p}(s,t)}{△t}.$$

$$p(s,t+△t)\in (0,2\pi )$$

$$p_t=\mathbf{e}\cdot \mathbf{\omega },$$

$$2{\mathbf{e}}_t=\mathbf{e}\times \mathbf{\omega }-\cot \left(\frac{p}{2}\right)\mathbf{e}\times (\mathbf{e}\times \mathbf{\omega }).$$

$$\mathbf{e}=A_1{\mathbf{e}}_1+A_2{\mathbf{e}}_2+A_3{\mathbf{e}}_3,\mathbf{\omega }=\omega {\mathbf{e}}_3,\omega :=\sqrt{{\omega }_1^2+{\omega }_2^2+{\omega }_3^2}.$$

$$2(A_1{)}_t=A_2\omega -\cot \left(\frac{p}{2}\right)A_1{A}_3\omega ,$$

$$2(A_2{)}_t=-A_1\omega -\cot \left(\frac{p}{2}\right)A_2{A}_3\omega ,$$

$$2(A_3{)}_t=-\cot \left(\frac{p}{2}\right)(A_3^2-1)\omega ,$$

$$p_t=A_3\omega .$$

$$\frac{2A_3(A_3{)}_t}{A_3^2-1}=-\cot \left(\frac{p}{2}\right)p_t,$$

$$(1-A_3^2){\sin }^2\left(\frac{p}{2}\right)=C,C_t=0.$$

$$4q_t^2={\omega }^2\left(1-q^2-C\right).$$

$$q=\sqrt{1-C}\sin \left(\frac{1}{2}\omega (C_1-t)\right),(C_1{)}_t=0.$$

$$A_1(s,t)=-\frac{\sqrt{C}\cdot sin(\frac{1}{2}\omega (C_2-t))}{\sqrt{{\omega }^2{\cos }^2(\frac{1}{2}\omega (C_1-t))+C{\sin }^2(\frac{1}{2}\omega (C_1-t))}},$$

$$A_2(s,t)=\frac{\sqrt{C}\cdot cos(\frac{1}{2}\omega (C_2-t))}{\sqrt{{\omega }^2{\cos }^2(\frac{1}{2}\omega (C_1-t))+C{\sin }^2(\frac{1}{2}\omega (C_1-t))}},$$

$$A_3(s,t)=\frac{\sqrt{{\omega }^2-C}\cdot cos(\frac{1}{2}\omega (C_1-t))}{\sqrt{{\omega }^2{\cos }^2(\frac{1}{2}\omega (C_1-t))+C{\sin }^2(\frac{1}{2}\omega (C_1-t))}},$$

$$p(s,t)=2\arccos \left(\frac{\sqrt{{\omega }^2-C}\sin (\frac{1}{2}\omega (C_1-t))}{\omega }\right),$$

$$C(s):={\omega }^2\left(1-A_3^2(s,t_0)\right){\sin }^2\left(\frac{p(s,t_0)}{2}\right),$$

Full text

1 Introduction

Deformable-body dynamics can be considered as a subarea of continuous mechanics that studies motion of deformable solids subject to the action of internal and external forces (cf. [11]). The equations describing the dynamics of such solids are nonlinear partial differential equations (PDEs) whose independent variables are three spatial coordinates and time. Given a particular deformable mechanical structure, to describe its dynamics, it is necessary to satisfy these equations at each point of the structure together with appropriate boundary conditions. For a mechanical structure having special geometric properties, it is worthwhile to exploit these properties to develop a simplified but geometrically exact mechanical model of the structure. Classical examples of such models are Cosserat theories of shells and rods; see e.g. [18] and references therein.

Rods are nearly one-dimensional structures whose dynamics can be described by the Cosserat theory of (elastic) rods (cf. [1], Ch. 8; [18], Ch. 5; and the original work [5]). This is a general and geometrically exact dynamical model that takes bending, extension, shear, and torsion into account, as well as rod deformations under external forces and torques. In this context, the dynamics of a rod is described by a governing system of twelve first-order nonlinear partial differential equations (PDEs) with a pair of independent variables $(s,t)$, where $s$ is the arc-length and $t$ the time parameter. In this PDE system, the two kinematic vector equations ((9a)–(9b) in [1], Ch. 8) are parameter free and represent the compatibility conditions for four vector functions $\boldsymbol{\kappa},\boldsymbol{\omega},\boldsymbol{\nu},\boldsymbol{v}$ in $(s,t)$. Whereas the first vector equation only contains two vector functions $\boldsymbol{\kappa}$ and $\boldsymbol{\omega}$, the second one contains all four vector functions $\boldsymbol{\kappa},\boldsymbol{\omega},\boldsymbol{\nu},\boldsymbol{v}$. The remaining two vector equations in the governing system are dynamical equations of motion and include two more dependent vector variables $\hat{\boldsymbol{m}}(s,t)$ and $\hat{\boldsymbol{n}}(s,t)$. Moreover, these dynamical equations contain parameters (or parametric functions of $s$) to characterize the rod and to include the external forces and torques. Studying the dynamics of Cosserat rods has various scientific and industrial applications, for example in civil and mechanical engineering (cf. [2]), microelectronics and robotics (cf. [3]), biophysics (cf. [7] and references therein), and visual computing (cf. [15]).

Because of its inherent stiffness caused by different deformation modes, the treatment of the underlying equations usually requires the application of specific solvers; see e.g. [16]. In order to reduce the computational overhead caused by the stiffness, we employed Lie symmetry based integration methods (cf. [9, 17]) and the theory of completion to involution (cf. [19]) to the two kinematic vector equations (cf. [1], Ch. 8, Eq. (9a)–(9b)) and constructed their general and analytical solution in [13, 14], which depends on two arbitrary vector functions in $(s,t)$.

In this contribution, we exploit the general analytic solution to the kinematic part of the governing Cosserat system constructed in [13, 14] and develop a symbolic-numeric approach to the integration of the dynamical part of the system. Our approach combines the ideas of numerical exponential integration (see e.g. [8, 12] and references therein) and symbolic integration of the intermediate system of nonlinear ordinary differential equations describing the dynamics of the arbitrary vector function in the general solution to the kinematic part in terms of the module of the twist vector function. The symbolic part of the integration is performed by means of Maple. We present an experimental comparison of the computational efficiency of our approach with that of the generalized $\alpha$-method for the numerical integration of problems in structural mechanics (see e.g. [20]).

This paper is organized as follows. In Section 2, we present the governing PDE system in the special Cosserat theory of rods and analytical solution to its kinematic part constructed in [13, 14]. In Section 3, we show first that the (naive) straightforward numerical integration of the dynamical part of Cosserat system has a severe obstacle caused by a singularity in the system. Then we describe a symbolic-numeric method to integrate the dynamical equations based on the ideas of exponential integration and the construction of a closed form analytical solution to the underlying nonlinear dynamical system. In doing so, we show that this symbolic-numeric method is free of the singularity problem. In Section 4, we present an experimental comparison of our method with the generalized $\alpha$-method. Some concluding remarks are given in Section 5.

2 Governing Cosserat Equations and the General Solution of their Kinematic Part

The governing PDE system in the special Cosserat theory of rods (cf. [1, 3, 5, 13, 14]) can be written in the following form:

[TABLE]

Here, the independent variable $t$ denotes the time and another independent variable $s$ the arc-length parameter identifying a material cross-section of the rod, which consists of all material points whose reference positions are on the plane perpendicular to the rod at $s$. The Darboux vector-function $\boldsymbol{\kappa}=\sum_{k=1}^{3}\kappa_{k}\boldsymbol{d}_{k}$ and the twist vector-function $\boldsymbol{\omega}=\sum_{k=1}^{3}\omega_{k}\boldsymbol{d}_{k}$ are determined by the kinematic relations

[TABLE]

where the vectors $\boldsymbol{d}_{1},\,\boldsymbol{d}_{2}$, and $\boldsymbol{d}_{3}:=\boldsymbol{d}_{1}\times\boldsymbol{d}_{2}$ form a right-handed orthonormal moving frame. These vectors are called directors. The use of the triple $(\boldsymbol{d}_{1},\,\boldsymbol{d}_{2},\,\boldsymbol{d}_{3})$ is natural for the intrinsic description of the rod deformation. Moreover, $\boldsymbol{r}$ describes the motion of the rod relative to the fixed frame $(\boldsymbol{e}_{1},\,\boldsymbol{e}_{2},\,\boldsymbol{e}_{3})$. This is illustrated in Figure 1.

In doing so, the motion of a rod is defined by the mapping

[TABLE]

Furthermore, the governing system (1a)–(1d) includes additional vector-valued dependent variables: linear strain $\boldsymbol{\nu}$ of the rod and the velocity $\boldsymbol{v}$ of the material cross-section:

[TABLE]

\contourlength

1.6pt

The components of the strain variables $\boldsymbol{\kappa}$ and $\boldsymbol{\nu}$ describe the deformation of the rod: the flexure with respect to the two major axes of the cross section $(\kappa_{1},\kappa_{2})$, torsion $(\kappa_{3})$, shear $(\nu_{1},\nu_{2})$, and extension $(\nu_{3})$.

The kinematic part of the governing Cosserat system consists of equations (1a)–(1b) ((9a)–(9b) in [1], Ch. 8). The remaining equations (1c)–(1d) ((9c)–(9d) in [1], Ch. 8) make up the dynamical part of the governing equations. For a rod density $\rho(s)$ and cross-section $A(s)$, these equations follow from Newton’s laws of motion:

[TABLE]

where $\boldsymbol{m}(s,t)=\sum_{k=1}^{3}m_{k}(s,t)\,\boldsymbol{d}_{k}(s,t)$ are the contact torques, $\boldsymbol{n}(s,t)=\sum_{k=1}^{3}n_{k}(s,t)\,\boldsymbol{d}_{k}(s,t)$ are the contact forces, $\boldsymbol{h}(s,t)=\sum_{k=1}^{3}h_{k}(s,t)\,\boldsymbol{d}_{k}(s,t)$ are the angular momenta, and $\boldsymbol{F}(s,t)$ and $\boldsymbol{L}(s,t)$ are the external forces and torque densities.

The contact torques $\boldsymbol{m}(s,t)$ and contact forces $\boldsymbol{n}(s,t)$ corresponding to the internal stresses, are related to the extension and shear strains $\boldsymbol{\nu}(s,t)$ as well as to the flexure and torsion strains $\boldsymbol{\kappa}(s,t)$ by the constitutive relations

[TABLE]

Under certain reasonable assumptions (cf. [1, 3, 13]) on the structure of the right-hand sides in (2), they take the form (1c)–(1d) in which $\boldsymbol{J}$ is the inertia tensor of the cross-section per unit length. Unlike the kinematic part, the dynamical part contains parameters characterizing the rod under consideration: $\rho,A$ and $\boldsymbol{J}$ together with the external force $\boldsymbol{F}$ and torque $\boldsymbol{L}$, whereas the kinematic part is parameter free.

In our previous papers [13, 14], by treating the kinematic Cosserat equations (1a)–(1b) with computer algebra aided methods of the modern Lie symmetry analysis (cf. [9, 17]) and the theory of completion of partial differential systems to involution (cf. [19]), we constructed the following closed form of an analytical solution to the kinematic part and proved its generality:

[TABLE]

where $\boldsymbol{p}(s,t)$ and $\boldsymbol{q}(s,t)$ are arbitrary analytic vector functions, and

[TABLE]

For the efficient numerical solving of the dynamical Cosserat equations (1c)–(1d), we use the following fact: the vector equation (3a) uniquely defines $\boldsymbol{p}_{t}$ in terms of $\boldsymbol{p}$ and $\boldsymbol{\omega}$. We formulate this fact as the following statement.

Proposition 1

The temporal derivative $\boldsymbol{p}_{t}$ of the vector function $\boldsymbol{p}$, as a solution of (3a), reads

[TABLE]

Proof

The vector function $\boldsymbol{p}_{t}$ occurs linearly in (3a). In the component form it is a linear system of three equations in three unknowns $(p_{1})_{t},(p_{2})_{t},(p_{3})_{t}$ whose matrix has a non-singular determinant (cf. formula (11) in [14])

[TABLE]

The vector form of $\boldsymbol{p}_{t}$ as a solution to equality (3a) is given by (4). This can be verified either by hand computation or by using the routines of the Maple package VectorCalculus after the substitution of (4) into the right-hand side of (3a) and simplification of the obtained expression to $\boldsymbol{\omega}$.

$\Box$

Instead of system (1a)–(1d) for unknowns $(\boldsymbol{\omega},\boldsymbol{\kappa},\boldsymbol{\nu},\boldsymbol{v})$, we are going to solve the equivalent system

[TABLE]

for unknown vector functions $(\boldsymbol{p},\boldsymbol{q},\boldsymbol{\omega},\boldsymbol{v})$, where $\boldsymbol{\kappa}$ and $\boldsymbol{\nu}$ are given by (3b)–(3c).

3 Symbolic-Numeric Integration Method

3.1 Naive Approach: Explicit Numerical Solving

Suppose we know the values of the vector-functions $\boldsymbol{\omega}$ and $\boldsymbol{p}$ on a time layer $t$. Then $\boldsymbol{p}_{t}$ in (6a) can be approximated by the forward Euler difference

[TABLE]

However, since equation (6a) has a singularity at $p=2\pi$ related with vanishing (5), there is a restriction to the time step $\Delta t$ caused by the condition

[TABLE]

to be held for all values of $s$. This is a severe problem for the numerical solving of the governing Cosserat system, because in the course of solving, one must control the time step at every value of $s$ to keep the values of $p(s,t)$ within the interval indicated in (7). This problem is resolved in the symbolic-numeric integration method described in the next subsection.

3.2 Advanced Approach based on Exponential Integration

To avoid the problem of controlling the condition (7) we use the differential equation (6a) for $\boldsymbol{p}(t)$ and rewrite it in terms of $p$ and the unit vector $\boldsymbol{e}$ where $\boldsymbol{p}=p\boldsymbol{e}$ . It leads to the following differential system:

[TABLE]

Now assume that the vector $\boldsymbol{\omega}$ is independent of $t$ on the time interval $\Delta t$ and choose the Cartesian coordinate system $\boldsymbol{e}_{1},\boldsymbol{e}_{2},\boldsymbol{e}_{3}$ such that $\boldsymbol{e}_{3}||\boldsymbol{\omega}$:

[TABLE]

Then, we obtain the following system of four first-order differential equations:

[TABLE]

From the equations (9c)–(9d) it follows

[TABLE]

and, hence,

[TABLE]

Equation (10) immediately implies the following statement providing fulfillment of (7).

Proposition 2

If $C\neq 0$, then $p(t)\in(0,2\pi)$ for all $t\geq t_{0}$ if $p(t_{0})\in(0,2\pi)$.

If one substitutes $A_{3}=p_{t}/\omega$ from (9c) and replaces $p$ with $q:=\cos\left(\frac{p}{2}\right)$, then (10) takes the form

[TABLE]

This equation is easily solvable by the Maple routine dsolve which outputs four solutions. These solutions can be unified into the general solution

[TABLE]

Then, the whole system (9a)–(9d) admits the following general analytical solution

[TABLE]

where $C,C_{1},C_{2}$ are functions of $s$. These functions are determined by the following initial data:

[TABLE]

Proposition 3

$C(s)\equiv 0$* if and only if $A_{3}(s,t)=\pm 1$ what corresponds to a degenerated solution

[TABLE]

Computationally, this solution is not of interest, since it is unstable: a small (e.g. numerical) deviation of $A_{3}(s,t_{0})$ from $\pm 1$ converts the solution into a generic one.

Proof

The solution

[TABLE]

is singular. Unlike the generic solution, the value of $|p(t)|$ in this solution may increase indefinitely. However, it is unstable, since any small perturbation $\epsilon>0$ to the initial value $A_{3}(s,t_{0})=\pm\left(1-\epsilon(s)\right)$ leads to the generic case when $p(t)$ remains bounded.

$\Box$

Our symbolic-numeric approach to the derivation of equations (9a)–(9d) and the construction of their explicit analytical solution (11a)–(11d) is in accord with the general principles of exponential integration (see e.g. [8]). The basic idea behind exponential integration is the identification of a prototypical differential system which has the stiffness properties similar to those in the original equation and which admits explicit solving.

In our case, the stiffness properties of the differential system (9a)–(9d) are similar to those in system (6a)–(6d). In doing so, the last system belongs to the second class of stiff problems (cf. [8], p. 210) whose stiffness is caused by the highly oscillatory behavior of their solutions. For such problems both explicit and implicit Euler schemes fail to provide the required stability unless the step size is strongly reduced to provide the resolution of all the oscillations in the solution. Thereby, the standard numerical treatment of the equations (1a)–(1d) and hence equations (6a)–(6d) is computationally inefficient. Just by this reason, special numerical solvers have been designed (cf. [10, 20]) for the equations (1a)–(1d).

Figure 2 illustrates the stiffness of the differential system (6a)–(6d). The behavior is shown, at $\omega=1$, of the functions $p(s_{0},t)$ and $A_{3}(s_{0},t)$ as solutions of (9c)–(9d) for the initial conditions $p(s_{0},0)=1$ and $A_{3}(s_{0},0)=0.99$. As illustrated, the solution oscillates and changes drastically over time. A numerical reconstruction of such a behavior is possible only for very small step sizes of difference approximations.

4 Numerical Comparison with the Generalized $\alpha$-Method

In 1993, Chung and Hulbert (cf. [4]) presented the generalized $\alpha$-method as a new integration algorithm for problems from structural mechanics. It is characterized primarily by a controllable numerical dissipation of high-frequency components in the numerical solution. These occur, for example in the context of finite element-based simulations, when the high-frequency states are too roughly resolved. Such methods usually improve the convergence behavior of iterative solving strategies for nonlinear problems.

The generalized $\alpha$-method is well-established in the field of structural mechanics and has the major advantage of unconditional stability as well as user controllable numerical damping. The idea of the introduction of a controllable numerical damping in the integration process is not new and found, among other things, realization in $\alpha$-HHT (cf. [6]) or in the WBC-$\alpha$-method (cf. [21]). The aim of the development of such methods is to maximize the attenuation of high-frequency components while preserving the important low-frequency components. Using the generalized $\alpha$-method, this ratio is optimal, i.e. for a given attenuation of the high-frequency components, the attenuation of the low-frequency components is minimized. A brief specification of this method for general problems in structural mechanics in given in Appendix 0.A.

Following [20], using a state vector

[TABLE]

describing the rod, we can rewrite its equations of motion (1a)–(1d) in terms of a system

[TABLE]

Here

[TABLE]

is the mass matrix and $\hat{\boldsymbol{\mathsf{K}}}=-\mathsf{adiag}(\boldsymbol{\mathsf{1}},\boldsymbol{\mathsf{1}},\boldsymbol{\mathsf{K}},\boldsymbol{\mathsf{1}})$ is the stiffness matrix with $\boldsymbol{\mathsf{K}}=\mathsf{diag}(EI_{1},EI_{2},G\mu)$, in which the bending stiffness into the direction of the principal components of the cross section $A$ is denoted by $EI_{1,2}$, and the torsional stiffness by $G\mu$. As above, the rod’s density is given by $\rho$, the Young’s modulus by $E$ and the shear modulus by $G$. The nonlinear terms are included in the nonlinearity $\boldsymbol{\Lambda}$. We do not explicitly write out the resulting equations here for brevity and refer to [20] for the explicit form of (12). In the generalized $\alpha$-method, the update schemes of positions and velocities at point $i$ in time correspond to those of the classical Newmark integrator. Its accurate and efficient application in the context of the simulation of elastic rods was demonstrated in [20].

According to [20], we consider four different test cases: (i) a sinus-like shaped rod which is released under gravity from a horizontal position (no damping); (ii) a highly damped helical rod subject to a time-varying end point load (the damping is obtained by setting the integration parameters (see Appendix 0.A) of the generalized $\alpha$-method to $\alpha:=\alpha_{m}=\alpha_{f}=0.4$, $\beta=0$, and $\gamma=1.0$); (iii) a straight rod (45 cm) subject to a time-varying torque; (iv) a helical rod with low damping that is excited by a force parallel to the axis of the helix and released after $0.1\,\text{s}$ showing the typical oscillating behavior of a steel-like coil spring.

We simulate these scenarios using the generalized $\alpha$-method and the symbolic-numeric integration scheme presented in this contribution. All fibers are discretized using $100$ individual segments. Using the symbolic-numeric method, damping is incorporated using the linear Rayleigh damping as described in [15]. We enforce a maximally tolerated relative $L^{2}$-error of 1% in the position and velocity space in order to ensure sufficient accuracy and measure the required computation time on a machine with an Intel(R) Xeon E5 with 3.5 GHz and 32 GB DDR-RAM without parallelization. For all test cases, we obtain significant speedups of the presented symbolic-numeric method (“snm”) compared to the generalized $\alpha$-method (“$\alpha$”), in particular:

1. (i)

speedup of a factor more than 20 ($\alpha$: $4.1\,\text{s}$; snm: $0.2\,\text{s}$); 2. (ii)

speedup of over $21\times$ ($\alpha$: $4.3\,\text{s}$; snm: $0.2\,\text{s}$); 3. (iii)

speedup of approx. $19\times$ ($\alpha$: $3.8\,\text{s}$; snm: $0.2\,\text{s}$); 4. (iv)

speedup of approx. $34\times$ ($\alpha$: $6.8\,\text{s}$; snm: $0.2\,\text{s}$).

Please note that the computation time of $0.2\,\text{s}$ of the presented symbolic-numeric method is constant for all test cases (with identical duration of $8\,\text{s}$).

5 Conclusison

Based on the closed form solution to the kinematic part (1a)–(1b) of the governing Cosserat system (1a)–(1d) and with assistance of Maple, we developed a new symbolic-numeric method for their integration. Our computational experiments demonstrate the superiority of the new method over the generalized $\alpha$-method for the accurate and efficient integration of Cosserat rods. Its application prevents from numerical instabilities and allows for highly accurate and efficient simulations. This clearly shows the usefulness of the constructed analytical solution to the kinematic equations and should enable more complex and realistic Cosserat rod-based scenarios to be explored in scientific computing without compromising efficiency.

Acknowledgements

The authors appreciate the insightful comments of the anonymous referees. This work has been partially supported by the King Abdullah University of Science and Technology (KAUST baseline funding), the Max Planck Center for Visual Computing and Communication (MPC-VCC) funded by Stanford University and the Federal Ministry of Education and Research of the Federal Republic of Germany (BMBF grants FKZ-01IMC01 and FKZ-01IM10001), the Russian Foundation for Basic Research (grant 16-01-00080) and the Ministry of Education and Science of the Russian Federation (agreement 02.a03.21.0008).

Appendix 0.A Generalized $\alpha$-Method

In this appendix, we briefly explain the application of the generalized $\alpha$-method for the common case of a system described by the standard equation from structural mechanics,

[TABLE]

in which $\boldsymbol{M}$,$\boldsymbol{D}$, and $\boldsymbol{K}$ denote the mass, damping, and stiffness matrices. The time-dependent displacement vector is given by $\boldsymbol{x}(t)$ and its first- and second-order temporal derivatives describe velocity and acceleration. The vector $\boldsymbol{\Lambda}(t)$ describes external forces acting on the system at time $t$. We are searching for functions $\boldsymbol{x}(t)$, $\boldsymbol{\upsilon}(t)=\dot{\boldsymbol{x}}(t)$, and $\boldsymbol{a}(t)=\ddot{\boldsymbol{x}}(t)$ satisfying (13) for all $t$ with initial conditions $\boldsymbol{x}(t_{0})=\boldsymbol{x}_{0}$ and $\boldsymbol{\upsilon}(t_{0})=\boldsymbol{\upsilon}_{0}$.

For the employment of the generalized $\alpha$-method, we can write the integration scheme with respect to (13) as follows:

[TABLE]

with the substitution rule $(\cdot)_{1-\alpha}:=(1-\alpha)(\cdot)_{i}+\alpha(\cdot)_{i-1}$ and the approximations

[TABLE]

The parameters $\alpha_{m}$, $\alpha_{f}$, $\gamma$, and $\beta$ are integration coefficients.