Isoparametric Elements Element Stiffness Matrices

1. Isoparametric ElementsElement Stiffness Matrices Structural Mechanics Displacement-based Formulations

2. General Approach – Specific Example • We will look at manipulation of the mechanics quantities (displacement, strain, stress) using shape functions • The approach is quite general, and is used to formulate a number of different elements • We will use a specific example to make the development more concrete (Q4) • We will start from the nodal displacement representation, work toward strain and stress, and finally element stiffness • There is a lot going on here, pay attention to both the overall themes and the detailed steps …

3. Master Element Mapping • Note: we will use a for x and b for h because I can’t remember, pronounce, or legibly write “xi” and “eta” master element actual element

4. Bilinear Quadrilateral (Q4) • Interpolation involves the summation of nodal values multiplied by corresponding shapes functions geometry interpolation field variable interpolation - where - nodal coordinates nodal displacements shape functions

5. Q4 - Displacements • Start with the element displacement field • We have to give it some functional form in order to work with it • Let it be defined over the element by our interpolation scheme {u} = displacements continuously defined (all components) over an element [N] = the element shape functions in master element coordinates {d} = the nodal (discrete) displacement values

6. Strain from {u} • Now calculate element strains from the displacement field • This is just the usual strain-displacement relationship written in compact form with an operator matrix

7. Q4 – Strain from {d} • Let’s now work toward an expression for element strain • We have a bit of a difficulty here with direct substitution • The shape functions (N1, N2, N3, N4) are defined in terms of the master element coordinates (a,b) • But we need to differentiate in terms of the global coordinates (x,y) this operation cannot be done directly - or -

8. Coordinate Transformation • Given any function of the master element coordinates (a,b): • We can find derivatives with respect to global (x,y) by using the chain rule: • We can combine and rearrange these relationships to get our derivatives:

9. The Jacobian • The Jacobian matrix is an important part of element formulation: • For the Q4 element this becomes: note the Jacobian matrix is a function of location within the master element local coordinate derivatives of the shape functions global coordinate locations of the element nodes

10. Jacobian Interpretation • The Jacobian contains information about element size and shape • The Jacobian determinant (j) is a scaling factor that relates the differential area of the actual element to the differential area of the master element • The Jacobian inverse (G) relates global coordinate system (x,y) function derivatives to master element coordinate system (a,b) function derivatives

11. Jacobian (determinant) Ratio • This is one measure of element quality (which affects element accuracy) • Ratio of the highest to lowest quadrature point Jacobian determinant • It is 1.0 for any square or rectangular element (same j throughout element) • It increases as element distortion increases

12. Strain/Displacement for Q4 • Start with the usual strain-displacement relationship in a slightly different form: • Now add the Jacobian approach to master/global coordinate derivative transformation:

13. Strain/Displacement cont. • Now represent the displacement field master element derivatives in terms of the shape functions:

14. All Together Now … [B] organization shape function derivatives, master coordinates - or - Jacobian inverse terms, master to global coordinate transformation nodal displacements, global coordinates

15. Stress • If we have strain, we can get to stress by bringing in material properties • We have to be a little careful here, this simple expression assumes: • No initial (residual, assembly) stresses present • Linear elastic behavior • The general form above does accommodate anisotropic behavior • If we further limit ourselves to 2D, isotropic, plane stress, we can write:

16. Element Stiffness Matrix • Recall where the element stiffness matrix fits into the finite element formulation: • Take it as a given for the present that the element stiffness matrix [k] is: • An integral over the element area in global coordinates (t = thickness) • Why is integration required? • Think about what [k] does • For displacements applied to the element nodes, it determines the required force • If one element is larger than another, the force required ought to be greater for the same nodal displacements • If an element has a rotated orientation, a coordinate axis displacement can produce forces with multiple coordinate components

17. Integration in Master Coordinates • It is not easy to integrate for the terms in [k] using the global coordinate system (elements are generally distorted and not aligned with global axes) • But we can do this instead (matrix dimensions for a Q4 element): • Integrate over the master element • It is undistorted and aligned with the coordinate system • Adjust for the change in coordinates by bringing in the Jacobian determinant j

18. Quadrature • Read “quadrature” as “numerical integration” • Why do we want to numerically integrate to establish [k]? • To integrate directly is still computationally expensive, even with the change to local coordinates • Quadrature involves sampling at discrete points, multiplying by a weighting factor, and summing to get an estimate of the integral this varies point-by-point too … these contain Jacobian inverse terms which vary point-by-point within the element

19. Gauss Points • Gauss quadrature is a method of numerical integration that has optimal characteristics when the underlying functions have polynomial form • The figure shows Gauss points for 2nd order and 3rd order quadrature • For (a), all four points have a weight of 1.0 (total = 4.0) • For (b): 1,3,7,9 weight = .3086; 2,4,6,8 weight = .4938; 5 weight = .7901 (total = 4.0) • Note: the quadrature rule is independent of element order (Q4, Q8, Q9)

20. Computational Procedure • Clear the array that will contain [k] • Loop over integration points in the a direction (i=1 to ni) • Set sampling point location ai and the weight factor Wi (i.e. +/- 0.577, 1.0) • Loop over integration points in thebdirection (j=1 to nk) • Set sampling point locationbkand the weight factor Wk (i.e.+/- 0.577, 1.0) • Call shape function subroutine, return [B] and j=det[J] at the point (ai, bk) • Calculate [B]T[E][B]tjWiWk, and add into [k] • End loop k • End loop i • Shape function subroutine • Input is the local Gauss point coordinates (a,b) • Access database for element global node coordinates (i.e. x1,y1,x2,y2,x3,y3,x4,y4) • Calculate [J], j = det[J], [g] = [J]-1 and [B]

21. Efficiency • One-time element calculations • The shape function local derivatives are the same for a given element type and integration scheme • They are only calculated once and stored each time the shape function subroutine is called • File vs. DRAM storage • If you had to access a hard drive file each time a node location was required, the shape function subroutine calculations could be slow • But you need enough DRAM to store the locations, along with all the other information required during the analysis

22. Element Distortion • One of the reasons a distorted element is less than ideal: • The integral is estimated by discrete sampling at specific locations within the element • If the element is not distorted, the sampled points are highly representative of the un-sampled near by regions of the element • If the element is highly distorted, the sampled points are not representative of the un-sampled regions of the element

23. Intra-Element Jacobian Variation • Here is a single Q4 element (highly-distorted, not recommended) • For integration purposes, element characteristics are sampled at the discrete Gauss points shown • How well does the point represent the region around it? approx area represented by the point Gauss point

24. Geometric Nonlinearity • The configuration of the element changes as soon as it deforms. • The global locations of the nodes is different at each load increment. • Since global node locations enters into the element stiffness calculation … • The element stiffness evolves during a high-deformation sequence • A single load step is often broken up into a series of increments … • And the element stiffness is updated at each increment. original element load increment 1 load increment 2 …