A stable finite element method for low inertia undulatory locomotion in three dimensions

TL;DR

This paper introduces a novel stable finite element method for simulating the low-inertia, three-dimensional undulatory locomotion of microswimmers, accurately capturing elastic and viscous deformations with preserved geometric properties.

Contribution

The authors develop a new finite element scheme that ensures stability and geometric fidelity for modeling low-inertia microswimmer dynamics in three dimensions.

Findings

The method is stable and preserves frame orthonormality up to machine precision.

Numerical experiments validate the method's accuracy and applicability.

The scheme effectively simulates undulatory locomotion in low-inertia regimes.

Abstract

We present and analyse a numerical method for understanding the low-inertia dynamics of an open, inextensible viscoelastic rod - a long and thin three dimensional object - representing the body of a long, thin microswimmer. Our model allows for both elastic and viscous, bending and twisting deformations and describes the evolution of the midline curve of the rod as well as an orthonormal frame which fully determines the rod's three dimensional geometry. The numerical method is based on using a combination of piecewise linear and piecewise constant functions based on a novel rewriting of the model equations. We derive a stability estimate for the semi-discrete scheme and show that at the fully discrete level that we have good control over the length element and preserve the frame orthonormality conditions up to machine precision. Numerical experiments demonstrate both the good properties of the method as well as the applicability of the method for simulating undulatory locomotion in the low-inertia regime.

Figures

Tables

Table 1. Table 1 . Summary of discretization choices for terms in model. Recall the V h subscript 𝑉 ℎ V_{h} is the space of piecewise linear functions and Q h subscript 𝑄 ℎ Q_{h} is the space of piecewise constant functions.

Variable	label	discrete space
Position	${𝒙}_h$	$V_h^3$
Lagrange multiplier	$p_h$	$Q_h$
Vector curvature	${𝜿}_h$	$V_{h,0}^3+{\kappa }_b$
Twist	${\gamma }_h$	$Q_h$
Tangential angular velocity	$m_h$	$V_h$
Normal moment	${𝒚}_h$	$V_{h,0}^3$
Tangential moment	$z_h$	$Q_h$
Frame vectors ($j=1,2$)	${𝒆}_h^j$	$V_{h,0}^3$

Table 2. Table 2 . Maximum error of local length constraint ℱ 1 ​ ( t n ) subscript ℱ 1 superscript 𝑡 𝑛 \mathcal{F}_{1}(t^{n}) for varying discretization parameters for relaxation test.

$\Delta t$	$N$	${\max }_n{ℱ}_1\left(t^n\right)$	$eoc\left(\Delta t\right)$
$1.00000$	16	$3.46788\cdot {10}^{-2}$	–
$2.50000\cdot {10}^{-1}$	32	$5.64486\cdot {10}^{-3}$	$1.30952$
$6.25000\cdot {10}^{-2}$	64	$4.89655\cdot {10}^{-4}$	$1.76355$
$1.56250\cdot {10}^{-2}$	128	$3.34948\cdot {10}^{-5}$	$1.93488$
$3.90625\cdot {10}^{-3}$	256	$2.14247\cdot {10}^{-6}$	$1.98330$
$9.76562\cdot {10}^{-4}$	512	$1.34687\cdot {10}^{-7}$	$1.99580$

Table 3. Table 3 . Maximum error of frame orthogonality constraint ℱ 2 ​ ( t n ) subscript ℱ 2 superscript 𝑡 𝑛 \mathcal{F}_{2}(t^{n}) for varying discretization parameters

$\Delta t$	$N$	${\max }_n{ℱ}_2\left(t^n\right)$	${\max }_n\left({ℱ}_2\left(t^n\right)-{ℱ}_2\left(t^{n-1}\right)\right)$
$1.00000$	16	$1.51801\cdot {10}^{-15}$	$2.46718\cdot {10}^{-16}$
$2.50000\cdot {10}^{-1}$	32	$5.09235\cdot {10}^{-15}$	$2.31888\cdot {10}^{-16}$
$6.25000\cdot {10}^{-2}$	64	$1.20420\cdot {10}^{-14}$	$2.11755\cdot {10}^{-16}$
$1.56250\cdot {10}^{-2}$	128	$3.95711\cdot {10}^{-14}$	$2.26700\cdot {10}^{-16}$
$3.90625\cdot {10}^{-3}$	256	$1.72853\cdot {10}^{-13}$	$2.35206\cdot {10}^{-16}$
$9.76562\cdot {10}^{-4}$	512	$8.70498\cdot {10}^{-13}$	$2.02678\cdot {10}^{-16}$

Table 4. Table 4 . Results for two dimensional test case.

$\Delta t$	$N$	${\max }_n{ℱ}_1\left(t^n\right)$	$eoc\left(\Delta t\right)$
$1.00000$	16	$6.93355\cdot {10}^{-1}$	–
$2.50000\cdot {10}^{-1}$	32	$4.71797\cdot {10}^{-2}$	$1.93868$
$6.25000\cdot {10}^{-2}$	64	$2.47289\cdot {10}^{-3}$	$2.12695$
$1.56250\cdot {10}^{-2}$	128	$1.58417\cdot {10}^{-4}$	$1.98220$
$3.90625\cdot {10}^{-3}$	256	$9.95399\cdot {10}^{-6}$	$1.99616$
$9.76562\cdot {10}^{-4}$	512	$6.23034\cdot {10}^{-7}$	$1.99895$

Table 5. Table 5 . Results for three dimensional test case.

$\Delta t$	$N$	${\max }_n{ℱ}_1\left(t^n\right)$	$eoc\left(\Delta t\right)$
$1.00000$	16	$2.43271$	–
$2.50000\cdot {10}^{-1}$	32	$9.14325\cdot {10}^{-2}$	$2.36686$
$6.25000\cdot {10}^{-2}$	64	$3.05806\cdot {10}^{-3}$	$2.45101$
$1.56250\cdot {10}^{-2}$	128	$1.98811\cdot {10}^{-4}$	$1.97157$
$3.90625\cdot {10}^{-3}$	256	$1.29005\cdot {10}^{-5}$	$1.97295$
$9.76562\cdot {10}^{-4}$	512	$8.11140\cdot {10}^{-7}$	$1.99566$

Table 6. Table 6 . Comparison of final centre of mass

$\Delta t$	$N$	Difference	Difference per time step
$1.00000$	$16$	$1.43005\cdot {10}^{-14}$	$5.72018\cdot {10}^{-16}$
$2.50000\cdot {10}^{-1}$	$32$	$1.57213\cdot {10}^{-14}$	$1.57213\cdot {10}^{-16}$
$6.25000\cdot {10}^{-2}$	$64$	$3.45280\cdot {10}^{-14}$	$8.63199\cdot {10}^{-17}$
$1.56250\cdot {10}^{-2}$	$128$	$6.04952\cdot {10}^{-13}$	$3.78095\cdot {10}^{-16}$
$3.90625\cdot {10}^{-3}$	$256$	$1.87147\cdot {10}^{-12}$	$2.92417\cdot {10}^{-16}$
$9.76562\cdot {10}^{-4}$	$512$	$4.87944\cdot {10}^{-13}$	$1.90603\cdot {10}^{-17}$

Table 7. Table 7 . Comparison of timings

$\Delta t$	$N$	Time (2d)	Time (3d)	Ratio
$1.00000$	$16$	$9.09700\cdot {10}^{-1}$ $\pm$ $1.87\%$	$9.25140\cdot {10}^{-1}$ $\pm$ $8.90\%$	$1.01697$
$2.50000\cdot {10}^{-1}$	$32$	$9.43000\cdot {10}^{-1}$ $\pm$ $2.94\%$	$9.89800\cdot {10}^{-1}$ $\pm$ $3.80\%$	$1.04963$
$6.25000\cdot {10}^{-2}$	$64$	$1.30110$ $\pm$ $3.09\%$	$1.74390$ $\pm$ $3.24\%$	$1.34033$
$1.56250\cdot {10}^{-2}$	$128$	$3.67120$ $\pm$ $5.04\%$	$8.17420$ $\pm$ $1.07\%$	$2.22657$
$3.90625\cdot {10}^{-3}$	$256$	$2.21780\cdot {10}^1$ $\pm$ $1.20\%$	$5.15220\cdot {10}^1$ $\pm$ $3.20\%$	$2.32311$
$9.76562\cdot {10}^{-4}$	$512$	$1.63716\cdot {10}^2$ $\pm$ $6.29\%$	$3.97100\cdot {10}^2$ $\pm$ $3.10\%$	$2.42554$

Equations

