A consistent dynamic finite element formulation for a pipe using Euler parameters

A pipe element developed earlier for the analysis of combined large bending and torsional deformations of blood vessels under static loading is extended to model behavior in the presence of large rotations and time-varying external forces. As in the case of the earlier element, the enhanced element supports ovalization and warping of its cross-section. The enhancements presented in this paper are comprised of a mass matrix and gyroscopic effects resulting from fast rotation rates and large deformations. The effectiveness of the element is demonstrated by two examples, which simulate the three-dimensional behavior of a highly flexible pipe under dynamic loading conditions.


Introduction
A straight pipe element that enables the efficient computation of large, three-dimensional deformations in blood vessels with initially circular cross-sections was presented in an earlier study by Jiang and Arabyan [5].The important features of that formulation are that it can account for both rigid-body and constant strain modes, and captures all stresses except the normal stress across the wall thickness (i.e., thin shell is assumed).Euler parameters are used to describe rotational rigid-body modes and are incorporated into the element's vector of degrees of freedom.Under general loading (axial, transverse, bending and torsion), the element (called the "RC element" in Jiang and Arabyan [5]) allows large ovalization and warping of its cross section and large, three-dimensional, angular changes in the orientation of its reference axis.The formulation used to derive the element incorporates the nonlinear coupling between torsional and bending deformations.The RC element differs from previous formulations by Bathe and Almeida [3], and Ohtsubo [6] in that it captures both constant strain and rigid-body modes, and can be applied to any cases involving geometrically nonlinear large deformation problems.
In [5] the variation in the performance of the RC element with grid refinement was shown in comparison to a pipe structure made up of four-node shell elements.It was shown that the stress resultants computed by the RC element converge to "exact" results (solutions obtained using very large number of four-node shell elements) with a smaller number of degrees of freedom than those obtained using four-node shell elements in ABAQUS [1].These results were shown using static loading examples.
In this paper, the RC element formulation is extended to include dynamic loads via a consistent formulation.In linear analysis, the simplest way to obtain the vector of inertia nodal forces is to take the inner product of a diagonal mass matrix with the nodal acceleration vector.This formulation and numerical procedure have been proven to be very effective in most applications.However, in the case of this pipe element, due to the special interpolation functions used, relations between the field variables and nodal displacements are nonlinear.Therefore the method used in linear analysis cannot account for all inertia effects.In order to capture all inertia effects and coupling among extension, torsion, and bending deformations, the dynamic formulation of the pipe element must be derived from the full geometrically non-linear shell theory by a consistent formulation.
This paper presents a consistent formulation for the dynamic analysis of the three-dimensional pipe elements in which the inertia forces due to large deformations are fully accounted by capturing all the inertia couplings.The assumption of a thin pipe wall is carried through from previous work.The nodal co-ordinates, Shock and Vibration 5 (1998) 111-117 ISSN 1070-9622 / $8.00 © 1998, IOS Press.All rights reserved incremental displacements and rotations are defined in terms of a fixed global co-ordinate system.

Review of the formulation for static analysis
Here, some basic concepts and equations that were developed in [5] are reviewed, these being essential for the dynamic formulation that follows.

Basic assumption
The following assumptions are made in the derivation of the non-linear behavior: 1.The pipe is assumed to be thin-walled, i.e., if h is the wall thickness, and R is the radius of the undeformed pipe, then h/R 1.As a result the behavior of the pipe can be described by a single surface which is chosen as the mid-surface (i.e., the surface exactly between the outer and inner surfaces of the pipe).2. The material is assumed to remain linearly elastic throughout the range of deformation of the pipe.3. Normals to the mid-surface remain straight after deformation.4. Stresses along the pipe thickness direction is negligible.

Co-ordinate systems
To describe the motion of the element, three reference frames (see Fig. 1) are defined: 1.A fixed global reference frame XY Z.
2. An element reference frame xyz in which the initial angular positions are defined.3. A floating reference frame moving with the undeformed pipe cross section.Ovalization and warping are measured relative to this frame.The angular position of this frame relative to xyz is described by a rotation vector φ.
In addition, two curvilinear coordinates, s and t, are defined such that t is a measure of length on the midsurface along a direction parallel to the pipe's longitudinal axis and s is a measure of length on the midsurface in a direction perpendicular to t.

Element interpolation functions
The general displacement of the pipe element is assumed to consist of three parts (see Fig. 1): 1.The displacement of a point on the pipe axis, denoted by the vector v which is a function of t only.Let a be the position vector defining a point on the pipe axis in the reference configuration.Then the position of a point on the mid-surface of the deformed configuration is

