Files
김경종 4cc312954f
Tests / Hermetic test suite (push) Has been cancelled
Tests / Skill frontmatter validation (push) Has been cancelled
add wiki
2026-05-28 17:16:48 +09:00

391 lines
27 KiB
Markdown
Raw Permalink Blame History

This file contains ambiguous Unicode characters
This file contains Unicode characters that might be confused with other characters. If you think that this is intentional, you can safely ignore this warning. Use the Escape button to reveal them.
<!-- source-page: 271 -->
of high enough order is used (see Section 5.5.5), and the errors in the solution of equations are normally small unless a very ill-conditioned set of equations is solved (see Section 8.2.6).
# 4.3.6 Calculation of Stresses and the Assessment of Error
We discussed above that for monotonic convergence to the exact results (“exact” within the mechanical, i.e., mathematical, assumptions made) the elements must be complete and compatible. Using compatible (or conforming) elements means that in the finite element representation of a $C^{m-1}$ variational problem, the displacements and their $(m - 1)$ st derivatives are continuous across the element boundaries. Hence, for example, in a plane stress analysis the u and v displacements are continuous, and in the analysis of a plate bending problem using the transverse displacement w as the only unknown variable, this displacement w and its derivatives, $\partial w/\partial x$ and $\partial w/\partial y$ , are continuous. However, this continuity does not mean that the element stresses are continuous across element boundaries.
The element stresses are calculated using derivatives of the displacements [see (4.11) and (4.12)], and the stresses obtained at an element edge (or face) when calculated in adjacent elements may differ substantially if a coarse finite element mesh is used. The stress differences at the element boundaries decrease as the finite element mesh is refined, and the rate at which this decrease occurs is of course determined by the order of the elements in the discretization.
For the same mathematical reason that the element stresses are, in general, not continuous across element boundaries, the element stresses at the surface of the structure that is modeled are, in general, not in equilibrium with the externally applied tractions. However, as for the stress jumps between elements, the difference between the externally applied tractions and the element stresses decreases as the number of elements used to model the structure increases.
The stress jumps across element boundaries and stress imbalances at the boundary of the body are of course a consequence of the fact that stress equilibrium is not accurately satisfied at the differential level unless a very fine finite element discretization is used: we recall the derivation of the principle of virtual work in Example 4.2. The development in this example shows that the differential equations of equilibrium are fulfilled only if the virtual work equation is satisfied for any arbitrary virtual displacements that are zero on the surface of the displacement boundary conditions. In the finite element analysis, the number of “real” and virtual displacement patterns is equal to the number of nodal degrees of freedom, and hence only an approximate solution in terms of satisfying the stress equilibrium at the differential level is obtained (while the compatibility and constitutive conditions are satisfied exactly). The error in the solution can therefore be measured by substituting the finite element solution for the stresses $\tau_{ij}^{h}$ into the basic equations of equilibrium to find that for each geometric domain represented by a finite element,
$$
\tau_ {i j, j} ^ {h} + f _ {i} ^ {B} \neq 0 \tag {4.106}
$$
$$
\tau_ {i j} ^ {h} n _ {j} - t _ {i} \neq 0 \tag {4.107}
$$
where $n_{j}$ represents the direction cosines of the normal to the element domain boundary and the $t_{i}$ are the components of the exact traction vector along that boundary (see Fig. 4.14). Of course, this traction vector of the exact solution is not known, and that the left-hand side of (4.107) is not zero simply shows that we must expect stress jumps between elements.
<!-- source-page: 272 -->
![](images/page-272_4a2efe20732e84dd5561f78b634424cd0b4464302c37ad2e2165295fe6ff7ed8.jpg)
<details>
<summary>text_image</summary>
y
n
t
τij
Domain of
finite element
x
</details>
(a) Exact solution to mathematical model
![](images/page-272_62ec4345e55fa2301a389426828e971eedf23c8f8ef6cd5528596b70ebd5fbd2.jpg)
<details>
<summary>text_image</summary>
y
n
t^h
τ^h_{ij}
t^h = \tau^h n
x
</details>
(b) Finite element solution
Figure 4.14 Finite element representing subdomain of continuum
It can be proven that for low-order elements the imbalance in (4.107) is larger than the imbalance in (4.106), and that for high-order elements the imbalance in (4.106) becomes predominant. In practice, (4.107) can be used to obtain an indication of the accuracy of the stress solution and is easily applied by using the isobands of stresses as proposed by T. Sussman and K. J. Bathe [A]. These isobands are constructed using the calculated stresses without stress smoothing as follows:
Choose a stress measure; typically, pressure or the effective (von Mises) stress is chosen, but of course any stress component may be selected.
Divide the entire range over which the stress measure varies into stress intervals, assign each interval a color (or use black and white shading or simply alternate black and white intervals).
A point in the mesh is given the color of the interval corresponding to the value of the stress measure at that point.
If all stresses are continuous across the element boundaries, then this procedure will yield unbroken isobands of stresses. However, in practice, stress discontinuities arise across the element boundaries, resulting in “breaks” in the bands. The magnitude of the intervals of the stress bands together with the severity of the breaks in the bands indicate directly the magnitude of stress discontinuities (see Fig. 4.15). Hence, the isobands represent an
<!-- source-page: 273 -->
![](images/page-273_f50adb269c07a3d2d1abd178231a4e4164089630129b99e0e2554627370e5933.jpg)
<details>
<summary>text_image</summary>
p = 30 MPa
25
20
15
10
5
(a)
</details>
![](images/page-273_80396a133ed062df8585c78dd5459f22b17097b667d349a18253b45d930128e8.jpg)
<details>
<summary>text_image</summary>
p = 30 MPa
25
20
15
10
5
(b)
</details>
![](images/page-273_e69a8d2e494a03d3fbe7e8dcc4608e26f398783ad63923d8026b46967c1cf1ca.jpg)
<details>
<summary>natural_image</summary>
Diagram of layered geological strata with shaded regions, labeled (c) at bottom (no text or symbols within diagram)
</details>
Figure 4.15 Schematic of estimating stress discontinuities using pressure bands, width of bands = 5 MPa; black and white intervals are used; (a) negligible discontinuities, $\Delta p \ll 5$ MPa; (b) visible discontinuities but bands still distinguishable, $\Delta p \simeq 2$ MPa; (c) visible discontinuities, bands not distinguishable, $\Delta p > 5$ MPa.
"eyeball norm" for the accuracy of the stress prediction $\tau_{ij}^{h}$ achieved with a given finite element mesh.
In linear analysis, the finite element stress values can be calculated using the relation $\tau^{h} = CB\hat{u}$ at any point in the element; however, this evaluation is relatively expensive and hardly possible in general nonlinear analysis (including material nonlinear effects). An adequate approach is to use the integration point values to bilinearly interpolate over the corresponding domain of the element. Figure 4.16 illustrates an example in two-dimensional analysis.
An alternative procedure for obtaining an indication of the error in the calculated stresses $\tau_{ij}^{h}$ is to first find some improved values $(\tau_{ij}^{h})_{\mathrm{impr.}}$ and then evaluate and display
$$
\Delta \tau_ {i j} = \tau_ {i j} ^ {h} - (\tau_ {i j} ^ {h}) _ {\text { impr. }} \tag {4.108}
$$
The display can again be achieved effectively using the isoband procedure discussed above.
Improved values might be found by simply averaging the stress values obtained at the nodes using the procedure indicated in Fig. 4.16 or by using a least squares fit over the integration point values of the elements (see E. Hinton and J. S. Campbell [A]). The least squares procedure might be applied over patches of adjacent elements or even globally over
<!-- source-page: 274 -->
![](images/page-274_61090bccc38bcc17b09197a2656c46983c717f717bbf7a5e5789d0776b47d397.jpg)
<details>
<summary>text_image</summary>
Domain over which stresses are interpolated bilinearly using
the four Gauss point values (3 × 3 Gauss integration is used)
Gauss point
b
b
2a
2a
b = \u221a{3/5} a
(see Section 5.5.3)
</details>
Figure 4.16 Interpolation of stresses from Gauss point stresses
a whole mesh. However, if the domain over which the least squares fit is applied involves many stress points, the solution will be expensive and, in addition, a large error in one part of the domain may affect rather strongly the least squares prediction in the other parts. Here, an important consideration is that when using the direct stress evaluation in (4.12), the stresses can be more accurate at the numerical integration points than at the nodal points, see J. Barlow [A] and J.-F. Hiller and K.J. Bathe [A]. Hence, for a least squares fit, it can be of value to use functions of order higher than that of the stress variations obtained from the assumed displacement functions because in this way improved values can be expected, see for example O. C. Zienkiewicz and J. Z. Zhu [A].
We demonstrate a least squares stress smoothing approach in the following example.
EXAMPLE 4.27: Consider the mesh of nine-node elements shown in Fig. E4.27. Propose reasonable schemes for improving the stress results by nodal point averaging and least squares fitting.
Let $\tau$ be a typical stress component. One simple and frequently effective way of improving the stress results is to bilinearly extrapolate the calculated stress components from the integration points of each element to node $i$ . In this way, for the situation and node $i$ in Fig. E4.27, four values for each stress component are obtained. The mean value, say $(\tau^h)_{\text{mean}}^i$ , of these four values is then taken as the value at nodal point $i$ . After performing similar calculations for each nodal point, the improved value of the stress component over a typical element is
$$
(\tau^ {h}) _ {\mathrm{impr.}} = \sum_ {i = 1} ^ {9} h _ {i} (\tau^ {h}) _ {\mathrm{mean}} ^ {i} \tag {a}
$$
where the $h_{i}$ are the displacement interpolation functions because the averaged nodal values are deemed to be more accurate than the values obtained simply from the derivatives of the displacements (which would imply that an interpolation of one order lower is more appropriate).
The key step in this scheme is the calculation of $(\tau^{h})_{\text{mean}}^{i}$ . Such an improved value can also be extracted by using a procedure based on least squares.
Consider the eight nodes closest to node $i$ , plus node $i$ , and the values of the stress component of interest at the 16 integration points closest to node $i$ (shown in Fig. E4.27). Let $(\tau^{h})_{\mathrm{integr.}}^{i}$ be the known values of the stress component at the integration points, $j = 1, \ldots, 16$ , and let $(\tau^{h})_{\mathrm{nodes}}^{k}$ be the unknown values at the nine nodes (of the domain corresponding to the
<!-- source-page: 275 -->
![](images/page-275_55c9b0260aeb5b9387b83f1755b0bfff5c9967b6f7fb2272b8519e6865d998e3.jpg)
<details>
<summary>text_image</summary>
Integration point
</details>
Figure E4.27 Mesh of nine-node elements. Integration points near node i are also shown.
integration points). We can use the least squares procedure (see Section 3.3.3) to calculate the values $(\tau^{h})_{\text{nodes}}^{k}$ by minimizing the errors between the given integration point values and the values calculated at the same points by interpolation from the nodal point values $(\tau^{h})_{\text{nodes}}^{k}$ ,
$$
\frac {\partial}{\partial \left(\tau^ {h}\right) _ {\text { nodes }} ^ {k}} \left[ \sum_ {j = 1} ^ {1 6} \left(\left(\tau^ {h}\right) _ {\text { integr. }} ^ {j} - \left(\tilde {\tau} ^ {h}\right) _ {\text { integr. }} ^ {j}\right) ^ {2} \right] = 0 \tag {b}
$$
$$
k = 1, \dots , 9
$$
where $(\hat{\tau}^{h})_{\mathrm{integr.}}^{j} = \sum_{k = 1}^{9}h_{k}\Bigg|_{\substack{\text{at integr.}\\ \text{point} j}}(\tau^{h})_{\mathrm{nodes}}^{k} \tag{c}$
Note that in (c) we evaluate the interpolation functions at the 16 integration stations shown in Fig. E4.27. The relations in (b) and (c) give nine equations for the values $(\tau^{h})_{\text{nodes}}^{k}, k = 1, \ldots, 9$ . We solve for these values but accept only the value at node i as the improved stress value, which is now our value for $(\tau^{h})_{\text{mean}}^{i}$ in (a). The same basic procedure is used for all nodes to arrive at nodal “mean” values, so that (a) can be used for all elements.
Of course, we presented in Fig. E4.27 a situation of four equal square nine-node elements. In practical analyses, the elements are generally distorted and fewer or more elements may couple into the node i. Also, element non-corner nodes and special mesh topologies at boundaries need to be considered.
However, in practice, frequently the use of low-order elements is preferred. A procedure to specifically enhance the stresses of low-order three-node elements in two-dimensional solutions and four-node elements in three-dimensional solutions has been presented in D. J. Payen and K. J. Bathe [A]. Here projection equations using the difference of higher-order assumed stresses and the element calculated stresses lead to stress solutions that converge quadratically, instead of just linearly.
Another procedure of improving the stress predictions when using low-order elements is to use interpolation covers. This approach was developed for the solid elements in J.H. Kim and
<!-- source-page: 276 -->
K. J. Bathe [A,B] and for a 3-node shell element in H.M. Jeon, P.S. Lee, and K.J. Bathe [A].
Using (4.108) with the improved stress values gives only an error indicator and, ideally, an actual error measure would be available. Much research effort has been devoted to establish effective error measures of finite element solutions, for an overview see M. Ainsworth and J. T. Oden [A] and T. Grätsch and K. J. Bathe [A]. The basic difficulty is to measure the finite element solution on the exact solution of the mathematical model, which is unknown and can be very complex, considering nonlinear or dynamic analyses. The error measure shall always be close to the exact error over the complete analysis domain, and always be conservative and inexpensive to compute. Such error measure is quite likely out of reach, and instead an error indicator based on (4.108) may need to be sufficient.
# 4.3.7 Exercises
4.25. Calculate the eight smallest eigenvalues of the four-node shell element stiffness matrix available in a finite element program and interpret each eigenvalue and corresponding eigenvector. (Hint: The eigenvalues of the element stiffness matrix can be obtained by carrying out a frequency solution with a mass matrix corresponding to unit masses for each degree of freedom.)
4.26. Show that the strain energy corresponding to the displacement error $e_{h}$ , where $e_{h} = u - u_{h}$ , is equal to the difference in the strain energies, corresponding to the exact displacement solution u and the finite element solution $u_{h}$ .
4.27. Consider the analysis problem in Example 4.6. Use a finite element program to perform the convergence study shown in Fig. 4.12 with the nine-node and four-node (Lagrangian) elements. That is, measure the rate of convergence in the energy norm and compare this rate with the theoretical results given in Section 4.3.5. Use N = 2, 4, 8, 16, 32; consider N = 32 to be the limit solution, and use uniform and graded meshes.
4.28. Perform an analysis of the cantilever problem shown using a finite element program. Use a two-dimensional plane stress element idealization to solve for the static response.
(a) Use meshes of four-node elements.
(b) Use meshes of nine-node elements.
In each case construct a sequence of meshes and identify the rate of convergence of strain energy.
![](images/page-276_49f74a1f860e3c48383c342febf6c80f20339df6e38c5641a91f3feed1326cf2.jpg)
<details>
<summary>text_image</summary>
P
6
</details>
E = 200,000
v = 0.3
![](images/page-276_5a2bf080c68a26f23902754d76e0499349d73051ea8a31f84523466eaa5d36de.jpg)
<details>
<summary>text_image</summary>
P
</details>
<!-- source-page: 277 -->
Also, compare your finite element solutions with the solutions using Bernoulli-Euler and Timoshenko beam theories (see S. H. Crandall, N. C. Dahl, and T. J. Lardner [A] and Section 5.4.1).
4.29. Consider the three-node bar element shown. Construct and plot the displacement functions of the element for the following two cases:
for case 1: $h_i = 1$ at node $i, i = 1, 2, 3$
$$
= 0 \text { at node } j \neq i
$$
for case 2: $h_i = 1$ at node $i, i = 1, 2$
$$
= 0 \text { at node } j \neq i, j = 1, 2
$$
$$
h _ {3} = 1 \text { at node } 3
$$
$$
h _ {3} = 0 \text { at node } 1, 2
$$
We note that the functions for case 1 and case 2 contain the same displacement variations, and hence correspond to the same displacement space. Also, the sets of functions are hierarchical because the three-node element contains the functions of the two-node element.
![](images/page-277_7e028b8ae672e9401023a5b6e17bfd9a2c6085f7b17e2b76c924d820ab661d3e.jpg)
<details>
<summary>text_image</summary>
1
-1
3
x
2
+1
</details>
4.30. Consider the eight-node element shown. Identify the terms of the Pascal triangle present in the element interpolations.
![](images/page-277_dacc444fecf48fa7ea5c70ef5b215e2c47a95806c439d45381c4d88f013800ba.jpg)
<details>
<summary>text_image</summary>
2
y
2
5
1
2
6
8
x
3
7
4
</details>
$$
h _ {1} = \frac {1}{4} (1 + x) (1 + y), h _ {2} = \frac {1}{4} (1 - x) (1 + y)
$$
$$
h _ {3} = \frac {1}{4} (1 - x) (1 - y), h _ {4} = \frac {1}{4} (1 + x) (1 - y)
$$
$$
h _ {5} = \frac {1}{2} (1 + y) \phi_ {2} (x), h _ {6} = \frac {1}{2} (1 - x) \phi_ {2} (y)
$$
$$
h _ {7} = \frac {1}{2} (1 - y) \phi_ {2} (x), h _ {8} = \frac {1}{2} (1 + x) \phi_ {2} (y)
$$
where $\phi_{2}$ is defined in Example 4.26.
<!-- source-page: 278 -->
4.31. A p-element of order p = 4 is obtained by using the following displacement functions.
$h_{i}, i = 1, 2, 3, 4,$ as for the basic four-node element (with corner nodes only; see Example 4.6). $h_{i}, i = 5, \ldots, 16$ to represent side modes.
side 1: $h_i^{(1)} = \frac{1}{2} (1 + y)\phi_j(x);$ $i = 5,9,13;j = 2,3,4$
side 2: $h_i^{(2)} = \frac{1}{2} (1 - x)\phi_j(y);$ $i = 6,10,14;j = 2,3,4$
side 3: $h_i^{(3)} = \frac{1}{2}(1 - y)\phi_j(x);$ $i = 7,11,15;j = 2,3,4$
side 4: $h_i^{(4)} = \frac{1}{2} (1 + x)\phi_j(y);$ $i = 8,12,16;j = 2,3,4$
where $\phi_{2}, \phi_{3}$ , and $\phi_{4}$ have been defined in Example 4.26. $h_{17}$ to represent an internal mode
$$
h _ {1 7} = (1 - x ^ {2}) (1 - y ^ {2})
$$
Identify the terms of the Pascal triangle present in the element interpolations.
![](images/page-278_0da72acfb229891c151acd4ab2f76d59d79c8b0bdeab20997a138e3f86f82f6a.jpg)
<details>
<summary>text_image</summary>
2
2
Side 1
1
2
Side 2
x
Side 4
3
4
y
</details>
4.32. Consider the analysis problem in Example 4.6. Use a finite element program to solve the problem with the meshes of nine-node elements in Exercise 4.27 and plot isobands of the von Mises stress and the pressure (without using stress smoothing). Hence, the isobands will display stress discontinuities between elements. Show how the bands converge to continuous stress bands over the cantilever plate.
# 4.4 INCOMPATIBLE AND MIXED FINITE ELEMENT MODELS
In the previous sections we considered the displacement-based finite element method, and the conditions imposed so far on the assumed displacement (or field) functions were completeness and compatibility. If these conditions are satisfied, the calculated solution converges in the strain energy monotonically (i.e., one-sided) to the exact solution. The completeness condition can, in general, be satisfied with relative ease. The compatibility condition can also be satisfied without major difficulties in $C^{0}$ problems, for example, in
<!-- source-page: 279 -->
plane stress and plane strain problems or in the analysis of three-dimensional solids such as dams. Yet, in the analysis of shell problems, and in complex analyses in which completely different finite elements must be used to idealize different regions of the structure, compatibility may be quite impossible to maintain. However, although the compatibility requirements are violated, experience shows that good results are frequently obtained.
Also, in the search for finite elements it was realized that for shell analysis and the analysis of incompressible media, the pure displacement-based method is not efficient. The difficulties in developing compatible displacement-based finite elements for these problems that are computationally effective, and the realization that by using variational approaches many more finite element discretizations can be developed, led to large research efforts. In these activities various classes of new types of elements have been proposed, and the amount of information available on these elements is voluminous. We shall not present the various formulations in detail but only briefly outline some of the major ideas that have been used and then concentrate upon a formulation for a large class of problems—the analysis of almost incompressible media. The analysis of plate and shell structures using many of the concepts outlined below is then further addressed in Chapter 5.
# 4.4.1 Incompatible Displacement-Based Models
In practice, a frequently made observation is that satisfactory finite element analysis results have been obtained although some continuity requirements between displacement-based elements in the mesh employed were violated. In some instances the nodal point layout was such that interelement continuity was not preserved, and in other cases elements were used that contained interelement incompatibilities (see Example 4.28). The final result was the same in either case, namely, that the displacements or their derivatives between elements were not continuous to the degree necessary to satisfy all compatibility conditions discussed in Section 4.3.2.
Since in finite element analysis using incompatible (nonconforming) elements the requirements presented in Section 4.3.2 are not satisfied, the calculated total potential energy is not necessarily an upper bound to the exact total potential energy of the system, and consequently, monotonic convergence is not ensured. However, having relaxed the objective of monotonic convergence in the analysis, we still need to establish conditions that will ensure at least a nonmonotonic convergence.
Referring to Section 4.3, the element completeness condition must always be satisfied, and it may be noted that this condition is not affected by the size of the finite element. We recall that an element is complete if it can represent the physical rigid body modes (but the element matrix has no spurious zero eigenvalues) and the constant strain states.
However, the compatibility condition can be relaxed somewhat at the expense of not obtaining a monotonically convergent solution, provided that when relaxing this requirement, the essential ingredients of the completeness condition are not lost. We recall that as the finite element mesh is refined (i.e., the size of the elements gets smaller), each element should approach a constant strain condition. Therefore, the second condition on convergence of an assemblage of incompatible finite elements, where the elements may again be of any size, is that the elements together can represent constant strain conditions. We should
<!-- source-page: 280 -->
note that this is not a condition on a single individual element but on an assemblage of elements. That is, although an individual element is able to represent all constant strain states, when the element is used in an assemblage, the incompatibilities between elements may prohibit constant strain states from being represented. We may call this condition the completeness condition on an element assemblage.
As a test to investigate whether an assemblage of nonconforming elements is complete, the patch test has been proposed (see B. M. Irons and A. Razzaque [A]). In this test a specific element is considered and a patch of elements is subjected to the minimum displacement boundary conditions to eliminate all rigid body modes and to the boundary nodal point forces that by an analysis should result in constant stress conditions. If for any patch of elements the element stresses actually represent the constant stress conditions and all nodal point displacements are correctly predicted, we say that the element passes the patch test. Since a patch may also consist of only a single element, this test ensures that the element itself is complete and that the completeness condition is also satisfied by any element assemblage.
The number of constant stress states in a patch test depends on course on the actual number of constant stress states that pertain to the mathematical model; for example, in plane stress analysis three constant stress states must be considered in the patch test, whereas in a fully three-dimensional analysis six constant stress states should be possible.
Fig. 4.17 shows a typical patch of elements used in numerical investigations for various problems. Here of course only one mesh with distorted elements is considered, whereas in fact any patch of distorted elements should be analyzed. This, however, requires an analytical solution. If in practice the element is complete and the specific analyses shown here produce the correct results, then it is quite likely that the element passes the patch test.
When considering displacement-based elements with incompatibilities, if the patch test is passed, convergence is ensured (although convergence may not be monotonic and convergence may be slow).
The patch test is used to assess incompatible finite element meshes, and we may note that when properly formulated displacement-based elements are used in compatible meshes, the patch test is automatically passed.
Figure 4.18(a) shows a patch of eight-node elements (which are discussed in detail in Section 5.2). The tractions corresponding to the plane stress patch test are also shown. The elements form a compatible mesh, and hence the patch test is passed.
However, if we next assign to nodes 1 to 8 individual degrees of freedom for the adjacent elements [e.g., at node 2 we assign two u and v degrees of freedom each for elements 1 and 2] such that the displacements are not tied together at these nodes (and therefore displacement incompatibilities exist along the edges), the patch test is not passed. Figure 4.18(b) gives some results of the solution.
The example in Fig. 4.18(b) uses, in essence, an element that was proposed by E. L. Wilson, R. L. Taylor, W. P. Doherty, and J. Ghaboussi [A]. Since the degrees of freedom of the midside nodes of an element are not connected to the adjacent elements, they can be statically condensed out at the element level (see Section 8.2.4) and a four-node element is obtained. However, as indicated in Fig. 4.18(b), this element does not pass the patch test. In the following example, we consider the element in more detail, first as a square element