$$\frac{1}{\left|{\mathbf{x}}_u\right|}{\mathbf{e}}_u^0=\alpha {\mathbf{e}}^1+\beta {\mathbf{e}}^2,\frac{1}{\left|{\mathbf{x}}_u\right|}{\mathbf{e}}_u^1=-\alpha {\mathbf{e}}^0+\gamma {\mathbf{e}}^2,\frac{1}{\left|{\mathbf{x}}_u\right|}{\mathbf{e}}_u^2=-\beta {\mathbf{e}}^0-\gamma {\mathbf{e}}^1,$$

$$\frac{1}{\left|{\mathbf{x}}_u\right|}{\mathbf{\tau }}_u=\mathbf{\Omega }\times \mathbf{\tau },\frac{1}{\left|{\mathbf{x}}_u\right|}{\mathbf{e}}_u^1=\mathbf{\Omega }\times {\mathbf{e}}^1,\frac{1}{\left|{\mathbf{x}}_u\right|}{\mathbf{e}}_u^2=\mathbf{\Omega }\times {\mathbf{e}}^2.$$

$${\mathbf{\tau }}_t=\mathbf{\omega }\times \mathbf{\tau },{\mathbf{e}}_t^1=\mathbf{\omega }\times {\mathbf{e}}^1,{\mathbf{e}}_t^2=\mathbf{\omega }\times {\mathbf{e}}^2.$$

$$\frac{1}{\left|{\mathbf{x}}_u\right|}{\mathbf{\omega }}_u-{\mathbf{\Omega }}_t=\mathbf{\omega }\times \Omega .$$

$$\frac{m_u}{\left|{\mathbf{x}}_u\right|}-{\gamma }_t=\left(\frac{{\mathbf{\omega }}_u}{\left|{\mathbf{x}}_u\right|}-{\mathbf{\Omega }}_t\right)\cdot \mathbf{\tau }+\mathbf{\omega }\cdot \frac{{\mathbf{\tau }}_u}{\left|{\mathbf{x}}_u\right|}+\mathbf{\Omega }\cdot {\mathbf{\tau }}_t.$$

$$\frac{m_u}{\left|{\mathbf{x}}_u\right|}-{\gamma }_t=\mathbf{\tau }\cdot (\mathbf{\omega }\times \mathbf{\Omega })+\mathbf{\omega }\cdot \mathbf{\kappa }+\mathbf{\tau }\cdot (\mathbf{\Omega }\times \mathbf{\omega })=\mathbf{\omega }\cdot \mathbf{\kappa }.$$

$$\frac{m_u}{\left|{\mathbf{x}}_u\right|}={\gamma }_t-\frac{{\mathbf{x}}_{tu}}{\left|{\mathbf{x}}_u\right|}\cdot \mathbf{\tau }\times \mathbf{\kappa }.$$

$$\mathbf{K}+\frac{1}{\left|{\mathbf{x}}_u\right|}{\mathbf{F}}_u=0,{\mathbf{K}}^{\mathrm{rot}}+\mathbf{\tau }\times \mathbf{F}+\frac{1}{\left|{\mathbf{x}}_u\right|}{\mathbf{M}}_u=0,$$

$$\boldsymbol{M}=\boldsymbol{\tau}\times\Bigl{\{}A\bigl{(}(\alpha-\alpha^{0})\boldsymbol{e}^{1}+(\beta-\beta^{0})\boldsymbol{e}^{2}\bigr{)}+B(\alpha_{t}\boldsymbol{e}^{1}+\beta_{t}\boldsymbol{e}^{2})\Bigr{\}}+C(\gamma-\gamma^{0})\boldsymbol{\tau}+D\gamma_{t}\boldsymbol{\tau}.$$

$$\mathbf{y}$$

$$z$$

$$\mathbf{M}=\mathbf{\tau }\times \mathbf{y}+z\mathbf{\tau }.$$

$${\left|{\mathbf{x}}_u\right|}_t=\mathbf{\tau }\cdot {\mathbf{x}}_{tu}=0.$$

$$\mathbf{K}=\mathbb{K}{\mathbf{x}}_t,{\mathbf{K}}^{\mathrm{rot}}=-K^{\mathrm{rot}}m\mathbf{\tau },$$

$$\frac{1}{\left|{\mathbf{x}}_u\right|}{\mathbf{M}}_u=\mathbf{\kappa }\times \mathbf{y}+\mathbf{\tau }\times \frac{1}{\left|{\mathbf{x}}_u\right|}{\mathbf{y}}_u+\frac{1}{\left|{\mathbf{x}}_u\right|}{\mathbf{z}}_u\mathbf{\tau }+z\mathbf{\kappa }.$$

$$\mathbf{\tau }\cdot \frac{1}{\left|{\mathbf{x}}_u\right|}{\mathbf{M}}_u$$

$$\mathbf{\tau }\times \frac{1}{\left|{\mathbf{x}}_u\right|}{\mathbf{M}}_u$$

$$\mathbb{K}{\mathbf{x}}_t+\frac{1}{\left|{\mathbf{x}}_u\right|}(p\mathbf{\tau }{)}_u+\frac{1}{\left|{\mathbf{x}}_u\right|}{\mathbf{f}}_u=0.$$

$$\frac{1}{\left|{\mathbf{x}}_u\right|}{\mathbf{f}}_u+\frac{1}{\left|{\mathbf{x}}_u\right|}{\left(\frac{1}{\left|{\mathbf{x}}_u\right|}\mathbf{\tau }\times {\mathbf{M}}_u\right)}_u$$

$$-K^{\mathrm{rot}}m+\frac{1}{\left|{\mathbf{x}}_u\right|}{\mathbf{M}}_u\cdot \mathbf{\tau }$$

$$\displaystyle\mathbb{K}\boldsymbol{x}_{t}+\frac{1}{|\boldsymbol{x}_{u}|}\bigl{(}p\boldsymbol{\tau}\bigr{)}_{u}+\frac{1}{|\boldsymbol{x}_{u}|}\bigl{(}(\mathbb{I}-\boldsymbol{\tau}\otimes\boldsymbol{\tau})\frac{\boldsymbol{y}_{u}}{|\boldsymbol{x}_{u}|}+z\boldsymbol{\tau}\times\boldsymbol{\kappa}\bigr{)}_{u}$$

$$-K^{\mathrm{rot}}m+\frac{z_u}{|{\mathbf{x}}_u|}+\mathbf{y}\cdot (\mathbf{\tau }\times \mathbf{\kappa })$$

$$\mathbf{\tau }\cdot {\mathbf{x}}_{tu}$$

$$p\mathbf{\tau }+(\mathbb{I}-\mathbf{\tau }\otimes \mathbf{\tau })\frac{{\mathbf{y}}_u}{|{\mathbf{x}}_u|}+z\mathbf{\tau }\times \mathbf{\kappa }$$

$$\mathbf{\tau }\times \mathbf{y}+z\mathbf{\tau }$$

$$\mathbf{\kappa }=\frac{1}{\left|{\mathbf{x}}_u\right|}{\mathbf{\tau }}_u=\frac{1}{\left|{\mathbf{x}}_u\right|}{\left(\frac{{\mathbf{x}}_u}{\left|{\mathbf{x}}_u\right|}\right)}_u.$$

$$\mathbf{\kappa }={\alpha }^0{\mathbf{e}}^1+{\beta }^0{\mathbf{e}}^2\text{mbox}atu=0,1.$$

$${\alpha }_t{\mathbf{e}}^1+{\beta }_t{\mathbf{e}}^2=(\mathbb{I}-\mathbf{\tau }\otimes \mathbf{\tau }){\mathbf{\kappa }}_t-m\mathbf{\tau }\times \mathbf{\kappa }.$$

$$\displaystyle\int_{0}^{1}\mathbb{K}\boldsymbol{x}_{t}\cdot\boldsymbol{\phi}|\boldsymbol{x}_{u}|-\int_{0}^{1}p\boldsymbol{\tau}\cdot\boldsymbol{\phi}_{u}-\int_{0}^{1}\bigl{(}(\mathbb{I}-\boldsymbol{\tau}\otimes\boldsymbol{\tau})\frac{1}{|\boldsymbol{x}_{u}|}\boldsymbol{y}_{u}+z\boldsymbol{\tau}\times\boldsymbol{\kappa}\bigr{)}\cdot\boldsymbol{\phi}_{u}$$

$$\displaystyle\int_{0}^{1}\bigl{(}\boldsymbol{y}-A(\boldsymbol{\kappa}-\alpha^{0}\boldsymbol{e}_{1}-\beta^{0}\boldsymbol{e}_{2})-B\bigl{(}(\mathbb{I}-{\boldsymbol{\tau}}\otimes{\boldsymbol{\tau}})\boldsymbol{\kappa}_{t}-m{\boldsymbol{\tau}}\times\boldsymbol{\kappa}\bigr{)}\bigr{)}\left|\boldsymbol{x}_{u}\right|$$

