김경종
22 KiB
Isoparametric Formulation
Introduction
In this chapter, we introduce the isoparametric formulation of the element sti¤ness matrices. After considering the linear-strain triangular element in Chapter 8, we can see that the development of element matrices and equations expressed in terms of a global coordinate system becomes an enormously di‰cult task (if even possible) except for the simplest of elements such as the constant-strain triangle of Chapter 6. Hence, the isoparametric formulation was developed [1]. The isoparametric method may appear somewhat tedious (and confusing initially), but it will lead to a simple computer program formulation, and it is generally applicable for two- and threedimensional stress analysis and for nonstructural problems. The isoparametric formulation allows elements to be created that are nonrectangular and have curved sides. Furthermore, numerous commercial computer programs (as described in Chapter 1) have adapted this formulation for their various libraries of elements.
We first illustrate the isoparametric formulation to develop the simple bar element sti¤ness matrix. Use of the bar element makes it relatively easy to understand the method because simple expressions result.
We then consider the development of the rectangular plane stress element sti¤- ness matrix in terms of a global-coordinate system that will be convenient for use with the element. These concepts will be useful in understanding some of the procedures used with the isoparametric formulation of the simple quadrilateral element sti¤ness matrix, which we will develop subsequently.
Next, we will introduce numerical integration methods for evaluating the quadrilateral element sti¤ness matrix and illustrate the adaptability of the isoparametric formulation to common numerical integration methods.
Finally, we will consider some higher-order elements and their associated shape functions.
10.1 Isoparametric Formulation of the Bar Element Stiffness Matrix
The term isoparametric is derived from the use of the same shape functions (or interpolation functions) ½N- to define the element’s geometric shape as are used to define the displacements within the element. Thus, when the shape function is u = a _ { 1 } + a _ { 2 } s for the displacement, we use x = a _ { 1 } + a _ { 2 } s for the description of the nodal coordinate of a point on the bar element and, hence, the physical shape of the element.
Isoparametric element equations are formulated using a natural (or intrinsic) coordinate system s that is defined by element geometry and not by the element orientation in the global-coordinate system. In other words, axial coordinate s is attached to the bar and remains directed along the axial length of the bar, regardless of how the bar is oriented in space. There is a relationship (called a transformation mapping) between the natural coordinate system s and the global coordinate system x for each element of a specific structure, and this relationship must be used in the element equation formulations.
We will now develop the isoparametric formulation of the sti¤ness matrix of a simple linear bar element [with two nodes as shown in Figure 10–1(a)].
text_image
x₁ 1 x, u L 2 x₂ (a)
text_image
s = -1 s = 0 | → s s = 1 1 L 2 (b)
Figure 10–1 Linear bar element in (a) a global coordinate system x and (b) a natural coordinate system s
Step 1 Select Element Type
First, the natural coordinate s is attached to the element, with the origin located at the center of the element, as shown in Figure 10–1(b). The s axis need not be parallel to the x axis—this is only for convenience.
We consider the bar element to have two degrees of freedom—axial displacements u _ { 1 } and u _ { 2 } at each node associated with the global x axis.
For the special case when the s and x axes are parallel to each other, the s and x coordinates can be related by
x = x _ {c} + \frac {L}{2} s \tag {10.1.1a}
where x _ { c } is the global coordinate of the element centroid.
Using the global coordinates x _ { 1 } and x _ { 2 } in Eq. (10.1.1a) with x _ { c } = ( x _ { 1 } + x _ { 2 } ) / 2 , we can express the natural coordinate s in terms of the global coordinates as
s = [ x - (x _ {1} + x _ {2}) / 2 ] (2 / (x _ {2} - x _ {1}) \tag {10.1.1b}
line
| s | N₁ |
|---|---|
| -1 | 1 |
| 0 | 0 |
| 1 | 0 |
(a)
line
| s | N₂ |
|---|---|
| -1 | 0 |
| 0 | 1 |
| 1 | 1 |
(b)
text_image
u = N₁u₁ + N₂u₂ Node 1 u₁ s = -1 s = 0 u s = 1 u₂ Node 2
（c）
Figure 10–2 Shape function variations with natural coordinates: (a) shape function N _ { 1 i } , (b) shape function N _ { 2 } , and (c) linear displacement field u plotted over element length
The shape functions used to define a position within the bar are found in a manner similar to that used in Chapter 3 to define displacement within a bar (Section 3.1). We begin by relating the natural coordinate to the global coordinate by
x = a _ {1} + a _ {2} s \tag {10.1.2}
where we note that s is such that - 1 \leqslant s \leqslant 1 . Solving for the \boldsymbol { a _ { i } } ^ { \prime } \mathbf { s } in terms of x _ { 1 } and x _ { 2 } , we obtain
x = \frac {1}{2} [ (1 - s) x _ {1} + (1 + s) x _ {2} ] \tag {10.1.3}
or, in matrix form, we can express Eq. (10.1.3) as
\{x \} = [ N _ {1} \quad N _ {2} ] \left\{ \begin{array}{l} x _ {1} \\ x _ {2} \end{array} \right\} \tag {10.1.4}
where the shape functions in Eq. (10.1.4) are
N _ {1} = \frac {1 - s}{2} \quad N _ {2} = \frac {1 + s}{2} \tag {10.1.5}
The linear shape functions in Eqs. (10.1.5) map the s coordinate of any point in the element to the x coordinate when used in Eq. (10.1.3). For instance, when we substitute s ¼ 1 into Eq. (10.1.3), we obtain x = x _ { 1 } . These shape functions are shown in Figure 10–2, where we can see that they have the same properties as defined for the interpolation functions of Section 3.1. Hence, N _ { 1 } represents the physical shape of the coordinate x when plotted over the length of the element for x _ { 1 } = 1 and x _ { 2 } = 0 , and N _ { 2 } represents the coordinate x when plotted over the length of the element for x _ { 2 } = 1 and x _ { 1 } = 0 . Again, we must have N _ { 1 } + N _ { 2 } = 1 .
These shape functions must also be continuous throughout the element domain and have finite first derivatives within the element.
Step 2 Select a Displacement Function
The displacement function within the bar is now defined by the same shape functions, Eqs. (10.1.5), as are used to define the element shape; that is,
\{u \} = \left[ \begin{array}{l l} N _ {1} & N _ {2} \end{array} \right] \left\{ \begin{array}{l} u _ {1} \\ u _ {2} \end{array} \right\} \tag {10.1.6}
When a particular coordinate s of the point of interest is substituted into [ N ] _ { : } , Eq. (10.1.6) yields the displacement of a point on the bar in terms of the nodal degrees of freedom u _ { 1 } and u2 as shown in Figure 10–2(c). Since u and x are defined by the same shape functions at the same nodes, comparing Eqs. (10.1.4) and (10.1.6), the element is called isoparametric.
Step 3 Define the Strain=Displacement and Stress=Strain Relationships
We now want to formulate element matrix ½B- to evaluate ½k-. We use the isoparametric formulation to illustrate its manipulations. For a simple bar element, no real advantage may appear evident. However, for higher-order elements, the advantage will become clear because relatively simple computer program formulations will result.
To construct the element sti¤ness matrix, we must determine the strain, which is defined in terms of the derivative of the displacement with respect to x. The displacement u, however, is now a function of s as given by Eq. (10.1.6). Therefore, we must apply the chain rule of di¤erentiation to the function u as follows:
\frac {d u}{d s} = \frac {d u}{d x} \frac {d x}{d s} \tag {10.1.7}
We can evaluate ( d u / d s ) and ( d x / d s ) using Eqs. (10.1.6) and (10.1.3). We seek ( d u / d x ) = \varepsilon _ { x } . Therefore, we solve Eq. (10.1.7) for ( d u / d x ) as
\frac {d u}{d x} = \frac {\left(\frac {d u}{d s}\right)}{\left(\frac {d x}{d s}\right)} \tag {10.1.8}
Using Eq. (10.1.6) for u, we obtain
\frac {d u}{d s} = \frac {u _ {2} - u _ {1}}{2} \tag {10.1.9a}
and using Eq. (10.1.3) for x, we have
\frac {d x}{d s} = \frac {x _ {2} - x _ {1}}{2} = \frac {L}{2} \tag {10.1.9b}
because x _ { 2 } - x _ { 1 } = L .
Using Eqs. (10.1.9a) and (10.1.9b) in Eq. (10.1.8), we obtain
\left\{\varepsilon_ {x} \right\} = \left[ - \frac {1}{L} \quad \frac {1}{L} \right] \left\{ \begin{array}{l} u _ {1} \\ u _ {2} \end{array} \right\} \tag {10.1.10}
Since \{ \varepsilon \} = [ B ] \{ d \} , the strain/displacement matrix ½B- is then given in Eq. (10.1.10) as
[ B ] = \left[ \begin{array}{c c} - \frac {1}{L} & \frac {1}{L} \end{array} \right] \tag {10.1.11}
We recall that use of linear shape functions results in a constant B matrix, and hence, in a constant strain within the element. For higher-order elements, such as the quadratic bar with three nodes, ½B- becomes a function of natural coordinate s (see Eq. (10.6.16).
The stress matrix is again given by Hooke’s law as
\underline {{\sigma}} = E \underline {{\varepsilon}} = E \underline {{B}} \underline {{d}}
Step 4 Derive the Element Stiffness Matrix and Equations
The sti¤ness matrix is
[ k ] = \int_ {0} ^ {L} [ B ] ^ {T} [ D ] [ B ] A d x \tag {10.1.12}
However, in general, we must transform the coordinate x to s because [ B ] is, in general, a function of s. This general type of transformation is given by References [4] and [5]
\int_ {0} ^ {L} f (x) d x = \int_ {- 1} ^ {1} f (s) | \underline {{J}} | d s \tag {10.1.13}
where J is called the Jacobian. In the one-dimensional case, we have | \underline { { J } } | = \underline { { J } } . . For the simple bar element, from Eq. (10.1.9b), we have
| \underline {{{J}}} | = \frac {d x}{d s} = \frac {L}{2} \tag {10.1.14}
Observe that in Eq. (10.1.14), the Jacobian relates an element length in the global-coordinate system to an element length in the natural-coordinate system. In general, | \underline { J } | is a function of s and depends on the numerical values of the nodal coordinates. This can be seen by looking at Eq. (10.3.22) for the quadrilateral element. (Section 10.3 further discusses the Jacobian.) Using Eqs. (10.1.13) and (10.1.14) in Eq. (10.1.12), we obtain the sti¤ness matrix in natural coordinates as
[ k ] = \frac {L}{2} \int_ {- 1} ^ {1} [ B ] ^ {T} E [ B ] A d s \tag {10.1.15}
where, for the one-dimensional case, we have used the modulus of elasticity E = [ D ] in Eq. (10.1.15). Substituting Eq. (10.1.11) in Eq. (10.1.15) and performing the simple integration, we obtain
[ k ] = \frac {A E}{L} \left[ \begin{array}{c c} 1 & - 1 \\ - 1 & 1 \end{array} \right] \tag {10.1.16}
which is the same as Eq. (3.1.14). For higher-order one-dimensional elements, the integration in closed form becomes di‰cult if not impossible (see Example 10.7). Even the simple rectangular element sti¤ness matrix is di‰cult to evaluate in closed form (Section 10.3). However, the use of numerical integration, as described in Section 10.4, illustrates the distinct advantage of the isoparametric formulation of the equations.
Body Forces
We will now determine the body-force matrix using the natural coordinate system s. Using Eq. (3.10.20b), the body-force matrix is
\{\hat {f} _ {b} \} = \iiint_ {V} [ N ] ^ {T} \{\hat {X} _ {b} \} d V \tag {10.1.17}
Letting d V = A d x . , we have
\{\hat {f} _ {b} \} = A \int_ {0} ^ {L} [ N ] ^ {T} \{\hat {X} _ {b} \} d x \tag {10.1.18}
Substituting Eqs. (10.1.5) for N _ { 1 } and N _ { 2 } into [ N ] and noting that by Eq. (10.1.9b), d x = ( L / 2 ) ds, we obtain
\{\hat {f} _ {b} \} = A \int_ {- 1} ^ {1} \left\{ \begin{array}{l} \frac {1 - s}{2} \\ \frac {1 + s}{2} \end{array} \right\} \{\hat {X} _ {b} \} \frac {L}{2} d s \tag {10.1.19}
On integrating Eq. (10.1.19), we obtain
\{\hat {f} _ {b} \} = \frac {A L \hat {X} _ {b}}{2} \left\{ \begin{array}{l} 1 \\ 1 \end{array} \right\} \tag {10.1.20}
The physical interpretation of the results for \{ \hat { f } _ { b } \} is that since AL represents the volume of the element and \hat { X } _ { b } the body force per unit volume, then A L \hat { X } _ { b } is the total body force acting on the element. The factor \textstyle { \frac { 1 } { 2 } } indicates that this body force is equally distributed to the two nodes of the element.
Surface Forces
Surface forces can be found using Eq. (3.10.20a) as
\{\hat {f} _ {s} \} = \iint_ {S} \left[ N _ {s} \right] ^ {T} \{\hat {T} _ {x} \} d S \tag {10.1.21}
Assuming the cross section is constant and the traction is uniform over the perimeter and along the length of the element, we obtain
\{\hat {f} _ {s} \} = \int_ {0} ^ {L} \left[ N _ {s} \right] ^ {T} \{\hat {T} _ {x} \} d x \tag {10.1.22}
where we now assume \hat { T } _ { x } is in units of force per unit length. Using the shape functions N _ { 1 } and N _ { 2 } from Eq. (10.1.5) in Eq. (10.1.22), we obtain
\{\hat {f} _ {s} \} = \int_ {- 1} ^ {1} \left\{ \begin{array}{l} \frac {1 - s}{2} \\ \frac {1 + s}{2} \end{array} \right\} \{\hat {T} _ {x} \} \frac {L}{2} d s \tag {10.1.23}
On integrating Eq. (10.1.23), we obtain
\{\hat {f} _ {s} \} = \hat {T} _ {x} \frac {L}{2} \left\{ \begin{array}{l} 1 \\ 1 \end{array} \right\} \tag {10.1.24}
The physical interpretation of Eq. (10.1.24) is that since \hat { T } _ { x } is in force-per-unit-length units, \hat { T } _ { x } L is now the total force. The 1 indicates that the uniform surface traction is equally distributed to the two nodes of the element. Note that if \hat { T } _ { x } were a function of x (or s), then the amounts of force allocated to each node would generally not be equal and would be found through integration as in Example 3.12.
10.2 Rectangular Plane Stress Element
We will now develop the rectangular plane stress element sti¤ness matrix. We will later refer to this element in the isoparametric formulation of a general quadrilateral element.
Two advantages of the rectangular element over the triangular element are ease of data input and simpler interpretation of output stresses. A disadvantage of the rectangular element is that the simple linear-displacement rectangle with its associated straight sides poorly approximates the real boundary condition edges.
The usual steps outlined in Chapter 1 will be followed to obtain the element sti¤- ness matrix and related equations.
Step 1 Select Element Type
Consider the rectangular element shown in Figure 10–3 (all interior angles are 90 ) with corner nodes 1–4 (again labeled counterclockwise) and base and height dimensions 2b and 2h, respectively.
The unknown nodal displacements are now given by
\{d \} = \left\{ \begin{array}{l} u _ {1} \\ v _ {1} \\ u _ {2} \\ v _ {2} \\ u _ {3} \\ v _ {3} \\ u _ {4} \\ v _ {4} \end{array} \right\} \tag {10.2.1}
text_image
v4 u4 4 b b v3 3 u3 h x, u h u1 1 2 u2 v1 v2
Figure 10–3 Basic four-node rectangular element with nodal degrees of freedom
Step 2 Select Displacement Functions
For a compatible displacement field, the element displacement functions u and v must be linear along each edge because only two points (the corner nodes) exist along each edge. We then select the linear displacement functions as
u (x, y) = a _ {1} + a _ {2} x + a _ {3} y + a _ {4} x y \tag {10.2.2}
v (x, y) = a _ {5} + a _ {6} x + a _ {7} y + a _ {8} x y
We can proceed in the usual manner to eliminate the \boldsymbol { a _ { i } } ^ { \prime } \boldsymbol { \mathfrak { s } } s from Eqs. (10.2.2) to obtain
\begin{array}{l} u (x, y) = \frac {1}{4 b h} [ (b - x) (h - y) u _ {1} + (b + x) (h - y) u _ {2} \\ \left. + (b + x) (h + y) u _ {3} + (b - x) (h + y) u _ {4} \right] \tag {10.2.3} \\ \end{array}
v (x, y) = \frac {1}{4 b h} [ (b - x) (h - y) v _ {1} + (b + x) (h - y) v _ {2}
\left. + (b + x) (h + y) v _ {3} + (b - x) (h + y) v _ {4} \right]
These displacement expressions, Eqs. (10.2.3), can be expressed equivalently in terms of the shape functions and unknown nodal displacements as
\{\psi \} = [ N ] \{d \} \tag {10.2.4}
where the shape functions are given by
N _ {1} = \frac {(b - x) (h - y)}{4 b h} \quad N _ {2} = \frac {(b + x) (h - y)}{4 b h} \tag {10.2.5}
N _ {3} = \frac {(b + x) (h + y)}{4 b h} \quad N _ {4} = \frac {(b - x) (h + y)}{4 b h}
and the \boldsymbol { N _ { i } } ^ { \prime } \boldsymbol { \mathrm { s } } are again such that N _ { 1 } = 1 at node 1 and N _ { 1 } = 0 at all the other nodes, with similar requirements for the other shape functions. In expanded form,
Eq. (10.2.4) becomes
\left\{ \begin{array}{l} u \\ v \end{array} \right\} = \left[ \begin{array}{c c c c c c c c} N _ {1} & 0 & N _ {2} & 0 & N _ {3} & 0 & N _ {4} & 0 \\ 0 & N _ {1} & 0 & N _ {2} & 0 & N _ {3} & 0 & N _ {4} \end{array} \right] \left\{ \begin{array}{l} u _ {1} \\ v _ {1} \\ u _ {2} \\ v _ {2} \\ u _ {3} \\ v _ {3} \\ u _ {4} \\ v _ {4} \end{array} \right\} \tag {10.2.6}
Step 3 Define the Strain= Displacement and Stress= Strain Relationships
Again the element strains for the two-dimensional stress state are given by
\left\{ \begin{array}{l} \varepsilon_ {x} \\ \varepsilon_ {y} \\ \gamma_ {x y} \end{array} \right\} = \left\{ \begin{array}{c} \frac {\partial u}{\partial x} \\ \frac {\partial v}{\partial y} \\ \frac {\partial u}{\partial y} + \frac {\partial v}{\partial x} \end{array} \right\} \tag {10.2.7}
Using Eq. (10.2.6) in Eq. (10.2.7) and taking the derivatives of u and v as indicated, we can express the strains in terms of the unknown nodal displacements as
\{\varepsilon \} = [ B ] \{d \} \tag {10.2.8}
where [ B ] = \frac { 1 } { 4 b h } \left[ \begin{array} { c c c c } { - ( h - y ) } & { 0 } & { ( h - y ) } & { 0 } \\ { 0 } & { - ( b - x ) } & { 0 } & { - ( b + x ) } \\ { - ( b - x ) } & { - ( h - y ) } & { - ( b + x ) } & { ( h - y ) } \end{array} \right]
\left. \begin{array}{c c c c} (h + y) & 0 & - (h + y) & 0 \\ 0 & (b + x) & 0 & (b - x) \\ (b + x) & (h + y) & (b - x) & - (h + y) \end{array} \right] \tag {10.2.9}
From Eqs. (10.2.8) and (10.2.9), we observe that \varepsilon _ { x } is a function of y , \varepsilon _ { y } is a function of x , and \gamma _ { x y } is a function of both x and y. The stresses are again given by the formulas in Eq. (6.2.36), where ½B- is now that of Eq. (10.2.9) and \{ d \} is that of Eq. (10.2.1).
Step 4 Derive the Element Stiffness Matrix and Equations
The sti¤ness matrix is determined by
[ k ] = \int_ {- h} ^ {h} \int_ {- b} ^ {b} [ B ] ^ {T} [ D ] [ B ] t d x d y \tag {10.2.10}
with ½D- again given by the usual plane stress or plane strain conditions, Eq. (6.1.8) or (6.1.10). Because the [ B ] matrix is a function of x and y , integration of Eq. (10.2.10) must be performed. The [ k ] matrix for the rectangular element is now of order 8 \times 8 .
The element force matrix is determined by Eq. (6.2.46) as
\{f \} = \iint_ {V} [ N ] ^ {T} \{X \} d V + \{P \} + \iint_ {S} \left[ N _ {s} \right] ^ {T} \{T \} d S \tag {10.2.11}
where ½N- is the rectangular matrix in Eq. (10.2.6), and N _ { 1 } through N _ { 4 } are given by Eqs. (10.2.5). The element equations are then given by
\{f \} = [ k ] \{d \} \tag {10.2.12}
Steps 5–7
Steps 5–7, which involve assembling the global sti¤ness matrix and equations, determining the unknown nodal displacements, and calculating the stress, are identical to those in Section 6.2 for the CST. However, the stresses within each element now vary in both the x and y directions.
As previously pointed out when describing defects for the CST in Chapter 6, the bilinear rectangle element described in this section also cannot provide pure bending. When this element is subjected to pure bending, it also develops false shear strain. This means that in a pure bending deformation, the bending moment needed to produce the deformation is predicted to be larger than the actual value when modeling with the rectangular element. Details of this shear locking are presented in [8].
To avoid the possibility of shear locking that occurs when the rectangular bilinear (four corner noded) element is subjected to bending, the higher order eightnoded quadratic rectangle has been developed and is described briefly in Section 10.6. This eight-noded element is created by adding mid-side nodes to the bilinear element.
10.3 Isoparametric Formulation of the Plane Element Stiffness Matrix
Recall that the term isoparametric is derived from the use of the same shape functions to define the element shape as are used to define the displacements within the element. Thus, when the shape function is u = a _ { 1 } + a _ { 2 } s + a _ { 3 } t + a _ { 4 } s t for the displacement, we use x = a _ { 1 } + a _ { 2 } s + a _ { 3 } t + a _ { 4 } s t for the description of a coordinate point in the plane element.
The natural-coordinate system s-t is defined by element geometry and not by the element orientation in the global-coordinate system x { - } y . Much as in the bar element example, there is a transformation mapping between the two coordinate systems for each element of a specific structure, and this relationship must be used in the element formulation.
We will now discuss the isoparametric formulation of the simple linear plane element sti¤ness matrix. This formulation is general enough to be applied to more