beam structures, a1, section beam, material behaviour, finite element formulation, torsion, beam, beam strains, reference point, constitutive equations, strains, formulation, ai, Lagrangian formulation, finite deformations, deformation, kinematic quantities, plastic strains, bi, Green's formula, Di, a2, a3, finite element interpolation, finite element, finite elements, stress, stress resultants, epl, finite element formulations
Content:
Published in: Int. J. Num. Meth. Engng., 48, 2000, 16751702. THEORY AND NUMERICS OF THREEDIMENSIONAL BEAMS WITH ELASTOPLASTIC MATERIAL BEHAVIOUR F. Gruttmann1, R. Sauer2, W. Wagner3 1 Institut fuЁr Statik, Technische UniversitaЁt Darmstadt, 64283 Darmstadt, Germany , 2 IngenieurbuЁro JЁager, 01445 Radebeul, Germany, 3 Institut fuЁr Baustatik, UniversitЁat Karlsruhe, 76128 Karlsruhe, Germany Keywords: Threedimensional beams, finite deformations, torsion warping deformation, arbitrary
cross sections, elastoplastic material behaviour, finite elements ABSTRACT A theory of space curved beams with arbitrary crosssections and an associated finite element formulation is presented. Within the present beam theory the reference point, the centroid, the center of shear and the loading point are arbitrary points of the cross section. The beam strains are based on a kinematic assumption where torsionwarping deformation is included. Each node of the derived finite element possesses seven
DEGREES OF FREEDOM. The update of the rotational parameters at the finite element nodes is achieved in an additive way. Applying the isoparametric concept the kinematic quantities are approximated using Lagrangian interpolation functions. Since the reference curve lies arbitrarily with respect to the centroid the developed element can be used to discretize eccentric stiffener of shells. Due to the implemented constitutive equations for elastoplastic material behaviour the element can be used to evaluate the load
carrying capacity of beam structures. 1 INTRODUCTION Threedimensional beamlike structures undergoing finite deformations occur in different areas of engineering practice. Here motions of flexible beams like helicopter blades, rotor blades, robot arms or beams in spacestructure technology are mentioned. Furthermore, the analysis of load carrying capacities requires the implementation of constitutive equations for inelastic material behaviour. Numerous papers have been published up to now using different approaches. Some of them are discussed in the following. Due to the noncommutativity of successive finite rotations about fixed axes, Argyris et al.1 introduced the socalled semitangential rotations to circumvent this difficulty. Bathe and Bolourchi2 developed an updated and a total Lagrangian formulation for beams with large displacements. Thin structures undergoing finite deformations are usually characterized by significant rigid body motions. This motivates the socalled corotational dedicated to F.G. Kollmann on the occasion of his 65th birthday 1
formulation where the rigid body motions are separated from the total deformation (e.g. see Belytschko and Hsieh,3 Crisfield4 or NourOmid5). Applying this procedure existing linear elements can be used for nonlinear computations. Several authors developed finite element formulations for threedimensional beams using beam strains derived from the internal virtual work, see Reissner.6 Here we mention the papers of Simo and VuQuoc7,8 and Cardona and Gґeradin.9 The authors in [7] show that the linearization of the virtual work principle yields a nonsymmetric geometric tangent stiffness matrix when applying a multiplicative update of the rotation tensor. The tangent matrix is only symmetric at a point of equilibrium assuming conservative loads. Ibrahimbegoviґc10 derived symmetric tangent matrices when using an additive update of the axial vector. Most of the finite element formulations are concerned with beams where centroid and center of shear coincide. The problem of coupled bending torsion deformation of beams has been studied theoretically e.g. by Reissner.11,12 In Ref. [8] torsionwarping deformation has been incorporated within the theory and the associated finite element formulation. In this paper we derive a threedimensional finite beamelement with arbitrary space curved reference axis. The novel aspects and essential features of the formulation are summarized as follows: (i) The introduced beam strains based on a kinematic assumption including torsion warping deformation are conform with the strains in Ref. [8]. Here, we additionally derive the relation to the GreenLagrangian strains and the associated variations. In our formulation the matrix which relates the variation of the GreenLagrangian strains to the variational beam strains depends on the deformation. Thus the complete nonlinearity within the underlying kinematic assumption is considered in the present theory. Due to this fact the numerical results show a good agreement with the results of higher valued models like shell discretizations for thinwalled structures. (ii) The constitutive equations for elastoplastic material behaviour are implemented. Within the present beam theory the
normal stress and the shear stresses enter into the von Mises yield condition and the associated flow rule. Thus interaction between the different stress components is included. The nonlinear stressstrain relations requires a
Numerical Integration over the crosssections. The developed element can be used to analyze the carrying capacity of beam structures. (iii) For linear elasticity and small strains the material matrix can be integrated analytically. For this purpose one has to apply the equations of SaintVenant torsion theory along with Green's formula to integrate the socalled warping coordinates. This yields a material matrix in terms of section quantities which describes the coupling effects. In contrast to Ref. [8] the reference point is an arbitrary point of the crosssection. (iv) The finite element formulation is presented using Lagrangian interpolation functions. Each node possesses seven degrees of freedom at the nodes. The update of the nodal rotation quantities is performed in an additive way. Within this procedure additional storage of the
rotation matrix with nine parameters can be avoided. The linearization yields a symmetric tangent matrix. However, due to the chosen finite 2
element approximation the expressions are essentially simpler compared with Ref. [10]. The external loading can be applied at an arbitrary point of the crosssection. The contribution of the loads to the stiffness matrix is derived. Loading with stress couple resultants leads to a symmetric load stiffness matrix. This is in contrast to the multiplicative procedure, see Ref. [7]. The contents of the paper is outlined as follows: In the next section we present the kinematics of a space curved beam. The beam strains are defined and the relation to the GreenLagrangian strains is shown. In section 3 the underlying variational formulation of the boundary value problem is given using a Lagrangian representation. The associated EulerLagrange equations are derived. Furthermore the linearized variational formulation is described. In section 4 the basic equations for elastoplastic material behaviour are given. For linear elasticity and small strains the material matrix which describes the torsion bending coupling is derived in terms of section quantities. The associated finite element formulation is given in section 5. The applicability of the developed formulation is demonstrated in section 6 with several examples. Comparisons with results obtained with shell discretizations are given.
2 KINEMATIC DESCRIPTION OF THE BEAM A beam with reference configuration denoted by B0 according to Fig. 1 is considered. Assuming arbitrary crosssections the centroid S and the center of shear M are independent of the reference point O. An orthogonal basis system Ai with local coordinates {1, 2, 3} is introduced. The axis of the beam is initially along A1 with the arclength parameter S = 1 [0, L] of the spatial curve. The crosssections of the beam therefore lie in planes described by the basis vectors {A2, A3}. Accordingly the frame ai is defined in the current configuration which is characterized by the time parameter t. Note, that within the underlying beam kinematic the vector a1 is not tangent vector of the deformed reference curve. The basis vectors Ai and ai follow from the orthogonal transformations
Ai(S) = R0(S) ei , ai(S, t) = R(S, t) ei with R0, R SO(3) . (1)
Several parametrizations of orthogonal tensor using Eulerian angles, Cardan angles, quaternions etc. have been discussed in the literature, see e.g. Gґeradin and Rixen,13 Betsch et al.14 The position vectors of the undeformed and deformed crosssections are given with the following kinematic assumption
X(2, 3, S) = X0(S) + 2A2(S) + 3A3(S)
(2)
x(2, 3, S, t) = x0(S, t) + 2a2(S, t) + 3a3(S, t) + (S, t) wЇ(2, 3) a1(t)
where a1 is assumed to be piecewise constant. The given warping function wЇ(2, 3) is defined within the SaintVenant torsion theory for bars. This is a basic assumption which restricts the present formulation to a certain class of problems. Our numerical investigations however showed good agreement of the results for thinwalled beam structures obtained with the present beam model compared with a higher valued shell model.
3
Ai
u
B
B0
x3
P
p3 ai
s3, m3
MS
X0
x0
ei
A3
0
A2
x 1 A1 m2 s2 p2
x2
Figure 1: Initial and current configuration of the beam
For linear elasticity and pure torsion the boundary value problem is given with the following Neumann problem, see e.g. Timoshenko and Goodier15
wЇ,22 +wЇ,33 = 0
in A
(3)
wЇ,2 n2 + wЇ,3 n3 = 3 n2  2n3 on C ,
where the outward normal vector n = [n2, n3]T is defined on the boundary C. Here, the notation (·), is used to denote partial derivatives with respect to the coordinates . The solution of (3) along with the normality conditions
w~ dA = 0 ,
w~ 2 dA = 0 ,
w~ 3 dA = 0
(4)
A
A
A
define the center of shear M with coordinates {m2, m3}. The warping function w~ refers
to M and reads
w~ = wЇ + m2Ї3  m3Ї2
(5)
where the coordinates Ї2 = 2  s2 and Ї3 = 3  s3 intersect at the centroid S.
Based on the kinematic assumption (2) the tangent vectors Gi = X,i and gi = x,i are
derived
G1 = X0 + 2A2 + 3A3
g1 = x0 + 2a2 + 3a3 + wЇa1
G2 = A2
g2 = a2 + wЇ,2 a1
(6)
G3 = A3
g3 = a3 + wЇ,3 a1
where (·) denotes the customary symbol for differentiation with respect to the arclength. The derivative of the orthogonal basis systems can be expressed using the vector products Ai = 0 Ч Ai and ai = Ч ai, where 0 and denote the socalled axial vectors. Pull back of the covariant basis vectors with the rotation tensors R0 = Ai ei and R = ai ei yields
F1 := RT0 G1 = 0 + 0 Ч d F2 := RT0 G2 = e2 F3 := RT0 G3 = e3
f1 := RT g1 = + Ч d + wЇe1
f2 := RT g2 = e2 + wЇ,2 e1
(7)
f3 := RT g3 = e3 + wЇ,3 e1
4
with d = 2e2 + 3e3. The strain vectors of the current configuration are defined with
:= RT x0
=
x0 · a1 x0 · a2
:= RT
=
a2 · a3 a3 · a1
.
(8)
x0 · a3
a1 · a2
The corresponding quantities of the reference configuration read 0 = RT0 X0 and 0 = RT0 0. Eqs. (6) are inserted into the GreenLagrangian strain tensor E = EijGi Gj, where the dual basis vectors Gi are defined in a standard way Gi · Gj = ij. The components which contribute to the virtual work read
E =
E11 2E12
=
1 2
(g11

G11)
g12  G12
(9)
2E13
g13  G13
with the metric coefficients
Gij = Gi · Gj = Gi · R0RT0 Gj = Fi · Fj
(10)
gij = gi · gj = gi · RRT gj = fi · fj
of the reference configuration and current configuration, respectively.
3 VARIATIONAL FORMULATION OF THE BOUNDARY VALUE PROBLEM In this section the virtual work of the stresses and the external forces are derived considering the beam kinematic. For this purpose stress resultants are defined. The associated EulerLagrange equations and the linearization of the virtual work expressions are given.
3.1 Internal Virtual Work and Definition of Stress Resultants Volume forces 0bЇ and applied surface loads Їt are acting on the considered body. Hence the equilibrium is given here in weak form
g(v, v) = S · E dV  0bЇ · u dV  Їt · u d = 0 .
(11)
B0
B0
B0
The independent kinematic quantities of the beam are v = [u, R, ]T , where u = x0  X0 and R = ai ei denote the displacement vector and the rotation tensor of the reference curve according to (2), respectively. The space of kinematically admissible variations is introduced by
V := {v = [u, w, ]T : [0, L]  R3 v = 0 on Su}
(12)
where Su describes the boundaries with prescribed displacements and rotations. Here, the axial vector w is defined by ai = RRT ai = w Ч ai. The GreenLagrangian strain tensor is work conjugate to the Second PiolaKirchhoff stress tensor S = Sij GiGj. Within the present beam theory the components S22, S33, S23
5
are neglected. Using matrix notation the vector of nonvanishing stress components is
defined
S = S11, S12, S13 T .
(13)
The variation of the work conjugate GreenLagrangian strains (9) yields with (7) (10)
E
=
E11 2E12
=
f1 · f1 f2 · f1 + f1 · f2
(14)
2E13
f3 · f1 + f1 · f3
where
f1 = + Ч d + wЇe1
f2 = wЇ,2 e1
(15)
f3 = wЇ,3 e1 .
Introducing
F^ = [f1, f2, f3]
0 Wd = skew d = 3
3 0
2 0
(16)
2 0 0
a = (f1 · e1) (wЇ,2 e2 + wЇ,3 e3)
one obtains
E = A E^ A = [F^ T , F^ T WdT , a, wЇ F^ T e1]
E^ =
.
(17)
The variation of the beam strains (8) yields
= RT x0 + RT x0 = RT (x0  w Ч x0)
(18)
= RT + RT = RT w .
Eq. (18)1 is evident, whereas a proof of (18)2 is given e.g. in [17]. Furthermore a representation of the components yields
=
a1 · u a2 · u
+ x0 · a1 + x0 · a2
a3 · u + x0 · a3
=
a3 · a2 + a2 · a3 a1 · a3 + a3 · a1
.
(19)
a2 · a1 + a1 · a2
Using (17) the internal virtual work in (11) can be rewritten with dV = dA dS as
gint(v, v) = ET S dV = E^ T S^ dS
B0
S
(20)
= (F · + M · + F w + M w ) dS
S
6
with the vector of stress resultants
F
T
S^ =
AT
S
dA
=
M Fw
=
L Tw
dA
(21)
A
Mw
A Lw
where T = [T 11, T 12, T 13]T = F^ S and L = Wd T. Finally a component representation of the vector of stress resultants yields
F1
T 11
S^ =
F2 F3 M1 M2 M3 Fw
= A
T 12 T 13 T 132  T 123 T 113 T 112 T~12wЇ,2 +T~13wЇ,3
dA
Mw
T 11wЇ
(22)
with T~1 = (f1 · e1)S1 and = 2, 3. Herein, F 1 is the normal force, F 2 and F 3 the shear forces, M 1 the torsion moment, M 2 and M 3 the
Bending Moments, F w the bishear and M w the bimoment, respectively.
3.2 External Virtual Work and EulerLagrange Equations In this paper the volume forces 0bЇ and the surface loads Їt according to (11) are considered with the load pЇ = pЇ(S) acting at the coordinates {p2, p3}, see Fig.1. The vector of the loading point in the current configuration reads xp = x0+rp with rp = p2a2+p3a3+ wЇp a1 and wЇp = wЇ(p2, p3). Hence the external virtual work reads
gext(v, v) =  pЇ · xp dS
S
=  pЇ · (x0 + p2a2 + p3a3 + wЇpa1 + wЇpa1) dS
(23)
S
=  (pЇ · x0 + mЇ · w + pЇ1) dS
S
with mЇ = rp Ч pЇ and pЇ1 = wЇp (pЇ · a1) .
The weak form of equilibrium (11) is now reformulated using the expressions (20), (18) and (23). This yields
g(v, v) = [f · (x0  w Ч x0) + m · w + F w + M w  (pЇ · x0 + mЇ · w + pЇ1)] dS = 0 S (24)
7
with f := RF and m := RM. Next,
Integration by parts yields with homogeneous stress
boundary conditions
g(v, v) =  [(f + pЇ) · x0 + (m + x0 Ч f + mЇ ) · w + (M w  F w + pЇ1) ] dS = 0 . S (25) Applying standard arguments from variational calculus we obtain
f + pЇ = 0
m + x0 Ч f + mЇ = 0
(26)
M w  F w + pЇ1 = 0 .
These are the EulerLagrange equations of the variational formulation. The first two equations in (26) are wellknown equilibrium equations of a threedimensional beam using a vector notation. The third equation is identical to that of the linear theory. Reissner12 derived (26)3 in the context of a secondorder geometrical nonlinear theory.
3.3 Linearization of the Virtual Work Expressions For the subsequent finite element formulation we need to derive the linearization of the weak form of equilibrium (11). This is formally achieved using the directional derivative
L [g(v, v)] = g + Dg · v ,
Dg
· v
=
d d
[G(v
+
v)]=0
(27)
where v = [u, w, ]T . The linearization of the internal virtual work yields a material part and a geometrical part
Dgint(v, v) · v = (E^ T D^ E^ + E^ T S^) dS .
(28)
S
Here the linearized beam strains E^ = [, , , ]T are obtained from (19) re
placing the operator by . Furthermore the linearized virtual strains are expressed by E^ = [, , 0, 0]T with
=
u u
· a1 + a1 · u · a2 + a2 · u
+ x0 · a1 + x0 · a2
u · a3 + a1 · u + x0 · a3
(29)
=
a2 · a3 + a3 · a2 + a3 · a2 + a2 · a3 a3 · a1 + a1 · a3 + a1 · a3 + a3 · a1
a1 · a2 + a2 · a1 + a2 · a1 + a1 · a2 .
Next the linearization of the stress resultants (21) applying the product rule leads to
D^
=
d S^ d E^
=
d d E^ [
AT S dA] =
(AT CT A + G^ ) dA .
(30)
A
A
The first part follows from differentiation of the stress vector S and the incremental GreenLagrangian strains dE = A dE^ . The tangential matrix CT := dS/dE is specified
8
for elasticplastic material behaviour in the next section. The second part follows from differentiation of the matrix product AT S where the stress vector S is held fixed. With
T L Tw
=
S11f1 + S12f2 + S13f3 Wd(S11f1 + S12f2 + S13f3) S~(f1 · e1)
Lw
S11wЇ(f1 · e1) + S~wЇ
(31)
where S~ = S12wЇ,2 +S13wЇ,3 and fi according to (7) one obtains
G^
=
T
L
(
T w
)T
T
L
(
T w
)T
T L T w
T S111
L T w
=
S 11 Wd S~eT1
S 11 WdT S 11 Wd2 S~d~ T
S~e1 S~d~ 0
S11wЇe1
S11wЇd~ S~wЇ
(
Lw
)T
(
Lw
)T
Lw
Lw
S11wЇeT1 S11wЇd~T S~wЇ S11wЇ2
(32)
where d~ = Wde1 = 3e2  2e3. The integration over the crosssection in (21) and (30) is performed numerically. The warping function wЇ is determined for arbitrary crosssections in a separate computing process using the finite element method as is discussed in [16]. The finite element mesh for the computation of wЇ is used here to perform the numerical Gauss integration. In case of small strains F^ 1 holds, A becomes independent of the displacement field and G^ vanishes. Finally we derive the linearization of the external virtual work as
Dgext(v, v) · v =  pЇ · xp dS
(33)
S
with xp = p2a2 + p3a3 + wЇpa1 + wЇpa1 + wЇpa1.
4 MATERIAL LAW AND STRESS RESULTANTS In this section we present the constitutive equations for elastoplastic material behaviour. We consider
metallic materials which can be described by the von Mises yield criterion with isotropic hardening and associated flow rule. For elasticity the section integrals are reformulated using the equations of SaintVenant's torsion theory and Green's theorem. This yields an elasticity matrix where for the general case torsion bending coupling occurs.
4.1 ElasticPlastic
stress analysis The GreenLagrangian strains (9) are decomposed in an elastic and plastic part
E = Eel + Epl
(34)
which holds for small elastic and plastic strains. The elastic part is described by the socalled St.VenantKirchhoff material law. Thus, a quadratic strain
energy function is postulated where the stresses are obtained by partial derivatives
Ws(Eel)
=
1 EelT CEel 2
,
S = Ws = C Eel. Eel
(35)
9
Since the stresses S22, S33 and S23 are neglected the constitutive matrix reads
E0 0
C = 0 G 0
(36)
0 0G .
Here, E and G denote Young's modulus and
shear modulus, respectively.
The plastic flow of metals can be described using the von Mises yield condition with
isotropic hardening
F (S, epl) = h(S)  y0 + r(epl)
(37)
where
h(S) = ST PS ,
100 P = 0 3 0
(38)
003 .
Assuming strain hardening the yield stress y is given here with the equivalent plastic
strain epl using a linear hardening function r(epl) = K epl, thus
y = y0 + K epl .
(39)
The initial yield stress y0 and the hardening parameter K are material constants. The rates of the plastic strains and of the equivalent plastic strain are described using an
associated flow rule
E pl = F
epl = F ,
(40)
S
r
where the gradient of the yield condition can be expressed as
F 1
F
= PS := N
= 1.
(41)
S h
r
In (40) the loadingunloading conditions
0 ,
F 0,
F = 0
(42)
must hold. Hence, a backward Euler
integration algorithm within a time step tn+1 = tn + t leads to
Epnl+1 = Epnl + Nn+1 epnl+1 = epnl +
(43)
where =
tn+1 tn
dt.
Inserting
(35)2
and
(43)1
into
(34)
at
time
tn+1
yields
Sn+1() = CЇ () Etr
Etr := En+1  Epnl
(44)
with
E
CЇ ()
=
(C1
+
y
P)1
=
1+E/y 0
0
0 G 1+3G/y 0
0
0
G
1+3G/y .
(45)
To express Sn+1 as an explicit function of we replaced h(S) by y() in (45) using F (S, epl) = 0. The consistency parameter is obtained with the solution of the yield
10
condition F (Sn+1, epnl+1) = 0. This is achieved iteratively within a Newton iteration pro
cedure
(l+1)
=
(l) 
F (l) dF (l)
d
F (l) = h Sn+1((l))  y epnl+1((l))
(46)
dF (l) d
=

