A Linearised $\mathit{\theta }$ Numerical Scheme for the Vibrations of Inextensible Beams

Theodosios K. Papathanasiou

Department of Civil and Environmental Engineering, Brunel University London, Uxbridge UB8 3PH, UK

E-mail: theodosios.papathanasiou@brunel.ac.uk

Received 26 April 2021; Accepted 22 June 2021; Publication 23 July 2021

Abstract

A linearised finite element numerical scheme for the vibration of inextensible beams is developed. The proposed scheme is based on the methodology introduced by S. Bartels [15] and satisfies a linearised form of the inextensibility constraint. The time m arching procedure is based on repeated use of the theta-parameter integration quadrature. Three parameters are introduced in total and appropriately selected such that the energy conservation features are improved compared to the Bartels algorithm while the inextensibility constraint is satisfied as accurately as possible. Cubic Hermite polynomials are employed for the spatial discretisation. The Bartels algorithm is retrieved as a special case. Several numerical experiments are presented demonstrating the theoretically predicted enhanced inextensibility mimicking and optimum values of the method parameters are identified.

Keywords: Flexural vibrations, inextensibility, finite element method, large deflections, time-integration, structural elements.

1 Introduction

Flexural deformations and flexural vibrations of beams is a core subject of structural mechanics with numerous applications. Models for the large deflection of beams and columns can be used to study instabilities of structures [1, 2], buckling related phenomena [3], biomimetic motion such as snake locomotion [4, 5] and compliant mechanisms in soft-robotics applications [6–8]. Inextensible beams form a large subcategory of this type of systems with equally important applications, such as modelling the mechanical behaviour of textiles [9], the response of flexible risers in the oil & gas industry [10, 11] and the properties of thin filaments in composite materials [12].

To model the large amplitude flexural vibrations of inextensible, incompressible beams, the structural component displacement can be simulated as a curve and be arc-length parameterised [13, 14]. For dynamic analysis let $T>0$, $a<b$ and introduce the arc-length parameter $s\in I=[a,b]$ that defines the position of each point along the beam during the evolution in a time interval $t\in [0,T]$. For $s\in I$, the vector valued function $u(s,t)={\left[\begin{array}{ccc}\hfill x(s,t)\hfill & \hfill y(s,t)\hfill & \hfill z(s,t)\hfill \end{array}\right]}^T$, defines the location in ${ℝ}^3$ of each point along the length of the inextensible beam. The curvature at each point $s$ is $\kappa (s)=\sqrt{u_{ss}\cdot u_{ss}}=|u_{ss}|$, where $u_s={\partial }_{s}u$. The osculating plane of the curve is defined by the unit vectors $\mathbf{\text{T}}=u_s/|u_s|$ and $\mathbf{\text{N}}={\kappa }^{-1}{u}_{ss}$, while the Frenet-Serret binormal vector, perpendicular to the osculating plane, is $\mathbf{\text{B}}={\kappa }^{-1}(u_s\times u_{ss})=\mathbf{\text{T}}\times \mathbf{\text{N}}$ (Figure 1).

Figure 1 Frenet vectors and osculating plane for the large displacement vibrations of inextensible beams.

The kinetic energy of the flexible beam is

$$E_K({\partial }_{t}u)=\frac{1}{2}{\int }_{I}m{|{\partial }_{t}u|}^2\mathit{d}s,$$

(1)

where $m$ is the mass distribution along the beam length. Adopting Euler-Bernoulli assumptions for bending of slender beams, the strain energy of the beam has the form

$$E_P(u)=\frac{1}{2}{\int }_{I}EI{|\kappa |}^2\mathit{d}s,$$

(2)

where $EI$ denotes the flexural rigidity. The inextensibility of the component is expressed through the no length change condition [13, 15]

$${|u_s|}^2-1=0.$$

(3)

For planar problems, the position vector has two components, while both the 2D and 3D position vector valued functions are denoted as $u(s,t):I\times [0,T]\to {ℝ}^d$, with $d=2$ or $d=3$.

In the following, boundary conditions which prescribe displacement, slope, or zero bending moment and shear force at the end points of the beam are considered. Periodic boundary conditions can be applied as well [15]. Introducing a characteristic length $\mathrm{\ell }$ and the nondimensional quantities $\tilde{s}=s/\mathrm{\ell }$, $\tilde{u}=u/\mathrm{\ell }$ and $\tilde{t}={\mathrm{\ell }}^{-2}t\sqrt{EI/m}$, the variational form of the beam vibration problem is (after dropping tildes)

$${\int }_I{\partial }_{tt}u\cdot \chi ds+{\int }_I{u}_{ss}\cdot {\chi }_{ss}\mathit{d}s+{\int }_I\mathrm{\Lambda }{u}_s\cdot {\chi }_s\mathit{d}s={\int }_{I}Q\cdot \chi \mathit{d}s,$$

(4)

for all admissible weight functions $\chi$, where $\mathrm{\Lambda }=\lambda {\mathrm{\ell }}^3/EI$ is the nondimensional Lagrange multiplier associated with the inextensibility constrain and $Q=q{L}^3/EI$ is the nondimansional expression of the externally applied load $q$. The spatial domain in the nondimensional setting is $I=[\alpha ,\beta ]=[a/\mathrm{\ell },b/\mathrm{\ell }]$. Variational form (4) applies for all functions $u$ such that ${|u_s|}^2-1=0$ in $I$. The Initial-Boundary Value problem features initial conditions on position and velocity at time $t=0$, of the form

$$u(s,0)=u_o(s),u_t(s,0)=v(s,0)=v_o(s),$$

(5)

To incorporate the inextensibility constraint, the procedure proposed in [15, 16] for the simulation of inextensible curves will be adopted. In particular, the inextensibility constraint will be used and the inextensibility constraint will be applied only at the mesh nodes. However, in the present study a variation of the approximation used in [15] will be considered. Differentiating the inextensibility constraint ${|u_s|}^2=1$ with respect to time, results to

