김경종
25 KiB
text_image
(a) Coordinate system and coordinate points y s 3 s = 1 6 5 s = 1/2 1 r = 0 4 r = 1/2 2 r = 1 s = 0 r x
(a) Coordinate system and nodal points
Include only if node i is defined
| i=4 | i=5 | i=6 | ||
| h1= | 1-r-s | -1/2h4 | -1/2h6 | |
| h2= | r | -1/2h4 | -1/2h5 | |
| h3= | s | -1/2h5 | -1/2h6 | |
| h4= | 4r(1-r-s) | |||
| h5= | 4rs | |||
| h6= | 4s(1-r-s) |
(b) Interpolation functions
Figure 5.11 Interpolation functions of three to six variable-number-nodes two-dimensional triangle
text_image
t 4 P(r, s, t) 1 3 s 2 r
Figure 5.12 Natural coordinate system of tetrahedral element
text_image
t 4 10 8 1 7 3 s 2 5 6 r
(a) Coordinate system and nodal points
Include only if node i is defined
| i=5 | i=6 | i=7 | i=8 | i=9 | i=10 | ||
| h1= | 1-r-s-t | -1/2h5 | -1/2h7 | -1/2h10 | |||
| h2= | r | -1/2h5 | -1/2h6 | -1/2h8 | |||
| h3= | s | -1/2h6 | -1/2h7 | -1/2h9 | |||
| h4= | t | -1/2h8 | -1/2h9 | -1/2h10 | |||
| h5= | 4r(1-r-s-t) | ||||||
| h6= | 4rs | ||||||
| h7= | 4s(1-r-s-t) | ||||||
| h8= | 4rt | ||||||
| h9= | 4st | ||||||
| h10= | 4t(1-r-s-t) | ||||||
| (b) Interpolation functions | |||||||
Figure 5.13 Interpolation functions of four to ten variable-number-nodes three-dimensional tetrahedral element
functions of a 4 to 10 three-dimensional variable-number-nodes element are given in Fig. 5.13.
To evaluate the element matrices, it is necessary to include the Jacobian transformation as given in (5.24) and to perform the r integrations from 0 to 1, the s integrations from 0 to (1 - r) , and the t integrations from 0 to (1 - r - s) . As for the quadrilateral elements, these integrations are carried out effectively in general analysis using numerical integration, but the integration rules employed are different from those used for quadrilateral elements (see Section 5.5.4).
EXAMPLE 5.19: The triangular element shown in Fig. E5.19 is subjected to the body force vector f^{B} per unit volume. Calculate the consistent nodal point load vector.
text_image
3 4 cm 6 4 cm 1 3 cm 5 cm 5 cm 4 3 cm 2 f^B = [f^B_x / f^B_y] = [10/20] N/cm^3 Element thickness t = 2.0 cm y x
Figure E5.19 Six-node triangular element
Let us use the displacement vector
\hat {\mathbf {u}} ^ {T} = \left[ \begin{array}{l l l l l l l l} u _ {1} & v _ {1} & u _ {2} & v _ {2} & u _ {3} & v _ {3} & \dots & v _ {6} \end{array} \right]
Hence, the load vector corresponding to the applied body force loading is
\mathbf {R} _ {B} = \int_ {V} \left[ \begin{array}{c c} h _ {1} & f _ {x} ^ {B} \\ h _ {1} & f _ {y} ^ {B} \\ \vdots \\ h _ {6} & f _ {y} ^ {B} \end{array} \right] d V
We have \mathbf{J} = \begin{bmatrix} 6 & 0 \\ 0 & 8 \end{bmatrix}
and note that the Jacobian matrix is diagonal and constant and that det J is equal to twice the area of the triangle. The integrations involve the following term for node i:
f _ {i} = \int_ {r = 0} ^ {1} \int_ {s = 0} ^ {1 - r} h _ {i} t \det \mathbf {J} d s d r
which gives f_{1} = f_{2} = f_{3} = 0 , whereas f_{4} = f_{5} = f_{6} = (t / 6) det \mathbf{J} . Hence, we obtain
\mathbf {R} _ {B} ^ {T} = \left[ \begin{array}{l l l l l l l l l} 0 & \dots & 0 & 1 6 0 & 3 2 0 & 1 6 0 & 3 2 0 & 1 6 0 & 3 2 0 \end{array} \right] \tag {a}
with the consistent nodal forces at the corner nodes equal to zero. Note that the total applied load is of course statically equivalent to the nodal forces listed in (a).
5.3.3 Convergence Considerations
We discussed in Section 4.3 the requirements for monotonic convergence of a finite element discretization. Since isoparametric elements are used very widely, let us address some key issues of convergence specifically for these elements.
Basic Requirements for Convergence
The two requirements for monotonic convergence are that the elements (or the mesh) must be compatible and complete.
To investigate the compatibility of an element assemblage, we need to consider each edge, or rather face, between adjacent elements. For compatibility it is necessary that the coordinates and the displacements of the elements at the common face be the same. We note that for the elements considered here, the coordinates and displacements on an element face are determined only by nodes and nodal degrees of freedom on that face. Therefore, compatibility is satisfied if the elements have the same nodes on the common face and the coordinates and displacements on the common face are in each element defined by the same interpolation functions.
Examples of adjacent elements that do and do not preserve compatibility are shown in Fig. 5.14.
Figure 5.14 Compatible and incompatible two-dimensional element assemblage
In practice, mesh grading is frequently necessary (see Section 4.3), and the isoparametric elements show particular flexibility in achieving compatible graded meshes (see Fig. 5.15).
Completeness requires that the rigid body displacements and constant strain states be possible. One way to investigate whether these criteria are satisfied for an isoparametric element is to follow the considerations given in Section 4.3.2. However, we now want to obtain more insight into the specific conditions that pertain to the isoparametric formulation of a continuum element. For this purpose we consider in the following discussion a three-dimensional continuum element because the one- and two-dimensional elements can be regarded as special cases of these three-dimensional considerations. For the rigid body and constant strain states to be possible, the following displacements defined in the local element coordinate system must be contained in the isoparametric formulation
\left. \begin{array}{l} u = a _ {1} + b _ {1} x + c _ {1} y + d _ {1} z \\ v = a _ {2} + b _ {2} x + c _ {2} y + d _ {2} z \\ w = a _ {3} + b _ {3} x + c _ {3} y + d _ {3} z \end{array} \right\} \tag {5.48}
text_image
4-node element 5-node element 8-node element 4-node element 4-node A element 4-node element C B L L C uC B vB uB A vA uA Constraint uA = (uC + uB) / 2 equations: vA = (vC + vB) / 2
(a) 4-node to 8-node element transition region
(b) 4-node to 4-node element transition; from one to two layers
natural_image
Pure geometric grid pattern with dots and lines, no text or symbols present
(c) 9-node to 9-node element transition region; from one to two layers
Figure 5.15 Some transitions with compatible element layouts
where the a_{j}, b_{j}, c_{j} , and d_{j}, j = 1, 2, 3 , are constants. The nodal point displacements corresponding to this displacement field are
\left. \begin{array}{l} u _ {i} = a _ {1} + b _ {1} x _ {i} + c _ {1} y _ {i} + d _ {1} z _ {i} \\ v _ {i} = a _ {2} + b _ {2} x _ {i} + c _ {2} y _ {i} + d _ {2} z _ {i} \\ w _ {i} = a _ {3} + b _ {3} x _ {i} + c _ {3} y _ {i} + d _ {3} z _ {i} \end{array} \right\} \tag {5.49}
where i = 1, \ldots, q and q = \text{number of nodes} .
The test for completeness is now as follows: show that the displacements in (5.48) are in fact obtained within the element when the element nodal point displacements are given by (5.49). In other words, we should find that with the nodal point displacements in (5.49), the displacements within the element are actually those given in (5.48).
In the isoparametric formulation we have the displacement interpolation
u = \sum_ {i = 1} ^ {q} h _ {i} u _ {i}; \quad v = \sum_ {i = 1} ^ {q} h _ {i} v _ {i}; \quad w = \sum_ {i = 1} ^ {q} h _ {i} w _ {i}
which, using (5.49), reduces to
\left. \begin{array}{l} u = a _ {1} \sum_ {i = 1} ^ {q} h _ {i} + b _ {1} \sum_ {i = 1} ^ {q} h _ {i} x _ {i} + c _ {1} \sum_ {i = 1} ^ {q} h _ {i} y _ {i} + d _ {1} \sum_ {i = 1} ^ {q} h _ {i} z _ {i} \\ v = a _ {2} \sum_ {i = 1} ^ {q} h _ {i} + b _ {2} \sum_ {i = 1} ^ {q} h _ {i} x _ {i} + c _ {2} \sum_ {i = 1} ^ {q} h _ {i} y _ {i} + d _ {2} \sum_ {i = 1} ^ {q} h _ {i} z _ {i} \\ w = a _ {3} \sum_ {i = 1} ^ {q} h _ {i} + b _ {3} \sum_ {i = 1} ^ {q} h _ {i} x _ {i} + c _ {3} \sum_ {i = 1} ^ {q} h _ {i} y _ {i} + d _ {3} \sum_ {i = 1} ^ {q} h _ {i} z _ {i} \end{array} \right\} \tag {5.50}
Since in the isoparametric formulation the coordinates are interpolated in the same way as the displacements, we can use (5.18) to obtain from (5.50),
\left. \begin{array}{l} u = a _ {1} \sum_ {i = 1} ^ {q} h _ {i} + b _ {1} x + c _ {1} y + d _ {1} z \\ v = a _ {2} \sum_ {i = 1} ^ {q} h _ {i} + b _ {2} x + c _ {2} y + d _ {2} z \\ w = a _ {3} \sum_ {i = 1} ^ {q} h _ {i} + b _ {3} x + c _ {3} y + d _ {3} z \end{array} \right\} \tag {5.51}
The displacements defined in (5.51), however, are the same as those given in (5.48), provided that for any point in the element,
\sum_ {i = 1} ^ {q} h _ {i} = 1 \tag {5.52}
The relation in (5.52) is the condition on the interpolation functions for the completeness requirements to be satisfied. We may note that (5.52) is certainly satisfied at the nodes of an element because the interpolation function h_{i} has been constructed to be unity at node i with all other interpolation functions h_{j}, j \neq i , being zero at that node; but in order that an isoparametric element be properly constructed, the condition must be satisfied for all points in the element.
In the preceding discussion, we considered a three-dimensional continuum element, but the conclusions are also directly applicable to the other isoparametric continuum element formulations. For the one- or two-dimensional continuum elements we simply include only the appropriate displacement and coordinate interpolations in the relations (5.48) to (5.52). We demonstrate the convergence considerations in the following example.
EXAMPLE 5.20: Investigate whether the requirements for monotonic convergence are satisfied for the variable-number-nodes elements in Figs. 5.4 and 5.5.
Compatibility is maintained between element edges in two-dimensional analysis and element faces in three-dimensional analysis, provided that the same number of nodes is used on connecting element edges and faces. A typical compatible element layout is shown in Fig. 5.14(a).
The second requirement for monotonic convergence is the completeness condition. Considering first the basic four-node two-dimensional element, we recognize that the arguments leading to the condition in (5.52) are directly applicable, provided that only the x and y coordinates and u and v displacements are considered.
Evaluating \Sigma_{i=1}^{4} h_i for the element, we find
\frac {1}{4} (1 + r) (1 + s) + \frac {1}{4} (1 - r) (1 + s) + \frac {1}{4} (1 - r) (1 - s) + \frac {1}{4} (1 + r) (1 - s) = 1
Hence, the basic four-node element is complete. We now study the interpolation functions given in Fig. 5.4 for the variable-number-nodes element and find that the total contribution that is added to the basic four-node interpolation functions is always zero for whichever additional node is included. Hence, any one of the possible elements defined by the variable number of nodes in Fig. 5.4 is complete. The proof for the three-dimensional elements in Fig. 5.5 is carried out in an analogous manner.
It follows therefore that the variable-number-nodes continuum elements satisfy the requirements for monotonic convergence.
Order of Convergence, the Effect of Element Distortions
The basic requirements for monotonic convergence, namely, compatibility and completeness, are satisfied by the isoparametric elements, as discussed above, when these elements are of general (but admissible) geometric shape. Therefore, the elements always have the capability to represent the rigid body modes and constant strain states, and convergence is guaranteed.
We discussed in Section 4.3.5 the rates of convergence of sequences of finite element discretizations with the assumptions that the elements are based on polynomial expansions and that uniform meshes of elements with characteristic dimension h are used. For the discussion we used the Pascal triangle to display which polynomial terms are present in various elements. The complete polynomial of highest order in the Pascal triangle determines the order of convergence. Let this degree (now for r, s, t) be k. Then if the exact solution u is sufficiently smooth and uniform meshes are used, the rate of convergence of the finite element solution u_{h} is given by [see (4.102)]
\left\| \mathbf {u} - \mathbf {u} _ {h} \right\| _ {1} \leq c h ^ {k} \tag {5.53}
where k is the order of convergence. The constant c is independent of h but depends on the exact solution of the mathematical model and the material properties.
In general practical finite element analysis, the exact solution of the mathematical model is not smooth (e.g., because of rapid load variations, changes in material properties), and with uniform meshes the order of convergence is much reduced. Therefore, mesh grading must be employed with fine discretizations in regions of nonsmooth stress distributions and coarse discretizations in the other regions. The meshes will therefore be nonuniform and based on geometrically distorted elements using, for example, in two-dimensional analysis general quadrilateral and triangular elements; see Fig. 5.16 for an example of a nine-node quadrilateral element.
The aim in such mesh constructions is then to use meshes in which the density of solution error is (nearly) constant over the domain considered and to use regular meshes. ^{3} When regular meshes are used, the rate of convergence is still given by a form such as (5.53) [see (4.101c)], namely,
\| \mathbf {u} - \mathbf {u} _ {h} \| _ {1} ^ {2} \leq c \sum_ {m} h _ {m} ^ {2 k} \| \mathbf {u} \| _ {k + 1, m} ^ {2}
where h_m denotes the largest dimension of element m (see Fig. 5.16). We note that, in essence, in this relation the interpolation errors over all elements are added to obtain the total interpolation error, which then gives us the usual bound on the actual error of the solution.
The (nearly) constant density of solution error can of course be achieved, in general, only by proper mesh grading and adaptive mesh refinement because the mesh to be used depends on the exact (and unknown) solution. In practice, a refinement of a mesh is constructed on the basis of local error estimates computed from the solution just obtained (with a coarser mesh).
text_image
s hₘ r
(a) Undistorted configuration used in uniform meshes
text_image
s h_m r ρ_m
(b) Aspect-ratio distortion
text_image
s h_m ρ_m r
(c) Parallelogram distortion
text_image
h_l s h_s h_m ρ_m r
(d) Angular distortion
text_image
s Δ hm r ρm
(e) Curved-edge distortion
text_image
s Δ r
(f) Midnode distortion
Figure 5.16 Classification of element distortions for the 9-node two-dimensional element; all midside and interior nodes are placed at their “natural” positions for cases (a) to (e). The value of \Delta should be smaller than h_{m}^{2} for (5.53) to be applicable. An actual distortion in practice would be a combination of those shown.
To introduce a measure of mesh regularity, the element geometric parameter \sigma_{m} , is used,
\sigma_ {m} = \frac {h _ {m}}{\rho_ {m}}
where h_{m} is the largest dimension and \rho_{m} is the diameter of the largest circle (or sphere) that can be inscribed in element m (see Fig. 5.16). A sequence of meshes is regular if \sigma_{m} \leq \sigma_{0} for all elements m and meshes used, where \sigma_{0} is a fixed positive value. In addition, when using meshes of quadrilateral elements in two-dimensional analysis and hexahedral ele-
ments in three-dimensional analysis, we also require that for each element the ratio of the largest to the smallest side lengths ( h_{i}/h_{s} in Fig. 5.16) is smaller than a reasonable positive number. These conditions prevent excessive aspect ratios and geometric distortions of the elements. Referring to Fig. 5.16, the elements in (b), (c), and in particularly (d) are used extensively in regular meshes. ^{4}
The above described mesh grading can in general be achieved with straight-sided elements, and the noncorner nodes can usually be placed at their natural positions (i.e., at the physical x, y, z locations in proportion to the r, s, t distances from the corner nodes); and the most typically used element in Fig. 5.16 is the quadrilateral in Fig. 5.16(d). However, when curved boundaries need to be modeled, the element sides will also be curved [see Fig. 5.16(e)], and we must ask what effect all these geometric distortions will have on the order of convergence.
Whereas the cases in Figs. 5.16(a) to (e) are used extensively in mesh designs, we note that the element distortion shown in Fig. 5.16(f) is avoided, unless specific stress effects need to be modeled, such as in fracture mechanics [where even larger distortions than those shown in Fig. 5.16(f) are used; see Fig. 5.9]. However, the distortion in Fig. 5.16(f) may also arise in geometrically nonlinear analysis.
P. G. Ciarlet and P. -A. Raviart [A] and P. G. Ciarlet [A] have shown in their mathematical analyses that the order of convergence of a sequence of regular meshes with straight element sides is still given as in (5.53) (even though, for example, in two-dimensional analysis, general straight-sided quadrilaterals are used instead of square elements) and that the order of convergence is also still given as in (5.53) with curved element sides and when the noncorner nodes are not placed at their natural positions provided these distortions are small, measured on the size of the element. For the element in Fig. 5.16 the distortions must be o(h^{2}) . The element distortions due to curved sides and due to interior nodes not placed at their natural positions must therefore be small, and in the refinement process the distortions must decrease much faster than the element size. The order of convergence in (5.53) is reached directly when the exact solution u is smooth, whereas when the exact solution is nonsmooth, mesh grading is necessary (to fulfill the requirement that the density of solution error be (almost) constant over the solution domain). We present some solutions to illustrate a few of these results in Section 5.5.5 (see Fig. 5.39). ^{5}
Of course, the actual accuracy attained with a given mesh is also determined by the constant c in (5.53). This constant depends on the specific elements used (all with the complete polynomial of degree k) and also on the geometric distortions of the elements. We should note that if the constant is large, the order of convergence may be of little interest because the h^{k} term may decrease the error sufficiently only at very small values of h.
These remarks pertain to the order of convergence reached when element sizes are small. However, interesting observations also pertain to a study of the predictive capability of elements when the element sizes are large. Namely, element geometric distortions can affect the general predictive capabilities to a significant degree.
^{4} In addition, we can also define a sequence of meshes that is quasi-uniform. In such sequence we also have, in addition to the regularity condition that the ratio of the maximum h_{m} encountered in the mesh over the minimum h_{m} encountered in that same mesh remains for all meshes below a reasonable positive number. Hence, whereas regularity allows that the ratio of element sizes becomes any value, quasi-uniformity restricts the relative sizes that are permitted. Therefore, the error measure in (4.101c) is also valid when a quasi-uniform sequence of meshes is used.
^{5} These solutions are given in Section 5.5.5 because the element matrices of the distorted elements are evaluated using numerical integration and the effect of the numerical integration error must also be considered.
text_image
100 10 M E = 200,000 v = 0.30 Unit thickness y x
text_image
-60
Analytical solution
\left(\tau_ {x x}\right) _ {\max} = 6 0
\tau_ {y y} = \tau_ {x y} = 0
text_image
Same stress distribution as analytical solution
other
| Dimension | Value | | --------- | ------ | | τ_xx | -4.0 | | τ_yy | 22.4 | | τ_xy | 2.0 | | τ_xx | -41.8 | | τ_yy | -38.3 | | τ_xy | -21.8 | | τ_xx | 28.5 | | τ_yy | 46.8 | | τ_xy | -7.2 | | τ_xx | 10.0 | | τ_yy | -13.6 | | τ_xy | -6.1 |
text_image
Same stress distribution as analytical solution
text_image
75 25 Same stress distribution as analytical solution
Figure 5.17 Example demonstrating effect of element distortions on predictive capability of elements
As an example to demonstrate a possible loss of predictive capability when an isoparametric element is geometrically distorted, consider the results given in Fig. 5.17. The single undistorted eight-node element gives the exact (beam theory) solution for the beam bending problem. However, when two elements of distorted shape are used, significant solution errors are obtained. On the other hand, when the same problem is analyzed with nine-node elements, the mesh of two distorted elements gives the correct result.
This example shows that in coarse meshes the stress predictive capability of certain elements can be significantly affected by element geometric distortions. Since, in practice, frequently rather coarse meshes are used and complete convergence studies are not performed, surely it is preferable to use those elements that are least sensitive to element geometric distortions.