h(1  y
y
dy depl
)
NT
CЇ N
+
dy depl
where l denotes the iteration number and dy/depl = K. As starting value we take (0) = 0. The equivalent plastic strain follows from (43)2. Note, that during the iteration y = h holds. The elastoplastic stresses at each integration point are evaluated using a
wellknown operator split method. The predictor step yields the socalled trial stresses
Str = C Etr. If above yield condition is not fulfilled the stresses are given with the
corrector step according to (44), thus
C Etr Sn+1 = CЇ Etr
if F (Str, epnl) 0 if F (Str, epnl) > 0 .
(47)
The consistency parameter depends on the strains E via eq. (46)2. This has to be considered when linearizing the stress vector (47). One obtains
dS
C
if F (Str, epnl) 0
dE
n+1
=
CЇ 
CЇ NNT CЇ NT CЇ N +
if F (Str, epnl) > 0
(48)
with = K/(1  K/y).
4.2 Elastic Stress Analysis
In case of small strains the transformation matrix F^ according to (16) is approximately the identity matrix, thus F^ 1. In this case A is given with (17)
1 0 0 0 3 2 0 wЇ
A = 0 1 0 3 0
0 wЇ,2 0
(49)
0 0 1 2 0 0 wЇ,3 0 .
One can see that A does not depend on the deformation and E = AE^ holds. For elasticity the stress vector is obtained from S = CAE^ . Now the matrix of the linearized stress resultants (30) can be reformulated using (48), (36), (49) and G^ = 0
E
D^ = A
0 G
0 0 G
0 G3 G2 G(22 + 32) sym
E3 0 0 0 E 32
E2 0 0 0 E23 E 22
0 GwЇ,2 GwЇ,3 G(2wЇ,3 3wЇ,2 ) 0 0 G(wЇ,22 +wЇ,23 )
EwЇ
0 0 0 EwЇ3 EwЇ2 0
dA .
(50)
EwЇ2
11
The centroid of the crosssection with coordinates {s2, s3} is denoted by S, see Fig. 1. Furthermore we denote the area of the crosssection by A, the moments of inertia relative to the centroid by IЇ22, IЇ33, IЇ23, the SaintVenant torsion modulus by IT and the warping constant by Iw~. The following definitions are given with , = 2, 3
S := dA = As
A IЇ := ЇЇ dA
I := dA = IЇ + Ass
A
A
I0 := I22 + I33
(51)
IT := [2(wЇ,3 +2)  3(wЇ,2 3)] dA = I0 + (2wЇ,3 3wЇ,2 ) dA
A
A
Iw~ := w~2 dA .
A
Using the equations of SaintVenants torsion theory (3) and application of Green's theorem yields after some algebra the missing quantities, see appendix A.1
wЇ,2 dA = As3
wЇ,3 dA = As2
A
A
(wЇ,22 +wЇ,23 ) dA = I0  IT
A
wЇ2 dA = IЇ22m3  IЇ23m2 := IwЇ2
wЇ3 dA = IЇ33m2 + IЇ23m3 := IwЇ3
A
A
wЇ2 dA = Iw~ + IЇ33m22 + IЇ22m23  2IЇ23m2m3 := IwЇ .
A (52)
A numerical computation of wЇ(2, 3) which fulfills A wЇ dA = 0 and above section quantities for arbitrary crosssections using the finite element method is discussed in e.g. [16].
In fact (52)4 and (52)5 are two equations which can be solved for the coordinates m2 and m3 if IwЇ2 and IwЇ3 are known, see Ref. [16]. Now, (50) can be expressed inserting the section quantities (51) and (52) as follows
EA D^ =
0 GA
0 0 GA
0 GAs3 GAs2 GI0 sym
EAs3 0 0 0 E I33
EAs2 0 0 0 EI23 E I22
0 GAs3 GAs2 G(I0  IT ) 0 0 G(I0  IT )
0
0 0 0 E IwЇ3 E IwЇ2 0
(53)
EIwЇ .
The elasticity matrix (53) is constant and symmetric. Hence, considering (30) the vector
of stress resultants follows from
 0
