Simulating Space-Curved Beams using Structure-Preserving Integration Techniques

(1)

Simulating Space-Curved Beams using Structure-Preserving

Integration Techniques

June 2019

Master's thesis

Muhammad Hamza Khalid

2019Muhammad Hamza Khalid NTNU Norwegian University of Science and Technology Faculty of Information Technology and Electrical Engineering Department of Mathematical Sciences

(2)

(3)

Simulating Space-Curved Beams using Structure-Preserving Integration

Techniques

Muhammad Hamza Khalid

Masters in Mathematical Sciences Submission date: June 2019

Supervisor: Professor Elena Celledoni Co-supervisor: Dr. Nikita Kopylov

Norwegian University of Science and Technology Department of Mathematical Sciences

(4)

(5)

Abstract

Spatial deformations in beams subjected to terminal twist and axial compression is a sig- nificant engineering problem. When compressed more than a limit, the geometry becomes complex and rods form a loop with itself. This unstable behavior is essential to understand post-buckling in rods. The mathematical model used for the analysis of nonlinear deformations is a boundary value problem with a system of ordinary differential equations in space.

In the present work, we use geometric integration schemes to study the deformations in space–curved beams subjected to twist and vertical compression. We present a numerical setup for the solution of boundary value problem and simulate commonly observed phenomena such as loop formation, axial shortening, and self–contact. We also study the conservation properties in space, and numerically compare the preservation performance of structure preserving schemes with non–preserving integration schemes, for a twist–shortening problem. A multi–symplectic formulation for space–curved beams is also derived and numerically implemented using RATTLE and symplectic Euler method.

The theoretical and numerical setup presented here will be useful to understand relatively complex dynamic problems.

(6)

Preface

This dissertation is submitted to fulfill the requirements for the International Masters in Mathematical Sciences, at the Department of Mathematical Sciences, NTNU. The pro- gram code is MSMNFMA.

I want to express gratitude to my supervisor Professor Elena Celledoni, and Dr. Nikita Kopylov for their engagement and guidance throughout the process. I want to dedicate this thesis to my parents, without their love, support and prayers it would never have been possible. In the end, I want to thank my family, and friends for the good time spent in Norway.

Muhammad Hamza Khalid, June 15, 2019.

(7)

List of Tables

3.1 Butcher tableau . . . 24

3.2 Butcher tableau RK4 . . . 24

3.3 Butcher tableau GL4 . . . 24

3.4 Butcher tableau EM . . . 26

3.5 Butcher tableau IM . . . 26

4.1 Convergence Rates . . . 32

4.2 Cross-sectional properties of rod . . . 34

4.3 Convergence Rates . . . 34

4.4 Max–norm of the error in the preservation of conserved quantities for sta- ble parameters{φ, ξ}={π,0.4}and unstable parameters{φ, ξ}={0,0.8}. 37 4.5 Convergence Rates . . . 51

(10)

(11)

List of Figures

2.1 Coordinate system representation for a Free Rigid Body . . . 10

2.2 Twist–Shortening Explained (Miyazaki and Kondo, 1997) . . . 20

3.1 Projection Explained . . . 29

4.2 Convergence Free Rigid Body . . . 32

4.1 Angular momentum calculated using IM.T ol= 10⁻¹⁰,h=T /2¹³and T = 1000. . . 32

4.3 Error in the preservation of invariants. T ol = 10⁻¹⁰,h = T /2¹³ and T = 1000. . . 33

4.4 Discretization of rod . . . 33

4.5 Convergence for{φ, ξ}={π,0.5}. . . 34

4.6 Deformations in line of centroid for varyingφandξ= 0.4. . . 35

4.7 Deformations in line of centroid for varyingξandφ= 0. . . 36

4.8 Effect of varyingξandφon strains. . . 36

4.9 Error in preservation ofP^TP.N = 200,T ol= 10⁻¹⁰. . . 38

4.10 Error in preservation ofP^TM.N = 200,T ol= 10⁻¹⁰ . . . 38

4.11 Error in preservation ofJ.N = 200,T ol= 10⁻¹⁰ . . . 38

4.12 Max(abs(P^T_i Pi−P^T₀P0)),i= 1,2, ..N, using EM for differentξandφ. . 39

4.13 Max(abs(P^T_i Mi−P^T₀M0)),i= 1,2, ..N, using EM for differentξandφ. 39 4.14 Max(abs(Ji− J0)),i= 1,2, ..N, using EM for differentξandφ. . . 40

4.17 Values of max(abs(J_i− J₀)) for varyingξandφ. . . 40

4.15 Values of max(abs(P^T_iPi−P^T₀P0)) for varyingξandφ. . . 41

4.16 Values of max(abs(P^T_iMi−P^T₀M0)) for varyingξandφ. . . 42

4.18 E(qi)^TPiplotted on sphere forh= 0.01and{φ, ξ}={π,0.4}. . . 43

4.19 E(qi)^TPiforh= 1and{φ, ξ}={π,0.4}. . . 44

4.20 E(qi)^TP_iforh= 0.5and{φ, ξ}={π,0.4}. . . 45

4.21 E(qi)^TPiforh= 0.1and{φ, ξ}={π,0.4}. . . 46

4.22 E(qi)^TPiforh= 0.5and{φ, ξ}={0,0.8}. . . 47

4.23 E(qi)^TPiforh= 0.05and{φ, ξ}={0,0.8}. . . 48

(12)

4.24 Center of mass (E(qi)^Te₃) for{φ, ξ}={π,0.4}. . . 49

4.25 Convergence{φ, ξ}={π,0.5} . . . 51

4.26 Solution for twist–shortening using multi–symplectic formulation. . . 51

4.27 Solution for twist–shortening using multi–symplectic formulation. . . 52

4.28 Center of mass (E(qi)^Te3) for{φ, ξ} ={π,0.4}using multi–symplectic for RM and SE. . . 53

A.1 Error in the preservation of invariants. T ol = 10⁻¹³,h = T /2¹³ and T = 1000. . . 61

A.2 Maximum of Error in preservation for varyingξandφusing IM. . . 62

A.3 Maximum of Error in preservation for varyingξandφusing KM. . . 62

A.4 Maximum of Error in preservation for varyingξandφusing RK4. . . 62

A.5 Maximum of Error in preservation for varyingξandφusing GL4. . . 63

(13)

Abbreviations

R : Real numbers

R³ : 3-dimensional vector space of real numbers SO(3) : Space of3×3orthogonal rotation matrices