$$u_s\cdot v_s=0,$$

(6)

where $v_s={\partial }_t{u}_s$. Assuming sufficient regularity, for time increments $\mathrm{\Delta }t\ll 1$ it is

$$(u_s(t)+\mathrm{\Delta }t{v}_s(t)+O(\mathrm{\Delta }{t}^2))\cdot v_s(t+\mathrm{\Delta }t)=0.$$

(7)

Equation (7) will be the basis for the implementation of the inextensibility constraint in a linearised $\theta$ scheme. The introduction of theta parameters in the approximation of integrals for time-marching procedures is expected typically to produce versatile methods with improved discrete energy conservation properties and accuracy [17].

In the following, the standard notation $H^2(I;{ℝ}^d)$ is used to denote the Sobolev (Hilbert) space $W^{2,2}(I;{ℝ}^d)$ of functions defined in $I$ with values in ${ℝ}^d$ and $(,),\parallel \parallel$ denote the underlying $L^2$ space inner product and induced norm respectively. For reasons of consistency, the presentation and notation follows in general the one introduced in [15].

2 The Linearised $\mathit{\theta }$ Scheme

For a given time-step $\tau =t_{n+1}-t_n>0$, integrating Equation (4) in $[t_n,t_{n+1}]$, results to

$$(v^{n+1},\chi )-(v^n,\chi )+{\int }_{t_n}^{t_{n+1}}(u_{ss},{\chi }_{ss})\mathit{d}t={\int }_{t_n}^{t_{n+1}}(Q,\chi )\mathit{d}t,$$

(8)

where the velocity at instant $n+1$ is denoted as $v^{n+1}$. Assuming that the time step is sufficiently small, the time integral on the left hand side is now approximated using a $\theta$ scheme as

$${\int }_{t_n}^{t_{n+1}}(u_{ss},{\chi }_{ss})\mathit{d}t\approx (1-{\theta }_1)\tau (u_{ss}^n,{\chi }_{ss})+{\theta }_1\tau (u_{ss}^{n+1},{\chi }_{ss}),$$

(9)

with the parameter ${\theta }_1\in [0,1]$. At the same time, the update for vector $u$ is assumed in the form

$$u^{n+1}=u^n+\tau [(1-{\theta }_2)v{}^n+{\theta }_{2}v{}^{n+1}],$$

(10)

where a second parameter ${\theta }_2\in [0,1]$ is introduced. Using now the approximation (10) for $u_{ss}^{n+1}$ in Equation (9) it is

$(1-{\theta }_1)\tau (u_{ss}^n,{\chi }_{ss})+{\theta }_1\tau (u_{ss}^{n+1},{\chi }_{ss})$

$\mathrm{\hspace{1em}\hspace{1em}}=\tau (u_{ss}^n,{\chi }_{ss})+{\tau }^2[{\theta }_1(1-{\theta }_2)(v_{ss}^n,{\chi }_{ss})+{\theta }_1{\theta }_2(v_{ss}^{n+1},{\chi }_{ss})].$

The inextensibility constraint at step $n+1$ is $v_s^{n+1}\cdot u_s^{n+1}=0$. The key idea at this point is that for a converging procedure, the value $u_s^n+{\theta }_3\tau v_s^n$ can be used as a good approximation for $u_s^{n+1}$. Regarding parameter ${\theta }_3\in [0,1]$, the value ${\theta }_3=0$ will lead to the linearised approximation introduced by S. Bartels in 2016 [15]. Selecting the value ${\theta }_3=1$ results in an expression of the form $u_s^{n+1}=u_s^n+\tau v_s^n$. For smooth evolution phenomena, the velocity derivative with respect to the arc length parameter is expected to satisfy an expansion of the form

$$u_s^{n+1}=u_s^n+\tau v_s^n+\frac{1}{2}{\tau }^2{a}_s^n+O({\tau }^3),$$

(12)

where $a^n$ is the acceleration at $t=t_n$. Therefore, if the velocity field is sufficiently smooth, the formula

$$u_s^{n+1}\approx u_s^n+\tau v_s^n,$$

(13)

is expected to provide a higher order approximation for $u_s^{n+1}$ with principal error term $O({\tau }^2)$. Hence, the proposed linearised form of the inextensibility constraint reads

$$v_s^{n+1}\cdot (u_s^n+{\theta }_3\tau v_s^n)=0\text{in }I.$$

(14)

The forcing integral in (8) can be evaluated analytically for simple forms of the $Q$ function or using numerical integration. For the fully discrete scheme, we define the function space to be used in spatial discretisation as

$V_{h,0}^d=\left\{\begin{array}{c}{v}^h\in H^2(I;{ℝ}^d):{v^h|}_{K_i}\in P^3(K_i),i=1,2,3,\mathrm{\dots },\mathrm{M}+1\hfill \\ \mathrm{and}{v}^h\mathrm{satisfies}\mathrm{homogeneous}\mathrm{essential}\mathrm{boundary}\mathrm{conditions}\hfill \\ \mathrm{at}x=a\mathrm{or}x=b,\mathrm{where}\mathrm{applicable}\hfill \end{array}\right\}.$ (15)

By definition it is $V_{h,0}^d\subset H^2(I;{ℝ}^d)$ and Hermite cubic polynomials, conforming to the function properties of the space defined in (2) are used as shape and weight functions in the fully discrete scheme.

Following the analysis in [15, 16, 18], the constraint equations are applied only at the mesh nodes. Then, the fully discrete algorithm is:

Given $u^{h,j},v^{h,j}\in V_{h,0}^d$ for $j=0,1,2,\mathrm{\dots },n$ $\tau >0$ and parameters ${\theta }_1,{\theta }_2,{\theta }_3\in [0,1]$, find $v^{h,n+1}\in V_{h,0}^d$ such that

