Theory Manual Volume 2

Theory Manual Volume 2LUSAS Version 14 : Issue 2 LUSAS Forge House, 66 High Street, Kingston upon Thames, Surrey, KT1 1HN, United Kingdom Tel: +44 (0)20 8541 1999 Fax +44 (0)20 8549 9399 Email: [email protected] http://www.lusas.com Distributors Worldwide Table of Contents Table of Contents 7 Element Formulations .................................................................................1 7.1 Bar Elements (BAR2, BAR3, BRS2, BRS3) ........................................... 1 7.1.1 Formulation.......................................................................................1 7.1.2 Evaluation and Output of Stresses/Forces....................................... 1 7.1.3 Nonlinear Formulation ......................................................................2 7.2 Beam Elements.......................................................................................7 7.2.1 2-D Straight Beam (BEAM) ............................................................10 7.2.2 2-D Straight Grillage (GRIL) ...........................................................12 7.2.3 2-D Ribbed Plate Beam (BRP2).....................................................14 7.2.4 3-D Straight Beam (BMS3).............................................................17 7.2.5 2-D Curved Thin Beam (BM3, BMX3) ............................................21 7.2.6 3-D Curved Thin Beam (BS3, BS4, BSX4) ....................................29 7.2.7 Semiloof Thin Beam (BSL3, BSL4, BXL4) .....................................37 7.2.8 3-D Straight Beam (BTS3)..............................................................45 7.3 Two-Dimensional Continuum Elements................................................55 7.3.1 Standard Isoparametric Elements ..................................................55 7.3.2 Enhanced Strain Elements (QPM4M, QPN4M, QAX4M)...............66 7.3.3 Incompatible Plane Membrane Element (PMI4).............................73 7.3.4 2D Explicit Dynamics Elements......................................................76 7.3.5 Two Phase Plane Strain Continuum Elements (TPN6P and QPN8P) ...................................................................................................87 7.3.6 Large-strain Mixed-type Elements (QPN4L, QAX4L).....................90 7.4 Three-Dimensional Continuum Elements .............................................95 7.4.1 Standard Isoparametric Elements (HX8, HX16, HX20, PN6, PN12, PN15, TH4, TH10)...................................................................................95 7.4.2 Enhanced Strain Element (HX8M) ...............................................102 7.4.3 3D Explicit Dynamics Elements (HX8E, PN6E, TH4E) ................ 105 7.4.4 Composite Solid Elements (HX8L,HX16L,PN6L PN12L).............113 7.4.5 Two Phase 3D Continuum Elements (TH10P, PN12P, PN15P, HX16P and HX20P)...............................................................................118 7.5 Space Membrane Elements................................................................122 7.5.1 Axisymmetric Membrane (BXM2, BXM3).....................................122 7.5.2 3-D Space Membrane (SMI4, TSM3)...........................................126 7.6 Plate Elements ....................................................................................128 7.6.1 Isoflex Thin Plate (QF4, QF8, TF3, TF6)......................................128 7.6.2 Isoflex Thick Plate (QSC4) ...........................................................135 7.6.3 Isoparametric Thick Mindlin Plate (QTF8, TTF6) .........................139 7.6.4 Ribbed Plate (RPI4, TRP3) ..........................................................144 7.7 Shell Elements ....................................................................................148 7.7.1 Axisymmetric Thin Shell (BXS3) ..................................................148 7.7.2 Flat Thin Shell (QSI4, TS3) ..........................................................155 7.7.3 Flat Thin Shell Box (SHI4) ............................................................159 i Table of Contents 7.7.4 Semiloof Thin Shell (QSL8, TSL6) ...............................................163 7.7.5 Thick Shells (TTS3, TTS6, QTS4, QTS8) ....................................176 7.8 Field Elements.....................................................................................189 7.8.1 Thermal Bar (BFD2, BFD3) ..........................................................189 7.8.2 Thermal Axisymmetric Bar (BFX2, BFX3) ....................................191 7.8.3 Thermal Link (LFD2, LFS2, LFX2) ...............................................192 7.8.4 Plane Field (QFD4, QFD8, TFD3, TFD6) .....................................194 7.8.5 Axisymmetric Field (QXF4, QXF8, TXF3, TXF6) .........................198 7.8.6 Solid Field (HF8, HF16, HF20, PF6, PF12, PF15, TF4, TF10) ....200 7.8.7 Solid Composite Field (HF8C, HF16C, PF6C, PF12C)................204 7.9 Joint Elements.....................................................................................208 7.9.1 Joints (JNT3, JPH3, JF3, JRP3, JNT4, JL43, JSH4, JL46, JSL4, JAX3, JXS3) ..........................................................................................208 7.9.2 Evaluation of Stresses/Forces......................................................209 7.9.3 Nonlinear Formulation ..................................................................209 7.9.4 Use of Joints With Higher Order Elements...................................212 7.10 Fourier Element Formulation (TAX3F, QAX4F, TAX6F, QAX8F) .....215 7.10.1 Global and Local Coordinate Systems .......................................215 7.10.2 Standard Isoparametric Elements ..............................................215 7.10.3 Strain-Displacement Relationships.............................................217 7.10.4 Constitutive Relationships ..........................................................218 7.10.5 Element Loading.........................................................................219 7.10.6 Inertial Loading ...........................................................................221 7.10.7 Evaluation of Stresses ................................................................225 7.11 Interface Elements (INT6, INT16) .....................................................225 7.11.1 Definition and interpolation .........................................................225 7.11.2 Internal force vector and stiffness matrix....................................226 Appendix A ..................................................................................................229 Quadrature Rules ......................................................................................229 Appendix B ..................................................................................................239 Restrictions On Element Topology............................................................239 Mid-Length and Mid-Side Nodes ...........................................................239 Warping of Flat Elements ......................................................................239 References...................................................................................................241 ii Notation Notation Standard matrix notation is used whenever possible throughout this manual and the expressions are defined as follows: Basic Expressions Vector Matrix or second order tensor Fourth order tensor : Matrix scalar product Determinant of a matrix Norm of a vector Trace of a matrix Transpose of a vector of matrix Inverse of a matrix Variation Virtual variation Rate Increment Summation | | || || bg bg bg db g δb g ej Δb g ∑b g tr T −1 • diag , Diagonal matrix with terms given <,> Orthogonality condition Variables defined in local axes ej ∧ iii y.Y. velocities and accelerations Component i Maximum value Normal component in slideline analyses Initial components (initial strains) Thermal components (thermal strains) Tangential component in slideline analyses Components in the local x.Z Z .z Cartesian system Components in the global X.y.Theory Manual 2 Subscripts cr g i max n o t x x. C2 Mooney-Rivlin constants c c Da Cohesion in friction based material models Wave speed Maximum distance between two adjacent contact nodes iv .Y.z X.Z Cartesian system Zienkiewicz constants Differentiation with respect to x (or other variable) Iteration I Local quantities (co-rotational for continuum elements) Values at time t Values at time t+Δt Rayleigh damping coefficient (multiplies mass matrix M) Bulk modulus Rayleigh damping coefficient (multiplies stiffness matrix K) neo-Hookean constant Scalars aR B bR C0 C1 .x Superscripts i l t t + Δt Critical value Ground displacements. I3 I1.Notation E fi Young’s Modulus Slideline interface force on contact node i Yield surface Initial gap for nonlinear joint models Shear Modulus Fracture energy in concrete model Transfer coefficient in field analysis First stress invariant Strain invariants Modified strain invariants Volume ratio (det F ) Second deviatoric stress invariant Third deviatoric stress invariant Interface stiffness coefficient Bulk modulus Thermal conductivity Spring stiffness when in contact for nonlinear joint models Spring stiffness after liftoff for nonlinear joint models Length of local contact segment Initial chord length of beam element Current chord length of beam element Moment Stress resultant Number of required eigenvalues Participation factor in spectral response analysis Axial force F g G Gf h I1 I1. I 2 J J2 J3 k k K Kc K1 l lo ln M N p P P v .I2. Theory Manual 2 q q Q Q hg Field variable flux in field analysis Number of starting iteration vectors for subspace iteration Rate of internal heat generation in field analysis Hourglass constant Constants for shock wave smoothing Contact zone radius Spectral displacement in special response analysis Spectral velocity in spectral response analysis Spectral acceleration in spectral response analysis Thickness of local contact segment Period of oscillation for transient and dynamic analysis Temperature in structural applications Torque Axial stretch Element volume Crack width in concrete model Work Normal penetration distance Radial overlap constant Coefficient of thermal expansion Softening Parameter in concrete model Constant used in dynamic recurrence algorithms Constant used in dynamic recurrence algorithms Shear retention factor in concrete model Constant used in dynamic recurrence algorithms Displacement norm used for convergence Q1.Q2 rz Sd Sv Sa t T T T u V w W X α α α α β β γ γd vi . αp ν ξ ρ ρ σ vii .Notation γψ γw γ ψ1 γ ψ2 εi η Residual norm used for convergence Work norm used for convergence Root mean square of residuals convergence criterion Maximum absolute residual convergence criterion Error estimate in subspace iteration Step length multiplier for line search Lode’s Angle Angle between old and new displacement vector in arch-length method Angle of orthotropy Angle defining crack directions in concrete model Local slope at a node Strain hardening parameter Load factor Plastic strain rate multiplier Eigenvalue (ith) Principal stretches Eigenvalue shift in subspace iteraction Friction coefficient Ogden constants Poisson’s ratio Modal Damping ratio Model coefficient for CQC combination in spectral response analysis Mass density Effective stress e e e e e κ λ λ λi λi μ μ μp. g .Theory Manual 2 φ φ φ Φ ψ ω Interface stiffness scale factor Friction angle in friction based models Structural damping in harmonic response analysis Field variable in field analysis Potential energy Circular frequency in transient and dynamic analysis Circular frequency of load in harmonic response analysis Local spin at the centroid Nodal displacement vector Unit vectors forming the co-rotated base axes Green-Lagrange strain vector Vector of master slideline surface forces Vector of nodal body forces ξ η ζ Ω Ω1 Vectors a ei E f f g . g Covariant base vectors Mm n P P Vector of master slideline surface mass Vector of unit segment normal Global internal force vector Local internal force vector Euler parameters External force vector Unit vectors defining the beam cross section at a gauss point Deviatoric Cauchy stress Second Piola-Kirchhoff stress vector Vector of surface tractions ~ q R ri s S t viii . 5.ω i ψ Ω1 Matrices/Tensors A Matrix of slopes B Strain-displacement matrix Linear strain-displacement matrix Displacement dependent strain-displacement matrix Damping matrix in dynamic analysis Matrix of constrain constants Green deformation tensor Compliance matrix of material moduli B B C C C C 0 1 ix . q x Y ε θ θ ψ i Unit vectors defining the beam cross section at a node Vector of unit segment tangent Generalised displacement vector in spectral response analysis Logarithmic strain vector Displacement gradient vector (geometric nonlinearity) Unscaled pseudovector (co-rotational formulation .Notation ti .1) Pseudovector of rotation Lagrangian multiplier vector Incompatible modes for enhanced elements Cauchy stress vector Jaumann variation of Cauchy stress Eigenvector Residual force vector Local spin at the centroid λ λi σ dσ J Φ i .section 3. Theory Manual 2 D D F G K Rate of deformation tensor Material modulus matrix Deformation gradient matrix Matrix of shape functions Stiffness matrix T K K Tangent stiffness matrix Stress stiffness matrix Mass matrix Shape function array σ M N N Q R i Principal directions of the Lagragian triad Vector of constraint constants Rotation tensor Skew symmetric matrix of the vector θ Matrix containing Second Piola-Kirchoff stresses Transformation matrix for co-rotational formulation Right stretch tensor Angular acceleration tensor Matrix of eigenvalues Local engineering strain tensor Density matrix Matrix containing Cauchy stresses Matrix containing Cauchy stresses Sθ $ S T U α Λ ε af ρ $ σ $ $ σ x . Notation σ Φ τ Ω Biot stress tensor Matrix of eigenvectors Kirchhoff stress tensor Angular velocity tensor xi . Theory Manual 2 xii . 1.1. The consistent and lumped mass matrices are evaluated using the procedures defined in (section 2.the axial force. BAR3. BRS2. BRS3) 7. 7.1 Bar Elements (BAR2.2 Evaluation and Output of Stresses/Forces The element output can be obtained at both the element nodes and Gauss points and consists of Fx . BAR3.the axial strain.7). Appropriate references are included when full details of the element derivation are not provided.1-1).7.Bar Elements (BAR2. (tension +ve) ∈x .1 Formulation The bar elements are 2-node and 3-node isoparametric elements that can only transmit longitudinal force (fig. BRS2. 7. The nodal variables are:BAR2 and BAR3 U and V BRS2 and BRS3 U. V and W The element strain-displacement relationship and thermal strain vector are defined in the local Cartesian system as ∈x = ∂u ∂x and b∈ g = αT o t The elastic constitutive relationship is defined as σ x = E ∈x A complete description of the element stiffness formulation is given in [B1]. BRS3) 7 Element Formulations This section of the Theory Manual covers the basic theoretical assumptions made for each element formulation. (tension +ve) 1 . For a curved element the local xy-plane is defined by the element nodes (fig.7. For a straight element not parallel to the global x-axis. The local y-axis forms a right-hand set with the local x and z-axes for all three cases.1-4).1-3). The geometric nonlinearity is a Total Lagrangian formulation which accounts for large displacements but small strains. BRS2 and BRS3 elements The local x-axis lies along the element axis in the direction in which the element nodes are specified (for a curved element it is tangent to the curve at the point concerned).Element Formulations The forces and strains are output in the local element coordinate system defined by BAR2 and BAR3 elements The element local x-axis lies along the element axis in the direction in which the element nodes are specified (fig. 7. Local y is perpendicular to local x and +ve on the convex side of the element.7.3 Nonlinear Formulation The bar elements can be employed in Materially nonlinear analysis utilising the elastoplastic constitutive laws [O1] (section 4.1-4). The local y and z axes form a right-hand set with the x-axis such that the y-axis lies in the global XY-plane and the z-axis is parallel to the global Z-axis (up out of page). Geometrically and materially nonlinear analysis utilising the nonlinear material laws specified in 1.1. Geometrically nonlinear analysis.2). Linear eigen-buckling analysis. Nonlinear dynamics utilising the nonlinear material laws specified in 1. The nonlinear strain-displacement relationship is defined by BAR2 and BAR3 ∈x = ∂u 1 ∂u + ∂x 2 ∂x LM OP N Q 2 + 1 ∂v 2 ∂x LM OP N Q 2 BRS2 and BRS3 2 .7. the local z-axis is defined by the unit vector z = j x x where j is a unit vector defining the Global Y-axis and x is a unit vector defining the local x-axis (fig. Note. the local z-axis is defined by the unit vector z = i x x where i is a unit vector defining the global X-axis (fig.7. For a straight element parallel to the global x-axis.1-4). BAR3. The loading is conservative. Y BAR2 2 V V 1 U V 1 U 2 U 3 U BAR3 X V V U (a) 2-D Bar Elements Y BRS2 W V U 1 W V W U V V U 2 V 3 W 2 U BRS3 X U 1 W Z (b) 3-D Bar Elements 3 . BRS3) ∈x = ∂u 1 ∂u + ∂x 2 ∂x LM OP N Q 2 + 1 ∂v 2 ∂x LM OP N Q 2 + 1 ∂w 2 ∂x LM OP N Q 2 with reference to the local x-axis. BRS2. referred to the undeformed configuration.Bar Elements (BAR2. The forces and strains output with the geometrically nonlinear analysis will be the 2nd Piola-Kirchhoff forces and Green-Lagrange strains respectively. 7.1-3 Local Cartesian System For BAR2 And BAR3 Elements 4 .1-2 Examples Illustrating Use Of BAR Elements 2 Y y 1 x 1 2 x y y x 3 x y X Fig 7.1-1 Nodal Freedoms For BAR Elements Struts represented with BAR2 elements Pressure Continuum elements 2-D Roof Truss Excavation Supports Fig.Element Formulations Fig.7. BRS2. BAR3. BRS3) Y x-y plane y y x y z 1 X x 3 z 2 z x Z (a) Curved Element Y z x 1 y 2 y X z x Z (b) Straight Element Parrallel With Global X-axis y 2 z y x 1 z X Y x Z (c) Arbitrarily Orientated Straight Element Fig. 7.Bar Elements (BAR2.1-4 Local Cartesian System For BRS2 And BRS3 Elements 5 . Element Formulations 6 . M y1.Beam Elements 7. w 2 . θ y 2 . M y 2 . Py 2 . M z1. v1. The nodal forces/moments and degrees of freedom (in local coordinates) for the 3D beam are F T = Px1. w1. Pz 2 . M x1. M z 2 a T = u1. Px 2 . Pz1. v2 . θ x1. u 2 . Py1. θ z1. M x 2 .2 Beam Elements The family of explicit straight beams are derived by restraining various degrees-offreedom of the full 3D beam. θ x 2 . The stiffness and mass matrices of these reduced elements may be obtained by deleting the appropriate rows and columns of the full stiffness and mass matrices. θ z 2 The corresponding stiffness and mass matrices are Element stiffness matrix K= LMK MNK 11 12 K 21 K 22 OP PQ OP PP PP PP PP PP b4 + Φ gEI PP Lb1 + Φ g d4 + Φ iEI PP 0 Ld1 + Φ i P Q Symmetric z y z y z y where submatrices are defined: K 11 LM EA MM L MM 0 MM 0 =M MM 0 MM 0 MM MM 0 N L3 1 + Φ y 0 0 0 12 EI z d i 12 EI y L 1 + Φz 0 −6 EI y L 1 + Φz 2 3 b g g GJ L 0 0 b L2 1 + Φ y d 6 EI z i 0 7 . θ y1. Element mass matrix M = ρAL LMM NMM 11 12 M 21 M 22 OP QP where submatrices are defined 8 .Element Formulations K 22 LM EA MM L MM 0 MM 0 =M MM 0 MM 0 MM MM 0 N L 1 + Φy 0 0 0 −6 EI z 3 12 EI z d i 12 EI y L3 1 + Φ z 0 L2 1 + Φ z 6 EI y b g g GJ L 0 0 b L 1 + Φy 2 d i 0 0 OP PP PP PP PP PP b4 + Φ gEI PP Lb1 + Φ g d4 + Φ iEI PP 0 Ld1 + Φ i P Q Symmetric z y z y z y K 12 = KT 21 LM −EA MM L MM 0 MM 0 =M MM 0 MM 0 MM MM 0 N 0 0 0 L 1 + Φy 0 0 0 3 −12 EI z d i 0 −12 EI y L3 1 + Φ z 0 −6 EI y L2 1 + Φ z b g g 0 − GJ L 0 0 b L2 1 + Φ y d 6 EI z i 12 EI y 0 OP PP −6 EI 0 L d1 + Φ i P PP 6 EI 0 PP L b1 + Φ g PP 0 0 PP b2 − Φ gEI 0 Lb1 + Φ g P d2 − Φ iEI PP 0 Ld1 + Φ i P Q 0 0 2 z y y 2 z z y z y z y and where Φy = 12 EI z GAs y L2 and Φz = GAs z L2 As y and As z are the cross-sectional areas effective in shear about the respective bending axis. a f ρAL 2 M 2.3 = a f ρAL 2 9 . M 11 = .Beam Elements M 11 LM 1 3 MM 0 MM M0 =M MM 0 MM 0 MM MN 0 LM 1 3 MM 0 MM MM 0 = MM 0 MM 0 MM MN 0 13 6I z + 35 5AL2 0 0 0 11L I + z 210 10 AL − 6I y 13 + 35 5AL2 0 Iy 11L − 210 10 AL 0 Jx 3A 0 0 2I y L2 + 105 15A 0 OP PP PP PP PP PP P 2I P L + 105 15A P Q Symmetric 2 z 13 6I z + 35 5AL2 0 0 0 − 11L I − z 210 10 AL M 22 6I y 13 + 35 5AL2 0 Iy 11L + 210 10 AL 0 Jx 3A 0 0 2I y L2 + 105 15A 0 OP PP PP PP PP PP P 2I P L + 105 15A P Q Symmetric 2 z M 21 = MT 12 LM 1 MM 6 MM 0 MM 0 MM 0 MM 0 MM MN 0 0 9 6I z − 70 5AL2 0 0 0 − 13L I + z 420 10 AL 0 0 6I y 9 − 70 5AL2 0 Iy 13L − 420 10 AL 0 0 0 0 Jx 6A 0 0 − 0 0 Iy 13L + 420 10 AL 0 − Iy L − 140 30 A 0 2 OP I P 13L − P 420 10 AL P PP 0 PP 0 PP 0 PP L I P − − 140 30 A P Q 0 z 2 z The lumped mass matrix contains terms only the following terms.2 = a f ρAL 2 M 3. +ve forces and moments are in the directions of the positive local Cartesian system.7. such that the y-axis lies in the global XY plane.2. 2-noded straight beam formulated by superimposing the bending. used in engineering beam theory.7. The nodal forces F are evaluated directly using F = Ka in the local Cartesian system. The stress resultant variations are constant axial force.1-1) U. The forces are output in the local Cartesian system which is defined as having its local x-axis along the element axis in the direction in which the element nodes are specified.5 = a f ρI y L 2 M 6.2. 10 .4 = a f ρJ x L 2 M 5. linear moments and linear shear forces.1 2-D Straight Beam (BEAM) Formulation This element is a 2-D. Evaluation of stresses/forces The element output obtained at the nodes consists of Fx . and are evaluated explicitly using R U R F U = | EAαaLΔTfO | SM V SEI α M ΔT P V T W | N dy Q | T W e x e z zz where ( ΔT )e and ( ΔT dz)e are average element values.Element Formulations M 4. linear rotation and cubic transverse displacements. shear and axial behaviour derived directly from the differential equations for beam displacements. The local y and z-axes form a right-hand set with the x-axis. The nodal degrees of freedom are (fig.6 = a f ρI z L 2 7. and the z-axis is parallel to the global Z-axis (up out of page) (fig. Fy and M z .1-3). The nodal forces due to the thermal strains are assumed to be constant within each element.2. See [P1] for further element details. V and θz at each node The displacement variations along the length of the beam are linear axial. 7.1-2 Examples Illustrating Use Of Beam Elements 11 .2.1-1 Nodal Freedoms For BEAM Element Load Load Cantilever Beam Plane Frame Fig. Y V θz U U 2 V θz 1 X Fig. but may be utilised in a nonlinear environment. These values are evaluated by combining the nodal values with the local element forces and moments calculated explicitly.7.2.Beam Elements The local Cartesian forces may also be output at eleven equally spaced points along the beam. The element cannot be employed for linear buckling analyses. Nonlinear formulation The element does not possess any nonlinear capability. 2.2. The stress resultant variations are constant axial.7. and torsional behaviour derived directly from the differential equations for beam displacements used in engineering beam theory. 2-noded straight beam formulated by superimposing the bending. M x and M y +ve forces and moments are in the directions of the positive local Cartesian system.2. and are evaluated explicitly using M y = EI yyα LM ΔT OP N dz Q e where ( ΔT dz)e is the average element value. 12 . The nodal forces due to the thermal strains are assumed to be constant within each element. Evaluation of stresses/forces The element output obtained at the nodes consists of Fz . See [P1] for further element details.Element Formulations Y y x 1 2 X Fig. and linear moment and linear shear. θ x and θ y The displacement variations along the length of the beam are linear axial. linear rotation and cubic transverse displacements.2 2-D Straight Grillage (GRIL) Formulation This element is a 2-D.2-1) W. shear.1-3 Local Cartesian System For BEAM Element 7. The nodal degrees of freedom are (fig.7. These values are evaluated by combining the nodal values with the local element forces and moments calculated explicitly. such that the y-axis lies in the global XY plane. Nonlinear formulation The element does not possess any nonlinear capability. The local y and z-axes form a right-hand set with the x-axis.2-3).2. The nodal forces F are evaluated directly using F = Ka in the local Cartesian system.Beam Elements The forces are output in the local Cartesian system which is defined as having its local x-axis along the element axis in the direction in which the element nodes are specified.7.7.2-1 Nodal Freedoms For GRIL Element 13 . The local Cartesian forces may also be output at eleven equally spaced points along the bar. but may be utilised in a nonlinear environment.2. The element cannot be employed for linear buckling analyses. Y θy θy w Z 1 θx θx w X Fig. and the z-axis is parallel to the global Z-axis (up out of page) (fig. 2-noded. torsional and axial behaviour derived directly from the differential equations for beam displacements used in engineering beam theory.3 2-D Ribbed Plate Beam (BRP2) Formulation This element is a 2-D. and are evaluated explicitly using 14 . W.7. straight eccentric beam formulated by superimposing the bending.2-3 Local Cartesian System For GRIL Element 7. V. The stress resultant variations are constant axial. θ X and θY at each node The displacement variations along the length of the beam are linear axial.2-2 Example Illustrating Use Of GRIL Elements Y y x 2 1 X Fig.2.7.2. linear moment and linear shear. linear rotation and cubic transverse displacements.3-1) U. The nodal degrees of freedom are (fig. shear.2.Element Formulations Z Y X Point Load Y X Problem Defintion Finite Element Mesh Fig. The nodal forces due to the thermal strains are assumed to be constant within each element.2.7. These values are evaluated by combining the nodal values with the local element forces and moments calculated explicitly. The nodal forces F are evaluated directly using F = Ka in the local Cartesian system.3-3). and the z-axis is parallel to the global Z-axis (up out of page) (fig. Evaluation of stresses/forces The element output obtained at the nodes consists of Fx .2. M x .7. 15 . Fy . M y +ve forces and moments are in the directions of the positive local cartesian system. The element cannot be employed for linear buckling analysis. such that the y-axis lies in the global XY-plane. The local y and z-axes form a right-hand set with the x-axis. See [P1] for further element details. The local Cartesian forces may also be output at eleven equally spaced points along the beam. The forces are output in the local Cartesian system which is defined as having its local x-axis along the element axis in the direction in which the element nodes are specified. Fz . but may be utilised in a nonlinear environment. Nonlinear formulation The element does not possess any nonlinear capability.Beam Elements R F U = R EAαaΔΔTTf U SM V |EI α L O | T W S MN dz PQ V | | T W e x e y yy where (ΔT) and (ΔT/dy) are average element values. 7.3-2 Ribbed Plate Illustrating Use Of BRP2 Element 16 .2.7.Element Formulations Y V θy U 2 θx V θy 1 U θx W W Z X Fig.3-1 Nodal Freedoms For BRP2 Element RPI4 elements Z Y X Y BRP2 elements X Problem Definition Finite Element Mesh Fig.2. shear.2.2. See [P1] for further element details. θ Y and θ Z at each node The displacement variations along the length of the beam are linear axial. and are evaluated explicitly using R | EAαaΔTf U | R F U | L ΔT O | |M | = |EI α | S V S MN dz PQ V |M | | T W | L ΔT O | |EI α MN dy PQ | | | | T W e x e y z yy e zz where ( ΔT )e and ( ΔT dz)e are average element values.7.7.4 3-D Straight Beam (BMS3) Formulation This element is a 3-D two noded straight beam formulated by superimposing the bending.Beam Elements Y y x 2 1 X Fig. The stress resultant variations are constant axial. V. torsional and axial behaviour derived directly from the differential equations for beam displacements used in engineering beam theory.3-3 Local Cartesian System For BRP2 Element 7. 17 .4-1) U.2. constant torsion and linear moment and linear shear. The nodal forces due to the thermal strains are assumed to be constant within each element. linear rotation and cubic transverse displacements. The nodal degrees of freedom are (fig. θ X . W. These values are evaluated by combining the nodal values with the local element forces and moments calculated explicitly.7. M x . The element cannot be utilised for linear buckling analysis. Nonlinear formulation The element does not possess any nonlinear capability. The local xy-plane is defined by the third element node and the element x-axis. The local y and z-axes form a right-hand set with the local x-axis (fig. Fz Forces in the local Cartesian system. Fy .4-3). but may be utilised in a nonlinear environment. The local x-axis lies along the element axis in the direction in which the element nodes are specified. The nodal forces F are evaluated directly using F = Ka in the local Cartesian system.4-1 Nodal Freedoms For BMS3 Element 18 . M z Moments in the local Cartesian system.2. M y .2.Element Formulations Evaluation of stresses/forces The element output obtained at the nodes consists of Fx .7. The local Cartesian forces may also be output at eleven equally spaced points along the bar. Y θY V θY 1 W θZ U θX X Z W 2 θZ θX Fig. 4-2 Local Cartesian System For BMS3 Element 19 .Beam Elements 3 y x 2 Y z 1 X Z Fig.2.7. 2.4-3 Examples Illustrating The Use Of BMS3 Elements 20 .7.Element Formulations (a) 3-D Frame Structure (a) 3-D Frame Structure Fig. 2.5 2-D Curved Thin Beam (BM3. curved. The global displacements and rotations are initially quadratic and are interpolated independently using linear Lagrangian shape functions for the end nodes and a hierarchical quadratic function for the central node. non-conforming beam elements formulated using the constraint technique. Δu at the mid-length node where Δu is the local axial relative (departure from linearity) displacement.7. Δv.2.7. V. The Kirchhoff condition of zero shear strain is applied at the two integration points. the initial degrees of freedom are (fig. BMX3) Formulation The BM3 and BMX3 elements are thin.5-1) U. Δj at the end nodes at the mid-length node. The infinitesimal strain-displacement relationship is defined in the local system as ∈x = ∂u ∂x ∂2 v ∂x 2 Cartesian ψz = − The elastic rigidity (resultant modulus) and modulus matrices are defined as Explicit Numerically Integrated EA $ D= EI z D= LM N EI z EI zz z LMN h OP Q OP Q Eb Eyb dy Eyb Ey 2 b The thermal strain vector is defined as 21 . by forcing ∂v ∂u ∂v + = − θz = 0 ∂x ∂z ∂x and eliminating the local transverse translational and rotational degrees of freedom at the central node.Beam Elements 7. θz at the end nodes. Therefore. j Δu. V.5-1) U. The final degrees of freedom for the element are (fig.2. 5-6).5-7).7.7. The true nodal moments for a beam element between supports is then obtained by adding the fixed end moments to the end node values.flexural strain The forces and strains are output in the local x-axis which lies along the element axis in the direction in which the element nodes are specified.2. The moments are +ve for tension in the top fibre of the element (hogging). such that the y-axis lies in the global XY plane. The consistent and lumped mass matrices are evaluated using the procedures defined in (section 2. Three options for interpreting the forces and moments within an element are available The axial force and moment are computed at the two Gauss points using numerical integration. Force and stress output may be obtained at either the nodes or element Gauss points.7.axial force (+ve tension) . The the fibre lies on the +ve local y side of the element. The local y and z-axes form a right-hand set with the x-axis. Greatest accuracy is obtained at the Gauss points. and the sagging moment to the mid-node value (fig. The values at the centre point are then interpolated from these end values and the values at the Gauss points assuming a cubic variation (fig. and the z-axis is parallel to the global Z-axis (up out of page) (fig.5-7).moment . The axial force and moment are computed at the two end nodes by using F end = T T Ka where T is the global-local transformation matrix. This is the default technique and must be used for nonlinear analyses.S1].axial strain . 22 .7). Note.2.Element Formulations R U | | αΔT | | eψ j = S d ΔT | a f Lα + ΔT dα OV | | dy MN dT P W Q| T 0 t A complete description of the element formulation is given in [M1.2. Evaluation of stresses/forces The element output obtained at the nodes or Gauss points consists of Fx Mz ex ψz . This method can only be used for linear analyses and is invoked via OPTION 136. The nonlinear strain-displacement relationship is defined by ∈x = ∂u 1 ∂u + ∂x 2 ∂x LM OP N Q 2 + 1 ∂v 2 ∂x LM OP N Q 2 ψz = − ∂ 2 v ∂u ∂ 2 v ∂ v ∂ 2 u − + ∂x 2 ∂x ∂ x 2 ∂ x ∂ x 2 with reference to the local element x-axis.5-8. The loading is conservative. The force and strain output with the geometrically nonlinear analysis will be the 2nd Piola-Kirchhoff stress resultants and Green-Lagrange strains respectively. BMX3 may be used with the concrete model and continuum-based plasticity models (section 4. provided that the rotations are small within a load increment.2).2. The BMX3 elements are valid for rotations (TL) or rotation increments (UL) greater than one radian. Geometrically and materially nonlinear analysis utilising the nonlinear material laws specified in 1. 23 . The axial force distribution from a simple problem is given in fig.5). Geometrically nonlinear analysis. The loading approximates to being non-conservative. ∂u / ∂x may no longer be interpreted as axial strain. An Updated Lagrangian formulation takes account of large displacements and large rotations but small strains. Notes BM3 and BMX3 may be used in conjunction with the stress resultant plasticity model (section 4. Nonlinear dynamics utilizing the nonlinear material laws specified in 1.7.2). Nonlinear formulation The beam elements can be employed in Materially nonlinear analysis utilising the elastoplastic constitutive laws [O1] (section 4. As rotations become large. The geometric nonlinearity may be either A Total Lagrangian formulation which accounts for large displacements but small strains. The output approximates to the true Cauchy stress resultants and logarithmic strains.Beam Elements This method is similar to (b) except that the stress resultants at the centre node are also computed by considering equilibrium and is invoked via OPTION 137. Linear eigen-buckling analysis. The initial assumptions used in deriving the BM3 and BMX3 elements limit the rotations to one radian in a Total Lagrangian analysis and rotation increments of one radian in an Updated Lagrangian analysis (section 3.2). referred to the undeformed configuration. Element Formulations V θZ U ΔV Δθ Z 2 3 ΔU ΔU θZ V U 3 2 Y θZ V Y θZ U V 1 1 U X Initial Variables X Final Variables Fig.7.7.2.5-1 Nodal Freedoms For BM3 And BMX3 Elements Quadrature Points Quadrature points coincide with frame joints Fig.2.5-2 Portal Frame Showing Locations Of Quadrature Points With A 3-Point Newton-Cotes Rule 24 . Beam Elements Y 2 3 1 4 Z Fig.5-3 Local Cartesian Axes For Cross-Section Y Element 1 or Quadrilateral 1 2 1 2 3 3 4 Element 2 or Quadrilateral 2 Z Element 3 or Quadrilateral 3 2 1 1 4 3 4 Fig.7.2.5-4 Cross-Section Of I-Beam Represented By Superimposing Three BMX3 Elements Or By Defining Three Quadrilaterals 25 .2.7. 7.5-5 Quadrature Rules For Cross-Section Integration Y y 2 x y 1 X x y x 3 Fig.2.5-6 Local Cartesian System For BM3 And BMX3 Elements 26 .Element Formulations Y Y Z Z 3-Point Newton-Coates 5-Point Newton-Coates Fig.7.2. 5-7 Interpretation Of Results Obtained Using BM3 And BMX3 Elements 27 .Beam Elements UDL Support True moment distribution LUSAS + fixed end moments wl/12 wl/24 wl/12 Values evaluated at Gauss points and extapolated to nodes (a) Adding Fixed End Moments Nodal values computed directly from F = K a Gauss points values (b) Cubic Fit Through Gauss and Nodal Values Nodal values computed directly from F = K a Mid-point moment evaluated using equilibrium (c) Quadratic Fit Through Nodal and Mid-length Values Fig.7.2. Element Formulations Load (a) Problem Definition BM3 Axial Force (b) Axial Force Distribution Fig.7.5-8 Axial Force Distributions Obtained For A Geometrically Nonlinear Analysis Of A Cantilever Beam 28 .2. 7. Δθ Z at the end nodes at the mid-side node The Kirchhoff condition of zero shear strain is applied at the two integration points. W. ΔW. BS4. The final degrees of freedom for the element are (fig.2. BS4. curved. ΔV. non-conforming beam elements formulated using the constraint technique.6-1) U. θ Y .6 3-D Curved Thin Beam (BS3. θ X . BSX4) Formulation The BS3. θ Y . Δθ Y . W. θ X .2. θ Z ΔU.7. The global displacements and rotations are initially quadratic and are independently interpolated using linear Lagrangian shape functions for the end nodes and a hierarchical quadratic function for the central node.2. V. This provides C(0) continuity of the in-plane displacement. The infinitesimal strain-displacement relationship is ∈x = ∂u ∂x ∂2 u ∂x 2 ∂2 v ∂x 2 ∂2 w ∂x∂y ψy = − ψz = − ψ xy = − 29 . by forcing ∂v ∂u ∂v + = − θz = 0 ∂x ∂y ∂x ∂w ∂u ∂w + = + θy = 0 ∂x ∂z ∂x and eliminating the local transverse translational and bending rotational freedoms at the central node. Δθ X . V. θ Z Δu and Δθ X at the end nodes at the mid-side node where Δu and Δθ X are the local relative (departure from linearity) axial displacement and torsional rotation of the central node. and BSX4 elements are 3-D thin.Beam Elements 7. The initial freedoms are (fig.6-1) U. 30 . The consistent and lumped mass matrices are evaluated using the procedures defined in (section 2. Numerically integrated $ D= zz h b LM E MMEy MMEz 0 MM 0 MN 0 Ey Ey 2 Eyz 0 0 0 Ez 0 Eyz 0 Ez 2 0 0 Gy 2 0 0 0 0 0 0 0 0 Gz 2 0 0 0 0 dydz 0 0 G OP PP PP PP PQ The thermal strain vector is defined as R daΔTf αΔT dα U | Lα + ΔT O| | dz MN dT P | | daΔTf Lα + ΔT dα QO| V eψ j = S dy M | N dT P | Q| | 0 | | | | 0 T W 0 t A description of the element formulation is given in [M2]. if K t has a non-zero value in the element geometric properties data section. for circular cross-sections K t = J .7). ψ xy + ψ xz = ψ z the total torsional strain The elastic rigidity (resultant modulus) and modulus matrices are defined as Explicit LM EA MMEI EI $ D=M MM 0 MNM 0 0 y z EI y EI yy EI yz 0 0 0 EI z EI yz EI zz 0 0 0 0 0 0 GI yy 0 0 0 0 0 0 GI zz 0 0 0 0 0 0 GA OP PP PP PP QP Alternatively.Element Formulations ψ xz = − ∂2w ∂x∂y Note. the polar second moment of area). the resultant torsional moduli GI yy and GI zz are replaced with GK t / 2 where K t is a torsional constant (typically. 2. the local z-axis is given the unit vector z = i x x (i is a unit vector along the global X-axis) (fig. the local z-axis is given by the unit vector z = j x x (j is a unit vector along the global Y-axis) (fig.6-7c) The local y-axis forms a right-hand set with the local x and z axes. Geometrically and materially nonlinear analysis utilising the nonlinear material laws specified in 1. The local y-axis is perpendicular to the local xaxis and +ve on the side of the element where the fourth node lies.7. Tz ∈x Ψy .2.2. BS4. The local y and z-axis form a right-hand set with the local x-axis (fig. 31 . The torques are +ve for anti-clockwise rotations at first node and clockwise rotations at third node. Ψz Ψxy . Greatest accuracy is obtained at the Gauss points. Ψxz axial force moments torques axial strain flexural strain torsional strain The forces and strains are output in the local Cartesian system which is defined by BS3 For a curved element the local xy-plane is defined by the three element nodes. BSX4 The local xy-plane is defined by all four element nodes which are assumed to be coplanar.6-6) Note.Beam Elements Evaluation of stresses/forces The element output obtained at the nodes or Gauss points consists of Fx My .7. Force and stress output may be obtained at either the nodes or element Gauss points. For a straight element parallel to the global X-axis. Local y is perpendicular to local x and +ve on the convex side of the element.2). Mz Ty .6-7a).6-7b) For a straight element not parallel to the global X-axis.2.7.7. The local y and z-axis form a right-hand set with the local x-axis (Fig. Nonlinear formulation The beam elements can be employed in Materially nonlinear analysis utilising the elasto-plastic constitutive laws [O1] (section 4. Geometrically nonlinear analysis. The force and strain output for a geometrically nonlinear analysis will be 2nd Piola-Kirchhoff stress resultants and Green-Lagrange strains respectively. Linear eigen-buckling analysis. The geometric nonlinearity utilises a Total Lagrangian formulation which accounts for large displacements but small strains.2). 32 . All continuum based nonlinear material models do not consider nonlinear torsional effects.2). The initial assumptions in deriving the BS3. BSX4 may be used with the concrete model and continuum based plasticity models (section 4. Notes BS3. and rotation increments of one radian in an Updated Lagrangian (UL) analysis (Section 3. referred to the undeformed configuration.Element Formulations Nonlinear dynamics utilising the nonlinear material laws specified in 1. BS4 and BSX4 elements limit the rotations to one radian in a Total Lagrangian (TL) analysis. The nonlinear straindisplacement relationship is defined by ∈x = ∂u 1 ∂u + ∂x 2 ∂x LM OP N Q 2 + 1 ∂v 2 ∂x LM OP N Q 2 ψy = − ψz = − ∂ 2 w ∂u ∂ 2 w ∂w ∂ 2 u ∂w ∂ 2 v − + + ∂x 2 ∂x ∂x 2 ∂x ∂x 2 ∂y ∂x 2 ∂ 2 v ∂u ∂ 2 v ∂ v ∂ 2 u ∂ w ∂ 2 w − + − ∂x 2 ∂x ∂x 2 ∂x ∂x 2 ∂y ∂x 2 ∂ 2 w ∂u ∂ 2 w ∂w ∂ 2 v − − ∂x∂y ∂x ∂x∂y ∂x ∂x 2 ∂2 w ∂u ∂ 2 w ∂v ∂2 v − − ∂x∂y ∂x ∂x∂y ∂x ∂x 2 ψ xz = − ψ xy = − γ yz = − ∂v ∂w ∂x ∂x with reference to the local element x-axis.5). The loading is conservative. BS4 and BSX4 may be used in conjunction with the stress resultant plasticity model (section 4. 6-1 Nodal Freedoms For BS3. BS4 And BSX4 Elements Quadrature Points Quadrature points coincide with frame joints Fig.2.6-2 Portal Frame Showing Locations Of Quadrature Points With A 3-Point Newton-Cotes Rule 33 .Beam Elements V θY W ΔV Δθ Y 2 Δθ Z ΔU Δθ X 3 U θX W θZ V θY 3 U θX θZ Δu ΔW 2 Δθ X V θY Y W θZ 1 X Z Z U θX Y V θY 1 X U θX W θZ Initial Variables Final Variables Fig.7.7.2. 6-4 Cross-Section Of An I-Beam Represented By Superimposing Three Bsx4 Elements Or By Defining Three Quadrilaterals 34 .7.2.6-3 Local Cartesian Axes For Cross-Section Y Element 1 or Quadrilateral 1 2 1 2 3 3 4 Element 2 or Quadrilateral 2 Z Element 3 or Quadrilateral 3 2 1 1 4 3 4 Fig.7.Element Formulations Y 2 3 1 4 Z Fig.2. 7.6-6 Local Cartesian System For Bs4 And BSX4 Elements 35 .6-5 Quadrature Rules For Cross-Section Integration Y x y 1 x y X 4 3 x z y 2 Z Fig.Beam Elements 3*3 Newton-Cotes 5*5 Newton-Cotes Fig.7.2.2. Element Formulations Y x-y plane y y x y z 1 X x 3 2 z x Z (a) Curved Element Y z x 1 y X 2 3 y z x Z (b) Straight Element Parrallel With Global X-axis y 3 z y x 1 z X 2 Y x Z (c) Arbitrarily Orientated Straight Element Fig.7.2.6-7 Local Cartesian System For The BS3 Element 36 . which are located at the quadrature points of the 2 point Gauss rule (fig.2.BSL4 and BXL4 elements are 3-D thin. θ y . W. the displacements and rotations are interpolated using quadratic and cubic shape functions respectively. The final degrees of freedom for the element are (fig. by forcing ∂v ∂u ∂v + = − θz = 0 ∂x ∂y ∂x ∂w ∂u ∂w + = + θy = 0 ∂x ∂z ∂x which provides four constraint equations and permits elimination of the two flexural degrees of freedoms at these positions. W at nodes 1 and 3 at node 2 and θX at nodes 4 and 5 Note. but are not relative rotations (departures from linearity) as with the other LUSAS beam elements based on Kirchhoff constraints. curved beam elements based on the Kirchhoff constraint technique. θz U. θ X . BXL4) Formulation The BSL3. Their formulation and nodal configuration has been specifically designed to provide an element compatible with the Semiloof shell element QSL8.7-1). Unlike the thick beam formulation presented by Irons [I1].7-1) U. V. BSL4. the present formulation utilises Kirchhoff constraints of zero shear strain applied at the 2-point Gauss quadrature locations.7.7. Initially. V.2.Beam Elements 7.7 Semiloof Thin Beam (BSL3. The infinitesimal strain-displacement relationship is ∈x = ∂u ∂x ∂2 w ∂x 2 ∂2 v ∂x 2 ψy = − ψz = − 37 .2. where the cubic variation is provided by the rotational degrees of freedoms of the 'loof' nodes. The rotations at the 'loof' nodes are local. 38 .I1. Numerically integrated D= zz h b LM E MME MME0 MM 0 MN 0 z y Ez 2 Eyz Eyz Ey 2 0 0 0 0 0 0 Ez Ey 0 0 0 Gz 2 0 0 0 0 0 0 Gy 2 0 OP PP P dydz 0P 0P P GP Q 0 0 0 The thermal strain vector is defined by LM daΔTf LαΔT dα OP MM dz MNα + ΔT dT OPQPP Δ eψ j = MM dadyTf LMα + ΔT dα OPPP MM N 0 dT QPP 0 QP NM 0 t A more detailed description of the element formulation is given in [A1. the resultant torsional moduli GI yy and GI zz are replaced with GK t / 2 where K t is a torsional constant (typically. ψ xy + ψ xz = ψ z the total torsional strain The elastic rigidity (resultant modulus) and modulus matrices are defined as Explicit LM EA MMEI EI $ D=M MM 0 MM 0 N0 y z EI y EI yy EI yz 0 0 0 EI z EI yz EI zz 0 0 0 OP PP P 0 0 P GeI + Ae j P 0 Ge I + Ae j 0 P 0 0 GA P Q 0 0 0 0 0 0 0 0 0 yy 2 z zz 2 y Alternatively if K t has a non-zero value in the element geometric properties data section. the polar second moment of area).Element Formulations ψ xy = − ψ xz = − ∂2 w ∂x∂y ∂2w ∂x∂y Note. for circular cross-sections K t = J .M1]. 7-7) Note. Tz ∈x ψy. The local y and z-axis form a right-hand set with the local x-axis (fig.7.7. For a straight element not parallel to the global X-axis. BSL4 The local xy-plane is defined by all four element nodes which are BXL4 assumed to be coplanar. Greatest accuracy is obtained at the Gauss points.2. Mz Ty .7). The local y-axis forms a right-hand set with the local x and z axes. The local y-axis is perpendicular to the local x-axis and +ve on the side of the element where the fourth node lies.Beam Elements The consistent and lumped mass matrices are evaluated using the procedures defined in (section 2. Force and stress output may be obtained at either the nodes or element Gauss points. Evaluation of stresses/forces The element output obtained at the nodes or Gauss points consists of Fx My .7-6a). the local z-axis is given the unit vector z = i∗ x (i is a unit vector along the global X-axis) (fig. ψ xz axial force moments torques axial strain flexural strain torsional strain The forces and strains are output in the local Cartesian system which is defined by BSL3 For a curved element the local xy-plane is defined by the three element nodes.7.7.7-6b).2.2. The torques are +ve for anti-clockwise rotations at the first node and clockwise rotations at the third node.2. 39 . For a straight element parallel to the global X-axis.7-6c). Local y is perpendicular to local x and +ve on the convex side of the element. the local z-axis is given by the unit vector z = j∗ x (j is a unit vector along the global Y-axis) (fig. The local y and z-axis form a right-hand set with the local x-axis (fig. ψz ψ xy . Nonlinear dynamics utilising the nonlinear material laws specified in 1. Geometrically nonlinear analysis.2). All continuum based nonlinear material models ignore nonlinear torsional effects.2). 40 .Lagrange strains respectively. BXL4 may be used with the concrete model and continuum based plasticity models (section 4. The nonlinear straindisplacement relationship is defined by ∈x = ∂u 1 ∂u + ∂x 2 ∂x LM OP N Q 2 + 1 ∂v 2 ∂x LM OP N Q 2 ψy = − ψz = − ∂ 2 w ∂u ∂ 2 w ∂w ∂ 2 u ∂w ∂ 2 v − + + ∂x 2 ∂x ∂x 2 ∂x ∂x 2 ∂y ∂x 2 ∂ 2 v ∂u ∂ 2 v ∂ v ∂ 2 u ∂ w ∂ 2 w − + − ∂x 2 ∂x ∂x 2 ∂x ∂x 2 ∂y ∂x 2 ∂ 2 w ∂u ∂ 2 w ∂w ∂ 2 v − − ∂x∂y ∂x ∂x∂y ∂x ∂x 2 ∂2 w ∂u ∂ 2 w ∂v ∂2 v − − ∂x∂y ∂x ∂x∂y ∂x ∂x 2 ψ xz = − ψ xy = − γ yz = − ∂v ∂w ∂x ∂x with reference to the local element x-axis. referred to the undeformed configuration. Notes BSL3. The force and strain output with the geometrically nonlinear analysis will be the 2nd Piola-Kirchhoff stress resultants and Green.2). The loading is conservative. Geometrically and materially nonlinear analysis utilising the nonlinear material laws specified in 1.Element Formulations Nonlinear formulation The beam elements can be employed in Materially nonlinear analysis utilising the elastoplastic constitutive laws [O1] (section 4. Linear eigen-buckling analysis. BSL4 and BXL4 may be used in conjunction with the stress resultant plasticity model (section 4. The geometric nonlinearity utilises a Total Lagrangian formulation which accounts for large displacements but small strains. BSL4 And BXL4 Elements QSL8 elements BSL3 elements Problem Definition Finite Element Mesh Fig.Beam Elements V θY W θZ 5 θZ θX θY 3 U θX V W 5 θX θZ V θY 3 U θX V W θY V θY Y W θZ 1 X Z θZ 4 2 U W 2 U V θX Y W θZ 1 X Z θY 4 U θX θX U θX Initial Variables Final Variables Fig.7.2.7-1 Nodal Freedoms For BSL3.7-2 Stiffened Shell Illustrating Use Of BSL3 Element 41 .7.2. 2.7.7-3 Local Cartesian Axes For Cross-Section Y Element 1 or Quadrilateral 1 2 1 2 3 3 4 Element 2 or Quadrilateral 2 Z Element 3 or Quadrilateral 3 2 1 1 4 3 4 Fig.Element Formulations Y 2 3 1 4 Z Fig.7-4 Cross-Section Of An I-Beam Represented By Superimposing Three BSL4 Elements Or By Defining Three Quadrilaterals 42 .7.2. 7.Beam Elements 3*3 Newton-Cotes 5*5 Newton-Cotes Fig.7-5 Quadrature Rules For Cross-Section Integration 43 .2. 2.7.7-6 Local Cartesian System For BSL3 Element 44 .Element Formulations Y x-y plane y y x y z 1 X x 3 2 z x Z (a) Curved Element Y z x 1 y 2 3 y X z x Z (b) Straight Element Parrallel With Global X-axis y 3 z y x 1 z X 2 Y x Z (c) Arbitrarily Orientated Straight Element Fig. Consistent and lumped mass matrices are available which are evaluated using the procedures defined in section 2.2. straight beam formulated using Timoshenko beam theory so that shear deformations are accounted for. W. (fig. The nodal degrees of freedom for BTS3 are identical to those of the BMS3 element.2.1. In essence.2. End releases may be applied to all the nodal freedoms.7. θ X . two noded.8 3-D Straight Beam (BTS3) Formulation This element is a 3-D. this element is formulated in a very straight forward manner. using linear shape functions and standard degrees of freedom.7.8.7-7 Local Cartesian System For BSL4 And BXL4 Elements 7.2. see section 7. All displacement and rotation variations along the length of the element are linear while all internal forces and moments are constant. The complexities in this formulation arise in the consistent derivation of the geometric tangent stiffness and in the treatment of the rotational degrees of freedom. θ Z at each node.4.Beam Elements Y x y 1 x y 4 3 x z y 2 X Z Fig.7. The nodal degrees of freedom are U. Evaluation of stresses/forces The element output consists of 45 . V. θ Y .8-1). Element Formulations Fx . These strains always relate to a local Cartesian system. Fy . the strains computed at the end of one increment do not depend on the strains computed at the 46 . These axes are consistent with those of the BMS3 element (fig. Linear eigen-buckling analysis. Nonlinear dynamics utilising the nonlinear material laws specified in 1. The local x-axis lies along the element axis in the direction in which the element nodes are specified.Moments in the local Cartesian system Element strains and curvatures are also available but nodal values are not output.2). The formulation is such that engineering strain measures are used in both linear and geometrically nonlinear applications.8-2). ∈ are the local strains and curvatures and D is the modulus matrix given by (terms not shown are zero) LM EA MM D=M MM MMEA e N xx EA xx e z GA sy GAsz G J xx + A xx e 2 z e j xx z E I yy + A xx e 2 z e j EI zz OP PP PP PP PQ Nonlinear formulation This element can be employed in Materially nonlinear analysis utilising the stress resultant plasticity model (section 4. My . Total local strains are computed using the current configuration and local frame. Geometric nonlinearity is accounted for using a co-rotational formulation.2. Geometrically nonlinear analysis. The local xy-plane is defined by the third element node and the element xaxis. The internal forces are computed using: P = D∈ Where P are the local internal forces. In this approach local strains are computed in a local Cartesian frame which is 'fixed' to the element and follows the element as it rotates in 3-D space.7. Fz Mx .Forces in the local Cartesian system . The local y and z-axes form a right-hand set with the local x-axis. Mz . In other words. Geometrically and materially nonlinear analysis utilising the nonlinear material laws specified in 1. X x . The current local gradients at the nodes are computed from 2 θ1 = t 2 T e3 − t 3T e 2 2 θ2 = t1T e 2 − t 2 T e1 2 θ3 = t1T e3 − t 3T e1 2 θ4 = q 2 Te 3 −q 3 Te 2 2 θ5 = q T e 2 − q 1 2 3 Te 1 2 θ6 = q T e3 − q 1 Te 1 Where e i are unit vectors defining the co-rotated base frame and t i . (fig. u is the axial stretch measured in the corotated frame. X y .2. These expressions may be thought of as being a means of computing an 'average' value for a local gradient at a node. t1T e 2 = − t 2 T e1 . 7. This is easily visualised in two dimensions where. γ y . lo is the initial element length.This is not true for the three dimensional case and a stricter derivation of the expressions for local gradients would involve the polar decomposition theorem. and θi . Solutions obtained using this element will not be load step size dependant. i=1. i=1.6 are the local gradients at the nodes or 'curvature producing' rotations relative to the co-rotated frame.3 are the i cartesian sets at nodes 1 and 2 respectively defining the orientation of the beam crosssection. The approach then taken is to decompose 47 .8-3). q . X z o t T where: ∈x = u / lo i e γ = − eθ + θ i 2 X = e θ − θ i 10 X = e θ − θ i 10 X = e θ − θ i 10 z 3 6 x y 4 1 2 5 z 3 6 γ y = − θ2 + θ5 2 Where. The local strains for the element are given by ∈ = ∈x .Beam Elements end of a previous increment. for example. γ z . i=1. (fig. The following expressions are used for defining e 2 and e3 . The local frame e i . The variation of these expressions is used in the virtual work equation to relate variations in local strains to variations in global nodal displacements.25 degrees for a local gradient of 15 degrees and 1. This lack of orthogonality has been shown to be 0.3 represent the 'average' of the nodal cartesian sets. These expressions have been used with a view to obtaining less costly derivatives in a consistent derivation of the tangent stiffness.Element Formulations the total rotation into a rigid body component and a local gradient. By defining the vectors x 21 = x 2 − x1 d 21 = d 2 − d1 and then by rearranging Pythagoras's theorem this may be expressed as u= 2 1 x 21 + d 21 ln + lo 2 R S T U V W T d 21 The vector x defines the nodes in the initial configuration while d 21 is the 'net' translational displacement vector. The expressions described above are the result of applying these principles.2. In three dimensions defining the local frame is more difficult. i=1. The variation of this expression reduces to 48 .8-4). The approach described by Crisfield [C7] has been used for this purpose. t i and q . The axial stretch may be taken as u = ln − lo where ln is taken as the current element length (or chord length).3 is easily established for a two dimensional problem. 7. e2 = r 2 − e3 = r 3 − r 2 T e1 2 r 3T e1 2 ne ne 1 + r1 + r1 s 1 s i The local frame is established at the centre of the element and the vectors r i . As these values actually represent the 'curvature producing' rotations in a single element the deformation would need to be very severe to reach these values. These expressions are approximations to the exact expressions for defining the 'smallest' rotation between vectors r1 and e1 .9 degrees for a gradient of 30 degrees [C7]. These axes. Notes This geometrically nonlinear formulation is consistently formulated and displays a quadratic rate of convergence in the limit. Using this expression in the virtual work equation allows global internal forces at the nodes to be expressed in terms of local internal forces as P = BT P The out of balance force vector is then given by ψ = BT P − R where R is the applied nodal loading. As explained in section 3. This element incorporates rotational degrees of freedom. Assuming conservative loading this gives d ψ = BT dP + dBT P d ψ = BT DBda + dBT P The first term on the right hand side of this equation may be recognised as the material or standard linear stiffness matrix. which are used in the computation of internal forces and the stiffness matrix. The total strains are computed from the current configuration and local frame only. are updated correctly 49 . Therefore. A set of Cartesian axes are established at each node to define the orientation of the beam cross section. To overcome this problem the rotation variables are never added to establish the current orientation of the element.Beam Elements δu = e1T δ d 21 Differentiation of the above equations relating to axial stretch and local gradients allows virtual variations of local strains to be related to virtual variations in global nodal displacements via a strain displacement matrix B δ∈= B δ a where a are the global nodal displacements. results obtained using this element are not load step size dependant. The variation of this equation gives the tangent stiffness matrix. The second term gives rise to the geometric stiffness.5. A consequence of this consistency is the ability of the element to cope with larger load increments. large rotations in three dimensions are non-vectorial in nature and therefore may not be summed as vectors. In view of the non-vectorial nature of these rotations it should be noted that the nodal rotation output represents approximate values which should be treated with caution. This arises as additional terms are added to the stiffness matrix to account for the variation in the load direction between iterations. In the deformed configuration the node is no longer completely shared and from (fig. with the gap in (fig. Full details of the derivation of these additional terms are given in [C6].8-5) drawn for illustrative purposes only.Element Formulations using the iterative increments in nodal rotations. m define the displacement vector and rotation matrix of the master 50 . the translational displacements and internal force output will be correct for problems involving arbitrary large nodal rotations. Prismatic (sliding).5.7. following conventional beam theory assumptions. spherical and cylindrical joints can be modelled by releasing the appropriate degrees of freedom at a node. Form an incremental Euler parameter from the iterative rotation increment. A nonsymmetric stiffness matrix will result if a follower load is specified. robots and rotating machinery.2. the origins of the vectors d m and d coincide. Examples include deployable space structures. and d and Q define the displacement vector and rotation matrix of the disconnected (at least partially) slave node. These freedoms relate to the local beam axes and a master-slave procedure has been adopted to model the release [J2]. one of which is not fully connected to the others. A more detailed derivation of this element formulation may be found in [C6] and [C8]. Consider a node which is initially shared by a number of elements. However. It should be noted that. section 3.8-5) the following relationships can be established: d = dm + ρ Q = Q* Q m where d m and Q node.7.2. Update the Euler parameter by manipulating the previous and incremental parameter using quarternion algebra. At present. The procedure for this operation is outlined as Extract the Euler parameters from the initial nodal Cartesian set. this facility is restricted to static analyses. revolute (hinges). End releases Many structures which are modelled with three-dimensional beam elements require joints at the nodes which follow the axes of the rotating system. Form the updated Cartesian set from the updated Euler parameter. d and Q .2. In a similar manner. equal to the vector of released displacements (fig. q . Q* and it is transformed to the master triad. between the master and slave variables is. d m and Q . ρ (with local components). some of the components of the displacement vectors. the master variables.2.q Q m3 and q . Depending on the type of joint. d m and d . when transformed into coordinates defined by the master triad. are generally not entirely independent from the slave variables. q .8-6) and (fig. if the rotational pseudovector β* . q (fig. slave and released freedoms.7. For translational joints. has zero components in non-released directions. In a geometrically nonlinear analysis these axes rotate together with the structure. and/or parameters of the rotation matrices Q and Q .7.7. m m When modelling different types of joints.2.8-7) illustrate a prismatic (sliding) and revolute (hinge) release.8-5): 1 2 3 m1 m and Q consist of orthonormal base vectors m = q . q m2 .8-6): s = QT ρ m m where the vector of released displacements.Beam Elements The columns of the rotation matrices Q q m1 .7. 51 .q m2 . Full details of this derivation can be found in [J2] while (fig. Different types of joints are defined by releasing displacements and/or rotations around chosen axes.8-7): ψ = Q T β* m Using these equations a relationship can be established between the variations of the master. the 'difference vector'. This relationship can then be used to derive a modified stiffness matrix and internal force vector which accounts for any released freedoms. q 1 2 3 The rotation matrix Q* is the matrix that defines the rotation of the master triad Q .7. s . the rotational pseudovector of released rotations is obtained (fig. on to the slave triad Q.2. is extracted from the 'rotation difference matrix'.2. q m3 Q= q . can be the same. 7.Element Formulations Y θY V θY 1 W θZ U θX X Z W 2 θZ θX Fig.8-2 Local Cartesian System For BTS3 Element 52 .2.2.7.8-1 Nodal Freedoms For BTS3 Element Y y x 2 3 z 1 X Z Fig. 8-4 Axial Stretch BTS3 Element 53 .7.8-3 Local Gradients BTS3 Element Final Configuration ln d2 lo d1 Y x1 X Z x21 x2 Initial Configuration Fig.Beam Elements q2 θ6 θ4 q1 q3 e2 e1 e3 t1 t2 Y θ3 θ1 t3 θ2 θ5 X Z Fig.2.7.2. 8-6 Prismatic (Sliding) Release 54 .Element Formulations e3 e2 e1 q3 d dm qm3 ρ qm2 q2 q1 qm1 Fig.8-5 General Displacements At A Node With Released Freedoms e3 e2 e1 qm3 ρ qm2 qm1 S Fig.2.2.7.7. Two-Dimensional Continuum Elements e3 e2 qm3 q3 e1 qm3= q1 q2 qm2 ϕ Fig.B2]. TPM3.8-7 Revolute (Hinged) Release 7. The consistent and lumped mass matrices are evaluated using the procedures defined in (section 2. TPM6.3.7.e. Fig. ηgUi i =1 n n geometry X= ∑ Ni bξ. TPK6) The plane stress elements are formulated by assuming that the variation of out of plane direct stress and shear stresses is zero. The nodal degrees of freedom are U and V. QPM8. b g 55 .5 on space membrane elements A complete description of their formulation is given in [H1. η is the element shape function for node i and n is the number of nodes. ηgXi i =1 where N i ξ. For 3-D plane membrane elements see section 7. i. displacement U= ∑ N i bξ.e.3 Two-Dimensional Continuum Elements 7. All the isoparametric elements described in this section must be defined using only X and Y coordinates. Plane stress (QPM4. QPK8. i.7.7).3.1 Standard Isoparametric Elements Isoparametric finite elements utilise the same shape functions to interpolate both the displacements and geometry.1-1 shows the nodal configurations available within LUSAS.2. 7. To obtain a valid material υ xy < E x / E y ∈z = − ∈z = − d i 1/ 2 υ σx + σy E d i for isotropic materials for orthotropic materials υyz υ xz σx − σy Ex Ey The thermal strain is defined by 56 .g.Element Formulations σ z = 0. The infinitesimal strain-displacement relationship is defined as ∈X = ∈Y = ∂U ∂X ∂V ∂Y ∂U ∂V + ∂Y ∂X γ XY = The isotropic and orthotropic elastic modulus matrices are Isotropic Orthotropic LM1 Mυ D= e1 − υ j MM0 N LM 1 / E D = M− υ / E MN 0 E 2 υ 1 0 OP 0 P a1 − υf PP 2 Q 0 x − υ xy / E x 1 / Ey 0 0 0 1 / G xy xy x OP PP Q −1 where υyx has been set to υxy E y / E x to maintain symmetry. thin plates subject to in-plane loading (fig.1-2).3. Note. σ yz = 0 The plane stress elements are suitable for analysing structures which are thin in the out of plane direction. The thickness of the material is defined at each node and may vary over the element. e. σ xz = 0. Note. Two-Dimensional Continuum Elements isotropic LM(1 − υ) E M υ D= (1 + υ)(1 − 2 υ) M MN 0 (∈0 )t = ΔT α x , α y , α xy T υ (1 − υ) 0 OP 0 P (1 − 2 υ) P P 2 Q 0 Orthotropic Plane strain (QPN4, QPN8, TPN3, TPN6, QNK8, TNK6) The plane strain elements are formulated by assuming that the variation of out of plane direct strain and shear strains is zero, i.e. ∈Z = 0, ∈YZ = 0, ∈XZ = 0 The plane strain elements are suitable for analysing structures which are thick in the out of plane direction, e.g. dams or thick cylinders (fig.7.3.1-3). The infinitesimal strain-displacement relationship is defined as ∈X = ∈Y = ∂U ∂X ∂V ∂Y ∂U ∂V + ∂Y ∂X γ XY = The isotropic and orthotropic elastic modulus matrices are Isotropic Orthotropic LM(1 − υ) E M υ D= (1 + υ)(1 − 2 υ) M MN 0 LM E − υ E E E MM −υ E − υ υ E D=M MM E E MM 0 N z 2 xz x x z xy z yz xz x x z υ (1 − υ) 0 OP 0 P (1 − 2 υ) P P 2 Q 0 − υ xy E z − υ xz υ yz E y E z − υ2 E y yz EyEz 0 EyEz 0 0 1 G xy OP PP PP PP PQ −1 where for symmetry E y υ xy E z + υ yz υ xz E x = E x υ xy E z + υ xz υ yz E y d i d i Note. To obtain a valid material 57 Element Formulations υ xy < E x / E y σz X d = υbσ i 1/ 2 υ xz < E x / E z b g 1/ 2 υ yz < E y / E z d i 1/ 2 + σY g for isotropic materials for orthotropic materials σ z = υ xz Ez E + υ yz z EX Ey The thermal strain is defined by Isotropic Orthotropic d∈ i = (1 + υ)ΔT α, α, 0 LM E υ α + α , E d∈ i = ΔT MM E υ α MN E 0 t T z x xz z x 0 t z y yz z + α y , α xy OP PP PQ T Axisymmetric (QAX4, QAX8, TAX3, TAX6, QXK8, TXK6) The axisymmetric elements are formulated by assuming that the variation of out of plane shear stresses is negligible, i.e. σ XZ = 0, σ YZ = 0 and the out of plane direct strain is defined as ∈Z = U R where R is the distance from the axis of symmetry. The axisymmetric elements are suitable for analysing solid structures which exhibit geometric symmetry about a given axis, e.g. thick cylinders or circular plates (fig.7.3.1-4). The elements are defined in the XY-plane and symmetry can be specified about either the X or Y axes. The infinitesimal strain-displacement relationship is defined as: ∈X = ∈Y = ∂U ∂X ∂V ∂Y ∂U ∂V + ∂Y ∂X γ XY = ∈Z = U R symmetry about the Y axis 58 Two-Dimensional Continuum Elements or ∈Z = V R symmetry about the X axis The isotropic and orthotropic linear elastic modulus matrices are defined as Isotropic Orthotropic LM(1 − υ) E M υ D= (1 − υ)(1 − 2 υ) M 0 MN υ LM 1 / E −υ / E −υ / E 1/ E D=M MM 0 0 −υ / E −υ / E NM x yx xy x y xz x yz υ (1 − υ) 0 υ 0 0 (1 − 2 υ) / 2 0 OP P 0 P P (1 − υ)Q υ υ y 0 0 1 / G xy 0 − υzx / E z − υzy / E z 0 1 / Ez y OP PP PP Q −1 where υyx , υzx and υzy are defined by υ yx = υ xy E y / E x υzx = υ xz E z / E x υzy = υ yz E z / E y to maintain symmetry. Note. To obtain a valid material υ xy < E x / E y d i 1/2 υ xz < E x / E z b g 1/2 υ yz < E y / E zy d i 1/ 2 The thermal strain vector is defined as Isotropic Orthotropic d∈ i = ΔT α, α, 0, α 0 t T d∈ i = ΔT α , α , α 0 t x y xy , αz T Evaluation of stresses The element output obtained at the element nodes and Gauss points consists of Stress Output σ X , σ Y , σ XY , σ Z σ max , σ min β σS σV the direct and shear stresses the maximum and minimum principal stresses the angle between the maximum principal stress and the positive X-axis the maximum shear stress von Mises equivalent stress the direct and shear strains Strain Output ∈X , ∈Y , γ XY , ∈z 59 Element Formulations ∈max , ∈min β ∈S ∈V the maximum and minimum principal strains the angle between the maximum principal strain and the positive X-axis the maximum shear strain von Mises equivalent strain Stress resultant output, which accounts for the thickness of the element, is available as an alternative to stress output for the plane stress elements, i.e Stress Resultant Output N X , N Y , N XY , N z the direct and shear stress resultants/unit length N max , N min the maximum and minimum principal stress resultants/unit length β the angle between the maximum principal stress resultant and the positive X-axis N S the maximum shear stress resultant/unit length N V von Mises equivalent stress resultant/unit length The sign convention for stress, stress resultants and strain output is shown in fig.7.3.1-6. The Gauss point stresses are usually more accurate than the nodal values. The nodal values of stress and strain are obtained using the extrapolation procedures detailed in section 6.1. Nonlinear formulation The 2-D isoparametric elements can be employed in:(Materially nonlinear analysis, utilising the elasto-plastic constitutive laws [O1] (section 4.2) and the concrete model (section.4.3) Geometrically nonlinear analysis. Geometrically and materially nonlinear analysis utilising the nonlinear material laws specified in 1. Nonlinear dynamics utilising the nonlinear material laws specified in 1. Linear eigen-buckling analysis. Notes The plane stress elements can be used with the nonlinear concrete model (section 4.3). The plane stress and plane strain elements may be used with the nonlinear interface model (section 4.2). The geometric nonlinearity may utilize: A Total Lagrangian formulation which accounts for large displacements but small strains. The nonlinear strain-displacement relationship is defined by 60 Two-Dimensional Continuum Elements Plane stress LM OP N Q ∂V 1 L ∂U O ∈ = + ∂X 2 M ∂Y P N Q ∈X = ∂U 1 ∂U + ∂X 2 ∂X Y 2 2 LM OP N Q 1 L ∂V O + M P 2 N ∂Y Q + 1 ∂V 2 ∂X 2 2 γ XY = ∂U ∂V ∂U ∂U ∂V ∂V + + + ∂Y ∂X ∂X ∂Y ∂X ∂Y 2 Plane Strain LM OP N Q ∂V 1 L ∂U O ∈ = + ∂X 2 M ∂Y P N Q ∈X = ∂U 1 ∂U + ∂X 2 ∂X Y 2 LM OP N Q 1 L ∂V O + M P 2 N ∂Y Q + 1 ∂V 2 ∂X 2 2 γ XY = ∂U ∂V ∂U ∂U ∂V ∂V + + + ∂Y ∂X ∂X ∂Y ∂X ∂Y 2 Axisymmetric ∈X = LM OP N Q ∂V 1 L ∂U O ∈ = + ∂X 2 M ∂Y P N Q ∂U 1 ∂U + ∂X 2 ∂X Y 2 LM OP N Q 1 L ∂V O + M P 2 N ∂Y Q + 1 ∂V 2 ∂X 2 2 γ XY = ∂U ∂V ∂U ∂U ∂V ∂V + + + ∂Y ∂X ∂X ∂Y ∂X ∂Y 2 or LM OP N Q V 1 LVO ∈ = + M P R 2 NRQ ∈Z = U 1 U + R 2 R Z symmetry about the Y axis symmetry about the X axis 2 The output is now in terms of the 2nd Piola-Kirchhoff stresses and Green-Lagrange strains referred to the undeformed configuration. The loading is conservative. An Updated Lagrangian formulation, which takes account of large displacements and moderately large strains provided that the strain increments are small. The output is now in terms of the true Cauchy stresses and the strains approximate to logarithmic strains. The loading approximates to being non-conservative. 61 The output is in terms of true Cauchy stresses and the strains approximate to logarithmic strains. 3 5 6 4 1 2 3 node triangle 1 2 3 6 node triangle 3 4 7 6 5 4 8 1 4 node quadrilateral 2 1 2 8 node quadrilateral 3 Fig. The loading is non-conservative.1-1 Nodal Configuration For Standard 2-D Isoparametric 62 .7.3. which takes account of large displacements and large strains.Element Formulations An Eulerian formulation. 7.Two-Dimensional Continuum Elements Problem Definition (a) Plate subject to Inplane Loading Finite Element Mesh Problem Definition Finite Element Mesh (b) Cantilever subject to a Point Loading Fig.3.1-2 Examples Illustrating The Use Of Plane Stress Elements 63 . 1-3 Examples Illustrating The Use Of Plane Strain 64 .3.7.Element Formulations Problem Definition QPN8 elements TPN6 elements Finite Element Mesh (a) Embankment Dam QPN4 elements Problem Definition (b) Thick Cylinder Finite Element Mesh Fig. Two-Dimensional Continuum Elements QAX4 elements r r Problem Definition (a) Thick Cylinder Finite Element Mesh QAX8 elements Problem Definition (b) Circular Plate Finite Element Mesh Fig.7.1-4 Examples Illustrating The Use Of Axisymmetric Solid Elements 65 .3. strains and displacements are represented by three independent functions in three separate vector spaces. particularly if bending predominates. QPN4M.Element Formulations Fig. σ Y +ve tension σ XY +ve into XY quadrant Fig.QPN4 and QAX4.3.3. The formulation provides for the following three conditions to be satisfied 66 . QAX4M) The lower order enhanced strain elements exhibit improved accuracy in coarse meshes when compared with their parent elements QPM4.3.2 Enhanced Strain Elements (QPM4M. The elements are based on a three-field mixed formulation [S8] in which stresses. In addition.1-6 Sign Convention For Stress/Strain Output 7.7.1-5 Deformed Mesh Illustrating Formation Of Spurious Mechanisms Y σY σ XY σX σ XY σY X σX σ X. these elements do not suffer from 'locking' in the nearly incompressible limit. These internal degrees of freedom are eliminated at the element level before assembly of the stiffness matrix for the structure.7. The formulation is based on the inclusion of an assumed 'enhanced' strain field which is related to internal degrees of freedom. compatibility and constitutive relationship may be expressed as z z z Ω δ∈c σ dΩ − R T δ a = 0 δσ ∈e dΩ = 0 δ∈ − σ + T T T Ω Ω LM N ∂W ∂∈ OP dΩ = 0 Q where R is the applied loading. these conditions also allow the stress field to be eliminated from the formulation. The weak form of the three field variational equations for equilibrium. By enforcing the so called L 2 orthogonality condition between stress and enhanced strain. Formulation The formulation requires that the total strain is expressed as the sum of a 'compatible' strain and an 'enhanced' strain ∈= ∈c + ∈e The compatible strain is directly related to the displacements of the element nodes in the standard manner. Capability of the element to model a constant state of stress after enforcing the orthogonality condition. W is the strain energy density.e. The compatible and enhanced strains are computed from ∈c = Ba ∈e = Gα e where G operates on the assumed strain parameters α e to provide the enhanced strains. The enhanced strain is related to internal degrees of freedom which are eliminated using static condensation at the element level. i. This allows the stress field to be eliminated from the formulation. 67 . a are nodal displacements and σ is the stress vector.Two-Dimensional Continuum Elements Independence of the enhanced and standard strain interpolation functions. L2 orthogonality of the stress and enhanced strains. The enhanced strains are therefore discontinuous between elements. In addition to ensuring that the element passes the patch test. terms involving σ T ∈e will disappear. requirement for passing the patch test. D is the modulus matrix at loadstep k. m el is the number element enhanced strain modes.Element Formulations Substitution of these expressions into the two remaining field variational equations yields z Ω δ∈c D Bda + Gd α e − R Tδa dΩ + T { } z Ω δ∈e D Bda + Gd α e T { } dΩ = 0 The following matrices are defined for use in discretising this equation Ka k f = H Γ k z af z af z = k Ω BT D B dΩ G T D G dΩ G T D B dΩ ( n el * n el matrix) ( m el * m el matrix) ( m el * n el matrix) Ω = Ω n el is the dimension of the element displacement field. Static condensation of this system of equations eliminates the equations included to enforce the orthogonality condition. The element stiffness and internal forces used to assemble the equations for the structure then become K a k fΔaa k +1f = ψ a k f a k f = Ka k f − Γ a k fT Ha k f −1 Γ a k f 68 where K . no iterations are necessary as h will be 0 and P will not be considered. The internal force vectors are given by Pa k f = ha k f = z z Ω Ba k fT σ a k f dΩ Ga k fT σ a k f dΩ Ω where h (k) is the internal force vector relating to the incompatible modes which is subsequently eliminated at the element level. Using standard finite element techniques for assembling the system of equations gives LMKa f MN Γa f k k Γ a k fT OR Δaak +1f U RR − Pakf U | | | | a k f PSΔαa k +1f V = S 0 − hak f V H P| | | | W W T T Q This nonlinear system of equations is solved using a Newton-Raphson iteration scheme. However. for the linear case. Standard transformations are applied and full details of this procedure are given in [S8]. U Rx a U Rx h U | = | | + ξ| | = g V S V S V N. 2 2 N2 η = a f 1 e1 − η j 2 2 and λ i represent the incompatible modes λ1 = u1. the enhanced strain parameters are updated as α a k +1f = α a k f − Ha k f Γ a k fΔ a a k f − ha k f The actual implementation of this formulation requires the orthogonality condition to be related to the isoparametric space. v1 l q. v2 l q T The covariant base vectors associated with the isoparametric space are g = ξ g ξ Rx | Sy | T Rx | =S |y T T N.Two-Dimensional Continuum Elements ψakf = R − P P akf a k f = P a k f − Γ a k f T H a k f −1 h a k f −1 In nonlinear analyses. T λ 2 = u 2 . Transformations are therefore required to assemble matrices and vectors that relate to covariant strains and contravariant stresses. Enhanced strain interpolation . QPN4M) The incompatible displacement field is given by U = N1 ξ λ1 + N 2 η λ 2 bg af where N1 ξ = b g 1 e1 − ξ j. It is postulated that the covariant enhanced strain field is given by ∈e = Eα e where E is the equivalent of G in the isoparametric space. | |y a | |y h | W T W T W ξ ξ T T 1 1 T T 0 1 + ηg η η T T 2 2 T T 0 2 + ξg where N = a 0 + a1ξ + a 2 η + hξη 69 . T T T U = Rx a U + ηRx hU = g | | | | | V S V Sy h V N. | |y a | T W W T W | | N.plane elements (QPM4M. i g =S | |du. The remaining terms satisfy the equilibrium equations. j | T T T ξ ξ η η T η ξ ξ T g η U | | V | | | W T 1 T 2 T 1 T 2 T 1 T 2 0 1 0 2 0 2 0 1 ξ 0 0 0 ξη = 0 η 0 0 0 0 0 ξ η ξ2 LM MM N R− λ g U |− λ g | | | 0 O| | P |− λ g | = E α ξηP S V −λ g | η P| Q| −λ g | | | | −λ g | T W i 2 i e The stress field for the element is derived from the linear uncoupled stress field [P2] σ = σy = xy Rσ | S |τ T x U | V | W LM1 ξ η MM N 1 ξ η OP PPβ 1 ξ ηQ The introduction of four internal degrees of freedom allows four of the nine stress parameters (β) to vanish. By basing the formulation on natural coordinates the element is less sensitive when distorted and possesses no zero-energy deformation modes. j g | | d u. The final contravariant stress field using five β parameters is defined as 70 . i g + eu. Full details of the elimination of the four stress parameters is described in [P2] for a hybrid element.Element Formulations and a0 = a1 = a2 = h= 1 T 1 1 1 1 4 1 −1 4 1 1 −1 1 1 T 1 −1 −1 4 1 −1 T 1 4 1 −1 T T T x = x1 x 2 x3 x 4 y = y1 y 2 y3 y 4 The initial enhanced strain field in isoparametric space is then given by ∈i R e u. and hence L 2 orthogonality. This is a requirement for passing the patch test [S8] and is implied in the sense that: zz 1 1 −1 − 1 E dξdη ≡ 0 Enhanced strain interpolation . The initial matrix is given by Ei LMξ 0 =M MM0 N0 1 0 η 0 0 0 0 ξ 0 0 0 0 0 η 0 0 ξη OP PP PQ For the axisymmetric case.η) will be included in the integrand for enforcing orthogonality zz 1 −1 − 1 r ∈T ∑ d ξd η ≡ 0 71 . the condition is satisfied if α 5 = −α 6 Forcing this equality. The final interpolation functions E also allow condition (III) to be satisfied. ∈ > L 2 zz 1 1 −1 − 1 r ∈T ∑ d ξd η ≡ 0 This condition is violated if the six initial enhanced strain parameters (α) are used. a factor r(ξ.axisymmetric element (QAX4M) The procedure for establishing the enhanced strain interpolation matrix for the axisymmetric element is similar to that used for the plane elements. gives the final enhanced strain interpolation matrix as LMξ E = M0 MN0 0 0 0 η 0 0 0 ξ η OP − ξη P ξ −η P Q ξη 2 2 This matrix is used in linear analyses but for nonlinear applications four enhanced strain parameters are used with the final column of E deleted [S8].Two-Dimensional Continuum Elements R∑ | ∑ = S∑ |∑ | T ξ η ξη U L1 | = M0 V M | NM0 | W 0 0 η 0 1 0 0 ξ β 0 1 0 0 OP PP Q To satisfy the L 2 orthogonality condition < ∑ . However. η= . The geometrically nonlinear performance of these elements is much improved in comparison with the standard elements. z3 .Element Formulations where ∑ = ∑ ξξ ∑ ηη ∑ ξη ∑ θθ T .1. Nonlinear formulation The comments made in section 7. z 2 .5. a 2 and h vectors have been defined for the plane elements. a1 . ξη = Ta Ta 3r 0 3r 0 9 rT a0 a 0 . ∈ = ∈ ∈ 2 ∈ξη ∈ ηη θθ ξξ r = rT N z = zT N r = r1. Simo and Rifai [S8] have derived interpolation functions which account for the factor r and satisfy this condition E = Ei − zz 1 1 1 −1 − 1 r dξdη zz 1 1 −1 −1 E i r dξdη LMξ − ξ 0 =M MM 0 MN 0 where ξ= 0 0 0 0 η− η 0 0 0 0 ξ − ξ η− η 0 0 0 0 ξη − ξη OP PP PP Q 1 r T a1 1 rT a2 1 rT h . z 4 Inclusion of the factor r( ∈. r4 z = z1.3. r2 . The nonlinear formulation for the enhanced strain elements involves enforcing orthogonality between assumed Green-Lagrange strains and 2nd Piola-Kirchhoff stresses. Evaluation of stresses The evaluation of stresses is identical to that described in section 7.6 regarding the nonlinear capability of the standard elements are also applicable to these elements. 72 .η) means that the orthogonality condition is violated using this interpolation matrix.3. r3 .1. η = 1 − η2 b g and a i are nodeless degrees of freedom which are condensed out before element i assembly. 4-noded. plane membrane element. ηgU i i =1 n is replaced with U= ∑ b g i =1 n N i ξ. i. i. The infinitesimal strain-displacement relationship is the same as QPM4. The nodal configuration and non-conforming shape functions are shown in fig.e.Two-Dimensional Continuum Elements 7.e.e. η U i + ∑ Pi bξ.3 Incompatible Plane Membrane Element (PMI4) Formulation This element is a high performance.3. The element passes the patch test (ensuring convergence as the mesh is refined) and the displacement field is approximately an order higher than the QPM4 element (i.3. non-conforming.3-1. quadratic displacement accuracy). ηga i i =1 2 where P1 ξ. It is formed by adding two non-conforming modes to the standard isoparametric formulation presented for QPM4.7. ∈X = ∈Y = ∂U ∂X ∂V ∂Y ∂U ∂V + ∂Y ∂X γ XY = The isotropic and orthotropic elastic modulus matrices are defined as Isotropic LM1 Mυ D= e1 − υ j MM0 N E 2 υ 1 0 OP 0 P a1 − υf PP 2 Q 0 73 . η = 1 − ξ 2 b g and P2 ξ. U= ∑ N i bξ. 7). ∈Y .Element Formulations Orthotropic L 1 D = M− υ / E MM N 0 xy x − υ yx / E y 1 / Ey 0 OP P 1 / Gxy P Q 0 0 −1 where υyx is set to υxy E x / E y to maintain symmetry. Note.7. Only a lumped mass matrix is evaluated using the procedure defined in (section 2. Strain Output ∈X . α T xy Full details of the formulation are presented in [T2. α. Evaluation of stresses The element output obtained at the element nodes consists of Stress Resultant Output N x . N min the maximum and minimum principal stress resultants/unit length β the angle between the maximum principal stress resultant and the positive X-axis NS the maximum shear stress resultant/unit length NV Von Mises equivalent stress resultant/unit length. N y . 0 0 t 0 t x y T d∈ i = ΔT α .W2]. γ XY ∈max . 74 .3. α . N xy the direct and shear stress resultants/unit length N max .3-4. ∈min β ∈S ∈V the direct and shear strains the maximum and minimum principal strains the angle between the maximum principal strain and the positive X-axis the maximum shear strain Von Mises equivalent strain The sign convention for stress resultant and strain output is shown in fig. For a valid material υxy < E x / E y The thermal strain is defined by Isotropic Orthotropic d i 1/ 2 d∈ i = ΔT α. The stress resultants are evaluated directly at the nodes. 3-2 Examples Illustrating The Use Of PMI4 Elements 75 . The element cannot be used for linear buckling analyses. Y.3.7.Two-Dimensional Continuum Elements Nonlinear formulation The element has no nonlinear capability. but may be utilised in a nonlinear environment.3-1 Nodal Configuration And Non-Conforming Shape Functions For The PMI4 Element Problem Definition (a) Plate subject to Inplane Loading Finite Element Mesh Problem Definition (b) Cantilever Plate subject to Point Loading Finite Element Mesh Fig.7.3.U (a) Nodal Configuration P1 = 1-ξ 2 P2 = 1-η 2 (b) Non-conforming shape functions Fig.V 3 4 2 1 X. 76 .3.7.3. σ Y +ve tension σ XY +ve into XY quadrant Fig.7.3-4 Sign Convention For Stress/Strain Output 7.3. A number of further advantages may also be obtained in explicit dynamic analyses The use of higher order shape functions creates difficulties at the contact interface in the form of uncontrolled overlap.3-3 LOCAL CARTESIAN SYSTEM FOR THE PMI4 ELEMENT Y σY σ XY σX σ XY σY X σX σ X.Element Formulations 3 Y 4 2 y x 1 X FIG. It has been shown that higher order continuum elements require a time step reduced from that of linear elements because of the greater mass associated with the interior nodes.4 2D Explicit Dynamics Elements Explicit time integration schemes have used simple linear elements rather than those of a higher order by virtue of their computational efficiency. 7. ηgXi i =1 where N i ξ. The combination of mass lumping with linear elements. i.e. η is the element shape function for node and n is the number of nodes. Fig. displacement U= ∑ N i bξ.3.g. All the explicit dynamics elements described in this section must be defined using only X and Y coordinates. e.3. thin plates subject to in-plane loading (fig. when applied in conjunction with the central difference operator. σ XZ = 0. They are for use only with the explicit central difference time integration scheme. A rate relationship is used to define the strain-displacement characteristics as t ∈x = & t ∈y = & & ∂t U t ∂X & ∂t V t ∂Y 77 . The linear explicit dynamics elements have been implemented to take advantage of these benefits. Note that the thickness of the material is defined at each node and may vary over the element. The explicit dynamics elements are based upon the isoparametric approach in which the same shape functions are used to interpolate both the displacements and the geometry. ηgU i i =1 n n geometry X= ∑ Ni bξ.e.Two-Dimensional Continuum Elements The mass lumping formulations for higher order elements are currently impractical for modelling shock wave propagation since the resulting numerical noise pollutes or destroys the solution.1-1 shows the nodal configurations available within LUSAS. Plane stress (QPM4E. i.7. σ YZ = 0 b g The plane stress elements are suitable for analysing structures which are thin in the out of plane direction. TPM3E) The plane stress elements are formulated by assuming that the variation of the out of plane direct stress and shear stresses is negligible. σ Z = 0.4-2). increases accuracy in solutions by virtue of their respective compensatory spectral errors. The nodal degrees of freedom are U and V. α T xy Plane strain (QPN4E.Element Formulations tγ & = & & ∂t U ∂t V + t t ∂Y ∂X XY t ∈z = − υ & & & LM ∂ U + ∂ V OP MN ∂ X ∂ Y PQ t t t t The isotropic and orthotropic elastic modulus matrices are Isotropic Orthotropic LM1 Mυ D= e1 − υ j MM0 N L 1 D = M− υ / E MM N 0 E 2 xy x υ 1 0 OP P a f PP Q 0 0 1− υ 2 0 − υ yx / E y 1 / Ey OP P 1 / Gxy P Q 0 0 −1 where υyx is set to υxy E y / E x to maintain symmetry. Note.3.4-3). ∈Z = 0. ∈YZ = 0. e. dams or thick cylinders (fig. i. α.7. To obtain a valid material υxy < E x / E y The initial thermal strain is defined by Isotropic Orthotropic d i 1/ 2 d∈ i = ΔT α. ∈XZ = 0 The plane strain elements are suitable for analysing structures which are thick in the out of plane direction. α . A rate relationship is used to define the strain-displacement characteristics as t& ∈X = & ∂t U t ∂X & ∂t V ∂tY t& ∈Y = 78 . 0 0 t 0 t x y T d∈ i = ΔT α . TPN3E) The plane strain elements are formulated by assuming that the variation of the out of plane direct strain and shear strains is negligible.e.g. α .e. TAX3E) The axisymmetric elements are formulated by assuming that the variation of out of plane shear stresses is negligible. 0 d∈ i = a1 + υfΔT α . α . σ XZ = 0.Two-Dimensional Continuum Elements t & & ∂t U ∂t V & γ XY = t + t ∂Y ∂X ∈Y = 0 t& The isotropic and orthotropic elastic modulus matrices are Isotropic LM(1 − υ) υ O 0 P E M P D= a1 − υf(1 − 2υ) MM υ (1 − υ) a1 −02υf PP 0 P 2 Q NM 0 LM E − υ E E E MM −υ E − υ υ E Orthotropic D=M MM E E MM 0 N z 2 xz x x z xy z yz xz x x z − υ xy E z − υ xz υyz E y E z − υ2 E y yz EyEz 0 EyEz 0 0 1 G xy OP PP PP PP PQ −1 where for symmetry E y υ xy E z + υ yz υ xz E x = E x υ xy E z + υ xz υ yz E y d i d i T x y The initial thermal strain is defined by Isotropic Orthotropic d∈ i = a1 + υfΔT α . 79 . σ YZ = 0 and the out of plane direct strain rate is defined as & ∈Z = & U R where R is the distance from the axis of symmetry. i. α 0 t 0 t x y T xy Axisymmetric (QAX4E. Element Formulations The axisymmetric elements are suitable for analysing solid structures which exhibit geometric symmetry about a given axis. For large strain axisymmetric analyses.e. e. These difficulties may be overcome by formulating the elements with the Petrov-Galerkin method [G2]. where the weighting functions are generally the element shape functions. The use of this particular formulation produces a time dependent mass matrix and as such must be computed each time. thick cylinders or circular plates (fig.4-4). in which the governing differential equation is utilised directly to form a weighted residual statement. the weighting functions are taken to be the product of the element shape functions and the inverse of the radius.g. Standard axisymmetric isoparametric elements are formulated with the Galerkin weighted residual method. This method is also a weighted residual method.7. however. The elements are defined in the XY-plane and symmetry can be specified about either the X or Y axes. i. the use of elements based on the Galerkin method leads to computational difficulties near the axis of symmetry. A rate relationship is used to define the strain-displacement characteristics as t& ∈X = ∈Y = & ∂t U t ∂X & ∂t V t ∂Y & & ∂t U ∂t V + t t ∂Y ∂X t& t& ∈XY = ∈Z = t t& & U R (symmetry about the Y axis) (symmetry about the X axis) or t& ∈Z = t & V R The isotropic and orthotropic linear elastic modulus matrices are defined as Isotropic E D= (1 − υ)(1 − 2 υ) LM(1 − υ) MM υ 0 MN υ υ (1 − υ) 0 υ 0 0 (1 − 2 υ) / 2 0 OP P 0 P P (1 − υ)Q υ υ Orthotropic 80 .3. eliminating the radial weighting in the governing equations. and integrate the stresses at the most accurate location. do not lock when incompressible behaviour is being modelled. plastic straining with von Mises plasticity.Two-Dimensional Continuum Elements LM 1 / E −υ / E D=M MM 0 MN− υ / E x xy xz x x 0 − υ yx / E y 1 / Ey 0 0 1 / G xy 0 − υyz / E y − υzx / E z − υzy / E z 0 1 / Ez OP PP PP Q −1 in which symmetry is maintained by defining υ yx = υ xy E y / E x υzx = υ xz E z / E x υzy = υ yz E z / E y Note.4-5).7. The effects of such modes are minimised by the viscous damping technique [H7].3. The rate of diagonal drifting is defined by the velocity at which the mid-points of the element are separating. The location of the integration point is given in Appendix A. This provides elements that are efficient. Element stabilisation The utilisation of one point Gauss quadrature has a limitation in that zero energy deformation or hourglass modes are generated (see fig. 0. The technique provides a damping force capable of preventing the formation of spurious modes but which has negligible influence on the true structural modes. α . α 0 t x y xy . This is possible since the spurious modes are orthogonal to the real deformations. αz T Integration rule for the elements A one point quadrature integration rule is utilised.g. This is utilised as the basis for hourglass detection.2f j =1 4 The viscous hourglassing forces are 81 . α 0 t x y T z d∈ i = ΔT α . α . giving the hourglass velocities as h ij = & ∑ xik Γjk ai = 1. e. To obtain a valid material υ xy < E x / E y d i 1/ 2 υ xz < E x / E z b g υ yz < E y / E zy d i 1/ 2 The initial thermal strain vector is defined as Isotropic Orthotropic d∈ i = ΔT α . and x ik is the nodal velocity of the kth node in the ith direction. This is achieved using an artificial bulk viscosity method. while c. D kk is the trace of the velocity strain tensor and L c is the characteristic length of the element which is related to the smallest element diagonal as Lc = 2A LD where L D = MAX 1 / 2 y 24 e 2 + 1 / 2 x 42 . Shock wave smoothing The shock discontinuities that occur in impact problems may promote numerical instabilities which must be smoothed out.j is given as 82 .1 1 .06 respectively. Q hg is a constant which is modified via the & SYSTEM command and is usually set to a value between 0. The salient characteristic of the method is the augmentation of element pressure with an artificial viscous term (q) prior to the evaluation of the element internal force. and may be modified as necessary via the SYSTEM command.05 and 0. 1 / 2 y31 2 2 + 1 / 2 x13 2 j in which the distance between any two nodal points i. The algorithm has the effect of spreading the shock front over a small number of elements.15. The exact form of artificial viscosity is somewhat arbitrary and the method used is based on the formulation originally proposed in [V1] q = ρ L c D kk Q1 L c D kk + Q 2 c where Q1 and Q 2 are dimensionless constants which default to 1. the material sound speed is defined from c2 = E 1− υ ρ 1 + υ 1 − 2υ a f a fa f The hourglass base vectors for the four node quadrilateral are defined as:Γi = 1 .1 T these viscous forces are included directly into the element force vector.5 and 0.Element Formulations fik = −1 / 4 Q hg ρ A1/2 c ∑ dh ij Γjk i j=1 4 in which A is the current element area. ρ is the current element density. This is zero in expanding elements and non-zero in contracting elements. however. Force calculations The direct stresses at time t+Δt are modified by the addition of the artificial viscosity pressure q as follows σ x = σ x + q and σ y = σ y + q The contribution to the force vector due to the element stresses is evaluated from the equilibrium equations of Timoshenko as Fx = ∂σ x / ∂x + ∂τ xy / ∂y + σ r − σ θ / r = 0 Fy = ∂σ y / ∂y + ∂τ xy / ∂x + τ xy / r = 0 b g Note that the terms σ r r and τ xy r from these two equations are not typically included in static analyses and occur as a result of the inertial effects. is included to control the small spurious oscillations following the shock waves in which the gradients are insufficient to make the quadratic term effective. Evaluation of stresses The element output obtained at the element nodes and Gauss points consists of Stress Output σ X .Two-Dimensional Continuum Elements x ij = x i − x j The quadratic term in strain rate is chosen to be small except in regions of very large gradients. In view of the abundance of excellent results. however. In converging geometries. σ Y . it is generally agreed that the effect is negligible. This occurs even though no shocks are generated and results in a nonphysical generation of pressure. σ min β the direct and shear stresses the maximum and minimum principal stresses the angle between the maximum principal stress and the positive X-axis 83 . σ Z σ max . The hourglass forces are included to give the final force vector. σ XY . The linear term. The mass matrix is computed as each node i as t M x = 1 / 4 t ρt A = 1 / 4 o ρt A t v / o v i tMy i e j = 1 / 4 ρ A = 1 / 4 ρ Ae v / vj t t o t t o where t v is the current volume and o v is the initial volume of an element. the centred strain rate term is negative and the q term is then non-zero. Care should be taken with the linear term since there is a danger of distorting the solution. N min Stress Resultant Output the direct and shear stress resultants/unit length the maximum and minimum principal stress max min resultants/unit length β the angle between the maximum principal stress resultant and the positive X-axis NS the maximum shear stress resultant/unit length NV von Mises equivalent stress resultant/unit length. The Green-Naghdi stress rate formulation is used to refer the constitutive variables to an unrotated 84 .2). Geometrically and materially nonlinear analysis utilising the nonlinear material laws specified in 1. Plain strain and axisymmetry are. Geometrically nonlinear dynamic analysis. All explicit dynamics elements may be used with nonlinear material models 61. i. N X . stress resultants and strain output is shown in fig. N XY . The nodal values of stress and strain are obtained using the extrapolation procedures detailed in section 6.e.4-5. 72. Nonlinear formulation The 2-D explicit dynamics elements can be employed in Materially nonlinear dynamic analysis utilising the elasto-plastic constitutive laws [O2] (section 4.Element Formulations σS σV the maximum shear stress von Mises equivalent stress the direct and shear strains the maximum and minimum principal strains the angle between the maximum principal strain and the positive X-axis the maximum shear strain von Mises equivalent strain Strain Output ∈X . supported. The Gauss point stress is usually more accurate than the nodal values. Eulerian geometric nonlinearity is always invoked with the use of the explicit elements in which the velocity strain measure is utilised. 64. γ XY . Notes The plane stress elements may not be used with nonlinear material model 75. however.3. N Y . ∈min β ∈S ∈V Stress resultant output which accounts for the thickness of the element is available as an alternative to stress output for the plane stress elements. ∈Z ∈max . N Z N max . ∈Y .1. The sign convention for stress.7. 4-2 Example Illustrating The Use Of Plane Stress Elements 85 .4-1 Nodal Configuration For 2d Explicit Dynamics Elements Problem Definition Plate subject to Inplane Loading Finite Element Mesh Fig.Two-Dimensional Continuum Elements configuration prior to the stress integration. The loading is non-conservative. 3 4 3 2 2 1 1 Fig.7. The output is in terms of true Cauchy stresses and the strains approximate to logarithmic strains.3.7.3. Element Formulations Problem Definition Thick Cylinder Finite Element Mesh Fig.3.7.4-4 Example Illustrating The Use Of Axisymmetric Solid Elements 86 .7.4-3 Example Illustrating The Use Of Plane Strain r r Problem Definition Thick Cylinder Finite Element Mesh Fig.3. 3.4-5 Deformed Mesh Illustrating Formation Of Spurious Mechanisms Y σ X.3.3.7.5 Two Phase Plane Strain Continuum Elements (TPN6P and QPN8P) Formulation These isoparametric finite elements utilise the same shape functions to interpolate the displacements and geometry. displacements U= ∑ N i bξ. ηgUi i =1 n n geometry X= ∑ Ni bξ. σ Y +ve tension σ XY +ve into XY quadrant σY σ XY σX σ XY σY σX X Fig.e. i.Two-Dimensional Continuum Elements Fig.7.4-6 Sign Convention For Stress/Strain Output 7. ηgXi i =1 87 . The nodal degrees of freedom are U.7.2.BT mN dv v z z v BT D' B dv 88 . In this instance the two phases comprise the soil skeleton and the pore water fluid. However. the pressures are only interpolated using the corner nodes pressures P= n corner i =1 b g ∑ N bξ. Undrained/fully drained conditions In this type of analysis no consolidation is assumed to take place and the coupled governing equations for static undrained conditions can be expressed as: LMK MNL T R U OPR U = | F − F | SΔP V S−L U − SPV S PT W | | Q T W L ΔU ext T int where the matrices K. ηgP i i where ncorner is the number of corner nodes. L and S are defined as: K= L = . taking into account the loading due to the pore pressure.3.5-1 shows the nodal configurations available within LUSAS. η is the element shape function for node i and n is the number of nodes. The plane strain assumptions and details of elastic modulus matrices applicable for these elements are described in section 7. These elements are used to model the behaviour of a two phase medium such as soil. The finite element method is used to solve the coupled equations in terms of nodal displacements and pore pressures. for consideration of stability.1. coupled by the interaction of the pore pressure and the soil deformation.7). V and P at the corner nodes and U and V at the midside nodes. The consistent and lumped mass matrices are evaluated using the procedures defined in (section 2. Fig. Two plane strain elements QPN8P (quadrilateral) and TPN6P (triangular) based on a mixed displacement-pressure formulation are available in LUSAS to solve these problems. The soil skeleton is analysed in terms of effective stress (total stress minus pore water pressure).2. whilst the pore fluid analysis takes account of the volumetric strain due to the soil skeleton deformation. Separate equations are derived for each phase.Element Formulations where N i ξ. 4) and D' the ‘effective’ soil modulus matrix. Drainage/consolidation process In the drainage/consolidation process. as: H = ∇ N T K ∇ N dv v p z p For nonlinear consolidation process.Two-Dimensional Continuum Elements S= - z v 1 T N N dv Ke K is the tangent stiffness matrix L is the coupling matrix S is the compressibility matrix.0 for backward Euler scheme) H the permeability matrix The permeability matrix H is defined in terms of the shape function derivatives and a permeability matrix of the soil.5.2. fluid flow in/out from the soil needs to be considered. where K e is the equivalent bulk modulus of the soil (see section 7. K . F and int F are external and internal forces Under static fully drained conditions the above coupled governing equations can be further simplified as ext LMK 0OPRΔ UU = R T W T N0 I Q S Δ P V S ext F − int F 0 U V W where I is a unit matrix block. the coupled governing equations can be written as 89 . For linear transient consolidation the coupled governing equations can be expressed as: LMK NML where: T OPRUU βδTH Q TP W PS V L = t + δt LMK NML T OPRUU + RΔ F U (1 − β)δTHQ TP W TΔ Q W PS V S V L t ΔF is the incremental load ΔQ the incremental flow β the time stepping scheme parameter (set to1. Geometrically and materially nonlinear drained/undrained/consolidation analysis utilising the nonlinear material laws specified in 1.3. Therefore the overall compressibility of the soil mass is approximated to be that of the pore fluid. The formulation utilises a nonlinear (spatial) Eulerian formulation. which overcomes the problems of near-incompressibility and effective incompressibility in standard plane-strain and axisymmetric elements.Element Formulations LM K NL T n +1 k k L F ext − n +1 F int ΔUk = Δ tβ H Δ P k − Δt n Q + βΔ Q − ΔtH n +1 P k − L T OPR U R S V | c | QT W S T h U | c U − U hV | W n +1 k n where the superscript on the left/right hand side represents the increment/iteration number. based on the logarithmic strain tensor.2). Geometrically and materially nonlinear dynamic drained/undrained analysis utilising the nonlinear material laws specified in 1. Geometrically nonlinear drained/undrained/consolidation analysis. η (1 − η) 1 η = + ≡ Ke Kf Ks Kf where: K e is the equivalent bulk modulus of the soil K f the bulk modulus of the pore fluid K s the bulk modulus of the solid soil particle η the porosity of the soil In practical geotechnical applications it is usually difficult to determine K f and Ks so a large value of the equivalent modulus K e is usually assumed. QAX4L) Nonlinear formulation These elements are based on a mixed displacement/pressure formulation. associated with the polar decomposition of the deformation gradient F = VR . 7. Material assumptions The bulk modulus of the soil particle Ks is very large compared to the bulk modulus of the pore fluid K f . where V is the left 90 .6 Large-strain Mixed-type Elements (QPN4L. Nonlinear formulation The two phase plane strain continuum elements can be employed in:Materially nonlinear drained/undrained/consolidation analysis utilising the elasto-plastic constitutive laws [O1] (section 4. 1012> K e >109. n 3 ] is the Eulerian triad (spatial orientation of the principal directions) and as τ=n ∑ μ p [λ p =1 N αp − α 1 tr (λ p )I]n T − pJ 3 for the Ogden material model. and as τi = ∑ μ p [λ i p =1 N αp α α 1 α − ( λ1 p + λ 2 p + λ 3 p )] + kJ J − 1 3 for the Ogden material model. k is the bulk modulus and λ i = λ i / 3 J are the “deviatoric” stretches. where N is the number of pairs of Ogden parameters μ p and α p . where λ is the diagonal matrix of deviatoric stretches and n = [ n1. The elements are currently available with Hencky and Ogden matrial models described in section 4.The deformation gradient is given as: F= ∂x ∂X where X and x denote the material and spatial position vector of a material particle. By introducing the independent pressure variable as p = −k J − 1 a f and by transforming τ i from the principal directions the Kirchhoff stress tensor τ is obtained as τ = 2Gn ln λ n T − pJ for the Hencky material model. while k and λ i have the same meaning as for the Hencky model. n 2 .10. so that the principal Kirchhoff stresses τ i = λ i obtained from the corresponding stored-energy function ψ as τ i = 2G ln λ i + kJ J − 1 ∂ψ [C16] are ∂λ i a f a f for the Hencky material model. where J = det F = λ1λ 2 λ 3 .Two-Dimensional Continuum Elements stretch tensor and R is the rotation of the axes of the stretches λ i . where G is the shear modulus. Element equilibrium (mixed formulation) The element equilibrium equations are given as 91 . The Kirchhoff (nominal) stress tensor τ is related to the (true) Cauchy stress σ via τ = J σ . in line with the adopted spatial approach. where R is the vector of applied loading and P is the vector of nodal internal forces. hence integration is still performed over the initial rather than the current volume.tangent stiffness matrix By expanding the element equilibrium into a Taylor’s series. Note that the formulation is defined in terms of the Kirchhoff and not the Cauchy stresses. and the entries in the tangent stiffness matrix are obtained by the consistent linearisation of the element equilibrium. By expressing the stress tensor in the vector form. the vector of nodal internal forces can be written as a f P= V0 z B x τdV0 bg where. In order to derive the subvector K12 and (in particular) submatrix K 11 it helps to regard the vector of nodal internal forces P as coming from the internal virtual power via & aT P = V0 z & ε: τdV0 da & 1 & is the time rate of the nodal displacements. and the second equation follows from p = − k J − 1 .Element Formulations g≡P−R=0 f≡− V0 z LMNa T J −1 + p f k OPQdV 0 =0 where the first equation is the conventional nodal equilibrium equation. Linearisation of the equilibrium . x is the spatial and not the material position vector. the following linearised equilibrium is obtained RδgU ≡ LMK K OPRδa U = −RgU Sδf V M T W MNK K PPQSδpV Sf V T W TW 11 T 12 12 22 where a is the vector of nodal displacements. and ε = ( L + L T ) is the where a = dt 2 strain-rate tensor with & & ∂u L=d= ∂x 92 . subvector K12 immediately follows as ∂τ ij K12 = − B ( x)iJdV0 V0 z T R1U R1U |1| || || with i = S1 V for the plane strain element QPN4L and i = S V for the axisymmetric |0 | |0 | TW |1| TW Piola-Kirchhoff stress S via τ = FSF T and bearing in mind that δ FF −1 = δ d we obtain ∇ ε τδε = F(∇ E Sδ E)F T + δ d τ + τδ d T = δτ + δ d τ + τδ d T . where the repeated indices indicate summation over the dimension of the space. by introducing standard FE ∂p & & matrix/vector notation whereby ε = B( x)a . or in indicial notation T element QAX4L. Also ε: τ = tr( τ T ε ) = ε ijτ ij . For & & configuration-independent loads. By noting the relationship between Kirchhoff stress τ and second ∂τ ij ∂ε kl T tTK δε kl = δτ T. = − Jδ ij so. for both material models. ij + δd ik τ kj + τ ik δd jk = D ijkl δε kl + δd ik τ kj + τ ik δd jk where δτ is called the Truesdell rate of Kirchhoff stress (which is often used in ratedependent constitutive models. here it is introduced because it enables a tTK straightforward formation of the material part of the stiffness matrix) and D ijkl is the tangent constitutive matrix relating the strain-rate tensor to the Truesdell rate of & & & Kirchhoff stress. By using ε ij = (d ij + d ji ) and noting the symmetry of the Kirchhoff & stress tensor τ ij = τ ji . hence & a T ( K δ a + K12δp) = 11 V0 z & & (δε ijτ ij + ε ij ∂τ ij ∂ε kl δε kl )dV0 + V0 z & ε ij ∂τ ij ∂p δpdV0 where. a Tδ g is equal to a Tδ P . the product ε ij & ε ij 1 2 ∂τ ij ∂ε kl δε kl can be written as ∂τ ij ∂ε kl & & & tTK δε kl = ε ijD ijkl δε kl + d ijδd ik τ kj + d jiδd ik τ kj 93 .Two-Dimensional Continuum Elements with L being the so-called velocity gradient and d = ∂u being only introduced for the ∂x & & & sake of convenience during the following derivation. which relates the strain-rate to the tTKE Truesdell rate of Kirchhoff stress can be defined in different ways. An easy way to . via tTK tTKE D ijkl = n ia n jb n kc n ld Dabcd where n ij denotes components of the Eulerian triad n . which relates the strain-rate with define it is by rotating the constitutive matrix D the Truesdell rate of Kirchhoff stress. the variation of which ∂X gives & ∂ ∂u & & δ X = δ dF + dδ F . after noting the symmetry of the Kirchhoff stress tensor. the “normal” components are defined as 94 . and by noting that the variation of the material ∂X ∂X & & & position vector δX is equal to zero. we obtain δ d = − dδ FF −1 = − dδ d . The components of the constitutive matrix D tTKE follow from the stretches and the principal Kirchhoff stresses. where both of these are given with components in the Eulerian frame. the product δε ijτ ij reduces to & & δε ijτ ij = − d ik δd kjτ ij The symmetry of the Kirchhoff stress tensor further implies & ε ij ∂τ ij ∂ε kl & & & tTK δε kl + δε ijτ ij = ε ijD ijkl δε kl + d ijτ kjδd ik so that eventually the submatrix K & a T K δa = 11 V0 11 z 11 follows from & & tTK (ε ijD ijkl δε kl + d ijτ kjδd ik )dV0 Following the standard FE notation.Element Formulations & & By noting that the above-mentioned FF −1 = d yields & ∂u & = dF . By dropping the summation convention. the submatrix K K 11 is then given as = V0 z [BT ( x)D tTK $ B( x) + G T ( x)τ G( x)]dV0 where the tangent constitutive matrix D tTK . which finally F I GH JK gives & δε ij = 1 & 1 & & & (δd ij + δd ji ) = − (d ik δd kj + d jk δd ki ) 2 2 & so that. Varying the second equilibrium equation gives T δf ≡ K12δa + K 22δp = − ∑ N μ pα p 2 λi αp − τ i for the p =1 V0 z LMN δJ + δp dV0 k OP Q where K12 has already been defined and K 22 = − V0 z dV0 k 7.Three-Dimensional Continuum Elements tTKE D iijj = λ j ∂τ i − 2 τ iδ ij ∂λ j where δ ij is the Kronecker symbol and λj ∂τ i 2 = − μ − pJ + 2μδ ij ∂λ j 3 for the Hencky model and λj ∂τ i = ∂λ j ∑ N μ pα p 3 p =1 α α α α α 1 α [3λ i p δ ij + ( λ1 p + λ 2 p + λ 3 p ) − λ i p − λ j p ] − pJ 3 for the Ogden model. HX16. 95 .4 Three-Dimensional Continuum Elements 7. HX20.1 Standard Isoparametric Elements (HX8. PN12.4. i. TH10) Three dimensional isoparametric finite elements utilise the same shape functions to interpolate both the displacements and geometry. TH4. in which case the shear components are given as tTKE tTKE D ijij = D ijji = D tTKE = D tTKE = jiji jiij λ i ∂τ i ∂τ i − τi − 2 ∂λ i ∂λ j F GH I JK which returns the result μ − τ i for the Hencky model and Ogden model. The “shear” components are defined as ( i ≠ j ) tTKE tTKE D ijij = D ijji = D tTKE = D tTKE = jiji jiij λ2j τ i − λ2 τ j i λ2 − λ2 i j unless λ i = λ j .e. PN6. PN15. ηgU i i =1 n n geometry X= ∑ Ni bξ.4. V and W at each node b g The infinitesimal strain-displacement relationship is fully 3-D and is defined as ∈X = ∈Y = ∈Z = ∂U ∂X ∂V ∂Y ∂U ∂Z ∂U ∂V + ∂Y ∂X ∂V ∂ W + ∂Z ∂Y ∂U ∂W + ∂Z ∂X γ XY = γ YZ = γ XZ = The isotropic and orthotropic elastic modulus matrices are Isotropic LMa1 − υf υ υ O 0 0 0 P MM P υ 0 0 0 P a1 − υf υ MM PP 0 0 P MM υ υ a1 − υf 0 E D= P a1 − υfa1 − 2υf M 0 0 0 a1 − 2υf 0 0 P MM P 2 a1 − 2υf 0 PP 0 MM 0 0 0 2 a1 − 2υf PP MM 0 0 0 0 0 2 P N Q 96 . η is the element shape function for node i and n is the number of i nodes.7. ηgXi i =1 where N i ξ. Fig.Element Formulations displacement U= ∑ N i bξ. The nodal degrees of freedom are U.1-1 shows the nodal configurations available within LUSAS. γ YZ . The consistent and lumped mass matrices are evaluated using the procedures defined in (section 2. α . 97 . ∈Z . σ Z . α. γ XZ Principal stresses and strains and the corresponding direction cosines may also be output. and υzy are defined by υ yx = υ xy E y / E x υzx = υ xz E z / E x υzy = υ yz E z / E y to maintain symmetry.Three-Dimensional Continuum Elements Orthotropic LM 1 / E E MM−υ // E −υ D=M MM 0 MM 0 N 0 x xy xz − υyx E y x x − υzx E z 0 0 0 1 / Ey − υyz / E y 0 0 0 − υzy / E z 0 0 0 1 / Ez 0 0 0 0 1 / G xy 0 0 0 0 1 / G yz 0 0 0 0 1 / G xz OP PP PP PP QP −1 where υyx . α xz T A complete description of their formulation is given in [H2. The sign convention for stress and strain output is shown in fig.4. 0.7). α. Evaluation of stresses The element output can be obtained at both the element nodes and Gauss points and consists of Stress Output σ X . σ XZ the direct and shear stresses the direct and shear strains Strain Output ∈X .B1].7. σ YZ . α 0 t x y z xy . 0 0 t T d∈ i = ΔT α . 0. γ XY . σ Y . ∈Y . To obtain a valid material υ xy < E x / E y d i 1/2 υ xz < E x / E z b g 1/2 υ yz < E y / E zy d i 1/ 2 The thermal strain is defined by Isotropic Orthotropic d∈ i = ΔT α. Note. α . σ XY .1-3. υzx . α yz . Notes The nonlinear interface model (section 4.Element Formulations The Gauss point stresses are usually more accurate than the nodal values. HX16.1. Geometrically and materially nonlinear analysis utilising the nonlinear material laws specified in 1. HX20. The geometric nonlinearity may utilise A Total Lagrangian formulation which accounts for large displacements but small strains. The nonlinear strain-displacement relationship is defined by LM OP + 1 LM ∂V OP + 1 LM ∂W OP N Q 2 N ∂X Q 2 N ∂X Q ∂V 1 L ∂U O 1 L ∂V O 1 L ∂W O ∈ = + ∂Y 2 M ∂Y P N Q + 2 MN ∂Y PQ + 2 MN ∂Y PQ ∂W 1 L ∂U O 1 L ∂V O 1 L ∂W O ∈ = + ∂Z 2 M ∂Z P N Q + 2 MN ∂Z PQ + 2 MN ∂Z PQ ∈X = ∂U 1 ∂U + ∂X 2 ∂X 2 2 2 2 Y 2 2 Z 2 2 2 γ XY = γ YZ = γ XZ = ∂U ∂V ∂U ∂U ∂V ∂V ∂W ∂W + + + + ∂Y ∂X ∂X ∂Y ∂X ∂Y ∂X ∂Y ∂V ∂W ∂U ∂U ∂V ∂V ∂W ∂W + + + + ∂Z ∂Y ∂Y ∂Z ∂Y ∂Z ∂Y ∂Z ∂U ∂W ∂U ∂U ∂V ∂V ∂W ∂W + + + + ∂Z ∂X ∂X ∂Z ∂X ∂Z ∂X ∂Z The output is now in terms of the 2nd Piola-Kirchhoff stresses and GreenLagrange strains referred to the undeformed configuration. which takes account of large displacements and moderately large strains provided that the strain increments 98 . Nonlinear formulation The 3-D isoparametric elements can be employed in Materially nonlinear analysis utilising the elastoplastic constitutive laws [O1] (section 4. PN6. An Updated Lagrangian formulation. PN12. Linear eigen-buckling analysis.2). The nodal values of stress and strain are obtained using the extrapolation procedures detailed in section 6.2) may be used with elements HX8. The loading is conservative. Nonlinear dynamics utilising the nonlinear material laws specified in 1. Geometrically nonlinear analysis. The output is in terms of true Cauchy stresses and the strains approximate to logarithmic strains. The loading approximates to being non-conservative. The output is now in terms of the true Cauchy stresses and the strains approximate to logarithmic strains. which takes account of large displacements and large strains.Three-Dimensional Continuum Elements are small. An Eulerian formulation. 99 . The loading is non-conservative. 1-1 Nodal Configuration For Solid Elements 100 .7.Element Formulations 8 5 6 1 2 19 18 16 12 7 6 10 2 3 12 7 8 9 1 2 3 4 10 7 3 1 TH4 HX8 7 4 3 9 16 15 14 12 11 8 7 4 3 6 5 HX16 13 10 1 2 20 6 17 11 5 4 HX20 4 5 1 2 15 10 11 7 1 8 2 3 9 5 4 2 3 TH10 13 9 1 14 15 8 3 PN6 11 10 6 4 PN12 14 13 12 9 6 4 PN15 5 5 8 6 1 2 Fig.4. 1-3 Sign Convention For Stress/Strain Output 101 .7.Three-Dimensional Continuum Elements Fig.7.1-2 Tractor Brake Component Arrows indicate +ve stress directions σY σ XY σ YZ Y σ YZ σ XZ σZ σ XZ σ XY σX X Z Fig.4.4. i. v3 l q.2. Formulation The general approach taken to formulate this element is identical to that described for the 2-D continuum elements in section 7. these conditions also allow the stress field to be eliminated from the formulation. N bζg = 1 e1 − ζ j 2 2 2 2 2 2 3 2 and λ i . particularly if bending predominates. N aηf = 1 e1 − η j. Enhanced strain interpolation The incompatible displacement field is given by u = N1 ξ λ1 + N 2 ξ λ 2 + N3 ζ λ 3 bg bg bg where N1 ξ = b g 1 e1 − ξ j. The formulation provides for the following three conditions to be satisfied Independence of the enhanced and standard strain interpolation functions. Capability of the element to model a constant state of stress after enforcing the orthogonality condition.Element Formulations 7.3. These internal degrees of freedom are eliminated at the element level before assembly of the stiffness matrix for the structure. In addition to ensuring that the element passes the patch test.e. L2 orthogonality of the stress and enhanced strains. represent the incompatible modes λ1 = u1. T 0 1 The covariant base vectors associated with the isoparametric space are g ξ Rx | = Sy |z T Ta 1 Ta 1 Ta 1 U Rx | + η|y V S | |z W T Th 1 Th 1 Th 1 U Rx | + ζ |y V S | |z W T U Rx | | h V + ηζ Sy h | W |z T 3 Tk Tk T U |=g V k| W + ηg1 + ζg1 + ηζ g 1 3 102 . The element is based on a three-field mixed formulation [S8] in which stresses. the element does not suffer from 'locking' in the nearly incompressible limit. T Th T T 3 3 λ 3 = u 3 .4. T λ 2 = u 2 . v2 l q. v1 l q. strains and displacements are represented by three independent functions in three separate vector spaces. requirement for passing the patch test. In addition.2 Enhanced Strain Element (HX8M) The low order enhanced strain element HX8M exhibits improved accuracy in coarse meshes when compared with the parent element HX8. The formulation is based on the inclusion of an assumed 'enhanced' strain field which is related to internal degrees of freedom. Three-Dimensional Continuum Elements g η Rx | = Sy |z T Rx | = Sy |z T Ta Ta Ta 2 2 2 U Rx | + ξ |y V S | |z W T Th 1 Th 1 Th 1 2 2 2 U Rx | + ζ |y V S | |z W T Th Th Th 2 2 2 3 3 U Rx | + ξζ|y V S | |z W T Tk Tk T U |=g V k| W 0 2 + ξ g1 + ζg1 + ξζg 1 2 Ta Ta T g ζ U Rx | + η|y V S a | T W |z 3 3 3 Th Th Th U Rx | + ξ |y V S | |z W T Th Th T U Rx | + ξη|y V S h | W |z T 3 Tk Tk T U |=g V k| W 0 3 + ξg1 + ξ g1 + ηξ g 2 3 where a1 = a2 = a3 = h1 = h2 = h3 = k= 1 −1 8 1 1 −1 − 1 1 1 1 −1 1 1 1 1 T 1 −1 −1 8 1 −1 − 1 1 1 T 1 − 1 −1 − 1 − 1 8 1 8 1 8 1 8 1 −1 1 1 −1 T 1 −1 1 −1 1 1 T 1 −1 −1 −1 −1 1 −1 1 1 T 1 −1 −1 1 −1 1 −1 1 −1 T T T 1 −1 8 1 −1 T x = x1 x 2 x3 x 4 x 5 x 6 x 7 x8 y = y1 y 2 y3 y 4 y 5 y 6 y 7 y8 z = z1 z 2 z3 z 4 z 5 z 6 z 7 z8 T The enhanced covariant strains are given by ∈= ∈ ξξ ∈ ηη ∈ ζζ 2∈ ξη 2∈ ηζ 2∈ ζξ T The enhanced strain field in isoparametric space can initially be expressed using 21-α parameter interpolation functions as follows 103 . i g | | | |eu. j g + du.Element Formulations U Reu. ∈ > L 2 zzz 1 1 1 −1 − 1 − 1 ∈T ∑ dξdηdζ ≡ 0 104 . j g V | | |eu. j g + e u. j g | | | |du. i g | | | | | |eu. The final enhanced strain interpolation matrix is assembled by enforcing the L 2 orthogonality condition < ∑ . j g | | | W T T ξ T ξ η T η i ζ ζ T T η T ξ ξ η T ζ T η η T ζ ξ ζ ζ ξ LMξ MM0 0 =M MM0 MM0 N0 0 0 0 0 0 0 0 0 ξη 0 η 0 0 0 0 0 0 0 0 ηζ 0 ζ 0 0 0 0 0 0 0 0 ξ 0 0 0 η 0 0 ζ 0 0 0 0 0 ζ 0 0 ξ = Ei α i e 0 0 0 0 0 η2 0 0 0 η 0 0 ξ2 0 0 ζξ 0 0 ζ2 ξζ 0 0 ηξ 0 0 0 ξ2 0 η2 0 0 0 0 ζη 0 ζ2 0 0 0 0 ξζ 0 0 0 ηζ 0 0 0 0 ξηζ 0 0 ξηζ 0 ξ 2ζ 0 η2 ξ 0 ηξ ζξ 0 η2 ξ 2 ξη 0 ζη ξ η 0 OP P ξηζ P Pα 0 P ζ ξP P ζ ηP Q 0 0 2 2 i e An element stress field with 12-β parameters is considered: LM1 MM0 0 ∑ = M0 MM MM0 N0 ∑ = 0 1 0 0 0 0 0 0 1 0 0 0 η 0 0 0 0 0 ζ 0 0 0 0 0 0 ξ 0 0 0 0 0 ζ 0 0 0 0 0 0 ξ 0 0 0 0 0 η 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 β 0 12 0 1 OP PP PP PP PQ where the contravariant stresses are defined as ∑ ξξ ∑ ηη ∑ ζζ ∑ ξη ∑ ηζ ∑ ζξ T This stress field is similar to the assumed five β stress field used by Pian [P2] for a hybrid stress quadrilateral element. i g + eu. The field satisfies both equilibrium and symmetry conditions. j g | | ∈ =S | |du. The geometrically nonlinear performance of this element is much improved in comparison with HX8.4.4.2. 7.3 regarding the nonlinear capability of the standard isoparametric element are also applicable to this element.Three-Dimensional Continuum Elements The final interpolation matrix involving eighteen β parameters is E 18 LMξ MM0 0 =M MM0 MM0 N0 0 η 0 0 0 0 0 0 ζ 0 0 0 0 0 0 ξ 0 0 0 0 0 0 η 0 0 0 0 0 0 ζ 0 0 0 η 0 0 0 0 0 0 ζ 0 0 ξη − ξη 0 0 0 0 ξ 2 − η2 0 0 ξ 0 0 ηζ − ηζ 0 2 η − ζ2 0 − ξζ 0 ξζ 0 0 ζ2 − ξ2 0 0 0 ξηζ 0 0 0 0 0 0 ξηζ 0 0 0 0 0 0 ξηζ ξζ ηζ 0 ξ 2ζ η2ζ 0 2 0 ηξ ζξ 0 η ξ ζ2ξ 2 ξη 0 ζη ξ η 0 ζ 2 η OP PP PP PP PQ A further enhanced strain interpolation matrix is also derived which is similar to an interpolation field defined in [S8] for planar elements.4. The nonlinear formulation for the enhanced strain element involves enforcing orthogonality between assumed Green-Lagrange strains and 2nd Piola-Kirchhoff stresses.1. A number of further advantages may also be obtained in explicit dynamic analyses 105 .3 3D Explicit Dynamics Elements (HX8E. TH4E) Explicit time integration schemes have used simple linear elements rather than those of a higher order by virtue of their computational efficiency. This is a requirement for passing the patch test [S8] and is implied in the sense that zzz 1 1 1 −1 − 1 − 1 E dξdηdζ ≡ 0 Evaluation of stresses The evaluation of stresses is identical to that described in section 7. Nonlinear formulation The comments made in section 7. PN6E. This matrix is based on nine a parameters and is also orthogonal to the twelve β stress field E 9 LMξ MM0 0 =M MM0 MM0 N0 0 η 0 0 0 0 0 0 ζ 0 0 0 0 0 0 ξ 0 0 0 0 0 η 0 0 0 0 0 0 η 0 0 0 0 0 ζ 0 0 0 0 0 0 ζ 0 0 0 0 0 ξ OP PP PP PP PQ and E 18 Both the final interpolation functions E 9 also allow condition III to be satisfied.1. It has been shown that higher order continuum elements require a time step reduced from that of linear elements. The combination of mass lumping with linear elements.7. The linear explicit dynamics elements have been implemented to take advantage of these benefits. The nodal degrees of freedom are U.e. A rate relationship is used to define the straindisplacement characteristics as t ∈x = & t ∈y = & t ∈z = & tγ & & ∂t U t ∂X & ∂t V t ∂Y & ∂t W t ∂Z b g XY = & & ∂t U ∂t V + t t ∂Y ∂X 106 . increases accuracy in solutions by virtue of their respective compensatory spectral errors. ηgU i i =1 n n geometry X= ∑ Ni bξ. They are for use only with the explicit central difference time integration scheme. The mass lumping formulations for higher order elements are currently impractical for modelling shock wave propagation since the resulting numerical noise pollutes or destroys the solution.3-1 shows the nodal configurations available within LUSAS. Evaluation of current strain increments The velocity strain rates e t+Dt/2are defined from the midpoint velocity ij gradients in the global axis system. The explicit dynamics elements are based upon the isoparametric approach in which the same shape functions are used to interpolate both the displacements and geometry. ηgXi i =1 where N i ξ. when applied in conjunction with the central difference operator.Element Formulations The use of higher order shape functions creates difficulties at the contact interface in the form of uncontrolled overlap. i. η is the element shape function for node i and n is the number of nodes. V and W at each node.4. Fig. displacement U= ∑ N i bξ. α. α.Three-Dimensional Continuum Elements tγ & tγ & = = & & ∂t V ∂t W + t t ∂Z ∂Y & & ∂t U ∂t W + t t ∂Z ∂X YZ XZ Evaluation of modulus matrices The isotropic and orthotropic elastic modulus matrices are as follows Isotropic LMa1 − υf υ υ O 0 0 0 P MM P υ 0 0 0 P a1 − υf υ MM PP υ υ 0 0 P a1 − υf 0 MM E D= P a1 − υfa1 − 2υf M 0 0 0 a1 − 2υf 0 0 P MM P 2 a1 − 2υf 0 PP 0 MM 0 0 0 2 a1 − 2υf PP MM 0 0 0 0 0 P 2 Q N Orthotropic LM 1 / E E MM−υ // E −υ D=M MM 0 MM 0 N 0 x xy xz − υ yx E y x x − υzx E z − υzy / E z 1 / Ez 0 0 0 0 0 0 1 / G xy 0 0 0 0 0 0 1 / G yz 0 0 0 0 0 0 1 / G xz 1 / Ey − υ yz / E y 0 0 0 OP PP PP PP QP −1 to maintain symmetry the following relations are utilised υ yx = υ xy E y / E x υzx = υ xz E z / E x υzy = υ yz E z / E y Note that a valid material is obtained only if υ xy < E x / E y d i 1/2 υ xz < E x / E z b g 1/2 υ yz < E y / E zy d i 1/ 2 The initial thermal strain is defined by Isotropic d∈ i = ΔT α. 0. 0. 0 0 t T 107 . α xz T The lumped mass matrix is computed as each node i as t t M x i = 1 / 8t ρt V M y i = 1 / 8t ρ t V where t v is the current volume of an element. This provides elements that are efficient and do not lock when incompressible behaviour is being modelled. ρ is the current element density.3f j =1 4 The viscous hourglassing forces are fik = − Q hg ρ v2/3 c / 4 e ∑ dh ij Γjk i M1 + 100 Qhg ∑ dh ij Γjk iP M P j=1 j =1 4 L N 4 O Q in which ve is the current element volume. The technique provides a damping force capable of preventing the formation of spurious modes but which has negligible influence on the true structural modes. and x ik is th th the nodal velocity of the k node in the i direction. α . This is possible since the spurious modes are orthogonal to the real deformations.g. α . Q hg is a constant which is modified via & the SYSTEM command and is usually set to a value between 0. while c. The location of the integration point is given in Appendix I. Element stabilisation The utilisation of one point Gauss quadrature has a limitation in that zero energy deformation or hourglass modes are generated (see Fig. giving the hourglass velocities as h ij = & ∑ xik Γjk ai = 1.05 and 0. The effects of such modes are minimised by the viscous damping technique [H7].3.15. Integration rule for the elements A one point quadrature integration rule is utilised.3-5). The rate of diagonal drifting is defined by the velocity at which the mid-points of the element are separating. The stresses are integrated at the most accurate location. α 0 t x y z xy .Element Formulations Orthotropic d∈ i = ΔT α . This is utilised as the basis for hourglass detection.7. e. the material sound speed is defined from 108 . α yz . plastic straining with von Mises plasticity. This is zero in expanding elements and non-zero in contracting elements. This is achieved using an artificial bulk viscosity method.0 a f ∂x ∂x * ∂ξ ∂η The surface Jacobian J may be evaluated from J 0. D kk is the trace of the velocity strain tensor and L c is the characteristic length of the element which is related to the smallest element diagonal as Lc = V Af where V is the current element volume and A f the current largest face area of the element.5 and 0. The exact form of artificial viscosity is somewhat arbitrary and the method used is based on the formulation originally proposed in [V1] q = ρ L c D kk Q1 L c D kk + Q 2 c where Q1 and Q 2 are dimensionless constants which default to 1. The salient characteristic of the method is the augmentation of element pressure with an artificial viscous term (q) prior to the evaluation of the element internal force. The algorithm has the effect of spreading the shock front over a small number of elements. and may be modified as necessary via the SYSTEM command.Three-Dimensional Continuum Elements c2 = E 1− υ ρ 1 + υ 1 − 2υ a f a fa f −1 1 −1 1 −1 1 −1 1 −1 −1 −1 −1 1 1 −1 −1 1 −1 1 1 −1 −1 1 −1 −1 1 −1 1 The hourglass base vectors Γij for the 8 node solid elements are given as Γij LM1 1 =M MM1 N1 OP PP PQ T these viscous forces are included directly into the element force vector. The face area is evaluated by considering each face in turn and using A f = 4 J 0. 109 .06 respectively. Shock wave smoothing The shock discontinuities that occur in impact problems may promote numerical instabilities which must be smoothed out.0 = a f in which the differentials are evaluated explicitly. The linear term. is included to control the small spurious oscillations following the shock waves in which the gradients are insufficient to make the quadratic term effective. 110 . Care should be taken with the linear term since there is a danger of distorting the solution. σ YZ .) Geometrically nonlinear dynamic analysis.Element Formulations The quadratic term in strain rate is chosen to be small except in regions of very large gradients. it is generally agreed that the effect is negligible. σ Y . Geometrically and materially nonlinear analysis utilising the nonlinear material laws specified in note I. σ XZ the direct and shear stresses the direct and shear strains Strain Output ∈X . In view of the abundance of excellent results. ∈Z .2). The output is in terms of true Cauchy stresses and the strains approximate to logarithmic strains. The loading is non-conservative. γ XZ Principal stresses and strains and the corresponding direction cosines may also be output. ∈Y . σ XY . however. This occurs even though no shocks are generated and results in a nonphysical generation of pressure. 72 and 75 (section 4. however. In converging geometries the centred strain rate term is negative and the q term is then non-zero. γ XY . The direct stresses at time t+Δt are modified by the addition of the artificial viscosity pressure q as follows σ ii = σ ii + q Nonlinear formulation The 3-D explicit dynamics elements can be employed in Materially nonlinear dynamic analysis (see note 1. Eulerian geometric nonlinearity is always invoked with the use of the explicit elements in which the velocity strain measure is utilised. σ Z . Notes The 3D explicit dynamics elements may be used with nonlinear material models 61 to 64. γ YZ . The Jaumann stress rate formulation is used to eliminate rigid motion prior to stress integration. Evaluation of stresses The element output can be obtained at both the element nodes and Gauss points and consists of Stress Output σ X . 7.4. The nodal values of stress and strain are obtained using the extrapolation procedures detailed in section 6. The Gauss point stress is usually more accurate than the nodal values.4.1.3-3.3-1 Nodal Configuration For 3d Explicit Dynamics Elements 111 .Three-Dimensional Continuum Elements The sign convention for stress and strain output is shown in fig. 8 5 6 1 2 HX8E 7 4 3 6 4 5 1 2 PN6E 3 4 3 1 2 TH4E Fig.7. 7.3-2 Compact Tension Fracture Specimen Arrows indicate +ve stress directions σY σ XY σ YZ Y σ YZ σ XZ σZ σ XZ σ XY σX X Z Fig.7.4.4.3-3 Sign Convention For Stress/Strain Output 112 .Element Formulations HX8E Elements PN6E Elements Fig. 4. or a single point (for PN6L). This limitation requires the calculation of only a 2x2 Jacobian matrix. In order to speed up the computation. are defined in terms of natural coordinates ξ and η. ζ + 1 N T top g ia ib 2 b g b g The in-plane and through-thickness shape functions can then be separated to give: N T = φ T + ζψ T br where φT = 1 T T N .4 Composite Solid Elements (HX8L. To form the complete shape functions for the brick element N br . and PN6L elements use the equivalent shape functions for 6. To overcome this difficulty layered brick elements were developed where several laminae are included in a single element. or a 3 corner point quadratic (for PN12L) gauss integration scheme outside the through thickness loop. The shape functions for the top and bottom surfaces of the composite elements can be considered to be single membrane element shape functions. 4 and 3 noded membranes.4. For these elements the three degrees of freedom per node are used to interpolate a displacement field that varies linearly over the thickness and quadratically in-plane for the higher order elements. see figure 7.N 2 i i 113 .HX16L. Each layer is specified by an orthotropic material stiffness matrix. The shape functions N i ( top) = N i ( bot ) = N i . linear interpolation is used between the functions for the top and bottom surfaces: NT = br 1 −ζ + 1 N T bot f . while the strain-displacement matrices by default are integrated using a plane 2x2 (for HX8L and HX16L).Three-Dimensional Continuum Elements 7.4-1.PN6L PN12L) If brick elements are used for an analysis of composite structures the number of degrees of freedom even for small laminate structures rapidly becomes very large leading to prohibitively excessive computer costs. for the HX16L element these are given by: Ni = Ni = 1 1 + ξ i ξ 1 + ηi η ξ i ξ + ηi η − 1 4 1 1 − η2 ξ 2 − ξ 2 η 1 + ξ i ξ + ηi η i i 2 b gb gb g for corner nodes for mid-side nodes e jb g The PN12L. the elements are restricted to reasonably constant layer thicknesses [H13]. For the integration of the element stiffness matrix. the material stiffnesses are summed layerwise through the thickness. HX8L. ... w ...... v ... ∂u + ∂w U S ∂x ∂y ∂z ∂y ∂x ∂z ∂y ∂z ∂x V T W The strain displacement relationship is given by: ∈= Ba where B is the strain displacement matrix...... can now be interpolated as: LMφ U=M MNM T + ζψ T 0T 0T 0T φ T + ζψ T 0T 0T 0T φ T + ζψ T OPR u U SV PP| v | TW QP|w| U = Ha with the displacement vectors in terms of the nodal degrees of freedom: l q v = lv ..... u 2 ...... ∂v + ∂w . ∂v .. U.Element Formulations ψT = 1 −NT ....... u n 1 2 n 1 2 n T T T The three-dimensional strain vector ∈ is defined as ∈T = R ∂u .. NT i i 2 The displacement field. ∂u + ∂v .. w q u = u1.... ∂w .... LM ∂φ + ζ ∂ψ MM ∂x ∂x MM 0 MM M 0 B= M MM ∂φ + ζ ∂ψ MM ∂y ∂y MM 0 MM 2 ψ N c T T T T T T T 0T ∂φ T ∂y +ζ 0T ∂φ T ∂x +ζ ∂ψ T ∂x ∂ψ T ∂y T 2 T ψ c 0T OP PP PP 0 PP 2 ψ PP c PP 0 ∂ψ P ∂φ P +ζ ∂y ∂y P P ∂ψ P ∂φ +ζ ∂x P ∂x Q 0T T T T T T T T B can be split into two matrices combining in-plane and through thickness terms: 114 . v q w = lw . Consequently for the transformation of the cartesian derivatives into the natural derivatives only a 2 dimensional Jacobian matrix is required. R ∂ U LM ∂x | ∂ξ | = M ∂ξ S ∂ V M ∂x | ∂η | M ∂η | | T W N or inverted ∂y ∂ξ ∂y ∂η OPR ∂ U PP| ∂∂x | S V | | PQ| ∂y | T W ∂ −1 ∂ −1 ∂ = J11 + J12 ∂η ∂x ∂ξ ∂ −1 ∂ −1 ∂ = J 21 + J 22 ∂η ∂y ∂ξ and an integration constant for the thickness is computed from: z= c ∂ 2 ∂ ζ→ = . The differential of the volume is given by dV = c J ζdξdη 2 115 .Three-Dimensional Continuum Elements B = B + ζB 1 2 where B 1 LM ∂φ MM ∂x MM 0 MM M0 =M MM ∂φ MM ∂y MM 0 MM 2 ψ Nc T T T T 0T ∂φ T ∂y 0T ∂φ T ∂x 2 T ψ c 0T 0T 2 T ψ c 0T ∂φ T ∂y ∂φ T ∂x T T 0T OP PP PP PP PP PP PP PP PP Q B 2 LM ∂ψ MM ∂x MM 0 MM M0 =M MM ∂ψ MM ∂y MM 0 MM 0 N T T T T T 0T ∂ψ T ∂y 0T ∂ψ T ∂x 0T 0T T OP PP 0 P PP 0 P PP 0 P P ∂ψ P P ∂y P P ∂ψ P ∂x P Q 0T T T T T T The restriction of constant layer thicknesses provides an uncoupling between the inplane coordinates and the through-thickness coordinate.4-1.4. 2 ∂z c ∂ζ where c is the depth of the element see figure 7. As the matrices B and B are independent of ζ.Element Formulations where |J| is the Jacobian determinant. only D varies from layer to layer. LM 1 / E E MM−υ // E −υ D=M MM 0 MM 0 N 0 x xy xz x x − υyx E y 1 / Ey − υ yz / E y 0 0 0 − υzx E z − υzy / E z 1 / Ez 0 0 0 0 0 0 1 / G xy 0 0 0 0 0 0 1 / G yz 0 0 0 0 0 0 1 / G xz OP−1 PP PP PP QP where υyx . The element stiffness matrix in basic form may be defined as K= z V BT DB dV where D is the modulus matrix for an orthotropic material. υzx and υzy are defined by υ yx = υ xy E y / E x υzx = υ xz E z / E x υzy = υ yz E z / E y to maintain symmetry. 1 2 Therefore the strain-displacement matrices can be left out of the integration through the thickness: K= zz 2 ξ η F GG B GG GG H T 1 LM MN ∑ nlay n =1 z ζlay D dζ B + B T n 1 1 nlay n =1 + BT 2 O LM ∑ ζ D dζ P B MN PQ LM O L ∑ ζ D dζ P B + B M ∑ MN PQ MN OP PQ nlay n =1 z ζlay n 2 z nlay n =1 ζlay n 1 T 2 z ζlay I JJ c JJ 2 J dηdξ O ζ D dζ P B J PQ JK 2 n 2 with B and B as: 1 116 . 2).Three-Dimensional Continuum Elements LMJ MM MM MM M B =M MMJ MM MM MM N LMJ MM MM MM M B =M MMJ MM MM MM N 1 2 −1 11 ∂φ T ∂ξ −1 + J12 ∂φ T ∂η −1 J 21 0T ∂φ T ∂ξ −1 + J 22 0T ∂φ T ∂η 0T 2 T ψ c 0T −1 J 21 −1 J11 0T 0T ∂φ T ∂ξ −1 + J 22 0T ∂φ T ∂η −1 J11 −1 21 ∂φ T 0T 2 T ψ c ∂ψ T ∂ξ −1 + J12 ∂ξ 2 T ψ c 0T −1 + J12 ∂φ T ∂η ∂φ T ∂ξ ∂φ T ∂ξ −1 + J 22 −1 + J12 OP PP PP PP PP PP P ∂φ P ∂η P P ∂φ P ∂η P Q T T −1 11 ∂ψ T ∂η −1 J 21 0T ∂ψ T ∂ξ −1 + J 22 0T ∂ψ T ∂η 0T 0T ∂ψ T ∂η −1 J 21 0T 0T −1 21 0T ∂ψ T ∂η −1 J11 ∂ψ T ∂ξ −1 + J 22 ∂ψ T ∂ξ −1 + J12 0T ∂ψ T −1 + J 22 −1 + J12 0T 0T 0T 0T J11 ∂ξ ∂ψ T −1 ∂ξ OP PP PP PP PP PP P ∂ψ P ∂η P P ∂ψ P ∂η P Q T T The through thickness dependency is condensed in the integration of the material modulus matrix which makes the assembly of the element stiffness matrix more efficient.2). 117 . The strain displacement matrices only have to be computed in-plane. Nonlinear formulation The 3-D solid composite elements can be employed in Materially nonlinear analysis utilising the elastoplastic constitutive laws [O1] (section 4. Geometrically and materially nonlinear analysis utilising the nonlinear material laws specified in 1.5. Geometrically nonlinear analysis utilising the corotational formulation (section 3. This is possible by restricting the element to a reasonably uniform thickness for a single layer. 4. Fig. ηgXi i =1 118 .Element Formulations Nonlinear dynamics utilising the nonlinear material laws specified in 1. displacements U= ∑ N i bξ.5 Two Phase 3D Continuum Elements (TH10P. ηgUi i =1 n n geometry X= ∑ Ni bξ.7. PN12P. Linear eigen-buckling analysis.4-1 Topology Of 3-D Layered Isoparametric Bricks 7. HX16P and HX20P) Formulation These isoparametric finite elements utilise the same shape functions to interpolate the displacements and geometry. PN15P.e. i.4. the pressures are only interpolated using the corner nodes pressures P= n corner i =1 b g ∑ N bξ. Separate equations are derived for each phase.5-1 shows the nodal configurations available within LUSAS. coupled by the interaction of the pore pressure and the soil deformation. In this instance the two phases comprise the soil skeleton and the pore water fluid. Fig. These elements are used to model the behaviour of a two phase medium such as soil. The soil skeleton is analysed in terms of effective stress (total stress minus pore water pressure). PN12P. taking into account the loading due to the pore pressure. HX16P and HX20P based on a mixed displacement-pressure formulation are available in LUSAS to solve these problems.4.7). The finite element method is used to solve the coupled equations in terms of nodal displacements and pore pressures. whilst the pore fluid analysis takes account of the volumetric strain due to the soil skeleton deformation. ηgP i i where ncorner is the number of corner nodes. PN15P.BT mN dv v z z v BT D' B dv 119 . W and P at the corner nodes and U.7. Undrained/fully drained conditions In this type of analysis no consolidation is assumed to take place and the coupled governing equations for static undrained conditions can be expressed as: R U LMK LOPRΔUU = | F − F | SΔP V S−L U − SPV | NML S PQT W | T W ext T int T where the matrices K. Five 3D elements TH10P. The consistent and lumped mass matrices are evaluated using the procedures defined in (section 2.3. However. for consideration of stability.1. V and W at the midside nodes. η is the element shape function for node i and n is the number of nodes. The nodal degrees of freedom are U. The details of elastic modulus matrices applicable for these elements are described in section 7. L and S are defined as: K= L = . V.Three-Dimensional Continuum Elements where N i ξ. where K e is the equivalent bulk modulus of the soil (see section 7.3. as: H = ∇ N T K ∇ N dv v p z p For nonlinear consolidation. Drainage/consolidation process In the drainage/consolidation process. the coupled governing equations can be written as 120 . K . fluid flow in/out from the soil needs to be considered. For linear transient consolidation the coupled governing equations can be expressed as: LMK NML where: T OPRUU βδTH Q TP W PS V L = t + δt LMK NML T OPRUU + RΔ F U (1 − β)δTHQ TP W TΔ Q W PS V S V L t ΔF is the incremental load ΔQ the incremental flow β the time stepping scheme parameter (set to 1. F and int F are external and internal equivalent nodal forces Under static fully drained conditions the coupled governing equations can be further simplified as ext LMK 0OPRΔ UU = R F− FU V T W T N0 I Q SΔ P V S 0 W ext int where I is a unit matrix block.5.Element Formulations S= - z v 1 T N N dv Ke K is the tangent stiffness matrix L is the coupling matrix S is the compressibility matrix.4) and D' the ‘effective’ soil modulus matrix.0 for backward Euler scheme) H the permeability matrix The permeability matrix H is defined in terms of the shape function derivatives and a permeability matrix of the soil. 121 . Material assumptions The bulk modulus of the soil particle Ks is very large compared to the bulk modulus of the pore fluid K f . Nonlinear formulation The two phase 3-D continuum elements can be employed in Materially nonlinear drained/undrained/consolidation analysis utilising the elastoplastic constitutive laws [O1] (section 4. Geometrically nonlinear drained/undrained/consolidation analysis. Geometrically and materially nonlinear drained/undrained/consolidation analysis utilising the nonlinear material laws specified in 1. 1012> K e >109.2). Therefore the overall compressibility of the soil mass is approximated to be that of the pore fluid. Geometrically and materially nonlinear dynamic drained/undrained analysis utilising the nonlinear material laws specified in 1. η (1 − η) 1 η = + ≡ Ke Kf Ks Kf where: K e is the equivalent bulk modulus of the soil K f the bulk modulus of the pore fluid K s the bulk modulus of the solid soil particle η the porosity of the soil In practical geotechnical applications it is usually difficult to determine K f and Ks so a large value of the equivalent modulus K e is usually assumed.Three-Dimensional Continuum Elements LM K NL T F OPRΔ U U = R S S Δ P V |− Δtc Q + βΔQh − ΔtH− Δ tβ H Q T | W T L k n +1 k ext k n k F int n +1 k P − LT n +1 U | c U − U hV | W n +1 k n where the superscript on the left/right hand side represents the increment/iteration number. 7.5.7).1-1) U and V at each node The infinitesimal strain-displacement relationship is defined in the local Cartesian system by ∈x = ∂u ∂x ∈z = U R The elastic modulus matrix is defined by D= 1 υ E 1 − υ2 υ 1 LM N OP Q The thermal strain is defined by b∈ g = ΔT α. The local y and z axes form a right-hand set with the x122 . α o t T The consistent and lumped mass matrices are evaluated using the procedures defined in (section 2. The nodal degrees of freedom are (fig.Element Formulations 7.5. BXM3) Formulation BXM2 and BXM3 elements are axisymmetric.1 Axisymmetric Membrane (BXM2. Evaluation of stresses The element output obtained at the element nodes and Gauss points consists of σx σz ∈x ∈z Stress Output Meridional stress (+ve tension) Circumferential stress (+ve tension) Strain Output Meridional strain (+ve tension) Meridional stress (+ve tension) The element local x-axis lies along the element axis in the direction in which the element nodes are specified. isoparametric membrane elements.5 Space Membrane Elements 7. They are defined in the XY-plane and symmetry may be specified about either the X or Y axes. 1. Nonlinear formulation The axisymmetric membrane elements can be employed in Materially nonlinear analysis utilising the elasto-plastic constitutive laws [O1] (section 4. Nonlinear dynamics utilising the nonlinear material laws specified in 1.V 1 2 1 2 BXM3 3 BXM2 X. Geometrically nonlinear analysis.Space Membrane Elements axis such that the y-axis lies in the global XY-plane and the z-axis is parallel to the global Z-axis (up out of page) (fig. Linear eigen-buckling analysis.1-4). The nonlinear straindisplacement relationship is defined by 2 + 1 ∂u 2 ∂x LM OP N Q 2 The output is now in terms of the 2nd Piola-Kirchhoff stresses and GreenLagrange strains referred to the undeformed configuration.5. Notes LM OP N Q U 1 LUO ∈= + M P R 2 NRQ ∈x = ∂u 1 ∂u + ∂x 2 ∂x 2 x The geometric nonlinearity utilises a Total Lagrangian formulation which accounts for large displacements but small strains.7. The nodal values of stress and strain are obtained using the extrapolation procedures detailed in section 6. The loading is conservative. The Gauss point stresses are usually more accurate than the nodal values.U 123 . Geometrically and materially nonlinear analysis utilising the nonlinear material laws specified in 1.2). Y. 1-2 'Stand-Alone' Applications For BXM2 And BXM3 Elements 124 .5.7.Element Formulations Fig.5.1-1 Nodal Configuration For BXM2 And BXM3 Elements Problem Definition (b) Circular Plate Finite Element Mesh Problem Definition (b) Circular Pipe Finite Element Mesh Fig.7. 7.5.7.5.Space Membrane Elements BXM3 elements QAX8 elements Problem Definition Finite Element Mesh Fig.1-4 Local Cartesian System For Bxm2 And Bxm3 Elements 125 .1-3 Fibre Reinforced Cylinder Illustrating Coupling Between QAX8 And BXM3 Elements y y Y 2 y x 1 y x 3 x y x 2 x 1 X Fig. 7).Element Formulations 7. Nxy the direct and shear stress resultants/unit length Nmax.Nmin the maximum and minimum principal stress resultants/ unit length b the angle between the maximum principal stress resultant and the positive X-axis Strain Output ∈x .7. Their formulations are exactly equivalent to their 2-D conterparts given in table 7. The nodal degrees of freedom are U. The stress resultants are evaluated directly at the nodes. ∈min the maximum and minimum principal strains β the angle between the maximum principal strain and the positive X-axis The sign convention for stress resultant and strain output is shown in fig. but may be utilized in a nonlinear environment.21 Space Membrane Plane Membrane SMI4 TSM3 PMI4 TPM3 TABLE 7.5.5. ∈y .2-4.2-1 Space Membrane Elements And Equivalent Plane Elements The nodal configurations are shown in fig. TSM3) Formulation SMI4 and TSM3 elements are membrane elements that function in 3-D. They are formulated in 2-D. Ny . V and W at each node Only a lumped mass matrix is evaluated using the procedure defined in (section 2.7.5. The element cannot be used for linear buckling analyses. by forming a local Cartesian system in the plane of the element (using a least squares fit through the element nodes).5. Once the element matrices have been formed they are then transformed to the global Cartesian basis.2-1. Evaluation of stresses The element output obtained at the element nodes consists of Stress Resultant Output Nx .5. γ xy the direct and shear strains ∈max . Nonlinear formulation The element has no nonlinear capability. 126 .2 3-D Space Membrane (SMI4. V 3 4 1 TSM3 2 1 SMI4 2 X.5.2-1 Nodal Configuration For SMI4 And TSM3 Elements Thin membrane SMI4 Elements Stiffening members QSI4 elements Problem definition Finite Element Mesh Fig. W Fig.7.2-2 Box Structure Illustrating The Use Of Space Membrane Elements 127 .5.7.Space Membrane Elements 3 Y. U Z. 2-3 Local Cartesian System For SMI4 And TSM3 Elements Y σ X.5.Element Formulations 4 Y z y 3 x 1 2 X Z Fig.6.6. σ Y +ve tension σ XY +ve into XY quadrant σY σ XY σX σ XY σY σX X Fig.7.5.1-1) 128 .7.1 Isoflex Thin Plate (QF4.2-4 Sign Convention For Stress/Strain Output 7. TF6) Formulation The Isoflex family of thin plate elements are formed by applying Kirchhoff constraints within elements formulated using Mindlin plate assumptions. TF3. QF8.6 Plate Elements 7.7. The displacements and rotations are considered independent and the unconstrained nodal configurations are (fig. These constraints are sufficient for the higher order elements and the extra constraints required for the lower order elements are provided by enforcing a linear variation of tangential rotation along the element sides. This provides 1 constraint suitable for removing the central translation of the quadrilaterals.1-3) w. Where γ t is the through-thickness shear strain tangential to the element edges. Δθ y at the corner at the mid side nodes of the quadrilateral.6. which are suitable for removing the rotations at the central node. This provides 6 and 8 constraints respectively for the triangles and quadrilaterals which are suitable for eliminating the mid-side translation and normal rotation.7.7. S γ n dS = 0 Where γ n is the transverse shear strain normal to the element sides and the integral is performed using 2-point quadrature along each side. z z A γ XZ dA = 0 . This removes 8 and 11 degrees of freedom for the 6 and 3 noded triangles and 11 and 15 degrees of freedom for the 8 and 4 noded quadrilaterals respectively.6. θ x . originally proposed by Irons for the Semiloof shell [I1] γt = ∂w − θy = 0 ∂x At the points shown in fig. These constraints provide extra equations that permit certain nodal degrees of freedom to be discarded.1-4. Δθ x . The final nodal configurations are (fig.7. This is achieved by using the following constraints. This provides 2 constraints for both the triangles and quadrilaterals. z A γ YZ dA = 0 Where the integral is performed using 2*2 Gauss quadrature. θ y Δθ at the corner at the mid side nodes where Δθ is the relative rotation about a tangent to the element edge.Plate Elements w. θ y Δw. at the central node of the triangle.1-2). θ x . where Δθ and θ re the relative (departure from linearity) and absolute rotations of the through-thickness normals after deformation.6. An element with thin plate performance is then produced by constraining the shear strains to zero at discrete points within the element. 129 . These rotations include the transverse shear deformations (fig. Δθ y Δθ x . ηg t i θY i =1 n i V=− z t n ∑ N i bξ. ηg Wi i =1 where t and t i are the thicknesses of the plate at the integration and nodal points respectively. and N(ξ. Therefore the discretised.η) are the element shape functions. so that ∈X = ∈Y = ∂U ∂X ∂V ∂Y ∂U ∂V + ∂Y ∂X γ XY = The continuum displacements for plates of varying thickness are related to the original degrees of freedom of the plate using U= z t ∑ Ni bξ. flexural strain-displacement relationship is LM ψ MM ψ Nψ X Y XY LM 1 ∂t ∂N X OP M 1t ∂∂X ∂∂N t = ∑M PP M t ∂Y ∂Y Q MM1 ∂t ∂N + 1 ∂t ∂N N t ∂Y ∂X t ∂X ∂Y i n i i =1 i 0 ti t t − i t − ∂N1 ∂Y ∂N1 ∂Y t i ∂N1 t ∂X i OP PPLM W OP 0 PPMMθ PP t ∂N N θ Q t ∂Y P Q i Xi Yi i 1 where the terms involving ∂t / ∂X and ∂t / ∂Y are the small strain contributions due to thickness variations.Element Formulations The infinitesimal strain-displacement relationship is derived from the 3-D continuum relationship [Section 7.4] by neglecting ∈Z which is zero in the Mindlin plate assumptions. ηg t i θX i =1 n i W= ∑ Ni bξ. generalised. For flat plates LM ψ MM ψ Nψ X Y XY LM − ∂ W OP ∂ OP MM ∂ XW PP PP = MM − ∂Y PP Q M−2 ∂ W P MN ∂X∂Y PQ 2 2 2 2 2 130 . and γ XZ and γ YZ which have been constrained to zero. ψ XY the flexural strains in the global Cartesian system. M XY the moments/unit width in the global Cartesian system. 131 . 0.7. α . 0. The nodal stress resultants are evaluated using the extrapolation procedures detailed in section 6. Nonlinear formulation The element has no nonlinear capability. α yz . 0 eψ j = daΔTf α . M Y . α . Evaluation of stresses/strains The element output obtained at the element nodes or Gauss points consists of Stress Resultant Output M X . Note. The element cannot be used for linear buckling analysis. For a valid material υxy < ( E x / E y )1/ 2 The thermal strain is defined by Isotropic Orthotropic Δ eψ j = dadzTf α.1. α xz T where υyx has been set to υxy E y / E x to maintain symmetry. 0 0 t T Δ eψ j = d(dzT) α . The Gauss point values are usually more accurate than the nodal values. α 0 t x y xy . α. α T 0 t 0 t T xy dz x y Full details of the element formulation are given in [L1]. Strain Output ψ X . The sign convention for stress resultant and strain output is shown in fig. Approximate shear forces evaluated by differentiating the moments may also be output.2-6.5. ψ Y .Plate Elements The isotropic and orthotropic elastic resultant modulus or rigidity matrices are Isotropic Orthotropic Δ eψ j = d(dzT) α. but may be utilised in a nonlinear environment.7). Both consistent and lumped mass matrices are available and are evaluated using the procedures defined in (section 2. α. Note. 1-1 Initial Nodal Configuration For Isoflex Plate Element ∂W / ∂X Z θY γ XZ Displacement of any point a distance z along normal is U = z θY where ∂W θY = − γ XZ ∂X θY X Fig.7.6.Element Formulations W θy θx W θy θx W θy θx 1 X Z Quadrilateral Element Triangular Element θx 2 8 θx W θy θx 3 θx 1 9 7 θx W θy 6 θx W θy 4 θx W θy Y W θy θx 2 W θy θx 3 θy 6 θx W θy W θy θx 4 W θy 5 W θy W θy θx 5 W Y Fig.7.6.1-2 Rotation Of The Through-Thickness Normal For A Thick Plate 132 . 6.6.1-3 Final Nodal Configuration For Isoflex Thin Plate Elements 1/ 3 1/ 3 1/ 3 1/ 3 1/ 3 2 2 1/ 3 2 1/ 3 1/ 3 2 1/ 3 1/ 3 2 Fig.1-4 Locations Where The Transverse Shear Strain Tangential To The Element Edge Is Constrained To Zero 133 .7.7.Plate Elements θx W 4 θx θy W 3 θx θy θx W θy W θy θx 1 QF4 W θy θx 2 W θy θx TF3 W θy θx W 4 θx 8 Δθ W Z θx 1 Y X θy 5 QF8 Δθ θy Δθ W 7 θx Δθ 6 W θy W 1 θx θy 4 Δθ Δθ 6 3 θy 3 W θy θx Δθ 5 W θy θx TF6 2 θx 2 Fig. 7.Element Formulations (a) Problem Definition QF4 elements Y X (b) Finite Element Mesh Fig.1-5 Thin Cantilever Plate Illustrating Use Of QF4 Element 134 .6. 6.7.7. θ X .7.6. The extra higher order degrees of freedom are condensed out before assembly so that the final nodal configuration is (fig. This is accomplished by first forming the constrained flexural strain-displacement relationship in exactly the same manner as for the QF4 element.6.1].Plate Elements MXY MY MX MXY MX MXY Z Y MY X MXY Fig.6.2-1. θ Y at the corner nodes The elastic resultant modulus or rigidity matrix is defined as 135 . ∈XZ = P1 ∈XZ1 + P3 ∈XZ 3 and ∈YZ = P2 ∈YZ 2 + P4 ∈YZ 4 where ∈XZi and ∈YZi are the transverse shear strains along the element sides and Pi are linear interpolation functions defined in fig.2 Isoflex Thick Plate (QSC4) Formulation The Isoflex thick plate element QSC4 is formed by imposing an assumed shear strain field on the isoflex thin plate element QF4 [section 7. and then imposing a bilinear shear strain field defined using nodeless degrees of freedom.e.2-2) w.1-6 Sign Convention For Stress Resultant Output 7.6. i. 0.2 0 LM N 0 G xz OP Q where υyx has been set to υxy E y / E x to maintain symmetry.4a1 + υf M0 1 P N Q Et 3 2 s and for orthotropic materials 1 / Ex 3 $ = t −υ / E D yx y b 12 0 LM MM N − υ xy / E x 1 / Ey 0 0 0 1 / G xy OP PP Q −1 and t G yz $ D = s 1. M XY the moments/unit width in the global Cartesian system. α. Evaluation of stresses/strains The element output obtained at the element nodes and Gauss points consists of Stress Resultant Output M X . 0. α 0 t 0 t T dz x y xy yz .7. Both consistent and lumped mass matrices are available and are evaluated using the procedures defined in section 2. α .Element Formulations $ LMD MN 0 0 $ D $ D= b s OP PQ where. α . 0 eψ j = daΔTf α . for isotropic materials LM1 υ 0 OP $ Mυ 1 0 PP D= 12e1 − υ j M MN0 0 a1 − υf PQ 2 Et L1 0 O $ D = 2. α xz T Full details of the element formulation are given in [C4]. For a valid material υxy < ( E x / E y )1/ 2 The thermal strain is defined by Isotropic Orthotropic Δ eψ j = dadzTf α. M Y . 136 . Note. Plate Elements S X . The nodal stress resultants are evaluated using the extrapolation procedures detailed in section 6.7.2-1 Interpolation Functions For Nodeless Freedoms Of The QSC4 Element 137 .7. 3 4 3 1 P5 2 1 P6 6 2 4 3 3 1 P7 2 1 P8 2 Fig.6. ψ XY the flexural strains in the global Cartesian system.6. The sign convention for stress resultant and strain output is shown in fig. Nonlinear formulation The element has no nonlinear capability but may be utilised in a nonlinear environment. The element cannot be used for linear buckling analysis. γ XZ the shear strains in the global Cartesian system. SY the shear forces/unit width in the global Cartesian system.2-4. Strain Output ψ X .1. The Gauss point values are usually more accurate than the nodal values. ψ Y . γ YZ . 2-3 Perferated Thick Plate Example Illustrating Use Of QSC4 Element 138 .7.6.2-2 Nodal Configuration For The QSC4 Element Fig.6.7.Element Formulations θx W 4 θx θy W 3 θx θy W θy Z θx Y X 1 W θy θx 2 Fig. 7. 'normals' to the mid-surface remain straight but not necessarily normal to the mid-surface after deformation (fig. The theory also permits treatment of lateral displacement and rotations as independent variables. so that 139 .7.6. θ Y at each node where θ X and θY are the rotation of the normals to the mid-surface and include the effects of shear deformations. Thus the elements account for the transverse shear effects associated with thicker plates and the elements are termed 'thick' plate elements.6. The infinitesimal. TTF6) Formulation The QTF8 and TTF6 elements are isoparametric plate elements formulated using Mindlin plate theory [M3]. generalized. θ X .3 Isoparametric Thick Mindlin Plate (QTF8. which assumes that Normal stress in the transverse stress is negligible in comparison with the in plane stresses. The nodal degrees of freedom are (fig.7.Plate Elements MXY MY MX MXY SY MX MXY SX Z Y SX MY SY X MXY Fig. producing elements which only require C(0) continuity.6.6.2-4 Sign Convention For Stress Resultant Output 7. flexural straindisplacement relationship is derived from the 3-D continuum strain-displacement relationship by neglecting the out of plane strain.3-1).3-2) W. 2 0 LM N 0 G xz OP Q where υyx has been set to υxy E y / E x to maintain symmetry. For a valid material υxy < ( E x / E y )1/ 2 140 .Element Formulations ψX = ∂θ Y ∂X ∂θ X ∂Y ∂θ Y ∂θ Y − ∂Y ∂X ∂W − θX ∂Y ∂W + θY ∂Y ψY = − ψ XY = γ YZ = γ XZ = The elastic resultant modulus or rigidity matrix is defined as $ D= $ LMD MN 0 b 0 $ D s OP PQ LM MM N 0 G xz 1 υ 0 υ 1 0 where.2 0 LM N OP Q − υ xy / E x 1 / Ey 0 0 0 1 / G xy and for orthotropic materials 1 / Ex 3 $ = t −υ / E D xy x b 12 0 LM MM N OP−1 PP Q and t G yz $ D = s 1. Note. for isotropic materials Et 3 $ D = b 12(1 − υ2 ) OP−1 0 PP (1 − υ) / 2 Q 0 and t G yz $ D = s 1. the flexural strains in the global Cartesian system. 0. ψ XY . SY . M XY . α.Plate Elements The thermal strain is defined by Isotropic Orthotropic Δ eψ j = dadzTf α. The Gauss point values are usually more accurate than the nodal values. 141 .3-4.6. Strain Output ψ X . α 0 t 0 t T dz x y xy yz . 0. The elements cannot be used for linear buckling analysis. α xz T Both consistent and lumped mass matrices are available and are evaluated using the procedures defined in section 2. 0 eψ j = daΔTf α .7.the bending strains in the global Cartesian system.7.the moments/unit width in the global Cartesian system. α . Evaluation of stresses/strains The element output obtained at the element nodes and Gauss points consists of Stress Resultant Output M X . Nonlinear formulation The element has no nonlinear capability but may be utilised in a nonlinear environment. The sign convention for stress resultant and strain output is shown in fig. α . γ YZ . ψ Y . The nodal stress resultants are evaluated using the extrapolation procedures detailed in section 6.the shear forces/unit width in the global Cartesian system. M Y . γ XZ . S X .1. Element Formulations Displacement of any point a distance z along normal is U = z θY where ∂W − γ XZ θY = ∂X ∂W / ∂X Z θY γ XZ θY X Fig.7.6.6.3-1 Rotation Of The Through-Thickness Normal For A Thick Plate W θy θx W θy θx W θy θx 1 X Z QTF8 θx 2 8 θx W θy θx 3 θx 1 W θy Y W θy θx 2 W θy θx 3 4 7 θx W θy 6 θx 5 W θy W θy θx θy 6 θx W θy W θy θx 4 W θy 5 W Y TTF6 Fig.7.3-2 Nodal Configuration For QTF8 And TTF6 Elements 142 . 6.6.7.3-4 Sign Convention For Stress Resultant Output 143 .7.3-3 Perforated Thick Plate Illustrating Use Of The TTF6 Element MXY MY MX MXY SY MX MXY SX Y Z SX MY SY X MXY Fig.Plate Elements Fig. 6.3 (in-plane) and section 7. σ max . and combined to give LMK MN membrane 0 K bending 0 OPRa | | PQS a T membrane bending U = RR | | V SR | | W T membrane bending U | V | W The component elements are listed in table 7.4-1) U. σ min the maximum and minimum principal membrane stresses. A lumped mass matrix is evaluated using the procedures presented in section 2.4-1 Element Membrane Bending RPI4 TRP3 PMI4 TPM3 QF4 TF3 Table 7. θ X . N xy 144 . For further details of the element formulation see Section 7.TPM3.TF3).PMI4) on the isoflex thin plate elements (QF4.3. N x . θ Y at each node The strain-displacement relationship. The membrane and bending stiffnesses are formed independently. β the angle between the maximum principal membrane stress and the local x-axis. resultant modulus matrix and thermal strains are defined in section 7.6 and [Z1.6. The final nodal variables are (fig. N y .6.4-1 Component Elements Used To Form Ribbed Plate Elements The element is formulated in a local Cartesian basis and then transformed to the global Cartesian system. V.4 Ribbed Plate (RPI4. Section 7.L1]. σ y . W.Element Formulations 7. Stress Resultant Output the membrane stress resultants/unit width in the local Cartesian system.6 (bending).7. TRP3) Formulation The 2-D flat ribbed plate elements are formulated by superimposing standard isoparametric plane stress elements (QPM4. Evaluation of stresses/strains The element output obtained at the element nodes consists of σ x .7. σ xy Stress Output direct and shear stresses in the local Cartesian system.6. 6. The local x-axis is defined as being a line joining the first and second element nodes.7.6. ψ min the maximum and minimum principal bending strains. M xy the moments/unit width in the local Cartesian system.7. The stress resultants are most easily interpreted if the local Cartesian axes are all parallel. ψ x .7. W W 4 U θx θy V θx U 3 θy V U 3 θx W θy V W Z U θx 1 θy V U RPI4 W θy θx 2 V U 1 W θy V U TRP3 W θy θx 2 V θx Y X Fig.6.Plate Elements M x . M y . ψ xy the flexural strains in the local Cartesian system. γ xy The sign convention for stress resultant and strain output is shown in fig. ψ max . The local y and z-axes are defined by a right hand screw rule (fig.4-4). The elements cannot be used for linear buckling analysis. Also. ψ y .1. ∈min the maximum and minimum principal membrane strains.4-1 Nodal Configuration For Ribbed Plate Elements 145 . Nonlinear formulation The element has no nonlinear capability but may be utilised in a nonlinear environment. ∈max .4-3. The xy-plane is defined by the third element node and the local x-axis. the presence of eccentricity requires that the forces and moments are examined at the mid-points of the element sides by averaging the nodal values. β the angle between the maximum principal bending strain and the local x-axis ∈x . ∈y . The nodal stress resultants are evaluated using the extrapolation procedures detailed in section 6. β the angle between the maximum principal membrane strain and the local x-axis. Strain Output the membrane strains in the local Cartesian system. 4-3 Sign Convention For Stress And Stress Resultant Output 146 .4-2 Ribbed Plate Illustrating Use Of RPI4 Element MXY MY Y σY σ X.6.6.Element Formulations RPI4 elements Y Z Y X X BRP2 elements Finite Element Mesh Fig.7. σ Y +ve tension σ XY +ve into XY quadrant MXY σ XY σX σ XY σY X σX Y Z MY MX MX MXY X MXY Stresses Stress Resultants Fig.7. 6.4-4 Local Cartesian System For Ribbed Plate Elements 147 .Plate Elements 4 3 1 x z y z y 1 Z Y x 2 2 3 X Fig.7. the initial degrees of freedom are (fig.7 Shell Elements 7.7. curved.7. at the mid-length node and Δu where Δu is the local axial relative (departure from linearity) displacement.1-1) U. non-conforming axisymmetric shell element formulated using the constraint technique.1-1) U. V.1 Axisymmetric Thin Shell (BXS3) Formulation The BXS3 element is a thin.7. Therefore. The final degrees of freedom for the element are (fig.7.Element Formulations 7. θ and Δu. Δv. The Kirchhoff condition of zero shear strain is applied at the two integration points by forcing ∂v ∂u ∂v + = − θz = 0 ∂x ∂z ∂x and eliminating the local transverse translational and rotational degrees of freedom at the central node. θz at the end nodes. V.7. The infinitesimal strain-displacement relationship is defined in the local cartesian system as ∈x = ∂u ∂x ∈z = U V cos φ − sin φ R R ∂2 v ∂x 2 ψx = − ψz = 1 ∂v cos φ R ∂x 148 . Δθ at the end nodes at the mid-length node. The global displacements and rotations are initially quadratic and are interpolated independently using linear Lagrangian shape functions for the end nodes and a hierarchical quadratic function for the central node. 7. The consistent and lumped mass matrices are evaluated using the procedures defined in section 2.Z1].C1. 149 .1-2) The elastic modulus and resultant modulus (or rigidity) matrices are defined as Explicit $ D= $ LMD MN 0 m 0 $ D b OP PQ $ D = b where Isotropic $ D = m Et 1 υ 1 − υ2 υ 1 t 1 − υ2 xz x Orthotropic Numerically Integrated E D= t 1 − υ2 $ D = m LM OP N Q LM E υ OP Nυ E Q xz z xz $ D b LM1 υOP 12e1 − υ j Nυ 1 Q LM E υ OP t = E Q 12e1 − υ j N υ Et 3 2 3 x xz z 2 xz xz z LM 1 MM υ MNυyy υ 1 υy y y υy y2 υy 2 OP P dy υy P P y Q υy y 2 2 The thermal strain vector is defined as Isotropic Orthotropic LM ΔTα OP MM daΔΔTTαf PP eψ j = M dy αP MM daΔTf αPP MN dy PQ LM ΔTα OP MM daΔTα PP Δ eψ j = M dyTf α P MM daΔTf PP α NM dy QP 0 t x z 0 t x z Further information on the element formulation is given in [S1.7.7.Shell Elements where R and φ are the radius and angle between the local and global Cartesian systems (fig. Element Formulations Evaluation of stresses/strains Element output is available at both the nodes and Gauss points and consists of (fig. ψ x . Linear and nonlinear buckling analysis. Geometrically and materially nonlinear analysis utilising the nonlinear material laws specified in 1.the meridional and circumferential membrane strains. Mz Strains ∈x . ∈z .7. Nonlinear formulation The axisymmetric shell element may be employed in Materially nonlinear analysis utilising the elastoplastic constitutive laws [O1] (section 4.1-5) Nx . the meridional and circumferential moments/unit width in the local Cartesian system.the meridional and circumferential bending strains. Mx . Nonlinear dynamics utilising the nonlinear material laws specified in 1. The forces and strains are output in the local Cartesian system. Note Layer stress output is also available when the nonlinear continuum plasticity model is utilised. Notes The BXS3 element may be used in conjunction with the stress resultant plasticity model (section 4. The top fibre lies on the +ve local y side of the element and +ve values define tension.7. Geometric nonlinearity utilises either 150 .7. defined as having its x-axis lying along the element axis in the direction in which the element nodes are specified. ψ z . Nz Stress Resultants the meridional and circumferential forces/unit width in the local Cartesian system.7. Geometrically nonlinear analysis.1-4). The forces have greatest accuracy at the Gauss points.2). such that the y-axis lies in the global XY plane. and the z-axis is parallel to the global Z-axis (up out of page) (fig. The local y and z-axes form a right-hand set with the x-axis.2). The nonlinear strain-displacement relationship is defined by ∈x = ∈z = ∂u 1 ∂u + ∂x 2 ∂x LM OP N Q 2 + 1 ∂v 2 ∂x LM OP N Q 2 u v u2 v2 uv cos φ − sin φ + cos2 φ + sin 2 φ − sin 2φ 2 2 R R 2R 2R 2R2 ∂ 2 v ∂u ∂ 2 v ∂v ∂ 2 u − + ∂x 2 ∂x ∂x 2 ∂x ∂x 2 ψx = − ψz = 1 ∂v u ∂v v ∂v cos φ − 2 cos2 φ + 2 cos φ sin φ R ∂x ∂x R R ∂x where R is the radius and φ is the angle between the local and global Cartesian systems.5). The output now approximates to the true Cauchy stresses and logarithmic strains. referred to the undeformed configuration. provided that the rotations are small within a load increment. The initial assumptions used in deriving the BXS3 element limit the rotations to one radian in a Total Lagrangian analysis and rotation increments of one radian in an Updated Lagrangian analysis (section 3. The loading approximates to being non-conservative.Shell Elements A Total Lagrangian formulation which accounts for large displacements but small strains. The loading is conservative. 151 . The forces and strains output with the geometrically nonlinear analysis will be the 2nd Piola-Kirchhoff forces and Green-Lagrange strains respectively. or An Updated Lagrangian formulation which takes account of large displacements and large rotations but small strains. Element Formulations V θZ U ΔV Δθ Z 2 3 ΔU ΔU 3 θZ U V 2 Y θZ V Y θZ U V 1 1 U X Initial Variables X Final Variables Fig.7.7.1-1 Nodal Configuration For The BXSs3 Element Axis of Revolution v, y φ u, x R Fig.7.7.1-2 Definition Of R And Φ For The Axisymmetric Shell Element 152 Shell Elements A A Plan Section A - A Finite Element Mesh (a) Spherical Shell Problem Definition Problem Definition (b) Circular Shell Finite Element Mesh Fig.7.7.1-3 Examples Illustrating The Use Of BXS3 Element 153 Element Formulations y x y x Y y x X Fig.7.7.1-4 Definition Of Local Cartesian System For BXS3 Element y z x Y z x y Mx Nz Mz Nx X Fig.7.7.1-5 Sign Convention For Stress Resultant Output 154 Shell Elements 7.7.2 Flat Thin Shell (QSI4, TS3) Formulation These flat shell elements are formulated in a local Cartesian system by superimposing standard isoparametric plane stress elements (QPM4,TPM3,PMI4) and the isoflex thin plate elements (QF4,TF3). The xy-plane of the local Cartesian system is evaluated using a least squares fit through the element nodes. The membrane and bending stiffnesses are then formed independently, and combined to give LMK MM MN membrane 0 K bending 0 0 K art 0 0 0 OPRa PP|a S PQ|θ T membrane bending z U RR |=| R V S | | M W T membrane bending z U | V | W where the component elements are listed in Table 7.7.2-1 Element Membrane Bending QSI4 TS3 PMI4 TPM3 QF4 TF3 Table 7.7.2-1 Primary Elements Used To Form Flat Thin Shell Elements Initially, the membrane stiffness is formed in terms of u and v, the in-plane displacements. An artificial in-plane rotational stiffness K is then added to prevent art singularities occurring when elements are co-planar. K art is defined as Triangles K art 1.0 −0.5 −0.5 1.0 −0.5 = k ip E x + E y tA −0.5 1.0 −0.5 −0.5 Quadrilaterals K art = k ip L i MM MN LM 1.0 dE + E itA MM−1 // 3 MN−1 / 3 −1 3 d x y OP PP Q OP PP PQ −1 / 3 −1 / 3 −1 / 3 1.0 −1 / 3 − 1 / 3 −1 / 3 1.0 −1 / 3 −1 / 3 −1 / 3 1.0 The in-plane stiffness parameter k ip has a default value of 0.02 which may be changed by using the SYSTEM command (variable STFINP). Once the local element matrices have been evaluated they are transformed to the global Cartesian system. The final nodal variables are (fig.7.7.2-1) U, V, θ x , θ y , θz at each node 155 Element Formulations The strain-displacement relationship is defined in section 7.3 (in-plane) and section 7.6 (bending). Note. The incompatible terms in the strain-displacement matrix are not used to evaluate nodal loads due to initial Gauss point stresses, e.g. thermal loading, initial stresses. For further details of the element formulation see section 7.3, section 7.6, [Z1,L1] A lumped mass matrix is evaluated using the procedures presented in section 2.7. Evaluation of stresses/strains The element output obtained at the element nodes consists of Stress Output σ x , σ y , σ xy direct and shear stresses in the local Cartesian system, σ max , σ min the maximum and minimum principal membrane stresses, β the angle between the maximum principal membrane stress and the local x-axis. Stress Resultant Output N x , N y , N xy the membrane stress resultants/unit width in the local Cartesian system, M x , M y , M xy the moments/unit width in the local Cartesian system. Strain Output the membrane strains in the local Cartesian system, ∈max , ∈min the maximum and minimum principal membrane strains, β the angle between the maximum principal membrane strain and the local x-axis, ψ x , ψ y , ψ xy the flexural strains in the local Cartesian system, ψ max , ψ min the maximum and minimum principal bending strains, β the angle between the maximum principal bending strain and the local x-axis The sign convention for stress resultant and strain output is shown in fig.7.7.2-3. The xy-plane of the local Cartesian system is evaluated using a least squares fit through the element nodes. The local x-axis is defined as being a line joining the first and second element nodes, and the local y and z-axes are defined by a right hand screw rule (fig.7.7.2-4) ∈x , ∈y , ∈xy 156 7. Average nodal stresses are in the global Cartesian system.2-2 Cylindrical Roof Example Illustrating Use Of Thin Flat Shell Elements 157 . The elements cannot be used for linear buckling analysis. The nodal stress is computed at each node directly.2-1 Nodal Configuration For Flat Thin Shell Elements Problem Description Finite Element Model Fig. Nonlinear formulation The elements have no nonlinear capability but may be utilised in a nonlinear environment. W W θz θy V W θx U 3 θz θy V θz θx U 3 V θy θx U 4 W θz Z U θx 1 θy V U QS4/QSI4 W θz θy V θz 1 U W V θy TS3 U W θz θy V θx 2 θx θx 2 Y X Fig.Shell Elements The nodal stress resultants are evaluated by extrapolating the strain-displacement relationship at the Gauss point to the nodes. The stress resultants are most easily interpreted if the local Cartesian axes are all parallel.7.7.7. 7.7.2-3 Sign Convention For Stress And Stress Resultant Output 4 3 1 x z y z y 1 Z Y x 2 2 3 X Fig.2-4 Local Cartesian System For Thin Flat Shell Elements 158 . σ Y +ve tension σ XY +ve into XY quadrant σY σ XY σX σ XY σY X σX Y Stresses Fig.7.7.Element Formulations MXY MY MX Stress Resultants MXY MX MXY Y Z MY X MXY σ X. Shell Elements 7.7.3 Flat Thin Shell Box (SHI4) Formulation The flat shell box element is formulated in a local Cartesian system by superimposing a non-standard isoparametric plane membrane element on the isoflex thin plate element. The xy-plane of the local Cartesian system is evaluated using a least squares fit through the element nodes. The membrane and bending stiffnesses are then formed independently and combined to give the total element stiffness in the local Cartesian system, i.e. LMK NM membrane 0 K bending 0 OPRa | S T QP| a membrane bending U = RR | | V SR | | W T membrane bending U | V | W The component bending stiffness and force vector used for this element is from the QF4 element [Section 7.6.1]. The elements use a non-standard plane membrane formulation which is more effective for modelling the in-plane bending in the web of box structures than the standard plane membrane formulation. The initial nodal configuration (fig.7.7.3-1) has 4 nodes with 3 in-plane degrees of freedom at each node u, v and ∂v / ∂x ξ where ∂v / ∂x ξ is the rotation of a line η = constant at each node and approximates to θz . In addition, incompatible displacement modes are utilised so that typically U= ∑ N i bξ, ηg Ui + ∑ Pi bξ, ηg a i i =1 i=1 n 2 where P1 ξ, η = 1 − ξ 2 and P2 ξ, η = 1 − ξ 2 b g b g and a i are nodeless degrees of freedom which are condensed out before element assembly. The extra incompatible modes are condensed out and the element matrices are then transformed to the global Cartesian system. This provides an element with the following nodal degrees of freedom (fig.7.7.3-2) U, V, W, θ x , θ y , θ z - at the corner nodes Δu - the relative (departure from linearity) local x-displacement for the mid-side nodes 159 Element Formulations The strain-displacement relationship is defined in section 7.3 (in-plane) and section 7.6 (bending). Note. No artificial in-plane rotational stiffnesses are required for this element. For further details of the element formulation see section 7.6, [L1,T2]. A lumped mass matrix is evaluated using the procedures presented in section 2.7. Evaluation of stresses/strains The element output obtained at the element nodes consists of Stress Output σ x , σ y , σ xy direct and shear stresses in the local Cartesian system, σ max , σ min the maximum and minimum principal membrane stresses, β the angle between the maximum principal membrane stress and the local x-axis. Stress Resultant Output N x , N y , N xy the membrane stress resultants/unit width in the local Cartesian system, M x , M y , M xy the moments/unit width in the local Cartesian system. Strain Output the membrane strains in the local Cartesian system, ∈max , ∈min the maximum and minimum principal membrane strains, β the angle between the maximum principal membrane strain and the local x-axis, ψ x , ψ y , ψ xy the flexural strains in the local Cartesian system, ψ max , ψ min the maximum and minimum principal bending strains, β the angle between the maximum principal bending strain and the local x-axis. The sign convention for stress resultant and strain output is shown in fig.7.7.3-4. The xy-plane of the local Cartesian system is evaluated using a least squares fit through the element nodes. The local x-axis is defined as being a line joining the first and second element nodes, the local y and z-axes are defined by a right hand screw rule (fig.7.7.3-5). The nodal stress resultants are evaluated by extrapolating the strain displacement relationship at the Gauss points to the nodes, and then computing the nodal stress at each node directly. The stress resultants are most easily interpreted if the local Cartesian axes are all parallel. ∈x , ∈y , ∈xy 160 Shell Elements Note. The averaged nodal stresses are output in the global Cartesian system. Nonlinear formulation The elements have no nonlinear capability but may be utilised in a nonlinear environment. The elements cannot be used for linear buckling analysis. V θz 4 U θz 3 V U V V Z θz 1 SHI4 U θz 2 U Y Fig.7.7.3-1 Initial In-Plane Nodal Configuration W θx U W θz V 4 θy θx U 3 θz V θy Z Y W W θz V U θx 1 θy SHI4 U θx 2 θz θy V X Fig.7.7.3-2 Final Nodal Configuration For Flat Thin Box Shell ElemENTS 161 Element Formulations Box Girder Box Girder Bridge Fig.7.7.3-3 Structures Suitable For Analysis With Flat Box Shell Elements MXY MY MX Stress Resultants MXY MX MXY Y Z MY X MXY σ X, σ Y +ve tension σ XY +ve into XY quadrant σY σ XY σX σ XY σY X σX Y Stresses Fig.7.7.3-4 Sign Convention For Stress And Stress Resultant Output 162 Shell Elements 4 3 z y 1 Z Y x 2 X Fig.7.7.3-5 Local Cartesian System For Thin Flat Box Shell Elements 7.7.4 Semiloof Thin Shell (QSL8, TSL6) Formulation The Semiloof shell element is a thin, doubly curved, isoparametric element formed by applying Kirchhoff constraints to a three dimensional degenerated thick shell element. The displacements and rotations are considered independent and the unconstrained nodal configurations are (fig.7.7.4-1) U, V, W θx , θy - at the corner and mid-side nodes, - at the loof nodes, - at the central node, and w, θ x , θ y where θ x and θ y are the rotations of the through-thickness normals. These rotations include transverse shear deformations. An element with thin shell performance is then produced by constraining the shear strains to zero at discrete points within the element, i.e. by ensuring that [I1] γt = ∂w − θy = 0 ∂x at the points shown in fig.7.7.4-2. Where γ t is the through-thickness shear strain tangential to the element edges. This provides 6 and 8 constraints respectively for the triangles and quadrilaterals which are suitable for eliminating the tangential rotations at the loof nodes. z A γ xz dA = 0, z A γ yz dA = 0 163 These constraints provide extra equations that permit certain nodal degrees of freedom to be discarded. Using the assumptions of thin shell theory. for Isotropic materials 164 . the strain-displacement relationship is written as ∈x = ∂u ∂x ∈y = ∂v ∂y ∂u ∂v + ∂y ∂x ∂2w ∂x 2 γ xy = ψx ≈ − ψy ≈ − ∂2 w ∂y 2 ∂2 w ∂x∂y ψ xy ≈ −2 The isotropic and orthotropic modulus and resultant modulus (rigidity) matrices are defined as Explicit $ D= $ LMD MN membrane 0 $ D bending 0 OP PQ where. V. .at the corner and mid-side nodes.at the loof nodes. This provides 1 constraint suitable for removing the central translation of the quadrilaterals. z S γ n dS = 0 where γ n is the transverse shear strain normal to the element sides and the integral is performed using 2-point quadrature along each side. This provides 2 constraints for both the triangles and quadrilaterals which are suitable for removing the rotations at the central node.Element Formulations where the integral is performed using 2*2 Gauss quadrature.7.7.4-3) U. W and θ . The final nodal configurations are (fig. Notes To obtain a valid material υxy < E x / E y d i 1/ 2 remove the γ YZ and γ XZ shear strain rows and columns.1 and is reduced to the plane stress modulus matrix in the following way: 165 .Shell Elements $ D membrane $ D bending LM1 υ 0 OP E M υ 1 0 P = 1− υ M a1 − υf PP MN0 0 2 Q LM1 υ 0 OP Et Mυ 1 0 PP = 12e1 − υ j M MN0 0 a1 − υf PQ 2 2 3 2 and for Orthotropic materials $ D membrane LM 1 / E = M− υ / E MN 0 x xy 3 x xy − υ xy / E x x 0 0 1 / G xy 0 1 / Ey 0 − υ xy / E x x OP PP Q −1 $ D bending L 1/ E t M = M− υ / E 12 M N 0 1 / Ey 0 0 1 / G xy OP PP Q −1 where υyx has been set to υxy E y / E x to maintain symmetry. invert the matrix so that the stress-strain relationship is obtained. remove the s Z row and column since this stress is assumed to be zero. This 6 by 6 modulus matrix is the same as that given in section 7.4. re-invert the matrix to obtain the stress-strain relationship (a 3 by 3 matrix). A three dimensional orthotropic modulus matrix may be specified by using the appropriate data input. the maximum and minimum principal membrane stresses. σ xy σ max . σ min β direct and shear stresses in the local Cartesian system. σ y .Element Formulations Numerically Integrated E D= t 1 − υ2 z LM 1 MM υ MM 0 MMυzz MM 0 N υ 1 0 υz z 0 a f a f 0 0 1− υ 2 0 0 1− υ z 2 z υz 0 z υz 2 0 2 υz z 0 υz 2 z2 0 OP a f PP PP PP a f PP Q 0 0 1− υ z 2 dz 0 0 1 − υ z2 2 The thermal strain vector is defined as αΔT OP LM αΔT MM PP 0 Δ eψ j = MM dadzTf LMNα + ΔT dα OPQPP dT P MM daΔTf L dα O P MM dz MNα + ΔT dT PQPP 0 Q N α ΔT OP LM α ΔT PP MM α ΔT MM daΔTf L dα P MNα + ΔT dT OPQ PP eψ j = MM dz α Δ MM dadzTf LMα + ΔT ddT OP PPP N MM daΔTf Lα + ΔT dα QOPP PP dT Q Q NM dz MN 0 t x y xy x y x y 0 t xy xy Isotropic Orthotropic Full details of the element formulation are given in [I1]. Both consistent and lumped mass matrices are available and are evaluated using the procedures defined in (section 2. the angle between the maximum principal membrane stress and the local x-axis. 166 .7). Evaluation of stresses/strains The element output obtained at the element nodes and Gauss points consists of Stress Output σ x . The nodal stress resultants are evaluated using the extrapolation procedures detailed in section 6. coincides with the curvilinear line ξ = constant (fig. with the +ve direction defined by the +ve ξ direction. ψ xy ψ max .1.7.4-8). Strain Output ∈x . ∈y . The sign convention for stress resultant and strain output is shown in fig. Geometrically nonlinear analysis.7. The average nodal stresses are in the global Cartesian system. the flexural strains in the local Cartesian system. The x-axis is then formed perpendicular to the y-axis and tangential to the shell mid-surface. the angle between the maximum principal membrane strain and the local x-axis. The local Cartesian system varies over the element for curved elements. 167 . N xy the membrane stress resultants/unit width in the local Cartesian system. The local z-axis forms a righthanded set with the x and y-axes.7.3). M x . the angle between the maximum principal bending strain and the local x-axis. ψ min β the membrane strains in the local Cartesian system.7. For the triangular element.4-9. For the quadrilateral element. The local x-axis is perpendicular to the local y-axis in the +ve ξ direction and is tangential to the shell mid-surface. the maximum and minimum principal membrane strains. ∈min β ψ x . the local Cartesian system is formed by orientating the local y-axis parallel to a line joining the mid-point of the first side with the 5th node. Geometrically and materially nonlinear analysis utilising the nonlinear material laws specified in 1. Nonlinear formulation The Semiloof shell element may be employed in Materially nonlinear analysis utilising the elastoplastic constitutive laws [O1] (section 4. ψ y . ∈xy ∈max .Shell Elements Stress Resultant Output N x . N y . Notes The Gauss point stresses are converted to the global Cartesian system before extrapolation.2) and the nonlinear concrete model (section 4. M y . the maximum and minimum principal bending strains. at any point within the element. The +ve z-axis defines the top surface. the local y-axis. M xy the moments/unit width in the local Cartesian system. which takes account of large displacements and moderately large strains provided that the strain increments are small. The loading approximates to being non-conservative. The nonlinear strain-displacement relationship is defined by LM OP N Q ∂v 1 L ∂u O ∈= + M P ∂y 2 N ∂y Q ∈x = ∂u 1 ∂u + ∂x 2 ∂x y 2 2 LM OP N Q 1 L ∂v O + M P 2 N ∂y Q + 1 ∂v 2 ∂x 2 2 LM OP N Q 1 L ∂w O + M P 2 N ∂y Q + 1 ∂w 2 ∂x 2 2 γ xy = ∂u ∂v ∂u ∂u ∂v ∂v ∂w ∂w + + + + ∂y ∂x ∂x ∂y ∂x ∂y ∂x ∂y ∂2w ∂x 2 ψx = − ψx = − ∂2w ∂y 2 ∂2 w ∂x∂y 2 ψ xy = −2 The output is now in terms of the 2nd Piola-Kirchhoff stresses and GreenLagrange strains referred to the undeformed configuration. Notes Geometric nonlinearity may be represented with either A Total Lagrangian formulation which accounts for large displacements but small rotations and strains. The initial assumptions used in deriving the shell elements limit the rotations to one radian in a Total Lagrangian analysis. or An Updated Lagrangian formulation. Linear and nonlinear buckling analysis. The loading is conservative. The output is now in terms of the True Cauchy stresses and the strains approximate to logarithmic strains. and rotation increments of one radian in an Updated Lagrangian analysis (section 3.Element Formulations Nonlinear dynamics utilising the nonlinear material laws specified in 1.5). 168 . Shell Elements V 7 W θx θy W V 8 W θx θy V Y U W 1 X Z θx θy θy V U W 4 U θx W θx W 2 θy V U θx W (a) QSL8 5 V U 3 θy θy θx V U U θx W 4 θy V U U θy θx 6 W V U θy θy θx θx W 5 V U V 6 θx θy V V Y W 1 X Z U θ x θy W 2 U θx W θy θy V θx W (b) TSL6 3 U Fig.7.7.4-1 Initial Nodal Configurations For QSL8 And TSL6 Elements 169 . Element Formulations 1/ 3 1/ 3 1/ 3 2 1/ 3 2 (a) QSL8 1/ 3 1/ 3 2 2 1/ 3 1/ 3 1/ 3 1/ 3 2 (b) TSL6 Fig.4-2 Locations Where Transverse Shear Strains Tangential To The Element Edge Are Constrained To Zero 170 .7.7. Shell Elements V 7 W V U θ2 θ1 W V 8 W U θx W 4 θy V U W 6 U θ1 θ2 W 5 V U θ2 V Y U W 1 X Z (a) QSL8 W θ1 V U 2 θ2 θ1 V U W V 5 W θ1 V θ2 U 3 V 6 W θ2 V V Y U W 1 X Z (b) TSL6 θ1 W 2 U U W 4 U θ1 V θ2 W U 3 Fig.4-3 Final Nodal Configuration For QSL8 And TSL6 Elements 171 .7.7. 7.7.4-4 Tubular Joint Example Illustrating Use Of QSL8 And TSL6 Elements Fig.Element Formulations Fig.7.7.4-5 Pressure Vessel Example Illustrating Coupling Of HX20 And QSL8 Elements 172 . 7.7.4-7 Bending Mechanism For QSL8 Element 173 .4-6 Stiffened Shell Illustrating Coupling Between QSL8 And BSL3 Elements Fig.Shell Elements QSL8 elements BSL3 elements Problem Definition Finite Element Mesh Fig.7.7. 7.7.4-8 Local Cartesian System 174 .Element Formulations 5 z η 6 y 3 x 4 ξ 7 2 8 1 (a) QSL8 Element z y 5 4 x 3 6 2 1 (b) TSL6 Element Fig. 7.7.4-9 Sign Convention For Stress And Stress Resultant Output 175 . Ν Y +ve tension Ν XY +ve into XY quadrant Y ΝY Ν XY ΝX Ν XY ΝY X ΝX MXY MY MX MXY MX MXY Y Z MY X MXY Stress Resultants Fig.Shell Elements σ X. σ Y +ve tension σ XY +ve into XY quadrant Y σY σ XY σX σ XY σY X σX Stresses Ν X. however.7. TTS6) are formulated using a standard isoparametric approach. QTS4. θ x . This definition of the rotations is used when a smooth surface configuration is to be modelled (fig.5-3. θz (fig. For the four noded quadrilateral (QTS4) the factors for interpolating from the sampling points to the gauss points are R1 = R2 = 1 1− η 2 a f a f OPa f Q 1 1+ η 2 while for the eight noded element (QTS8) the factors are R1 = ξ 1 1 1− 1 − η − R5 a 4 4 LM N 176 .5 Thick Shells (TTS3.7. These rotations include transverse shear deformations and relate to a set of 'local' axes set up at each node.5-1).QTS8) adopt an assumed strain field for interpolation of the transverse shear strains. The triangular elements (TTS3.7. The displacements and rotations are considered independent and the nodal degrees of freedom are (fig. the direction of these axes is dictated by the direction of the nodal normal. The quadrilateral elements (QTS4.7.7.5-1) U. The location of the transverse shear sampling points for defining the assumed strain fields are shown in fig. TTS6. the axis chosen corresponds with the smallest component of the nodal vector.5-1).7.Element Formulations 7. θα and θβ are the rotations of the through-thickness normals.at all nodes. The cross product of this axis and the nodal vector defines the second axis of rotation for θβ (fig. W. θα . or a branched shell junction. V.7.7. The director is subsequently referred to as the normal.5-2). QTS8) Formulation The formulation for this family of thick shell elements is based on the degeneration of a three dimensional continuum. The normal is considered to remain straight during deformation for computation of displacements through the element thickness. In this approach.7. θβ . θ y . these rotations are transformed to relate to global axes. the displacements at any point in the shell are defined by the translation of the reference surface together with the rotation of a director. The inclusion of an assumed strain field prevents the element from 'shear locking' when used as a thin shell. In the event of a discontinuity. To avoid singularities.7. the director need not be initially normal to the reference surface. One of the global axes is chosen to define the θα rotation.7. connection with a beam element. Approximation of the Kirchhoff-Love thin shell hypothesis. ∈i are the transverse shears at the sampling points. The local axes are set up using $ e1 = G ξ / G ξ $ $ $ e3 = e1 x G η / e1 x G η e j 177 . It is necessary to express the transverse shear strains in terms of covariant components so that interpolation can be carried out using the isoparametric map.Shell Elements R2 = R3 R4 LM OPa f N Q 1 L ξO 1 = M1 + Pa1 − ηf − R 4 N aQ 4 1 L ξO 1 = M1 + Pa1 + ηf − R 4 N aQ 4 ξ 1 1− a 4 ξ 1 1 1− 1 + η − R5 a 4 4 5 5 R5 = LM L O OP 1 − η MN MN PQ PQe j 2 2 where a =1/ 3 and Si ( η. ξ) = R i (ξ. η) The covariant transverse shear strains at the gauss points are then given by ∈ξζ = i ∑ R i bξ. ηg ∈iηζ i =1 where ∈ξζ and ∈ηζ are the covariant transverse shear strains at the gauss points and i ∈ξζ . No spurious zero energy modes using full numerical integration. The stress and strain terms are ultimately transformed to relate to a local orthogonal set of axes at each gauss point. ηg ∈ξζ i =1 n n ∈ηζ = ∑ Si bξ. ηζ Using this representation of shear strains allows Correct representation of the six rigid body modes. Element Formulations $ $ $ e 2 = e3 x e1 where G ξ and G η are the covariant base vectors at a gauss point. Strains in the curvilinear system ∈lm may then be transformed to strains in the $ orthogonal local system ∈ij by using the contravariant base vectors $ d∈ i =∈ ij lm d i dG ⋅ e$ ieG ⋅ e$ j 1 i m j b g The elements are formulated using the plane stress hypothesis so that σ zz in the thickness direction is set to zero. The continuum strains are evaluated at integration points through the thickness, and for the geometrically linear case these strains are given by $ ∈xx = ∂u ∂x $ ∈yy = ∂v ∂y ∂u ∂v $ γ xy = + ∂y ∂x ∂v ∂w $ γ yz = + ∂z ∂y ∂u ∂w $ γ yz = + ∂z ∂x Material properties are specified in the local orthogonal axes. For a thick shell the modulus matrix is condensed so that the plane stress hypothesis is observed. The isotropic modulus matrix is given by [Z1] D= E 1 − υ2 LM 1 MM MM MM MNSymm. υ 1 0 0 1− υ 2 0 0 0 1− υ 2.4 OP P 0 P PP 0 P 1− υP P 2.4 Q 0 0 If orthotropic properties are specified the modulus matrix becomes 178 Shell Elements E υ LM E MM d1 − υ υ i d1 − υ υ i E MM d1 − υ υ i D=M MM MM MN Symm. x x yx xy yx xy yx y xy yx 0 0 G xy 0 0 0 G yz 1.2 OP PP 0 P P 0 P PP 0 P G P P 1.2 Q 0 xz Factors of 5/6 have been included in the transverse shear terms to take account of a parabolic distribution through the thickness. As the material properties are specified in local element directions and the element formulation is based on covariant components of strain, the modulus matrix must be transformed. The required transformation of the modulus matrix is C ijkl = G i ⋅ ea G j ⋅ e b G k ⋅ e c G l ⋅ e d D d id id id i abcd where G m m = ξ, η, ζ are the contravariant base vectors. Full details of the element formulations may be found in [D4],[H9] and [S7]. Both consistent and lumped mass matrices are available and are evaluated using the procedures defined in (section 2.7). Evaluation of stresses/strains The element output obtained at the element nodes and Gauss points consists of Stress Output σ x , σ y , σ xy σ yz , σ xz , σ e direct and shear stresses in the local Cartesian system, together with von Mises equivalent stress Three dimensional principal stresses and the corresponding direction cosines may also be output Stress Resultant Output N x , N y , N xy the membrane stress resultants/unit width in the local Cartesian system, M x , M y , M xy the moments/unit width in the local Cartesian system, Sx ,Sy the shear stress resultants/unit width in the local Cartesian system 179 Element Formulations Strain Output ∈x , ∈y , ∈xy , ∈yz , ∈xz , ∈e the direct and shear strains in the local Cartesian system, together with von Mises equivalent strain The local cartesian systems are set up at the element reference surface. For curved elements, the local Cartesian system will vary over the reference surface. The local xaxis, at any point within the element, coincides with the curvilinear line η = constant in the direction of increasing ξ (fig.7.7.5-4). The direction of the local z-axis is defined by the vector product of the local x-axis and the curvilinear line ξ = constant (in the direction of increasing η). The local y-axis is defined by the vector product of the local z and local x-axes. The +ve z-axis defines the element top surface. The position of the origin of the curvilinear system for each element together with the directions of increasing values are shown in (fig.7.7.5-5). The sign convention for stress and strain output is shown in fig.7.7.5-6 and fig.7.7.57. The nodal stresses and strains are evaluated using the extrapolation procedures detailed in section 6.1. Notes The Gauss point stresses are converted to the global Cartesian system before extrapolation. The average nodal stresses are in the global Cartesian system. Nonlinear formulation The thick shell elements may be employed in Materially nonlinear analysis utilising the elastoplastic constitutive laws [O1] (section 4.2) and the nonlinear concrete model (section 4.3). Geometrically nonlinear analysis using a Total Lagrangian formulation. Geometrically and materially nonlinear analysis utilising the nonlinear material laws specified in 1. Nonlinear dynamics utilising the nonlinear material laws specified in 1. Linear and nonlinear buckling analysis. Creep analysis Note. The Total Lagrangian formulation used for these elements is valid for both large displacements and large rotations. However, the formulation is only valid for small strains. The nonlinear strain-displacement relationship is defined by ∈x = ∂u 1 ∂u + ∂x 2 ∂x LM OP N Q 2 + 1 ∂v 2 ∂x LM OP N Q 2 + 1 ∂w 2 ∂x LM OP N Q 180 2 Shell Elements ∈y = ∂v 1 ∂u + ∂y 2 ∂y LM OP N Q 2 + 1 ∂v 2 ∂y LM OP N Q 2 + 1 ∂w 2 ∂y LM OP N Q 2 γ xy = γ yz = γ xz = ∂u ∂v ∂u ∂u ∂v ∂v ∂w ∂w + + + + ∂y ∂x ∂x ∂y ∂x ∂y ∂x ∂y ∂v ∂w ∂u ∂u ∂v ∂v ∂w ∂w + + + + ∂z ∂y ∂z ∂y ∂z ∂y ∂z ∂y ∂u ∂w ∂u ∂u ∂v ∂v ∂w ∂w + + + + ∂z ∂x ∂z ∂x ∂z ∂x ∂z ∂x The output is now in terms of the 2nd Piola-Kirchhoff stresses and Green-Lagrange strains referred to the undeformed configuration. The loading is conservative. 181 Element Formulations θβ V θα Z,w θβ (a) 5 degrees of freedom Definition of nodal rotations when global X defines θ α . X,u θα V Y,v Z,w θz θy Y,v θx X,u (b) 6 degrees of freedom Fig.7.7.5-1 Nodal Variables For Thick Shell Elements 182 7.7.5-2 Smooth And Discontinuous Surface 183 .Shell Elements Default angle < 20o Averaged nodal vector V2 V1 2 1 Element 1 Element 2 3 V3 (a) SMOOTH SURFACE (5 degrees of freedom) Default angle > 20o Separate nodal vector V 21 V 22 V1 2 V3 1 Element 1 Element 2 3 (b) DISCONTINUOUS SURFACE (5 degrees of freedom) Fig. 5-3 Shear Sampling Points 184 .7.Element Formulations 2 4 3 4 3 η 1 ξ 2 η ξ 1 2 1 1 Shear γ ξζ 2 Shear γ ηζ (a) QTS4 3 5 7 7 6 6 a a η 4 5 3 a 8 a 1 5 η ξ a a 4 8 5 ξ 4 2 a a = 3-1/2 a 1 2 Shear γ ηζ 3 1 1 2 Shear γ ξζ 3 2 (b) QTS8 Fig.7. 5-4 Local Cartesian Systems 185 .Shell Elements z y 5 4 3 ξ = constant x 2 6 η = constant 1 (a) TTS6 (TTS3 axes coincide when element is flat) 5 ξ = constant 6 y η = constant z x 4 7 3 8 2 1 (a) QTS8 (QTS4 axes coincide when element is flat) Fig.7.7. 5-5 Curvilinear Coordinates 186 .Element Formulations 2 3 4 ζ 3 ξ ζ 5 ξ η 2 η 6 1 (a) TTS3 3 6 1 (b) TTS6 5 4 η ζ ξ 7 η ζ ξ 4 2 8 2 3 1 (c) QTS4 1 (d) QTS8 Fig.7.7. 5-6 Sign Convention For Continuum Stress Output 187 .7.Shell Elements y σy σ xy σx σ xy σy x σx z z σ xz σ yz σ xz x σ yz y Direct stress (+ve) Tension Shear stress (+ve) Shear into XY.7. YZ and ZX quadrants Note: Positive values shown in figure Fig. 5-7 Sign Convention For Stress Resultant Output 188 .7.Element Formulations Y Y Mxy My Mx Mxy Sy Sx Mxy Mx Mxy Sy Sx My Y X Nx Nxy Ny X Nxy Nxy Ny Nxy Nx X Membrane stress Flexural stress Shear stress Note: (+ve) Direct tension (+ve) In-plane shear into XY quadrant (+ve) Hogging moment (producing +ve stresses on the element top surface) (+ve) In-plane shear into YZ and XZ quadrants Positive values shown in figure Fig.7. The field gradient is related to the flow by q x = k gx Element Output The element output consists of the gradients of the field variables and flow in the local element axis system at either the Gauss or nodal points.e. convection. radiation or applied heat flux. The Gauss point values are generally more accurate. The cross sectional area of the material is defined at each node and may vary over the element. The bar is assumed to be perfectly insulated along its length and may transfer heat across its end areas via conduction.8-1) has either 2 or 3 nodes and may transfer heat along its length.8 Field Elements 7.Field Elements 7. BFD3) The thermal bar element (fig.8. i.7. The nodal degree of freedom is the field variable Φ. 189 . The gradient-field variable relationship is defined as gx = ∂Φ ∂x where x represents the local x-direction of the element.1 Thermal Bar (BFD2. gx qx field gradient flow The local x-axis lies along the element axis in the direction in which the element nodes are specified. 7.7.8-1 Nodal Configuration For Bar Field Elements Struts represented with BFD2 elements Pressure Continuum elements Heat transfer between members Fig.U Fig.V 1 2 1 BFD2 2 BFD3 3 X.Element Formulations Y.8-2 Example Illustrating The Use Of Bar Field Elements 190 . Y.8-3 Nodal Configuration For Axisymmetric Bar Field Elements 191 . i. convection. radiation or applied heat flux.8.2 Thermal Axisymmetric Bar (BFX2.e. gx qx field gradient flow The local x-axis lies along the element axis in the direction in which the element nodes are specified. The Gauss point values are generally more accurate. The nodal degree of freedom is the field variable Φ. BFX3) The axisymmetric thermal bar element (fig. The gradient-field variable relationship is defined as gx = ∂Φ ∂x The field variable is related to the flow by q x = k gx Element output The element output consists of the gradients of the field variables and flows in the local element axis system at either the Gauss or nodal points.U Fig.8-3) has either 2 or 3 nodes and may transfer heat along its length.Field Elements 7.7. The bar is assumed to be perfectly insulated along its length and may transfer heat across its end areas via conduction.7.V 1 2 1 BFX2 2 BFX3 3 X. convection or radiation.8-5) has 2 nodes and may transfer heat between two nodal points by either conduction. The heat flow is positive in the direction of the local x-axis and is defined as Conduction Convection Radiation q = K(Φ1 − Φ 2 ) q = h c (Φ1 − Φ 2 ) 4 4 q = h r (Φ1 − Φ 2 ) where K is the gap conductance.7.3 Thermal Link (LFD2.8-4 Spinning Cylinder Example Illustrating The Use Of Axisymmetric Bar Field Elements 7.Element Formulations Problem Definition Finite Element Mesh Fig. The cross sectional area of the material is defined at each node and may vary over the element. LFS2. LFX2) The thermal link element (fig.8. The nodal degree of freedom is the field variable Φ. 192 . h c and h r are the convective and radiative heat transfer coefficients.7. Element output The element output consists of the gradients of the field variables and flows in the local element axis system at either the Gauss or nodal points. i. The element stiffness matrix may be derived from the nodal flows which are defined as Q1 = K t A Φ1 − Φ 2 = − Q 2 b g where K t is the combined heat transfer coefficient for the element.Field Elements Note that the material properties may be dependent on both the gap distance and the temperature evaluated at the centre of the element. The Gauss point values are generally more accurate. 193 .e. The stiffness matrix is then rewritten as K=K +K o o 1 where K is the linear contribution to the stiffness matrix defined as follows K = KtA o LM 1 −1OP N −1 1 Q t 1 2 and K1 is the nonlinear contribution to the stiffness matrix defined as K 1 LM ∂K AbΦ − Φ g ∂Φ =M MM− ∂K AbΦ − Φ g N ∂Φ 1 t 1 1 2 ∂K t A Φ1 − Φ 2 ∂Φ 2 ∂K t − A Φ1 − Φ 2 ∂Φ 2 b b g OP P gPPQ The nodal flows are evaluated directly using Q = KΦ where K is evaluated using the current properties in temperature dependent analyses. gx qx field gradient flow The local x-axis lies along the element axis in the direction in which the element nodes are specified. QFD8. The gradientfield variable relationship is defined as 194 .8-5 Nodal Configuration For Link Field Elements LFX2 elements QAX8 elements Problem Definition Finite Element Mesh Fig.U Fig.7.V 1 LFD2 LFX2 LFS2 2 X.8-7) is defined in the global XY-plane.8.4 Plane Field (QFD4. TFD6) The plane field element (fig.Element Formulations Y.7.7. TFD3.8-6 Thermal Analysis Of An Interface Fit Illustrating The Use Of Link Field Elements 7. 195 . qY LMk 0 OP N0 k Q Lk 0 OP k=M N0 k Q k= x y field gradient flow The Gauss point values are generally more accurate. gX . gY qX .e.Field Elements gX = gY = ∂Φ ∂X ∂Φ ∂Y The isotropic and orthotropic thermal conductivity modulus matrices are defined as follows Isotropic Orthotropic Element output The element output consists of the gradients of the field variables or flows in the global axis system at either the Gauss or nodal points. i. Element Formulations 3 6 7 5 4 8 2 1 Y QFD4 1 1 2 QFD8 5 4 6 2 3 TFD3 2 TFD6 3 4 3 1 X Fig.8-7 Nodal Configuration For Plane Field Elements 196 .7. 8-8 Cofferdam Example Illustrating The Use Of Plane Field Elements 197 .Field Elements (a) Problem Definition (b) Finite Element Mesh Fig.7. i. TXF3. QXF8.5 Axisymmetric Field (QXF4. TXF6) The axisymmetric field elements (fig. gY qX . The gradient-field variable relationship is defined as gX = gY = ∂Φ ∂X ∂Φ ∂Y The isotropic and orthotropic thermal conductivity modulus matrices are defined as follows Isotropic Orthotropic Element output The element output consists of the gradients of the field variables or flows in the global axis system at either the Gauss or nodal points.7. gX . 198 .8-9) are formulated using the axisymmetric quasi-harmonic equation (section 2.e.10) and are defined in the global XY-plane. qY LMk 0 OP N0 k Q Lk 0 OP k=M N0 k Q k= x y field gradient flow The Gauss point values are generally more accurate.Element Formulations 7.8. Field Elements 3 6 7 5 4 8 2 1 Y QXF4 1 1 2 QXF8 5 4 6 2 3 TXF3 2 TXF6 3 4 3 1 X Fig.7.8-9 Nodal Configuration For Axisymmetric Field Elements 199 . PF15.7. TF4. The nodal degree of freedom is the field variable Φ. PF6. PF12.Element Formulations Well Flow (a) Problem Definition Point sink to represent well (b) Finite Element Mesh Fig.7.8.10). TF10) The solid field elements (fig. The gradient-field variable relationship is defined as gX = gY = gZ = ∂Φ ∂X ∂Φ ∂Y ∂Φ ∂Z 200 .6 Solid Field (HF8. HF16.8-11) are formulated using the 3-D quasi-harmonic equation (section 2.8-10 Groundwater Flow Example Illustrating The Use Of Axisymmetric Field Elements 7. HF20. Field Elements The isotropic and orthotropic thermal conductivity modulus matrices are defined as follows Isotropic Orthotropic Element output LMk MM N0 LMk k=M0 MN 0 k= 0 0 0 0 ky 0 0 k 0 OP PP kQ OP 0P k P Q 0 z x The element output consists of the gradients of the field variables or flows in the global axis system at either the Gauss or nodal points. i. g Z qX. qZ field gradient flow The Gauss point values are generally more accurate.e. qY. g X . gY . 201 . Element Formulations 8 5 6 1 2 19 18 16 12 7 6 10 2 3 12 7 8 9 1 2 3 4 10 7 3 1 TF4 HF8 7 4 3 9 16 15 14 12 11 8 7 4 3 6 5 HF16 13 10 1 2 20 6 17 11 5 4 HF20 4 5 1 2 15 10 11 7 1 8 2 3 9 5 4 2 3 TF10 13 9 1 14 15 8 3 PF6 11 10 6 4 PF12 14 13 12 9 6 4 PF15 5 5 8 6 1 2 Fig.7.8-11 Nodal Configuration For Solid Field Elements 202 . Field Elements Beryllia heat sinks Silicon Chips Stainless steel side walls Copper Base Fig.8-12 Thermal Analysis Of A Hybrid Power Assembly Illustrating Use Of Solid Field Elements 203 .7. HX8L and HX16L structural elements. or a 3 corner point quadratic (for PF12C) gauss integration scheme outside the through thickness loop. For these elements the field variable at a node (temperature) is used to interpolate a temperature field that varies linearly over the thickness of the whole element and quadratically in-plane for the higher order elements. For the integration of the element thermal stiffness matrix. As with the PN6L. while the shape function derivative matrices are integrated using a plane 2x2 (for HF8C and HF16C). HF8C. are defined in terms of natural coordinates ξ and η. or a single point (for PF6C). in order to speed up the computation the elements are restricted to constant layer thicknesses [H13].7 Solid Composite Field (HF8C. The shape functions N i ( top) = N i ( bot ) = N i .Element Formulations 7. To overcome this difficulty layered brick elements were developed where several laminae are included in a single element. PN12L. HF16C. for the HF16C element these are given by: Ni = Ni = 1 1 + ξ i ξ 1 + ηi η ξ i ξ + ηi η − 1 4 1 1 − η2 ξ 2 − ξ 2 η 1 + ξ i ξ + ηi η i i 2 b gb gb g for corner nodes for mid-side nodes e jb g The PF12C. To form the complete shape functions for the brick element N br . 4 and 3 noded membranes. see figure 7. the thermal conductivity matrices are summed layerwise through the thickness. PF12C) If brick elements are used for a structural analysis of composites the number of degrees of freedom even for small laminate structures rapidly becomes very large leading to prohibitively excessive computer costs. PF6C. and PF6C elements use the equivalent shape functions for 6.8-13. The shape functions for the top and bottom surfaces of the composite elements can be considered to be single membrane element shape functions. Composite field elements are designed to compliment the structural elements in a thermo-mechanical analysis.8. linear interpolation is used between the functions for the top and bottom surfaces: NT = br 1 −ζ + 1 N T bot f . This indicates that the (in plane) shape function derivative matrix is independent of the non-dimensional coordinate through the thickness. This limitation requires the calculation of only a 2x2 Jacobian matrix. ζ + 1 N T top g ia ib 2 b g b g The in-plane and through-thickness shape functions can then be separated to give: N T = φ T + ζψ T br where 204 . . 205 ... LM ∂φ + ζ ∂ψ OP MM ∂ x ∂ x PP ∂ψ ∂φ B=M +ζ MM ∂ y 2 ∂ y PPP MM c ψ PP Q N T T T T T and ‘c’ is the overall depth of the element shown in figure 7....Field Elements φT = ψT = 1 T T N .variable relationship is given by: g = BΦ where B is the field gradient .N 2 i i 1 −NT . Φ 2 .8-13. ∂ Φ ....variable matrix.. Φ.. Φ n l q T The field gradient vector g is defined as g = T R∂ Φ . ∂ Φ U S∂ x ∂ y ∂ z V T W The field gradient . can now be interpolated as: Φ = φ T + ζψ T Φ Φ = HΦ lq with the field variable in terms of the nodal variables: Φ = Φ1 .... NT i i 2 The scalar temperature field. 8-13. Consequently for the transformation of the cartesian derivatives into the natural derivatives only a 2 dimensional Jacobian matrix is required. The element thermal stiffness matrix in basic form may be defined as 206 . R ∂ U LM ∂x | ∂ξ | = M ∂ξ S ∂ V M ∂x | ∂η | M ∂η | | N T W or inverted ∂y ∂ξ ∂y ∂η OPR ∂ U PP| ∂∂x | S V | | PQ| ∂y | T W ∂ −1 ∂ −1 ∂ = J11 + J12 ∂x ∂ξ ∂η ∂ −1 ∂ −1 ∂ = J 21 + J 22 ∂y ∂ξ ∂η and an integration constant for the thickness is computed from: z= c ∂ 2 ∂ ζ→ = . The differential of the volume is given by dV = c J ζdξdη 2 where |J| is the 2x2 Jacobian determinant.Element Formulations B can be split into two matrices combining in-plane and through thickness terms: B = B + ζB 1 2 where LM ∂φ OP MM ∂ x PP ∂φ B =M MM 2∂ y PPP MM c ψ PP N Q T T 1 B2 T LM ∂ψ OP ∂ MM ∂ψx PP =M MM ∂y PPP 0 MN PQ T T T The restriction of constant layer thicknesses provides an uncoupling between the inplane coordinates and the through-thickness coordinate. 2 ∂z c ∂ζ where c is the depth of the element see figure 7. variable matrices can be left out of the integration through the thickness: As the matrices B1 and B 2 are independent of ζ.Field Elements K= z V B T k B dV where k is the principal thermal conductivity matrix for an orthotropic material. Therefore the strain-displacement matrices can be left out of the integration through the thickness: F B L∑ k dζ O B + B L∑ ζ k dζ O B I PQ MN z PQ GG MN z JJ c K = zz GG LM∑ ζ k dζ OP B + B LM∑ ζ k dζ OP B JJ 2 J dηdξ +B H N z Q N z Q K T 1 nlay n =1 ζlay n 1 T 1 nlay n =1 ζlay n 2 ξ η T 2 nlay n =1 ζlay n 1 T 2 nlay n =1 2 ζlay n 2 with B1 and B 2 as: LM J MM B = MJ MM MNM 1 −1 11 ∂φ T ∂φ T −1 + J12 ∂ξ ∂η T ∂φ ∂φ T −1 −1 + J22 21 ∂ξ ∂η 2 T ψ c OP PP PP PP QP B2 LM J MM = MJ MM MM N −1 11 ∂ψ T ∂ψ T −1 + J12 ∂ξ ∂η T ∂ψ ∂ψ T −1 −1 + J22 21 ∂ξ ∂η 0 T OP PP PP PP QP 207 . Therefore the field gradient . only k varies from layer to layer. only k varies from layer to layer. LMk k =M0 MN 0 x 0 ky 0 0 0 kz OP PP Q As the matrices B1 and B 2 are independent of ζ. variable matrices only have to be computed inplane. JNT4.9. e.9-1). two translational springs for plane elements or one translation and two rotations for plate bending elements (fig. JXS3) Joint elements are composed of translational and rotational springs that may be used to connect two nodes on adjacent finite elements (fig.g. JSH4. JF3. The field gradient . JPH3. JL43.7.7. JL46.1 Joints (JNT3.8-13 Nodal Configuration For Solid Composite Field Elements 7.9 Joint Elements 7.7. JRP3. JSL4. The number of springs utilised in each element type is chosen to be compatible with a corresponding finite element type.9-2).Element Formulations The through thickness dependency is condensed in the integration of the thermal conductivity matrix which makes the assembly of the element thermal stiffness matrix more efficient. Fig. 208 . JAX3. This is possible by restricting the element to a reasonably uniform thickness for a single layer. 4). however they may be used as nonlinear support conditions in geometrically nonlinear analysis.2 Evaluation of Stresses/Forces Element forces f are evaluated directly in the local Cartesian system from f = ka' where k and a are stiffness matrix and displacement vector in the local Cartesian system.9-5).9. Note.3 Nonlinear Formulation The joint elements may be employed in Materially nonlinear analysis using the nonlinear joint models (section 4. element output may be obtained at either nodal or Gauss points.9.7. 7. The joint element possess no geometrically nonlinear terms. The element mass matrix is lumped and formed directly from user input masses.Joint Elements The local element stiffness matrix is formulated directly from user input stiffness coefficients and is then transformed to the global Cartesian system.9-6). 4 noded element The local x-axis is defined by a line joining the first and third element nodes. For convenience.7. The local xy-plane is defined by the fourth element node and the element x-axis. Element strain output is evaluated in the local element system as and R∈ U R u |∈ | = | v S V S |∈ | |w T W T Rψ U Rθ |ψ | = |θ S V S |ψ | |θ T W T x y z x y z 2 2 2 U | V −w | W − u1 − v1 1 translational strain x2 y2 z2 − θ x1 − θ y1 − θz1 U | V | W rotational strain The element local axes are defined by 3 noded element The local x-axis is defined by a line joining the first and third element nodes. The local y-axis forms a right handed set with the xaxis such that the local z-axis is upwards (out of page) (fig. The local y and z-axes form a right handed set with the local x-axis (fig. 7. 209 . 7.Element Formulations 2 Y 1 Springs X Fig.9-1 Springs Connecting Two Nodes In 2-D (No Rotational Stiffness) Rotational springs y 3 1 2 Translational spring x z Fig.7.9-2 Joint Element For Plate Bending Elements (One Translation And Two Rotations) 210 . Joint Elements z y x Z Y θ y is given zero stiffness to allow rotation X Fig.9-3 Modelling A Hinged Connection Between Shell Elements Fig.7.7.9-4 Excavation Example Illustrating Use Of Joint Elements For Nonlinear Support Conditions 211 . as depicted in figure 7.9-7.9-6 Local Cartesian System For 4-Noded Joint Elements 7. 2 x X Fig.4 Use of Joints With Higher Order Elements To illustrate the behaviour of joints with higher order elements. For convenience.9. 2 x X Z Fig.Element Formulations y Y 3 1. the coordinate system is centred at the mid-node of the element.7. the elastic foundation must apply an equal and opposite stress per unit length. The beam is modelled by a single quadratic element of length two and the elastic foundation by three joint elements. consider the problem of the beam on the elastic foundation shown in figure 7. For equilibrium. The equivalent nodal loads are calculated using the principle of virtual work as 212 .9-8.9-5 Local Cartesian System For 3-Noded Joint Elements z Y y 4 3 1.7. A constant stress per unit length q along the top face is transmitted through the beam to the elastic foundation beneath. the nodal loads calculated must be identical to the nodal loads developed by the discrete joint elements connected to the nodes on the lower surface. q may be substituted in each of the nodal equilibrium equations resulting with discrete spring stiffnesses K of (see fig.Joint Elements z −1 qδu dx = q δu1 1 LM N LM 1 δu =q N3 z 11 −1 2 x x − 1 dx + δu 2 a f 1 + 4 1 δu 2 + δu3 3 3 OP Q ze 1 −1 1 − x 2 dx + δu3 j z 11 −1 2 x x + 1 dx a f OPQ where the virtual work of the load is calculated from the perturbation of shape functions particular to each node. then the spring stiffness at connecting nodes must be summed (fig.7.7. Note. the computer program may be used to calculate the ratio of joint stiffnesses. it is recommended that for nonlinear contact problems. and a unit face load applied.9-10). then the stress per unit length is related to k by q = ku For a constant deflection u along the lower face. there is a corresponding difference in the internal strain energy associated with each virtual perturbation. If the appropriate boundary nodes are restrained. If the elastic foundation has a stiffness per unit length of k. However. the nodal loads are different since the virtual work arising from the perturbation of the mid-node is larger than that of the side nodes. the resulting nodal reactions will correspond to the integrated shape functions. these are also the ratio of spring stiffnesses to be used. 213 . From equilibrium.9-9) At node 1 and 3 At node 2 K= K= 1 k 3 4 k 3 If the beam is modelled by two or more elements. linear elements should be used if possible as higher order elements poorly represent the discontinuities in the boundary conditions. This may result in either poor convergence or divergence of the solution. For non-central midside nodes. With a little ingenuity a variety of spring boundary conditions can be evaluated using appropriate loading and the program to calculate equivalent reactions. Element Formulations q k Fig.9-7 Beam On An Elastic Foundation l l q 3 4 2 1 8 y x 7 6 5 Fig.7.7.9-8 Finite Element Discretisation q/3 4q/3 q/3 y x k/3 4k/3 k/3 Fig.7.9-9 Discrete Joint Stiffnesses 214 . The tangential displacement w is positive in the direction of increasing j. u and v are positive in the direction of increasing x and y respectively and may be either axial or radial displacements depending on the definition of the axis of symmetry. unless CBF loading is applied.2 Standard Isoparametric Elements The geometry of the body is defined using the shape functions 215 . The element axes are defined in the cylindrical coordinate system xyz. QAX4F.v. 7. 7.10.1 Global and Local Coordinate Systems The structural mesh is specified in the global XY-plane and. QAX8F) These are axisymmetric solid elements that may be subjected to non-axisymmetric loading (fig.Joint Elements q/3 4q/3 q/3 q/3 4q/3 q/3 k/3 4k/3 k/3 k/3 2k/3 4k/3 k/3 Fig. TAX6F.7.w.10. with associated displacements u. where θ is the positive rotation defined by the right-hand coordinate system about the axis of symmetry. together with the global Z-coordinate.10-1). in which case the structure is restricted to be axisymmetric about the X-axis. the global coordinate axes form a right-hand coordinate system.7.10 Fourier Element Formulation (TAX3F.9-10 Summation Of Joint Stiffnesses 7. In general the structure may be axisymmetric about either the X or Y axis. V X. The nodal degrees of freedom of the element are u.v.U Figure 7. η are standard linear or quadratic isoparametric element shape functions for node i and m is the number of nodes. ηgYi i =1 m where N i ξ. η X i and Y= ∑ Ni bξ.Element Formulations X= ∑ b g i =1 m N i ξ.w 4 b g at each node in the cylindrical coordinate system 3 6 7 5 8 2 1 QAX4F 1 2 QAX8F 4 3 5 4 6 2 TAX6F 3 1 1 2 3 TAX3F Y.10-1 Nodal Configurations For Fourier Elements 216 . u= n=0 m ∑ m us cos nθ + n ∑ uan sin nθ n =1 m n =1 m m v= n=0 m ∑ vsn cos nθ + ∑ van sin nθ n =0 w= ∑ ws sin nθ − n n=0 ∑ wan cos nθ n=0. the discretised displacement is defined as u = N' a where a is the nodal displacement vector and for any node i a i = as .m represents the range of harmonics considered and superscripts 's' and 'a' denote the symmetric and asymmetric components. vs .10.v. u a ..w are given by the Fourier expansions.Fourier Elements where u.. aa i i T = u s .1. 7. w s ..3 Strain-Displacement Relationships The infinitesimal strain-displacement relationships for structures symmetric about the global X-axis are given by ∈x = ∂u ∂x ∈y = ∂v ∂y 1 ∂w v + Y ∂θ Y ∈θ =∈z = γ XY = ∂u ∂v + ∂Y ∂X 217 . For each harmonic. w a i i i i i i T and N' is the shape function matrix defined as LMN cos nθ N' = M 0 MN 0 i 0 N i cos nθ 0 0 0 N i sin nθ N i sin nθ 0 0 0 N i sin nθ 0 0 0 − N i cos nθ OP PP Q where N i are the standard isoparametric shape functions. v a . LM 1 / E MM−υ 0 / E D=M MM 0 MM 0 N 0 x xy x − υyx / E y 1 / Ey 0 0 0 0 0 0 0 0 0 0 0 1 / G xy 0 0 0 0 0 0 0 0 0 0 0 0 0 OP PP P 0P 0P P 0P Q 218 .Element Formulations γ Yθ = γ YZ = γ Xθ = γ XZ = ∂w W 1 ∂v − + ∂Y Y Y ∂θ 1 ∂u ∂w + Y ∂θ ∂X 7. The constitutive relationship is. a plane stress material model may also be utilised. υ xz υzx υ yz υzy = . which for isotropic and orthotropic elasticity is defined as Isotropic 0 0 0 O LMa1 − υf υ υ P υ 0 0 0 P a1 − υf υ MM 0 P MM υ υ a1 − υf a1 −02υf 0 P E 0 0 0 0 0 P D= 2 a1 + υfa1 − 2υf MM a1 − 2υf 0 PPP 0 0 0 0 MM 2 a1 − 2υf PP MM 0 0 0 0 0 P 2 Q N Orthotropic LM 1 / E E MM−υ // E −υ D=M MM 0 MM 0 N 0 x xy xz − υ yx / E y x x − υzx / E z − υzy / E z 1 / Ez 0 0 0 0 0 0 1 / G xy 0 0 0 0 0 0 1 / G yz 0 0 0 0 0 0 1 / G xz 1 / Ey − υ yz / E y 0 0 0 OP PP PP PP QP −1 where for symmetry υ xy Ex = υ yx Ey . = Ex Ez Ey Ez In addition to the solid material definition of the element.4 Constitutive Relationships The modulus matrix D.10. 10-2a. the global loads are related to the local loads via LMP MMP NP x1 y1 z1 OP LM1 PP = MM0 Q N0 0 cos θ − sin θ OPLMP sin θ P M P cos θP M P QN 0 XG YG ZG OP PP Q OPLMP sin θ cos nθ P M P cos θ sin nθ P M P QN 0 0 Applying the appropriate virtual perturbations gives For symmetric contributions LMP MMP NP x1 cos nθ y1 z1 OP LMcos nθ cos nθP = M 0 sin nθ P M 0 Q N OP LMsin nθ sin nθP = M 0 cos nθP M 0 Q N 0 cos θ cos nθ − sin θ sin nθ XG YG ZG OP PP Q OP PP Q For asymmetric contributions LM P MM P N− P x1 sin nθ y1 0 cos θ sin nθ sin θ cos nθ z1 OPLMP sin θ sin nθ P M P − cos θ cos nθP M P QN XG YG ZG 219 . initial strains and thermal loading are applied in the local xyz directions. 7. σ xy . σ y and σ xy The plane stress material option is intended to allow the modelling of fan blades for which the use of the full modulus matrix is inappropriate.10. surface tractions. Note that concentrated loads (and nodal reactions) are applied as forces per unit length of the structural surface. σ z .5 Element Loading Concentrated loads.7. initial stresses. in contrast. A complete description of the element formulation is given in [C1]. (i) Axisymmetric about x-axis From fig. σ y . Two possible axes of symmetry must be considered. The resolution of the global loads into the local tangential and radial directions is performed using the following matrices. σ yz and σ xz to σ x . constant body forces and body force potentials are applied in the global XYZ directions. Note that elements using this material model should be adequately restrained in the tangential w direction.Fourier Elements The use of this material model results in a reduction of the active stresses from σ x . Element Formulations (ii) Axisymmetric about y-axis From fig. 220 . the global loads are related to the local loads via LMP MMP NP x1 y1 z1 OP LM cos θ PP = MM 0 Q N− sin θ − sin θ PXG 1 0 PYG 0 − cos θ PZG 0 OPLM PPMM QN OP PP Q OPLM PPMM QN OP PP Q Applying the appropriate virtual perturbations gives For symmetric contributions LMP MMP NP x1 cos nθ y1 z1 OP LM cos θ cos nθ 0 cos nθP = M sin nθ P M− sin θ sin nθ Q N OP LMcos θ sin nθ 0 sin nθP = M cos nθP Msin θ cos nθ Q N − sin θ cos nθ PXG 0 cos nθ PYG − cos θ sin nθ PZG 0 0 For asymmetric contributions LM P MM P N− P x1 sin nθ y1 z1 − sin θ sin nθ PXG 0 sin nθ PYG 0 cos θ cos nθ PZG 0 OPLM PPMM QN OP PP Q For dynamic and harmonic response analyses where the automatic evaluation of Fourier loads is not available.7.10-2b. the global loads for the 'nth' harmonic must be converted to local loads using the above expressions. 10.6 Inertial Loading The inertial loading due to angular accelerations and rotations require explicit evaluation and are directly applied to the element. 221 .10-2 Local And Global Loads 7.Fourier Elements Z Pz PZ Py PY θ x Y (a) Axi-symmetric about the global X-axis Pz Px PX θ Y PZ X Z (b) Axi-symmetric about the global Y-axis Fig 7. 10-3). Y'. X'. To apply the resulting inertial forces in the cylindrical coordinate system. Z define the shift in the global coordinate system. the accelerations are resolved (fig. Y. Substituting for r gives the accelerations with respect to the fixed system. Z) with respect to a fixed system r' ( X' . Y' and Z'.7.Element Formulations Assuming an arbitrary origin for the XYZ coordinate system (as shown in fig. Z') axes and && && && && = Xi + Y j + Zk r are the linear accelerations in the global X'.7. about which the angular accelerations and velocities are to be applied. Y' . results in the following definition for the displacement vector r r = X + x i + Y + r cos θ j + Z + r sin θ k d i d i d i where X.10-4) && = &&' x x && L = &&' cos θ + &&' sin θ y y z &&L = − &&' sin θ + &&' cos θ z y z 222 . Y' and Z' directions. Y. Definition of the acceleration vector The instantaneous acceleration of a particle in a rotating coordinate system r( X. Z' ) is [S4] &&' = && + 2Ωxr + Ωx Ωxr + αxr & r r b g where a dot signifies the derivative with respect to time and the vectors Ω = Ω X i + ΩY j + Ωz k α = αX i + αY j + αz k are the angular velocities and accelerations about the r'(X'. α X Ω Y.10-4 Resolution Of Global To Local Accelerations 223 .Fourier Elements z zl Z yl Y θ x y Figure 7.10-3 Definition Of Rotation Axes Origin Z Ω Z.α Z z θx x r θ y o Z Ω X.α Y Y X X Y Figure 7. Element Formulations which gives the following acceleration terms in cylindrical coordinates &&' = c1 + c2 cos θ + c3 sin θ + c 4 cos 2θ + c5 sin 2θ r where c1 c2 c3 c4 c5 LM&& − eΩ + Ω jdX + xi + bΩ Ω − α gY + bΩ Ω + α gZOP x 1 M PP & =M − F Ω + Ω + 2eΩ + θj I r K 2H MM PP θ e&& + α jr N Q OP LM & eΩ eΩ + 2θj − α j r && = MY + bΩ Ω + α gd X + xi − eΩ + Ω jY + bΩ Ω − α g Z P MM && PP NMZ + bΩ Ω − α gdX + xi + bΩ Ω + α gY − eΩ + Ω jZ QP OP LM & eΩ eΩ + 2θj + α j r && = M Z + bΩ Ω − α gd X + xi + bΩ Ω + α gY − eΩ + Ω jZ P MM && P − Y − bΩ Ω + α gd X + xi + eΩ + Ω jY − bΩ Ω − α gZ P MN PQ L 0 OP 1 M = r MΩ − Ω P 2 MN 2Ω Ω PQ L 0 O 1 = r M 2Ω Ω P 2 M MNΩ − Ω PPQ 2 Y 2 Z X Y Z X Z Y 2 Y 2 Z 2 X X Y X Z X Y Z 2 X 2 Z Y Z X X Z Y Y Z X 2 X 2 Y Z X Y X Z Y Y Z X 2 X 2 Y X Y Z 2 X 2 Z Y Z X 2 Y 2 Z Y Z 2 Z Y Z 2 Y Using d'Alembert's principle. the inertia force may be included as part of the load vector Rb = ze j e v N' T f − ρ&& dv x j = R b( f ) + R b(a ) where R b( f ) is the element nodal and body forces and R b( a ) is the inertial force vector which includes the effects of angular velocities and accelerations R b( a ) = − zzzbg + 1 +1 2 π −1 − 1 0 N' T ρ && Y J dξ dη x 224 . σ y are the coefficients of the Fourier series expansion of the stresses and strains.10. Centripetal forces are proportional to the angular rotation squared and the lever arm of the mass from the centre of rotation. ∈y . contrary. but there is no nonlinear stress stiffening contribution. which increases the applied load. INT16) 7. ∈z . actually decreases the stiffness of the structure. The centripetal load stiffness is applied by default. 7. They may be inserted at planes of potential delamination to model inter-laminar failure. The centripetal load stiffening matrix. σ yz . 7. 7. γ xz where σ x . This may.7 Evaluation of Stresses The stresses are evaluated at the element Gauss points and are extrapolated to the nodal points. and crack initiation and propagation.11. The principal stresses and strains are evaluated at θ = 0 only. 15 14 13 zero 6 5 3 4 1 2 zero 16 7 8 10 6 11 12 4 5 9 1 2 3 225 . but it may be omitted by setting OPTION 102. conversely. the lever arm is lengthened by positive displacements. These elements have no geometric properties and are assumed to have no thickness (see Fig. σ xz direct and shear stresses direct and shear strains Strain Output ∈x .1 Definition and interpolation Interface elements INT6 (6-noded line element for 2D analyses) and INT16 (16-node surface element for 3D analyses) are used to model composite delamination in an incremental non-linear analysis. to its name. As the body spins.11 Interface Elements (INT6.11-1). γ xy .Fourier Elements Centripetal load stiffening Centripetal load stiffening has been applied to the n = 0 term. The output consists of Stress Output σ x . be thought of as reducing the stiffness. σ y . σ xy . σ z . γ yz . 11.11-1 Interface Elements INT 6 and INT16 The displacement field u for these elements contains the bottom displacement u b and the top displacement u t (number of components in u b and u t is two for INT6 and three for INT16) so that u T = u T u T . The vector of nodal internal forces can be written in a standard form as 226 . The bottom and the top displacement are b t interpolated as u b = Hp b b and u t = Hp t where p and p are the vectors of the bottom and the top nodal displacements t (number of components in these vectors is six for INT6 and twenty-four for INT16) and H is the matrix of shape functions of the type INT6: Lh H=M NM0 T T 0T hT OP QP or INT16: LMh H = M0 MN0 T T T 0T hT 0T 0T 0T hT OP PP Q where ith component in the vector h (i=3 for INT6 and i=8 for INT16) is the value of the ith shape function at a particular point. The actual constitution of the interface element is defined in terms of the relative displacements between the bottom and the top surfaces ε = ut − ub = B Rp U = Bp | | Sp V | | T W b t where p is the vector of the nodal displacements (which has twelve components for INT6 and forty-eight components for INT16) and the matrix B follows as LM−h MN 0 LM−h B= M 0 MN 0 B= T hT 0T hT 0T 0T 0T −hT 0T −hT 0T 0T hT 0T hT 0T T T T T OP for INT6 and PQ 0 0 O P 0 0 P for INT16.2 Internal force vector and stiffness matrix The element equilibrium equation in a vector form is given as g≡P−R=0 where R is the vector of applied loading and P is the vector of nodal internal forces.Element Formulations Fig.7. −h h P Q T T T T T T 7. Interface Elements P = B σdA A z T where the integration domain A is a line for INT6 and an area for INT16. the internal force vector and the stiffness matrix are integrated using a Newton-Cotes integration rule rather than a reduced or full Gauss integration rule. The 3-point Newton-Cotes scheme is utilised for INT6 and the 3*3-point NewtonCotes scheme is utilised for INT16. For a given arbitrary constitutive relationship (note that ε is a relative displacement between the element’s surfaces rather than a strain measure) σ=σ ε bg T where σ is a stress vector. The interface elements INT6 and INT16 can currently be used only with the delamination damage model (non-linear material model 25).ij = t ∂σ i for a given ∂ε j material model. In order to eliminate spurious oscillations of the stress field along the element [H13]. which follows from D t. 227 . the stiffness matrix in a geometrically linear analysis follows as K = B D BdA A t z where D is a tangent constitutive matrix. . BEAMS.000000000 TABLE A-2 SAMPLING POINTS AND WEIGHTS FOR 5-POINT RULE FOR 2D QUADRILATERALS AND SHELLS 229 .A-1 to fig.20000000 0.8611363116 ± 0.7745966692 2.0000000000 ± 0.592348877 ± 0.0000000000 0.000000000 0.592348877 0. ORDER LOCATION ξi WEIGHT WI 1 2 3 4 0.5555555555 0.0000000000 1.QUADRILATERAL 2-D SOLIDS. SHELLS AND 3-D HEXAHEDRA AND PENTAHEDRA RULE ει LOCATION ηι WEIGHT 5 point ± 0.3478548454 0.6521451549 0.A-7.Appendix A Appendix A Quadrature Rules The locations and weights of the quadrature points used in integrating the element matrices are listed in table A-1 to table A-7 and are shown in fig.3399810436 TABLE A-1 SAMPLING POINTS AND WEIGHTS FOR BARS. PLATES.5773502692 ± 0.8888888888 0.00000000000 ± 0.95000000 0. 50000000 0.3333333333 0.5208333333 0.1259391805 TABLE A-3 SAMPLING POINTS AND WEIGHTS FOR TRIANGULAR 2-D SOLIDS.25000000 0.00000000 0.3333333333 0.0000000000 0.3333333333 0.0000000000 0.6000000000 0. PLATES.4701420641 0.3333333333 0.25000000 0.2000000000 0.0000000000 0.25000000 0.0597158717 0.3333333333 -0. SHELLS AND 3-D PENTAHEDRA RULE A1 LOCATION A2 A3 WEIGHT 3-Point 1.58541020 0.2000000000 0.25000000 0.3333333333 0.3333333333 TABLE A-4 SAMPLING POINTS AND WEIGHTS FOR TRIANGULAR SEMILOOF SHELL RULE V1 V2 LOCATION V3 V4 WEIGHT 1-Point 2-Point 3-Point 0.Appendix A RULE A1 LOCATION A2 A3 WEIGHT 1-point 3-point 4-point 7-point 0.0000000000 0.0000000000 0.13819660 0.25000000 0.1012865073 0.5625000000 0.16666666 TABLE A-5 SAMPLING POINTS AND WEIGHTS FOR 3-D TETRAHEDRA 230 .50000000 0.00000000 1.3333333333 0.00000000 0.13819660 0.4701420641 0.0000000000 0.5000000000 0.7974269853 0.13819660 0.3333333333 0.1012865073 1.3333333333 0.1323941527 0.2250000000 0.3333333333 0. 333333333 ± 1.758786911 -0.Appendix A RULE ξi LOCATION ηi ζi WEIGHT 13-Point 0.000000000 ± 0.000000000 -0.333333333 0.000000000 -0.79562143 0.335180055 0.750000000 0.250000000 0.00000000 ± 0.507644216 0.795822426 0.500000000 0.758786911 0.758786911 -0.000000000 ± 1.025293237 1.000000000 2.000000000 ± 0.758786911 0.54498736 0.49584802 ± 0.000000000 ± 1.684210565 0.00000000 ± -0.88030430 ± 0.SAMPLING POINTS AND LOCATIONS FOR NEWTON-COTES RULES 231 .711111111 0.758786911 -0.79562143 0.000000000 0.758786911 0.758786911 0.266666667 0.49584802 ± 0.000000000 ± 1.355555556 0.758786911 TABLE A-6 .758786911 -0.166666667 1.000000000 TABLE A-7 .335180055 14-Point ± 0.758786911 0.335180055 0.SAMPLING POINTS AND WEIGHTS FOR SPECIAL RULES FOR 3-D SOLIDS RULE LOCATION WEIGHT 1-Point 2-Point 3-Point 4-Point 5-Point 0.155555556 0.758786911 0.000000000 1.00000000 ± -0.758786911 0.335180055 0. BEAM AND AXISYMMETRIC SHELL ELEMENTS 232 .A-1 GAUSS QUADRATURE RULES FOR BAR.Appendix A 1 1 2 (a) 1-Point Rule 1 1 2 2 (b) 2-Point Rule 1 1 2 3 2 (c) 3-Point Rule 1 1 2 3 4 2 (d) 4-Point Rule FIG. A-2 GAUSS QUADRATURE RULES FOR QUADRILATERAL 2-D CONTINUUM. PLATE AND SHELL ELEMENTS 233 .Appendix A 4 4 3 6 8 1 2 1 2 3 3 7 9 5 6 8 4 4 7 5 1 (a) 2*2 Rule 2 1 2 (b) 3*3 Rule 3 7 13 9 8 5 1 1 6 2 6 14 10 15 11 7 3 2 (a) 4*4 Rule 16 12 5 7 4 6 3 5 5 4 8 8 1 2 4 4 3 1 2 (b) 5-Point Rule 3 FIG. A-4 SPECIAL 3-POINT RULE FOR TRIANGULAR SEMILOOF SHELL ELEMENT 234 . PLATE AND SHELL ELEMENTS 3 3 1 1 2 2 FIG.Appendix A 3 3 3 1 1 1 (a) 1-Point Rule 2 1 2 (b) 3-Point Rule 2 3 2 1 3 5 7 4 1 2 (b) 7-Point Rule 4 3 1 (a) 4-Point Rule 2 2 3 6 1 FIG.A-3 GAUSS QUADRATURE RULES FOR TRIANGULAR 2-D CONTINUUM. Appendix A 4 1 1 3 2 (a) 1-Point Rule 4 3 4 1 1 2 3 2 (b) 4-Point Rule 4 6 3 5 1 4 2 (c) 6-Point Rule 2 3 1 FIG.A-5 GAUSS QUADRATURE RULES FOR SOLID TETRAHEDRA ELEMENTS 235 . A-6 QUADRATURE RULES FOR SOLID PENTAHEDRA AND HEXAHEDRA ELEMENTS 236 .Appendix A 8 5 8 3 5 6 4 4 5 5 6 6 7 7 4 2 3 1 2 1 2 1 1 2 3 3 (a) 3*2 Rule (b) 2*2*2 Rule 8 16 5 10 13 11 7 4 1 2 5 3 14 12 8 6 6 17 18 15 7 8 25 5 19 6 18 26 27 7 9 3 10 4 1 2 11 5 3 2 (d) 3*3*3 Rule 12 6 9 3 1 (c) 3*3*2 Rule 2 1 FIG. Appendix A 1 (a) 1-Point Rule 1 2 (b) 2-Point Rule 1 2 3 (c) 3-Point Rule 1 2 3 4 (d) 4-Point Rule (e) 5-Point Rule FIG.A-7 NEWTON-COTES RULESAppendix B 237 . Appendix A 238 . a small amount of warping is permitted provided that z < 0.B-1.02 where a. However.05 (a + b)/c < 1.B-1 DEFINITION OF PARAMETERS FOR CURVATURE LIMITS Warping of Flat Elements The four nodes defining a flat quadrilateral element in 3-D should be coplanar. b and c are defined in fig.01 a where and z is the distance of the out of plane node from the plane a is the length of the side between the first and second nodes. b a c  Fig.Appendix B Appendix B Restrictions On Element Topology Mid-Length and Mid-Side Nodes The mid-length and mid-side nodes of elements should be equidistant from the two end nodes.b|/(a + b) < 0. and the element curvature must satisfy the following requirements (i) (ii) |a . 239 . Appendix B 240 . 241 .References References Contact [email protected] for details of all references stated in this manual. References 242 .