Sensitivity computations for a geometrically nonlinear 3D beam model.

Pierre Aubert^{∗}  Bernard Rousselet^{†}

June 6, 2000

Abstract

A geometrically nonlinear three dimensional initialy curved beam model including flexion, shear and torsion is described and studied for sensitivity analysis. The sensitivity with respect to the thickness of the cross-section and to the mid-line of the beam is presented.

Key words: sensitivity analysis, geometrically nonlinear curved beam.

AMS subject classifications: 73 K 40, 49 Q 10.

1 Introduction

In [1], we adressed the thicknes and shape sensitivity of an initially curved plane beam; the finite displacement model used has been proposed by Fried [2]; it does not include any transverse shear measure. Here we consider the model proposed by Reissner [3] Antman and Kenney [4] which includes transverse shear; it is three dimensional. Numerical implementation of this model has been adressed by several authors: Simo and Vu-Quoc [5], Bourgat, Le Tallec & Mani [6], Ibrahim, Frey & Kozar [7]. The shape sensitivity has been considered in [8] and for small loads, proved in [9]. Here we consider a material version of the equations in order to make easier the implementaton of the shape sensitivity; we use a three parameter parametrization of the rotations by the so called rotation vector through the use of the Rodrigues formula as proposed by Ibrahim, Frey & Kozar [7]. Implementation of thickness and shape sensitivity has been realized with the use of a formal computation system; numerical results will be presented and discussed. The case of sensitivity with instabilities is not considered here (see [AR98] for the case of optimization with instabilities of the plane model proposed by Fried [2]).

2 A three dimensional beam model with transverse shear

2.1 Description of a 3D beam

Following many authors[4, 10, 11], a beam deformed configuration may be described by two applications ϕ and Γ; ϕ describes a curve in ℝ^{3} which is the line of centroids of the beam and Γ is an orthogonal frame attached to the cross-section of the beam. The underlying mechanical hypothesis is that the cross-section remains plane but its normal is not necessarily colinear with ϕ′. More precisely, let (e⃗_{1}, e⃗_{2}, e⃗_{3}) be the Euclidean reference frame, we consider ϕ_{0} : [0, L] → ℝ^{3}; S ↦ ϕ_{0}(S) = (x_{0}, y_{0}, z_{0}) and ϕ : [0, L] → ℝ^{3}; S ↦ ϕ(S) = (x, y, z). where ϕ_{0} is a map which describes the line of centroids of the undeformed configuration and ϕ describes the line of centroids of the deformed (or loaded) configuration; S is the arc lenth of the line of centroids of the undeformed configuration. If the line of centroids of the undeformed configuration is defined by a map, we have

∗maply, umr n^{o}5585, batiment 401, 20, avenue Albert Einstein, insa, Villeurbanne, France, pierre.aubert@insa-lyon.fr.

†Laboratoire J.-A. Dieudonné, umr n^{o}6621, Université de Nice-Sophia Antipolis, F-06108 Nice Cedex 2, br@math.unice.fr.

This space left blank for copyright notice.

S(ξ) = ∫_{0}^{ξ} ∥ϕ′_{0}(t)∥dt  and  L = ∫_{0}^{1} ∥ϕ′_{0}(t)∥dt. (2.1)

The physical properties of the beam are given by L, the length of the beam and e, b : [0, L] → ℝ_{+}^{∗} two functions that describe the thickness and the width of the rectangular cross-section of the beam.

We call E⃗_{1} the unit vector normal to the cross-section of the undeformed beam, and we complete it with E⃗_{2} and E⃗_{3} (co-linear to the axes of principal inertia) such that (E⃗_{1}, E⃗_{2}, E⃗_{3}) is an orthonormal basis of the undeformed configuration of the beam called the material basis. It should be noted that E⃗_{1} and ϕ′ = \frac{dϕ_{0}}{dS} are not necessarily co-linear: the unloaded beam may already have transverse shear or torsion produced for example by a pre-stress.