$${\int }_0^1\mathbf{\kappa }\cdot \mathbf{\psi }|{\mathbf{x}}_u|+\frac{{\mathbf{x}}_u}{|{\mathbf{x}}_u|}\cdot {\mathbf{\psi }}_u$$

$${\int }_0^1-K^{\mathrm{rot}}mv|{\mathbf{x}}_u|+{\int }_0^1\mathbf{y}\cdot (\mathbf{\tau }\times \mathbf{\kappa })v|{\mathbf{x}}_u|-{\int }_0^{1}z{v}_u$$

$${\int }_0^1(z-C(\gamma -{\gamma }^0)-D{\gamma }_t)q|{\mathbf{x}}_u|$$

$${\int }_0^1{\gamma }_{t}q|{\mathbf{x}}_u|-{\int }_0^1{m}_{u}q+{\int }_0^1\mathbf{\tau }\times \mathbf{\kappa }\cdot {\mathbf{x}}_{tu}q$$

$${\int }_0^{1}q\mathbf{\tau }\cdot {\mathbf{x}}_{tu}$$

$$\displaystyle\int_{0}^{1}\bigl{(}\boldsymbol{e}^{j}_{t}-\bigl{(}\boldsymbol{\tau}\times\boldsymbol{\tau}_{t}+m\boldsymbol{\tau}\bigr{)}\times\boldsymbol{e}^{j}\bigr{)}\cdot\boldsymbol{\phi}\left|\boldsymbol{x}_{u}\right|=0,\quad\mbox{ for }j=1,2,$$

$$\mathbf{x}(\cdot ,0)={\mathbf{x}}^0,{\mathbf{e}}^1(\cdot ,0)={\mathbf{e}}^{1,0},{\mathbf{e}}^2(\cdot ,0)={\mathbf{e}}^{2,0},$$

Full text

Funding. This work was supported by a Leverhulme Trust Early Career Fellowship.

Keywords. Kirchhoff rod; viscoelastic materials; finite element methods; biomechanics; undulatory locomotion

Subject class. 74S05, 74K10, 74D05, 65M60, 74L15, 92C10

1. Introduction

1.1. Background

The dynamics of active microswimmers has long captured the interest of physicists, mathematicians and engineers, not to mention researchers in various biological fields. Undulation, passing bending waves down the body, is a common strategy used across many orders of magnitude [40] ranging from bacteria [61, 57] and spermatozoa [37] to larger fish and whales [38]. It is especially common in the low inertia regime, where viscous forces of the surrounding media dominate over inertial forces and the scallop theorem prohibits self-propelled locomotion for time-reversal symmetric sequences of body postures [68, 79]. For additional reading, we refer the reader to a large body of excellent reviews on animal locomotion [44, 58, 19], fluid dynamics at low Reynolds number [18, 60, 61, 16, 68, 83, 37], and the biophysics and biology of cell motility [8, 47, 48, 14, 7, 15]. There are many computational studies in this area which have focussed on solving a fully detailed three dimensional formulation capturing many aspects of the problem including the dynamics and its interaction with the body (e.g. [76, 66, 74, 34, 65, 73]). These large scale computational studies provide many interesting results giving insights on both complex fluids and animal locomotion (see, e.g., the review [58] and references therein). An alternative is to choose a reduced, lower dimensional (in our case one dimensional) description of the body and represent all forces on this reduced description. This allows orders of magnitude faster simulations by considering only a one dimensional problem which significantly reduces the numbers of computational degrees of freedom.

In this work we develop, analyse and show experiments using a novel computational approach to the evolution of an open viscoelastic rod, representing the body of a long, thin microswimmer - a three dimensional object with one axis much longer than the other two. We assume the body is embedded in three dimensions which can undergo bending and twisting deformations. The resulting problem can be described as a system of one dimensional partial differential equations for a midline curve and an orthonormal frame which describes the conformation of cross sections to the midline. The model is a natural generalisation of the locomotion model [46] to three spatial dimensions. The elastic terms we consider arise from the classical Kirchhoff-Love model for an elastic rod [52] which are combined in a simple linear Voigt viscoelastic model [1, 53]. Our model avoids considering shearing or extensional deformations from the full Cosserat or Timoshenko models. In general this nonlinear system can not be solved analytically and requires the use of computer methods. The reduction to a one dimensional object allows for significant reduction in the complexity of the resulting mathematical model in particular with respect to the computational effort required to solve the problem.

Our new method tackles three key challenges. First when considering active locomotion one must be able to map directly between anatomical detail of the organism under consideration and the geometry defined by the model. For example, we should be able to clearly identify where the muscles are located and in which direction do they apply force. Our approach to capture this results in equations for the midline of the rod coupled to equations for the orientation of the cross section of the rod (see the discussion in [56] and Sec. 4.2 for more details). These equations must be carefully coupled to ensure that we have an accurate and robust representation of the geometry which allows us to apply the biological forces appropriately. Second we have a moving geometry which is a priori unknown. Our scheme captures this with a parameterization which is equivalent to a moving mesh. As is typical in this type of problem we must ensure that the moving mesh does not become too distorted (see, for example, [33, 4] and references therein). Finally, in many biological systems various parameters and problem data may be unknown or poorly characterized. Any computational method should be both cheap enough to run so that parameterization studies may be run and also robust to input parameters so that a wide variety of behaviours can be observed. See [75, 49] for example. We will demonstrate through analysis of our method and numerical experiments that we can address all three challenges.

Similar models to those considered in this paper arise in many areas of natural science and engineering. For example, similar models have been considered for elastic ribbons and filaments [10, 9, 39], tangled hair and fibres [12, 81, 24, 30], plants [43, 42], and woven cloth [50]. A historical overview of the model used in this paper is given in [29].

We are interested in the locomotion of low inertia microswimmers. We explore the applicability of our method through the case study of the 1mm long nematode Caenorhabditis elegans. This animal has become a model organism studied by physicists, computer scientists and engineers partly due to its simple, undulating, periodic gait but also its experimental tractability and simple neuroanatomy [21]. In its natural environment, C. elegans grows mostly in rotten vegetation: a highly complex, heterogeneous three dimensional environment. It is cultured on the surface of agar gels for extensive use in laboratories. On agar, C. elegans move by propagating bending waves from head to tail to generate forward thrust [45, 23, 67, 11, 36, 59] where the bending arises from body wall muscles which line sides of the body. A two dimensional version of our viscoelastic model [46], developed originally for capturing worms or snakes moving across land, has previously been applied to C. elegans locomotion both on surfaces and in a range of both Newtonian and non-Newtonian media [77, 78, 36, 20, 28]. However, previous formulations have either linearized the equations, which means although postures are recovered trajectories are not, or viscosity is neglected, which limits the applicability of the model in less resistive environments. Recent experiments [13, 72] have shown that C. elegans achieves a different gait in three dimensions. Here we demonstrate that our computational tool is capable of accurately and robustly capturing such behaviours. The simulation results are meant to be indicative of the capabilities of the method rather than demonstrating any properties of the underlying model which is left to future work. The numerical method used in this section builds on the unpublished work [20] (see also [28]).

1.2. Model

In our model, the rod is described through a smooth, time-dependent parameterization of the midline $\boldsymbol{x}\colon[0,1]\times[0,T]\to\mathbb{R}^{3}$ and an orthonormal frame $\boldsymbol{e}^{0},\boldsymbol{e}^{1},\boldsymbol{e}^{2}\colon[0,1]\times[0,T]\to\mathbb{R}^{3}$ up to some final time $T>0$. We call $u\in[0,1]$ the material parameter. Since we do not allow shear deformations, we set $\boldsymbol{e}^{0}\equiv\boldsymbol{\tau}:=\boldsymbol{x}_{u}/\left|\boldsymbol{x}_{u}\right|$ the unit tangent to the midline. The other two coordinates $\boldsymbol{e}^{1}$ and $\boldsymbol{e}^{2}$ form an orthonormal basis of the cross section of the rod to the midline. See Fig. 1. We can form the skew system:

[TABLE]

where $\alpha$ and $\beta$ are smooth fields denoting the curvature of the midline in directions $\boldsymbol{e}^{1}$ and $\boldsymbol{e}^{2}$ and $\gamma$ is a smooth field denoting the twist of the orthonormal frame about the midline. At the core of the model is a moment which is the sum of elastic and viscous contributions. The elastic terms are proportional to differences between the fields $\alpha,\beta,\gamma$ and some desired values $\alpha^{0},\beta^{0},\gamma^{0}$. The viscous contribution is proportional to the time derivatives $\alpha_{t},\beta_{t},\gamma_{t}$.

