김경종
446 lines
25 KiB
Markdown
446 lines
25 KiB
Markdown
<!-- source-page: 481 -->
|
||
|
||
# 10.4 Gaussian and Newton-Cotes Quadrature (Numerical Integration)
|
||
|
||
In this section, we will describe Gauss’s method, one of the many schemes for numerical evaluation of definite integrals, because it has proved most useful for finite element work.
|
||
|
||
For completion sake, we will also describe the more common numerical integration method of Newton-Cotes. The Newton-Cotes methods for one and two intervals of integration are the well-known trapezoid and Simpson’s one-third rule, respectively. After describing both methods, we will then understand why the Gaussian quadrature method is used in finite element work.
|
||
|
||
# Gaussian Quadrature:
|
||
|
||
To evaluate the integral
|
||
|
||
$$
|
||
I = \int_ {- 1} ^ {1} y d x \tag {10.4.1a}
|
||
$$
|
||
|
||
where $y = y ( x )$ , we might choose (sample or evaluate) y at the midpoint $y ( 0 ) = y _ { 1 }$ and multiply by the length of the interval, as shown in Figure 10–8, to arrive at $I = 2 y _ { 1 }$ , a result that is exact if the curve happens to be a straight line. This is an example of what is called one-point Gaussian quadrature because only one sampling point was used. Therefore,
|
||
|
||
$$
|
||
I = \int_ {- 1} ^ {1} y (x) d x \cong 2 y (0) \tag {10.4.1b}
|
||
$$
|
||
|
||
which is the familiar midpoint rule. Generalization of the formula [Eq. (10.4.1b)] leads to
|
||
|
||
$$
|
||
I = \int_ {- 1} ^ {1} y d x = \sum_ {i = 1} ^ {n} W _ {i} y _ {i} \tag {10.4.2}
|
||
$$
|
||
|
||
That is, to approximate the integral, we evaluate the function at several sampling points n, multiply each value $y _ { i }$ by the appropriate weight $W _ { i } ,$ , and add the terms. Gauss’s method chooses the sampling points so that for a given number of points, the best possible accuracy is obtained. Sampling points are located symmetrically with respect to the center of the interval. Symmetrically paired points are given the
|
||
|
||
|
||
|
||
<details>
|
||
<summary>text_image</summary>
|
||
|
||
Approximate area = 2y₁
|
||
y₁
|
||
-1
|
||
0
|
||
1
|
||
x
|
||
y
|
||
</details>
|
||
|
||
Figure 10–8 Gaussian quadrature using one sampling point
|
||
|
||
<!-- source-page: 482 -->
|
||
|
||
same weight $W _ { i } .$ . Table 10–1 gives appropriate sampling points and weighting coe‰- cients for the first three orders—that is, one, two, or three sampling points (see Reference [2] for more complete tables). For example, using two points (Figure 10–9), we simply have $I = y _ { 1 } + y _ { 2 }$ because $W _ { 1 } = W _ { 2 } = 1 . 0 0 0$ . This is the exact result if $y = f ( x )$ is a polynomial containing terms up to and including $x ^ { 3 }$ . In general, Gaussian quadrature using n points (Gauss points) is exact if the integrand is a polynomial of degree $2 n - 1$ or less. In using n points, we e¤ectively replace the given function $y = f ( x )$ by a polynomial of degree $2 n - 1$ . The accuracy of the numerical integration depends on how well the polynomial fits the given curve.
|
||
|
||
If the function $f ( x )$ is not a polynomial, Gaussian quadrature is inexact, but it becomes more accurate as more Gauss points are used. Also, it is important to understand that the ratio of two polynomials is, in general, not a polynomial; therefore, Gaussian quadrature will not yield exact integration of the ratio.
|
||
|
||
Table 10–1 Table for Gauss points for integration from minus one to one, $\int _ { - 1 } ^ { 1 } y ( x ) d x = \sum _ { i = 1 } ^ { n } W _ { i } y _ { i }$ Wiyi
|
||
|
||
<table><tr><td>Number of Points</td><td>Locations, $x_i$ </td><td>Associated Weights, $W_i$ </td></tr><tr><td>1</td><td> $x_1 = 0.000 \ldots$ </td><td>2.000</td></tr><tr><td>2</td><td> $x_1, x_2 = \pm 0.57735026918962$ </td><td>1.000</td></tr><tr><td>3</td><td> $x_1, x_3 = \pm 0.77459666924148$ </td><td> $\frac{5}{9} = 0.555 \ldots$ </td></tr><tr><td></td><td> $x_2 = 0.000 \ldots$ </td><td> $\frac{8}{9} = 0.888 \ldots$ </td></tr><tr><td>4</td><td> $x_1, x_4 = \pm 0.8611363116$ </td><td>0.3478548451</td></tr><tr><td></td><td> $x_2, x_3 = \pm 0.3399810436$ </td><td>0.6521451549</td></tr></table>
|
||
|
||
|
||
|
||
<details>
|
||
<summary>line</summary>
|
||
|
||
| x | y |
|
||
| ---- | ----- |
|
||
| -1 | y2 |
|
||
| x₁ | y₁ |
|
||
| 1 | y1 |
|
||
</details>
|
||
|
||
Figure 10–9 Gaussian quadrature using two sampling points
|
||
|
||
# Two-Point Formula
|
||
|
||
To illustrate the derivation of a two-point ðn ¼ 2Þ Gauss formula based on Eq. (10.4.2), we have
|
||
|
||
$$
|
||
I = \int_ {- 1} ^ {1} y d x = W _ {1} y _ {1} + W _ {2} y _ {2} = W _ {1} y (x _ {1}) + W _ {2} y (x _ {2}) \tag {10.4.3}
|
||
$$
|
||
|
||
<!-- source-page: 483 -->
|
||
|
||
There are four unknown parameters to determine: $W _ { 1 } , W _ { 2 } , x _ { 1 }$ , and $x _ { 2 }$ . Therefore, we assume a cubic function for y as follows:
|
||
|
||
$$
|
||
y = C _ {0} + C _ {1} x + C _ {2} x ^ {2} + C _ {3} x ^ {3} \tag {10.4.4}
|
||
$$
|
||
|
||
In general, with four parameters in the two-point formula, we would expect the Gauss formula to exactly predict the area under the curve. That is,
|
||
|
||
$$
|
||
A = \int_ {- 1} ^ {1} \left(C _ {0} + C _ {1} x + C _ {2} x ^ {2} + C _ {3} x ^ {3}\right) d x = 2 C _ {0} + \frac {2 C _ {2}}{3} \tag {10.4.5}
|
||
$$
|
||
|
||
However, we will assume, based on Gauss’s method, that $W _ { 1 } = W _ { 2 }$ and $x _ { 1 } = x _ { 2 }$ as we use two symmetrically located Gauss points at $x = \pm a$ with equal weights. The area predicted by Gauss’s formula is
|
||
|
||
$$
|
||
A _ {G} = W y (- a) + W y (a) = 2 W (C _ {0} + C _ {2} a ^ {2}) \tag {10.4.6}
|
||
$$
|
||
|
||
where $y ( - a )$ and $y ( a )$ are evaluated using $\operatorname { E q . }$ (10.4.4). If the error, $e = A - A _ { G } ,$ is to vanish for any $C _ { 0 }$ and $C _ { 2 }$ , we must have, using Eqs. (10.4.5) and (10.4.6) in the error expression,
|
||
|
||
$$
|
||
\frac {\partial e}{\partial C _ {0}} = 0 = 2 - 2 W \quad \text { or } \quad W = 1 \tag {10.4.7a}
|
||
$$
|
||
|
||
and ${ \frac { \hat { c } e } { \hat { o } C _ { 2 } } } = 0 = { \frac { 2 } { 3 } } - 2 a ^ { 2 } W \qquad { \mathrm { o r } } \qquad a = \sqrt { \frac { 1 } { 3 } } = 0 . 5 7 7 3 \dots$ ð10:4:7bÞ
|
||
|
||
Now $W = 1$ and $a = 0 . 5 7 7 3 \dots$ . are the $W _ { i } { ' } \mathrm { s }$ and $\boldsymbol { a } _ { i } ^ { \cdot } \mathbf { s } \left( \boldsymbol { x } _ { i } ^ { \cdot } \mathbf { s } \right)$ for the two-point Gaussian quadrature given in Table 10–1.
|
||
|
||
# Example 10.2
|
||
|
||
Evaluate the integrals ðaÞ $I = \int _ { - 1 } ^ { 1 } [ x ^ { 2 } + \cos ( x / 2 ) ]$ dx and (b) $\begin{array} { r } { I = \int _ { - 1 } ^ { 1 } ( 3 ^ { x } - x ) d x } \end{array}$ using three-point Gaussian quadrature.
|
||
|
||
(a) Using Table 10–1 for the three Gauss points and weights, we have $x _ { 1 } = x _ { 3 } =$ $\pm 0 . 7 7 4 5 9 \ldots , x _ { 2 } = 0 . 0 0 0 \ldots , W _ { 1 } = W _ { 3 } = \frac { 5 } { 9 } ,$ and $\begin{array} { r } { W _ { 2 } = \frac { 8 } { 9 } . } \end{array}$ . The integral then becomes
|
||
|
||
$$
|
||
\begin{array}{l} I = \left[ (- 0. 7 7 4 5 9) ^ {2} + \cos \left(- \frac {0 . 7 7 4 5 9}{2} \mathrm{rad}\right) \right] \frac {5}{9} + \left[ 0 ^ {2} + \cos \frac {0}{2} \right] \frac {8}{9} \\ + \left[ (0. 7 7 4 5 9) ^ {2} + \cos \left(\frac {0 . 7 7 4 5 9}{2} \mathrm{rad}\right) \right] \frac {5}{9} \\ = 1. 9 1 8 + 0. 6 6 7 = 2. 5 8 5 \\ \end{array}
|
||
$$
|
||
|
||
Compared to the exact solution, we have $I _ { \mathrm { e x a c t } } = 2 . 5 8 5$ .
|
||
|
||
In this example, three-point Gaussian quadrature yields the exact answer to four significant figures.
|
||
|
||
<!-- source-page: 484 -->
|
||
|
||
(b) Using Table 10–1 for the three Gauss points and weights as in part (a), the integral then becomes
|
||
|
||
$$
|
||
\begin{array}{l} I = \left[ 3 ^ {(- 0. 7 7 4 5 9)} - (- 0. 7 7 4 5 9) \right] \frac {5}{9} + \left[ 3 ^ {0} - 0 \right] \frac {8}{9} + \left[ 3 ^ {(0. 7 7 4 5 9)} - (0. 7 7 4 5 9) \right] \frac {5}{9} \\ = 0. 6 6 7 5 5 + 0. 8 8 8 8 9 + 0. 8 6 0 6 5 = 2. 4 2 2 9 (2. 4 2 3 \text { to four significant figures }) \\ \end{array}
|
||
$$
|
||
|
||
Compared to the exact solution, we have $I _ { \mathrm { e x a c t } } = 2 . 4 2 7$ . The error is $2 . 4 2 7 - 2 . 4 2 3 =$ 0:004.
|
||
|
||
In two dimensions, we obtain the quadrature formula by integrating first with respect to one coordinate and then with respect to the other as
|
||
|
||
$$
|
||
\begin{array}{l} I = \int_ {- 1} ^ {1} \int_ {- 1} ^ {1} f (s, t) d s d t = \int_ {- 1} ^ {1} \left[ \sum_ {i} W _ {i} f \left(s _ {i}, t\right) \right] d t \\ = \sum_ {j} W _ {j} \left[ \sum_ {i} W _ {i} f (s _ {i}, t _ {j}) \right] = \sum_ {i} \sum_ {j} W _ {i} W _ {j} f (s _ {i}, t _ {j}) \tag {10.4.8} \\ \end{array}
|
||
$$
|
||
|
||
In Eq. (10.4.8), we need not use the same number of Gauss points in each direction (that is, i does not have to equal $j )$ , but this is usually done. Thus, for example, a four-point Gauss rule (often described as a 2
|
||
2 rule) is shown in Figure 10–10. Equation (10.4.8) with i ¼ 1; 2 and j ¼ 1; 2 yields
|
||
|
||
$$
|
||
I = W _ {1} W _ {1} f (s _ {1}, t _ {1}) + W _ {1} W _ {2} f (s _ {1}, t _ {2}) + W _ {2} W _ {1} f (s _ {2}, t _ {1}) + W _ {2} W _ {2} f (s _ {2}, t _ {2}) \tag {10.4.9}
|
||
$$
|
||
|
||
where the four sampling points are at $s _ { i } , t _ { i } = \pm 0 . 5 7 7 3 \ldots = \pm 1 / \sqrt { 3 } ,$ and the weights are all 1.000. Hence, the double summation in Eq. (10.4.8) can really be interpreted as a single summation over the four points for the rectangle.
|
||
|
||
In general, in three dimensions, we have
|
||
|
||
$$
|
||
I = \int_ {- 1} ^ {1} \int_ {- 1} ^ {1} \int_ {- 1} ^ {1} f (s, t, z) d s d t d z = \sum_ {i} \sum_ {j} \sum_ {k} W _ {i} W _ {j} W _ {k} f (s _ {i}, t _ {j}, z _ {k}) \tag {10.4.10}
|
||
$$
|
||
|
||
|
||
|
||
<details>
|
||
<summary>text_image</summary>
|
||
|
||
s = -0.5773 ... (i = 1)
|
||
(s₁, t₂)
|
||
2
|
||
(s₂, t₂)
|
||
4
|
||
t = 0.5773 ... (j = 2)
|
||
(s₁, t₁)
|
||
1
|
||
(s₂, t₁)
|
||
3
|
||
t = -0.5773 ... (j = 1)
|
||
s = 0.5773 ... (i = 2)
|
||
</details>
|
||
|
||
Figure 10–10 Four-point Gaussian quadrature in two dimensions
|
||
|
||
<!-- source-page: 485 -->
|
||
|
||
# Newton-Cotes Numerical Integration:
|
||
|
||
We now describe the common numerical integration method called the Newton-Cotes method for evaluation of definite integrals. However, the method does not yield as accurate of results as the Gaussian quadrature method and so is not normally used in finite element method evaluations, such as to evaluate the sti¤ness matrix.
|
||
|
||
To evaluate the integral
|
||
|
||
$$
|
||
I = \int_ {- 1} ^ {1} y d x
|
||
$$
|
||
|
||
we assume the sampling points of $y ( x )$ are spaced at equal intervals. Since the limits of integration are from 1 to 1 using the isoparametric formulation, the Newton-Cotes formula is given by
|
||
|
||
$$
|
||
I = \int_ {- 1} ^ {1} y d x = h \sum_ {i = 0} ^ {n} C _ {i} y _ {i} = h \left[ C _ {0} y _ {0} + C _ {1} y _ {1} + C _ {2} y _ {2} + C _ {3} y _ {3} + \dots + C _ {n} y _ {n} \right] \tag {10.4.11}
|
||
$$
|
||
|
||
where the $C _ { i }$ are the Newton-Cotes constants for numerical integration with i intervals (the number of intervals will be one less than the number of sampling points, n) and h is the interval between the limits of integration (for limits of integration between 1 and 1 this makes $h = 2 )$ ). The Newton-Cotes constants have been published and are summarized in Table 10–2 for $i = 1$ to 6. The case $i = 1$ corresponds to the wellknown trapezoid rule illustrated by Figure 10–11. The case $i = 2$ corresponds to the well-known Simpson one-third rule. It is shown [9] that the formulas for $i = 3$ and $i = 5$ have the same accuracy as the formulas for $i = 2$ and $i = 4$ , respectively. Therefore, it is recommended that the even formulas with $i = 2$ and $i = 4$ be used in practice. To obtain greater accuracy one can then use a smaller interval (include more evaluations of the function to be integrated). This can be accomplished by using a higher-order Newton-Cotes formula, thus increasing the number of intervals i.
|
||
|
||
It is shown [9] that we need to use n equally spaced sampling points to integrate exactly a polynomial of order at most $n - 1$ . On the other hand, using Gaussian quadrature
|
||
|
||
Table 10–2 Table for Newton-Cotes intervals and points for integration,
|
||
$\int _ { - 1 } ^ { 1 } y ( x ) d x = h \sum _ { i = 0 } ^ { n } C _ { i } y _ { i }$
|
||
<table><tr><td>Intervals, i</td><td>No. of Points, n</td><td> $C_0$ </td><td> $C_1$ </td><td> $C_2$ </td><td> $C_3$ </td><td> $C_4$ </td><td> $C_5$ </td><td> $C_6$ </td></tr><tr><td>1</td><td>2</td><td>1/2</td><td>1/2</td><td></td><td></td><td colspan="3">(trapezoid rule)</td></tr><tr><td>2</td><td>3</td><td>1/6</td><td>4/6</td><td>1/6</td><td></td><td colspan="3">(Simpson's 1/3 rule)</td></tr><tr><td>3</td><td>4</td><td>1/8</td><td>3/8</td><td>3/8</td><td>1/8</td><td colspan="3">(Simpson's 3/8 rule)</td></tr><tr><td>4</td><td>5</td><td>7/90</td><td>32/90</td><td>12/90</td><td>32/90</td><td colspan="3">7/90</td></tr><tr><td>5</td><td>6</td><td>19/288</td><td>75/288</td><td>50/288</td><td>50/288</td><td>75/288</td><td>19/288</td><td></td></tr><tr><td>6</td><td>7</td><td>41/840</td><td>216/840</td><td>27/840</td><td>272/840</td><td>27/840</td><td>216/840</td><td>41/840</td></tr></table>
|
||
|
||
<!-- source-page: 486 -->
|
||
|
||
|
||
|
||
<details>
|
||
<summary>line</summary>
|
||
|
||
| x | y |
|
||
| --- | ----- |
|
||
| -1 | y₀ |
|
||
| 1 | y₁ |
|
||
</details>
|
||
|
||
Figure 10–11 Approximation of numerical integration (approximate area under curve) using i ¼ 1 interval, $n = 2$ sampling points (trapezoid rule), for
|
||
|
||
$$
|
||
I = \int_ {- 1} ^ {1} y (x) d x = h \sum_ {i = 0} ^ {2} C _ {i} y _ {i}
|
||
$$
|
||
|
||
we have previously stated that we use unequally spaced sampling points n and integrate exactly a polynomial of order at most $2 n - 1$ . For instance, using the Newton-Cotes formula with $n = 2$ sampling points, the highest order polynomial we can integrate exactly is a linear one. However, using Gaussian quadrature, we can integrate a cubic polynomial exactly. Gaussian quadrature is then more accurate with fewer sampling points than Newton-Cotes quadrature. This is because Gaussian quadrature is based on optimizing the position of the sampling points (not making them equally spaced as in the Newton-Cotes method) and also optimizing the weights $W _ { i }$ given in Table 10–1. After the function is evaluated at the sampling points, the corresponding weights are multiplied by these evaluated functions as was illustrated in Example 10.2.
|
||
|
||
Example 10.3 is used to illustrate the Newton-Cotes method and compare its accuracy to that of the Gaussian quadrature method previously described.
|
||
|
||
# Example 10.3
|
||
|
||
Solve Example 10.2 using the Newton-Cotes method with $i = 2$ intervals $( n = 3$ sampling points). That is, evaluate the integrals (a) $\begin{array} { r } { I = \int _ { - 1 } ^ { 1 } [ x ^ { 2 } + \cos ( x / 2 ) ] d x } \end{array}$ and (b) $\begin{array} { r } { I = \int _ { - 1 } ^ { 1 } \left( 3 ^ { x } - \overset { \cdot } { x } \right) } \end{array}$ dx using the Newton-Cotes method.
|
||
|
||
Using Table 10–2 with three sampling points means we evaluate the function inside the integrand at $x = - 1 , x = 0$ , and $x = 1$ , and multiply each evaluated function by the respective Newton-Cotes numbers, $1 / 6 , 4 / 6 .$ , and $1 / 6 .$ . We then add these three products together and finally multiply this sum by the interval of integration $( h = 2 )$ as follows:
|
||
|
||
$$
|
||
I = 2 \left[ \frac {1}{6} y _ {0} + \frac {4}{6} y _ {1} + \frac {1}{6} y _ {2} \right] \tag {10.4.12}
|
||
$$
|
||
|
||
<!-- source-page: 487 -->
|
||
|
||
(a): Using Eq. (10.4.12), we obtain
|
||
|
||
$$
|
||
y _ {0} = x ^ {2} + \cos (x / 2) \text { evaluated at } x = - 1, \text { etc. as follows: }
|
||
$$
|
||
|
||
$$
|
||
y _ {0} = (- 1) ^ {2} + \cos (- 1 / 2 \mathrm{rad}) = 1. 8 7 7 5 8 2 6
|
||
$$
|
||
|
||
$$
|
||
y _ {1} = (0) ^ {2} + \cos (0 / 2) = 1 \tag {10.4.13}
|
||
$$
|
||
|
||
$$
|
||
y _ {2} = (1) ^ {2} + \cos (1 / 2 \mathrm{rad}) = 1. 8 7 7 5 8 2 6
|
||
$$
|
||
|
||
Substituting $y _ { 0 } - y _ { 2 }$ from Eq. (10.4.13) into Eq. (10.4.12), we obtain the evaluation of the integral as
|
||
|
||
$$
|
||
I = 2 \left[ \frac {1}{6} (1. 8 7 7 5 8 2 6) + \frac {4}{6} (1) + \frac {1}{6} (1. 8 7 7 5 8 2 6) \right] = 2. 5 8 5
|
||
$$
|
||
|
||
This solution compares exactly to the evaluation performed using Gaussian quadrature and to the exact solution. However, for higher-order functions the Gaussian quadrature method yields more accurate results than the Newton-Cotes method as illustrated by part (b) as follows:
|
||
|
||
(b): Using Eq. (10.4.12), we obtain
|
||
|
||
$$
|
||
y _ {0} = 3 ^ {(- 1)} - (- 1) = \frac {4}{3}
|
||
$$
|
||
|
||
$$
|
||
y _ {1} = 3 ^ {0} - 0 = 1
|
||
$$
|
||
|
||
$$
|
||
y _ {2} = 3 ^ {1} - (1) = 2
|
||
$$
|
||
|
||
Substituting $y _ { 0 } - y _ { 2 }$ into Eq. (10.4.12) we obtain I as
|
||
|
||
$$
|
||
I = 2 \left[ \frac {1}{6} \left(\frac {4}{3}\right) + \frac {4}{6} (1) + \frac {1}{6} (2) \right] = 2. 4 4 4
|
||
$$
|
||
|
||
The error is 2:444 2:427 ¼ 0:017. This error is larger than that found using Gaussian quadrature (see Example 10.2 (b)).
|
||
|
||
# 10.5 Evaluation of the Stiffness Matrix and Stress Matrix by Gaussian Quadrature
|
||
|
||
Evaluation of the Stiffness Matrix
|
||
|
||
For the two-dimensional element, we have shown in previous chapters that
|
||
|
||
$$
|
||
\underline {{k}} = \iint_ {A} \underline {{B}} ^ {T} (x, y) \underline {{D}} \underline {{B}} (x, y) h d x d y \tag {10.5.1}
|
||
$$
|
||
|
||
where, in general, the integrand is a function of x and y and nodal coordinate values.
|
||
|
||
<!-- source-page: 488 -->
|
||
|
||
|
||
|
||
<details>
|
||
<summary>flowchart</summary>
|
||
|
||
```mermaid
|
||
graph TD
|
||
A["Read in four Gauss points and weight functions s_i, t_i = ±0.5773 ... ; W_i, W_{j=i} = 1., 1."] --> B["Zero k̃^(e)"]
|
||
B --> C["DO i = 1, 4"]
|
||
C --> D["Let s = s_i, t = t_i"]
|
||
D --> E["Compute |J(s, t)|, B(s, t), D"]
|
||
E --> F["Compute k̃ = B^T DB |J|h"]
|
||
F --> G["k̃^(e) = k̃^(e) + kW_i W_j"]
|
||
G --> C
|
||
```
|
||
</details>
|
||
|
||
Figure 10–12 Flowchart to evaluate $\underline { { k } } ^ { ( e ) }$ by four-point Gaussian quadrature
|
||
|
||
We have shown in Section 10.3 that k for a quadrilateral element can be evaluated in terms of a local set of coordinates s-t, with limits from minus one to one within the element, and in terms of global nodal coordinates as given by Eq. (10.3.27). We repeat Eq. (10.3.27) here for convenience as
|
||
|
||
$$
|
||
\underline {{k}} = \int_ {- 1} ^ {1} \int_ {- 1} ^ {1} \underline {{B}} ^ {T} (s, t) \underline {{D}} \underline {{B}} (s, t) | \underline {{J}} | h d s d t \tag {10.5.2}
|
||
$$
|
||
|
||
where $| \underline { J } |$ is defined by Eq. (10.3.22) and $\underline { { B } }$ is defined by Eq. (10.3.18). In Eq. (10.5.2), each coe‰cient of the integrand $\underline { { B } } ^ { T } \underline { { D } } \underline { { B } } | \underline { { J } } |$ must be evaluated by numerical integration in the same manner as $f ( s , t )$ was integrated in Eq. (10.4.9).
|
||
|
||
A flowchart to evaluate k of Eq. (10.5.2) for an element using four-point Gaussian quadrature is given in Figure 10–12. The four-point Gaussian quadrature rule is relatively easy to use. Also, it has been shown to yield good results [7]. In Figure 10–12, in explicit form for four-point Gaussian quadrature (now using the single summation notation with $i = 1 , 2 , 3 , 4 )$ , we have
|
||
|
||
$$
|
||
\begin{array}{l} \underline {{k}} = \underline {{B}} ^ {T} (s _ {1}, t _ {1}) \underline {{D}} \underline {{B}} (s _ {1}, t _ {1}) | \underline {{J}} (s _ {1}, t _ {1}) | h W _ {1} W _ {1} \\ + \underline {{{B}}} ^ {T} (s _ {2}, t _ {2}) \underline {{{D}}} \underline {{{B}}} (s _ {2}, t _ {2}) | \underline {{{J}}} (s _ {2}, t _ {2}) | h W _ {2} W _ {2} \\ + \underline {{{B}}} ^ {T} (s _ {3}, t _ {3}) \underline {{{D}}} \underline {{{B}}} (s _ {3}, t _ {3}) | \underline {{{J}}} (s _ {3}, t _ {3}) | h W _ {3} W _ {3} \\ + \underline {{B}} ^ {T} (s _ {4}, t _ {4}) \underline {{D}} \underline {{B}} (s _ {4}, t _ {4}) | \underline {{J}} (s _ {4}, t _ {4}) | h W _ {4} W _ {4} \tag {10.5.3} \\ \end{array}
|
||
$$
|
||
|
||
where $s _ { 1 } = t _ { 1 } = - 0 . 5 7 7 3 , s _ { 2 } = - 0 . 5 7 7 3 , t _ { 2 } = 0 . 5 7 7 3 , s _ { 3 } = 0 . 5 7 7 3 , t _ { 3 } = - 0 . 5 7 7 3$ , and $s _ { 4 } = t _ { 4 } = 0 . 5 7 7 3$ as shown in Figure 10–9, and $W _ { 1 } = W _ { 2 } = W _ { 3 } = W _ { 4 } = 1 . 0 0 0$ .
|
||
|
||
<!-- source-page: 489 -->
|
||
|
||
Evaluate the stiffness matrix for the quadrilateral element shown in Figure 10–13 using the four-point Gaussian quadrature rule. Let $E = 30 \times 10^{6}$ psi and v = 0.25. The global coordinates are shown in inches. Assume h = 1 in.
|
||
|
||
Using Eq. (10.5.3), we evaluate the $\underline{k}$ matrix. Using the four-point rule, the four points are
|
||
|
||
$$
|
||
(s _ {1}, t _ {1}) = (- 0. 5 7 7 3, - 0. 5 7 7 3)
|
||
$$
|
||
|
||
$$
|
||
\begin{array}{l} \begin{array}{l} \left(s _ {2}, t _ {2}\right) = (- 0. 5 7 7 3, 0. 5 7 7 3) \\ \left(s _ {2}, t _ {2}\right) = (0. 5 7 7 3, 0. 5 7 7 3) \end{array} \tag {10.5.4a} \\ (s _ {3}, t _ {3}) = (0. 5 7 7 3, - 0. 5 7 7 3) \\ \end{array}
|
||
$$
|
||
|
||
$$
|
||
(s _ {4}, t _ {4}) = (0. 5 7 7 3, 0. 5 7 7 3)
|
||
$$
|
||
|
||
with weights $W_{1} = W_{2} = W_{3} = W_{4} = 1.000$ .
|
||
|
||
Therefore, by Eq. (10.5.3), we have
|
||
|
||
$$
|
||
\begin{array}{l} \underline {{k}} = \underline {{B}} ^ {T} (- 0. 5 7 7 3, - 0. 5 7 7 3) \underline {{D}} \underline {{B}} (- 0. 5 7 7 3, - 0. 5 7 7 3) \\ \times | \underline {{{J}}} (- 0. 5 7 7 3, - 0. 5 7 7 3) | (1) (1. 0 0 0) (1. 0 0 0) \\ + \underline {{{B}}} ^ {T} (- 0. 5 7 7 3, 0. 5 7 7 3) \underline {{{D}}} \underline {{{B}}} (- 0. 5 7 7 3, 0. 5 7 7 3) \\ \times | \underline {{{J}}} (- 0. 5 7 7 3, 0. 5 7 7 3) | (1) (1. 0 0 0) (1. 0 0 0) \\ + \underline {{{B}}} ^ {T} (0. 5 7 7 3, - 0. 5 7 7 3) \underline {{{D}}} \underline {{{B}}} (0. 5 7 7 3, - 0. 5 7 7 3) \\ \times | \underline {{{J}}} (0. 5 7 7 3, - 0. 5 7 7 3) | (1) (1. 0 0 0) (1. 0 0 0) \\ + \underline {{{B}}} ^ {T} (0. 5 7 7 3, 0. 5 7 7 3) \underline {{{D}}} \underline {{{B}}} (0. 5 7 7 3, 0. 5 7 7 3) \\ \times | \underline {{{J}}} (0. 5 7 7 3, 0. 5 7 7 3) | (1) (1. 0 0 0) (1. 0 0 0) \tag {10.5.4b} \\ \end{array}
|
||
$$
|
||
|
||
To evaluate $\underline{k}$ , we first evaluate $|\underline{J}|$ at each Gauss point by using Eq. (10.3.22). For instance, one part of $|\underline{J}|$ is given by
|
||
|
||
$$
|
||
\begin{array}{l} | \underline {{{J}}} (- 0. 5 7 7 3, - 0. 5 7 7 3) | = \frac {1}{8} [ 3 \quad 5 \quad 5 \quad 3 ] \\ \times \left[ \begin{array}{c c c c} 0 & 1 - (- 0. 5 7 7 3) & - 0. 5 7 7 3 - (- 0. 5 7 7 3) & - 0. 5 7 7 3 - 1 \\ - 0. 5 7 7 3 - 1 & 0 & - 0. 5 7 7 3 + 1 & - 0. 5 7 7 3 - (- 0. 5 7 7 3) \\ - 0. 5 7 7 3 - (- 0. 5 7 7 3) & - 0. 5 7 7 3 - 1 & 0 & - 0. 5 7 7 3 + 1 \\ 1 - (- 0. 5 7 7 3) & - 0. 5 7 7 3 + (- 0. 5 7 7 3) & - 0. 5 7 7 3 - 1 & 0 \end{array} \right] \\ \times \left\{ \begin{array}{l} 2 \\ 2 \\ 4 \\ 4 \end{array} \right\} = 1. 0 0 0 \tag {10.5.4c} \\ \end{array}
|
||
$$
|
||
|
||
|
||
|
||
<details>
|
||
<summary>text_image</summary>
|
||
|
||
(3, 4) (5, 4)
|
||
4 3
|
||
1 2
|
||
(3, 2) (5, 2)
|
||
x
|
||
</details>
|
||
|
||
Figure 10–13 Quadrilateral element for stiffness evaluation
|
||
|
||
<!-- source-page: 490 -->
|
||
|
||
Similarly,
|
||
|
||
$$
|
||
| \underline {{J}} (- 0. 5 7 7 3, 0. 5 7 7 3) | = 1. 0 0 0
|
||
$$
|
||
|
||
$$
|
||
| \underline {{J}} (0. 5 7 7 3, - 0. 5 7 7 3) | = 1. 0 0 0 \tag {10.5.4d}
|
||
$$
|
||
|
||
$$
|
||
| \underline {{{J}}} (0. 5 7 7 3, 0. 5 7 7 3) | = 1. 0 0 0
|
||
$$
|
||
|
||
Even though $|\underline{J}| = 1$ in this example, in general, $|\underline{J}| \neq 1$ and varies in space.
|
||
|
||
Then, using Eqs. (10.3.18) and (10.3.19), we evaluate $\underline{B}$ . For instance, one part of $\underline{B}$ is
|
||
|
||
$$
|
||
\underline {{{B}}} (- 0. 5 7 7 3, - 0. 5 7 7 3) = \frac {1}{| \underline {{{J}}} (- 0 . 5 7 7 3 , - 0 . 5 7 7 3) |} \left[ \begin{array}{l l l l} \underline {{{B}}} _ {1} & \underline {{{B}}} _ {2} & \underline {{{B}}} _ {3} & \underline {{{B}}} _ {4} \end{array} \right]
|
||
$$
|
||
|
||
where, by Eq. (10.3.19),
|
||
|
||
$$
|
||
\underline {{B}} _ {1} = \left[ \begin{array}{c c} a N _ {1, s} - b N _ {1, t} & 0 \\ 0 & c N _ {1, t} - d N _ {1, s} \\ c N _ {1, t} - d N _ {1, s} & a N _ {1, s} - b N _ {1, t} \end{array} \right] \tag {10.5.4e}
|
||
$$
|
||
|
||
and by Eqs. (10.3.20) and (10.3.21), $a, b, c, d, N_{1,s}$ , and $N_{1,t}$ are evaluated. For instance,
|
||
|
||
$$
|
||
a = \frac {1}{4} \left[ y _ {1} (s - 1) + y _ {2} (- 1 - s) + y _ {3} (1 + s) + y _ {4} (1 - s) \right]
|
||
$$
|
||
|
||
$$
|
||
= \frac {1}{4} \left\{2 (- 0. 5 7 7 3 - 1) + 2 [ - 1 - (0. 5 7 7 3) ] \right\} + 4 [ 1 + (- 0. 5 7 7 3) ] + 4 [ 1 - (- 0. 5 7 7 3) ]
|
||
$$
|
||
|
||
$$
|
||
= 1. 0 0 \tag {10.5.4f}
|
||
$$
|
||
|
||
with similar computations used to obtain b, c, and d. Also,
|
||
|
||
$$
|
||
N _ {1, s} = \frac {1}{4} (t - 1) = \frac {1}{4} (- 0. 5 7 7 3 - 1) = - 0. 3 9 4 3 \tag {10.5.4g}
|
||
$$
|
||
|
||
$$
|
||
N _ {1, t} = \frac {1}{4} (s - 1) = \frac {1}{4} (- 0. 5 7 7 3 - 1) = - 0. 3 9 4 3
|
||
$$
|
||
|
||
Similarly, $\underline{B}_2$ , $\underline{B}_3$ , and $\underline{B}_4$ must be evaluated like $\underline{B}_1$ , at $(-0.5773, -0.5773)$ . We then repeat the calculations to evaluate $\underline{B}$ at the other Gauss points [Eq. (10.5.4a)].
|
||
|
||
Using a computer program written specifically to evaluate $\underline{B}$ at each Gauss point and then $\underline{k}$ , we obtain the final form of $\underline{B}(-0.5773, -0.5773)$ as
|
||
|
||
$$
|
||
\underline {{{B}}} (- 0. 5 7 7 3, - 0. 5 7 7 3) =
|
||
$$
|
||
|
||
$$
|
||
\left[ \begin{array}{c c c c c c c c} - 0. 1 0 5 7 & 0 & 0. 1 0 5 7 & 0 & 0 & - 0. 1 0 5 7 & 0 & - 0. 3 9 4 3 \\ - 0. 1 0 5 7 & - 0. 1 0 5 7 & - 0. 3 9 4 3 & 0. 1 0 5 7 & 0. 3 9 4 3 & 0 & - 0. 3 9 4 3 & 0 \\ 0 & 0. 3 9 4 3 & 0 & 0. 1 0 5 7 & 0. 3 9 4 3 & 0. 3 9 4 3 & 0. 1 0 5 7 & - 0. 3 9 4 3 \end{array} \right] \tag {10.5.4h}
|
||
$$
|
||
|
||
with similar expressions for $\underline{B}(-0.5773,0.5773)$ , and so on.
|