In the same way, we define an orthonormal basis (t⃗_{1}, t⃗_{2}, t⃗_{3}) attached to the deformed configuration of the beam such that t⃗_{1} is normal to the cross-section of the beam and t⃗_{2}, t⃗_{3} co-linear to the the axis of principal inertia; this basis is called the spatial basis. These basis define two rotation transformations such that : Γ_{0} : [0, L] → SO_{3}; S ↦ Γ_{0}(S) and Γ : [0, L] → SO_{3}; S ↦ Γ(S) where SO_{3} is the special orthogonal group and

∀I = 1, 2, 3,  E⃗_{I} = Γ_{0}e⃗_{I},  t⃗_{I} = ΓE⃗_{I}. (2.2)

We introduce ω̃_{0}, ω̃, κ̃_{0} and κ̃ the skew-symmetric tensors defined by :

(2.3)  ω̃_{0} = Γ′_{0}Γ_{0}^{T} ,  ω̃ = Γ′Γ^{T} ,  κ̃_{0} = Γ_{0}^{T} Γ′_{0},  κ̃ = Γ^{T} Γ′.

By differentiating the material and spatial basis with respect to S, we get :

E⃗′_{i} = ω̃_{0}E⃗_{i} and t⃗′_{i} = (ΓΓ_{0})′e⃗_{i} = (Γ′Γ_{0} + ΓΓ′_{0})(Γ_{0}^{T} Γ^{T})t⃗_{i} = (ω̃ + Γω̃_{0}Γ^{T})t⃗_{i}.

An example of a set of possible configuration of the rod is given by :

C_{1} = {(ϕ, Γ) ∈ (H^{1}([0, L]))^{3} × L([0, L], SO_{3}), ϕ(0) = ϕ_{0}(0), Γ(0) = Γ_{0}(0)} . (2.4)

The set describes a rod which is clamped at one end and free at the other; other cases have been dealt with by [12, 4].

2.2 Parametrisation of finite rotation

The sets C_{1} is a nonlinear differentiable manifold. This nonlinearity introduces some complex behaviour[5], which depends on the parametrisation of finite rotations. In particular, we should use tangent spaces which are useless in the case of a vector space. As discussed in [13], if we use quaternions in order to parametrise finite rotations and a linear interpolation of variables, the Hessian of the energy function may not be symmetric and a quadratic rate of convergence for the Newton-Raphson process used in the resolution is difficult to achieve. If we use quaternions and a nonlinear interpolation scheme, we recover the quadratic rate of convergence but this leads to more difficult implementation and a special update procedure is required.

An alternative approach has been proposed by [14]: an element of SO_{3} has 9 entities but due to the relation ΓΓ^{T} = I_{3} there are only 3 independent parameters. Following [14], we parametrise the two-point tensor Γ with its axial vector θ⃗. Denoting θ = ∥θ⃗∥, the Rodrigues formula gives [15] :

(2.5) Γ(θ⃗) = exp(θ̃) = cos θI_{3} + \frac{sin θ}{θ} θ̃ + \frac{1 − cos θ}{θ^{2}} θ⃗ ⊗ θ⃗.

The advantage of this parametrization is that we work in ℝ^{3}; the trouble is that this parametrization is singular for θ = 2kπ. Let δψ⃗ be an infinitesimal rotation, i.e. the axial vector associated with ψ̃ ∈ so_{3}. Due to the exponential mapping (see [16]), the Gateaux derivative of (2.5) in direction δψ̃ reads :

(2.6) δΓ[δψ̃] = lim_{η→0} [\frac{d}{dη} Γ(θ⃗) exp(ηδψ̃)] = Γ(θ⃗)δψ̃.

As we want to use the parameter θ⃗, we need to relate the variation δ⃗ψ to the variation of θ⃗. We rewrite Γ(θ⃗ + δθ⃗) = exp(θ̃+ ηδθ̃) thusexp(ηδψ̃) = exp(−θ̃) exp(θ̃+ ηδθ̃+ o(η)). By differentiation, and setting η = 0, we get after a lenthy calculation (see [7]) :