so(3) : Lie Algebra whose elements are3×3skew-symmetric matrices S³ : 3-sphere,S³={q= (q0,q)∈R×R³:q²₀+||q||²= 1}

s³ : Lie algebra of 3-sphere i.es³=TeS³={u= (0,u)∈ {0} ×R³ :u∈R³}

du

dt : u˙

du

ds : u⁰

P : Stress Resultant

M : Stress Couple

{γi}i∈I : Force Strains {κi}i∈I : Moment Strains

EM : Explicit Midpoint Method IM : Implicit Midpoint Method

KM : Kahan Method

SM : Spherical Midpoint Method ERK4 : 4th order Runge-Kutta Method

GL4 : 4th order collocation GaussLegendre Method SE : Symplectic Euler Method

RM : RATTLE Method

(14)

(15)

Chapter 1 Introduction

In the present work, we consider a three–dimensional elastic rod model formulated in Reissner (1981). The rod model discusses deformations in space–curved beams and is a static rod model; however, Simo (1985) presents an extension to the generalized dynamic rod model. This model has equal principal stiffness and can experience flexure, torsion, extension, and shear, in contrast to the Kirchhoff–Love model, (Love, 1944) that neglected extension and shear strain. Strain caused due to extension is an important aspect, when investigating the formation of knots and Kirchhoff–Love’s model is not feasible to analyze such phenomenon. However, with Reissner (1981) model, it is possible to set up a mathematical model that can experience all four strains. The resulting boundary value formulation can be used for various applications, for instance, twist-shortening problem or a contact problem (Miyazaki and Kondo, 1997). The former is a relatively simplified application, with boundary conditions on the start and end point of the rod, whereas the later has few additional conditions and assumptions which make it a relatively complex problem.

An analytical solution and possible applications for this model were presented in Miyazaki and Kondo (1997); however, the solutions are expressed in integral forms and require appropriate numerical integration techniques. Choice of coordinates for the configuration manifold is a crucial issue for the numerical simulation, and the first numerical discretization of this model presented by Simo and Vu-Quoc (1986) used rotation matrices to represent the configuration manifold. Miyazaki and Kondo (1997) choose Euler angles to represent rotations; however, we have opted to use unit quaternions (Euler parameters).

We use the configuration assumptions on the rod from Miyazaki and Kondo (1997) to write the equilibrium equations for space–curved beams in a compact form using vectors. The resulting model is a boundary value problem with a system of ordinary differential equations in space that can be numerically solved using geometric integrators. The problem was investigated for a twist–shortening problem discussed in subsection 2.2.4. Due to the nature of the problem, it is challenging to use integration techniques that preserve geometry, and also satisfy boundary conditions. To overcome this problem, we use the shooting method combined with schemes based on Runge–Kutta methods that treat a boundary

(16)

Chapter 1. Introduction

value problem as an initial value problem and uses root finding techniques for the solution. The advantage of the shooting method lies in the usability of time-marching schemes for initial value problems and will be discussed in detail in section 3.1.

One simple approach to the numerical discretization of Hamiltonian PDEs is to semi–

discretize them in such a way that the resulting ODE is a Hamiltonian system. This enables the use of symplectic methods or integral preserving methods in time. Another approach is to rewrite the Hamiltonian PDE into a multi–symplectic form (Bridges, 1997) and then provide a space–time discretization of the problem satisfying a discrete local conservation law of multi–symplecticity. As the problem is time–independent, the kinetic energy is zero, and the resulting Hamiltonian formulation has conjugate variables zeros. There- fore it is insignificant to investigate a Hamiltonian formulation. We adopt the work of Celledoni and S¨afstr¨om (2010), to write the multi–symplectic formulation for Reissner (1981) static model and discuss the conserved quantities in the problem. The simplified formulation presented here will be helpful to understand the multi–symplectic formula- tions for complex dynamical problems. Ringheim (2013) presented a framework for the multi–symplectic formulation on a static rod model, but no numerical experiments were considered. We do a successful implementation of the multi–symplectic formulation using symplectic integration techniques discussed in chapter 3. The goal of this thesis, therefore, is to numerically simulate deformations in space–curved beams when subjected to torsion and terminal thrust using geometric integration.

The sections below primarily focuses on the definitions of the mathematical terms and gives a brief insight into the requisite concepts for this thesis.

1.1 Lie Group and Lie Algebra

Definition 1.1.1(Lie Group). Anr-parameterLie group is a groupGcarrying a structure of an r-dimensional smooth manifold, such that the multiplication mapm:G×G→G,

m(g, h) =g.h, g, h∈G, and the inversion mapi:G→G,

i(g) =g⁻¹, g∈G, are smooth maps between manifolds (Olver, 2012).

Definition 1.1.2(Lie Bracket or Commutator). Letf: M →Rbe smooth function, and vandwvector fields onM, then their Lie Bracket[v,w](Olver, 2012) is the unique vector field satisfying

[v,w](f) =v(w(f))−w(v(f)).

Proposition 1. The Lie Bracket has the following properties, (Olver, 2012) (a) Bilinearity

[cv+c⁰v⁰,w] =c[v,w] +c⁰[v⁰,w], wherec, c⁰are constants.

(17)

1.2 Multi-Symplectic PDEs (b) Skew-Symmetry

[v,w] =−[w,v], (c) Jacobi Identity

[u,[v,w]] + [[w,u],v] + [v,[w,u]] = 0.

Definition 1.1.3(Lie Algebra). A Lie Algebra (Olver, 2012) over a fieldKis a vector spacegwith a bilinear map,[,] :g×g→g, such that the properties mentioned in Propo- sition 1 are fulfilled, i.e Lie Bracket is bilinear, skew-symmetric and Jacobi Identity is satisfied.

1.2 Multi-Symplectic PDEs

Definition 1.2.1(Multi-Symplectic PDE). If a PDE can be written in the form of a linear system of ordinary differentiable equations of type,

Mzt+Kzs=∇_zS(z), (1.1) whereM ∈R^d×dandK∈R^d×dare skew-symmetric matrices,z∈R^dandS:R^d→R is a smooth map, then the PDE is said to be multi-symplectic.

Define two-forms,ωandκ,

ω:=dz∧Mdz, κ:=dz∧Kdz, (1.2)

such thatω is associated with time direction and defines a symplectic structure onR^m (m = rankM), whereasκis associated with space direction and defines a symplectic structure onR^k (k=rankK)(Bridges and Reich, 2001). Any solution of the variational equation associated with Equation (1.1) satisfies the multi–symplectic conservation law,

