김경종
27 KiB
text_image
(-1, +1) (+1, +1) E = 200,000 ν = 0.30 x, u (-1, -1) (+1, -1)
N elements per side, N = 2, 4, 8,...
(a) Square domain considered
\begin{array}{l} u = c _ {1} (1 - x ^ {2}) (1 - y ^ {2}) e ^ {k y} \cos k x \\ v = c _ {1} (1 - x ^ {2}) (1 - y ^ {2}) e ^ {k y} \sin k x \\ c _ {1} = \text { constant }; k = 5 \\ \end{array}
(b) Exact in-plane displacements
Obtain the finite element solution for the body loads f_{x}^{B} and f_{y}^{B} , where
f _ {x} ^ {B} = - \left(\frac {\partial \tau_ {x x}}{\partial x} + \frac {\partial \tau_ {x y}}{\partial y}\right)
f _ {\gamma} ^ {B} = - \left(\frac {\partial \tau_ {\gamma \gamma}}{\partial \gamma} + \frac {\partial \tau_ {\gamma x}}{\partial x}\right)
and \tau_{xx}, \tau_{yy}, \tau_{xy} are the stresses corresponding to the exact in-plane displacements given in (b).
(c) Test problem
Figure 4.12 Ad-hoc test problem for plane stress (or plane strain, axisymmetric) elements. We use, for h small, E - E_{h} = c h^{\alpha} and hence \log (E - E_{h}) = \log c + \alpha \log h (see also (4.101)). The numerical solutions give \alpha = 3.91 .
4.3.5 Rate of Convergence
In the previous sections we considered the conditions required for monotonic convergence of the finite element analysis results and discussed how in general convergence is reached, but we did not mention the rate at which convergence occurs.
As must be expected, the rate of convergence depends on the order of the polynomials used in the displacement assumptions. In this context the notion of complete polynomials is useful.
Figure 4.13 shows the polynomial terms that should be included to have complete polynomials in x and y for two-dimensional analysis. It is seen that all possible terms of the form x^{\alpha}y^{\beta} are present, where \alpha + \beta = k and k is the degree to which the polynomial is complete. For example, we may note that the element investigated in Example 4.6 uses a
line
| log10 h | log10 (E - Eh) |
|---|---|
| -2 | 0 |
| -1 | 2 |
| 0 | 5 |
(d) Solution for plane stress problem
Figure 4.12 (continued)
polynomial displacement that is complete to degree 1 only. Figure 4.13 also shows important notation for polynomial spaces. The spaces P_{k} correspond to the complete polynomials up to degree k. They can also be thought of as the basis functions of triangular elements: the functions in P_{1} correspond to the functions of the linear displacement (constant strain) triangle (see Example 4.17); the functions in P_{2} correspond to the functions of the parabolic displacement (linear strain) triangle (see Section 5.3.2); and so on.
In addition, Fig. 4.13 shows the polynomial spaces Q_{k}, k = 1, 2, 3 , which correspond to the 4-node, 9-node, and 16-node elements, referred to as Lagrangian elements because the displacement functions of these elements are Lagrangian functions (see also Section 5.5.1).
In considering three-dimensional analysis of course a figure analogous to Fig. 4.13 could be drawn in which the variable z would be included.
Let us think about a sequence of uniform meshes idealizing the complete volume of the body being considered. A mesh of a sequence of uniform meshes consists of elements of equal size—square elements when the polynomial spaces Q_{k} are used. Hence, the parameter h can be taken to be a typical length of an element side. The sequence is obtained by taking a starting mesh of elements and subdividing each element with a natural pattern to obtain the next mesh, and then repeating this process. We did this in solving the ad hoc test problem in Fig. 4.12. However, considering an additional analysis problem, for example, the problem in Example 4.6, we would in Fig. 4.11 subdivide each four-node element into four equal new four-node elements to obtain the first refined mesh; then we would
flowchart
Mathematical diagram illustrating polynomial degree and corresponding space for polynomial degrees 1 to 5, with labeled nodes and lines.
Figure 4.13 Polynomial terms in two-dimensional analysis, Pascal triangle
subdivide each element of the first refined mesh into four equal new four-node elements to obtain the second refined mesh; and so on. The continuation of this subdivision process would give the complete sequence of meshes.
To obtain an expression for the rate of convergence, we would ideally use a formula giving d(\mathbf{u}, V_{h}) in (4.95) as a function of h. However, such a formula is difficult to obtain, and it is more convenient to use interpolation theory and work with an upper bound on d(\mathbf{u}, V_{h}) .
Let us assume that we employ elements with complete polynomials of degree k and that the exact solution u to our elasticity problem is “smooth” in the sense that the solution satisfies the relation ^{14}
\begin{array}{l} \| \mathbf {u} \| _ {k + 1} = \left\{\int_ {\text { Vol }} \left[ \sum_ {i = 1} ^ {3} (u _ {i}) ^ {2} + \sum_ {i = 1} ^ {3} \sum_ {j = 1} ^ {3} \left(\frac {\partial u _ {i}}{\partial x _ {j}}\right) ^ {2} \right. \right. \\ \left. + \sum_ {n = 2} ^ {k + 1} \sum_ {i = 1} ^ {3} \sum_ {r + s + t = n} \left(\frac {\partial^ {n} u _ {i}}{\partial x _ {1} ^ {r} \partial x _ {2} ^ {s} \partial x _ {3} ^ {t}}\right) ^ {2} \right] d \text { Vol } \Bigg \} ^ {1 / 2} < \infty \tag {4.99} \\ \end{array}
where of course k \geq 1 .
Therefore, we assume that all derivatives of the exact solution up to order (k + 1) in (4.99) can be calculated.
A basic result of interpolation theory is that there exists an interpolation function \mathbf{u}_I \in V_h such that
\| \mathbf {u} - \mathbf {u} _ {j} \| _ {1} \leq \hat {c} h ^ {k} \| \mathbf {u} \| _ {k + 1} \tag {4.100}
where h is the mesh size parameter indicating the “size” of the elements and \hat{c} is a constant independent of h. Typically, h is taken to be the length of the side of a generic element or the diameter of a circle encompassing that element. Note that u_{l} is not the finite element solution in V_{h} but merely an element in V_{h} that geometrically corresponds to a function
close to u. Frequently, as here, we let u_{j} , at the finite element nodes, take the value of the exact solution u.
Using (4.100) and Property 3 discussed in Section 4.3.4 [see (4.91)], we can now show that the rate of convergence of the finite element solution u_{h} to the exact theory of elasticity solution u is given by the error estimate
\left\| \mathbf {u} - \mathbf {u} _ {h} \right\| _ {1} \leq c h ^ {k} \left\| \mathbf {u} \right\| _ {k + 1} \tag {4.101}
where c is a constant independent of h but dependent on the material properties. Namely, using (4.95) and (4.100), we have
\begin{array}{l} \| \mathbf {u} - \mathbf {u} _ {h} \| _ {1} \leq c d (\mathbf {u}, V _ {h}) \\ \leq c \left\| \mathbf {u} - \mathbf {u} _ {I} \right\| _ {1} \tag {4.101a} \\ \leq c \hat {c} h ^ {k} \| \mathbf {u} \| _ {k + 1} \\ \end{array}
which gives (4.101) with a new constant c . For (4.101), we say that the rate of convergence is given by the complete right-hand-side expression, and we say that the order of convergence is k or, equivalently, that we have o(h^k) convergence.
Another way to look at the derivation of (4.101)—which is of course closely related to the previous derivation—is to use (4.79) and (4.91). Then we have
\begin{array}{l} \left\| \mathbf {u} - \mathbf {u} _ {h} \right\| _ {1} \leq \frac {1}{c _ {1}} \left[ a \left(\mathbf {u} - \mathbf {u} _ {h}, \mathbf {u} - \mathbf {u} _ {h}\right) \right] ^ {1 / 2} \\ \leq \frac {1}{c _ {1}} \left[ a \left(\mathbf {u} - \mathbf {u} _ {I}, \mathbf {u} - \mathbf {u} _ {I}\right) \right] ^ {1 / 2} \tag {4.101b} \\ \leq \frac {c _ {2}}{c _ {1}} \| \mathbf {u} - \mathbf {u} _ {I} \| _ {1} \\ \leq c h ^ {k} \| \mathbf {u} \| _ {k + 1} \\ \end{array}
Hence, we see directly that to obtain the rate of convergence, we really only expressed the distance d(\mathbf{u}, V_h) in terms of an upper bound given by (4.100).
In practice, we frequently simply write (4.101) as
\left\| \mathbf {u} - \mathbf {u} _ {h} \right\| _ {1} \leq c h ^ {k} \tag {4.102}
and we now recognize that the constant c used here is independent of h but depends on the solution and the material properties [because c in (4.101a) and c_{2} , c_{1} in (4.101b) depend on the material properties]. This dependence on the material properties is detrimental when (almost) incompressible material conditions are considered because the constant then becomes very large and the order of convergence k results in good accuracy only at very small (impractical) values of h. For this reason we need in that case the property (4.95) with the constant independent of the material properties, and this requirement leads to the condition (4.156) (see Section 4.5).
The constant c also depends on the kind of element used. While we have assumed that the element is based on a complete polynomial of order k, different kinds of elements within that class in general display a different constant c for the same analysis problem (e.g., triangular and quadrilateral elements). Hence, the actual magnitude of the error may be considerably different for a given h, while the order with which the error decreases as the mesh is refined is the same. Clearly, the magnitude of the constant c can be crucial in practical analysis because it largely determines how small h actually has to be in order to reach an acceptable error.
These derivations of course represent theoretical results, and we may question in how far these results are applicable in practice. Experience shows that the theoretical results indeed closely represent the actual convergence behavior of the finite element discretizations being considered. Indeed, to measure the order of convergence, we may simply consider the equal sign in (4.102) to obtain
\log \left(\left\| \mathbf {u} - \mathbf {u} _ {h} \right\| _ {1}\right) = \log c + k \log h \tag {4.103}
Then, if we plot from our computed results the graph of \log\left(\|\mathbf{u}-\mathbf{u}_{h}\|_{1}\right) versus \log h , we find that the resulting curve indeed has the approximate slope k when h is sufficiently small.
Evaluating the Sobolev norm may require considerable effort, and in practice, we may use the equivalence of the energy norm with the 1-norm. Namely, because of (4.79), we see that (4.101) also holds for the energy norm on the left side, and this norm can frequently be evaluated more easily [see (4.97)]. Figure 4.12 shows an application. Note that the error in strain energy can be evaluated simply by subtracting the current strain energy from the strain energy of the limit solution (or, if known, the exact solution) [see (4.90)]. In the solution in Fig. 4.12 we obtained an order of convergence (of the numerical results) of 3.91, which compares very well with the theoretical value of 4 (here k = 2 and the strain energy is the energy norm squared). Further results of convergence for this ad hoc problem are given in Fig. 5.39 (where distorted elements with numerically integrated stiffness matrices are considered).
The relation in (4.101) gives, in essence, an error estimate for the displacement gradient, hence for the strains and stresses, because the primary contribution in the 1-norm will be due to the error in the derivatives of the displacements. We will primarily use (4.101) and (4.102) but also note that the error in the displacements is given by
\left\| \mathbf {u} - \mathbf {u} _ {h} \right\| _ {0} \leq c h ^ {k + 1} \left\| \mathbf {u} \right\| _ {k + 1} \tag {4.104}
Hence, the order of convergence for the displacements is one order higher than for the strains.
These results are intuitively reasonable. Namely, let us think in terms of a Taylor series analysis. Then, since a finite element of “dimension h” with a complete displacement expansion of order k can represent displacement variations up to that order exactly, the local error in representing arbitrary displacements with a uniform mesh should be o(h^{k+1}) . Also, for a C^{m-1} problem the stresses are calculated by differentiating the displacements m times, and therefore the error in the stresses is o(h^{k+1-m}) . For the theory of elasticity problem considered above, m = 1, and hence the relations in (4.101) and (4.104) are what we might expect.
EXAMPLE 4.25: Consider the problem shown in Figure E4.25. Estimate the error of the finite element solution if linear two-node finite elements are used.
text_image
L x → f^B(x)
Constant cross-sectional area A Young's modulus E
(a) Bar subjected to load per unit length f^B(x) = ax
line
| x | u(x) | u_h(x) |
|---|---|---|
| x | 0 | 0 |
| L | ∞ | ∞ |
(b) Solutions (for finite element solution three elements are used)
Figure E4.25 Analysis of bar
The finite element problem in this case is to calculate u_{h} \in V_{h} such that
(E A u _ {h} ^ {\prime}, v _ {h} ^ {\prime}) = (f ^ {B}, v _ {h}) \quad \forall v _ {h} \in V _ {h}
with V_{h} = \left\{v_{h}\mid v_{h}\in L^{2}(\mathrm{Vol}),\frac{\partial v_{h}}{\partial x}\in L^{2}(\mathrm{Vol}),v_{h}\big|_{x = 0} = 0\right\}
To estimate the error we use (4.91) and can directly say for this simple problem
\int_ {0} ^ {L} (u ^ {\prime} - u _ {h} ^ {\prime}) ^ {2} d x \leq \int_ {0} ^ {L} (u ^ {\prime} - u _ {l} ^ {\prime}) ^ {2} d x \tag {a}
where u is the exact solution, u_{n} is the finite element solution, and u_{l} is the interpolant, meaning that u_{l} is considered to be equal to u at the nodal points. Hence, our aim is now to obtain an upper bound on \int_0^L (u' - u_l')^2 dx .
Consider an arbitrary element with end points x_{i} and x_{i+1} in the mesh. Then we can say that for the exact solution u(x) and x_{i} \leq x \leq x_{i+1} ,
u ^ {\prime} (x) = u ^ {\prime} \big | _ {x _ {c}} + (x - x _ {c}) u ^ {\prime \prime} \big | _ {x = \overline {{{x}}}}
where x = x_{c} denotes a chosen point in the element and \overline{x} is also a point in the element. Let us choose an x_{c} where u^{\prime}|_{x_{c}} = u^{\prime} , which can always be done because
u _ {I} (x _ {i}) = u (x _ {i}), u _ {I} (x _ {i + 1}) = u (x _ {i + 1})
Then we have for the element
\left| u ^ {\prime} (x) - u _ {I} ^ {\prime} \right| \leq h \left(\max _ {0 \leq x \leq L} \left| u ^ {\prime \prime} \right|\right) \tag {b}
where we have introduced the largest absolute value of the second derivative of the exact solution to obtain an upper bound.
With (b) we have
\int_ {0} ^ {L} \left(u ^ {\prime} - u _ {I} ^ {\prime}\right) ^ {2} d x \leq L h ^ {2} \left(\max _ {0 \leq x \leq L} \left| u ^ {\prime \prime} \right|\right) ^ {2}
and hence \left(\int_0^L (u' - u_h')^2 dx\right)^{1 / 2}\leq ch (c)
where the constant c depends on A, E, L , and f^B but is independent of h .
We should recognize that this analysis is quite general but assumes that the exact solution is smooth so that its second derivative can be calculated (in this example given by -f^{B}/EA ). Of course the result in (c) is just the error estimate (4.102).
An interesting additional result is that the nodal point displacements of the finite element solution are for two reasons the exact displacements. First, the exact solution at the nodes due to the distributed loading is the same as that due to the equivalent concentrated loading (the "equivalent" loading calculated by the principle of virtual work). Second, the finite element space V_h contains the exact solution corresponding to the equivalent concentrated loading. Of course, this nice result is a special property of the solution of one-dimensional problems and does not exist in general two- or three-dimensional analysis.
In the above convergence study it is assumed that uniform discretizations are used (that, for example, in two-dimensional analysis the elements are square and of equal size) and that the exact solution is smooth. Also, implicitly, the degree of the element polynomial displacement expansions is not varied. In practice, these conditions are generally not encountered, and we need to ask what the consequences might be.
If the solution is not smooth—for example, because of sudden changes in the geometry, in loads, or in material properties or thicknesses—and the uniform mesh subdivision is used, the order of convergence decreases; hence, the exponent of h in (4.102) is not k but a smaller value dependent on the degree of “loss of smoothness.”
In practice of course graded meshes are used in such analyses, with small elements in the areas of high stress variation and larger elements away from these regions. The order of convergence of the solutions is then still given by (4.101) but rewritten as
\left\| \mathbf {u} - \mathbf {u} _ {h} \right\| _ {1} ^ {2} \leq c \sum_ {m} h _ {m} ^ {2 k} \left\| \mathbf {u} \right\| _ {k + 1, m} ^ {2} \tag {4.101c}
where m denotes an individual element and h_{m} is a measure of the size of the element. Hence the total error is now estimated by summing the local contributions in (4.101) from each element. A good grading of elements means that the error density in each element is about the same.
In practice when mesh grading is employed, geometrically distorted elements are invariably used. Hence, for example, general quadrilateral elements are very frequently encountered in two-dimensional analyses. We discuss elements of general geometric shapes in Chapter 5 and point out in Section 5.3.3 that the same orders of convergence are applicable to these elements so long as the magnitude of the geometric distortions is reasonable.
In the above sequence of meshes the same kind of elements are used and the element sizes are uniformly decreased. This approach is referred to as the h-method of analysis. Alternatively, an initial mesh of relatively large and low-order elements may be chosen, and then the polynomial displacement expansions in the elements may be successively increased. For example, a mesh of elements with a bilinear displacement assumption may be used (here k = 1), and then the degree of the polynomial expansion is increased to order 2, 3, . . . p, where p may be 10 or even higher. This approach is referred to as the p-method of analysis. To achieve this increase in element polynomial order efficiently, special interpolation functions have been proposed that allow the calculation of the element stiffness matrix corresponding to a higher interpolation by using the previously calculated stiffness matrix and simply amending this matrix, and that have valuable orthogonality properties (see B. Szabó and I. Babuška [A]). However, unfortunately, these functions lack the internal element displacement variations which are important when elements are geometrically distorted (see K. Kato, N. S. Lee, and K. J. Bathe [A] and Section 5.3.3). We demonstrate the use of these functions in the following example.
EXAMPLE 4.26: Consider the one-dimensional bar element shown in Fig. E4.26. Let (\mathbf{K}_{p}) be the stiffness matrix corresponding to the order of displacement interpolation p, where p = 1, 2, 3, \ldots , and let the interpolation functions corresponding to p = 1 be
h _ {1} = \frac {1}{2} (1 - x); \quad h _ {2} = \frac {1}{2} (1 + x) \tag {a}
text_image
f^B = f(x) 1 x = -1 x 2 x = +1 Young's modulus E Cross-sectional area A
Figure E4.26 Bar element subjected to varying body force
For the higher-order interpolations use
h _ {i} (x) = \phi_ {i - 1} (x) \quad i = 3, 4, \dots \tag {b}
where \phi_j = \frac{1}{[2(2j - 1)]^{1/2}} [P_j(x) - P_{j-2}(x)] (c)
and the P_{j} are the Legendre polynomials (see, for example, E. Kreyszig [A]),
P _ {0} = 1
P _ {1} = x
P _ {2} = \frac {1}{2} (3 x ^ {2} - 1)
P _ {3} = \frac {1}{2} (5 x ^ {3} - 3 x)
P _ {4} = \frac {1}{8} (3 5 x ^ {4} - 3 0 x ^ {2} + 3)
• • •
(n + 1) P _ {n + 1} = (2 n + 1) x P _ {n} - n P _ {n - 1}
Calculate the stiffness matrix (\mathbf{K})_{p} and corresponding load vector of the element for p \geq 1 .
Let us first note that these interpolation functions fulfill the requirements of monotonic convergence: the displacement continuity between elements is enforced, and the functions are complete (they can represent the rigid body mode and the constant strain state). This follows because the functions in (a) fulfill these requirements and the functions in (b) merely add higher-order displacement variations within the element with h_{i}=0 at x=\pm1 , i\geq3 .
The stiffness matrix and load vector of the element are obtained using (4.19) and (4.20). Hence, typical elements of the stiffness matrix and load vector are
K _ {i j} = \int_ {- 1} ^ {+ 1} A E \frac {d h _ {i}}{d x} \frac {d h _ {j}}{d x} d x \tag {d}
R _ {i} ^ {B} = \int_ {- 1} ^ {+ 1} f (x) h _ {i} d x
The evaluation of (d) gives
(\mathbf {K}) _ {p} = \frac {A E}{2} \left[ \begin{array}{c c c c c} 1 & - 1 & & \text { zero } \\ - 1 & 1 & & \text { entries } \\ & & 2 \cdot & & \\ & & & \ddots & 2 \cdot \\ & & & & 2 \end{array} \right] (p + 1) \times (p + 1) \tag {e}
where we note that, in essence, the usual 2 \times 2 stiffness matrix corresponding to the interpolation functions (a) has been amended by diagonal entries corresponding to the internal element displacement modes (b). In this specific case, each such entry is uncoupled from all other entries because of the orthogonality properties of the Legendre functions. Hence, as the order of the element is increased, additional diagonal entries are simply computed and all other stiffness coefficients are unchanged.
This structure of the matrix (\mathbf{K})_{p} makes the solution of the governing equations of an element assemblage simple, and the conditioning of the coefficient matrix is always good irrespective of how high an order of element matrices is used. Note also that if the finite element solution is known for elements with a given order of interpolation, then the solution for an increased order of interpolation within the elements is obtained simply by calculating and adding the additional displacements due to the additional internal element modes.
Since the sets of displacement functions corresponding to the matrix (\mathbf{K})_{p+1} contain the sets of functions corresponding to the matrix (\mathbf{K})_{p} , we refer to the displacement functions and the stiffness matrices as hierarchical functions and matrices. This hierarchical property is generally available when the interpolation order is increased (see Exercise 4.29 and Section 5.2).
The concept given in Example 4.26 is also used to establish the displacement functions for higher-order two- and three-dimensional elements. For example, in the two-dimensional case, the basic functions are h_{i}, i = 1, 2, 3, 4 , used in Example 4.6, and the additional functions are due to side modes and internal modes (see Exercises 4.30 and 4.31).
We noted that in the analysis of a bar structure idealized by elements of the kind discussed in Example 4.26, the coupling between elements is due only to the nodal point displacements with the functions h_{1} and h_{2} , and this leads to the very efficient solution. However, in the two- and three-dimensional cases this computational efficiency is not present because the element side modes couple the displacements of adjacent elements and the governing equations of the finite element assemblage have, in fact, a large bandwidth (see Section 8.2.3).
A very high rate of convergence in the solution of general stress conditions can be obtained if we increase the number of elements and at the same time increase the order of displacement variations in the elements. This approach of mesh/element refinement is referred to as the h/p method and can yield an exponential rate of convergence of the form (see B. Szabó and I. Babuška [A])
\left\| \mathbf {u} - \mathbf {u} _ {h} \right\| _ {1} \leq \frac {c}{\exp [ \beta (N) ^ {\gamma} ]} \tag {4.105}
where c, \beta , and \gamma are positive constants and N is the number of nodes in the mesh. If for comparison with (4.105) we write (4.101) in the same form, we obtain for the h method the algebraic rate of convergence
\| \mathbf {u} - \mathbf {u} _ {h} \| _ {1} \leq \frac {c}{(N) ^ {k / d}}
where d = 1, 2, 3 , respectively, in one-, two-, and three-dimensional problems. The effectiveness of the h/p method lies in that it combines the two attractive properties of the h and p methods: using the p method, an exponential rate of convergence is obtained when the exact solution is smooth, and using the h method, the optimal rate of convergence is maintained by proper mesh grading independent of the smoothness of the exact solution.
While the rate of convergence can be very high in the h/p solution approach, of course, whether the solution procedure is effective depends on the total computational effort expended to reach a specified error (which also depends on the constant c).
A key feature of a finite element solution using the h, p, or h/p methods must therefore be the “proper” mesh grading. The above expressions indicate a priori how convergence to the exact solution will be obtained as the density of elements and the order of interpolations are increased, but the meshes used in the successive solutions must be properly graded. By this we mean that the local error density in each element should be about constant. We discuss the evaluation of errors in the next section.
We also assumed in the above discussion on convergence—considering the linear static model problem—that the finite element matrices are calculated exactly and that the governing equilibrium equations are solved without error. In practice, numerical integration is employed in the evaluation of the element matrices (see Section 5.5), and finite precision arithmetic is used to solve the governing equilibrium equations (see Section 8.2.6); hence some error will clearly be introduced in the solution steps. However, the numerical integration errors will not reduce the order of convergence, provided a reliable integration scheme