$(v^{h,n+1},{\chi }^h)-(v^{h,n},{\chi }^h)+\tau (u_{ss}^{h,n},{\chi }_{ss}^h)$
$\mathrm{\hspace{1em}\hspace{1em}}\mathit{\hspace{1em}}+{\tau }^2[{\theta }_1(1-{\theta }_2)(v_{ss}^{h,n},{\chi }_{ss}^h)+{\theta }_1{\theta }_2(v_{ss}^{h,n+1},{\chi }_{ss}^h)]$
$\mathrm{\hspace{1em}\hspace{1em}}={\displaystyle {\int }_{t_n}^{t_{n+1}}}(Q,{\chi }^h)\mathit{d}t,\text{and}$	(16)

$v_s^{h,n+1}\cdot (u_s^{h,n}+{\theta }_3\tau v_s^{h,n})=0$ in $I$, for all ${\chi }^h\in V_{h,0}^d$, with ${\chi }_s^h\cdot (u_s^{h,n}+{\theta }_3$ $\tau v_s^{h,n})=0$ in $I$. Update $u$, using $u^{h,n+1}=u^{h,n}+\tau [(1-{\theta }_2)v{}^{h,n}+{\theta }_{2}v{}^{h,n+1}]$.

The method introduced by S. Bartels (2016) is retrieved using the parameter values ${\theta }_1={\theta }_2=1$ and ${\theta }_3=0$. Analysis of the energy conservation characteristics for the modified $\theta$ method indicates that for ${\theta }_2=1/2$ the energy conservation properties of the method will improve as ${\theta }_1\to 1/2$. This is expected to some extend due to the trapezoidal integration rule employed and also verified by the result in the following section. Selecting the value ${\theta }_3=1$ for the constraint equation is expected to lead to better approximation for the $u_s^{n+1}$ value when smooth fields are considered since for $\tau \to 0$ it should be $u_s^{n+1}=u_s^n+\tau v_s^n+O({\tau }^2)$. Analysis of the error introduced due to the inextensibility constraint is included in the following subsection. In the following, two cases will be considered as summarised in the following table (Table 1).

Table 1 Parameter values for the S. Bartels method (2016) and the linearised $\theta$ scheme

Parameter Values	Method
${\theta }_1={\theta }_2=1$, ${\theta }_3=0$	S. Bartels Method 2016
${\theta }_1\in (1/2,1)$, ${\theta }_2=1/2$, ${\theta }_3=1$	Linearised $\theta$ method

3 Properties of the $\mathit{\theta }$ Scheme

In this section the discrete energy conservation and inextensibility violation properties of the linearized $\theta$ method will be studied.

3.1 Energy Conservation Properties

The discrete kinetic and strain energy values at time step $n$ are

$$E_K^{h,n}=\frac{1}{2}{\parallel v^{h,n}\parallel }^2\mathit{\hspace{1em}}\text{and}\mathit{\hspace{1em}}{E}_P^{h,n}=\frac{1}{2}{\parallel u_{ss}^{h,n}\parallel }^2.$$

(17)

The values ${\theta }_1\in (1/2,1)$, ${\theta }_2=1/2$ and ${\theta }_3=1$ are used for the modified method. Setting ${\chi }^h=\frac{\tau }{2}(v{}^{h,n}+v{}^{h,n+1})$ in Equation (16) and using the update definition in the form $v{}^{h,n+1}=\frac{2}{\tau }(u^{h,n+1}-u^{h,n})-v{}^{h,n}$ the discrete energy balance for $Q=0$ is

$$E_K^{h,n+1}+E_P^{h,n+1}=E_K^{h,n}+E_P^{h,n}+\frac{1-2{\theta }_1}{2}{\parallel u_{ss}^{h,n+1}-u_{ss}^{h,n}\parallel }^2.$$

(18)

Using now again the update definition (10) it is

$$E_K^{h,n+1}+E_P^{h,n+1}=E_K^{h,n}+E_P^{h,n}+\frac{(1-2{\theta }_1){\tau }^2}{8}{\parallel v_{ss}^{h,n+1}+v_{ss}^{h,n}\parallel }^2.$$

(19)

Since ${\theta }_1\in (1/2,1)$ it is $1-2{\theta }_1=-|2{\theta }_1-1|$ and Equation (19) indicates that numerical dissipation of energy occurs. Iterating Equation (19) results in

$E_K^{h,n+1}+E_P^{h,n+1}=E_K^{h,0}+E_P^{h,0}-{\displaystyle \frac{|2{\theta }_1-1|{\tau }^2}{8}}{\displaystyle \sum _{j=0}^n}{\parallel v_{ss}^{h,j+1}+v_{ss}^{h,j}\parallel }^2.$

(20)

From Equation (20) it follows that the selection ${\theta }_1=1/2$ leads to exact conservation of the discrete energy. This was expected due to the use of the trapezoidal integration rule. However, extensive numerical experiments indicate that the method is unstable for realistic selection of time step $\tau$ when ${\theta }_1\to 1/2$. As the parameter increases and tends to become one, the energy conservation attribute of the method deteriorates while at the same time the violation in the inextensibility constraint becomes smaller. The conducted numerical experiments indicate that a reasonable choice is ${\theta }_1\in (0.55,0.95)$. The method performance in terms of the inextensibility constraint violation is examined in Subsection 3.2.

3.2 Inextensibility Constraint Violation

In this section the inextensibility constraint approximation will be analysed. It will be shown that for sufficiently small time step the error terms are proportional to ${\tau }^2$ and that furthermore their multiplicative factors attain small values. According to the update $u^{n+1}=u^n+\frac{\tau }{2}(v{}^n+v{}^{n+1})$ it is