S^ = D^ E^
E^ =
 0
(54)
.
12
As can be seen torsionbending coupling occurs if the reference point O and the shear center M with coordinates {m2, m3} do not coincide. If the coordinates {2, 3} define principal axes of the crosssection and if S = M = O all offdiagonal terms ex cept G(I0  IT ) are zero. Furthermore, if the contribution of the bimoment and of the bishear to the strain energy is neglectable a formulation with six stress resultants can be derived, see.17 In this case one obtains the wellknown elasticity matrix D^ = diag [EA, GA, GA, GIT , EI33, EI22] .
5 FINITE ELEMENT FORMULATION
According to the isoparametric concept, the following kinematic variables are interpo lated using standard Lagrangian shape functions NI() where [1, +1]. Within a typical element the position vector of the reference curve S0h, the curve Sh of the current configuration and the parameter are interpolated by
nel Xh0 = NI ()XI , I =1
nel xh0 = NI ()(XI + uI ) , I =1
nel h = NI ()I . I =1
(55)
Here nel denotes the number of nodes at the element. For nel 3 the reference curve of
a spacecurved beam is approximated by
polynomial functions.
Furthermore the basis systems of the reference configuration and the current configuration
are approximated using the same interpolation functions
nel Ahm = NI () AmI , I =1
nel ahm = NI () amI . I =1
(56)
Thus the orthogonality condition of the basis systems Ahm and ahm is only fulfilled at the nodes. Numerical investigations however show that no significant loss of accuracy follows from this approximation. The initial basis system AmI is generated within the input of the finite element mesh whereas the current basis system at the finite element nodes is computed using the socalled Rodrigues formula
RI
=
amI
em
=
1
+
sin I I
I
+
1
 cos I I2