Inertial terms are ignored and external forces are simply modelled as linear drag terms [51]. This model can be seen as a quasi-static or low inertia approximation of the full rod dynamics. Although this very simple approach to modelling external forces for undulatory locomotion has been used previously in undulatory locomotion models (see the references above), a more accurate model would include non-local terms [51, 60] or capture the external fluid using dynamic equations (e.g. [22, 70, 64, 62]). Any model for the surrounding fluid must be coupled to a model of the body. The focus in this work is to accurately capture the body mechanics combined with a simple model of surrounding fluid. Study of a more detailed model of the environment coupled to the model of body mechanics presented here is left to future work.

Our model can be seen as a simplification of the model presented in [55, 54] or a three dimensional version of [46]. Full details of the model are given in Sec. 2.

1.3. Computational method

Our approach is to discretize an appropriate formulation of the continuous equations directly. This allows us to use the structure of the equations to recover a robust numerical scheme through careful discretization choices. The discretization extends the approach of [63] to the case of non-constant twist and open curves.

Key to our numerical method is a well suited formulation of the continuous partial differential equation system. We start from the balance laws for linear and angular momentum and the linear viscoelastic constitutive law. This is combined with geometric constraints so that the solution variables are the position of the midline, line tension (a Lagrange multiplier for enforcing the length constraint), the curvature of the midline, the twist and angular velocity of the frame and two auxiliary variables which describe the bending and twisting moments. The continuous system of equations is discretized in space by a mixed finite element method where we use a mix of piecewise linear and piecewise constant approximation spaces. We use a Lagrange multiplier to enforce inextensibility (i.e. that the length of the curve is locally fixed) in a similar approach to [5]. Finally, we discretize in time using a semi-implicit method which results in a linear system of equations to solve at each time step, inspired by the approach in [31, 32, 63] (see also [4]), together with an update formula for the frame. The use of lower order finite element approximations, along with mass lumping [80], allows us to derive identities for various geometric quantities at the mesh vertices. The frame is updated using a Rodrigues formula where the rotation is specified by the frame’s angular moment which can be derived from the solution variables. We use the angular momentum and Rodrigues formula as an alternative to using Euler angles [71] or quaternions [55, 54]. Although direct use of the Rodrigues formula is not recommended for numerical applications in general [6], we will show that our formulation avoids problems for any angles. More details of the numerical scheme are given in Sec. 3.

We will demonstrate three key properties of our scheme:

• •

a semi-discrete stability result (coupled with computational evidence of fully discrete stability) which shows that we recover a discrete Lyapunov functional for our scheme;

• •

control over the length element which ensures that vertices in the moving mesh do not collide;

• •

preservation of the frame orthogonality conditions to ensure that the frame really does remain orthonormal over time.

We see these three properties as allowing our method to provide a computational tool both for the understanding of viscoelastic rods and also for domain experts working on undulatory locomotion in the low inertia regime. Further work is required to link this model for undulatory locomotion with a more detailed model of the surrounding fluid which would give a more accurate model and allow a greater variety of behaviours to be captured. Where such links are to be used a balance must be struck between the computational efficiency of a one dimensional approach against the greater detail provided by a fully three dimensional model.

In our previous unpublished work [20], we presented a similar scheme (also used in [28]). This paper builds on that work extending the scheme to three dimensions including twisting as well as bending contributions and generalising the material law to include viscous terms.

There are several very successful methods from discrete differential geometry which solve either the problem we consider or a generalisation. In this approach the rod is first discretized and then equations are derived to evolve those discrete quantities. We mention in particular the approach of [10] (extended to viscous threads in [9, 2]) which uses a discrete set of vertices (in our notations a piecewise linear curve) with a material frame on each edge between the vertices (piecewise constant in our notation) to represent a rod which was stored using a reduced curve-angle formulation [56]. The approach uses curvature and twist defined as integrated quantities based at vertices to define a Kirchhoff-Love energy and then applies discrete parallel transport and variation of holonomy to derive the update equations. This work was generalized to Cosserat rods by [41] who revert away from the reduced curvature-angle formulation storing the full frame. The super-helix and super-clothoid approaches [12, 17] use a piecewise constant or piecewise linear approximation of generalized curvature of a rod and then recover the geometry of the rod (position of midline and frame) using analytic expressions. The scheme results in a smooth curve with well defined curvatures and twist.

In the context of these schemes our model can be seen as using a set of vertices with a frame at each vertex to define the discrete geometry of the rod. The equations are straight discretization of the continuous scheme using a finite element approach and discrete equations to define curvature and relate twist with tangential angular velocity. Our choice to discretize continuum equations allows us to use standard finite element tools to both analyse and implement the method in a single unified framework. This approach is well suited to the parabolic nature of the equations we consider. It is our particular choice of both geometric discretization and direct discretization of the equations that allows us to demonstrate the properties of our scheme.

1.4. Outline

In Sec. 2, we present the continuous model we use and the discretization is shown in Sec. 3. Numerical experiments to demonstrate the efficacy of the method are shown in Sec. 4. The restriction of our scheme to a two dimensional problem is given in App. A.

2. Governing equations

2.1. Geometry

We consider a smooth, inextensible, unshearable rod embedded in $\mathbb{R}^{3}$ over a time interval $[0,T]$ for some $0<T<\infty$. The rod can be described by a (non-arc length) parameterization of the centre line $\boldsymbol{x}\colon[0,1]\times[0,T]\to\mathbb{R}^{3}$ and an oriented frame of reference $\boldsymbol{Q}\colon[0,1]\times[0,T]\to SO(3)$. For a discussion of rod representations see [56]. We will write derivatives with respect to the first coordinate, the material coordinate $u$, and the second coordinate, time $t$, with subscripts $(\cdot)_{u}$ and $(\cdot)_{t}$, respectively. Our assumptions imply $\left|\boldsymbol{x}_{u}\right|_{t}=0$ so that the length of the midline is fixed. We call the length of the midline curve $L$.

Rather than use the tensor $\boldsymbol{Q}$, we will use the equivalent orthonormal triad of unit vectors $\boldsymbol{Q}=\{\boldsymbol{e}^{0},\boldsymbol{e}^{1},\boldsymbol{e}^{2}\}$. We will use the convention for an unshearable rod that $\boldsymbol{e}^{0}\equiv\boldsymbol{\tau}$ the unit tangent to the centre line given by $\boldsymbol{\tau}=\boldsymbol{x}_{u}/|\boldsymbol{x}_{u}|$. See Fig. 1 for an example rod conformation.

We can recover the generalized curvature (sometimes called the Darboux vector) $\boldsymbol{\Omega}\colon[0,1]\times[0,T]\to\mathbb{R}^{3}$ which satisfies

[TABLE]

We decompose $\boldsymbol{\Omega}=\boldsymbol{\tau}\times\boldsymbol{\kappa}+\gamma\boldsymbol{\tau}$. We call $\boldsymbol{\kappa}=\boldsymbol{\tau}_{u}/\left|\boldsymbol{x}_{u}\right|$ the vector curvature with decomposition $\boldsymbol{\kappa}=\alpha\boldsymbol{e}_{1}+\beta\boldsymbol{e}_{2}$ for fields $\alpha=\boldsymbol{\kappa}\cdot\boldsymbol{e}_{1},\beta=\boldsymbol{\kappa}\cdot\boldsymbol{e}_{2}$ and call $\gamma$ the twist. Further, if the orthonormal triad transforms smoothly in time, then we recover the angular velocity of the frame $\boldsymbol{\omega}$ given by

[TABLE]

Again, we decompose $\boldsymbol{\omega}=\boldsymbol{\tau}\times\boldsymbol{\tau}_{t}+m\boldsymbol{\tau}$. The field $m$ denotes the angular velocity of the frame about the tangent vector field which must be computed separately to the centre line velocity.

There is a geometric relation between the strain vector $\boldsymbol{\Omega}$ and the angular velocity $\boldsymbol{\omega}$ [82]. Since the $u$ and $t$ derivatives commute, using the inextensibility of the parameterization (i.e. $\left|\boldsymbol{x}_{u}\right|_{t}=0$), we can compute that

[TABLE]

Next we subtract the time derivative of $\gamma=\boldsymbol{\Omega}\cdot\boldsymbol{\tau}$ from the scaled $u$ derivative (i.e. arc-length derivative) of $m=\boldsymbol{\omega}\cdot\boldsymbol{\tau}$ to see

[TABLE]

Applying 2.3 and 2.2, the definition of $\boldsymbol{\kappa}$ and rearranging the vector triple product gives