∂tω+∂sκ= 0. (1.3)

Another important property associated with Equation (1.1) is the local conservation of energy and momentum, i.e

∂te(z) +∂sf(z) = 0, ∂ti(z) +∂sg(z) = 0, (1.4) where,

e(z) =S(z)−1

2z^T_sK^Tz, f(z) = 1 2z^T_tK^Tz, g(z) =S(z)−1

2z^T_tM^Tz, i(z) =1

2z^T_sM^Tz.

Global conservation quantities for energyE(z)and momentumI(z)can be obtained by the integration off(z)andi(z)over the spatial domain,

E(z) = Z L

0

e(z)ds, I(z) = Z L

0

i(z)ds, (1.5)

such that both are conserved, i.e(d/dt)E(z) = (d/dt)I(z) = 0.

(18)

1.3 First Integrals

Given a differential equation, z˙ = f(z(t)), the solution traces a curve in R^d and for a given function I(z) : R^d 7→ R the solution can be evaluated asI(z(t)). For some special functions,I(z(t))may be constant along the solutions and such special functions are known as first integrals.

Definition 1.3.1(First Integral). A non-constant functionI(z)is called afirst integral, if for any solutionz(t)of a linear system following holds. (Hairer et al., 2006),

I⁰(z)f(z(t)) = 0, =⇒ I(z(t)) =Const,

They are also known as invariants, and are crucial in geometric integration. The importance of first integrals varies for mechanical systems. They either help in defining the problem, with constant values having physical significance, or help in confining the solution on a bounded region. They also serve as a comparable measure for evaluating the performance of different numerical methods.

For a constant vectord,I(z) =d^Tzis calledlinear invariantifd^Tf(z) = 0for allz, and for a symmetricm×mmatrixC,I(z) =z^TCzis aquadratic invariant.

Theorem 1.3.1. For a systemm˙ = A(m)m, If the matrix A is skew-symmetric for all m, then the quadratic functionI(m) =m^Tmis an invariant (Marsden and Ratiu, 2013).

For example, in Equations (1.4) and (1.5) energy and momentum are invariant quantities that are preserved either locally or globally. Global conservation (1.5) is a consequence of local conservation (1.4), therefore, having local conservation is a stronger property. In case of time-independent problems, the conserved quantities hold locally in space (in fact point-wise), hence their integrals over the space domain also remain constant. Their importance for our problem, will be discussed in chapter 2.

1.4 Euler Angles

The orientation of a body in space can be specified using three successive rotations with respect to set of some coordinate axes fixed in the body, provided no two successive axes of the sequence are same. The specific sequence of rotations (Leimkuhler and Reich, 2004, page 228) used by Euler is as follows:

1. Counterclockwise rotationαabout the z-axis.

2. Counterclockwise rotationβabout the x-axis.

3. Counterclockwise rotationγabout the z-axis.

The rotation matrix∧can be expressed using the product of the following three planar rotations, i.e∧=ABCwhere,

A=





cosγ sinγ 0

−sinγ cosγ 0

0 0 1



, B=





1 0 0

0 cosβ sinβ 0 −sinβ cosβ



, C =





cosα sinα 0

−sinα cosα 0

0 0 1



,

(19)

1.5 Quaternions Elementary rotationsA, B, Cdo not commute, so there are 6 possible conventions for the representation of rotation matrix and the choice of axes for rotations is not unique.

Euler angles used for solving equations of motion of a rigid body do not cover all possible orientations, and singularities in the equations complicate the numerical integration of the resulting equations of motion, not only at the singularity but also in the neighbor- hood. One way to counter this problem is choosing a different sequence of axes (with new angle variables) whenever integration proceeds in the vicinity of a singular point, but this procedure will make the algorithm cumbersome. Because of this problem with the Euler angle, an alternative set of parameters based on Hamilton’s quaternions is often used (Leimkuhler and Reich, 2004).

1.5 Quaternions

A quaternion is quadruple of Euler parametersq = (q0,q)whereq0 is a scalar andq = (q₁, q₂, q₃)is a vector in three dimensional space. Quaternions (with unit length) are points on the unit three dimensional,S³described by the constraint,{q ∈R⁴: ||q|| = 1}. In other words,q∈R×R³is equipped with a certain Lie group structure and the product (non–commutative) “⊗” is defined as

(p₀,p)⊗(q₀,q) := (p₀q₀−p^Tq, p₀q+q₀p+p×q). (1.6) Product of quaternions can also be expressed as a matrix–vector product, i.e forp,q∈R⁴, p⊗q =R(q)p=L(p)q, whereL(p)andR(q)are orthogonal commutative matrices defined as

L(q) =

q0 −q^T q (q01+ ˆq)

andR(q) =

q0 −q^T q (q01−ˆq)

, where1∈R³is an identity matrix. Their sum is defined as

q+p= (q₀+p₀,q+p).

Forq6= (0,0), there exists an inverse,q⁻¹, q⁻¹= ¯q/||q||, ||q||=

q

q₀²+||q||²₂,

where the conjugate of quaternion is q¯ = (q0,−q)and the identity element is qid = (1,0). A unit quaternion with scalar as zero, i.eq = (0,q), is called apure quaternion.

Quaternions are often used to find rotation matrices to avoid dealing with singularities.

A unit quaternion q = (q0,q)can be transformed to rotation matrices using an Euler–

Rodrigues mapE:S³→SO(3)defined as

E(q) =I3+ 2q0bq+ 2bq², (1.7) where, “b” is the hat-map.

(20)

Definition 1.5.1(Hat-map). A mapb:R³→so(3)is called hat-map, if forv∈R³,

v=



 v1

v2

v3



→ bv=





0 −v3 v2

v3 0 −v1

−v2 v1 0



, and it satisfies,vub =v×u.

E(q)is surjective submersion, i.e its differential is surjective everywhere for a differential map between differential manifolds. However, it is not injective asE(q) =E(−q)that implies each rotation matrix will have two pre-images. For this reason, it is insignificant to consider either of one.E(q)is also a group homomorphism, sinceE(pq) =E(p)E(q).

Explicitly writing the Equation (1.7) gives the structure for the rotation matrix E(q) =





1−2(q²₂+q₃²) −2q0q₃+ 2q₁q₂ 2q₀q₂+ 2q₁q₃ 2q₀q₃+ 2q₁q₂ 1−2(q₁²+q₃²) −2q0q₁+ 2q₂q₃

