A NOTE ON RELATIONSHIP BETWEEN FIXED-POLE AND MOVING ...

2 downloads 0 Views 211KB Size Report
Sep 14, 2012 - ... beam theory provided by Simo in 1985 [1] makes one of the mile- ..... tum, defined with respect to the fixed-pole rather than the reference line ...
European Congress on Computational Methods in Applied Sciences and Engineering (ECCOMAS 2012) J. Eberhardsteiner et.al. (eds.) Vienna, Austria, September 10-14, 2012

A NOTE ON RELATIONSHIP BETWEEN FIXED-POLE AND MOVING-POLE APPROACHES IN STATIC AND DYNAMIC ANALYSIS OF NON-LINEAR SPATIAL BEAM STRUCTURES Gordan Jeleni´c1 , Maja Ga´ceˇsa1 , and Miran Saje2 1 University

of Rijeka, Faculty of Civil Engineering Radmile Matejˇci´c 3, 51000 Rijeka, Croatia e-mail: {gordan.jelenic,maja.gacesa}@gradri.hr

2

University of Ljubljana, Faculty of Civil and Geodetic Engineering Jamova cesta 2, 1001 Ljubljana, Slovenia e-mail: [email protected]

Keywords: 3D beams, non-linear analysis, configuration-dependent interpolation, fixed-pole approach Abstract. Fixed-pole approach in the non-linear analysis of beam structures offers a different and promising framework for development of new element types with particular potential in the area of conservative time integrators. The finite elements formulated using this approach, however, utilise non-standard degrees of freedom which makes them difficult to combine with standard finite elements in existing codes. A modified approach which combines the benefits of the fixed-pole approach with the use of standard kinematic unknowns from the existing (movingpole) approaches is presented in this work along with a family of solution procedures for static analysis of geometrically non-linear spatial beam structures. Such a hybrid formulation opens a number of questions which are specifically addressed.

Gordan Jeleni´c, Maja Ga´ceˇsa, and Miran Saje

1

INTRODUCTION

1.1

Notation

L, x Ei I, 0 r, Λ u N, M n, m ¯ m ¯ n, Γ, K γ, κ ¯, κ ¯ γ ne , me ¯ e, m ¯e n F 0 , F L ,M 0 , M L ¯ 0, F ¯ L, M ¯ 0, M ¯L F CN , CM E, G A1 , A2 , A3 J, I2 , I3 A Jρ ρ k, π δr, δϑ N I i (x) G (if not in C N ) g i , q ii , q im , q ie K ij

initial length of the beam, co-ordinate of a point along the centroid axis of the beam orthonormal base vectors of the material co-ordinate system (i = 1, 2, 3) 3 × 3 unit and null matrix, respectively position vector of a point along the beam, rotation matrix of a triad rigidly attached to the cross-section displacement vector vectors of stress and stress-couple resultants in material co-ordinate system vectors of stress and stress-couple resultants in spatial co-ordinate system fixed-pole vectors of stress and stress-couple resultants in spatial co-ordinate system translational and rotational strain measures in material co-ordinate system translational and rotational strain measures in spatial co-ordinate system fixed-pole translational and rotational strain measures in spatial co-ordinate system vectors of applied distributed forces and torques fixed-pole vectors of applied distributed forces and torques applied concentrated forces at ends of the beam fixed pole applied concentrated forces at ends of the beam translational and rotational constitutive matrices Young’s (elastic) modulus, shear modulus cross-sectional and shear areas torsional and cross-sectional moments of inertia diag [A1 , A2 , A3 ] diag [J, I2 , I3 ] material density vectors of specific translational momentum and angular momentum variation of the position vector, spin variable number of interpolation points (nodes) of a finite element i-th Lagrangian interpolation polynomial of order N − 1 weak form residual and nodal vectors of internal, inertial and external forces at node i stiffness matrix relating the trial functions at node j to the residual at node i

˙ (•)0 , (•) c (•)t (•), (•)i , (•)j δ(•), ∆(•)

derivatives of (•) with respect to x, and time t cross-product operator (•)×, transposed quantity (•) (•) at nodes i and j variation of (•), incremental change of (•)