[TABLE]

Then applying the definition of $\boldsymbol{\omega}$ and rearranging the resulting vector triple product, we see that

[TABLE]

This final equation is important for providing a closing relation between the angular velocity and twist of the frame.

2.2. Model derivation

For an inextensible, unshearable rod we can write down the conservation of linear and angular momentum as [53]:

[TABLE]

where $\boldsymbol{K}$ is the external force, $\boldsymbol{F}$ is the internal force resultant, $\boldsymbol{K}^{\mathrm{rot}}$ is the external moment and $\boldsymbol{M}$ is the internal moment.

We assume that the rod is viscoelastic with preferred curvatures and preferred twist so that the internal moment is given by a linear Voigt model:

[TABLE]

Here $\alpha^{0},\beta^{0}$ and $\gamma^{0}$ are given fields which we allow to depend on the parameter $u$ and time $t$. We call $\alpha^{0}$ and $\beta^{0}$ preferred curvatures and $\gamma^{0}$ a preferred twist. The material parameters, which we allow to depend on the parameter $u$ but not time $t$, are $A$ the bending modulus, $C$ the twisting modulus, $B$ the bending viscosity and $D$ the twisting viscosity. Note that the material parameters will depend on the precise geometry of the cross section of the rod [54]. We assume that $A=A(u)\geq A_{0}>0$, $B=B(u)\geq 0$, $C=C(u)>C_{0}>0$ and $D=D(u)\geq 0$. We introduce the variables $\boldsymbol{y}$, the bending moment, and $z$ the twisting moment given by

[TABLE]

so that the moment is given by

[TABLE]

We assume that the tangential forces, $p\boldsymbol{\tau}$, act as a Lagrange multiplier to enforce the inextensibility constraint:

[TABLE]

We assume a linear drag response from the environment onto the rod by

[TABLE]

with strictly positive definite matrix $\mathbb{K}$ and strictly positive scalar coefficient $K^{\mathrm{rot}}$. Our model of drag is inspired by resistive force theory [51].

We combine the above model components in a way that is well suited to numerical computation. We start by calculating from 2.8 that

[TABLE]

Noting that $\boldsymbol{\tau}\cdot\boldsymbol{\kappa}=\boldsymbol{\tau}\cdot\boldsymbol{y}=\boldsymbol{\tau}\times(\boldsymbol{\kappa}\times\boldsymbol{y})=0$, we infer

[TABLE]

Here $\mathbb{I}$ is the $3\times 3$ identity matrix and $\otimes$ is the outer product given by $(\boldsymbol{a}\otimes\boldsymbol{b})_{ij}=\boldsymbol{a}_{i}\boldsymbol{b}_{j}$ for $i,j=1,2,3$, $\boldsymbol{a},\boldsymbol{b}\in\mathbb{R}^{3}$ which satisfies $\boldsymbol{\tau}\times(\boldsymbol{\tau}\times\boldsymbol{a})=(\mathbb{I}-\boldsymbol{\tau}\otimes\boldsymbol{\tau})\boldsymbol{a}$ for all $a\in\mathbb{R}^{3}$.

We decompose the internal force resultant into tangential and normal components by $\boldsymbol{F}=p\boldsymbol{\tau}+\boldsymbol{f}$. We call $p$ the pressure and $\boldsymbol{f}$ the normal force resultant. We then use the force decomposition and 2.10 in the linear momentum balance 2.5:

[TABLE]

From taking cross product of $\boldsymbol{\tau}$ with the angular momentum balance 2.5 and then the arc-length derivative (the scaled $u$ derivative), we infer:

[TABLE]

Next, we take the scalar product of $\boldsymbol{\tau}$ with the angular momentum balance 2.5 using 2.10 to see:

[TABLE]

Finally we apply the expressions for derivatives of $\boldsymbol{M}$ from 2.11 in 2.12, 2.13 and 2.14 and combine with 2.9 to give our model:

[TABLE]

Remark 1*.*

The system of partial differential equations 2.15 can be seen as a nonlinear fourth order parabolic equation for the parameterization $x$, subject to the nonlinear inextensibility constraint, coupled to a nonlinear second order parabolic equation for the twist $\gamma$.

2.3. Weak form

We will write down a weak form which we will use for our finite element method in the next section. When writing down the weak formulation, we combine equations for the conservation laws 2.15a and 2.15b, the constitutive laws 2.7, the geometric relation 2.4, and the inextensibility constraint 2.15c. When writing down the constitute equation we add the Laplace-Beltrami identity for $\boldsymbol{\kappa}$ the vector curvature:

[TABLE]

We must also impose boundary conditions for $\boldsymbol{\kappa}$, since the curvature is not well defined at $u=0,1$, which we set to be equal to the prescribed curvatures here:

[TABLE]

We also note that for the bending viscosity terms we can write

[TABLE]

We derive the weak form of the problem by multiplying by appropriate test functions and integrating over the centre line.

We let $Q=L^{2}(0,1)$ denote the space of square integrable functions on $(0,1)$, $V=H^{1}(0,1)$ the Sobolev space of functions in $L^{2}(0,1)$ with a weak derivative in $L^{2}(0,1)$ and $V_{0}$ the space of functions in $V$ with zero trace [35]. Unless otherwise stated, integrals are with respect to the measure $\mathrm{d}u$.

Problem 2**.**

Given preferred curvatures, $\alpha^{0},\beta^{0}$, a preferred twist, $\gamma^{0}$, and initial conditions for the parameterization $\boldsymbol{x}^{0}$ and frame $\boldsymbol{e}^{1,0},\boldsymbol{e}^{2,0}$, find $\boldsymbol{x},\boldsymbol{y},\boldsymbol{\kappa}\colon[0,1]\times[0,T)\to\mathbb{R}^{3}$ (with the conditions 2.15d, 2.15e and 2.17 at the boundaries), $m,z,\gamma\colon[0,1]\times[0,T)\to\mathbb{R}$, and $\boldsymbol{e}^{1},\boldsymbol{e}^{2}\colon[0,1]\times[0,T)\to\mathbb{R}^{3}$ such that, for almost every $t\in(0,T)$:

[TABLE]

3. Numerical method

3.1. Finite element spaces

We take a partition of $[0,1]$ by $N$ points $u_{1}=0<u_{2}<\ldots<u_{N}=1$. We call $\{u_{1},\ldots,u_{N}\}$ the mesh. We will use a combination of piecewise linear and piecewise constant functions. We introduce the spaces:

[TABLE]

We will denote by $V_{h,0}$ the space of finite element functions in $v_{h}\in V_{h}$ such that $v_{h}(0)=v_{h}(1)=0$. We will use the subscript $h$, defined to be the maximum mesh spacing, to denote discrete quantities. Temporal and spatial derivatives of discrete functions will be denoted by subscript $u$ and $t$ with a comma separating from the subscript $h$: e.g. $\eta_{h,u}$ denote the spatial derivative of $\eta_{h}$.

For a discrete parameterization $\boldsymbol{x}_{h}\in V_{h}^{3}$, we introduce two different tangent vector fields. First, $\boldsymbol{\tau}_{h}\in Q_{h}^{3}$ as the piecewise constant normalized derivative of $\boldsymbol{x}_{h}$:

[TABLE]

We will also require a piecewise linear approximation of $\boldsymbol{\tau}_{h}$ written $\tilde{\boldsymbol{\tau}}_{h}\in V_{h}^{3}$ with vertex values given by

[TABLE]

where $\boldsymbol{\tau}_{h}(u_{i}^{\pm},\cdot)$ is $\boldsymbol{\tau}_{h}$ evaluated on the left (or right) element to the vertex $u_{i}$.

We will apply mass lumping [80] using the notations:

[TABLE]

where $I_{h}$ is the Lagrangian interpolation operator $C([0,1])\to V_{h}$.

Finally we denote by $V_{h,0}^{3}+\boldsymbol{\kappa}_{b}(\cdot,t)$ the space of finite element functions $\boldsymbol{v}_{h}\in V_{h}^{3}$ which match the boundary conditions for $\boldsymbol{\kappa}_{h}$:

[TABLE]

where $\boldsymbol{e}^{1}_{h},\boldsymbol{e}^{2}_{h}\in V_{h}^{3}$ will denote components of the orthonormal frame that we will solve for as part of the method. The space $V_{h,0}^{3}+\boldsymbol{\kappa}_{b}$ will in general be time dependent.

3.2. Semi-discrete problem

We directly discretize the weak form 2.18. The choice of piecewise linear or piecewise constant approximation spaces for the different functions is determined by the properties we will show in Lem. 5.