−2q₀q₂+ 2q₁q₃ 2q₀q₁+ 2q₂q₃ 1−2(q₁²+q²₂)



.

Euler parameters are defined in terms of angle-axis parameter. Ifφis the angle of rotation around the axisk∈R³such that,||k||= 1, thenq= (cos(^φ₂),ksin(^φ₂))(Schwab, 2002;

Coutsias and Romero, 2004; Leimkuhler and Reich, 2004).

It is also possible to obtain Euler–parameters from the rotation matrix. This is done using the method presented by Shepperd (1978). As the map (1.7) is not injective, hence for each rotation matrix there are two possible choices of quaternions. The appropriate choice is decided using criterion such as continuity orq0>0. For the ease of evaluation we assume

z=





 z0

z1

z2

z3





 := 2





 q0

q1

q2

q3







, and T :=q11+q22+q33=Tr(E(q)),

where{qij}i,j=1,2,3 represents the(i, j)entry of the rotation matrix andq00 := T. On comparison with Euler–Rodrigues map we can write following symmetric equations for the diagonal entries ofE(q),

z₀²= 1 + 2q₀₀−T, z²₂= 1 + 2q₂₂−T,

z₁²= 1 + 2q11−T, z²₃= 1 + 2q33−T, (1.8) Similarly the off–diagonal entries gives the following equations

z0z1=q32−q23, z2z3=q32+q23, z0z2=q13−q31, z3z1=q13+q31, z₀z₃=q₂₁−q₁₂, z₁z₂=q₂₁+q₁₂.

(1.9)

The above information can be summarized into the following algorithm, that will be used to find Euler-parameters from a rotation matrix. (Shepperd, 1978; Egeland and Gravdahl, 2002).

(21)

1.5 Quaternions Algorithm

1. Computep_ii:=max{qii} 2. Compute|z_i|=√

1 + 2p_ii−T

3. Find the sign ofz_iusing some criteria, like continuity, orq₀>0.

4. Find remainingzjfrom the system of equations (1.8) and (1.9).

5. Computeq:= (q0,q)fromqi=zi/2fori∈ I.

(22)

(23)

Chapter 2 Space–Curved Beams

To understand a mechanical system that involves deformations in rods, an easy example to start with is a rigid body problem. Rigid body models, such as in Celledoni et al.

(2018, 2008); Celledoni and S¨afstr¨om (2010), were used to investigate applications like offshore pipe-lay operations, free rigid body and molecular dynamics, and they can serve as a good starting point to understand the mechanics of related problems. For convenience, we start with a free rigid body problem. We discuss the possible invariants and their importance and then write the Hamiltonian formulation using unit quaternions. We then write the mathematical model for the equilibrium equations of space–curved beams (Reissner, 1981; Miyazaki and Kondo, 1997) in a compact form using vectorial notation, and discuss the spatially conserved quantities. Using the mathematical model, we derive the multi–symplectic formulation of the rod model and then explain the idea behind twist–shortening. Time derivatives in this thesis will be represented by dot( ˙u)and spatial derivatives with prime(u⁰).

2.1 Free Rigid Body

2.1.1 Equations of Motion

Consider a rigid body (as shown in Figure 2.1) with center of mass fixed at origin in a con- vected coordinate system, with{n_i(s, τ)}_i∈I(I={1,2,3}) representing the orthonormal basis of the body frame, and{e_i(τ)}_i∈Ithe orthonormal basis of a stationary frame. Pa- rameters sandτ stands for space and time respectively. The configuration of the rigid body can then be determined using a rotation matrixQ∈SO(3), which transforms coordinate vectors from stationary frame to the body coordinates, in particular usingni =Qei

fori= 1,2,3. As is conventional, we denote stationary coordinate vectors by lower case letters, whereas body coordinate vectors by upper case letters.

IfM = (M₁, M₂, M₃)represents the angular momenta andI = diag(I_xx, I_yy, I_zz) the inertia tensor in body coordinates, where(Ixx, Iyy, Izz)are principal moments of inertia, then in the absence of any external force the equations of motion can be written as

(24)

Chapter 2. Space–Curved Beams

followed:

M˙ =I\⁻¹M M, (2.1)

Q˙ =QI\⁻¹M. (2.2)

Equation (2.1) is written for the angular momenta rather than angular velocityΩ=I⁻¹M

x-axis y-axis

z-axis

x-axis y-axis

z-axis

x’-axis y’-axis

z’-axis

Static Frame Body Frame

x(t)

O

Figure 2.1:Coordinate system representation for a Free Rigid Body

and is known as the Euler equation whereas Equation (2.2) is known as the Arnold equation. (Marsden and Ratiu, 2013; Hairer et al., 2006; Celledoni et al., 2008; Leimkuhler and Reich, 2004) provides more details on rigid bodies. In the next section, we formulate the equations of motion using unit quaternions.

2.1.2 Quaternionic representation of Equations of Motion

In this section, we will write the Equations (2.1) and (2.2) using quaternions, following (Egeland and Gravdahl, 2002). Before finding the quaternionic representation, we will write some helpful results for quaternions. Consider an arbitraryq∈S³. Differentiating q⊗q¯ =qidwith respect tot, we get

q˙ ⊗q¯+q⊗q˙¯=0, =⇒ q˙¯=−¯q⊗q˙ ⊗q¯. (2.3) Now considerq˙ ⊗q¯,

q˙0

˙ q

⊗ q0

−q

=

q˙0q0+ ˙q^Tq

−q˙0q+q0q˙+q×q˙

= 0

v

. (2.4)

wherev=−q˙0q+q0q˙+q×q, and we have used the following identity.˙ 1

2 d

dt(q²₀+q^Tq) = 0.

Now we have gathered the results and we can write a proposition that will give quaternionic representation for Equation (2.2).

(25)

2.1 Free Rigid Body Proposition 2. The quaternionic representation for Equation (2.2) is given by,

q˙ = 1

2ωe⊗q= 1

2q⊗Ω,e (2.5)

whereωe = (0, ω)andΩ = (0,e Ω), are unit quaternions representing angular velocity in stationary and body frame, respectively.

Proof. Consider an arbitrary vectoru= (0,u) ∈ s³ then the coordinate transformation using quaternion product is,

0 Qu

=q⊗ 0

u

⊗q¯, (2.6)

whereQ∈SO(3)is the rotation matrix. Differentiating with respect to time, 0

Qu˙

+ 0

Qu˙

= ˙q⊗ 0

u

⊗q¯+q⊗ 0