2I
(57)
with I = I. Note, that formula (57) is singularity free for 0 I < 2. The singularity at n2 (n = 1, 2, 3, ...) can be overcome by a multiplicative update of the
rotation tensor after a certain number of load steps. The skew
symmetric tensor I
follows from the independent rotational parameters by I = skew I. The basic equation reads I h = I Ч h for all h R3, thus the components are given by
1I I = 2I
I =
0 3I 2I
3I
0
1I
(58)
3I ,
2I 1I 0 .
The finite element approximations of the virtual displacements, basis vectors and the associated linearization are expressed as follows
nel
uh =
NI ()uI ,
I =1
nel
ahm =
NI ()amI ,
I =1
nel
h =
NI ()I ,
I =1
nel
ahm =
NI ()amI .
I =1
(59)
13
The derivative of the shape function NI() with respect to the arclength S is obtained using the chain rule NI() = NI(), /Xh0 , . Hence, the tangential vectors X0 and x0, the derivatives of the basis vectors Am and am and associated variations and linearizations are given replacing NI by NI in (55), (56) and (59). The variation of the orthogonal basis system amI yields for all h R3
amI = wI Ч amI = WmT I wI ,
(60)
h · amI = bmI (h) · wI , bmI (h) = amI Ч h
and nel h · am = NI bmI (h) · wI I =1
nel h · am = NI bmI (h) · wI . (61) I =1
Thus, the virtual beam strains E^ considering (19) can be expressed as follows
nel E^ h = BI vI , I =1
BI
=
NI RT 0 0
0
NI BTI NI BTI + NI BTI 0 0
0
0 NI
NI
(62)
with the virtual nodal displacement vector vI = [uI, wI, I]T and
R := [a1, a2, a3] BI := [b2I (a3), b3I (a1), b1I (a2)]
BI := [b1I (x0), b2I (x0), b3I (x0)] (63) BI := [b3I (a2), b1I (a3), b2I (a1)] .
Next the finite element interpolation (55) (63) is inserted into the linearized boundary value problem (27)
Anumel nel nel
L[g(vh, vh)] = e=1
vIT (fIe + KeIK vK ) . I=1 K=1
(64)
A Here, denotes the standard assembly operator, numel the total number of elements to discretize the problem and vK = [uK, wK, K]T the incremental displacement vector. Furthermore fIe and KeIK denote the sum of the
internal and external nodal forces of node I and the tangential stiffness matrix of element e related to nodes I and K,
respectively. Considering (62) one obtains
fIe = (BTI S^  NI qЇ) dS ,
KeIK = (BTI D^ BK + GIK + PIK ) dS
(65)
S
S
with
qЇ = [pЇ, mЇ I , pЇ1I ]T
mЇ I = rI Ч pЇ rI = p2 a2I + p3 a3I + wЇp a1I
(66)
pЇ1I = wЇp pЇ · a1I .
For linear elasticity the vector of stress resultants follows from the constitutive equation (54). For elastoplastic material behaviour S^ and D^ are obtained by numerical integration according to (22) and (30), respectively.
14
Next the geometric matrix GIK is derived with the linearized virtual strains (29). For this purpose the second variation of the current basis system is derived as
h · amI = wI · M(amI , h)wI
M(amI , h)
=
1 2 (amI
h
+
h
amI )
+
1 2 (tmI
I
+
I
tmI )
+
c10
1
(67)
where tmI and c10 are specified in appendix A.3. This yields
nel
h · am = NI wI · M(amI , h) wI
I =1 nel
(68)
h · am =
NI wI · M(amI , h) wI
I =1
for all h R3. As can be seen, the linearization of the basis system leads to a symmetric
bilinear form. One obtains
0
GIK = NI NK WfI
NI NK WfTK GwI Kw
0 0
(69)
0
0
0
,
where GwIKw = NI NK W^ I1K + NI NK W^ I2K + IK [M(a1I , h1I ) + M(a2I , h2I ) + M(a3I , h3I )] h1I = NI F 1x0 + NI M 2a3 + NI M 3a2 h2I = NI F 2x0 + NI M 3a1 + NI M 1a3 h3I = NI F 3x0 + NI M 1a2 + NI M 2a1 W^ I1K = M 1W2I W3TK + M 2W3I W1TK + M 3W1I W2TK W^ I2K = M 1W3I W2TK + M 2W1I W2TK + M 3W2I W3TK WfI = skew (F 1a1I + F 2a2I + F 3a3I ) . (70)
The external loading yields a contribution to the stiffness matrix if rI = 0. Considering
(33) one obtains
0
0
0
PIK =  0 NI IK M(rI , pЇ) NI NK mЇ 1I
(71)
0 NI NK mЇ T1K
0
where mЇ 1I = wЇp a1I Ч pЇ. Application of the chain rule yields the differential arclength dS = Xh0,  d. The inte
gration of the element residual and element stiffness matrix with respect to the coordinate
S is performed numerically. To avoid shear locking uniform reduced integration is applied
to all quantities. Finally the transformations of the axial vectors are introduced as follows
wI = HI I
wK = HK K
(72)
where the tensor H is specified in appendix A.2. This leads to
L[g(vh, vh)] v~I
Anumel nel nel
=
e=1
v~IT (~fIe + K~ eIK v~K ) I=1 K=1
= [uI , I , I ]T v~K = [uK, K, K]T
(73)
~fIe = TTI fIe
K~ eIK = TTI KeIK TK
TI = diag[1, HI, 1] .
15
From eq. (73) one obtains a linear system of equations for the incremental nodal degrees of freedom. It is emphasized, that the update of the nodal displacements as well as of the rotational parameters is performed in an additive way. Thus, the
equilibrium configuration is computed iteratively within the Newton iteration procedure.
Remark: Within an alternative procedure one keeps wI as unknown incremental rotational parameters. Thus the transformation (73) is not necessary anymore. However I must be known for the finite element interpolation (57) (71). It can be obtained within the following simple update procedure
(Kn+1) = (Kn) + K
K = HK1wK
HK1
=
1