At this stage of discretization the choices are between which discrete function spaces each solution variable should live in and how to implement boundary conditions. We choose piecewise linear approximations of position $\boldsymbol{x}$, bending moment $\boldsymbol{y}$, curvature $\boldsymbol{\kappa}$, angular momentum $m$ and frame $\boldsymbol{e}_{1}$ and $\boldsymbol{e}_{2}$ and piecewise constant approximations of twisting moment $z$, twist $\gamma$ and the Lagrange multiplier $p$. We choose to enforce boundary conditions for the bending moment $\boldsymbol{y}$ and curvature $\boldsymbol{\kappa}$ in the function spaces but boundary conditions for the twisting moment $z$ and twist $\gamma$ arise as natural boundary conditions. Imposing natural conditions means that the boundary conditions are only achieved exactly in the limit of small mesh spacing (numerical confirmation not shown). We will see our choices naturally lead to the key properties of our scheme. A summary of discretization choices in given in Tab. 1.

Problem 3**.**

Given preferred curvatures $\alpha^{0},\beta^{0}$, a preferred twist $\gamma^{0}$, and initial conditions for the parameterization $\boldsymbol{x}_{h}^{0}$ and frame $\boldsymbol{e}_{h}^{1,0},\boldsymbol{e}_{h}^{2,0}$, for $t\in[0,T)$, find $\boldsymbol{x}_{h}(\cdot,t)\in V_{h}^{3},\boldsymbol{y}_{h}(\cdot,t)\in V_{h,0}^{3},\boldsymbol{\kappa}_{h}(\cdot,t)\in V_{h,0}^{3}+\boldsymbol{\kappa}_{b}(\cdot,t)$, $m_{h}(\cdot,t)\in V_{h},z_{h}(\cdot,t),\gamma_{h}(\cdot,t),p_{h}(\cdot,t)\in Q_{h}$, $\boldsymbol{e}^{1}_{h}(\cdot,t),\boldsymbol{e}^{2}_{h}(\cdot,t)\in V_{h}^{3}$ such that, for all $t\in(0,T)$:

[TABLE]

Remark 4*.*

Decoupling variables for curvature $\boldsymbol{w}_{h}$ and position $\boldsymbol{x}_{h}$ should be interpreted as a tool for solving the partial differential equation system. This is widely used when solving geometric partial differential equations (see e.g. [32] with the convergence results in [26] and the review [27]). Since we compute with piecewise linear curves we cannot formulate an exact curvature of the discrete curve. This is a key difference between the approach presented here and super-helix and super-clothoid approaches [12, 17].

Using $\boldsymbol{e}^{0}_{h}\equiv\tilde{\boldsymbol{\tau}}_{h}$, we will see $\{\boldsymbol{e}^{0}_{h},\boldsymbol{e}^{1}_{h},\boldsymbol{e}^{2}_{h}\}$ is a vertex-wise orthonormal frame. Indeed, we note that 3.3h implies we recover the vertex-wise relations

[TABLE]

where $\boldsymbol{\omega}_{h}=\tilde{\boldsymbol{\tau}}_{h}\times\tilde{\boldsymbol{\tau}}_{h,t}+m_{h}\tilde{\boldsymbol{\tau}}_{h}\in V_{h}^{3}$. This implies the following vertex-wise identities hold:

[TABLE]

so long as the initial values satisfy corresponding versions of these identities. In other words we have the $(\tilde{\boldsymbol{\tau}}_{h},\boldsymbol{e}^{1}_{h},\boldsymbol{e}^{2}_{h})$ form an orthonormal frame at each vertex.

Next, we note that we can write 3.3c as:

[TABLE]

So that, using the fact that $\left|\boldsymbol{\tau}_{h}\right|=1$, we can infer that $\boldsymbol{\kappa}_{h}$ and $\tilde{\boldsymbol{\tau}}_{h}$ are orthogonal at the vertices. Indeed, for all $i=2,\ldots,N-1$, we have

[TABLE]

Further, at $i=1$ and $i=N$, the boundary conditions give us that $\boldsymbol{\kappa}_{h}$ and $\tilde{\boldsymbol{\tau}}_{h}$ are orthogonal directly. This implies we can create a decomposition of $\boldsymbol{\kappa}_{h}$, at the vertices, into fields $\alpha_{h},\beta_{h}\in V_{h}$ given by

[TABLE]

Similarly it can be shown that $\boldsymbol{y}_{h}(u_{i},\cdot)\cdot\tilde{\boldsymbol{\tau}}_{h}(u_{i},\cdot)=0$ for $i=1,\ldots,N$.

Lemma 5**.**

If $\alpha^{0},\beta^{0},\gamma^{0}$ are independent of time, any solution to the above problem satisfies:

[TABLE]

Proof.

First we see that (3.3g) implies that $\left|\boldsymbol{x}_{h,u}\right|_{t}=0$.

We test 3.3a with $\boldsymbol{x}_{h,t}$, 3.3d with $-m_{h}$ and 3.3g with $p_{h}$. Adding the resulting equations results in

[TABLE]

Next, we take a time derivative of 3.3c and test the result with $\boldsymbol{y}_{h}$.

[TABLE]

From the result we subtract 3.3b tested with $\boldsymbol{\kappa}_{h,t}$ to see

[TABLE]

Using 3.6, the frame equations 3.4 and some simple vector identities gives

[TABLE]

We can further reduce the right hand side using the definition of $\boldsymbol{\omega}_{h}$ and the fact that $\boldsymbol{\kappa}_{h}\cdot\tilde{\boldsymbol{\tau}}_{h}=0$:

[TABLE]

Similarly, we see that

[TABLE]

Combining 3.10, 3.11 and 3.12 with 3.9 gives

[TABLE]

We identify the last two terms of the right-hand side with terms in 3.3b tested with the test function $\boldsymbol{\psi}_{h}$ given by

[TABLE]

Noting that $\boldsymbol{y}_{h}\cdot\boldsymbol{\psi}_{h}=\boldsymbol{y}_{h}\cdot(m_{h}\tilde{\boldsymbol{\tau}}_{h}\times\boldsymbol{\kappa}_{h})$, we infer that:

[TABLE]

We sum the result of testing 3.3e with $-\gamma_{h,t}$ and 3.3f with $z_{h}$ and rearrange:

[TABLE]

Adding 3.9, 3.13 and 3.14 gives the desired result. ∎

3.3. Fully discrete problem

To discretize in time we use a uniform partition of the time interval $[0,T]$ into time steps $0=t_{0}<t_{1}<\ldots<t_{M}=T$ where $t_{n}=n\Delta t$ for $n=0,\ldots,M$. We denote discrete variables at a time step $t_{n}$ with a superscript $n$ and the frame vectors will be denoted by $\boldsymbol{e}^{1,n}_{h}$. Our approach is a first order semi-implicit scheme which results in a linear problem to solve at each time step. During the numerical experiments, we demonstrate that by choosing to take certain terms implicitly we recover the semi-discrete stability result. For a variable $\eta$ defined at each time step $0<n<M$, we denote the backward difference $\bar{\partial}\eta^{n}:=(\eta^{n}-\eta^{n-1})/\Delta t$. We define $\boldsymbol{\kappa}_{b}^{n}$ by

[TABLE]

As well as choosing whether to take terms implicitly or explicitly, we also integrate the constraint equation 3.3g forwards in time and write the frame update equation as an algebraic relation which preserves the nature of the angular velocity vector $\boldsymbol{\omega}_{h}$.

Problem 6**.**

Given preferred curvatures $\alpha^{0},\beta^{0}$, a preferred twist $\gamma^{0}$, and initial conditions for the parametrization $\boldsymbol{x}_{h}^{0}$ and the frame $\boldsymbol{e}^{1,0}_{h},\boldsymbol{e}^{2,0}_{h}$, first define $\boldsymbol{\kappa}_{h}^{0}\in V_{h,0}^{3}$ and $\gamma_{h}^{0}\in V_{h}$ as the solutions of

[TABLE]

Then for $n=1,\ldots,M$ find $\boldsymbol{x}_{h}^{n}\in V_{h}^{3},\boldsymbol{y}_{h}^{n}\in V_{h,0}^{3},\boldsymbol{\kappa}_{h}^{n}\in V_{h,0}^{3}+\boldsymbol{\kappa}_{b}^{n}$, $m_{h}^{n}\in V_{h},z_{h}^{n},\gamma_{h}^{n},p_{h}^{n}\in Q_{h}$, $\boldsymbol{e}^{1,n}_{h},\boldsymbol{e}^{2,n}_{h}\in V_{h}^{3}$ such that

[TABLE]

The scheme results in a linear system of equations at each time step followed by an algebraic update formula for the frame. The linear system can be solved using a direct sparse solver. In this work, the numerical results are computed using the UMFPACK library [25].