$${|u_s^{n+1}|}^2={|u_s^n+\tau v_s^n+\frac{\tau }{2}(v_s^{n+1}-v_s^n)|}^2.$$

(21)

The orthogonality condition $(u_s^n+\tau v_s^n)\cdot v_s^{k+1}=0$, implies that

${|u_s^{n+1}|}^2$	$={|u_s^n+\tau v_s^n|}^2-\tau (u_s^n+\tau v_s^n)\cdot v_s^n+{\displaystyle \frac{{\tau }^2}{4}}{|v_s^{n+1}-v_s^n|}^2,\text{or}$
${|u_s^{n+1}|}^2$	$={|u_s^n|}^2+\tau u_s^n\cdot v_s^n+{\displaystyle \frac{{\tau }^2}{4}}{|v_s^{n+1}-v_s^n|}^2.$	(22)

For $n\ge 1$, using again the orthogonality condition it is

$u_s^n\cdot v_s^n=\left(u_s^{n-1}+\tau v_s^{n-1}+{\displaystyle \frac{\tau }{2}}(v_s^n-v_s^{n-1})\right)\cdot v_s^n={\displaystyle \frac{\tau }{2}}{v}_s^n\cdot (v_s^n-v_s^{n-1}),$(23)

and therefore

$${|u_s^{n+1}|}^2={|u_s^n|}^2+\frac{{\tau }^2}{2}{v}_s^n\cdot (v_s^n-v_s^{n-1})+\frac{{\tau }^2}{4}{|v_s^{n+1}-v_s^n|}^2.$$

(24)

REMARK: The terms responsible for altering the initial length of the curve are proportional to ${\tau }^2$ as was also the case in the method proposed by S. Bartels [15]. Hence, for bounded values of the velocity, the error in the inextensibility constraint is expected to be small as $\tau \to 0_{+}$. In addition to that, the present formulation has an extra important feature. If the evolution is smooth and the time step is sufficiently small, the velocity will not undergo significant changes between successive steps. If this is indeed the case then for $n\ge 1$ $v_s^{j+1}\approx v_s^j\approx v_s^{j-1}$. Noting that the terms proportional to ${\tau }^2$ are multiplied by the differences $v_s^{j+1}-v_s^j,v_s^j-v_s^{j-1}$, the error terms in (24) will become even smaller. It is therefore expected that for smooth vibration phenomena the proposed modified scheme will provide a very good approximation of the inextensibility property. If the velocity changes rapidly between successive steps or if there appear reversals in the velocity sign, then the violation of the inextensibility constraint is expected to be severe. However, for small time steps, the terms responsible for the constraint violation will remain small since they are proportional to ${\tau }^2$.

4 Implementation and Numerical Results

In the following, M denotes the mass matrix produced by applying the shape functions to the inner product $({\chi }^h,{\chi }^h)$ and K the stiffness matrix produced by applying the shape functions to the inner product $({\chi }_{ss}^h,{\chi }_{ss}^h)$. The force vector results from the integral ${\int }_{t_n}^{t_{n+1}}(Q,{\chi }^h)\mathit{d}t$, using again the shape functions as test functions and is denoted as $F_{n\to n+1}^h$.

The method is programmed in MATLAB§and appropriate selection of the $\theta$ parameters is used to reproduce the two schemes in Table 1. Let ${\mathbf{\text{B}}}_n$ be the matrix that corresponds to the nodal inextensibility constraint conditions and appropriate boundary conditions. Denote $U^n,V^n$ the vectors of nodal unknowns and their time derivatives respectively and ${\mathrm{\Lambda }}^n$ the vector of Laplace multipliers. Then, the general implementation form is:

Given $U^n,V^n$ for $j=0,1,2,\mathrm{\dots },n$, $\tau >0$ and parameters ${\theta }_1,{\theta }_2,{\theta }_3\in [0,1]$ set $\mathbf{\text{A}}=\mathbf{\text{M}}+{\theta }_1{\theta }_2{\tau }^2\mathbf{\text{K}}$ and $\mathbf{\text{D}}=\mathbf{\text{M}}-{\theta }_1(1-{\theta }_2){\tau }^2\mathbf{\text{K}}$. Find admissible solution vectors $V^{n+1}$ such that

$\left[\begin{array}{cc}\hfill \mathbf{\text{A}}\mathit{\hspace{1em}}\hfill & \hfill {\mathbf{\text{B}}}_n^T\hfill \\ \hfill {\mathbf{\text{B}}}_n\mathit{\hspace{1em}}\hfill & \hfill O\hfill \end{array}\right]\left[\begin{array}{c}\hfill V^{n+1}\hfill \\ \hfill {\mathrm{\Lambda }}^{n+1}\hfill \end{array}\right]=\left[\begin{array}{c}\hfill \mathbf{\text{D}}{V}^n-\tau \mathbf{\text{K}}{U}^n+{\mathrm{F}}_{n\to n+1}^h\hfill \\ \hfill \mathbf{\text{0}}\hfill \end{array}\right],\text{and}$

(25)

update $U^n$ using $U^{n+1}=U^n+\tau [(1-{\theta }_2)V{}^n+{\theta }_{2}V{}^{n+1}]$.

REMARK: The Schur complement with respect to A of the partitioned matrix in (25) is $\mathbf{\text{W}}=-{\mathbf{\text{B}}}_n{\mathbf{\text{A}}}^{-1}{\mathbf{\text{B}}}_n^T$. Then, the updated solution can be represented as