1 2
(Kn)
+
1 2
c3
2K(n)
(74)
where n denotes the index of the Newton iteration procedure. One can easily show that above expression fulfills HK1HK = 1. Again, the advantage of (74) compared with (73) is, that the transformation wI = HI I drops out. However, the algorithm (74) requires additional storage of the parameters (Kn).
6 NUMERICAL EXAMPLES The element scheme has been implemented in an enhanced version of the program FEAP documented in Ref. [18]. In this section three examples with finite deformations and elasticplastic material behaviour are presented. With the first example we investigate the stability behaviour of a beam assuming linear elastic material behaviour. The last two examples are concerned with the load carrying capacity of beam structures. The second example is a channelsection beam, where centroid, center of shear and loading point are not identical. With the last example we solve a coupled beamshell problem. For comparison the investigated thinwalled beam structures are discretized using four noded shell elements. These elements possess six nodal degrees of freedom which are identical to the first six degrees of freedom of the beam.
6.1 Lateral Torsional Buckling of a Single Span Girder The first example is a single span girder with length L = 150 cm where an
axial force F is applied at the centroid. Fig. 2 shows the crosssection of the channel section beam. The ratio of the width of the flange to the height of the web is relatively large. This type of crosssection is fairly sensitive against torsional buckling. The geometrical and material data are given in (75), where the moments of inertia are denoted according to the definitions in (51). The section quantities are determined with the finite element program as is described in [16].
16
x3
M h
S
x2
s
t
b m2 Figure 2: Crosssection of a single span girder
b= h= s= t= m2 =
10.0 cm 10.0 cm 0.2 cm 0.2 cm 7.55 cm
A = 5.92 cm2
I33 = 110.8 cm4 I22 = 64.49 cm4
IT = 0.0792 cm4 Iw~ = 1108.2 cm6
(75)
E = 21000 kN/cm2 G = 8077 kN/cm2
The
boundary conditions are chosen as follows: The displacements at both supports
except the axial displacement at the loaded crosssection are fixed. Furthermore the
torsion angle is fixed at both ends. The other degrees of freedom are not restrained. The
theoretical solution based on a secondorder geometrical nonlinear Bernoulli beam theory including torsionwarping deformation may be found e.g. in Kollbrunner and Meister.19
For this purpose we compute the following reference quantities with n = 1 for the lowest
eigenvalue
F1
:=
n2
2 E I22 L2
=
594.0 kN
F2
:=
n2
2 E I33 L2
=
1020.0 kN
(76)
F3
:=
1 i2M
(GIT
+
n2
2 E Iw~ L2
)
=
125.1 kN
.
Here, i2M = i2P + m22 follows with i2P = (I22 + I33)/A. Furthermore F1 describes the critical load for buckling in the symmetry plane. The theoretical critical load is given according
to
Fcr = 2
11 + F2 F3
+
11
2 +
4
F2 F3
F2 F3
1
m2
2
= 115.4 kN .
iM
(77)
17
As the result shows Fcr is close to the pure torsional buckling load F3. Thus one obtains a significant reduction of the critical load perpendicular to the symmetry plane due to
torsional buckling.
For the numerical solution we discretize the beam with 4, 8 and 16 threenoded beam
elements. The boundary conditions are set as described above. We solve the eigenvalue
problem
(KT  21) = 0
(78)
where the contribution of the element nodes to the tangential stiffness matrix KT are given in (65). The critical loads Fcr are characterized by zeroeigenvalues i with associated eigenvectors i. The
minimum values for the different discretizations are given in table 1. As can be seen the agreement of our results with the solution of the secondorder geometrical nonlinear Bernoulli theory is very good.
Table 1: Critical loads of the axially compressed single span girder
numel 4 8 16 analytical
Fcr in kN 115.3 115.0 115.0 115.4
18
6.2 Channelsection beam A channelsection beam clamped at one end and subjected to a tip force at the free end is investigated next, see Fig. 3. We assume linear elastic and ideal plastic material behaviour with constants according to Fig. 3. The yield stress is y0 = 36 kN/cm2. The developed beam model is compared with a shell model. The discretization is performed with 30 twonoded
beam elements and in the second case with 360 fournoded shell elements. The shell discretization consists of 36 elements along the length direction, 6 elements along the web and 2 elements for each flange. In the following the vertical displacement w of point O at the cantilever tip is computed. For the elasticplastic case we load up to a tip displacement w = 250 cm and then unload the structure. The results for both models agree very good in the total range of the computed load deflection curve, see Fig. 4. This holds for elastic as well as for inelastic material behaviour. Fig. 5 shows a plot of the von Mises stresses according to eq. (38) for the ultimate state and the unloaded state. In Fig. 6 the stress distribution of the crosssection at a distance of 225 cm from the clamped end is given again for the ultimate state and the unloaded state. In the first case the crosssection is completely plastified, whereas in the second case
residual stresses can be seen. P
P L
L = 900 cm
h = 30 cm
h
0
b = 10 cm
s = 1:0 cm
s t
t = 1:6 cm
E = 21000 kN=cm2
= 0:30
b
Figure 3: Channelsection beam with geometrical and material data
19
load P in kN
20 18 16 14 12 10 8 6 4 2 0 0
elastic
beam shell
plastic
25 50 75 100 125 150 175 200 225 250 vertical displacement w in cm
Figure 4: Load deflection curves of the channelsection beam
0.000E+00 min 6.000E+00 1.200E+01 1.800E+01 2.400E+01 3.000E+01 3.600E+01 max Figure 5: Von Mises stresses in kN/cm2 for the ultimate state and unloaded state 20
0.000E+00 min 6.000E+00 1.200E+01 1.800E+01 2.400E+01 3.000E+01 3.600E+01 max Figure 6: Von Mises stresses in kN/cm2 of the crosssection at a distance of 225 cm from the clamped end for the ultimate state and unloaded state 21
6.3 Stiffened plate
With the last example coupling of the beam element with shell elements is investigated. Fig. 7 shows the crosssection of a rectangular plate of length and width b with an Lshaped stiffener. The structure is loaded by a single force, see Fig. 7. The geometrical and material data are given as follows
= 300 cm
E = 21000 kN/cm2
b = 80 cm t = 0.5 cm
= y0 =
0.3 36 kN/cm2
(79)
K=
0
.
For comparison we apply two types of different discretizations. In both cases the plate is modeled using 20 fournoded shell elements in length direction and 10 elements in transverse direction. In the first case the L shaped stiffener is discretized with 20 shell elements in length direction, 3 elements for the vertical leg and 2 elements for the horizontal leg. In the second case 20 twonoded beam elements are used. A computation without stiffener shows that the beam leads to a significant reduction of the deformation. The vertical displacement of the loading point is depicted in Fig. 8 for elastic and elasticplastic material behaviour. Both models are in good agreement in large ranges of the computed load deflection curve. Plots of the von Mises stresses and of the equivalent plastic strains are depicted in Fig. 9, 10 and 11. The maximum plastic strain is about 5%. Fig. 11 shows that the structure undergoes large deformations.
l
P
t
L 5.0 x 4.0 x 0.5
b/2
b/2
Figure 7: Clamped rectangular plate with stiffener
22
load in kN
16
14
12
elastic
10
shell  beam
shell  shell
8
6
4 plastic 2
0
0
50
100
150
200
vertical deflection of the loading point in cm
Figure 8: Load deflection curves of the stiffened plate
.000E+00 min 6.000E+00 1.200E+01 1.800E+01 2.400E+01 3.000E+01 3.600E+01 max Figure 9: Von Mises stresses in kN/cm2 for the stiffened plate at w = 25 cm 23
.000E+00 min 6.000E+00 1.200E+01 1.800E+01 2.400E+01 3.000E+01 3.600E+01 max Figure 10: Von Mises stresses in kN/cm2 for the stiffener at the clamped crosssection at w = 25 cm .000E+00 min 1.000E02 2.000E02 3.000E02 4.000E02 5.029E02 max Figure 11: Equivalent plastic strains of the stiffened plate 24
7 CONCLUSIONS In this paper a theory for threedimensional beams with arbitrary crosssections and an associated finite element formulation is developed. The beam strains are derived from the GreenLagrangian strain tensor. Elastoplastic material behaviour applying the von Mises yield condition and associated flow rule is considered. Due to the nonlinear stressstrain relations the stress resultants and associated linearizations are obtained by numerical integration over the crosssections. A finite beam element is developed using Lagrangian interpolation functions to approximate the kinematic quantities. Each node of the element possesses seven degrees of freedom. The basis systems at the nodes are evaluated using orthogonal transformations. Due to the chosen finite element approximation the linearization yields relative simple expressions for the symmetric tangent matrix. Examples show the applicability of the developed beam element applied to geometrical and physical nonlinear problems. Alternative discretizations of thinwalled crosssections with shell elements show good agreement between the different models. Thus, the derived element can effectively be used to analyze the loadcarrying capacities of spatial beam structures. 25
A APPENDIX A.1 Proof of integrals (52) Using (3)1 the following integral is reformulated such that Green's formula can be applied. Hence inserting boundary condition (3)2 yields
wЇ,2 dA = [(2wЇ,2 ),2 +(2wЇ,3 ),3 ] dA = [2(wЇ,2 n2 + wЇ,3 n3)] dC
A
A
C
(80)
= [2(3n2  2n3)] dC .
C
Again application of Green's formula and considering (51)1 leads to
[2(3n2  2n3)] dC = [(23),2 22,3 ] dA = 3 dA = A s3
(81)
C
A
A
which proves (52)1. Proceeding in an analogous way leads to (52)2.
The following integral is reformulated in the same way as done above. Thus, we apply (3)1, Green's formula and boundary condition (3)2
(wЇ,22 +wЇ,23 ) dA = [(wЇ wЇ,2 ),2 +(wЇ wЇ,3 ),3 ] dA
A
A
(82)
= [wЇ(wЇ,2 n2 + wЇ,3 n3)] dC = wЇ (3n2  2n3) dC .
C
C
Again application of Green's formula and considering the definition of SaintVenant torsion modulus (51)4 yields
wЇ(3n2  2n3) dC = (wЇ,2 3  wЇ,3 2) dA = I0  IT
(83)
C
A
which proves (52)3.
Next using (4) and (5) we get
wЇ2 dA = (Ї2 + s2)(w~  m2Ї3 + m3Ї2) dA = IЇ22m3  IЇ23m2 = IwЇ2
(84)
A
A
which yields (52)4. The expression for IwЇ3 is obtained in an analogous way. Finally eq. (52)6 can be derived considering the orthogonality conditions (4) and definition (5).
A.2 First Variation of the Orthogonal Basis System
The current orthogonal basis system can be written using the Rodrigues formula (57)
ai = R ei
R
=
1
+
sin
+
1
 cos 2