The geometrically exact beam theory provided by Simo in 1985 [1] makes one of the milestones in the development of the non-linear finite-element method in introducing a non-linear manifold as the solution space and the related differential geometry tools in order to implement a robust solution procedure [2]. This approach, however, turns out to be incompatible with algorithmic preservation of strain invariance with respect to a rigid-body motion, or simultaneous conservation of energy and the momentum vectors during a free motion of an unloaded beam. Borri and Bottasso [3, 4] propose the fixed-pole approach in which they introduce new stresscouple resultant and angular momentum - defined with respect to a unique point for the whole structure - the fixed pole. Different implementations of this general concept result in the algo2

Gordan Jeleni´c, Maja Ga´ceˇsa, and Miran Saje

rithms which naturally inherit the strain-invariance of the underlying formulation with respect to a rigid body motion [3] or are capable of simultaneous conservation of both energy and momentum vectors [4]. The main problem with the fixed-pole finite element implementation is that it uses nonstandard problem unknowns, and therefore is not easily combinable with existing finite element codes. For this reason, we derive and present a formulation which uses standard kinematic unknowns but is inherently based on the fixed-pole approach, along with with the algorithmic benefits of this approach. 2

KINEMATIC AND CONSTITUTIVE EQUATIONS

Kinematic equations follow from the Reissner-Simo theory [1]. The material translational and rotational strain measures are given as follows: Γ = Λt γ = Λ t r 0 − E 1 ˆ = Λt κ b Λ = Λt Λ 0 , K

(1)

while the spatial translational and angular velocity vectors, v and w, are defined as v = r˙ and ˙ Varying (1) leads to b = Λt Λ. w δγ = δr 0 + r 0 × δϑ (2) δκ = δϑ0 . Let us consider a linear elastic material with the material cross-sectional stress and stresscouple resultants defined as      N CN 0 Γ = , (3) M 0 CM K with CN 3

  EA 0 0 0 , =  0 GA1 0 0 GA2

CM

  GIt 0 0 =  0 EI1 0  . 0 0 EI2

(4)

THE STANDARD APPROACH

The standard approach [2] given by Simo and Vu-Quoc follows naturally from the principle of virtual work: Vi + Vm − Ve = 0 ,

(5)

where Vi , Vm and Ve are the virtual works of internal, inertial and external forces, respectively, which are defined as follows:      d Z L Z L I dx 0 n n t t t t dx , (6) Vi = hδγ δκ i dx = hδr δϑ i 0 I d b m m − r 0 0 dx   Z L k˙ t t Vm = hδr δϑ i dx , (7) π˙ 0       Z L ne F0 FL t t t t t t Ve = hδr δϑ i dx + hδr 0 δϑ0 i + hδr L δϑL i . (8) me M0 ML 0 3

Gordan Jeleni´c, Maja Ga´ceˇsa, and Miran Saje

A certain irregularity in the integrals should be noted – while the virtual quantities simply dot-multiply the specific vectors of inertial and external forces, an additional operator occurs between the virtual quantities and the stress resultants:   d I dx 0 . (9) d −rb0 I dx After interpolating virtual quantities δr, δϑ using Lagrangian polynomials I i (x): δr =

N X

i

I (x)δr i ,

δϑ =

i=1

N X

I i (x)δϑi ,

(10)

i=1

we have the approximate weak form: h

G ≡

N X

hδr i t δϑi t ig i = 0 ,

(11)

i=1

from which, because of the arbitrariness of δr i and δϑi , follows the vector of nodal dynamic residuals: g i = q ii + q im − q ie = 0 ,

(12)

with q ii , q im and q ie as the nodal vectors of internal, inertial and external forces, respectively, given as follows:   Z L  i0 I I 0 n i qi = dx , (13) i0 i b0 m −I r I I 0 Z L   k˙ i dx , (14) Ii qm = π˙ 0     Z L   ne F0 FL i i i i qe = I dx + δ1 + δN . (15) me M0 ML 0 4

THE FIXED-POLE APPROACH

