TL bar element

. 8 The TL Bar Element 8–1 Chapter 8: THE TL BAR ELEMENT 8–2 TABLE OF CONTENTS Page §8.1. Introduction 8–3 §8.2. The Two-Dimensional Bar Element 8–3 §8.2.1. Element Kinematics . . . . . . . . . . . . . . . 8–3 §8.2.2. Strain Measure . . . . . . . . . . . . . . . . 8–4 §8.2.3. Stress Measure . . . . . . . . . . . . . . . . . 8–7 §8.2.4. Total Potential Energy and Residual Force Equations . . . . 8–8 §8.2.5. The Tangent Stiffness Matrix . . . . . . . . . . . . 8–8 §8.3. FEM Coding Using Mathematica 8–10 §8.3.1. The Example Structure . . . . . . . . . . . . . . 8–10 §8.3.2. Forming the Internal Force . . . . . . . . . . . . . 8–11 §8.3.3. The Equilibrium Equations . . . . . . . . . . . . 8–14 §8.3.4. Plotting the Equilibrium Paths . . . . . . . . . . . . 8–15 §8.3.5. Having Some Fun: Snapshots and Animation . . . . . . 8–15 §8.3.6. Forming the Tangent Stiffness Matrix . . . . . . . . . 8–18 §8.3.7. Critical Point Study . . . . . . . . . . . . . . . 8–18 §8. Exercises . . . . . . . . . . . . . . . . . . . . . . 8–25 8–2 8–3 §8.2 THE TWO-DIMENSIONAL BAR ELEMENT §8.1. Introduction In the present Chapter the key concepts of nonlinear continuum mechanics reviewed in Chapter 7 are applied to the development of governing equations of bar (truss) elements based on the Total Lagrangian (TL) kinematic description. For brevity these will be referred to as “TL bar elements.” There are two ways to construct TL elements: 1. The Standard Formulation (SF) 2. The Core Congruential Formulation (CCF). The first method is easier to describe and will be presented in this Chapter through examples. The second one is more flexible and powerful but it is more difficult to teach because it proceeds in stages. As such it will be relegated to Chapters 10-11. §8.2. The Two-Dimensional Bar Element The element developed in this Chapter is a prismatic bar element that can be used to model pin- jointed plane truss structures of the type sketched in Figure 8.1. These structures undergo large displacements and rotations but their strains are assumed to remain small so that the material behavior stays in the linear elastic range. These assumptions allows us to consider only geometric nonlinear effects. A two-node bar element appropriate to model members of such truss structures is shown in Figure 8.2. The element moves in the (X, Y ) plane. In the reference (base) configuration the element has cross section area A0 (constant along the element) and length L0. In the current configuration the cross section area and length become A and L , respectively. The material has an elastic modulus E that links the axial-stress and axial-strain measures defined below. Because thisChapter deals primarilywith the formulation of an individual element, the identification superscript (e) will be omitted to reduce clutter until assemblies are considered. The element has four node displacements and associated node forces. These quantities are collected in the vectors u =   uX1 uY1 uX2 uY2   , f =   fX1 fY1 fX2 fY2   , (8.1) The loads acting on the nodes will be assumed to be conservative. §8.2.1. Element Kinematics In accordance with bar theory, to describe the element motion it is sufficient to consider a generic point of coordinates X located on the longitudinal axis of the reference configuration C0. That point maps to point x in the current configuration C. The bar remains straight in any configuration. These coordinates can be parametrically interpolated from the end nodes as X(ξ) = [ X (ξ) Y (ξ) ] = [ 1 2 (1 − ξ)X1 + 12 (1 + ξ)X2 1 2 (1 − ξ)Y1 + 12 (1 + ξ)Y2 ] x(ξ) = [ x(ξ) y(ξ) ] = [ 1 2 (1 − ξ)x1 + 12 (1 + ξ)x2 1 2 (1 − ξ)y1 + 12 (1 + ξ)y2 ] (8.2) 8–3 Chapter 8: THE TL BAR ELEMENT 8–4 Motion Reference configuration� (same as base) Current configuration Figure 8.1. A plane truss structure undergoing large displacements while its material stays in the linear elastic range. Here ξ is the usual isoparametric coordinate that varies from −1 at node 1 to +1 at node 2. The displacement field is obtained by subtracting the foregoing position vectors: u(ξ) = x(ξ) − X(ξ) = [ uX (ξ) uY (ξ) ] = [ 1 2 (1 − ξ)uX1 + 12 (1 + ξ)uX2 1 2 (1 − ξ)uY1 + 12 (1 + ξ)uY2 ] , (8.3) which expressed in matrix form is u(ξ) = [ uX (ξ) uY (ξ) ] = [ 1 2 (1 − ξ) 0 12 (1 + ξ) 0 0 12 (1 − ξ) 0 12 (1 + ξ) ]  uX1 uY1 uX2 uY2   = N(ξ) u. (8.4) The element kinematic defined by these equations is depicted in Figure 8.3. §8.2.2. Strain Measure As discussed in Chapter 7, in the Total Lagrangian (TL) description the Green-Lagrange (GL) strains and the second Piola-Kirchhoff (PK2) stresses are frequently used as conjugate measures in the formulation of the internal energy. The only GL strain that appears in the energy expression is 8–4 8–5 §8.2 THE TWO-DIMENSIONAL BAR ELEMENT X, x Y, y X1 u u u u Y1 X2 Y2 f X1 fY1 f X2 fY2 1(X1, Y1) 2(X2, Y2) 1(x1, y1) 2(x2, y2) Reference configuration Area A , length L C0 PK2 stress ,s0 GL strain e0 = 0 Current configuration C PK2 stress ,s GL strain e 0 0 Area A, length L � Figure 8.2. The geometrically nonlinear, two-node, two-dimensional bar element in Total Lagrangian description. This element may be used to model members of a plane truss such as that shown in Figure 8.1. the axial strain e1 ≡ e which is most expediciously defined using the length change as e = L 2 − L20 2L20 , (8.5) rather than through displacement gradients. Because of the linear displacement assumptions (8.4) the strain e is constant over the element.1 This expression can bemaneuvered into a matrix function of the node displacements as follows. Let X21 = X2−X1, Y21 = Y2−Y1, aX = X21/L0, aY = Y21/L0, uX21 = uX2−uX1, uY21 = uY2−uY1, umX = (uX2 − uX1)/2 and umY = (uY2 − uY1)/2. Some of these quantities can be geometrically interpreted as illustrated in Figure 8.4. Then 1 This is in fact the only use of the displacement interpolation (8.4) in the following derivations. 8–5 Chapter 8: THE TL BAR ELEMENT 8–6 X X , x Y, y 1(X1, Y1) 2(X2, Y2) 1(x1, y1) 2(x2, y2) X(ξ) = { X (ξ) Y (ξ) } = { 1 2 (1 − ξ)X1 + 12 (1 + ξ)X2 1 (1 − ξ)Y + 1 (1 + ξ)Y } 2 1 2 2 x(ξ) = { x(ξ) y(ξ) } = { 1 2 (1 − ξ)x1 + 12 (1 + ξ)x2 1 2 (1 − ξ)y1 + 12 (1 + ξ)y2 } u(ξ) = x(ξ)−X(ξ) = { uX (ξ) uY (ξ) } = { 1 2 (1 − ξ)vX1 + 12 (1 + ξ)vX2 1 2 (1 − ξ)vY1 + 12 (1 + ξ)vY2 } Bar longitudinal axis C0 C _ Figure 8.3. The definition of displacement field for the two-dimensional TL bar element. L2 = (X21 + uX21)2 + (Y21 + uY21)2, e = L 2 − L20 2L20 = 1 L0 (aXuX21 + aY uY21) + 1 2L20 ( u2X21 + u2Y21 ) = 1 L0 [−aX −aY aX aY ]u + 1 L20 [−umX −umY umX umY ]u = (Bl + Bn(u))u. (8.6) Observe that the GL strain e has been separated into two parts: e = el + en , where 1. el = Blu, where Bl is constant, depends linearly in the node displacements u. This is called the linear part of the strain. 2. en = Bnu, in which Bn is a function of the node displacements, depends quadratically on the displacements. This is called the nonlinear part of the strain, because Bn(u) vanishes if u → 0. Note that both matrices are constant over the element. 8–6 8–7 §8.2 THE TWO-DIMENSIONAL BAR ELEMENT X, x Y, y 1(X1, Y1) 2(X2, Y2) C0 C X21 Y21 L0 L ψ ψ0 aX = X21 L0 = cosψ0 aY = Y21 L = sinψ0 0 ay = y21 L0 = L L0 sinψ1(x1 = X1+ u uX1, y1 = Y1+ Y1) u u2(x2 = X2+ X2, y2 = Y2+ Y2) ux21 = X21 + X21 uy21 = Y21 + Y21 ax = x21 L0 = L L0 cosψ X _ Figure 8.4. Geometric interpretation of quantities used in the study of element kinematics. The variation of e induced by variations δ u in the node displacements is δe = Bl δu + δ(Bnu) = B δu, (8.7) The matrix B that links δe to δu is B = ∂(B� + Bnu) ∂u = 1 L0 [−ax −ay ax ay ] (8.8) in which ax = aX + uX21 L0 = x2 − x1 L0 = x21 L0 , ay = aY + uY21 L0 = y2 − y1 L0 = y21 L0 . (8.9) For a geometric interpretation of ax and ay see Figure 8.4. 8–7 Chapter 8: THE TL BAR ELEMENT 8–8 §8.2.3. Stress Measure The stress measure conjugate to GL strains is the second Piola-Kirchhoff (PK2) stress tensor. The only component that appears in the internal energy is the axial stress s, which is related to e through the constitutive equation s = s0 + Ee, (8.10) where s0 is the axial stress in the reference configuration, and E is the elastic modulus. The axial force based on this stress is N = A0s (8.11) Note that this is not the true axial force in the current configuration C , which would be Ntrue = Aσ (8.12) where σ is the true or Cauchy stress in C and A is the actual area. §8.2.4. Total Potential Energy and Residual Force Equations In what follows it is assumed that the element is subjected only to node forces f that are conservative and proportional, so that f = λq. The Total Potential Energy (TPE) of the element in the current configuration is � = U − P = ∫ V0 (s0e + 12 Ee2) dV0 − fT u = ∫ L0 A0(s0 e + 12 Ee2) d X¯ − λqT u. (8.13) where X¯ is directed along the bar longitudinal axis in C0, as shown in Figure 8.4. This energy expression is separable because the internal energy depends only on u through e and not on λ. The finite element residual equations are obtained by differentiating this equation or, equivalently, making (8.13) stationary with respect to virtual displacement variations δ u: δ� = δU − δP = ( p − f)T δu = 0. (8.14) The first variation of U is given by δU = p T δu = ∫ L0 A0(s0 δe + Eeδe) d X¯ = ∫ L0 A0 s δe d X¯ = ∫ L0 N0B δu d X¯ = NL0B δu (8.15) where N = A0s is the axial force in the current configuration measured per area of the reference configuration.2 It follows that the internal force vector is p = N0L0BT = N   −ax −ay ax ay   . (8.16) This equation admits of a simple geometric interpretation; see Figure 8.5. The relation between N0 = A0s and the true axial force N = Aσ can be worked out from inspection of this diagram. The load potential variation δP simply generates f = λq as can be expected. 2 This is the PK2 axial force; cf (8.11). 8–8 8–9 §8.2 THE TWO-DIMENSIONAL BAR ELEMENT 1 2 N0 fX1 = −N0ax fY1 = −N0ay fX2 = N0ax fY2 = N0ay Figure 8.5. Geometrical interpretation of the internal force vector. The axial force N0 = A0s would be positive as shown. §8.2.5. The Tangent Stiffness Matrix Because the residual equations are separable the tangent stiffness matrix is obtained simply by differentiating the internal force with respect to the node displacements u: K = ∂p ∂u = ∂(NL0B T ) ∂u = A0L0BT ∂s ∂u + A0L0s ∂B T ∂u = KM + KG . (8.17) The above expression shows that K splits naturally into two parts: KM and KG , which are called the material stiffness matrix and geometric stiffness matrix, respectively, in the FEM literature. To get KM note that ∂s ∂u = ∂(s0 + Ee) ∂u = E ∂e ∂u = EB. (8.18) Consequently KM = E A0L0BT B. (8.19) Inserting the expression (8.8) for B yields KM = E A0L0   a2x axay −a2x −axay axay a 2 y −axay −a2y −a2x −axay a2x axay −axay −a2y axay a2y   (8.20) This component of K looks formally similar to the stiffness matrix of a linear bar element,3 except that B now depends on u. The dependence of KM on the material properties (here the elastic modulus E) explains the name “material stiffness” given in the FEM literature. 3 To which it reduces if u = 0. In that case ax and ay become the sine and cosine of the angle ψ0 shown in Figure 8.4. 8–9 Chapter 8: THE TL BAR ELEMENT 8–10 The other component can be obtained by differentiating B with respect to the node displacements, the result being a constant 4 × 4 matrix: ∂BT ∂u = 1 L20   1 0 −1 0 0 1 0 −1 −1 0 1 0 0 −1 0 1   . (8.21) Inserting this into (8.17) one gets KG = NL0   1 0 −1 0 0 1 0 −1 −1 0 1 0 0 −1 0 1   . (8.22) This component ofK depends only on the stress state in the current configuration, because N = A0s. No material properties appear. Thus the name “geometric stiffness” applied to KG .4 Remark 8.1. Assuming that E , A0 and L0 are nonzero, the rank of KM is obviosly one because B is a 1 × 4 matrix. On the other hand the rank of the numerical matrix in (8.21) is 2 (because its eigenvalues are 2, 2, 0, 0). Consequently KG has rank 2 if s is nonzero and 0 otherwise. Combining these results it can be shown that the rank of K = KM +KG is 1 if the configuration is stressed and 2 otherwise. In other words, the rank deficiency is 3 and 2, respectively. The implications of this property in the analysis of stability are studied later. Remark 8.2. The addition of KG increases the bar stiffness if the current configuration is in tension (s > 0), but it reduces it if the current configuration is in compression (s < 0). This is in accord with physical intuition. The main effect of this stiffness is on the rotational rigid-body motions of the bar about the Z axis. §8.3. FEM Coding Using Mathematica This section illustrates the use of Mathematica as a quick way to do some simple problems with the TL bar element. The use of a computer symbolic system as prototyping language has several advantages: (i) The code is very compact, because of higher lever array notation. Because of its compactness, it can be easily translated into other high-level languages such as, for example, Matlab or Maple. (ii) Debugging is often quicker than with programming languages because of the interpretative nature of the language front end. This is of course platform dependent; the implementation of Mathematica on a Macintosh or PC is friendler than on a Unix workstation. (iii) Graphic output is part of the language. On the Macintosh, graphics can be transported into documents such as this one via PostScript and Adobe Illustrator. (iv) The language can manipulate both algebraic (symbolic) and numerical expressions. The later may be done exactly or in floating-point arithmetic. This flexibility is convenient for many situations, as the examples below illustrate. 4 In the pre-1970 FEM literature, the name “initial stress stiffness” was used for KG by some authors. 8–10 8–11 §8.3 FEM CODING USING MATHEMATICA �� �� �� �� α S H X, x Y, y (2)(1) 1 2 3 uX uY fY = λ E, A0E, A0 Here λ > 0 means crown� load is applied upwards Figure 8.6. An arch as a two-bar FE model. This model has two degrees of freedom: uX and uY , which are the displacements of node 2. This structure will be used as example for use of Mathematica in §8.3. §8.3.1. The Example Structure The two-spring arch studied in Exercise 6.3 will be used as a structure that illustrates the use of Mathematica in conjunction with finite elements. The structure is shown in Figure 8.6. It is modeled with two TL bar elements labeled (1) and (2), which join nodes 1-2 and 2-3, respectively. The modulus of elasticity and the bar areas are unity. The arch geometry is defined by the span S and height H . The reference configuration is assumed stress free. Because nodes 1 and 3 are pinned to the ground, the only degrees of freedom are the horizontal and vertical displacements of node 2. To avoid clutter these will be denoted in the sequel as uX and uY rather than uX2 and uY2, respectively. The vertical motion uY is positive upwards. The node-2 components of the applied and internal force vector will be also denoted simply by ( fX , fY ) and (pX ,pY ), respectively. Displacement and force components for nodes 1 and 3 will be explicitly removed from the governing equations. The arch is loaded by a vertical force fY = λ applied at node 2, positive upwards. The applied horizontal force fX will be always zero. Recall that in Exercise 6.3 only symmetric motions of the arch under a vertical load were considered. That constraint effectively reduces the model to one degree of freedom: uY , and allows only limit point behavior. When two degrees of freedoms are considered, however, a richer set of possibilities emerges; notably lateral bifurcation under a vertical load. §8.3.2. Forming the Internal Force The Mathematicamodule FormIntForce2DTwoNodeBar listed in Figure 8.7 forms the internal force vector of an individual 2D bar finite element handled by the TL description. The arguments of this module are XY1 Coordinates {X1, Y1} of end node 1 in C0. XY1 Coordinates {X2, Y2} of end node 2 in C0. uX1 Displacements {uX1, uY1} of end node 1. 8–11 Chapter 8: THE TL BAR ELEMENT 8–12 FormIntForce2DTwoNodeBar[XY1_,XY2_,uXY1_,uXY2_,Em_,A0_,s0_]:= Module[{X1,Y1,X2,Y2,X21,Y21,uX1,uY1,uX2,uY2,uX21,uY21,L0,L, e,s,ax,ay,ifv}, {X1,Y1}=XY1; {X2,Y2}=XY2; X21=X2-X1; Y21=Y2-Y1; {uX1,uY1}=uXY1; {uX2,uY2}=uXY2; uX21=uX2-uX1; uY21=uY2-uY1; L0=Sqrt[X21^2+Y21^2]; L=Sqrt[(X21+uX21)^2+(Y21+uY21)^2]; e=(L-L0)*(L+L0)/(2*L0^2); s=s0+Em*e; ax=(X21+uX21)/L0; ay=(Y21+uY21)/L0; ifv=A0*s*{{-ax}, {-ay}, {ax}, {ay}}; Return[Simplify[ifv]] ]; pe=FormIntForce2DTwoNodeBar[{2,3},{5,7},{1,0},{1,0}, 20,12, 0]; Print["Internal elem force pe=",pe]; pe=FormIntForce2DTwoNodeBar[{2,3},{5,7},{1,0},{1,0}, 20,12, s0]; Print["Internal elem force pe=",pe]; Internal elem force pe={{0}, {0}, {0}, {0}} -36 s0 -48 s0 36 s0 48 s0 Internal elem force pe={{------}, {------}, {-----}, {-----}} 5 5 5 5 Figure 8.7. A Mathematica module that computes the internal force vector p of an individual 2-node 2D bar element. Output from two test statements is shown. uX2 Displacements {uX2, uY2} of end node 2. Em Elastic modulus E . A0 Cross sectional area A0. s0 PK2 stress in C0 Two test statements, one with purely numeric arguments and the other one with one symbolic argument (s0 for the bar prestress), are shown in the figure along with their output. Figure 8.8 shows two modules, AssembleMasterIntForceOfTwoBarArch and MergeElemIntoMasterIntForce that, together with FormIntForce2DTwoNodeBar, form the internal force vector p for the two-bar example structure. The arguments of the former are the element internal force vector eif, the element freedom table eftab that maps element-level freedom to master-level freedoms and the master internal force vector mif before the element is merged. The function returns the updated internal force vector. The arguments of the second module are: span Span S height Height H . Em Elastic modulus E . A0 Cross sectional area A0. uX,uY Crown displacements uX , uY . 8–12 8–13 §8.3 FEM CODING USING MATHEMATICA MergeElemIntoMasterIntForce[pe_,eftab_,p_]:= Module[{i,ii,neldof,fmaster}, fmaster=p; neldof=Dimensions[eftab][[1]]; For[i=1, i<=neldof, i++, ii=eftab[[i]]; If [ii>0,fmaster[[ii,1]]+=pe[[i,1]]] ]; Return[fmaster] ]; AssembleMasterIntForceOfTwoBarArch[span_,height_,Em_,A0_,uX_,uY_]:= Module[{f1,f2,fmaster}, fmaster=Table[0,{2},{1}]; f1=FormIntForce2DTwoNodeBar[{-span/2,0},{0,height}, {0,0},{uX,uY},Em,A0,0]; fmaster=MergeElemIntoMasterIntForce[f1,{0,0,1,2},fmaster]; f2=FormIntForce2DTwoNodeBar[{0,height},{span/2,0}, {uX,uY},{0,0},Em,A0,0]; fmaster=MergeElemIntoMasterIntForce[f2,{1,2,0,0},fmaster]; Return[Simplify[fmaster]] ]; ClearAll[Em,A0,S,H,uX,uY,f]; p=AssembleMasterIntForceOfTwoBarArch[2,2.5,10,.75,-0.4,0.25]; Print["Master Int Force p= ",p]; p=AssembleMasterIntForceOfTwoBarArch[S,H,Em,A0,uX,uY]; Print["Master Int Force p= ",p]; Master Int Force p= {{-0.5336499821957073}, {1.555758744891104}} Master Int Force p= {{(4*A0*Em*uX* (S^2 + 2*uX^2 + 4*H*uY + 2*uY^2))/(4*H^2 + S^2)^(3/2)}, {(8*A0*Em*(H + uY)*(uX^2 + 2*H*uY + uY^2))/ (4*H^2 + S^2)^(3/2)}} Figure 8.8. Two Mathematica modules that assemble the internal force vector p for the two-bar arch structure of Figure 8.6. Two test statements, one numeric and the other symbolic, are shown along with output. Two test statements are also shown in the figure. The first one is purely numeric and sets the following properties for the structure: S = 2, H = 2.5, E = 10, A0 = 0.75, uX = 0.4, uY = −0.25 (8.23) and the module AssembleMasterIntForceOfTwoBarArch returns p = [ pX pY ] = [−0.53365 1.55576 ] (8.24) The second one specifies symbolic values for all of the arguments: S for span S, H for the height H , etc. The result (dis