2
(85)
where =  and = skew . To alleviate the notation the node index is omitted. The first variation, denoted here by the symbol , yields
ai = Rei = RRT ai = w Ч ai
(86)
26
with
sin
1  cos
R =
+
2
( + )
+
cos  2
sin
+
sin
+ 2 cos 3

2 2
(87)
RT
=
1

sin
+
1
 cos 2
2
with = ( · )/. Inserting the identities
2 =  21 3 = 2 =  ( · )1
(88)
yields after some lengthy algebraic manipulation the skewsymmetric tensor
RRT = (1  c22) + c1 (  ) + c2 ( · )
1  cos
 sin
(89)
c1 =
2 ,
c2 = 3 .
The associated axial vector reads
w = H ,
H = 1 + c1 + c2 2 .
(90)
A.3 Second Variation of the Orthogonal Basis System
The second variation of the orthogonal basis system is denoted by . One obtains for all h R3 h · ai = h · (w Ч ai)
= h · [w Ч (w Ч ai) + H Ч ai]
(91)
= w · [ai h  (ai · h)1] w + bi · H
where bi = ai Ч h. The second part in (91) yields with = ( · )/
bi · H
=
bi
·
[( c1
+
c2
2)
+
c1
+
c2(
+
)]
= · {[(c4 + c52)bi ] + c1bi Ч + c2( + )bi} (92)
with
sin + 2 cos  2
c3 =
2(cos  1)
c4
=
1
c1
=
c1 c3
,
c5
=
1
c2
=
c6