u˙

⊗q¯+q⊗ 0

u

⊗q˙¯, 0

Qu˙

= ˙q⊗ 0

u

⊗q¯+q⊗ 0

u

⊗q˙¯,

Using(2.6)foru˙ 0

Qu˙

= ˙q⊗q¯⊗q⊗ 0

u

⊗q¯−q⊗ 0

u

⊗q¯⊗q˙ ⊗q¯,

using(2.3) 0

Qu˙

= ( ˙q⊗q¯)⊗ 0

Qu

− 0

Qu

⊗( ˙q⊗q¯), 0

Qu˙

= 0

v

⊗ 0

Qu

− 0

Qu

⊗ 0

v

,

using(2.4) , 0

Qu˙

= 2 0

v×Qu

,

Comparing it with (2.2) gives,2v=ω. Using this we can write the quaternionic representation as followed

q˙ ⊗q¯ =1

2eω, =⇒ q˙ = 1

2ωe⊗q, (2.7)

whereeω= (0, ω). Using (2.7), we can write the following proposition.

Equation (2.5) can also be written using matrix–vector multiplication as following, q˙ =1

2R(q)ωe= 1

2L(q)eΩ. (2.8)

At each time-step, Equation (2.5) or (2.8) gives quaternion q ∈ S³, and using Euler–

Rodrigues map (1.7), we can find the rotation matrix. For a general overview of quaternion representation for kinematics equations see Egeland and Gravdahl (2002); Celledoni et al.

(2008); Schwab (2002); Coutsias and Romero (2004).

(26)

2.1.3 Quadratic Invariants

The two quadratic invariants for a free rigid body problem are as following, (Hairer et al., 2006; Marsden and Ratiu, 2013)

M^TM=M₁²+M₂²+M₃², (2.9)

H(M) =1

2hM, I⁻¹Mi= 1 2

M₁² I₁ +M₂²

I₂ +M₃² I₃

, (2.10)

The first quadratic invariant (2.9) confines the solution to a sphere, whereas the second quadratic invariant (2.10) confines it to the ellipsoid; therefore the solution is constrained to the intersection of sphere and ellipsoid. Preservation of either of these will guarantee the boundedness of solution ast→ ∞. These two quadratic invariants will later be used to distinguish the performance of different integrators in chapter 4.

2.1.4 Hamiltonian Formulation

Free rigid body problem can also be solved by writing the problem in terms of canonical coordinatesr = (p,q), using the Hamiltonian formulation. In this section we present a Hamiltonian formulation for the free rigid body problem in unit quaternions, using the idea presented in Maciejewski (1985). For angular velocityωe = (0, ω)∈s³, the kinetic energy (2.11) can be written as

L=1

2hω, I ωi=1

2 I₁ω₁²+I₂ω₂²+I₃ω²₃

, (2.11)

L= 1

2hω,e Ieωie =1

2heΩ,IeΩie = 2hq˙, L(q)I L(e q^c) ˙qi, where

Ie=

α 0^T 0 I

,

whereα6= 0is of no physical importance, but is chosen to ensure the regularity of solution.

Using Legendre transformation, the expression for the conjugate momenta can be written as,

p=∂L

∂q˙ = 4L(q)I L(e q^c) ˙q =⇒ q˙ = 1

4L(q)Ie⁻¹L(q^c)p. Therefore, the Hamiltonian formulation for a free rigid body can be written as,

p˙ =−1

4L(p)Ie⁻¹L(p^c)q, (2.12) q˙ = 1

4L(q)Ie⁻¹L(q^c)p. (2.13) In the next section, we will derive the equilibrium equations for space–curved beams and their multi–symplectic formulation in quaternions, and discuss the important conserved quantities, as we did for the free rigid body.

(27)

2.2 Space–Curved Beam

Spatial elastica is the mathematical model used to analyze deformations in space–curved beams. The static rod model under consideration was formulated by Reissner (1981), where the rod can experience flexure, torsion, shear, and extension and has equal principal stiffness. The model, however, was a static problem, and a generalized dynamic problem was later presented by Simo (1985). In this section, we write the equilibrium equations from Miyazaki and Kondo (1997) in a more compact form using vectors and state them as a boundary value problem. We then present a multi–symplectic formulation and discuss the conserved quantities for the model.

2.2.1 Equilibrium Equations in compact form using vectors

Consider a rod of lengthL, vertically clamped from both ends with the following features, (Miyazaki and Kondo, 1997)

• Equal principle stiffness and uniform cross-section.

• Direct proportionality of stress resultants and stress couples (irrespective of large displacements) on the force and moment strains, respectively.

• No deformations within the cross-section.

• No distributed load along the rod, and only terminal load.

The rod is made of a hyperelastic material, and there are no deformations within the cross- section, therefore its configuration can be fully described by position of its line of centroids by means of mapx: [0, L]7→R³, and the orientation of cross-section ats∈[0, L](Simo, 1985). Suppose{ni(s, τ)}i∈I(I ={1,2,3}) and{ei(τ)}i∈Irepresents the orthogonal basis vector for the body frame and the stationary frame, respectively where sdenotes space andτ time parameter. The origin of the axis is fixed at the centroid of the cross- section, and the equation for the centroid of rod is given by,

x⁰ =γ₁^bn₁+γ₂^bn₂+ (1 +γ₃^b)n₃, (2.14) where{γ_i^b}_i∈Irepresents force strains withγ₁^bandγ₂^bmeasuring shear aboutn1andn2, while γ₃^b measures extension. For convenience, “b” and “r” in superscript will denote corresponding terms in body frame and stationary frame, respectively. If{Pi(s)}_i∈I is the stress resultant and{Mi(s)}i∈Istress couple in body coordinates, then the equilibrium equations can be written as (Reissner, 1981)

P₁⁰−P₂κ^b₃+P₃κ^b₂= 0, M₁⁰ −M₂κ^b₂+M₃κ^b₂−P₂(1 +γ₃^b) +P₃γ₂^b= 0, P₂⁰−P3κ^b₁+P1κ^b₃= 0, M₂⁰ −M3κ^b₁+M1κ^b₃−P3γ₁^b+P1(1 +γ₃^b) = 0, P₃⁰−P1κ^b₂+P2κ^b₁= 0, M₃⁰−M1κ^b₂+M2κ^b₁−P1γ^b₂+P2γ^b₁= 0,

(2.15)

