Mixed methods and flexibility approaches for nonlinear frame analysis

Journal of Constructional Steel Research 58 (2002) 967–993 www.elsevier.com/locate/jcsr

Mixed methods and flexibility approaches for nonlinear frame analysis K.D. Hjelmstad a,*, E. Taciroglu b a

b

Department of Civil and Environmental Engineering, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA Center for Simulation of Advanced Rockets, University of Illinois at Urbana-Champaign, Urbana, IL 61801, USA Received 1 June 2001; received in revised form 16 August 2001; accepted 29 November 2001

Abstract The principle of virtual forces has long been an attractive tool for linear matrix structural analysis because it provides a means to compute the exact flexibility matrix for a linear, elastic, non-prismatic Bernoulli–Euler beam element while the classical stiffness method, based on the principle of virtual displacements, does not. Nonlinear flexibility methods have recently appeared in the literature as a means of improving the accuracy of frame analysis when the curvature fields are more complex than the moment fields (e.g., when inelastic response is important). This paper surveys the basic methods of formulating nonlinear Bernoulli–Euler beam elements based on classical and mixed variational principles. We show that the variational structure provided by the Hellinger–Reissner and Hu–Washizu functionals gives a framework that reaps the benefits of the so-called nonlinear flexibility methods while allowing a standard treatment of non-linear problems by Newton’s method.  2002 Elsevier Science Ltd. All rights reserved. Keywords: Beam finite elements; Mixed-enhanced elements; Mixed finite element methods; Hu–Washizu; Hellinger–Reissner; Variational principles; Nonlinear flexibility methods

* Corresponding author. E-mail address: [email protected] (K.D. Hjelmstad). 0143-974X/02/$ - see front matter  2002 Elsevier Science Ltd. All rights reserved. PII: S 0 1 4 3 - 9 7 4 X ( 0 1 ) 0 0 1 0 0 - 6

968

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

1. Introduction The analysis of framed structures is almost exclusively done within the context of Bernoulli–Euler beam theory. This theory is easily extended to problems in which inelasticity manifests using the theory of plasticity in stress resultant space through postulation of axial force and flexural moment interaction relationships and associative flow rules. However, when the constitutive equations are nonlinear, many of the advantages of the classical displacement-based structural analysis procedures vanish. In particular, it is no longer possible to exactly solve the governing differential equations of motion. Methods based upon exact stiffness coefficients are only valid for ‘plastic hinge’ approaches where all of the inelasticity is considered to manifest in deformations at a certain point (i.e., the location of the ‘hinge’). Plastic hinge methods often do not represent frame response well, particularly if strain hardening is important. For these problems ‘plastic zone’ models that are capable of representing the distributed inelasticity are better. Displacement-based methods (usually called ‘stiffness’ methods), which are formulated using the principle of virtual work, provide an approximate approach to solving problems in which the distribution of curvature is different from the distribution of bending moment. These approximations do, of course, converge in the limit as more elements are taken to discretize the beam. However, mesh refinement in framed structures is not generally well-understood (or practiced) by those who use frame analysis. Furthermore, the amount of refinement needed to accurately capture regions of steep curvature often makes the problem difficult and tedious. Accurate estimates of response in these regions are essential to design. This paper examines alternative methods of estimating beam response. For nonlinear problems the difference between internal resistance and applied load is a nonlinear function of the deformed configuration, r(u). The solution of this nonlinear problem by Newton’s method involves iteratively solving the system of equations Ki⌬ui ⫽ ⫺r(ui), where Ki ⬅ ∂r / ∂u is the tangent stiffness evaluated at the configuration ui, and updating the configuration as ui ⫹ 1 ⫽ ui ⫹ ⌬ui until ||r(ui)|| ⬍ tol. The role of the stiffness matrix in nonlinear analysis is clearly different from linear analysis. If the tangent stiffness matrix is consistently derived from the residual then Newton’s method converges quadratically. The classical methods of structural analysis benefit from the variational structure of the principles of virtual work. Displacement-based methods have been more popular for structural analysis because the continuity of the displacement field can be trivially satisfied by selecting shape functions whose coefficients are the values of the displacements at the ends of the member and by simply assigning a common displacement variable to those nodes that connect the members. The familiar notion of assembly of the stiffness matrix and force vector from element contributions takes care of establishing the proper continuity. Because the classical flexibility method is based upon force interpolation, the interelement continuity of the displacements is difficult to enforce. On the other hand, the force-based interpolation has the advantage that, at least for linear problems, it can generate the exact flexibility matrix for the element, even for nonprismatic members.

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

969

Recently, Ciampi and Carlesimo [1], Spacone et al. [2,3], Neuenhofer and Filippou [4], and others have proposed algorithms that attempt to bring the benefits of the flexibility method (i.e., the ability to compute an exact flexibility matrix) to problems with nonlinear material response. The motivation for these efforts is very clear. For a typical structural element in a frame the maximum forces tend to occur at the ends. Hence, the ends of the elements are the zones where the inelastic deformations accumulate. The bending strain (curvature) field along the length of an element can be highly nonlinear while the moment field (in absence of applied loads) remains linear. These methods, nominally based on the Hellinger–Reissner variational principle, are generally referred to as ‘nonlinear flexibility methods.’ The difficulty with these methods is that they lack a complete variational structure. Furthermore, there seems to be no natural way to extend them to account for geometric nonlinearities. Petrangeli and Ciampi [5] discuss mixed variational principles as the more natural setting for understanding the nonlinear flexibility methods and make a crucial connection with the work of Simo and Rifai [6] who discuss ‘assumed strain’ methods in the context of two- and three-dimensional solids. The purpose of the present paper is to summarize, clarify, and extend what is known about the application of mixed variational principles to Bernoulli–Euler beam theory, particularly as it applies to the analysis of framed structures. By providing a survey that starts with the classical principles of virtual work we hope to clarify and extend the ideas presented by Petrangeli and Ciampi [5] and others. One of our main goals is to make these ideas accessible to anyone with a classical education in the methods of structural analysis. As such, this paper is partly a state-of-the-art survey, partly a transfer of technology from two- and three-dimensional finite elements back to nonlinear analysis or reticulated structures, and partly an exposition of new ideas. To keep the discussion and notation simple, we will focus on the bending problem and ignore axial effects. This limitation will make the theoretical developments much easier to follow but will technically restrict us to continuous beam structures, not arbitrary plane frames. We shall amend these limitations in a forthcoming paper. In this paper we remind the reader of the basic features of the two dual principles of virtual work and their manifestations when implemented for the analysis of reticulated structures. We demonstrate that mixed finite element methods, based upon the Hellinger–Reissner and Hu–Washizu variational principles, have the advantages of both the flexibility approach (exact stiffness matrices for linear problems) and the stiffness approach (natural enforcement of displacement continuity). Further, we suggest that the so-called nonlinear flexibility methods do not have a consistent variational structure.

2. Bernoulli–Euler beam theory and the principles of virtual work Let us consider the geometrically linear, planar Bernoulli–Euler beam theory (i.e., shearing and axial deformations are constrained to be zero) as a model problem. The classical equilibrium and strain–displacement equations of this theory are

970

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

m⬙ ⫽ q, ␬ ⫽ w⬙

(1)

where m is the bending moment field, q is the applied transverse load per unit length, ␬ is the curvature field, and w is the transverse displacement field. All of these fields are functions of the axial coordinate x which is measured from the left end of the beam, as shown in Fig. 1 which also indicates the convention for positive values of the variables. A prime designates differentiation with respect to the axial coordinate. Throughout this paper a lowercase italic character will represent a field variable that depends upon the axial coordinate whereas a bold character will represent a vector of parameters that do not depend upon the spatial coordinate x. In Bernoulli–Euler beam theory the bending moment m plays the role of the internal stress resultant, while the curvature ␬ is the corresponding strain resultant. These two variables are related through a constitutive hypothesis, which we write alternately as m⫽m ˆ (␬), ␬ ⫽ ␬ˆ (m)

(2)