c2 c3
,
c6
=
c3  c2 2
.
(93)
Considering
bi = (bi Ч ) Ч
bi = [(bi · )1  bi ]
Di := c1Bi + c2Ci
(94)
Bi := h ai  ai h
Ci := bi  bi
si := (c2 1  c4 + c5 2) bi
27
we get
bi · H = · [si + Di + c2(bi · )1] .
(95)
Inserting this result into (91) yields
h · ai = w · Mi w (96) Mi = ai h  (ai · h)1 + HT 1 [si + Di + c2(bi · )1] H1 .
Within a multiplicative update procedure as is applied in [7] or [8] the last term in (96) is missing and one obtains a nonsymmetric tangent operator. It follows if w = is chosen, thus H = 1. However, this is justified for a multiplicative procedure with unknown rotation increments since H() approaches the unit tensor at an equilibrium configuration. The symmetry of Mi is shown in the following. For this purpose Mi is split in a symmetric and a skewsymmetric part, thus Mi = MSi +MAi where we show that the skewsymmetric part cancels out. We introduce
wi := HT 1si = c3bi + c6 (bi · ) HT 1 =
(97)
where the first eq. follows immediately by multiplying (97)1 with HT , comparing the coefficients and considering (93). The second equation is evident with (74) and = 0. One obtains
MAi
=
1 2
(Mi

MTi
)
=
1 2 (ai
h

h
ai)