x(s, t) = a(t) + v(t) + B(t)(r + r(s, t)).
( The director field t 2 of the deformed configuration is defined in terms of Euler parameters e, as: where T 2 (s) is the unit normal of the mid-surface in the undeformed configuration, A(e, t) is the rotation matrix, and e is the vector of four Euler parameters which are functions of the rotational vector φ defined in [5].Note that A(e, t) is defined independently from B in order to include transverse shears strain.

Generalized elasticity forces
The virtual work method is employed to obtain the generalized elasticity forces.The virtual work done by internal stress resultants is given by where N , M , and Q are the stress resultants as described by Cook and Malkus [4] shown in Fig. 2, and ε, κ, γ are membrane strains, bending strains, and transverse shear strains, respectively, as defined in [5].Using this formulation, the strain variations are where U T = {q T p T }, q is the vector of translational degrees of freedom and p is the vector of rotational degrees of freedom.B ε , B κ , and B γ are the material constitutive matrices.By substituting the above relations into Eq.( 3), the virtual work of the elasticity force can be obtained as: where Ψ e is the generalized elastic force defined as

Dynamic equations
The virtual work of the inertia force can be expressed as where δx is the variation of the position vector of a point on the deformed configuration of the midsurface, referred to an inertial frame of reference.The quantities δx and its time derivatives are expressed explicitly as: and where are the Lagrangian interpolation functions defined in [5], and R Bi (i = 1, 2 for two axial-node element) is defined by Substituting Eqs (8) and (10) into Eq.( 13), one obtains the variational expression for the virtual work done by the inertia forces: where Ψ d is the generalized inertia force, and In Eq. ( 14) C is a damping matrix and M is a consistent mass matrix.
In the development of the inertia force, the entire mass of the pipe is assumed to be concentrated on the mid-surface.Normally, this approximation leads to adequate precision for thin walled pipes.

Equations of motion
Now the elastic force, the inertia force, and external force T can be added together to obtain the equations of motion: Using the mass and damping coefficient matrices of Eq. ( 14), we obtain:

Numerical method
The equilibrium relation in Eq. ( 16) must be satisfied throughout the history of load application.Due to its nonlinearity, Eq. ( 16) must be solved iteratively for U(τ ) for any time τ .
The basic approach described by Bathe [2], in an incremental step-by-step solution is to assume that the solution U at time τ is known.In other words it is assumed that is satisfied.The solution for the discrete time τ + ∆τ is required, where ∆τ is a suitably chosen time increment.Hence, considering Eq. ( 16) at time τ + ∆τ we want: where the left superscript denotes "at time τ + ∆τ ".Since the nodal forces τ +∆τ Ψ e + τ +∆τ Ψ d depend nonlinearly on nodal point displacements, τ +∆τ U must be found by an iterative scheme.In this work, the generalized Newton-Raphson method is used for this iteration.The relations used in the Newton iteration are where τ +∆τ K (i−1) is the tangential stiffness matrix in iteration i−1.This matrix is assumed to be nonsingular and is given by Bathe [2].
where α is an integration parameter (see below).
The basic iteration formulae (19) provides a solution for nodal displacements for dynamic equilibrium at τ + ∆τ .This solution must then be used to compute nodal velocities and accelerations.The Newmark method is used for that purpose, using the general formula in Bathe [2]: where α and δ are parameters that are chosen to obtain integration accuracy and stability.When δ = 1/2 and α = 1/6, relations (21) correspond to the linear acceleration method (also called trapezoidal rule).In this work, these parameters were picked as δ = 1/2 and α = 1/4, making the integration unconditionally stable with O(∆τ ) accuracy.The relations employed in the integration scheme are It is evident that the iterative scheme used in dynamic analysis is the same as that used in Jiang and Arabyan [5], except that the tangential stiffness matrix K now contains contributions from the inertia of the system.
All the numerical results presented in the following sections have been obtained using this numerical scheme.For most dynamic analysis examples, the time step sizes were of the order of 10 −3 -10 −2 s.

Examples
Two examples are given below to demonstrate the effectiveness of the formulation described above in dynamic analysis.

Flexible pipe with end rotation
In this example, a pipe rotating vertically about a horizontally axis passing through one end is considered.Figure 3    The time step size was set to 0.001 and the relative error tolerance was set to 10 −3 .The deformation of the pipe at several points in time are shown in Fig. 4.

Fast rotating pipes with end lateral force
The second example is a pipe rotating about a fixed horizontal axis passing through one end as illustrated in Fig. 5.The pipe has the following properties: E = 200 5 N/m 2 , ν = 0.3, ρ = 1.0 3 kg/m 2 , L = 3 m, R = 0.05 m, and h = 0.01 m.In this simulation, 15 elements were used.Around the pipe cross section two Fourier terms were included for ovalization and warping deformations.This pipe was subjected to a time varying moment at one end, and a lateral force was applied at the other end (see Fig. 5).The moment was increased linearly with time beginning with zero and reaching a final value of 376 N m at τ = 0.2365 s.The final lateral force was 2 N (linearly increasing).The time step size is set to 0.0005 second and the relative error tolerance was set to 10 −3 .The deformation at several points in time are shown in Figs 6-7.For a better view of the deformation, the radius of the pipe was magnified by a factor in Figs 6-7.The computation time this simulation was approximately 3 hours on a Sun Sparc Station 5.

Concluding remarks
The pipe element described by Jiang and Arabyan [5] and this paper exhibits properties that may be very advantageous in the analysis of pipe like structures that undergo large and fast rotations and large deformations.Example of such structure range from umbilical cords used by astronautics and undersea divers to blood vessels in biosystems.
The principal feature of the new element is its ability to capture both rigid-body and constant-strain modes.This feature enables the element to be used in circumstances of combined bending, torsional and pressure loads, and large deformations including large cross-sectional ovalization and warping.This is a capability that other pipe elements (Bathe and Almeida [3], Ohtsubo [6], and ABAQUS theory manual [1]) do not have.
In addition, the new element is substantially more cost efficient computationally than other elements because a pipe-like structure can be modeled with few elements.This advantage was demonstrated earlier in Jiang and Arabyan [5].
In its present form the new element cannot model large strains (by virtue of the strain definition and the total Lagrangian formulation used).The element can be easily modified to model large strains by replacing the strain definitions and the constitutive laws with the appropriate relationships.Such modifications would also have a significant impact on Eqs (13) and ( 14) and make the dynamic formulation substantially more complex, but remove the thin-shell restriction.
shows the geometry of the pipe, which has the following properties: E = 200 5 N/m 2 , ν = 0.3, ρ = 1.0 3 kg/m 2 , L = 2 m, R = 0.05 m, and h = 0.01 m.The pipe is subjected to a timevarying moment at the cantilevered end, and the other is free.The external moment increases linearly with time, beginning with zero and reaching a final value of 692 N m.The pipe is discretized using ten elements.