In addition to the usual time discretization, we have chosen to integrate the constraint equation forwards in time. This gives us more control over the length element $|\boldsymbol{x}_{h,u}^{k}|$ as shown in the following lemma.

Lemma 7**.**

If there exists a solution such that $\left|\boldsymbol{\tau}_{h}^{n}-\boldsymbol{\tau}_{h}^{n-1}\right|^{2}<2$, then

[TABLE]

Proof.

Testing equation 3.15g with $q_{h}=\chi_{[u_{i},u_{i+1}]}$, the characteristic function of the interval $[u_{i},u_{i-1}]$, gives the element-wise identity

[TABLE]

Then we have

[TABLE]

Since $\frac{1}{2}\left|\boldsymbol{x}_{h,u}^{n}\right|\left|\boldsymbol{\tau}_{h}^{n}-\boldsymbol{\tau}_{h}^{n-1}\right|^{2}\geq 0$, we have $\left|\boldsymbol{x}_{h,u}^{n}\right|\geq\left|\boldsymbol{x}_{h,u}^{0}\right|$. Furthermore if it holds that $\left|\boldsymbol{\tau}_{h}^{n}-\boldsymbol{\tau}_{h}^{n-1}\right|^{2}<2$, this equation can be rearranged to see the desired result. ∎

The first rotation 3.15i maps $\tilde{\boldsymbol{\tau}}_{h}^{n-1}$ to $\tilde{\boldsymbol{\tau}}_{h}^{n}$ and the second rotation, 3.15j, rotates the frame about the new $\tilde{\boldsymbol{\tau}}_{h}^{n}$ (leaving $\tilde{\boldsymbol{\tau}}_{h}^{n}$ unaffected). Since we apply the same rotations to the two frame vectors, and these rotations map $\tilde{\boldsymbol{\tau}}_{h}^{n-1}$ to $\tilde{\boldsymbol{\tau}}_{h}^{n}$, this update procedure results preserves the orthogonality of the frame vectors at each vertex. In practical computations we will see an accumulation of floating point errors (see Fig. 7, for example). If the errors become too large we may renormalize the frame and continue the computation.

4. Results

We provide three test cases for our numerical scheme. In the first we relax a straight rod to one with prescribed curvatures and twist and in the other two we demonstrate the applicability of the method to C. elegans locomotion.

4.1. Relaxation test

We take as initial configuration of the rod a unit length straight midline curve and constant frame. We then simulate to $T=25$ with

[TABLE]

We take material parameters all equal to 1: $L=K_{\mathrm{rot}}=A=B=C=D=1$, $\mathbb{K}=\mathbb{I}$. An example of the final configuration is shown in Fig. 2. We will use this example to show how the stability result (Lem. 5) translates to the discrete case. We also explore the errors in the length element (Lem. 7) and failure to preserve exact orthogonality of the frame. To show the properties of the scheme, we first simulate with $\Delta t=1$ and $N=16$ and repeat with the time step $\Delta t$ reduced by a factor of four and doubling $N$: We simulate $\Delta t=4^{-l}$ and $N=2^{4+l}$ for refinement levels $l=0,1,\ldots,5$.

To investigate the fully discrete stability of the scheme, the elastic energy $\mathcal{E}(t^{n})$ is shown at each time step $t^{n}$ in Fig. 3 . We define $\mathcal{E}(t^{n})$ by

[TABLE]

We see that across all configurations the energy decreases across all our simulation results apart from the very coarsest time steps.

Next, we look at the error in the length element. We follow the refinement procedure detailed above and show the results in Tab. 2 and Fig. 4. The error shown is

[TABLE]

The experimental order of convergence, $eoc(\Delta t)$, is

[TABLE]

We observe that the error in length element decrease in time after the increase from the first initially perfect time step. This matches with the analysis of Lem. 7 that the error only depends on the change in tangent from one time step to the next. Since the scheme converges to a stable solution the change in tangent vector reduces in time which results in the reduction of error. Moreover, we see that the error reduces to second order in the time step which is an order higher than the expected error in the scheme overall.

Next, we test for errors in the frame orthogonality conditions with results shown in Tab. 3 and Fig. 4. Here, we look at the maximum over time of the $L^{2}$-norm of the errors in the orthogonality conditions:

[TABLE]

We observe that the errors are very small across all simulations although the error does increase as we refine in space and time. Furthermore we see that the errors increase over time which we attribute to accumulation of rounding errors. The final column in Tab. 3 shows the maximum change in this error over time. This value is close to machine precision epsilon which indicates that the increase in errors, both in time and as we refine the time step, is due to an increase in the number of time steps.

4.2. Application to nematode locomotion in two and three spatial dimensions

We augment the method detailed above by changing the linear drag term to a resistive force term [51]:

[TABLE]

which we approximate in the fully discrete scheme by

[TABLE]

We demonstrate that we can use the method to simulate C. elegans locomotion in two and three dimensions. We set $L=1$ and restrict our considerations to a stiff environment with $K=40,K_{\mathrm{rot}}=1$ which should correspond to a crawling behaviour. We model the C. elegans body as an elastic tapered cylinder. We assume that the internal viscous forces do not play a role in the stiff environment [11] so we consider ($B=D=0$). For material parameters, we take $A=C=8((\varepsilon+u)(\varepsilon+1-u))^{3/2}/(1+2\varepsilon)^{3}$, for $\varepsilon>0$ small, which corresponds to a uniform elasticity across the shell of a tapered body shape. We assume that the frame directions correspond to physically meaningful directions within the worm. We assume that $\boldsymbol{e}^{0}$ follows the midline of the body pointing head to tail, $\boldsymbol{e}^{1}$ points in the ventral-dorsal plane - the usual bending direction when considering two dimensional locomotion - and that $\boldsymbol{e}^{2}$ points in the left-right direction. Muscle contractions generate bending in either the $\boldsymbol{e}^{1}$ or $\boldsymbol{e}^{2}$ directions.

In the usual two dimensional scenario, C. elegans generates bending waves in the dorsal-ventral plane. This will be our first test case:

[TABLE]

We have previously seen [20] that the first condition can be seen to recreate realistic looking C. elegans locomotion postures. We explore here how well our updated numerical method captures this behaviour. Further results for this test case are shown in App. A.

It is assumed that C. elegans generates undulations in the dorsal-ventral plane due to symmetries in its neural control. However these symmetries do not exist in the head and neck regions (see the discussion in [13]). Therefore, we propose an alternative control strategy which results in three dimensional, non-planar body configurations. Using the notation $\chi_{[0,1/3]}$ for the characteristic function of the interval $[0,1/3]$, we simulate with

[TABLE]

For both cases, as initial condition we start with an initially straight rod with constant frame and simulate with $\alpha^{0}(u,0),\beta^{0}(u,0),\gamma^{0}(u,0)$ until $t=5$ and use the resulting curve as initial condition for our simulation. This means we have an initial condition where the curvatures and twist match the initial preferred curvatures and twist exactly however the frame orthogonality conditions will not hold exactly.

We show some characteristic body positions in Fig. 5(a) and head trajectories in Fig. 5(b). We note that simulations for the two dimensional case we have that the third component is zero and the twist is exactly zero (Fig. 5(c)) whereas there is some twist in the three dimensional scenario even though the preferred twist is zero (Fig. 5(d)). Further the three dimensional test case demonstrates a non-planer body position and trajectory (see Fig. 5(b) and B). The second test case demonstrates a possible strategy for generating three dimensional postures. We check the errors in the length element and frame orthogonality conditions with results shown in Fig. 6, 6(a) and 6(b) for the two dimensional case and Fig. 7, 7(a) and 7(b) for the three dimensional case. We observe similar results to the relaxation case. We see the same second order convergence in the error in the length element although now this error increases and decreases periodically depending on the periodic undulations. The error is higher overall since the midline continues to move throughout the simulation. We see that the frame mismatch is again very small across all simulations. This error is initially higher since the initial conditions are derived from simulations.

5. Discussion

We have presented a new finite element scheme for viscoelastic rods suitable for simulating undulatory locomotion in three dimensions. We have shown analytically that the semi-discrete problem preserves the energy structure of the continuous problem and that the scheme preserves geometric constraints exactly (up to machine precision in practical computations). Further our numerical experiments show, first, that the analytic properties are realized in a practical fully discrete scheme and that we can capture simple two and three dimensional undulatory behaviours. Although our method is not the first to tackle this problem, we believe the theoretical grounding and our numerical experiments demonstrate that our method is very well suited for simulations of undulatory locomotion.