1 2 c3(bi

bi)
+
HT 1
Di
H1
(98)
=

1 2
Bi
+
1 2 c3Ci
+
HT 1
Di
H1
.
Considering (74) we get
HT 1 Di H1
=
Di
+
1 2
(Di

Di)
+
1 2
c3(2Di
+
Di2)
+
1 4
c3(Di2

2Di)

1 4
Di
+
1 4
c232Di
2
with
Di = di  (di · )1 Di = (di · ) Di2 = 2Di = (di · ) 2 2Di + Di2 = (di · )  2Di 2Di2 = 2(di · )
thus
HT 1
Di
H1
=
c7
Di
+
1 2
(Di

Di)
+
c8(di
·
)
c7
=
1

1 2
c32
c8
=
1 4

1 2
c3
+
1 4
c23
2
with di = c1 bi + c2 bi Ч . Using di · = c1bi · , (94)3 and
1 2 (Di

Di)
=

1 2
c2(bi
·
)
+
1 2
c22
Bi

1 2 c1
Ci
(99) (100) (101) (102)
28
one obtains the final version of the skewsymmetric part as
MAi
=
(c1c8

1 2 c2)
(bi
·
)
+
( 1 2
+
c1c7
+
1 2
c22)
Bi
+
1 ( 2 c3
+
c2c7

1 2 c1)
Ci
,
(103)
where one can easily show that the coefficients vanish, and thus MAi 0. The symmetric part of Mi reads
MSi
=
1 2
(Mi + MTi )
=
1 2
(ai
h
+
h
ai)

(ai
·
h)1
+
1 2
(wi
+
wi)
+
c2(bi
·
)
HT 1
H1
(104)
and with H1 according to (74) and (88)
HT 1H1 = 1 + c92
c9
=
1 4
+
c3

1 4
2
c23
(105)
we may express (104) as follows
MSi
=
1 2
(ai
h
+
h
ai)
+
1 2
(ti
+
ti)
+
c101
ti = c3 bi + (c6 + c2c9) (bi · )
c10 = c2(1  c9 2) (bi · )  (ai · h) .
(106)
Applying a multiplicative update procedure the rotational parameters are replaced by which vanish in the Newton iteration process. Thus all terms in in eq. (106) cancel out. This again shows that M becomes symmetric at an equilibrium configuration within the multiplicative procedure.7
29
References 1. J.H. Argyris, O. Hilpert, G.A. Malejannakis and D.W. Scharpf, On the geometrical stiffness of a beam in space A consistent v.w. approach, Comp. Meth. Appl. Mech. Engrg. 20 (1979) 105131. 2. K.J. Bathe and S. Bolourchi, Large displacement analysis of threedimensional beam structures, Int. J. Num. Meth. Engng. 14 (1979) 961986. 3. T. Belytschko and B.J. Hsieh, Nonlinear finite element analysis with convected coordinates, Int. J. Num. Meth. Engng. 7 (1973) 255271. 4. M.A. Crisfield, A consistent corotational formulation for nonlinear three dimensional beam elements, Comp. Meth. Appl. Mech. Engrg. 81 (1990) 131150. 5. B. NourOmid and C.C. Rankin, Finite rotation analysis and consistent linearization using projectors, Comp. Meth. Appl. Mech. Engrg. 93 (1991) 353384. 6. E. Reissner, On finite deformations of spacecurved beams,
J. Appl. Math. Phys. (ZAMP) 32 (1981) 734744. 7. J.C. Simo and L. VuQuoc, A threedimensional finitestrain rod model. Part II: Computational aspects, Comp. Meth. Appl. Mech. Engrg. 58(1) (1986) 79116. 8. J.C. Simo and L. VuQuoc, A geometricallyexact rod model incorporating shear and torsionwarping deformation, Int. J. Solids Structures 27(3) (1991) 371393. 9. A. Cardona and M. Gґeradin, A beam finite element nonlinear theory with finite rotations, Int. J. Num. Meth. Engng. 26 (1988) 24032438. 10. A. Ibrahimbegoviґc, Computational aspects of vectorlike parametrization of three dimensional finite rotations, Int. J. Num. Meth. Engng. 38 (1995) 36533673. 11. E. Reissner, Some considerations on the problem of torsion and flexure of prismatical beams, Int. J. Solids Structures 15 (1979) 4153. 12. E. Reissner, On a simple variational analysis of small finite deformations of prismatical beams, J. Appl. Math. Phys. (ZAMP) 34 (1983) 642648. 13. M. Gґeradin and D. Rixen, Parametrization of finite rotations in computational dynamics: a review, Revue europґeenne des ґelґements finis, 4 (1995) 497553 14. P. Betsch, A. Menzel, E. Stein, On the parametrization of finite rotations in computational mechanics A classification of concepts with application to smooth shells, Comp. Meth. Appl. Mech. Engrg. 155 (1998) 273305. 15. S.P. Timoshenko and J.N. Goodier, Theory of Elasticity,
3rd Edition, McGrawHill International Book Company, 1984. 16. F. Gruttmann, W. Wagner and R. Sauer, Zur Berechnung von WoЁlbfunktion und Torsionskennwerten beliebiger Stabquerschnitte mit der Methode der finiten Elemente, Bauingenieur 73(3) (1998) 138143. 30
17. F. Gruttmann, R. Sauer and W. Wagner W, A geometrical nonlinear eccentric 3D beam element with arbitrary crosssections, Comp. Meth. Appl. Mech. Engrg. 160 (1998) 383400. 18. O.C. Zienkiewicz and R.L. Taylor, The Finite Element Method,
4th Edition, Volume 1 (
McGraw Hill, London, 1988). 19. C.F. Kollbrunner and M. Meister, Knicken, Biegedrillknicken, Kippen : Theorie und Berechnung von KnickstЁaben; Knickvorschriften,
2nd Edition (Springer, Berlin, 1961). 31
F Gruttmann, R Sauer