ˆ (␬) indicates that the moment is a function of the curvature. SimiThe notation m larly, ␬ˆ (m) indicates that the curvature is a function of the moment. The functions ˆ and ␬ˆ are duals of each other, in the sense that they are meant to represent exactly m the same constitutive model. These two representations can have significantly different ramifications in the numerical implementation of the model. For linear, elastic ˆ (␬) ⫽ D␬, where D is the bending beams the constitutive function takes the form m stiffness of the cross section (usually taken as Young’s modulus times the second moment of the cross sectional area about the centroid). For the linear model ␬ˆ (m) ⫽ Cm, where C ⫽ 1 / D is the compliance. For non-prismatic elastic members the moduli C and D vary with x but do not evolve with the deformation. In nonlinear elasticity and inelasticity we generally think of the moduli as the current tangent to the moment–curvature response function and so the idea of modulus can only be understood in an incremental sense. Remark 1. For nonlinear problems, we will think of the relationship between m and ␬ as a prescribed nonlinear function. For inelastic structures it is common to express the constitutive equations in rate form rather than as a direct functional dependence. Once the inelastic rate equations are integrated numerically they can be expressed as m ⫽ m ˆ (␬,h) where h is a vector of internal variables that may arise

Fig. 1. Planar Bernoulli–Euler beam (a) definition of position and loading, (b) definition of positive transverse displacement and bending moment.

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

971

in the description of the constitutive model (e.g., the plastic strain in an elasto⫺plastic problem). Throughout these developments we will assume that the internal variables are treated in a classical sense. Therefore, it will suffice for our purposes here to think of the moment as being a function of the curvature alone (and vice versa).䊐 2.1. The principles of virtual work The principles of virtual work are based on the fundamental theorem of the calculus of variations in conjunction with the notion of the work done by forces moving through displacements (see, for example, [7]). In the principle of virtual displacements, the forces are taken to be the real forces acting on the system. The virtual work is the work done by those forces when subjected to imagined (virtual) displacements. In the principle of virtual forces, the displacements are taken to be the real displacements of the system. The virtual work is the work done when imagined (virtual) forces act through those displacements. The principle of virtual displacements leads to an alternative statement of equilibrium while the principle of virtual forces leads to an alternative statement of compatibility. The principle of virtual displacements uses the real loads q and the real moments m and posits a virtual displacement field w¯ . The difference between the internal and external virtual work gives the functional S(m,q,w¯ ) ⬅

冕

ᐉ

(mw¯ ⬙⫺qw¯ )dx

(3)

0

The principle of virtual displacements is as follows: If S(m,q,w¯ ) ⫽ 0 for all w¯ that satisfy the essential boundary conditions (i.e., prescribed displacements and rotations and the ends of the beam), then the system is in equilibrium (i.e., m⬙ ⫽ q). One can prove this assertion by integrating the first term by parts twice and invoking the fundamental theorem of the calculus of variations [8]. The principle of virtual forces uses the real displacements w and real curvatures ␬ and posits a virtual force field m ¯ that is in equilibrium with an applied virtual load q¯ (i.e., the virtual moment field must satisfy m ¯ ⬙ ⫽ q¯ ). The difference between internal and external complementary virtual work gives the functional F(␬,w,m ¯) ⬅

冕

ᐉ

(␬ m ¯ ⫺wm ¯ ⬙)dx

(4)

0

The principle of virtual forces is as follows: If F(␬,w,m ¯ ) ⫽ 0 for all m ¯ , then the curvatures are compatible with the displacements, (i.e., ␬ ⫽ w⬙). The character of the second term in the integral in Eq. (4) as external work can be clearly seen through its equivalent expression wq¯ . One of the best known applications of the principle of virtual forces is the dummy unit load method for computing the deflections in structures, for which the strain field is known, by specializing to a loading q¯ ⫽ d(x0), where d(x0) is a delta distribution centered at x ⫽ x0. It is trivial to prove the assertion of the principle of virtual forces in the same way as the principle of virtual displacements above.

972

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

3. The classical stiffness and flexibility methods The classical implementation of the stiffness and flexibility methods have formed the backbone of linear structural analysis for a long time (see, for example, [9]). The two classical methods, based on the principle of virtual displacements and force, respectively, really need no introduction, but are summarized here to illuminate the features that naturally appear in each so that we can set up the notation and provide the formulations that will connect with the nonlinear flexibility methods and mixed methods later in the paper. 3.1. The classical stiffness method In the displacement-based formulation the moment m that appears in Eq. (3) is viewed as being a (constitutive) function of the curvature ␬, which, according to the classical equations, is the second derivative of the displacement field ␬ ⫽ w⬙. Let us enforce both the constitutive equation and the strain–displacement equation in a classical sense. (For the linear theory, for example, such an implementation is trivˆ (␬(w)) ⫽ Dw⬙). With these constraints, the virtual work functional ially realized as m involves only the displacement field w (and its virtual counterpart). Hence, we write the functional as

冕

S(w,w¯ ) ⬅

ᐉ