where{κ^b_i}_i∈Iare the moment strains, for torsion and bending. Using the assumption of linear dependence of stress and strain i.ePi∝γ_i^bandMi∝κ^b_i, expressions forκ^b, γ^bcan be calculated.

(28)

P1=Kγ₁^b, P2=Kγ₂^b, P3=K3γ₃^b, =⇒



 γ^b₁ γ^b₂ γ^b₃



=





1

K 0 0

0 _K¹ 0 0 0 _K¹

3







 P1

P2

P3





M1=Aκ^b₁, M2=Aκ^b₂, M3=A3κ^b₃, =⇒



 κ^b₁ κ^b₂ κ^b₃



=





1

A 0 0

0 _A¹ 0 0 0 _A¹

3







 M1

M2

M3





(2.16)

whereA, A₃, K, K₃are constants representing flexural, torsional, shear, and extensional rigidity respectively. IfP= [P₁, P₂, P₃]^T andM = [M₁, M₂, M₃]^T, then (2.15) can be written in compact form as

P⁰=−C\_M⁻¹MP,

M⁰=−C\_M⁻¹MM−C\_P⁻¹PP−eb3P,

(2.17) where

C_M =diag[A, A, A₃], C_P =diag[K, K, K₃].

Following Simo (1985), we take the derivative of the body frame with respect to space, to find the equation for the rotation matrix.

n⁰_i=k×ni, where k=κ^b₁n1+κ^b₂n2+κ^b₃n3=Qκ^b =QC_M⁻¹M,

whereQ = [n1, n2, n3] ∈SO(3)is the orthogonal rotation matrix, that will specify the orientation of rod. Similarly, we can write the Equation (2.14) as followed

Q⁰ =QC\_M⁻¹M. (2.18)

Similarly, the equation for the position vector can be re-written in compact form as follows, x⁰=QC_P⁻¹P+Qe₃=QC_P⁻¹P+n3. (2.19) Combining Equations (2.17), (2.18) and (2.19), results in a boundary value problem that can be used to study deformations in space–curved beams. Now we will discuss the Hamil- tonian structure and associated energy for the space–curved beams that will be helpful in the derivation of multi–symplectic formulation in subsection 2.2.3.

Following Simo et al. (1995, 1988), Equations (2.17) can also be calculated by differentiating the free-energy function with respect to the strain measures. If the free-energy is a quadratic function given by the strain measures, then we can write the expression for energy density in space–curved beams as

J(γ^r, κ^r) :=1 2

hγ^r, DPγ^ri+hκ^r, DMκ^ri

, (2.20)

where,

DP =QCPQ^T, DM =QCMQ^T, γ^r=x⁰−Qe3,

(29)

2.2 Space–Curved Beam and,

γ^r=Qγ^b, κ^r=Qκ^b, p=QP, m=QM. (2.21) Differentiating potential energy (2.20) with respect toγ^randκ^r, we get

p= ∂

∂γ^rJ =DPγ^r =⇒ γ^r=D⁻¹_P p, m= ∂

∂κ^rJ =D_Mκ^r.

Using (2.21), the above expressions can be written in body coordinates as followed P=Q^TDP(Qγ^b) =CPγ^b, whereγ^b=Q^Tx⁰−e3,

M=Q^TDM(Qκ^b) =CMκ^b.

The general form for the balance of linear and angular momentum in space are given by (Reissner, 1981)

∂_sp+ep= 0, (2.22)

∂sm+ (∂sx)×p+me = 0, (2.23)

whereep,m, are the external forces ande xis the position vector for line of centroid. For derivations of (2.22) and (2.23), see Green and Zerna (1992). The Hamiltonian, i.e the total energy of the system, is given by the sum of kinetic energy and the potential energy, but since the kinetic energy is zero, therefore its Hamiltonian is

U =1 2

Z L 0

hγ^r, DPγ^ri+hκ^r, DMκ^rids. (2.24) Proposition 3. With no external forces acting on the system equilibrium equations (2.22) and (2.23) gives a system equivalent to system (2.17).

Proof. Consider Equations (2.22) and (2.23) in body frame withep=me = 0, (2.22) =⇒ p⁰ = (QP)⁰=Q⁰P+QP⁰,

QP⁰ =−QC\_M⁻¹MP, =⇒ P⁰ =−C\_M⁻¹MP.

(2.23) =⇒ m⁰ =−x⁰×p,

(QM)⁰ =−(QC_P⁻¹P+Qe3)×QP, Q⁰M+QM⁰ =−QC_P⁻¹P×QP−Qe3×QP,

QM⁰ =−QC\_M⁻¹MM−(QC_P⁻¹P)×(QP)−Qe3×QP, QM⁰ =−QC\_M⁻¹MM+ (QP)×(QC_P⁻¹P)−Qe3×QP,

M⁰ =−C\_M⁻¹MM+P×(C_P⁻¹P)−e3×P, M⁰ =−C\_M⁻¹MM−C\_P⁻¹PP−be3P.

Hence, we conclude that both systems give identical equations.

(30)

2.2.2 Conserved Quantities

As seen in the time–dependent problem, described in subsection 2.1.3, preservation of quadratic first integrals such as energy, momentum, over time is important. Likewise, we consider here scalar functions of solution which are invariant with respect to space variables, i.e the derivative with respect to spacesis zero. We refer to them as conserved quantities. The first conserved quantity for the problem, as remarked in Miyazaki and Kondo (1997) is

κ^b₁P1+κ^b₂P2+A3

A κ^b₃P3=Const, 1

AM₁P₁+ 1

AM₂P₂+ A₃

A3AM₃P₃=Const,

using Equation(2.16)

=⇒ 1

A(M1P1+M2P2+M3P3) =Const.

Similarly, the other conserved quantities are the applied terminal force (defined as the integral of stress resultants over spatial domain) and curvatureκ3. Using compact vector notation the conserved quantities are proved in the propositions below.

Proposition 4. P^TPand P^TMare the two conserved quantities for the boundary value problem of space–curved beams.

Proof. Consider the first conservedP^TP, and differentiate with respect tos, P^0TP+P^TP⁰ = (−C\_M⁻¹MP)^TP+P^T(−C\_M⁻¹MP),

= (P×C_M⁻¹M)^TP+P^T(P×C_M⁻¹M).

Using property of scalar triple productu^T(v×w) =v^T(w×u) =w^T(u×v)and commutative propertyu^Tv=v^Tu,

P^0TP+P^TP⁰=C_M⁻¹M(P×P) +C_M⁻¹M(P×P) = 0,

=⇒ P^TP=Const.