In [3, 4] Borri and Bottasso thoroughly investigated the idea of replacing the stress-couple resultant m and the specific angular momentum π, which are defined with respect to the beam ¯ and specific angular reference axis at a cross-section, with another stress-couple resultant m ¯ which are to be defined with respect to a unique point for the whole structure— momentum π, the fixed pole—naturally taken to be the origin of the spatial frame, ie: ¯ = r × n + m, m ¯ = r×k+π. π

(16) (17)

¯ = k. In order to obtain a relationship between standard ¯ = n and k We also note that n and fixed-pole strain measures, we request the strain energy density to remain invariant to the change of virtual quantities: 1 1 ¯ +κ ¯ · m) ¯ = (γ · n + κ · m) (¯ γ·n 2 2 4

(18)

Gordan Jeleni´c, Maja Ga´ceˇsa, and Miran Saje

By substituting the relationship (16) between standard and fixed pole stress resultants, we obtain the relationship between standard and fixed-pole spatial strain measures: ¯ =γ + r × κ, γ ¯ =κ. κ

(19) (20)

It is also useful to note the relationship between the time derivatives of standard and fixedpole specific momenta: ¯˙ =k˙ , k π ¯˙ = r × k} + r × k˙ + π˙ , |˙ {z

(21)

˙ ˙ r×Aρ r=0

and considering virtual work of the inertial forces (7) to be an objective quantity gives us the fixed-pole kinematics which is virtual-work conjugate to the time-derivatives of the specific fixed-pole momenta as δ¯ r = δr + r × δϑ , (22) ¯ δ ϑ = δϑ , from where their respective derivatives follow as δ¯ r 0 = δr 0 + r 0 × δϑ + r × δϑ0 , ¯ 0 = δϑ0 . δϑ

(23)

Analogously to (16), the fixed-pole distributed moment is introduced: ¯ e = me + r × ne , m

(24)

¯ e = ne . We also introduce the relationship between standard and fixed-pole noting that n concentrated moments at ends of the beam: ¯ 0 = M 0 + r0 × F 0 , M ¯ L = M L + rL × F L , M

(25)

¯ 0 = F 0 and F ¯ L = F L. noting that F 4.1

Original formulation

We rewrite the expression for the virtual work of internal forces (6) and substitute (23) and (16) into it:   Z L Z L n t 0t 0 0 Vi = h(δr + r × δϑ) δϑ i dx = ((δr 0 + r 0 × δϑ) · n + δϑ0 · m) dx = m 0 0 Z L ¯ − r × n)) dx = ((δ¯ r 0 − r 0 × δϑ − r × δϑ0 + r 0 × δϑ) · n + δϑ0 · (m 0 Z L  ¯0 · m ¯ + δϑ ¯ dx , = δ¯ r0 · n 0

(26)

5

Gordan Jeleni´c, Maja Ga´ceˇsa, and Miran Saje

and obtain the virtual work of internal forces using the fixed-pole virtual quantities and stress and stress-couple resultants:   Z L  d ¯ t t ¯i n hδ¯ Vi = dx . (27) r δϑ ¯ m dx 0 Note that there is no irregularity (9) in this term! In a like manner, we rewrite the expression for virtual work of the inertial forces (7) and substitute (22) and (21) into it: ( ) Z L ˙ ¯ k h(δ¯ dx = Vm = r − r × δϑ)t δϑt i π ¯˙ − r × k˙ 0 Z L Z L   ˙ ¯˙ + δ ϑ ¯ ·π ˙ ¯ δ¯ r·k ¯˙ dx , (28) (δ¯ r − r × δϑ) · k + δϑ · π ¯˙ − r × k dx = = 0

0

to get the fixed-pole virtual work of inertial forces as   Z L ¯˙ t k t ¯ Vm = hδ¯ dx . r δϑ i π ¯˙ 0

(29)

Finally, the virtual work of external forces (8) may be treated in the same way, so that after substituting (25) and (22) into it:   Z L ¯e n t t Ve = h(δ¯ dx+ r − r × δϑ) δϑ i ¯ e − r × ne m 0   ¯0 F t t +h(δ¯ + r 0 − r 0 × δϑ0 ) δϑ0 i ¯ M 0 − r0 × F 0   ¯L F t t , (30) +h(δ¯ r L − r L × δϑL ) δϑL i ¯ M L − rL × F L we arrive at the following fixed-pole virtual work of external forces:       Z L ¯0 ¯L ¯e F F n t t t t t t ¯ ¯ ¯ hδ¯ dx + hδ¯ + hδ¯ . Ve = r δϑ i r 0 δ ϑ0 i ¯ r L δ ϑL i ¯ ¯e m M0 ML 0 ¯ using Lagrangian polynomials Deciding to interpolate the fixed-pole virtual quantities δ¯ r, δϑ I i (x), δ¯ r=

N X

¯= δϑ

i

I (x)δ¯ ri ,

i=1

N X

¯i , I i (x)δ ϑ

(31)

i=1

we obtain the approximate weak form: Gh ≡

N X i ¯ t i¯ hδ¯ r ti δ ϑ i g = 0,

(32)

i=1

¯ i, from which follows, because of the arbitrariness of δ¯ r i and δ ϑ g¯ i = q¯ ii + q¯ im − q¯ ie = 0 , 6

(33)

Gordan Jeleni´c, Maja Ga´ceˇsa, and Miran Saje

with g¯ i as the fixed-pole vector of nodal dynamic residuals and q¯ ii , q¯ im and q¯ ie as the fixed-pole nodal vectors of internal, inertial and external forces, respectively, given as Z L   ¯ n 0 i Ii dx , (34) q¯ i = ¯ m 0 Z L ˙  ¯ k i Ii dx , (35) q¯ m = π ¯˙ 0     Z L   ¯L ¯0 ¯e F n F i i i i . (36) + δN ¯ q¯ e = I dx + δ1 ¯ ¯e ML M0 m 0 4.2

Modified formulation

As shown in (32), the fixed-pole nodal dynamic residual g¯ i is virtual-work conjugate to the fixed-pole nodal virtual displacements and rotations which are non-standard and prevent the elements from being combined with standard elements in a finite-element mesh. This problem, however, is easily overcome by noting that from (22) the nodal fixed-pole virtual quantities are related to the standard nodal virtual quantities via:        δ¯ ri δr i + r i × δϑi I rbi δr i = = , (37) ¯ δ ϑi δϑi 0 I δϑi so that the fixed-pole weak form (32) becomes h

G ≡

N X

hδr ti

δϑti i



 N I 0 i X t gi = 0 , g¯ = hδr i δϑti i˜ −rbi I

(38)

i=1

i=1

with g˜ i = q˜ ii + q˜ im − q˜ ie ,

(39)

and Z

L



I I = d r − ri 0 Z L  I Ii d q˜ im = r − ri 0 Z L  I q˜ ie = Ii d r − ri 0

q˜ ii

i0

 n dx , m   0 k˙ dx , I π˙       0 ne F0 FL i i dx + δ1 + δN . me M0 ML I

0 I



(40) (41) (42)

After substituting (40), (41) and (42) in g˜ i and then back to (32), it can be shown that, using a non-linear interpolation of the type: δr =

δϑ =

N X i=1 N X

I i [δr i + (r i − r) × δϑi ] (43) I i δϑi

i=1

is equivalent to the standard (Lagrangian) interpolation of the fixed-pole virtual quantities (31). 7

Gordan Jeleni´c, Maja Ga´ceˇsa, and Miran Saje

5

SOLUTION PROCEDURE

The non-linear equation (39) may be solved using the Newton-Raphson solution procedure. In order to concentrate on the issues that arise when attempting to utilise a fixed-pole approach in the hybrid manner presented here, in which the standard displacement and rotation quantities have been kept as the problem unknowns, we proceed by considering only the d and static case. When linearising (39), account has to be taken of the fact that ∆Λ = ∆ϑΛ t 0 0 ∆Γ = Λ (∆r + r × ∆ϑ) as well as [1]    t    0 b Λold b Λold = cnew = Λt Λ0 = exp ∆ϑ K exp ∆ ϑ new new   t   0   t   b b b b Λ0 = Λtold exp ∆ϑ exp ∆ϑ Λold + Λtold exp ∆ϑ exp ∆ϑ old {z } | I

c+K cold , = ∆K

(44)

where   b + 1 − cos ∆ϑ ∆ϑ b2 , b = I + sin ∆ϑ ∆ϑ exp ∆ϑ ∆ϑ ∆ϑ2

(45)

∆K = Λtold H∆ϑ0 ,

(46)

and

with H=

∆ϑ − sin ∆ϑ sin ∆ϑ 1 − cos ∆ϑ b ∆ϑ∆ϑt + I− ∆ϑ . 3 ∆ϑ ∆ϑ ∆ϑ2

(47)

The strain measure updates for K and Γ then follow as K new = K old + ∆K , Γnew := Λtnew r 0new − E 1 .

(48) (49)

After expanding equation (39) into a Taylor series around a known configuration and omitting higher-order terms we have g ˜i + ∆˜ gi = 0 ,

(50)

∆˜ g i = ∆˜ q ii − ∆˜ q ie .

(51)

where

After the linearisation of the modified nodal internal force vector (40) we have

∆˜ q ii

  Z L  0 0 ∆r i − ∆r I 0 = I dx + Ii b 0 ∆ϑi − ∆ϑ b − rbi n r 0 0    t Z L  I 0 Λ 0 CN 0 Λ 0 + Ii b b r − r I 0 Λ 0 C 0 i M 0 Z

L

i0



8

   0 0 −b n ∆r dx+ I 0 −c m ∆ϑ  0  0 ∆r + r 0 × ∆ϑ dx . Λt ∆ϑ0 (52)

Gordan Jeleni´c, Maja Ga´ceˇsa, and Miran Saje

Linearizing the modified nodal external force vector (42) we have:   Z L  0 ∆r i − ∆r i i 0 dx . I ∆˜ qe = ce 0 ∆ϑi − ∆ϑ n 0

(53)

An interesting phenomenon occurs when linearising the external force vector, which is a direct consequence of the use of the fixed-pole virtual quantities and their conversion back to standard virtual quantities: even though the external loading is configuration-independent, the modified nodal external loading vector (42) is not. As a result of this, the stiffness matrix will have an additional part as a result of the linearisation of q˜ei . In order to complete the solution procedure, the unknown functions and/or their iterative changes must be interpolated in some way. Within this paper we consider three different interpolation options, which arise as a consequence of the fact that the virtual position vector has been interpolated in a non-linear manner, dependent on the actual position vector not only at the nodal points, but also at x (i.e. the integration point). 5.1

Interpolation option 1

Within this option ∆r and ∆ϑ are interpolated in the same way as δr and δϑ (43): ∆r =

N X

I j [∆r j + (r j − r) × ∆ϑj ]

j=1

∆ϑ =

N X

(54) I j ∆ϑj ,

j=1

while the unknown displacement functions are interpolated as: r=

N X

I k rk .

(55)

k=1

It should be noted that this interpolation is in contradiction with the interpolation for ∆r in (54), which obviously does not follow as a linearisation of (55). The above contradiction necessarily results in loss of quadratic convergence of the Newton-Raphson solution process, and has been introduced as a simplest means of providing r(x) in (54). After introducing (54) into (52) and (53) we have:   ∆r j ij i ∆˜ q i = K int , (56) ∆ϑj   ∆r j ij i ∆˜ q e = K ext , (57) ∆ϑj where K ij int

Z = 0

L

   Z L 0 0 −b n i0 j 0 I dx + I I dx b (r jd c (δij − I j )b n −I j n − r) 0 −(r d − r i )b n−m 0     t   Z L b I 0 Λ 0 CN 0 Λ 0 I rbj − r i0 j 0 + I I dx , b − rbi I 0 Λ r 0 CM 0 Λt 0 I 0 (58) i0



9

Gordan Jeleni´c, Maja Ga´ceˇsa, and Miran Saje

and K ij ext

Z

L

=−

I

i



0

 0 0 dx , ce (r jd (I j − δij )c ne I j n − r)

(59)

and finally ij K ij = K ij int − K ext .

5.2

(60)

Interpolation option 2

Within this option, r is interpolated using (55). Interpolation of the displacement increments follows consistently by linearising r and ∆ϑ: ∆r =

N X

j

I ∆r j ,

∆ϑ =

j=1

N X

I j ∆ϑj .

(61)

j=1

After introducing (61) into (52) and (53) we have   ∆r j ij i , ∆˜ q i = K int ∆ϑj   ∆r j ij i ∆˜ q e = K ext , ∆ϑj

(62) (63)

where K ij int

L

   Z L 0 0 I 0 0 j i0 j = I (δij − I ) dx + I I b 0 b − rbi I 0 r n 0 0    t  0 Z L  0 Λ 0 Ij I I 0 Λ 0 CN i0 + I b − rbi I 0 Λ r 0 CM 0 Λt 0 0 Z

i0



 −b n dx+ −c m  I j rb0 dx , 0 Ij I

and K ij ext

Z =− 0

L

 0 0 dx , I (I − δij ) ce 0 n i

j



(64)

and finally ij K ij = K ij int − K ext .

(65)

In this interpolation option, different interpolations have been used for the test functions and the trial functions, which is bound to make the tangent stiffness matrix strongly non-symmetric. 5.3

Interpolation option 3

Within this option r is not interpolated, but only updated as follows: r new (x) = r old (x) + ∆r(x) ,

(66)

where r old (x) is the last known value for r(x), not necessarily associated with an equilibrium state and ∆r(x) following interpolation of the Newton-Raphson changes given in (54). This is a consistent choice which yields a tangent stiffness matrix. In this case the values of r at both integration and nodal points must be saved. The earlier results (56)–(60) can still be used, but now noting that: r = r g = r(xg ) . 10

(67)

Gordan Jeleni´c, Maja Ga´ceˇsa, and Miran Saje

6

CONCLUDING REMARKS

The fixed-pole approach, first introduced by Borri and Bottasso [3] has been analysed and discussed. The approach uses a new set of stress-couple resultants and specific angular momentum, defined with respect to the fixed-pole rather than the reference line at the cross-section, which in turn implies that r(x), the standard position vector, ceases to be the problem unknown. The presented modified fixed-pole approach makes it easy to combine the standard choice of the problem unknowns with the benefits of the original fixed-pole approach. The solution procedure for a spatial beam under static loading has been presented and the complexities of the modified approach discussed. In particular, the approach is based on a non-linear interpolation of the virtual displacements, which is dependent on the position vector field itself, which in turn opens a number of combinations in interpolating the position vector field and its iterative change. The solution procedure for a spatial beam under dynamic loading will be derived and the numerical examples shown in future works. ACKNOWLEDGEMENT Research resulting with this paper was made within the scientific project No 114-00000003025: “Improved accuracy in non-linear beam elements with finite 3D rotations” financially supported by the Ministry of Science, Education and Sports of the Republic of Croatia. The second author additionally acknowledges the support of the Croatian Science Foundation which funded the bilateral project No 03.01/129: “Preservation of mechanical constants in numerical time integration of non-linear beam equations of motion”. REFERENCES [1] J.C. Simo. A finite strain beam formulation. The three-dimensional dynamic problem. Part I. Computer Methods in Applied Mechanics and Engineering, 49:55 – 70, 1985. [2] J.C. Simo and L. Vu-Quoc. On the dynamics in space of rods undergoing large motions — A geometrically exact approach. Computer Methods in Applied Mechanics and Engineering, 66:125–161, 1986. [3] M. Borri and C. Bottasso. An intrinsic beam model based on a helicoidal approximation— Part I: Formulation. International Journal for Numerical Methods in Engineering, 37(13):2267–2289, 1994. [4] C. Bottasso and M. Borri. Integrating finite rotations. Computer Methods in Applied Mechanics and Engineering, 164(3-4):307–331, 1998.

11