$V^{n+1}$	$=(\mathbf{\text{I}}+{\mathbf{\text{A}}}^{-1}{\mathbf{\text{B}}}_n^T{\mathbf{\text{W}}}^{-1}{\mathbf{\text{B}}}_n){\mathbf{\text{A}}}^{-1}(\mathbf{\text{D}}{V}^n-\tau \mathbf{\text{K}}{U}^n+{\mathrm{F}}_{n\to n+1}^h),$	(26a)
${\mathrm{\Lambda }}^{n+1}$	$=-{\mathbf{\text{W}}}^{-1}{\mathbf{\text{B}}}_n{\mathbf{\text{A}}}^{-1}(\mathbf{\text{D}}{V}^n-\tau \mathbf{\text{K}}{U}^n+{\mathrm{F}}_{n\to n+1}^h).$	(26b)

In (26a) the term ${\mathbf{\text{A}}}^{-1}(\mathbf{\text{D}}{V}^n-\tau \mathbf{\text{K}}{U}^n+{\mathrm{F}}_{n\to n+1}^h)$ represents the updated state without accounting for the constraint equations, while $\mathbf{\text{I}}+{\mathbf{\text{A}}}^{-1}{\mathbf{\text{B}}}_n^T{\mathbf{\text{W}}}^{-1}{\mathbf{\text{B}}}_n$ is the factor related to constraint effect. In the following subsections several numerical examples are presented and the effect of the $\theta$ parameters is studied.

4.1 Unwinding Helix 1 (S. Bartels 2016)

In this first set of numerical experiments the unwinding helix configuration introduced by S. Bartels [13] is reproduced as a benchmark example. The following parameters, boundary and initial conditions are selected:

$$\mathrm{\ell }=1,I=[0,2\pi ],u(0,t)={\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 1\hfill & \hfill 0\hfill \end{array}\right]}^T,u_s(0,t)={\left[\begin{array}{ccc}\hfill \beta \hfill & \hfill 0\hfill & \hfill \gamma \hfill \end{array}\right]}^T.$$

$$u(0,s)={\left[\begin{array}{ccc}\hfill \sin (\beta s)\hfill & \hfill \cos (\beta s)\hfill & \hfill \gamma s\hfill \end{array}\right]}^T,v(0,s)={\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]}^T,$$

$${\beta }^2=\frac{99}{100},{\gamma }^2=\frac{1}{100},$$

For the first set of results the parameter values for the modified method are set to ${\theta }_1=\theta =0.75$, ${\theta }_2=1/2$, ${\theta }_3=1$. The S. Bartels method is retrieved using ${\theta }_1={\theta }_2=1$, ${\theta }_3=0$. The example is reproduced using 2 different discretization levels. First, a total of $\mathrm{N}=40$ elements with an element size over time-step ratio of $h/\tau =2\pi$. Then, the number of elements used is increase to $\mathrm{N}=80$ while the element time-step ratio is again $h/\tau =2\pi$.

Figure 2 shows the evolution of the discrete system for the mesh with $\mathrm{N}=40$ elements. The linearised $\theta$ scheme corresponds to the continuous black line in Figure 2. The S. Bartels method results are also plotted using a dashed blue line. Small variations between the two solutions are observed as time values increase. The discrete energy (upper row) and change in length (lower row) due to the approximate character of the methods is depicted in Figure 3 for the two discretization levels. The results in Figure 3 demonstrate that the $\theta$ scheme with an appropriate selection of the parameters is more efficient as for the same discretization level it improves the discrete energy conservation properties of the system and leads to significantly smaller violation of the inextensibility constraint. The latter attribute is of course related to the relatively smooth evolution dynamics of the specific example.

To demonstrate the improved behaviour of the proposed method, the performance with the parameters $\mathrm{N}=40,\tau =1/40,\theta =0.75$ is compared to that of the S. Bartels scheme resulting from using the same parameters but also from using the refined mesh $\mathrm{N}=80,\tau =1/80$. The comparison is plotted in Figure 4. It is evident from the right hand side subplot in Figure 4 that the proposed $\theta$ scheme with a less fine mesh produces results comparable to the more refined mesh S. Bartels method output.

Figure 2 Evolution of the unwinding helix in example 4.1 at specific time instances for SB method (blue dashed line) and the proposed modified method with ${\theta }_1=\theta =0.75$ (black line). The total number of elements is $\mathrm{N}=40$ and the time-step $\tau =1/40$.

The performance of the linearised $\theta$ scheme for different values of the parameter in the interval ${\theta }_1=\theta \in (0.55,1)$ is evaluated in Figure 5. The energy conservation ability of the method deteriorates as the parameter increases while at the same time the violation in the inextensibility constrain decreases. Therefore, an optimum balance on the effect on these two criteria needs to be considered. As the time step decreases and since the evolution is relatively smooth, the proposed scheme leads to almost perfect reproduction of the inextensibility (results in the right hand side column of Figure 5) apart from the first few time steps. This can be attributed to the fact that in this example the velocity increases faster during the first stages of the time evolution, while at later stages the velocity between successive steps does not feature significant differences. This result verifies the validity of Equation (24). The method of S. Bartels outperforms the proposed scheme for these few first time-steps. This finding suggests that both methods can be used in conjunction in order to minimise the violation in the inextensibility constraint depending on the magnitude of the acceleration during specific phases of the beam dynamic response.

Figure 3 Discrete energy evolution (upper row) and change in the beam length for example 4.1. S. Bartels method (blue dashed line) and the proposed modified method with $\theta =0.75$ are considered in a relatively coarse and a finer mesh.

REMARK: To suppress the violation in the inextensibility constraint that appears during the first time instances when the modified method is used, the first time step can be performed using ${\theta }_1={\theta }_2=1$ and ${\theta }_3=0$, that is the S. Bartels 2016 method. This option is explored in example 4.2.

Figure 4 Snapshots of the unwinding helix in example 4.1 using $\theta =0.75$ (continuous black line). Evaluation of the proposed scheme using a relatively coarse mesh against the S. Bartels method (blue dashed line) results with a coarse (left) and a finer (right) mesh.