δψ̃ = T_{m}(θ̃)δθ̃ with T_{m}(θ̃)δθ̃ = c_{1}δθ̃ − c_{2}[θ̃, δθ̃] − 2c_{3}tr(θ̃δθ̃)θ̃ (2.7)

where c_{1}(θ) = (sin θ)/θ, c_{2}(θ) = (1 − cos θ)/θ^{2} and c_{3}(θ) = (θ − sin θ)/θ^{3}. An equivalent relation in the other basis is given by δψ⃗ = T_{m}(θ⃗)δθ⃗ with T_{m}(θ⃗) = c_{1} I_{3} − c_{2} θ̃ + c_{3}θ⃗ ⊗ θ⃗.

Lemme 2.1. The variation of a rotation matrix may be expressed as :

δΓ(θ̃)[δθ̃] = Γδψ̃ = Γ(θ̃)T_{m}(θ̃)δθ̃ = Γ(θ⃗)({ \atop T_{m}(θ⃗)δθ⃗}. (2.8)

2.3 Definition of strains

The following expressions of spatial strain measures ϵ^{ϕ} and κ^{ϕ} are introduced by [10] :

ϵ⃗^{ϕ}(ϕ, Γ) = ϕ′ − t⃗_{1} and κ̃^{ϕ}(ϕ, Γ) = ω̃ + Γω̃_{0}Γ^{T} .

As we have in mind the differentiation (see §3.3) with respect to the shape of reference configuration we shall use the expressions of material strain measures (2.3) (denoted by bold symbol) :

(2.9)  ϵ⃗(ϕ, Γ) = Γ^{T} ϵ⃗^{ϕ} = Γ^{T} ϕ′ − E⃗_{1}
(2.10)  κ̃(ϕ, Γ) = Γ^{T} κ̃^{ϕ}Γ = Γ^{T} ω̃Γ + ω̃_{0} = Γ^{T} Γ′ + ω̃_{0} or κ⃗(ϕ, Γ) = Γ^{T} ω⃗ + ω⃗_{0}.

If we parametrise the orthogonal tensor Γ resp (Γ_{0}) with its rotational vector θ⃗ (resp. θ⃗_{0}) then strain expressions (2.9) and (2.10) become by using (2.3) and (2.7) : ϵ(ϕ, θ⃗) = Γ^{T} ϕ′−E⃗_{1} and κ̃(ϕ, θ) = T_{m}(θ̃)θ̃′+ω̃_{0} = T_{m}(θ̃)θ̃′+Γ_{0}T_{m}(θ̃_{0})θ̃′_{0}Γ_{0}^{T} or κ⃗(ϕ, θ⃗) = T_{m}(θ⃗)θ⃗′ + ω⃗_{0} = T_{m}(θ⃗)θ⃗′ + Γ_{0}T_{m}(θ⃗_{0})θ⃗′_{0}. and their virtual variations can be formally computed by : δϵ(ϕ, θ⃗)[δϕ, δθ⃗] = Γ^{T} δϕ′ + (Γ^{T} ϕ′) × T_{m}(θ⃗)δθ⃗ and δκ⃗(ϕ, θ⃗)[δϕ, δθ⃗] = δT_{m}(θ⃗)[δθ⃗]θ⃗′ + T_{m}(θ⃗)δθ⃗′. As shown by [17] this definition of strains amounts to elimination of rigid body motions.

2.4 Definition of stress

We assume an hyper-elastic material behaviour defined with a strain energy density function w and to make things easier we assume a Saint-Venant, Kirchoff class. For simplicity, the material elasticity tensor is taken diagonal and written in two parts, C_{M}(e, b) = diag(EA, GA_{2}, GA_{3}) and C_{R} = diag(GJ, EI_{2}, EI_{3}) where E is the Young modulus, G is the shear modulus, A, A_{2}, A_{3}, I, I_{2}, I_{3}, J are mechanical properties of the cross-section which depend only of the shape of the cross-section[18]. Expressions for a square cross-section are given in section III.

The material description of internal stress resultant N⃗ and the internal stress couple M⃗ are related to the strain measures : N⃗(ϕ, θ⃗) = C_{M}(e, b)ϵ(ϕ, θ⃗) and M⃗(ϕ, θ⃗) = C_{R}(e, b)κ(ϕ, θ⃗). The strain density energy w is then given by :

w(ϕ, θ⃗) = \frac{1}{2}ϵ(ϕ, θ⃗)^{T} N(ϕ, θ⃗) + \frac{1}{2}κ(ϕ, θ⃗)^{T} M(ϕ, θ⃗) = \frac{1}{2}ϵ(ϕ, θ⃗)^{T} C_{M}(e, b)ϵ(ϕ, θ⃗) + \frac{1}{2}κ(ϕ, θ⃗)^{T} C_{R}(e, b)κ(ϕ, θ⃗) (2.11)

and the potential energy π is given by

π(ϕ, θ⃗) = w(ϕ, θ⃗) − (N̄(ϕ, θ⃗)δϕ + M̄(ϕ, θ⃗)T_{m}(θ⃗)δθ⃗) (2.12)

where N̄ and M̄ are the density of external forces and moment applied to the beam. The total potential energy is then given by integration of π along the line of centroids : Π(ϕ, θ) = ∫_{0}^{L} π(ϕ, θ⃗)dS.

2.5 Equilibrium equation

As we have parametrised Γ by its axial vector, the spaces defined by (2.4) become :

C_{1} = {(ϕ, θ⃗) ∈ (H^{1}([0, L]))^{3} × (H_{1}([0, L], B(0, π))^{3}, ϕ(0) = ϕ_{0}(0), θ⃗(0) = θ⃗_{0}(0)} .

H^{1}([0, L]) is a Sobolev space [19] and B(0, π) the ball of center 0 and radius π in ℝ^{3}. We note that δψ̃ is an element of so_{3} but θ⃗ is an element of B(0, π)^{3} and δθ⃗ ∈ ℝ^{3}. We restrict the norm of the rotation because the matrix T_{m} is not invertible at θ = 2kπ, k ≥ 1.

These spaces are linear manifolds contrary to those defined by (2.4). The tangent spaces are respectively given by [6] : TC_{1} = {(δϕ, δθ⃗) ∈ (H^{1}([0, L]))^{3} × H_{1}([0, L])^{3}, δϕ(0) = 0, δθ⃗(0) = 0}. The principle of virtual work may be written : find (ϕ, θ⃗) ∈ C_{1,2,3} such that ∀(δϕ, δθ⃗) ∈ TC_{1,2,3}

G(ϕ, θ⃗)[δϕ, δθ⃗] = 0 (2.13) with

G(ϕ, θ⃗)[δϕ, δθ⃗] = ∫_{0}^{L} δϵ(ϕ, θ⃗)^{T} C_{M}(e, b)ϵ(ϕ, θ⃗) dS + ∫_{0}^{L} δκ(ϕ, θ⃗)^{T} C_{R}(e, b)κ(ϕ, θ⃗) dS
− ∫_{0}^{L} (N̄(ϕ, θ⃗)δϕ + M̄(ϕ, θ⃗)T_{m}(θ⃗δθ⃗) dS.

(2.14)

Theorem 2.1. Provided that the external forces do not depend on ϕ and θ (dead loading), the problem (2.13) has at least one solution. A proof is given in [6].

3 Sensitivity analysis

3.1 Principle

The weak form of the balance equations are given by (2.14). We rewrite them as : Let F be the space of all admissible shapes, let η_{0} ∈ F and let u ∈ C_{1,2,3} such that : ∀δv ∈ TC_{1,2,3}, G(η_{0}; u)[δv] = 0 where the design variable η_{0} has been added in the variables of G. If we consider a thickness optimisation problem, η_{0} will be a parametrisation of the cross-section of the beam for example: η_{0} = (e, b). If we investigate shape optimisation then η_{0} will be a parametrisation of the unloaded configuration of the beam, η_{0} = (ϕ_{0}, θ⃗_{0}). As mentioned in the previous section, the unloaded beam may already have transverse shear or torsion produced for example by a prestress; the design variable should take this fact into account; this is in contrast with [8] where ϕ_{0} is the only design variable; E⃗_{i} is the Frénet frame associated to the line of centroids so that Γ_{0} is a function of ϕ_{0}. Remark 3.1. As noted by [8], the tangent space TC_{1} do not depend on η_{0}.

We introduce the cost functional j(η_{0}) = J(η_{0}; u) and in order to minimise this functional, we compute a sensitivity gradient of j, i.e. we differentiate the solution of (2.14) with respect to η_{0}. If this solution is smooth enough then :

\frac{dj(η_{0})}{dη_{0}}Δη_{0} = \frac{∂J}{∂η_{0}} Δη_{0} + \frac{∂J}{∂u} \frac{∂u}{∂η_{0}} Δη_{0}. (3.1)

To get rid of \frac{∂u}{∂η_{0}} we introduce an adjoint equation : find a ∈ T C_{1,2,3} such that :

(3.2) ∀v ∈ T C_{1,2,3}, \frac{∂G(η_{0}; a, v)}{∂u}dv = −\frac{∂J(η_{0}; u)}{∂u}dv

With the use of (3.2), the sensitivity gradient can be expressed as :

\frac{dj(η_{0})}{dη_{0}}Δη_{0} = \frac{∂J}{∂η_{0}} Δη_{0} + \frac{∂G(η_{0}; u, a)}{∂η_{0}}Δη_{0} (3.3)

where a is a solution of the adjoint equation.

A precise proof of (3.3) can be found in [9] in the case where (3.2) is invertible. This is the case if we are at regular path but other equations have to be derived for limit points and bifurcation points[20], [21].

3.2 Thickness sensitivity

In this part we consider the case where the design variable η_{0} = e is a parametrisation of the thickness of the beam. We suppose that the cross-section is rectangular with a thickness e : [0, L] → ℝ_{∗}^{+} and a constant width b. The formula of C_{m} and C_{r} can be written explicitly as [18] :

A = eb, A_{2} = be^{2}/4, A_{3} = eb^{2}/4, I_{2} = be^{3}/12 I_{3} = eb^{3}/12. (3.4)

The torsional stiffness has a more complex expression and is given in closed form[22] by :

J = be^{3}/3 − \frac{64e^{4}}{π^{5}} ∑_{k=0}^{+∞} \frac{tanh \frac{(2k+1)πb}{2e}}{(2k + 1)^{5}}

A differentiation with respect to e is easily performed. The sensitivity gradient can be written :

(3.5) ∀v ∈ T C_{1,2,3}, \frac{dG(e, u)}{de}v = ∫_{0}^{L}(γ(u)δv)^{T} C′(e)γ(u) dS

where C′(e) is the derivative of the elasticity tensor with respect to e.

3.3 Shape sensitivity

In this part we considere the case where η_{0} = (ϕ_{0}, θ⃗_{0}); we recall that the undeformed beam may have an initial strain so that so that we have to take into account the position of the cross section θ⃗_{0}. The solution of equilibrium equation (2.14) becomes : find (ϕ, θ⃗) ∈ C_{1,2,3} such that ∀(δϕ, δθ⃗) ∈ TC_{1,2,3}, G(ϕ_{0}, θ⃗_{0}; ϕ, θ⃗)[δϕ, δθ⃗] = 0, where G = G_{i} − G_{e}. The internal and external parts have the following form :

G_{i}(ϕ_{0}, θ⃗_{0}; ϕ, θ⃗)[δϕ, δθ⃗] = ∫_{0}^{L} δϵ(ϕ_{0}, θ⃗_{0}; ϕ, θ⃗)^{T} C_{M}(e, b)ϵ(ϕ_{0}, θ⃗_{0}; ϕ, θ⃗) dS
+ ∫_{0}^{L} δκ(ϕ_{0}, θ⃗_{0}; ϕ, θ⃗)^{T} C_{R}(e, b)κ(ϕ_{0}, θ⃗_{0}; ϕ, θ⃗) dS
G_{e}(ϕ_{0}, θ⃗_{0}; ϕ, θ⃗)[δϕ, δθ⃗] = ∫_{0}^{L} (N̄(ϕ, θ⃗)δϕ + M̄(ϕ, θ⃗)T_{m}(θ⃗)δθ⃗) dS.

(3.6)

We seek the equation satisfied by the derivative of the solution (ϕ, θ⃗) with respect to η_{0} = (ϕ_{0}, θ⃗_{0}). We note Δ(.) the derivative of (.) with respect to shape variables in contrast to δ(.) which denotes derivative with respect to state variables.

We should be cautious about functions like ϕ, θ⃗: they do not depend explicitely on ϕ_{0}, θ⃗_{0} but their derivative with respect to S, the arc lenth along the line of centoids of the unloaded beam does depend on this line; we get easily :

(3.7) Δ\frac{dΦ}{dS} = −\frac{dΦ}{dS}(\frac{dΦ_{0}}{dξ} . Δ\frac{dΦ_{0}}{dξ})\frac{1}{∥\frac{dΦ_{0}}{dξ} ∥^{2}}

and a similar equation for Δ\frac{dθ⃗}{dS}. By differentiation wrt arc length of the line of centroids of the unloaded beam, we get :

Δ dS = \frac{dΦ_{0}}{dξ} .Δ\frac{dΦ_{0}}{dξ} \frac{dξ}{∥\frac{dΦ_{0}}{dξ} ∥} . (3.8)

Thus the change of length is given by :

Δϵ = Γ^{T} Δϕ′ − ΔE⃗_{1} with ΔE_{1} = (Γ_{0}T_{m}(θ⃗_{0})Δθ⃗_{0}) × E⃗_{1}.

For the change of curvature κ⃗, we get :

Δκ⃗(ϕ, Γ) = T_{m}(θ⃗)Δθ⃗′ + ΔΓ(θ⃗_{0})T_{m}(θ⃗_{0})θ⃗′_{0} + Γ(θ⃗_{0})ΔT_{m}(θ⃗_{0})θ⃗′_{0} + Γ(θ⃗_{0})T_{m}(θ⃗_{0})Δθ⃗′_{0}.

We get similar expressions for Δδϵ and Δδκ. We note that C_{M} and C_{R} do not depend on ϕ_{0} and θ⃗_{0}. We are know ready to differentiate (3.6) with respect to η_{0} : forall Δη_{0} = (Δϕ_{0}, Δθ⃗_{0}), the shape gradient is given by :

ΔG_{i} = \frac{dG_{i}}{dη_{0}} (η_{0}; u, δu)Δη_{0}

= ∫_{0}^{L}(δϵ^{T} C_{M}ϵ + δκ^{T} C_{M}κ) ΔdS + ∫_{0}^{L}(Δδϵ^{T} C_{M}ϵ + Δδκ^{T} C_{M}κ) dS + ∫_{0}^{L}(δϵ^{T} C_{M}Δϵ + δκ^{T} C_{M}Δκ) dS.

The variation of the work of external basis should be considered in a case study.

4 Conclusion

We have recalled a large displacement beam model; for this model we have presented a complete design sensitivity. Although not presented here, a finite element aproximation has been completed; C or fortran code has been generated with some helps from a formal computation system. Numerical results will be presented in another paper.

References

[1] P. Aubert and B. Rousselet. ‘Sensitivity computations and shape optimization for a nonlinear arch model with simple instabilities.’ International Journal for numerical methods in engineering, 42:15–48 (1998).

[2] I. Fried. ‘Nonlinear finite element computation of the equilibrium, stability and motion of the extensional beam and ring.’ Computer Methods in applied mechanics and engineering, 38(1):29–44 (1983).

[3] E. Reissner. ‘On a one dimensional finite strain beam theory: the plane problem.’ J. Appl. Math. Phys. (1972).

[4] S. S. Antman and C. S. Kenney. ‘Large buckled states of nonlinearly elastic rods under torsion, thrust, and gravity.’ Archive for rational mechanics and analysis, 289–338 (1981).

[5] J. Simo and L. Vu-Quoc. ‘A three dimensional finite-strain rod model. part II computational aspects.’ Computer Methods in applied mechanics and engineering, 58:79–116 (1986).

[6] J. Bourgat et al. ‘Modélisation et calcul des grands déplacements de tuyaux élastiques en flexion-torsion.’ International Journal of Bifurcation and Chaos, 7(4):379–408 (1988).

[7] A. Ibrahimbegović et al. ‘Computational aspects of vector-like parametrisation of three-dimensional finite rotations.’ Tech. Rep. 6, École Polytechnique Fédéral de Lausanne, Laboratoire de mécanique des structures et milieux continus (1994).

[8] B. Rousselet. ‘A finite strain rod model and its design sensibility.’ Mech. Struct.& Mach., 20(4):415–432 (1992).

[9] B. Rousselet et al. ‘Design sensitivity for a hyperelastic rod in large displacement with respect to its midcurve shape.’ J. Optim. Theory Appl., 96(3):683–708 (1998).

[10] J. Simo. ‘A finite strain beam formulation. the three dimensional dynamic problem. part I.’ Computer Methods in applied mechanics and engineering, 49:55–70 (1984).

[11] A. Ibrahimbegović. ‘On the choice of finite rotation parameters.’ Computer Methods in applied mechanics and engineering, (149):49–71 (1997).

[12] J. C. Alexander and S. S. Antman. ‘Global and local behavior of bifurcating multidimensional continua of solutions for multi parameter nonlinear eigenvalue problems.’ Archive for rational mechanics and analysis, 339–354 (1989).

[13] A. Ibrahimbegović et al. ‘On the role of geometrically exact and second order theories in buckling and postbuckling analysis of three-dimensional beam.’ Tech. rep., Department of mechanical engineering, Université de technologie de Compiègne, Lab. Lg2mS, Div. MNM, BP-649 Compiègne, France (1996).

[14] A. Ibrahimbegović and M. A. Mikdad. ‘On dynamics of finite rotations of 3d beams.’ In Computational methods in applied sciences, 447–453. Eccomas, Wiley (september 1996).

[15] J. Stocker. Differential Geometry. Wiley (1989).

[16] M. D. Carmo. Differential geometry of curves and surfaces. Prentice Hall, Englewood Cliffs, New Jersey (1976).

[17] A. Ibrahimbegović. ‘Stress resultant geometrically nonlinear shell theory with drilling rotations, part I. a consistent formulation.’ Computer Methods in applied mechanics and engineering, 118:264–284 (1994).

[18] I. Sokolnikoff. Mathematical theory of elasticity. Mc Graw-Hill, 2nd edn. (1956).

[19] R. Adams. Sobolev spaces. Academic Press (1978).

[20] Z. Mróz and J. Piekarsky. ‘First and second order design sensitivity at a bifurcation point.’ In J. Herskovits, ed., Advances in Structural Optimization, vol. 25 of Solid Mechanics and its Applications, 349–380. Kluwer Academic Publishers, Dordrecht, The Netherlands (1995).

[21] P. Aubert. ‘Thickness optimization of a geometrically non linear arch at limit points.’ Applied Mathematics and Computer Science, 6(2):101–117 (1996).

[22] J. P. Smith. ‘Highly accurate beam torsion solutions.’ Tech. Rep. 3608, Nasa/Johnson space center (1996).