Now, consider the second conserved quantity,P^TMand differentiate again with respect to s

P^0TM+P^TM⁰= (−C\_M⁻¹MP)^TM+P^T(−C\_M⁻¹MM−C\_P⁻¹PP−ce₃P),

= (P×C_M⁻¹M)^TM+P^T(M×C_M⁻¹M+P×C_P⁻¹P−e₃×P),

= (P×C_M⁻¹M)^TM+P^T(M×C_M⁻¹M) +P^T(P×C_P⁻¹P)

| {z }

=0

−P^T(e3×P)

| {z }

=0

,

=M^T(P×C_M⁻¹M) +P^T(M×C_M⁻¹M),

Using property of scalar triple product we can write, M^T(P×C_M⁻¹M) = −P^T(M× C_M⁻¹M),

P^0TM+P^TM⁰=−P^T(M×C_M⁻¹M) +P^T(M×C_M⁻¹M) = 0,

=⇒ P^TM=Const.

(31)

2.2 Space–Curved Beam To verify local the conservation of energy (1.4) over the spatial domain, we re-write equations (2.17) in the formy⁰=B(y)∇_yJ(y),

P⁰ M⁰

= 0 ˆP

Pˆ Mˆ

C_P⁻¹P+e₃ C_M⁻¹M

, (2.25)

whereB(y)is a skew-symmetric matrix andJ(y)energy (2.20) in body coordinates. i.e, J = 1

2hP, C_P⁻¹P+ 2e3i+1

2hM, C_M⁻¹Mi. (2.26) Proposition 5. Energy density (2.26) is constant over the spatial domain.

Proof. Differentiating Equation (2.26) with respect to space “s”.

J⁰(y) = (C_P⁻¹P)^TP⁰+e^T₃P⁰+ (C_M⁻¹M)^TM⁰,

= (C_P⁻¹P)^T(P×C_M⁻¹M) +e^T₃(P×C_M⁻¹M)

+ (C_M⁻¹M)^T(M×C_M⁻¹M+P×C_P⁻¹P+P×e₃),

= (C_P⁻¹P)^T(P×C_M⁻¹M) +e^T₃(P×C_M⁻¹M)

+ (C_M⁻¹M)^T(P×C_P⁻¹P) + (C_M⁻¹M)^T(P×e3),

=−(C_M⁻¹M)^T(P×C_P⁻¹P) +e^T₃(P×C_M⁻¹M)

+ (C_M⁻¹M)^T(P×C_P⁻¹P)−e^T₃(P×C_M⁻¹M),

= 0.

FromJ⁰(y) = 0, we conclude (2.26) is constant over spatial domain.

The idea of conservation of energy is followed from local conservation laws (1.4) and (1.5) for multi–symplectic PDEs discussed in section 1.2. However, since the problem is time–independent, the energy density differentiated with respect to time in (1.4) disap- pears. We are then left with a derivative of energy density due to strain, i.eJ with respect tos, and in Proposition 5 we have shown it is equal to zero, therefore we say energy is locally (in fact point-wise) conserved. Since having local conservation is a stronger property and global conservation is its consequence, therefore we sayJ is conserved globally as well, i.e (1.5) is satisfied.

2.2.3 Multi–Symplectic Formulation

Using the idea presented in subsection 2.1.2, the equilibrium equation (2.18) for the position matrix can also be written in unit quaternions as followed

q⁰ =q⊗eκ^r=eκ^b⊗q,

forκe^r= (0, κ^r)∈s³andeκ^b= (0, κ^b)∈s³. Or using matrix–vector multiplication, q⁰= 1

2L(q)eκ^b= 1

2R(q)κe^r.

(32)

Similarly, Equation (2.19) in unit quaternions is as followed x⁰=E(q)C_P⁻¹P+n3,

whereE(q)is the rotation matrix calculated using Euler–Rodrigues map (1.7). As discussed in chapter 1, it is insignificant to consider Hamiltonian formulation, because of zero conjugate variables. We therefore derive a multi–symplectic formulation for space–

curved beams, that will be applied to twist–shortening problem in chapter 4.

Following Celledoni and S¨afstr¨om (2010), consider a smooth mapS defined onu= (x,q)for conjugate variablesv= (vx,v)defined as

S(u,v) :=hv,us(v)i − L(u,us), (2.27) whereLis the Lagrangian. For convenience,CM ∈R^3×3is extended to invertibleCeM ∈ R^4×4so that the new Lagrangian becomes regular on the cotangent bundle.

DeM =L(q)R(q^c)CeML(q^c)R(q), whereCeM =

α 0^T 0 CM

, α6= 0.

The Lagrangian in general is given by the difference of the potential energy and the kinetic energy, but since the kinetic energy is zero, we can use the potential energy (2.20) to write the Lagrangian in quaternions with a holonomic constraintg(q) :=||q||²−1 = 0.

L(u,us) =−1 2

hhγ^r, DPγ^ri+ 4hq⁰, R(q)DeMR(q^c)q⁰ii

−λ

||q||²−1

, (2.28) whereγ^r = x⁰ −n₃. Analogous to the canonical system written for a free rigid body, we can write a similar system for conjugate variables by taking partial derivatives of the Lagrangian with respect to the space derivatives

vx:= ∂L

∂x⁰ =−DPγ^r, =⇒ x⁰(q,vx) =−D_P⁻¹vx+(q)e3, (2.29) v:= ∂L

∂q⁰ =−4R(q)DeMR(q^c)q⁰ ∈T^∗S³, =⇒ q⁰(q,v) =−1

4R(q)De⁻¹_MR(q^c)v, (2.30) Using these conjugate variables, we rewrite the Equations (2.28) and (2.27),

(2.28) =⇒ L(u,us) =−1 2

hhγ^r, DPγ^ri+ 4hq⁰, R(q)DeMR(q^c)q⁰ii

−λ

||q||²−1 , L(u,us) =−1

2

hh−D_P⁻¹vx, DP(−D_P⁻¹vx)i+1

4hv, R(q)De⁻¹_M R(q^c)vii

−λ

||q||²−1 , L(u,us) =−1

2

hhvx, D⁻¹_P vxi+1

4hv, R(q)De⁻¹_MR(q^c)vii

−λ

||q||²−1 . Calculatinghv,u_s(v)iin (2.27),

hv,us(v)i=hvx,x⁰i+hv,q⁰i,

=hvx,−D_P⁻¹vx+(q)e3i+1

4hv, R(q)De_M⁻¹R(q^c)vi,