Figure 5 Discrete energy conservation (upper row) and change in total length (lower row) of the unwinding helix in example 4.1 for different values of the parameter ${\theta }_1=\theta$. The number of elements is $\mathrm{N}=40$ and the number of time-steps is doubled for the results in the right hand side column graphs.

4.2 Unwinding Helix 2

To remedy the problematic behaviour of the modified $\theta$ method during the initial time instances, it is proposed that the method is initiated using the S. Bartels scheme for the first step and then switch to the more general $\theta$ scheme. To demonstrate the effectiveness of this approach a second example with the unwinding helix is considered. The geometry of the configuration and initial conditions are the same as in example 4.1 only now the beam is not clamped at the lower end. Both ends are considered to be traction free, i.e. the selected conditions are

$$\mathrm{\ell }=1,I=[0,2\pi ],u(0,s)={\left[\begin{array}{ccc}\hfill \sin (\beta s)\hfill & \hfill \cos (\beta s)\hfill & \hfill \gamma s\hfill \end{array}\right]}^T,$$

$$v(0,s)={\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]}^T,{\beta }^2=\frac{99}{100},{\gamma }^2=\frac{1}{100},$$

In this example, the response will be analysed in the time interval $[0,60]$. For the modified scheme the parameters are ${\theta }_1=\theta =0.6$, ${\theta }_2=1/2$, ${\theta }_3=1$. Figure 6 displays the energy balance (left) and inextensibility constraint violation (right) when the modified scheme is applied. It is again observed that the S. Bartels method outperforms the modified one in the inextensibility property during the first time steps. This is due to the rapid changes in the velocity field that occur initially. In Figure 7, the modified method is initiated with the S. Bartels scheme. It is now observed that the modified method improves significantly with respect to the constraint violation, while retaining to some extend its improved energy preserving attributes. Snapshots of the evolution are presented in Figure 8.

Figure 6 Discrete energy (left) and change in the beam length (right) for example 4.2.

Figure 7 Discrete energy (left) and change in the beam length (right) for example 4.2.

Figure 8 Snapshots of the evolution when the modified method is initiated with the S. Bartels scheme in example 4.2.

4.3 Faulted Ring

In this third example, a circular ring with a discontinuity is analysed. An inextensible thin beam is initially bended in the shape of a circular ring such that the two end points are coinciding, at the lower point of the configuration, but are not fixed one to another. The initial velocity is zero. The ring shaped beam is released from this position at time instant $t=0$. The ring unfolds and encapsulates repeatedly, undergoing flexural vibrations. The parameters selected are

$$I=[-\pi /2,\mathrm{\hspace{0.17em}3}\pi /2],\mathrm{\ell }=1,d=2\text{and}T=110.$$

The boundary conditions are set to zero bending moment and shear force, i.e.

$u_s(-\pi /2,t)$	$=u_s(3\pi /2,t)={\left[\begin{array}{cc}\hfill 0\hfill & \hfill 0\hfill \end{array}\right]}^T,$
$u_{ss}(-\pi /2,t)$	$=u_{ss}(3\pi /2,t)={\left[\begin{array}{cc}\hfill 0\hfill & \hfill 0\hfill \end{array}\right]}^T.$

and the initial conditions are

$$u(s,0)={\left[\begin{array}{cc}\hfill \cos (s)\hfill & \hfill -\sin (s)\hfill \end{array}\right]}^T,v(s,0)={\left[\begin{array}{cc}\hfill 0\hfill & \hfill 0\hfill \end{array}\right]}^T.$$

The spatial domain is discretised using $\mathrm{N}=100$ elements and it therefore is $h=2\pi /100$. A total of 2200 time steps are used for the simulation, defining a time-step of $\tau =110/2200=0.05$. Initially, the ring centre of mass is located at point $(x,y)=(0,0)$. Since there are no externally applied forces and moments to this system, the centre of mass should not undergo translation along the $x$ and $y$ axes during the flexural vibrations. However, due to the discretization errors introduced, the numerical solution will feature such motions of the centre of mass. Its location at each time instant, compared to point $(x,y)=(0,0)$ can be used as another indicator of the method’s effectiveness for this example. The centre of mass location can be calculated as

$$x_m^{h,n}=\frac{{\int }_a^b{x}^{h,n}|u_s^{h,n}|\mathit{d}s}{{\int }_a^b|u_s^{h,n}|\mathit{d}s},y_m^{h,n}=\frac{{\int }_a^b{y}^{h,n}|u_s^{h,n}|\mathit{d}s}{{\int }_a^b|u_s^{h,n}|\mathit{d}s}.$$

(27)

Figure 9 presents snapshots of the faulted ring during the first few unfolding-refolding cycles. The black line represents the modified method with ${\theta }_1=\theta =0.6$. The blue dashed line represents the S. Bartels method result. Due to the improved energy conservation properties of the modified method for parameter values close to 0.5, it can be seen that the folding-unfolding sequence is more accurately simulated by the modified method. This is more evident as the numerical simulation proceeds, where significant amount of energy is lost due to numerical dissipation. However, the effect is manifested even in the first cycle as indicated from the second row graphs in Figure 10.

Figure 9 Evolution of the faulted ring in example 4.3 at specific time instances during the first stages of the phenomenon for the SB method (blue dashed line) and the proposed modified method with ${\theta }_1=\theta =0.6$ (black line). The total number of elements is $\mathrm{N}=100$ and the time-step $\tau =0.05$.

Figure 10 Evolution of the faulted ring in example 4.3 during the first unfolding – refolding cycle for proposed modified method with ${\theta }_1=\theta =0.6$ (upper row) and the S. Bartels method (second row). The total number of elements is $\mathrm{N}=100$ and the time-step $\tau =0.05$.