We hope in future work to be able to extend the analysis to show stability at the fully discrete level or other convergence properties. Stability of similar schemes for second order geometric flow curve shortening flow based on parametric representations has been shown at the fully discrete level often by using a non-degenerate re-parameterization of the underlying partial differential equation (e.g. [31, 3, 33] and references therein). The basis of our numerical scheme are the papers [32, 63] for which there exists numerical analysis of the semi-discrete problem [26] but not the fully discrete problem. There is however further analysis available when alternative parametric schemes are used. The author of [5] uses a local length constraint for the elastic curve problem in a finite element scheme using Bézier curves. The resulting scheme is unconditionally stable at the discrete level.

In order to apply the method to derive insights for undulatory locomotion, especially for C. elegans, we must try to understand both material parameters and control mechanisms. The discussion in [20] of material parameters for a two-dimensional model still holds true and in three dimensions there are further parameters to explore. We must also try to understand what gaits are possible and try to find their underlying mechanisms. For example, the authors of [13] propose a roll manoeuvre in C. elegans called a “torsion turn”. Work is required to understand what control mechanisms would allow such a behaviour.

There are also many details which we have chosen not to consider in our model. The largest simplification is the use of resistive force theory (RFT). RFT is a simple first order approximation to the full fluid dynamics which could appear in our model. The drag coefficients have been experimentally verified for C. elegans simple forwards locomotion in a variety of media [11], although other studies have suggested nonlinear, but still local, drag coefficients give a more accurate representation of the underlying physics [69]. Alternative approaches include using slender body theory [51, 60] which gives a second order non-local correction to RFT or solving the full fluid structure interaction problem using embedded filaments (e.g. [22, 70], a boundary element method [64] or an immersed boundary method [62]. However the RFT approach is crucial in allowing us to solve such a simple one-dimensional system of equations but the model does not hold when non-local interactions become important such as in turning manoeuvres, mating or undulations near walls. The model and numerical scheme would need to be extended in order to capture these situations.

Acknowledgements

The author would like to thank Netta Cohen and Felix Salfelder for discussions which have improved the writing of the manuscript and the computer implementation of the method. I would also like to thank the anonymous reviewers for their careful consideration and helpful comments.

Declarations of interest

None

Appendix A A related two dimensional problem

The numerical method presented in the main part of this paper may also be used to understand a related two dimensional scenario. We consider the problem of a three dimensional rod embedded in a plane without twist - this is the relevant model for nematode locomotion when restricted to two dimensions where the nematode lies on a plate on its left or right side. This can be used to develop a simpler numerical scheme for these scenarios as we no longer need to solve for the frame and can simply compute the evolution of the midline only. Indeed we demonstrate numerically that we recover the same results either using the numerical scheme in the main paper or the restricted version here. The restricted two dimensional scheme here also has the advantage of reducing the numbers of degrees of freedom. In practice this leads to a speed up of simulations by at least a factor of 2 for fine discretizations.

This two dimensional method extends the method previously presented in [20]. Here, we include the viscous as well as elastic contribution to the moment.

A.1. The model

We consider the same model as in Sec. 2 except now we impose that the parameterization of the midline lies in the plane $\boldsymbol{x}\in\mathbb{R}^{2}\times\{0\}$ and the frame has zero twist ($\gamma\equiv 0$). For the orthonormal frame, we can use $\boldsymbol{e}^{0}\equiv\boldsymbol{\tau},\boldsymbol{e}^{1}\equiv\boldsymbol{\nu}:=\boldsymbol{\tau}^{\perp},\boldsymbol{e}^{2}:=(0,0,1)$, where $(\cdot)^{\perp}$ denotes rotation by $\pi/2$ about $\boldsymbol{e}^{2}$. This orthonormal frame has zero twist ($\gamma\equiv 0$) and we have zero curvature in the $\boldsymbol{e}_{2}$ direction $(\beta\equiv 0)$. Therefore we can derive an appropriate model by replacing 2.6 by

[TABLE]

More compactly, we can write

[TABLE]

Then 2.15 may be replaced by a two dimensional version:

[TABLE]

A.2. Weak form

The weak form of our two dimensional problem is:

Problem 8**.**

Given a preferred curvature $\alpha^{0}$ and an initial condition for the parameterization $\boldsymbol{x}^{0}$, find $\boldsymbol{x},\boldsymbol{y},\boldsymbol{\kappa}\colon[0,1]\times[0,T]\to\mathbb{R}^{2}$ (with the conditions 2.15d, 2.15e and 2.17 at the boundaries), $p\colon[0,1]\times[0,T]\to\mathbb{R}$, such that, for almost every $t\in(0,T)$:

[TABLE]

subject to the initial condition

[TABLE]

and initial equation

[TABLE]

A.3. Numerical method

A.3.1. Semi-implicit scheme

Our approach follows the same steps as Sec. 3. We again use the spaces of piecewise linear and piecewise constant spaces $V_{h}$ and $Q_{h}$ and the two approximations of the tangent vector 3.1 and 3.2. We use the second of these definitions to define an approximation to the normal vector $\boldsymbol{\nu}_{h}:=\tilde{\boldsymbol{\tau}}_{h}^{\perp}$. We make the same choice of approximation spaces as before and arrive at the scheme:

Problem 9**.**

Given a preferred curvature $\alpha^{0}$, and an initial condition for parameterization $\boldsymbol{x}_{h}$, for $t\in[0,T]$, find $\boldsymbol{x}_{h}(\cdot,t)\in V_{h}^{2},\boldsymbol{y}_{h}(\cdot,t)\in V_{h,0}^{2},\boldsymbol{\kappa}_{h}(\cdot,t)\in V_{h,0}^{2}+\boldsymbol{\kappa}_{b}(\cdot,t)$, $p_{h}(\cdot,t)\in Q_{h}$, $\boldsymbol{\nu}_{h}(\cdot,t)=\tilde{\boldsymbol{\tau}}_{h}^{\perp}\in V_{h}^{2}$ such that, for all $t\in(0,T)$:

[TABLE]

subject to the initial condition:

[TABLE]

and initial equation:

[TABLE]

The stability result can be shown in the same way. We do not show the proof here.

Lemma 10**.**

If $\alpha^{0}$ is independent of time, any solution to the above problem satisfies:

[TABLE]

A.3.2. Fully discrete scheme

Problem 11**.**

Given a preferred curvature $\alpha^{0}$, and an initial condition for the parameterization $\boldsymbol{x}_{h}^{0}$, first define $\boldsymbol{\kappa}_{h}^{0}\in V_{h,0}^{2}$ and as the solutions of

[TABLE]

Then, for $n=1,\ldots,M$, find $\boldsymbol{x}_{h}^{n}\in V_{h}^{2},\boldsymbol{y}_{h}^{n}\in V_{h,0}^{2},\boldsymbol{\kappa}_{h}^{n}\in V_{h,0}^{2}+\boldsymbol{\kappa}_{b}^{n}$, $p_{h}^{n}\in Q_{h}$ such that

[TABLE]

We see that the treatment of the length constraint has not changed so we immediately recover the result of Lem. 7. Furthermore, we can find the frame without solving any more equations so that we preserve the frame orthogonality exactly.

A.4. Numerical results

We demonstrate the similarities between the two and three dimensional versions of the models by simulating the two dimensional test case from the main paper. We use the same refinement strategy as the main paper also. First we compute the centre of mass at the final time step for the two dimensional scheme A.8 and the three dimensional scheme 3.15 and show the norm of the difference in Tab. 6. We see that the error at the final time is very small and is approximately machine precision epsilon per time step. This demonstrates that the schemes are exactly the same up to rounding errors. Second, we compare the how long each simulation takes to run taking the average run time and range across ten simulations. This includes the 5 time unit initial simulation in order to generate the initial condition. We see that for reasonable levels of resolution the two dimensional version of the code runs in less than half the time of the three dimensional version.

Appendix B Simulation videos

We also attach three videos of simulations presented in the main paper:

**Video 1.: **

Conformation of the rod during the relaxation test with $\Delta t=10^{-3}$ and $N=128$ for times in $[0,25]$. The colouring is the same as Fig. 2. https://www.youtube.com/watch?v=aih4-yPn0lk

**Video 2.: **

Conformation of the body during the two dimensional locomotion test with $\Delta t=10^{-3}$ and $N=128$ for times in $[0,25]$. The colouring is the same as Fig. 5(a)(i). https://www.youtube.com/watch?v=mWtLNb93RGY

**Video 3.: **

Conformation of the body during the three dimensional locomotion test with $\Delta t=10^{-3}$ and $N=128$ for times in $[0,25]$. The colouring is the same as Fig. 5(a)(ii). https://www.youtube.com/watch?v=FO8opOIykLs