김경종
23 KiB
On studying the cause of loss of predictive capability, we find that this effect is due to the elements no longer being able to represent the same order of polynomials in the physical coordinates x, y, z after the geometric distortion as they did without the distortion. For example, the general quadrilateral nine-node element, shown in Fig. 5.16(d), is able to represent the x^2, xy, y^2 displacement variations exactly, whereas the corresponding eight-node quadrilateral element is not able to do so. Hence, the general quadrilateral eight-node element does not contain the quadratic terms in the Pascal triangle of the physical coordinates.
This observation explains the results in Fig. 5.17, and an investigation into this phenomenon for widely used elements and common distortions is of general interest. In such a study, we can measure the loss of predictive capability by identifying which terms in the physical coordinates of the Pascal triangle can no longer be represented exactly, (see N. S. Lee and K. J. Bathe [A]).
Let us consider the two-dimensional element in Fig. 5.16 as an example. For elements with undistorted configurations or with aspect ratio distortions only, the physical coordinates x, y are linearly related to the isoparametric coordinates r, s ; i.e., we have x = c_1 r , y = c_2 s , where c_1 and c_2 are constants. Hence, the Pascal triangle terms in physical coordinates are simply the r, s terms obtained from the interpolation functions h_i replaced by x and y , respectively.
The effects of the parallelogram, general angular, and curved edge distortions shown in Figs. 5.16(c) to (e) can be studied by establishing the physical coordinate variations for these specific cases with the coordinate interpolations (5.18) and then asking what polynomial terms in x and y are contained in the r, s polynomial terms of the displacement variations given in (5.19) (see Example 5.21).
Table 5.1 summarizes the results obtained in such a study for two-dimensional quadrilateral elements (see N. S. Lee and K. J. Bathe [A]). The first column in Table 5.1 gives the terms in the Pascal triangle when the element is undistorted or is subjected to an aspect ratio or parallelogram distortion only. The terms below the dashed line are present only when the element is undistorted, or subjected only to an aspect ratio distortion, and also unrotated. Table 5.1 in particular shows that a general angular distortion significantly affects the predictive capability of 8- and 12-node elements; i.e., with such distortions the elements can represent only linear displacement variations in x and y exactly, whereas the 9- and 16-node elements can in distorted form still represent the parabolic and cubic displacement fields exactly.
On the other hand, curved edge distortions reduce the order of displacement polynomials that can be represented exactly for all the elements considered in Table 5.1, and indeed with such distortions only the biquartic 25-node element can still represent the parabolic displacement field exactly.
While the information given in Table 5.1 shows clearly that the Lagrangian elements are preferable to the 8- and 12-node elements in terms of predictive capability, of course, we also need to recall that the Lagrangian elements are computationally slightly more expensive, and for fine meshes the order of convergence is identical [although the constant c in (5.53) is different].
We demonstrate the procedure of analysis to obtain the information given in Table 5.1 in the following example.
TABLE 5.1 Polynomial displacement fields in physical coordinates that can be solved exactly by various elements in their undistorted and distorted configurations ^{ \dagger }
| Type of element | Fields for undistorted configuration, aspect ratio and/or parallelogram distortions | ||
| Additional fields if also unrotated | Angular distortion | Quadratic curved-edge distortion | |
| 8-node element | 1 | 1 | 1 |
| $x \ y$ | $x \ y$ | $x \ y$ | |
| $x^{2} \ xy \ y^{2}$ $x^{2}y \ xy^{2}$ | |||
| 12-node element | 1 | 1 | 1 |
| $x \ y$ | $x \ y$ | $x \ y$ | |
| $x^{2} \ xy \ y^{2}$ $x^{3} \ x^{2}y \ xy^{2} \ y^{3}$ $x^{3}y \ xy^{3}$ | |||
| 9-node Lagrangian element | 1 | 1 | 1 |
| $x \ y$ | $x \ y$ | $x \ y$ | |
| $x^{2} \ xy \ y^{2}$ $x^{2}y \ xy^{2}$ $x^{2}y^{2}$ | $x^{2} \ xy \ y^{2}$ | ||
| 16-node Lagrangian element | 1 | 1 | 1 |
| $x \ y$ | $x \ y$ | $x \ y$ | |
| $x^{2} \ xy \ y^{2}$ $x^{3} \ x^{2}y \ xy^{2} \ y^{3}$ $x^{3}y \ x^{2}y^{2} \ xy^{3}$ $x^{3}y^{2} \ x^{2}y^{3}$ $x^{3}y^{3}$ | $x^{2} \ xy \ y^{2}$ $x^{3} \ x^{2}y \ xy^{2} \ y^{3}$ | ||
| 25-node Lagrangian element | 1 | 1 | 1 |
| $x \ y$ $x^{2} \ xy \ y^{2}$ $x^{3} \ x^{2}y \ xy^{2} \ y^{3}$ $x^{4} \ x^{3}y \ x^{2}y^{2} \ xy^{3} \ y^{4}$ $x^{4}y \ x^{3}y^{2} \ x^{2}y^{3} \ xy^{4}$ $x^{4}y^{2} \ x^{3}y^{3} \ x^{2}y^{4}$ $x^{4}y^{3} \ x^{3}y^{4}$ $x^{4}y^{4}$ | $x \ y$ $x^{2} \ xy \ y^{2}$ $x^{3} \ x^{2}y \ xy^{2} \ y^{3}$ $x^{4} \ x^{3}y \ x^{2}y^{2} \ xy^{3} \ y^{4}$ | $x \ y$ $x^{2} \ xy \ y^{2}$ | |
^{\dagger} Two-dimensional quadrilateral elements are considered.
EXAMPLE 5.21: Consider the general angularly distorted eight-node element in Fig. E5.21. Evaluate the Pascal triangle terms in x, y for this element.
The physical coordinate variations are obtained by using the interpolation functions in Fig. 5.4, which give for this element with its midside nodes placed halfway between the corner nodes,
x = \gamma_ {1} + \gamma_ {2} r + \gamma_ {3} s + \gamma_ {4} r s \tag {a}
y = \delta_ {1} + \delta_ {2} r + \delta_ {3} s + \delta_ {4} r s \tag {b}
text_image
s y x r
Figure E5.21 Eight-node isoparametric element, with angular distortion
with \gamma_{1},\ldots ,\gamma_{4} and \delta_1,\ldots ,\delta_4 constants. We use (a) and (b) to identify which x and y terms are contained in the displacement interpolations
u = \sum_ {i = 1} ^ {8} h _ {i} u _ {i} \tag {c}
v = \sum_ {i = 1} ^ {8} h _ {i} v _ {i} \tag {d}
where the h_i are again those in Fig. 5.4.
Consider the u-displacement interpolation. The constant and x and y terms in (a) and (b) are clearly contained in (c) because (c) interpolates u in terms of the functions (1, r, s, r^{2}, rs, s^{2}, r^{2}s, rs^{2}) multiplied by constants. We discussed this fact earlier when considering the requirements for convergence.
However, if we next consider the term x^2 , we notice that the term r^2 s^2 [obtained by squaring the right-hand side of (a)] is not present in (c). Similarly, the terms xy, y^2, x^2y , and xy^2 are not present in the displacement interpolation (c).
The analysis for the v-displacement interpolation is of course identical. Hence, when an eight-node isoparametric element is subjected to a general angular distortion, the predictive capability is diminished in that quadratic displacement variations in x and y can no longer be represented exactly (see Table 5.1).
This analysis also shows that the quadratic displacement variations in x and y are retained when the nine-node displacement-based element is subjected to the same angular distortions. These conclusions explain the results shown in Fig. 5.16.
5.3.4 Element Matrices in Global Coordinate System
So far we have considered the calculation of isoparametric element matrices that correspond to local element degrees of freedom. In the evaluation we used local element coordinates x, y, and z, whichever were applicable, and local element degrees of freedom u_{i}, v_{i} , and w_{i} . However, we may note that for the two-dimensional element considered in Examples 5.5 to 5.7 the element matrices could have been evaluated using the global coordinate variables X and Y, and the global nodal point displacements U_{i} and V_{i} . Indeed, in the calculations presented, the x and y local coordinates and u and v local displacement components simply needed to be replaced by the X and Y global coordinates and U and V global displacement components, respectively. In such cases the matrices then would correspond directly to the global displacement components.
In general, the calculation of the element matrices should be carried out in the global coordinate system, using global displacement components if the number of natural coordinate variables is equal to the number of global variables. Typical examples are two-dimensional elements defined in a global plane and the three-dimensional element in
Fig. 5.5. In these cases the Jacobian operator in (5.24) is a square matrix, which can be inverted as required in (5.25), and the element matrices correspond directly to the global displacement components.
In cases where the order of the global coordinate system is higher than the order of the natural coordinate system, it is usually most expedient to calculate first the element matrices in the local coordinate system and corresponding to local displacement components. Afterward, the matrices must be transformed in the usual manner to the global displacement system. Examples are the truss element or the plane stress element when they are oriented arbitrarily in three-dimensional space. However, alternatively, we may include the transformation to the global displacement components directly in the formulation. This is accomplished by introducing a transformation that expresses in the displacement interpolation the local nodal point displacements in terms of the global components.
EXAMPLE 5.22: Evaluate the element stiffness matrix of the truss element in Fig. E5.22 using directly global nodal point displacements.
text_image
Y, V Y2 V2 U2 Node 2 V1 α r, u Y1 Node 1 U1 L X1 X2 X, U Young's modulus E Cross section A
Figure E5.22 Truss element in global coordinate system
The stiffness matrix of the element is given in (5.27); i.e.,
\mathbf {K} = \int_ {V} \mathbf {B} ^ {T} \mathbf {C} \mathbf {B} d V
where B is the strain-displacement matrix and C is the stress-strain matrix. For the truss element considered we have
u = \left[ \cos \alpha \quad \sin \alpha \right] \left[ \begin{array}{l} \frac {1}{2} (1 - r) U _ {1} + \frac {1}{2} (1 + r) U _ {2} \\ \frac {1}{2} (1 - r) V _ {1} + \frac {1}{2} (1 + r) V _ {2} \end{array} \right]
Then, using \epsilon = \partial u / \partial x , expressed in the natural coordinate system as \epsilon = (2 / L)\partial u / \partial r (see Section 5.2), we can write the strain-displacement transformation corresponding to the displacement vector \mathbf{U}^T = [U_1\quad V_1\quad U_2\quad V_2] as
\mathbf {B} = \frac {1}{L} [ \cos \alpha \quad \sin \alpha \quad \cos \alpha \quad \sin \alpha ] \left[ \begin{array}{c c c c} - 1 & & & \\ & - 1 & & \text { zeros } \\ & & & 1 \\ \text { zeros } & & & 1 \end{array} \right]
Also, as given in Section 5.2, we have
d V = \frac {A L}{2} d r \quad \text { and } \quad \mathbf {C} = E
Substituting the relations for B, C, and dV and evaluating the integral, we obtain
\mathbf {K} = \frac {A E}{L} \left[ \begin{array}{c c c c} \cos^ {2} \alpha & \cos \alpha \sin \alpha & - \cos^ {2} \alpha & - \cos \alpha \sin \alpha \\ \sin \alpha \cos \alpha & \sin^ {2} \alpha & - \sin \alpha \cos \alpha & - \sin^ {2} \alpha \\ - \cos^ {2} \alpha & - \cos \alpha \sin \alpha & \cos^ {2} \alpha & \cos \alpha \sin \alpha \\ - \sin \alpha \cos \alpha & - \sin^ {2} \alpha & \sin \alpha \cos \alpha & \sin^ {2} \alpha \end{array} \right]
5.3.5 Displacement/Pressure Based Elements for Incompressible Media
We discussed in Section 4.4.3 the fact that pure displacement-based elements are not effective for the analysis of incompressible (or almost incompressible) media and introduced two displacement/pressure formulations. In the u/p formulation, the pressure is interpolated individually for each element and can (in the almost incompressible case) be statically condensed out prior to the element matrix assemblage, whereas in the u/p-c formulation the element pressures are defined by nodal variables which, as for the displacements, pertain to adjacent elements. Various effective elements of these formulations were given (see Tables 4.6 and 4.7) and discussed (see Section 4.5).
As for the pure displacement-based elements, we assumed in Chapter 4 that the displacement and pressure interpolation matrices are constructed using the generalized coordinate approach. However, it is now apparent that these matrices can be obtained directly using the isoparametric formulation.
In the u/p formulation, we use the same coordinate and displacement interpolations for an element as in the pure displacement formulation [see (5.18) and (5.19)], and we interpolate the pressure using
p = p _ {0} + p _ {1} r + p _ {2} s + p _ {3} t + \dots \tag {5.54}
where p_0, p_1, p_2, p_3, \ldots are the pressure parameters to be calculated and r, s , and t are the isoparametric coordinates. Of course, as an alternative, we could also interpolate the pressure using
p = p _ {0} + p _ {1} x + p _ {2} y + p _ {3} z + \dots
where x, y , and z are the usual Cartesian coordinates.
In the u/p-c formulation we also use the displacement and coordinate interpolations as in the pure displacement formulation, and
p = \sum_ {i = 1} ^ {q _ {p}} \bar {h} _ {i} \hat {p} _ {i} \tag {5.55}
where the \tilde{h}_i, i = 1, \ldots, q_p are the nodal point pressure interpolation functions and the \hat{p}_i are the unknown nodal pressures. We note that the \tilde{h}_i are different from the h_i which are used for the displacement and coordinate interpolations. For example, for the 9/4-c two-dimensional element, the displacement and coordinate interpolations are the functions
corresponding to the nine element nodes in Table 5.4, whereas the \tilde{h}_i are the functions corresponding to the four element corner nodes.
In practice, the isoparametric formulation of the u/p and u/p-c elements is effective because the generality of nonrectangular and curved elements is then also available (see Fig. 4.21 and T. Sussman and K. J. Bathe [A, B]).
5.3.6 Exercises
5.1. Use the procedure in Example 5.1 to prove that the functions in Fig. 5.3 for the four-node truss element are correct.
5.2. Use the procedure in Example 5.1 to prove that the functions in Fig. 5.4 for the two-dimensional element are correct.
5.3. Use the functions in Fig. 5.4 to construct the interpolation functions of the six-node element shown. Plot the interpolation functions in an oblique/aerial view (as in Example 5.4).
text_image
s 2 5 1 6 r 3 4
5.4. Prove that the construction of the interpolation functions in Fig. 5.5 gives the correct functions for the three-dimensional element.
5.5. Determine the interpolation function h_i for the element shown for use in a compatible finite element mesh.
text_image
2/3 2/3 2/3 s 2 6 5 1 2/3 2/3 2/3 r 7 8 3 4
5.6. In the computation of isoparametric element matrices, the integration is performed over the natural coordinates r, s, t, which requires the transformation (5.28). Derive this transformation using the elementary volume dV = (\mathbf{r} \, dr) \times (\mathbf{s} \, ds) \cdot (\mathbf{t} \, dt) shown.
text_image
r, s constant t, r constant z y x t dt s ds r dr s, t constant dV = (r dr) × (s ds) • (t dt)
Here the vectors r, s, t are given by
\mathbf {r} = \left[ \begin{array}{l} \frac {\partial x}{\partial r} \\ \frac {\partial y}{\partial r} \\ \frac {\partial z}{\partial r} \end{array} \right]; \qquad \mathbf {s} = \left[ \begin{array}{l} \frac {\partial x}{\partial s} \\ \frac {\partial y}{\partial s} \\ \frac {\partial z}{\partial s} \end{array} \right]; \qquad \mathbf {t} = \left[ \begin{array}{l} \frac {\partial x}{\partial t} \\ \frac {\partial y}{\partial t} \\ \frac {\partial z}{\partial t} \end{array} \right]
5.7. Evaluate the Jacobian matrices for the following four-node elements.
text_image
6 4 y x Element 1 6 4 30° Element 2 6 4 30° Element 3 60°
Show explicitly that the Jacobian matrices of elements 2 and 3 contain a rotation matrix representing a 30-degree rotation.
5.8. Calculate the Jacobian matrix of the element shown for all r, s. Identify the values of r, s for which the Jacobian matrix is singular.
scatter
| Point | x | y | |---|---|---| | 1 | 7 | 9 | | 2 | 6 | 4 | | 3 | 3 | 2 | | 4 | 9 | 2 |
text_image
s 2 1 r 3 4
5.9. Evaluate the Jacobian matrices J for the following four-node elements.
text_image
y (3, 5) (3, 2) (8, 6) (8, 2) Element 1 x
text_image
y (1 + 5√3/2, 9/2 + 2√3) (3/2, 2 + 3√3/2) (3 + 5√3/2, 9/2) 30° Element 2 x
Show how the Jacobian matrix J of element 2 can be obtained by applying a rotation matrix to the Jacobian matrix J of element 1. Give this rotation matrix.
5.10. Consider the isoparametric elements given by
(a) Case 1:
x = \sum_ {i = 1} ^ {8} h _ {i} x _ {i}; \quad x _ {1} = 1 2, x _ {2} = 4, x _ {3} = 4, x _ {4} = 1 2, x _ {5} = 9, x _ {6} = 5, x _ {7} = 8, x _ {8} = 1 1
y = \sum_ {i = 1} ^ {8} h _ {i} y _ {i}; \quad y _ {1} = 1 2, y _ {2} = 8, y _ {3} = 2, y _ {4} = 2, y _ {5} = 8, y _ {6} = 5, y _ {7} = 1, y _ {8} = 7
(b) Case 2:
x = \sum_ {i = 1} ^ {6} h _ {i} ^ {*} x _ {i}; \quad x _ {1} = 8, x _ {2} = 2, x _ {3} = 1, x _ {4} = 9, x _ {5} = 5, x _ {6} = 5
y = \sum_ {i = 1} ^ {6} h _ {i} ^ {*} y _ {i}; \quad y _ {1} = 1 0, y _ {2} = 8, y _ {3} = 3, y _ {4} = 1, y _ {5} = 9, y _ {6} = 2
Draw the elements accurately on graph paper and show the physical locations of the lines r = \frac{1}{2} , r = -\frac{1}{4} , s = \frac{3}{4} , and s = -\frac{1}{3} for each case. (You may also write a small program to perform this task.)
5.11. Consider the isoparametric finite elements shown. Sketch the following for each element.
(a) The lines, s as the variable and constant r = -\frac{2}{3}, -\frac{1}{3}, 0, \frac{1}{3}, \frac{2}{3} .
(b) The lines, r as the variable and constant s = -\frac{2}{3}, -\frac{1}{3}, 0, \frac{1}{3}, \frac{2}{3} .
(c) The determinant of the Jacobian over element 1 (in an oblique aerial view).
text_image
y v₁ u₁ (3, 4) 2 1 (7, 5) 3 (2, 2) 4 (4, 2) Element 1 x
line
| Point | x | y | |---|---|---| | 1 | 6 | 6 | | 2 | 3 | 3 | | 3 | 2 | 2 | | 4 | 2 | 2 |5.12. Prove that for any parallelogram-shaped isoparametric element the Jacobian determinant is constant. Also, prove that the Jacobian determinant always varies with r and/or s whenever the element is not square, rectangular, or parallelogram-shaped.
5.13. Consider the isoparametric element shown. Calculate the coordinates x, y of any point in the element as a function of r, s and establish the Jacobian matrix.
text_image
y 2 s 1 0.5 5 6 9 8 r 4 1 1 3 2 2 4 x
5.14. Calculate the nodal point forces corresponding to the surface loading on the axisymmetric element shown (consider 1 radian).
text_image
1 radian 2 5 1 6 9 8 3 3 7 4 2 p 2 p y x p C
5.15. Consider the five-node plane strain isoparametric element shown.
(a) Evaluate appropriate interpolation functions h_i , i = 1 to 5.
(b) Evaluate the column in the strain-displacement matrix corresponding to the displacement u_{1} at the point x = 2.5 , y = 2.5 .
line
| Point | x | y | v₁ |
|---|---|---|---|
| 1 | 4.0 | 4.0 | u₁ |
| 2 | 1.0 | 3.0 | |
| 3 | 1.0 | 2.0 | |
| 4 | 4.0 | 1.0 | |
| 5 | 2.5 | 2.5 |
5.16. Consider the isoparametric axisymmetric two-dimensional finite element shown.
(a) Construct the Jacobian matrix J.
(b) Give an analytical expression of the column in the strain-displacement matrix \mathbf{B}(r,s) that corresponds to the displacement u_{1} .
text_image
Axis of revolution y 1 u₁ (-3, 3) 2 (5, 5) (-2, -1) 3 x 4 (3, -3) 7
5.17. The eight-node isoparametric element shown has all its nodal point displacements constrained to zero except for u_{1} . The element is subjected to a concentrated load P into u_{1} .
(a) Calculate and sketch the displacements corresponding to P.
(b) Sketch all element stresses corresponding to the deformed configuration. Use an oblique/aerial view for your sketches.
text_image
s u₁ P 2 r 4 y x Plane stress condition
(unit thickness)
Poisson's ratio \nu = 0.25
Young's modulus E
5.18. The eight-node element assemblage shown is used in a finite element analysis. Calculate the diagonal elements of the stiffness matrix and consistent mass matrix corresponding to the degree of freedom U_{100} .
text_image
6 3 U₁₀₀ All Yo Po Pla Ma
All elements of equal size
Young's modulus E
Poisson's ratio \nu = 0.3
Plane stress analysis, thickness = 0.5
Mass density ρ