=−hvx, D_P⁻¹vx−(q)e3i+1

4hv, R(q)De_M⁻¹R(q^c)vi.

(33)

2.2 Space–Curved Beam (2.27) =⇒ S(u,v) =hv,u_s(v)i − L(u,u_s),

=−hvx, D_P⁻¹vx−(q)e₃i+1

4hv, R(q)De_M⁻¹R(q^c)vi +1

2

hhvx, D⁻¹_P vx)i+1

4hv, R(q)De⁻¹_MR(q^c)vii +λ

||q||²−1 ,

=−1 2

hhv_x, D_P⁻¹v_x−2(q)e₃i+1

||q||²−1 . Summarizing (2.27) as

S(u,v) =−1 2

hhvx, D⁻¹_P vx−2(q)e3i+1

||q||²−1 . (2.31) The partial derivatives ofSwith respect toxandqare given as

∂S

∂x = 0, (2.32)

and

∂S

∂q =1

4R(q)L(v)L(q^c)De_M⁻¹R(q^c)v+R(q)

L(vx)−R(vx)

0 D⁻¹_P v_x−(q)e₃

+ 2h((q)−1)vx, D⁻¹_P vx−(q)e3i+ 2λq. (2.33) We can then write equations for the multi–symplectic system (1.1) as follows (Celledoni and S¨afstr¨om, 2010):

∂S

∂u =−∂sv, (2.34)

∂S

∂v =∂_su, (2.35)

0 =||q||²−1, (2.36)

Or in matrix–vector notation,





0 −1 0 1 0 0

0 0 0







 u⁰ v⁰ λ⁰



=





∂uS

∂vS

∂λS



 (2.37)

Hence, forz= (u,v, λ), and a skew-symmetric matrixK, the multi–symplectic formulation can be written as

Kz_s=∇zS, (2.38)

The multi–symplectic formulation calculated in this section corresponds to the formulation for a dynamic problem presented in Celledoni and S¨afstr¨om (2010) with time variables

(34)

zero. Constraintg(q) = 0ensures the solution stays on the surface of sphere, but since v∈T^∗S³we need another constraint to ensurevstays on the co-tangent bundle. To find the hidden constraint, differentiatingg(q) = 0once, we get

g⁰(q) = 2hq,q⁰i= 2hq,−1

4R(q)De_M⁻¹R(q^c)vi=−α

2hqid, R(q^c)vi=−α 2hq,vi,

=⇒ g⁰(q) =hq,vi= 0.

ForG(q) =∂g/∂q, this can be re-written as

g⁰(q) =G(q)v= 0, (2.39)

The second constraint will ensure v stays on the co-tangent bundle. We can find the expression for the Lagrange multiplier, by differentiating it again

g⁰⁰(q) =G⁰(q)v+G(q)v⁰

=hq⁰,vi+hq,v⁰i

=λ+h((q)−1)vx, D⁻¹_P vx−(q)e3i= 0,

=⇒ λ=−h((q)−1)vx, D_P⁻¹vx−(q)e3i. (2.40)

2.2.4 Twist–Shortening Problem

ɸ

m

l

e

e e

3

1

ξl

Figure 2.2:Twist–Shortening Explained (Miyazaki and Kondo, 1997)

As discussed in the chapter 1, when a vertically clamped rod is subjected to axial and moment stresses, it can experience deformation in its shape, depending on the configuration of the rod. Suppose such a rod of lengthl is twisted from one end by an angleφin the positive z-direction and then vertically compressed down by an amountξl, then it will

(35)

2.2 Space–Curved Beam experience deformation in its shape, as shown in Figure 2.2. The boundary conditions for this problem will then be

Q0=





1 0 0 0 1 0 0 0 1



, Ql=





cos(φ) −sin(φ) 0 sin(φ) cos(φ) 0

0 0 1



,

x₀=



 0 0 0



,x_l=



 0 0 (1−ξ)l



.

(2.41)

The conditions have been represented using Euler angles but we will use algorithm (1.5) with conditionq0>0, to find the corresponding Euler parameters in the numerical imple- mentations chapter 4.

(36)

(37)

Chapter 3 Structure Preserving Geometric Schemes

In the previous chapter, we derived explained the mechanical systems that are considered for the present study. In this chapter, we highlight the essential properties of geometric integrators that are suitable for initial value problems and posses additional geometric structures, e.g preservation of first integrals, energy, and symplecticity. We will later use these methods in combination with shooting techniques (see chapter 4) to solve boundary value problems. Consider a non-autonomous differential equation,

˙

z(t) =f(t, z(t)), z(t₀) =z₀. (3.1) The solution is approximated at equally spaced specified points using numerical scheme, suitable for the differential equation. For the present work, we investigate structure preserving geometric integrators based on the class of Runge–Kutta methods. They are used widely for the integration of equations of type (3.1), however many of them do not ex- hibit properties necessary for the solution. In the next section, we present the numerical schemes for initial value problems that were used in the numerical experiments.

3.1 Initial Value Problem

3.1.1 Runge–Kutta Methods

For an ordinary differential equation (3.1), Runge–Kutta methods take the form Ki=f(tn+cih, zn+h

i−1

X

j=1

aijKj), for i= 1,2, . . . , s

=⇒zn+1=zn+h

s

X

i=1

biKi,

Simulating Space-Curved Beams using Structure-Preserving Integration Techniques

Simulating Space-Curved Beams using Structure-Preserving

Integration Techniques

Master's thesis

Muhammad Hamza Khalid

Simulating Space-Curved Beams using Structure-Preserving Integration

Techniques

Muhammad Hamza Khalid

Abstract

Preface

Table of Contents

List of Tables

List of Figures

Abbreviations

Chapter 1

Introduction

1.1 Lie Group and Lie Algebra

1.2 Multi-Symplectic PDEs

1.3 First Integrals

1.4 Euler Angles

1.5 Quaternions

Chapter 2

Space–Curved Beams

2.1 Free Rigid Body

2.1.1 Equations of Motion

2.1.2 Quaternionic representation of Equations of Motion

2.1.3 Quadratic Invariants

2.1.4 Hamiltonian Formulation

2.2 Space–Curved Beam

2.2.1 Equilibrium Equations in compact form using vectors

2.2.2 Conserved Quantities

2.2.3 Multi–Symplectic Formulation

2.2.4 Twist–Shortening Problem

Chapter 3

Structure Preserving Geometric Schemes

3.1 Initial Value Problem

3.1.1 Runge–Kutta Methods