ˆ (␬)w¯ ⬙⫺qw¯ )dx (m

(5)

0

As a notational convenience we designate the argument of the constitutive function as the curvature with the understanding that the curvature is a function of the disˆ (␬(w)). The placement field. Hence, we understand the functional dependence as m abbreviated notation will be used throughout, with an indication in each case of the exact functional dependence. The dependence is important for linearization of the functional. Let us assume that we will solve the nonlinear problem incrementally by Newton’s method. Thus, we must linearize the functional S with respect to the displacement w. Let wi represent a configuration of the beam (not necessarily an equilibrium configuration). We can linearize the constitutive equation about this configuration as ˆ (␬i) ⫹ Di⌬w⬙ ˆ ( ␬) ⫽ m m

(6)

where m ˆ (␬i) is the moment field associated with the curvature ␬i ⫽ wi⬙. The tangent ˆ / d␬ is evaluated at the state wi and ⌬w ⬅ w⫺wi is the incremental modulus Di ⬅ dm displacement field. The linear part of virtual work functional is ᏸ[S] ⬅

冕

ᐉ 0

ˆ (␬i)w¯ ⬙⫺qw¯ )dx ⫹ (m

冕

ᐉ

w¯ ⬙Di⌬w⬙dx

(7)

0

In Newton’s method we set ᏸ[S] ⫽ 0 for all w¯ to provide a means of estimating an incremental state ⌬w that, when added to the estimate wi, will come closer to

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

973

satisfying the nonlinear principle of virtual work S(w,w¯ ) ⫽ 0. The process can be repeated to convergence. In the classical stiffness method, we interpolate the displacement field, the incremental displacement field, and the virtual displacement field as a linear combination of known base functions hT ⬅ {h1(x),…,hN(x)}, where N is the number of terms in the interpolation. The interpolation of the three fields can be written as w ⫽ hTw, w¯ ⫽ hTw ¯,

⌬w ⫽ hT⌬w

(8)

where w, w ¯ , and ⌬w are the discrete displacement, virtual displacement, and incremental displacement parameters, respectively. Substituting these interpolations into the linearized virtual work functional and insisting that ᏸ[S] ⫽ 0 for all w ¯ gives the discrete stiffness equations ˆ h(wi) Ki⌬wi ⫽ q⫺m

(9)

ˆ h(w) where the tangent stiffness Ki, the equivalent load q, and the internal force m are given from the general expressions Ki ⬅

冕

ᐉ

h⬙Dih⬙Tdx, q ⬅

0

冕

ᐉ

ˆ h(w) ⬅ qhdx, m

0

冕

ᐉ

ˆ (␬)h⬙dx m

(10)

0

Here we introduce a notational convention that will be used throughout the paper ˆ h(w) in different forms. We indicate the field w as the argument of the function m ˆ h(w) to remind us that w is the unknown field sought in the problem and that m involves an integral of that field. We recognize, however, that the field w is interpolated with base functions. Thus, the actual unknowns are the generalized displaceˆ h(w) is selected to remind us which base ment parameters w. The subscript on m functions are involved in the integral. Equilibrium is established by iteratively solving Eq. (9) and updating the estimate of the displacements with the equation wi ⫹ 1 ⫽ wi ⫹ ⌬wi, until the force residual ˆ h(wi)|| ⬍ tol. The iterative computation must be started is small enough, i.e., ||q⫺m with an estimate w0 and the tolerance tol on the satisfaction of equilibrium must be specified a priori. Remark 2. In the linear case, if the functions h are capable of representing the exact solution to the classical equations, then K is exact. For example, the cubic Hermitian polynomials h1(x) ⫽ 1⫺3x2 ⫹ 2x3 h2(x) ⫽ 3x2⫺2x3 h3(x) ⫽ (x⫺2x2 ⫹ x3)ᐉ

(11)

h4(x) ⫽ (⫺x2 ⫹ x3)ᐉ

where x ⬅ x / ᐉ, are exact for constant modulus. The nonlinear case is more like the non⫺prismatic case because the tangent modulus varies along the length of the beam. Hence, K is only an approximation of the true tangent stiffness of the beam for nonlinear material response.䊐

974

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

3.2. The classical flexibility method The principle of virtual forces forms the variational basis of the classical flexibility method. As a practical simplification, the virtual load can be taken as q¯ ⫽ m ¯⬙ ⫽ 0 to eliminate the term involving w from the functional in Eq. (4). If we view the curvature as being a constitutive function of the moment field, i.e., ␬ ⫽ ␬ˆ (m), then the functional involves only the moment field m (and, of course, its virtual counterpart). We write the functional as

冕

ᐉ

F(m,m ¯) ⬅

␬ˆ (m)m ¯ dx

(12)

0

Let us again assume that we will solve the nonlinear problem incrementally by Newton’s method. Hence, we must linearize the functional F with respect to the moment m. Let mi represent a moment field (not necessarily in equilibrium with the applied loads). We can linearize the constitutive equation about this configuration as

␬ˆ (m) ⫽ ␬ˆ (mi) ⫹ Ci⌬m

(13)

where ␬ˆ (mi) is the curvature field associated with the moment mi, Ci ⬅ d␬ˆ / dm is the tangent compliance evaluated at the state mi, and ⌬m ⬅ m⫺mi is the incremental moment field. The linearization of the virtual work functional about the state mi is ᏸ[F] ⬅

冕

ᐉ

0

␬ˆ (mi)m ¯ dx ⫹

冕

ᐉ

m ¯ Ci⌬mdx

(14)

0

Again, we try to find the incremental moment field that satisfies ᏸ[F] ⫽ 0 for all m ¯ to provide an improved estimate of the moment field that satisfies F(m,m ¯ ) ⫽ 0. The process is repeated until the solution converges. In the classical flexibility method we interpolate the moment, virtual moment, and incremental moment fields as a linear combination of known base functions bT ⬅ {b1(x),…,bM(x)}, where M is the number of terms in the interpolation. The interpolations can be written as m ⫽ bTm ⫹ mp, m ¯ , ⌬m ⫽ bT⌬m ¯ ⫽ bTm

(15)

where m, m ¯ and ⌬m are the (constant) coefficients of the real, virtual, and incremental bending moment fields, respectively, and mp(x) is the particular moment field that equilibrates the applied loads q when m ⫽ 0. If m represents the two end moments, for example, then mp is simply the bending moment field of a (statically determinate) beam with pinned ends. The freedom in selecting base functions is limited by the fact that the virtual bending moment field m ¯ must equilibrate the applied virtual loads. Since q¯ ⫽ 0, the virtual bending moment field must be linear, and the base functions must be some linear combination of the functions {1,x}. This choice establishes the so-called ‘essential’ or statical degrees of freedom of the element. Substituting these interpolations into the linearized virtual work functional and insisting that ᏸ[F] ⫽ 0 for all m ¯ gives the discrete equations

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

Fi⌬mi ⫹ ␬ˆ b(mi) ⫽ 0

975

(16)

where the tangent flexibility Fi and the residual ␬ˆ b(m) can be computed as Fi ⬅

冕

ᐉ

bCibTdx, ␬ˆ b(m) ⬅

0

冕

ᐉ

␬ˆ (m)bdx

(17)

0

where Ci is the tangent compliance d␬ˆ / dm, evaluated at the state mi. Again, note that we use the notation ␬ˆ b(m) to remind us that this quantity depends upon the moment field. The moment field is, of course, computed through the interpolation given in Eq. (15). For a single beam element, the problem can be solved in complete analogy with the stiffness approach by iteratively solving Fi⌬mi ⫹ ␬ˆ b(mi) ⫽ 0 and updating the estimate of the moments with the equation mi ⫹ 1 ⫽ mi ⫹ ⌬mi until the residual is small enough, i.e., ||␬ˆ b(mi)|| ⬍ tol. The iterative computation must be started with an estimate m0 and the tolerance tol must be specified a priori. Remark 3. If the functions b are capable of representing the exact solution to the classical homogenous equations, then F is exact. For linear problems one can use the inverse of this flexibility matrix to compute an exact stiffness matrix for use in a global displacement based solution procedure. The stiffness matrix can be computed from the flexibility as K ⬅ LF⫺1LT

(18)

where L is a statical transformation matrix relating the element forces m to the complete set of element end forces associated with the kinematic degrees of freedom. This stiffness matrix can be used to carry out matrix structural analysis in exactly the same way as the stiffness matrix that arises from the principle of virtual displacements or the one obtained by using the classical solution to the differential equation. The merit of this stiffness matrix is that it yields exact nodal displacements for any variation of cross sectional modulus because the (linear) force interpolations are exact.䊐 3.3. Global structural analysis by the displacement method The preceding developments were for a single element having any of a variety of boundary conditions at the two ends. In the stiffness approach the boundary conditions are imposed by having the interpolation functions h satisfy the displacement boundary conditions. In the flexibility approach the boundary conditions are implicit in the particular moment field mp and the selection of interpolation functions b. Structural analysis is about finding the response of structures composed of many elements. Hence, we must extend the formulation to account for the many elements and the communication of information between them. The global virtual work functional for the principle of virtual displacements can be constructed by summing element virtual work functionals and subtracting the work of the applied nodal loads. Hence, let us define the global functional

976

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

冘冕 E

SG ⬅

e⫽1

ᐉe

冘 N

(mew¯ ⬙e⫺qew¯ e)dxe⫺

0

¯ n ⫹ QnW ¯ n) (Mn⍜

(19)

n⫽1

¯ n} are the ¯ n,⍜ where (Qn, Mn) are the nodal applied force and moment at node n, {W nodal virtual displacement and rotation at node n, E is the number of elements, and N is the number of nodes in the structure. Equilibrium of the structure is guaranteed if SG ⫽ 0 for all virtual displacement fields w¯ (a function defined on the entire struc¯ n and first derivature that takes on values in each element of w¯ e(xe) and has value W ¯ n at node n). tive ⍜ To see why the above assertion of global equilibrium is true, one can integrate the first term in Eq. (19) by parts twice to get the equivalent expression

冘冕 E

SG ⫽

e⫽1

ᐉe

0

冘 E

(m⬙e⫺qe)w¯ edxe ⫹

e⫽1

冘 N

sTew ¯ e⫺

¯n RTnW

(20)

n⫽1

where se ⬅ {⫺m(0), m⬘(0), m(ᐉ), ⫺m⬘(ᐉ)} and w ¯ e ⬅ {w¯ (0), w¯ ⬘(0), w¯ (ᐉ), w¯ ⬘(ᐉ)} are the element end forces and virtual displacements that appear as a result of the inte¯ n} are the nodal forces and ¯ n, ⍜ ¯ n ⬅ {W gration by parts, and Rn ⬅ {Qn, Mn} and W virtual nodal displacements, as shown in Fig. 2. Note that ⫺m⬘ is the shear force, in accord with the classical equations of equilibrium. It is not difficult to show that, for any structure, the sum over all E elements of the virtual work done by the element end forces can be reorganized as a sum over nodes of a sum of all elements framing into each node. Let us write se ⬅ {sie, sje}, where sie is the force and moment acting on the end of the element attached to global node i, and sje, is the force acting on end of the element attached to global node j. Similarly, w ¯ ie, w ¯ je}, are the element end displacements and rotations at global ¯ e ⬅ {w nodes i and j. We can recast Eq. (20) as

Fig. 2.

Definition of nodal displacements and element end forces and end displacements.

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

冘冕

ᐉe

E

SG ⫽

e⫽1

冘 N

(m⬙e⫺qe)w¯ edxe⫺

0

¯n⫹ [RTnW

n⫽1

冘

snT ¯ ne] e w

977

(21)

e苸Kn

where Kn is the set of element numbers that frame into global node n. If the elements share global displacement degrees of freedom at the nodes, then, by construction, ¯ n, and the the virtual displacement field is continuous at the nodes. Hence, w ¯ ne ⫽ W functional can be written as

冘冕

ᐉe

E

SG ⬅

e⫽1

冘 N

(me⬙⫺qe)w¯ edxe⫺

0

¯ Tn[Rn ⫹ W

n⫽1

冘

sne]

(22)

e苸Kn

The fundamental theorem of the calculus of variations thus guarantees that the moment field in each element equilibrates the applied loads in that element and that the element end forces are in equilibrium with the applied nodal forces at each node. The power of the stiffness approach derives from the simplicity of enforcing the continuity of the displacement field at the nodes. It is easy to see that any functional with the term mw¯ ⬙ will have this feature. In particular, the mixed methods presented in the sequel share this feature with the stiffness approach, and hence these methods are ideally suited to matrix structural analysis. 3.4. Global structural analysis by the force method The principle of virtual forces can be treated in a manner similar to the principle of virtual displacements, with some very different consequences. The global virtual work functional for the principle of virtual forces can be stated as

冘冕

ᐉe

E

FG ⬅

e⫽1

冘 N

(␬em ¯ e⫺wem ¯ e⬙)dxe⫺

0

¯ TnWn R

(23)

n⫽1

¯ n is the vector of virtual nodal forces and Wn, is the vector of nodal displacewhere R ments at node n. The principle of virtual forces for the structure as a whole then emanates from the condition that FG ⫽ 0 for all virtual moment fields m ¯ that are in equilibrium with the applied virtual transverse load q¯ ⫽ m ¯ ⬙ in each element and the ¯ n. virtual nodal loads R Following an argument similar to the principle of virtual displacements leads to the equivalent form of the global functional FG ⫽

冘冕 E

ᐉe

e⫽1

0

冘 N

(␬e⫺w⬙e)m ¯ edxe⫺

n⫽1

¯ TnWn ⫹ [R

冘

n s¯nT e we ]

(24)

e苸Kn

where Kn is the set of element numbers that frame into global node n. In order to invoke the fundamental theorem of the calculus of variations we must have a relationship between the virtual end forces and the virtual nodal forces. To wit, if ¯n ⫽ ⫺ R

冘

e苸Kn

s¯ne

(25)

978

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

then we can write Eq. (24) as FG ⫽

冘冕 E

ᐉe

e⫽1

0

冘冘 N

(␬e⫺w⬙e)m ¯ edxe ⫹

n s¯nT e [Wn⫺we ]

(26)

n ⫽ 1e苸Kn

The fundamental theorem of the calculus of variations then guarantees that the curvatures are equal to the second derivative of displacement in each element and the displacement at the end of each element equals the displacement of the node. Unfortunately, it is very difficult to satisfy Eq. (25) for statically indeterminate structures. The classical way around this dilemma is to introduce a set of redundant ¯ and their associated conjugate forces X and corresponding redundant virtual forces X displacements ⌬ (which, of course, are zero), and augment the virtual work functional ¯ . If X ⫽ 0, the resulting structure is statically determinate and to be F∗G ⬅ FG ⌬TX it is possible to write explicit relationships between the member end forces and nodal loads. In fact, it is possible to express the generalized moments in member e as a linear combination of loads and redundant forces as me ⫽ Be0R ⫹ BexX, where B0 and Bx, are the (constant) equilibrium transformation matrices of the structure. An identical relationship holds for the virtual generalized moments. Substituting these relationships, along with the interpolation m ⫽ bTm ⫹ mp in each element leads to the classical flexibility method. Clearly, with this extension, the flexibility method can be extended to nonlinear problems, but with the same disadvantage that it has for linear problems—it is difficult to automate the selection of redundant forces and thereby establish the matrices Be0 and Bex. An automatic approach to the flexibility method that resolves these difficulties has been described by Heath et al. [10].

4. Nonlinear flexibility methods Motivated by the observation that the curvature distribution in an inelastic beam can have steep gradients (for example, in the plastic hinge region) but that the distribution of moments is always linear for q ⫽ 0, in accord with the equations of equilibrium, researchers have proposed that the principle of virtual forces, or the flexibility method, will give more accurate results than a standard displacement-based finite element method [1–4]. These approaches are nominally based upon the Hellinger– Reissner variational principle (see, for example, [2]), but appear not to be completely variationally consistent in the sense described by Simo and Hughes [11]. The algorithm of Neuenhofer and Filippou [4] is exemplary of the general class of nonlinear flexibility methods. When recast in the present notation, the computation goes as follows. Assume that the state at step i⫺1 is known and the displacement increment ⌬wi has been computed from the global stiffness equations. The estimate of the moment can be computed alternatively from the curvature, which can be computed from the displacement field, through the constitutive relationship m⫽m ˆ (␬). The moment field should also satisfy the interpolation m ⫽ bTm with generalized moments m. The following algorithm gives a means of computing the generalized moments m and the curvature field ␬, and does so by introducing an

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

979

auxiliary field m ˜ (x) defined as the difference between the moment field determined from the interpolation and the moment field determined from the constitutive function. Let ⌬mi ⬅ F⫺1 i⫺1L⌬wi⌬mi, where L is the statical transformation described earlier and F is the flexibility matrix. The curvature field at the new state i is estimated from the expression

␬i ⫽ ␬i⫺1 ⫹ Ci⫺1[bT⌬mi ⫹ m ˜ i⫺1]

(27)

where m ˜ (x) is called the residual section moment. The section moment is interpolated from the updated generalized moment mi as mi ⫽ bTmi, where

冕

ᐉ

ˆ (␬i)⫺m ˆ (␬i⫺1)⫺m Ci[m mi ⫽ mi⫺1 ⫹ F⫺1 ˜ i⫺1]bdx i

(28)

0

where Fi is computed from Eq. (17). Finally, the residual section moment is updated as m ˆ (␬i) ˜ i ⫽ bTmi⫺m

(29)

The iteration can start with m ˜ 0 ⫽ 0. The integral in Eq. (28) can be carried out numerically. Hence, one must store values of (␬i, m ˜ i) at each integration point. In ) and m must be stored from the previous iteration. At each addition, Fi (and/or F⫺1 i i iteration the global stiffness matrix is assembled from element stiffness matrices that are computed from element flexibility matrices computed with Eq. (17) using Ci ˆ (␬i). Neuenhofer and Filippou [4] report that this formulation gives with mi ⫽ m excellent accuracy both in displacement and in strain compared to a displacementbased finite element method with a similar number of elements. The basis of the global analysis aspects of the various nonlinear flexibility methods reported in the literature are not very clear. Most authors set up the global equations KT⌬W ⫽ R that yields incremental nodal displacements ⌬W from the residual force R and the ‘tangent stiffness’ KT. In these approaches the element stiffness is generally computed as the inverse of the element flexibility matrix mentioned in Remark 3. Since the force interpolation and the moment field given by the constitutive function are identical, the global residual can be computed simply as the nodal force residual. As we shall show in Section 5, the two-field variational formulation is central to the connection between the element-level flexibility computations and the displacement-based global equations. However, the nonlinear flexibility methods employ a state determination algorithm that is not variationally consistent with the linearization of a two-field variational statement based on the original nonlinear equations. Hence, it is unlikely that the global stiffness matrix computed with these algorithms is the consistent tangent stiffness associated with the nodal force residual. In fact, the basic merit of the flexibility-based methods (i.e., the computation of an ‘exact’ stiffness matrix) may be irrelevant in an iterative computation where obtaining a quadratic (asymptotic) convergence rate depends upon having a tangent stiffness matrix that is consistently derived from the nonlinear residual function. None of the papers on

980

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

the nonlinear flexibility method have reported on the asymptotic rate of convergence of the global iterations.

5. Mixed variational principles and mixed finite element methods The Hellinger–Reissner and Hu–Washizu energy principles represent multi-field variational principles and provide a natural setting for the formulation of mixed finite element methods. In a mixed method, the opportunity to independently approximate all of the fields that appear in the functional exists and gives rise to interesting methods of analysis. In this section we show, in the context of the model problem, how these methods arise from variational considerations, how computations can be organized for nonlinear problems, and how these methods relate to the nonlinear flexibility methods. 5.1. Method based on the Hellinger–Reissner type functional Consider the following two-field (Hellinger–Reissner type) functional for the Bernoulli–Euler beam H(w,m,w¯ ,m ¯) ⬅

冕

ᐉ

[mw¯ ⬙⫺qw¯ ⫺(␬ˆ (m)⫺w⬙)m ¯ ]dx

(30)

0

This functional depends upon the two fields m and w (and their virtual counterparts). Note that the constitutive equation ␬ ⫽ ␬ˆ (m) is implemented in a classical sense. Integration of the first term by parts twice reveals that H ⫽ 0 for all m ¯ and w¯ (that satisfy the essential boundary conditions) implies that m⬙ ⫽ q and w⬙ ⫽ ␬ˆ (m). Remark 4. We have referred to the functional H as a Hellinger–Reissner type functional. In fact, the Hellinger–Reissner functional is an energy functional, not a virtual work functional, and it exists only under certain conditions. If the constitutive ˆ (m) such equation is hyperelastic, that is if there exists a stress energy function U ˆ / dm, then Vainberg’s theorem guarantees the existence of the that ␬ˆ (m) ⫽ dU energy functional H⬅

冕

ᐉ

ˆ (m)]dx [mw⬙⫺qw⫺U

(31)

0

that has the property that DH ·(w¯ ,m ¯ ) ⫽ H, i.e., the directional derivative of the energy H in the direction of the virtual fields (w¯ ,m ¯ ) is H. Thus, the virtual work functional H is simply a two⫺field variational principle that has a close variational relationship with the Hellinger–Reissner functional H.䊐 If we substitute our interpolations of displacement and moment, i.e., w ⫽ hTw and m ⫽ bTm and their virtual counterparts, then the discrete form of the functional H can be written as

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

H⫽w ¯ T[BTw⫺␬ˆ b(m)] ¯ T[Bmⴚq]⫺m

981

(32)

where the equivalent load q, defined in Eq. (10), is the same one that appears in the stiffness approach. The matrix B and vector ␬ˆ b(m) are defined as B⬅

冕

ᐉ

h⬙bTdx, ␬ˆ b(m) ⬅

0

冕

ᐉ

␬ˆ (m)bdx

(33)

0

Setting H ⫽ 0 for all w¯ and m ¯ implies that the displacement w and moment m are determined by the following set of equations Bm ⫽ q, BTw ⫽ ␬ˆ b(m)

(34)

where we have used the notation ␬ˆ b(m) to remind us that the residual depends upon the moment field which, of course, is determined from the moment parameters m through the interpolation. Observe that the first of these equations is linear in m. The matrix B that appears in both equations depends only on the interpolation functions and not on the state of deformation. Remark 5. It is interesting to note that Eq. (34)2 represents the classical moment– area equations when the moment interpolation b is linear. Recall that, with the linear interpolation the matrix B is identical to the statical transformation L, the transpose of which converts end displacements to relative end deformations as LTw.䊐 The nonlinear function ␬ˆ b(m) linearizes to ␬ˆ b(m)⬇␬ˆ b(mi) ⫹ Fi⌬mi about the state with moment field mi, where Fi is the flexibility matrix given by Eq. (17). The linearized system of equations that must be solved is B⌬mi ⫽ q⫺Bmi, BT⌬wi⫺Fi⌬mi ⫽ ␬ˆ b(mi)⫺BTwi

(35)

Provided that Fi is invertible, the incremental generalized moment can be eliminated from these equations to produce and incremental equations for the displacement alone as ¯ iwi ¯ i⌬wi ⫽ q⫺B[mi ⫹ F⫺1 K ˆ b(mi)]⫺K i ␬

(36)

¯ i ⬅ BF B is the tangent stiffness matrix. Note that, due to the linearity of Eq. K (34)2 with respect to the nodal displacements w, Eq. (36) can also be expressed as ⫺1 i

T

¯ iwi+1 ⫽ q⫺B[mi ⫹ F⫺1 K ˆ b(mi)] i ␬

(37)

where, as usual, wi ⫹ 1 ⫽ wi ⫹ ⌬wi is the updated displacement field. The increment to the moment field can then be computed by solving Fi⌬mi ⫽ BTwi ⫹ 1⫺␬ˆ b(mi) and updating the moment parameters as mi ⫹ 1 ⫽ mi ⫹ ⌬mi. The formulation based on the Hellinger–Reissner type functional appears to have the variational structure that is missing from the nonlinear flexibility methods. Because the Hellinger–Reissner type functional, like the standard displacement-based approach, has the feature that the equilibrium equations are satisfied in weak form, the global structural analysis involves only the standard continuity of the displacement field. Thus, the formulation can be carried out in the context of a standard (displacement-based) finite element program.

982

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

Unlike the displacement-based formulation, the constitutive equations are expressed in complementary form (strains are a function of stresses). This format is not the familiar ‘strain-driven’ format wherein the constitutive equations are satisfied locally assuming that the strains within each iteration are fixed (from the global displacement state) and stresses are found that satisfy the constitutive equations for those given strains. The strain-driven constitutive equations fit the framework of inelasticity computations much better than the stress-driven format of the Hellinger– Reissner type functional for multiaxial (interacting) states of stress and strain. However, for uniaxial states of stress and strain there is no particular disadvantage to the stress-driven format. Simo et al. [12] describe a variational approach to the problem of stress-driven inelasticity. On the other hand, the behavior of Newton’s method can be affected by the shape of the response function. It is well-known, for example that Newton’s method can experience problems for displacement-based formulations for problems in which the material stiffens with increasing deformation. One would naturally expect this same phenomenon to manifest for the Hellinger–Reissner formulation for problems in which the material softens with increasing deformation. Material softening is much more common in engineering materials than is material stiffening. Remark 6. Note that if the response is linear then ␬ˆ (m) ⫽ Cm and the function ␬ˆ b(m) reduces to the linear expression ␬ˆ b(m) ⫽ Fm. The necessary condition H ⫽ 0 for all m ¯ and w¯ implies the discrete equations Bm ⫽ q, BTw ⫽ Fm

(38)

The flexibility matrix F is guaranteed to be positive definite so we can solve the second of these equations for the moment as m ⫽ F⫺1BTw. We can substitute this result into the first equation to get [BF⫺1BT]w ⫽ q

(39)

It is straightforward to show that if b represents a linear interpolation, which is consistent with m⬙ ⫽ 0 for the homogenous linear problem, and h are the standard cubic Hermitian beam shape functions, then B ⫽ L, the statical transformation of Eq. (18). Therefore, the stiffness matrix that comes from the Hellinger–Reissner variational functional is exactly the element stiffness matrix computed by the flexibility method and enjoys the feature of being exact for nonprismatic members. 䊐 5.2. Methods based on the Hu–Washizu type functional Consider the following three-field (Hu–Washizu type) functional for the Bernoulli–Euler beam J(w,m,␬,w¯ ,m ¯ ,␬¯ ) ⬅

冕

ᐉ

ˆ (␬))␬¯ ]dx [mw¯ ⬙⫺qw¯ ⫺(␬⫺w⬙)m ¯ ⫺(m⫺m

(40)

0

This functional depends upon the three fields m, ␬, and w (and their virtual

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

983

counterparts). Note that none of the governing equations is implemented in a classical sense. Integration of the first term by parts twice reveals that J ⫽ 0 for all m ¯ , ␬¯ and w¯ ˆ (␬). implies that m⬙ ⫽ q, w⬙ ⫽ ␬, and m ⫽ m Remark 7. Like the Hellinger–Reissner functional, the Hu–Washizu functional is an energy functional, not a virtual work functional, and it exists only under certain conditions. If the constitutive equation is hyperelastic, that is if there exists a strain energy function Vˆ (␬) such that m ˆ (␬) ⫽ dVˆ / d␬, then Vainberg’s theorem guarantees the existence of the energy functional I⬅

冕

ᐉ

[mw⬙⫺qw⫺␬m ⫹ Vˆ (␬)]dx

(41)

0

that has the property that DI·(w¯ ,m ¯ ,␬¯ ) ⫽ J, i.e., the directional derivative of the energy I in the direction of the virtual fields (w¯ ,m ¯ ,x¯ ) is J. Thus, the virtual work functional J is simply a three⫺field variational principle that has this particular relationship to the Hu–Washizu functional I.䊐 5.3. Classical formulation based on the Hu–Washizu type functional If we substitute interpolations of displacement, moment, and curvature i.e., w ⫽ hTw, m ⫽ bTm, ␬ ⫽ cT␬ (where the interpolation functions c and the generalized curvatures ␬ can be understood in the same sense as the other interpolations) then the discrete form of the functional J can be written as ˆ c(␬)] J⫽w ¯ T[C␬⫺BTw]⫺␬¯ T[CTm⫺m ¯ T[Bm⫺q]⫺m

(42)

where B and q are defined in Eq. (33) and (10) and the new terms are defined as C⬅

冕

ᐉ

ˆ c(␬) ⬅ bcTdx, m

0

冕

ᐉ

ˆ (␬)cdx m

(43)

0

ˆ c(␬) to remind us that the residual depends upon where we have used the notation m the curvature field which, of course, is determined from the curvature parameters ␬ through the interpolation. The necessary condition J ⫽ 0 for all m ¯ ,␬¯ , and w¯ implies the discrete equations ˆ c (␬) ⫽ 0 Bm ⫽ q, BTw⫺C␬ ⫽ 0, ⫺CTm ⫹ m

(44)

Again, it is interesting to note that only the last of these three equations in nonlinear. ˆ c(␬)⬇m ˆ c(␬i) ⫹ Di⌬␬i about the state ˆ c(␬) linearizes to m The nonlinear function m with curvature field ␬i, where Di is given by Di ⬅

冕

ᐉ

cDicTdx

0

ˆ / d␬ is the tangent modulus evaluated at the state ␬i. where Di ⬅ dm

(45)

984

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

Linearization of these equations is very similar to the two-field mixed formulation described in the last section. The incremental equations can be put in the form = ⫺1 ¯ i wi ˆ c(␬i)]⫺K (46) Ki⌬wi+1 ⫽ q ⫹ BA⫺1 i C[␬i⫺Di m = T ⫺1 T where Ki ⬅ BA⫺1 i B is the effective stiffness matrix with Ai ⬅ CDi C . The question of invertibility of the matrix Ai is central to the success of the mixed method (see, for example, [13]). It is obvious, for example, that there must be more terms in c than there are in b (i.e., the curvature interpolation must be of higher order than the moment interpolation). However, we will not dwell on this generic discretization of the Hu–Washizu type functional, but rather turn to the so-called enhanced strain methods. 5.4. Enhanced strain methods based on the Hu–Washizu type functional The method that issues from the Hu–Washizu type functional leaves considerable room for choice in the interpolation functions. A closer inspection of the problems in selecting the interpolation functions leads to improved methods. The basic ideas behind making these choices evolved from the methods of incompatible modes that goes back at least to Wilson et al. [14]. Simo and Hughes [11] describe the formal variational structure that applies to these methods and discuss the idea of variationally consistent recovery of the stresses in these methods. Simo and Rifai [6] further elaborate this class of methods and show the connections among the various historical contributions to the subject. Let us take the three-field functional of Eq. (40) as a point of departure. The mixed, enhanced strain methods make the observation that the curvature field of the beam can be expressed as the sum of a compatible part (i.e., the part given by the appropriate derivatives of the displacement field) and an enhanced part as follows

␬ ⬅ w⬙ ⫹ ␬e

(47)

where ␬ is the ‘enhanced’ curvature field. We will also assume that the virtual curvature field admits the same decomposition. To wit, ␬¯ ⬅ w¯ ⬙ ⫹ ␬¯ e. If we substitute these expressions into the basic functional of Eq. (40) we find the new functional e

E(w,m,␬e,w¯ ,m ¯ ,␬¯ e) ⬅

冕

ᐉ

ˆ (␬)[w¯ ⬙ ⫹ ␬¯ e]⫺qw¯ ⫺(␬em [m ¯ ⫹ m␬¯ e)]dx

(48)

0

We can now observe that if the real and virtual enhanced strain fields are chosen to be orthogonal (in some sense) to the real and virtual moment fields then the last term in the functional will vanish. That is, we require that

冕

ᐉ

(␬em ¯ ⫹ m␬¯ e)dx ⫽ 0

(49)

0

for all ␬¯ e and m ¯ . This rather enigmatic requirement was proposed by Simo and Hughes [11] and is discussed in the context of Bernoulli–Euler beams by Petrangeli and Ciampi [5].

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

985

Now let us interpolate the fields as w ⫽ hTw, m ⫽ bTm, ␬e ⫽ cT␬e (where the interpolation functions c and the generalized enhanced curvatures ␬e are understood as the suitable modifications of those introduced in the previous section) then the discrete form of the functional E can be written as ˆ h(w,␬e)⫺q]⫺␬¯ eTm ˆ c(w,␬e) E⫽w ¯ T[m

(50)

where ˆ h(w,␬e) ⬅ m

冕

ᐉ

ˆ c(w,␬e) ⬅ ˆ (␬)h⬙dx, m m

0

冕

ᐉ

ˆ (␬)cdx m

(51)

0

It should be understood that the curvature that appears as the argument of the ˆ (␬) in these integrals is evaluated from Eq. (47) with the appropriate function m ˆ h(w,␬e) and m ˆ c,(w,␬e) interpolations as ␬ ⫽ h⬙Tw ⫹ cT␬e. Hence the notation m reminds us that these matrices depend upon the displacement and enhanced curvature fields. Eq. (50) implies the nonlinear equations ˆ c(w,␬e) ⫽ 0 ˆ h(w,␬e) ⫽ q, m m

(52)

which linearize about the state {wi,␬ } to e i

ˆ h(wi,␬ei ) Ki⌬wi ⫹ Hi⌬␬ei ⫽ q⫺m

(53)

ˆ c(wi,␬ei ) HTi⌬wi ⫹ Di⌬␬ei ⫽ ⫺m where a subscript i indicates that the matrix is evaluated at the state {wi, ␬ei ). The matrix Ki is the ordinary tangent stiffness matrix of the stiffness approach, defined in Eq. (10). The coefficient matrix Di is defined in Eq. (45). The matrix Hi is defined as Hi ⬅

冕

ᐉ

h⬙DicTdx

(54)

0

ˆ / d␬ is the material tangent stiffness of the element where, in each case, Di ⬅ dm evaluated at the curvature ␬i ⫽ h⬙Twi ⫹ cT␬ei . From the second equation we can find the generalized enhanced strains in terms of the displacements as T ⫺1 ˆ c(wi,␬ei ) ⌬␬ei ⫽ ⫺D⫺1 i Hi ⌬wi⫺Di m

(55)

This result can be substituted into the first equation to give T ˆ h(wi,␬ei ) ⫹ HiD⫺1 ˆ c(wi,␬ei ) [Ki⫺HiD⫺1 i Hi ]⌬wi ⫽ q⫺m i m

(56)

It is interesting to note that the orthogonality condition, Eq. (49), can be expressed in discrete form as m ¯ TC␬e ⫹ ␬¯ eTCTm ⫽ 0

(57)

where the matrix C was defined in Eq. (43) to be the outer product of b and c. It is fairly obvious that the only way to guarantee the orthogonality condition implied by Eq. (49) is to have the functions b and c orthogonal to each other on the domain

986

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

[0, ᐉ]. If these functions are chosen to be orthogonal, then C ⫽ 0 and Eq. (57) is satisfied automatically. One of the ramifications of forcing C ⫽ 0 is that the moment field does not enter into the equations governing the response of the beam, except through its influence on the choice of interpolation of the other two fields. However, we generally need to know the moment field and thus it must be recovered from the enhanced curvatures and displacements. One approach to this recovery is simply to compute the curvature field from the interpolation of displacements and enhanced curvatures as ␬ ⫽ h⬙Tw ⫹ cT␬e and then use this curvature in the constitutive equation m ⫽ ˆ (␬). As pointed out by Simo and Hughes [11], this approach is not variationally m consistent. In fact, it ignores the requirement given by Eq. (49). Simply stated, the moment m ˆ (␬) does not lie within the interpolation m ⫽ bTm for any set of genˆ (␬) is not a eralized moments m (e.g., if b represents a linear interpolation and m linear function). One solution to this problem is to find the values m that give the closest projection ˆ (␬) on the basis b. Specifically, we wish to satisfy the overdetermined set of of m equations ˆ ( ␬) bTm ⫽ m

(58)

The simplest solution to this problem is to find the least-squares projection. Premultiplying Eq. (58) by b and integrating over the length of the beam results in the following least-squares projection ˆ b(wi,␬ei ) Pm ⫽ m

(59)

where P⬅

冕

ᐉ

ˆ b(w,␬e) ⬅ bbTdx, m

0

冕

ᐉ

ˆ (␬)bdx m

(60)

0

ˆ ⫽ D␬. If we Remark 8. It is interesting to note that for the linear problem m premultiply Eq. (58) by bD⫺1 and integrate over the length, then the moment field is recovered as m ⫽ F⫺1BTw

(61)

where F is the flexibility matrix of Eq. (17). This estimate of the moment parameters is exactly the same expression that appears in the Hellinger–Reissner type formulation. This approach was proposed by Simo and Rifai [6] and is applicable only in the linear problem.䊐 5.5. Mixed-enhanced strain methods As a modification of the enhanced strain methods proposed by Simo and co-workers, Kasper and Taylor [15] have proposed a class of mixed-enhanced strain methods based on the Hu–Washizu type functional. These methods were proposed for three dimensional low-order finite element interpolations with the objective of avoiding

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

987

volumetric locking and to increase the performance of these elements in bending dominated problems. The ideas can be easily etched onto the Bernoulli–Euler beam considered here. Take the three-field functional, Eq. (40), as a point of departure. Let the displacement and moment fields be interpolated as w ⫽ hTm and m ⫽ bTm, respectively, and let the curvature field be interpolated by two parts as ␬ ⫽ bT␬o ⫹ cT␬e. The main difference with the previous approach is that there is no ‘compatible part’ of the curvature field. Hence, the enhanced curvature parameters do not necessarily represent curvature in addition to that given by the appropriate derivatives of the displacement field. Let the virtual fields be interpolated in exactly the same way as the real fields. Substituting these interpolations into Eq. (40), gives the discrete version of the functional ˆ b(␬o,␬e)] M⫽w ¯ T[P␬o ⫹ C␬e⫺BTw]⫺␬¯ oT[Pm⫺m ¯ T[Bm⫺q]⫺m

(62)

ˆ c(␬ ,␬ )] ⫺␬¯ [C m⫺m eT

T

o

e

ˆ c, are defined as ˆ b and m where the vector functions m ˆ b(␬o,␬e) ⬅ m

冕

ᐉ

ˆ c(␬o,␬e) ⬅ ˆ (␬)bdx, m m

0

冕

ᐉ

ˆ (␬)cdx m

(63)

0

The matrices B, C, and P have been defined previously in Eqs. (33, 43) and (60), respectively. Setting M ⫽ 0 for all w ¯ ,␬¯ o, and ␬¯ e, and applying the fundamental ¯ ,m theorem of the calculus of variation, we get the following governing equations ˆ b(␬o,␬e), CTm ⫽ m ˆ c(␬o,␬e) Bm ⫽ q, P␬o ⫹ C␬e ⫽ BTw, Pm ⫽ m

(64)

For convenience we can take the functions b and c to be orthogonal (this time it is a choice, not a requirement) to give C ⫽ 0. The first and second equations are linear. The second equation is satisfied by ␬o ⬅ ⌫Tw

(65)

where ⌫ ⬅ BP⫺1 (note that P is symmetric). Thus, the curvature ␬o can be thought of as being a function of the displacement w, as one might expect from the developments of the previous section. The moment parameters m can be eliminated by substituting the third equation into the first. Thus, the nonlinear equations that must be solved are ˆ c(␬o,␬e) ⫽ 0 ˆ b(␬o,␬e) ⫽ q, m ⌫m

(66)

These equations can be linearized to give ˆ b(␬oi,␬ei ) [⌫Ji⌫T]⌬wi ⫹ [⌫Gi]⌬␬ei ⫽ q⫺⌫m

(67)

ˆ c ( ␬ ,␬ ) [G ⌫ ]⌬wi ⫹ Di⌬␬ ⫽ ⫺m T i

T

e i

o i

e i

where we have substituted ⌬␬oi ⬅ ⌫t⌬wi, and the coefficient matrices are defined as Ji ⬅

冕

ᐉ 0

bDibTdx, Gi ⬅

冕

ᐉ

0

bDicTdx

(68)

988

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

For the system of Eq. (67) one again has the choice of condensing the enhanced curvature parameters ␬e at the element level to generate a tangent stiffness and right hand side for displacements only T T ˆ i ⬅ ⌫[Ji⫺GiD⫺1 K i Gi ]⌫ (69) ˆ b(␬o,␬e)⫺GiD⫺1 ˆ c(␬o,␬e)] rˆ i ⬅ q⫺⌫[m i m ˆ i⌬wi ⫽ rˆ i, for the displaceThe Newton iteration can be carried out by solving K ment increment ⌬wi, updating the displacement as wi ⫹ 1 ⫽ wi ⫹ ⌬wi. The enhanced curvature increment ⌬␬ei can then be found by solving Eq. (67)2 and the enhanced curvature can be updated as ␬ei ⫹ 1 ⫽ ␬ei ⫹ ⌬␬ei . This process is repeated until the displacement and enhanced curvatures converge. With these converged values, the moment parameters can be computed directly from Eq. (64)3 as mi ⫹ 1 ⫽ ˆ b(␬oi ⫹ 1,␬i ⫹ 1e). Note that this estimate of the moments turns out to be exactly P⫺1m the least-squares estimate proposed by Simo and Rifai [6]. Remark 9. It is interesting to note that if D is constant, then J ⫽ DP and G ⫽ 0 (by orthogonality of b and c). The flexibility matrix is F ⫽ P / D. Thus, ⌫J⌫T ⫽ [BP⫺1][DP][P⫺1BT] ⫽ BF⫺1BT

(70)

which is the flexibility matrix that appeared naturally from the Hellinger–Reissner type formulation. As mentioned previously, if b is linear and h is cubic, then ˆ ⫽ ⌫J⌫T is the exact stiffness computed by the classical flexibility formulation.䊐 K Remark 10. Often, the interpolation represented by b and h⬙ span the same space. For example, if h represents the cubic hermitian beam functions and b is linear, then both b and h⬙ span the space of linear functions. One can then write h⬙ ⫽ ⌳b for some constant matrix ⌳. Using this relationship for h⬙ in the definition of B, we find that B ⫽ ⌳P. Therefore, we can observe the equality ⌫ ⫽ BP⫺1 ⫽ ⌳. Thus, the matrix ⌫, which appears in the mixed⫺enhanced formulation, can be formed without the evaluation of any integrals. Furthermore, this observation implies that the variationally consistent stress recovery of Kasper and Taylor [15] is exactly equivalent to the least⫺squares projection proposed by Simo and Rifai [6], at least for this problem. 䊐 5.6. Global analysis methods The methods based upon the Hellinger–Reissner and Hu–Washizu type functionals can be extended to the analysis of structures as assemblies of elements. The global functionals of the Hellinger–Reissner and Hu–Washizu types are given as

冘 冘 E

HG ⬅

¯ He(w,m,w¯ ,m ¯ )⫺RTW

e⫽1 E

JG ⬅

¯ Je(w,m,␬,w¯ ,m ¯ ,␬¯ )⫺RTW

e⫽1

(71)

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

989

where the functional He and Je for each element are given by Eq. (30) and Eq. (40), ¯ are the concatenation of the nodal force and virtual disrespectively, and R and W placement components. Because these functionals both include the terms mw¯ ⬙⫺qw¯ , they share the feature with the displacement-based approach that nodal equilibrium is strongly enforced simply by selecting displacement interpolations that are continuous at the nodes (i.e., by sharing global degrees of freedom). Hence, global structural analysis based on these functionals includes direct assembly of the tangent stiffness matrix and residual identical to the pure displacement formulation. Therefore, these methods can be easily implemented in a computer program designed for carrying out displacement-based calculations. There are some important observations to make about the mixed formulations. In accord with the fundamental theorem of the calculus of variations, only the displacement field needs to have any interelement continuity. The moment and curvature fields do not, in fact often should not, be continuous between elements. Two simple examples illustrate. First, for a node where more than two elements meet, it makes no sense to even suggest the notion of continuity of the moment field. The end forces must satisfy equilibrium. Second, for a node where there is an abrupt change in cross section, the curvatures should not be continuous (in this case the moments should be continuous, but equilibrium implies the continuity). Because the moments and curvatures do not need to meet continuity requirements, they can be condensed out at the element level, as described by Eqs. (36) and (46). Condensation at the element level involves the inversion of matrices, which could add significantly to the computational expense of element stiffness and residual matrices. One always has the option of retaining the generalized moments and curvatures as global variables and assembling the complete system of equations for each element.

6. Example Ultimately, the distinction between methods—and hence the decision as to which method is best—mainly comes down to accuracy, computational efficiency, and ease of implementation. These features of an algorithm are not so straightforward to compare and some authors do better than others in making these comparisons. There is not enough room in this paper to make a reasonable comparison of estimates so we simply provide a benchmark example to demonstrate the features of the different approaches. Consider the inverted cantilever beam of unit length subjected to a unit tip load, as shown in Fig. 3 The linearly varying depth gives a bending modulus that varies as D(x) ⫽ (1 ⫹ x)3. The moment field is m(x) ⫽ 1⫺x, while the curvature field has the expression ␬(x) ⫽ (1⫺x) / (1 ⫹ x)3. This example, although irritating as a structural form, is an excellent simple benchmark problem for the methods described in this paper because of the great differences among the displacement, moment, and curvature fields. Fig. 4 shows the results of analyzing this beam with several different approaches. The displacement errors shown are the difference between the computed and exact

990

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

Fig. 3. Simple benchmark example that has linear moment field but a highly nonlinear curvature field even for a linearly elastic material.

values, divided by the exact end displacement. The moment and curvature errors are the difference between computed and exact values, divided by the exact value at the fixed end. Fig. 4(a–c) gives the results for the classical stiffness approach with fourth, fifth, sixth, and seventh order polynomial interpolation. The difficulties of the classical stiffness approach for low-order interpolations are clearly manifested in the struggle of the same basis try to fit a linear moment and a nonlinear curvature. It is evident that with higher-order interpolation the classical stiffness approach gives good results. Fig. 4(d–i) show the results for the two-field and three-field approaches, respectively. There are two cases done for the two-field approach. The 6-function results have a cubic displacement field and a linear moment field. The 8-function results have a quintic displacement field and a quadratic moment field. As expected, the moment field is exact, but the computed curvature is less accurate. There are also two cases done for the three-field approach. The 7-function result has a cubic displacement field, a linear moment field, and a single quadratic term for the enhanced curvature field. The 8-function result adds a cubic term to the enhanced curvature field. Note that the displacement error does not go down for the 8-function approximation because the displacement field remains cubic, unlike the 8-function twofield approach. The problems associated with comparing these methods should be evident from this example. The total number of parameters is one measure of comparison. It is evident in this example that the classical stiffness approach competes very favorably with the two- and three-field approximations. The amount of computational effort for each method is quite different, however. A more detailed study of these tradeoffs is left to a future paper.

7. Summary The Hellinger–Reissner and Hu–Washizu type functionals provide the natural context for developing mixed methods of structural analysis that capitalize on the benefits of the classical methods based on the principle of virtual forces without the

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

991

Fig. 4. Distribution of displacement, moment, and curvature errors for the cantilever benchmark example problem. (a–c) Classical Stiffness [䊐, 䊊,+ and 왕 for 4, 5, 6 and 7 base functions, respectively], (d–f) Hellinger–Reissner [䊐, 䊊 for 6 and 8 base functions, respectively], (g–i) Hu–Washizu [䊐, 䊊 for 7 and 8 base functions, respectively].

problems associated with the selection of redundant forces. These methods are attractive for beam problems in which the curvatures and moments have vastly different gradients. These problems commonly arise for inelastic materials in which plastic zones with high curvature gradients develop. Different moment and curvature gradients also occur for nonprismatic members, even with elastic response. In this paper we have presented the theoretical underpinnings of the variational structure associated with mixed methods for the particular case of Bernoulli–Euler beam theory along with some historical perspective from the classical methods. We

992

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

have demonstrated how different methods of nonlinear structural analysis result from applying simple interpolation schemes to the various field variables that appear in the multifield functionals. Through linearization procedures, we have shown how to formulate the element stiffness matrices and element residuals that are appropriate to each functional and have examined some of the connections among those quantities. We have also compared the variational approaches with the so-called nonlinear flexibility methods that have recently been reported in the literature and conclude that these approaches, while certainly having merit, are not variationally consistent. In particular, they are only loosely related to the Hellinger–Reissner type two-field functional. Questions of implementation and performance of the various methods described herein still remain. The tradeoffs between the classical stiffness method and the mixed methods presented here are many. It remains unclear if it is more efficient to compute with a classical stiffness approach with h- or p-refinement (e.g., the PEP element described by Chan and Zhou [16]) or to use a mixed method that introduces additional variables for the moment and curvature approximations. These questions are addressed in a forthcoming paper [17], which presents the formulation of the mixed methods in a more general setting and presents numerical examples to illustrate the trade-offs. References [1] Ciampi V, Carlesimo L. A nonlinear beam element for seismic analysis of structures. In: Proceedings of the 8th European Conference on Earthquake Engineering. Laboratorio National de Engenharia Civil, Lisbon (Portugal): 1986. [2] Spacone E, Ciampi V, Filippou FC. Mixed formulation of nonlinear beam finite element. Comp Struct 1996;58:71–83. [3] Spacone E, Filippou FC, Tauter FF. Fibre beam-column model for nonlinear analysis of RC frames: formulation. Earthquake Engng Struct Dyn 1996;25:711–25. [4] Neuenhofer A, Filippou FC. Evaluation of nonlinear frame finite elements models. J Struct Engng ASCE 1997;123(7):958–66. [5] Petrangeli M, Ciampi V. Equilibrium based iterative solutions for the non-linear beam problem. Int J Numer Meth Engng 1997;40:423–37. [6] Simo JC, Rifai MS. A class of mixed assumed strain methods and the method of incompatible modes. Int J Numer Meth Engng 1990;29:1595–638. [7] Hjelmstad KD. Fundamentals for structural mechanics. Upper Saddle River (NJ): Prentice Hall, 1997. [8] Gelfand IM, Fomin SV. Calculus of variations. Englewood Cliffs (NJ): Prentice Hall, 1963. [9] Argyris JH, Kelsey S. Energy theorems and structural analysis. London: Butterworths, 1960. [10] Heath MT, Plemmons RJ, Ward RC. Sparse orthogonal schemes for structural optimization using the force method. SIAM J Sci Stat Comput 1984;5(9):514–32. [11] Simo JC, Hughes TJR. On the variational foundations of assumed strain methods. J Appl Mech ASME 1986;53:51–4. [12] Simo JC, Kennedy JG, Taylor RL. Complementary mixed finite element formulations for elastoplasticity. Comp Meth Appl Mech Engng 1989;74:177–206. [13] Zienkiewicz OC, Taylor RL. The finite element method, 1. London: McGraw-Hill, 1989. [14] Wilson EL, Taylor RL, Doherty WP, Ghaboussi J. Incopatible displacement models. In: Fenves SJ, editor. Numerical and computer models in structural mechanics. New York: Academic Press; 1973. [15] Kasper EP, Taylor RL. A mixed-enhanced strain method: part 1—geometrically linear problems. Comp. Struct. 2000;75:237–50.

K.D. Hjelmstad, E. Taciroglu / Journal of Constructional Steel Research 58 (2002) 967–993

993

[16] Chan SL, Zhou ZH. Pointwise equilibrating polynomial element for nonlinear analysis of frames. J Struct Engng ASCE 1994;120(6):1703–17. [17] Hjelmstad KD, Taciroglu E. Mixed variational methods for finite element analysis of geometrically nonlinear, inelastic Bernoulli–Euler beams. Comm Num Meth Engng (in press).