The change in length of the beam and the translation of the centre of mass due to numerical errors are depicted in Figure 11. The upper row subplot depicts the change in length. The modified method performs excellent with respect to that, as the evolution is relatively smooth. Apart from some elongation during the first steps, the length of the beam does not change, verifying the prediction of Equation (24). The second subplot shows the horizontal location of the centre of mass, $x_m^{h,n}$. In reality this value should be equal to zero at for time instances due to the absence of external loads. However, due to numerical errors, the location is displaced towards the negative $x$ axis. This displacement manifests as an oscillation of the position with reducing amplitude. The modified method performs better in the sense that the $x_m^{h,n}$ value remains closer to zero. The vertical translation of the centre of mass $y_m^{h,n}$ is depicted in the third subplot of Figure 11. In this case the position oscillates around the zero coordinate assuming both positive and negative values. The amplitude of the oscillation is significantly lower for the modified method. In the case of the S. Bartels method, the amplitude of this spurious oscillation initially increases until approximately $t=55$ and then decreases due to numerical dissipation.

Figure 11 Faulted ring. Violation in the inextensibility constraint and spurious oscillation of the centre of mass around point $(x,y)=(0,0)$ for proposed modified method and S. Bartels method. The total number of elements is $\mathrm{N}=100$ and the time-step $\tau =0.05$.

The total simulation lasts for about twenty and a half cycles. Snapshots of the unfolding ring during the last stages of the simulation are depicted in Figure 12. The increased length in the case of the S. Bartels method, accompanied by a significant reduction in the vibration amplitude can be seen in the last snapshots when compared with the modified method.

Figure 12 Evolution of the faulted ring in example 4.3 at specific time instances during the last stages of the phenomenon for the SB method (blue dashed line) and the proposed modified method with ${\theta }_1=\theta =0.6$ (black line). The total number of elements is $\mathrm{N}=100$ and the time-step $\tau =0.05$.

Finally, the effect of parameter ${\theta }_1=\theta$ on the inextensibility constraint violation and the discrete energy conservation is studied in Figure 13. As expected, the energy (left hand side graph) conservation properties of the methods deteriorate as $\theta$ increases, and at the same time the violation in the inextensibility constraint decreases as well. The effect on the inextensibility constraint is much smaller though.

Figure 13 Effect of parameter ${\theta }_1=\theta$ on the inextensibility violation and discrete energy conservation properties of the modified method. Quantities calculated at $T=110$.

4.4 Swinging Rope (S. Bartels 2016)

In this last example, the model of the swinging rope presented in [15] will be solved using the proposed linearised $\theta$ scheme. The parameters selected for the domain are $I=[0,1]$, $d=3$ and the total duration of the phenomenon is $T=4$. All parameters are identical to the values used for the same example in [13] and no further nondimensional quantities are introduced. Hence, the flexural rigidity is set to $1/20$ and the gravitational force is $g={\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 40\hfill \end{array}\right]}^T$. Initial conditions for the position and velocity are

$$u(s,0)={\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill -s\hfill \end{array}\right]}^T,v(s,0)={\left[\begin{array}{ccc}\hfill 4s^2\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]}^T.$$

The upper end-point of the rope is clamped implying the boundary conditions $u(t,0)={\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]}^T$ and $u_s(t,0)={\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill -1\hfill \end{array}\right]}^T$ at $(x,y,z)=(0,0,0)$ for $t\in [0,T]$. The spatial discretization involves $40$elements, while two different numbers of time steps are considered, namely ${\mathrm{N}}_t=6400$ (same as in Bartels 2016) and ${\mathrm{N}}_t=9600$.

Figure 14 Evolution of the swinging rope in example 4.4 for the SB method and the proposed modified method with ${\theta }_1=\theta =0.75$. The total number of elements is $\mathrm{N}=40$ and time-steps are ${\mathrm{N}}_t=6400$. The modified method produces significant violations in the inextensibility constraint.

Figure 15 Evolution of the swinging rope in example 4.4 for the SB method and the proposed modified method with ${\theta }_1=\theta =0.75$ (black line). The total number of elements is $\mathrm{N}=40$ and time-steps are ${\mathrm{N}}_t=9600$. Virtually zero change in the beam length is produced by the modified scheme.

Figure 14 depicts the evolution of the swinging rope configuration and the change in length when the number of time steps is ${\mathrm{N}}_t=6400$. The modified method performs poorly in this case producing unacceptable violation in the inextensibility constraint. The S. Bartels method is very robust featuring good very small length changes. The same example is solved with $\mathrm{N}{}_t=9600$ and the results are depicted in Figure 15. In this case the modified method produces a virtually no changes in the swinging rope length and outperforms the S. Bartels method. This can be attributed to the very small time step selected and the effect of result in Equation (24) since due to the fine time discretisation the velocity profiles at successive time steps are rendered almost identical. It is worth noting that the S. Bartels method in this example is very robust even for less fine discretisation.

5 Conclusions

A linearized $\theta$ scheme for the vibrations of inextensible beams is introduced. The method is a modification of the S. Bartels technique [15] and has been shown to feature enhanced energy conservation characteristics and improved inextensibility attributes when the velocity evolution is relatively smooth. In certain cases the method benefits when initiated (first time step) by the S. Bartels algorithm. By appropriately selecting the $\theta$ parameters introduced, the energy conservation properties or the inextensibility attributes can be improved. The method becomes unstable for realistic selection of time steps if the free ${\theta }_1=\theta$ parameter approaches the value 1/2. Feasible choices of this parameter, to achieve acceptable balance between energy conservation and minimum inextensibility constraint violation, are in the range $(0.6,0.9)$. Due to its linearized character, the method is very fast and can be used as prediction step for fully nonlinear algorithms as well.

Acknowledgements

This work was supported by the Brunel Research Initiative & Enterprise Fund BRIEF 2018/2019-award number 1069.

References

[1] Maddocks JH. (1984) Stability of nonlinearly elastic rods. Archive for Rational Mechanics and Analysis, 85:311–54.

[2] Bosi F, Misseroni D, Dal Corso F, Bigoni D. (2014) An elastica arm scale. Proc. R. Soc. A 470,20140232, http://dx.doi.org/10.1098/rspa.2014.0232

[3] Bigoni D, Bosi F, Dal Corso F, Misseroni D. (2014) Instability of a penetrating blade. J. Mech. Phys.Solids 64, 411–425, http://dx.doi.org/10.1016/j.jmps.2013.12.008

[4] Xiao J, Chen X. (2013) Buckling morphology of an elastic beam between two parallel lateral constraints: implication for a snake crawling between walls. J. R. Soc. Interface 10, 20130399, https://doi.org/10.1098/rsif.2013.0399

[5] Dal Corso, F., Misseroni, D., Pugno, N.M., Movchan, A.B., Movchan, N.V., Bigoni, D. (2017) Serpentine locomotion through elastic energy release. J. R. Soc. Interface, 14, 20170055, doi: 10.1098/rsif.2017.0055

[6] Rus D, Tolley, MT. (2015) Design, fabrication and control of soft robots. Nature 521, 467–475, https://doi.org/10.1038/nature14543

[7] Laschi C, Cianchetti M, Mazzolai B, Margheri L, Follador M, Dario P. (2012) Soft robot arm inspired by the octopus. Adv. Robotics, 26, 709–727, https://doi.org/10.1163/156855312X626343

[8] Armanini, C., Dal Corso, F., Misseroni, D., Bigoni, D. (2017). From the elastica compass to the elastica catapult: an essay on the mechanics of soft robot arm. Proc. R. Soc. A, 473, 20160870, https://doi.org/10.1098/rspa.2016.0870

[9] Grothaus, M. & Marheineke, N. (2015) On a nonlinear partial differential algebraic system arising in technical textile industry: analysis and numerics. IMA J. Numer. Anal. 36, 4, 1783–1803, https://doi.org/10.1093/imanum/drv056

[10] Vlahopoulos N. and Bernitsas, M. M., (1988) Three dimensional nonlinear dynamics of nonintegral riser bundle, Journal of Ship Research, 35, 1, 40–57, ISSN 0022-4502.

[11] Vlahopoulos, N. and Bernitsas, M. M., (1988) Static Dynamic and Eigen Analysis of Non-Integral Production Risers, Applied Ocean Research, 10, 3, 144–154, https://doi.org/10.1016/S0141-1187(88)80014-2

[12] Liwen He, Jia Lou, Jianke Du, Huaping Wu (2018) Voltage-driven nonuniform axisymmetric torsion of a tubular dielectric elastomer actuator reinforced with one family of inextensible fibers. Eur J Mech A Solids, 71, 386–393, https://doi.org/10.1016/j.euromechsol.2018.06.004

[13] Caflisch, R. E. & Maddocks, J. H. (1984) Nonlinear dynamical theory of the elastica. Proc. Roy. Soc. Edinb. Sect. A, 99, 1–23, https://doi.org/10.1017/S0308210500025920

[14] Manning, R. S. (2014) A catalogue of stable equilibria of planar extensible or inextensible elasticords for all possible Dirichlet boundary conditions. J. Elasticity, 115, 105–130, https://doi.org/10.1007/s10659-013-9449-y

[15] Bartels, S. (2016) A simple scheme for the approximation of elastic vibrations of inextensible curves, IMA J. Numer. Anal, 36, 3, 1051–1071, https://doi.org/10.1093/imanum/drv054

[16] Bartels, S. (2013) A simple scheme for the approximation of the elastic flow of inextensible curves. IMA J. Numer. Anal., 33, 4, 1115–1125, https://doi.org/10.1093/imanum/drs041

[17] Lee P., Kim S., (2020) A variable-$\theta$ method for parabolic probems of nonsmooth data, Computers & Mathematics with applications, 79(4), 962–981, https://doi.org/10.1016/j.camwa.2019.08.006

[18] Bartels S., Reiter Ph., and Riege J. (2018). A simple scheme for the approximation of selfavoiding inextensible curves. IMA J. Numer. Anal., 38, 2, 543–565, https://doi.org/10.1093/imanum/drx021

Biography

Theodosios K. Papathanasiou is Lecturer (Structures Engineering) at the department of Civil and Environmental Engineering, Brunel University London. He specialises in Computational Mathematics/Mechanics and in particular the development of Finite Element schemes for Coupled-Field problems (e.g. Hydroelasticty, Thermoelasticity). Dr Papathanasiou is a Mechanical Engineer (National Technical University of Athens) and holds two MSc degrees (Computational Mechanics/Applied Mechanics) and a PhD from the NTUA school of Applied Mathematical and Physical Sciences. He has co-authored more than 30 publications in peer reviewed journals, more than 30 publications in international scientific conferences, and has participated in numerous funded research projects. Dr Papathanasiou was a Marie Curie Research Fellow (IAPP) at the University of Trento (Italy), where he conducted research on the thermomechanical response of refractory ceramics.

European Journal of Computational Mechanics, Vol. 30_1, 121–144.
doi: 10.13052/ejcm1779-7179.3015
© 2021 River Publishers

Abstract

1 Introduction

2 The Linearised $\mathit{\theta }$ Scheme

3 Properties of the $\mathit{\theta }$ Scheme

3.1 Energy Conservation Properties

3.2 Inextensibility Constraint Violation

4 Implementation and Numerical Results

4.1 Unwinding Helix 1 (S. Bartels 2016)

4.2 Unwinding Helix 2

4.3 Faulted Ring

4.4 Swinging Rope (S. Bartels 2016)

5 Conclusions

Acknowledgements

References

Biography