MultiPhysicsVault/.raw/AFirstCourseInTheFiniteElementMethod/AFirstCourseInTheFiniteElementMethod_053.md

as follows:

\frac {\partial f}{\partial x} = \frac {\left| \begin{array}{c c c} \frac {\partial f}{\partial s} & \frac {\partial y}{\partial s} & \frac {\partial z}{\partial s} \\ \frac {\partial f}{\partial t} & \frac {\partial y}{\partial t} & \frac {\partial z}{\partial t} \\ \frac {\partial f}{\partial z ^ {\prime}} & \frac {\partial y}{\partial z ^ {\prime}} & \frac {\partial z}{\partial z ^ {\prime}} \end{array} \right|}{| \underline {{{J}}} |} \quad \frac {\partial f}{\partial y} = \frac {\left| \begin{array}{c c c} \frac {\partial x}{\partial s} & \frac {\partial f}{\partial s} & \frac {\partial z}{\partial s} \\ \frac {\partial x}{\partial t} & \frac {\partial f}{\partial t} & \frac {\partial z}{\partial t} \\ \frac {\partial x}{\partial z ^ {\prime}} & \frac {\partial f}{\partial z ^ {\prime}} & \frac {\partial z}{\partial z ^ {\prime}} \end{array} \right|}{| \underline {{{J}}} |} \tag {11.3.8}

\frac {\partial f}{\partial z} = \frac {\left| \begin{array}{c c c} \frac {\partial x}{\partial s} & \frac {\partial y}{\partial s} & \frac {\partial f}{\partial s} \\ \frac {\partial x}{\partial t} & \frac {\partial y}{\partial t} & \frac {\partial f}{\partial t} \\ \frac {\partial x}{\partial z ^ {\prime}} & \frac {\partial y}{\partial z ^ {\prime}} & \frac {\partial f}{\partial z ^ {\prime}} \end{array} \right|}{| \underline {{J}} |}

Using Eqs. (11.3.8) by substituting u, v, and then w for f and using the definitions of the strains, we can express the strains in terms of natural coordinates (s, t, z') to obtain an equation similar to Eq. (10.3.14). In compact form, we can again express the strains in terms of the shape functions and global nodal coordinates similar to Eq. (10.3.15). The matrix B, given by a form similar to Eq. (10.3.17), is now a function of s, t, and z' and is of order 6 \times 24 .

The 24 \times 24 stiffness matrix is now given by

[ k ] = \int_ {- 1} ^ {1} \int_ {- 1} ^ {1} \int_ {- 1} ^ {1} [ B ] ^ {T} [ D ] [ B ] | \underline {{{J}}} | d s d t d z ^ {\prime} \tag {11.3.9}

Again, it is best to evaluate [k] by numerical integration (also see Section 10.4); that is, we evaluate (integrate) the eight-node hexahedral element stiffness matrix using a 2 \times 2 \times 2 rule (or two-point rule). Actually, eight points defined in Table 11–1 are used to evaluate k as

\underline {{k}} = \sum_ {i = 1} ^ {8} \underline {{B}} ^ {T} (s _ {i}, t _ {i}, z _ {i} ^ {\prime}) \underline {{D}} \underline {{B}} (s _ {i}, t _ {i}, z _ {i} ^ {\prime}) | \underline {{J}} (s _ {i}, t _ {i}, z _ {i} ^ {\prime}) | W _ {i} W _ {j} W _ {k}

where W_{i} = W_{j} = W_{k} for the two-point rule.

As is true with the bilinear quadrilateral element described in Section 10.3, the eight-noded linear hexahedral element cannot model beam-bending action well because the element sides remain straight during the element deformation. During the bending process, the elements will be stretched and can shear lock. This concept of shear locking is described in more detail in [12] along with ways to remedy it. However, the quadratic hexahedral element described subsequently remedies the shear locking problem.

Table 11–1 Table of Gauss points for linear hexahedral element with associated weightsa

Points, i	$s_i$	$t_i$	$z_i'$	Weight, $W_i$
1	$-1/\sqrt{3}$	$-1/\sqrt{3}$	$1/\sqrt{3}$	1
2	$1/\sqrt{3}$	$-1/\sqrt{3}$	$1/\sqrt{3}$	1
3	$1/\sqrt{3}$	$1/\sqrt{3}$	$1/\sqrt{3}$	1
4	$-1/\sqrt{3}$	$1/\sqrt{3}$	$1/\sqrt{3}$	1
5	$-1/\sqrt{3}$	$-1/\sqrt{3}$	$-1/\sqrt{3}$	1
6	$1/\sqrt{3}$	$-1/\sqrt{3}$	$-1/\sqrt{3}$	1
7	$1/\sqrt{3}$	$1/\sqrt{3}$	$-1/\sqrt{3}$	1
8	$-1/\sqrt{3}$	$1/\sqrt{3}$	$-1/\sqrt{3}$	1

^ a 1 / { \sqrt { 3 } } = 0 . 5 7 7 3 5 .

Quadratic Hexahedral Element

For the quadratic hexahedral element shown in Figure 11–6, we have a total of 20 nodes with the inclusion of a total of 12 midside nodes.

The function describing the element geometry for x in terms of the 20 \boldsymbol { a _ { i } } ^ { \prime } \mathbf { s } is

\begin{array}{l} x = a _ {1} + a _ {2} s + a _ {3} t + a _ {4} z ^ {\prime} + a _ {5} s t + a _ {6} t z ^ {\prime} + a _ {7} z ^ {\prime} s + a _ {8} s ^ {2} + a _ {9} t ^ {2} \\ + a _ {1 0} z ^ {\prime 2} + a _ {1 1} s ^ {2} t + a _ {1 2} s t ^ {2} + a _ {1 3} t ^ {2} z ^ {\prime} + a _ {1 4} t z ^ {\prime 2} + a _ {1 5} z ^ {\prime 2} s \\ + a _ {1 6} z ^ {\prime} s ^ {2} + a _ {1 7} s t z ^ {\prime} + a _ {1 8} s ^ {2} t z ^ {\prime} + a _ {1 9} s t ^ {2} z ^ {\prime} + a _ {2 0} s t z ^ {\prime 2} \tag {11.3.10} \\ \end{array}

Similar expressions describe the y and z coordinates.

The development of the stiffness matrix follows the same steps we outlined before for the linear hexahedral element, where the shape functions now take on new forms. Again, letting s _ { i } , t _ { i } , z _ { i } ^ { \prime } = \pm 1 , we have for the corner nodes ( i = 1 , 2 , \dots , 8 ) ,

N _ {i} = \frac {(1 + s s _ {i}) (1 + t t _ {i}) (1 + z ^ {\prime} z _ {i} ^ {\prime})}{8} (s s _ {i} + t t _ {i} + z ^ {\prime} z _ {i} ^ {\prime} - 2) \tag {11.3.11}

For the midside nodes at s _ { i } = 0 , t _ { i } = \pm 1 , z _ { i } ^ { \prime } = \pm 1 \ ( i = 1 7 , 1 8 , 1 9 , 2 0 ) , we have

N _ {i} = \frac {(1 - s ^ {2}) (1 + t t _ {i}) (1 + z ^ {\prime} z _ {i} ^ {\prime})}{4} \tag {11.3.12}

text_image

3 11 10 4 12 2 9 1 1 17 15 8 16 13 5 18 13 6 7 t s z'

Figure 11–6 Quadratic hexahedral isoparametric element

natural_image

3D wireframe model of a mechanical part with no visible text or symbols

Figure 11–7 Finite element model of a forging using linear and quadratic solid elements

For the midside nodes at s _ { i } = \pm 1 , t _ { i } = 0 , z _ { i } ^ { \prime } = \pm 1 \ ( i = 1 0 , 1 2 , 1 4 , 1 6 ) , we have

N _ {i} = \frac {(1 + s s _ {i}) (1 - t ^ {2}) (1 + z ^ {\prime} z _ {i} ^ {\prime})}{4} \tag {11.3.13}

Finally, for the midside nodes at s _ { i } = \pm 1 , t _ { i } = \pm 1 , z _ { i } ^ { \prime } = 0 \ ( i = 9 , 1 1 , 1 3 , 1 5 ) , we have

N _ {i} = \frac {(1 + s s _ {i}) (1 + t t _ {i}) (1 - z ^ {\prime 2})}{4} \tag {11.3.14}

Note that the stiffness matrix for the quadratic solid element is of order 6 0 \times 6 0 because it has 20 nodes and 3 degrees of freedom per node.

The stiffness matrix for this 20-node quadratic solid element can be evaluated using a 3 \times 3 \times 3 rule (27 points). However, a special 14-point rule may be a better choice [9, 10].

As with the eight-noded plane element of Section 10.6 (Figure 10–16), the 20-node solid element is also called a serendipity element.

Figures 1–7 and 11–7 show applications of the use of linear and quadratic (curved sides) solid elements to model three-dimensional solids.

Finally, commercial computer programs, such as [11] (also see references [46–56] of Chapter 1), are available to solve three-dimensional problems. Figures 11–8, 11–9, and 11–10 show a steel foot pedal, a hollow cast-iron member, and an alternator bracket solved using a computer program [11]. We emphasize that these problems have been solved using the three-dimensional element as opposed to using a twodimensional element, such as described in Chapters 6 and 8, as these problems have a three-dimensional stress state occurring in them. That is, the three normal and three shear stresses are of similar order of magnitude in some parts of the foot pedal, the cast-iron member, and the alternator bracket. The most accurate results will then occur when modeling these problems using the three-dimensional brick or tetrahedral elements (or a combination of both).

For the foot pedal, the largest principal stress was 4111 psi and the largest von Mises stress was 4023 psi, both located at the interior corner of the elbow.

text_image

0.5 in. 2 in. D' D A' A x 10 in. B' B C C' E F F 1 in. 0.25 in. 5 lb G' G H' H z y z

Figure 11–8 Three-dimensional steel foot pedal

text_image

18,000 N 18,000 N 0.275 m 0.025 m 0.05 m 0.025 m 0.275 m 0.025 m 0.4 m 0.2 m 0.6 m x y z

Figure 11–9 Cast iron hollow member, E = 1 6 5 \mathsf { G P a } with opening on frontside

natural_image

3D rendered mechanical component with mesh texture, no visible text or symbols

Figure 11–10 Meshed model of an alternator bracket (Courtesy of Andrew Heckman, Design Engineer, Seagrave Fire Apparatus, LLC)

The maximum displacement was 0.01234 in. down at the free end point H in Figure 11–8. The model had 889 nodes and 648 brick elements.

For the cast-iron member, the maximum principal stress was 19 MPa and the maximum von Mises stress was 23.2 MPa. Both results occurred at the top near the fixed wall support. The free end vertical displacement was 0.300 mm.

Table 11–2 Table comparing results for cantilever beam modeled using 4-noded-tetrahedral, 10-noded tetrahedral, 8-noded brick, and 20-noded brick element

Solid Element Used	Number of Nodes	Number of Degrees of Freedom	Number of Elements	Free End Displ., in.	Principal Stress, psi
4-noded tet	30	90	61	0.0053	562
4-noded tet	415	1245	1549	0.0282	2357
4-noded tet	896	2688	3729	0.0420	3284
4-noded tet	1658	4974	7268	0.0548	4056
10-noded tet	144	432	61	0.1172	6601
10-noded tet	2584	7752	1549	0.1277	7970
8-noded brick	64	192	27	0.1190	5893
8-noded brick	343	1029	216	0.1253	6507
8-noded brick	1331	3993	1000	0.1277	6836
20-noded brick	208	624	27	0.1250	7899
20-noded brick	1225	3675	216	0.1285	8350
20-noded brick	4961	14,883	1000	0.1297	8323
Classical solution				0.1286	6940

(Mr. William Gobeli for creating the results for Table 11–2)

For the alternator bracket made of ASTM-A36 hot-rolled steel, the model consisted of 13,298 solid brick elements and 10,425 nodes. A total load of 1000 lb was applied downward to the flat front face piece. The bracket back side was constrained against displacement. The largest von Mises stress was 11,538 psi located at the top surface near the center (narrowest) section of the bracket. The largest vertical deflection was 0.01623 in. at the front tip of the outer edge of the alternator bracket.

It has been shown [3] that use of the simple eight-noded hexahedral element yields better results than use of the constant-strain tetrahedral discussed in Section 11.1. Table 11.2 also illustrates the comparison between the corner-noded (constant-strain) tetrahedral, the linear-strain tetrahedral (mid-edge nodes added), the 8-noded brick, and the 20-noded brick models for a three-dimensional cantilever beam of length 100 in., base 6 in., and height 12 in. The beam has an end load of 10,000 lb acting upward and is made of steel ( E = 3 0 \times 1 0 ^ { 6 } psi). A typical 8-noded brick model with the principal stress plot is shown in Figure 11–11. The classical beam theory solution for the vertical displacement and bending stress is also included for comparison. We can observe that the constant-strain tetrahedral gives very poor results, whereas the linear

bar

Bricks	Stress (bcf/(in^2))
1	7899.803
2	6776.506
3	5653.207
4	4529.908
5	3406.61
6	2283.312
7	1160.014
8	36.71533
9	-1086.583
10	-2209.881
11	-3333.179

Figure 11–11 Eight-noded brick model (27 Bricks) showing principal stress plot

tetrahedral gives much better results. This is because the linear-strain model predicts the beam-bending behavior much better. The 8-noded and 20-noded brick models yield similar but accurate results compared to the classical beam theory results.

In summary, the use of the three-dimensional elements results in a large number of equations to be solved simultaneously. For instance, a model using a simple cube with, say, 20 by 20 by 20 nodes (¼ 8000 total nodes) for a region requires 8000 times 3 degrees of freedom per node (¼ 24; 000) simultaneous equations.

References [4–7] report on early three-dimensional programs and analysis procedures using solid elements such as a family of subparametric curvilinear elements, linear tetrahedral elements, and 8-noded linear and 20-noded quadratic isoparametric elements.

d References

[1] Martin, H. C., ‘‘Plane Elasticity Problems and the Direct Stiffness Method.’’ The Trend in Engineering, Vol. 13, pp. 5–19, Jan. 1961.
[2] Gallagher, R. H., Padlog, J., and Bijlaard, P. P., ‘‘Stress Analysis of Heated Complex Shapes,’’ Journal of the American Rocket Society, pp. 700–707, May 1962.
[3] Melosh, R. J., ‘‘Structural Analysis of Solids,’’ Journal of the Structural Division, American Society of Civil Engineers, pp. 205–223, Aug. 1963.
[4] Chacour, S., ‘‘DANUTA, a Three-Dimensional Finite Element Program Used in the Analysis of Turbo-Machinery,’’ Transactions of the American Society of Mechanical Engineers, Journal of Basic Engineering, March 1972.
[5] Rashid, Y. R., ‘‘Three-Dimensional Analysis of Elastic Solids-I: Analysis Procedure,’’ International Journal of Solids and Structures, Vol. 5, pp. 1311–1331, 1969.
[6] Rashid, Y. R., ‘‘Three-Dimensional Analysis of Elastic Solids-II: The Computational Problem,’’ International Journal of Solids and Structures, Vol. 6, pp. 195–207, 1970.

[7] Three-Dimensional Continuum Computer Programs for Structural Analysis, Cruse, T. A., and Griffin, D. S., eds., American Society of Mechanical Engineers, 1972.
[8] Zienkiewicz, O. C., The Finite Element Method, 3rd ed., McGraw-Hill, London, 1977.
[9] Irons, B. M., ‘‘Quadrature Rules for Brick Based Finite Elements,’’ International Journal for Numerical Methods in Engineering, Vol. 3, No. 2, pp. 293–294, 1971.
[10] Hellen, T. K., ‘‘Effective Quadrature Rules for Quadratic Solid Isoparametric Finite Elements,’’ International Journal for Numerical Methods in Engineering, Vol. 4, No. 4, pp. 597–599, 1972.
[11] Linear Stress and Dynamics Reference Division, Docutech On-line Documentation, Algor, Inc., Pittsburgh, PA.
[12] Cook, R. D., Malkus, D. S., Plesha, M. E., and Witt, R. J., Concepts and Applications of Finite Element Analysis, 4th ed., Wiley, New York, 2002.

d Problems

11.1 Evaluate the matrix B for the tetrahedral solid element shown in Figure P11–1.

text_image

y 2 (0, 2, 0) 4 (2, 0, 0) 3 x 1 (0, 0, 2) z

text_image

y 4 (0, 2, 0) (1, 0, 0) (3, 0, 0) 1 3 x 2 (0, 0, 2) z

(b)
Figure P11–1

11.2 Evaluate the stiffness matrix for the elements shown in Figure P11–1. Let E = 3 0 \times 1 0 ^ { 6 } psi and \nu = 0 . 3 :
11.3 For the element shown in Figure P11–1, assume the nodal displacements have been determined to be

u _ {1} = 0. 0 0 5 \text {   in.   } \quad v _ {1} = 0. 0 \quad w _ {1} = 0. 0

u _ {2} = 0. 0 0 1 \text {   in.   } \quad v _ {2} = 0. 0 \quad w _ {2} = 0. 0 0 1 \text {   in.   }

u _ {3} = 0. 0 0 5 \text {   in.   } \quad v _ {3} = 0. 0 \quad w _ {3} = 0. 0

u _ {4} = - 0. 0 0 1 \text {   in.   } \quad v _ {4} = 0. 0 \quad w _ {4} = 0. 0 0 5 \text {   in.   }

Determine the strains and then the stresses in the element. Let E = 3 0 \times 1 0 ^ { 6 } psi and \nu = 0 . 3 .

11.4 What is special about the strains and stresses in the tetrahedral element?

11.5 Show that for constant body force Z _ { b } acting on an element ( X _ { b } = 0 and Y _ { b } = 0 ) ,

\{f _ {b i} \} = \frac {V}{4} \left\{ \begin{array}{l} 0 \\ 0 \\ Z _ {b} \end{array} \right\}

where \{ f _ { b i } \} represents the body forces at node i of the element with volume V.

11.6 Evaluate the B matrix for the tetrahedral solid element shown in Figure P11–6. The coordinates are in units of millimeters.

text_image

y (25, 25, 0) 2 (10, 0, 0) 1 4 (40, 0, 0) x z 3 (25, 0, 25)

text_image

(10, 7, 0) 2 (4, 2, 0) 1 4 (12, 2, 0) 3 (10, 2, 5) x y z

(b)
Figure P6–6

11.7 For the element shown in Figure P11–6, assume the nodal displacements have been determined to be

u _ {1} = 0. 0 \quad v _ {1} = 0. 0 \quad w _ {1} = 0. 0

u _ {2} = 0. 0 1 \mathrm{mm} \quad v _ {2} = 0. 0 2 \mathrm{mm} \quad w _ {2} = 0. 0 1 \mathrm{mm}

u _ {3} = 0. 0 2 \mathrm{mm} \quad v _ {3} = 0. 0 1 \mathrm{mm} \quad w _ {3} = 0. 0 0 5 \mathrm{mm}

u _ {4} = 0. 0 \quad v _ {4} = 0. 0 1 \mathrm{mm} \quad w _ {4} = 0. 0 1 \mathrm{mm}

Determine the strains and then the stresses in the element. Let E = 2 1 0 \ \mathrm { G P a } and \nu = 0 . 3 .

11.8 For the linear strain tetrahedral element shown in Figure P11–8, (a) express the displacement fields u, v, and w in the x, y and z directions, respectively. Hint: There are

natural_image

Geometric diagram of a polyhedron with solid and dashed edges (no text or labels)

Figure P11–8

10 nodes each with three translational degrees of freedom, u _ { i } , v _ { i } , , and w _ { i } . Also look at the linear strain triangle given by Eq. (8.1.2) or the expansion of Eqs. (11.2.2).

11.9 Figure P11–9 shows how solid and plane elements may be connected. What restriction must be placed on the externally applied loads for this connection to be acceptable?

text_image

z x Solid elements Plane elements y

Figure P11–9

11.10 Express the explicit shape functions N _ { 2 } through N _ { 8 } . , similar to N _ { 1 } given by Eq. (11.3.4), for the linear hexahedral element shown in Figure 11–5 on page 501.
11.11 Express the explicit shape functions for the corner nodes of the quadratic hexahedral element shown in Figure 11–6 on page 504.
11.12 Write a computer program to evaluate \underline { { k } } of Eq. (11.3.9) using a 2 \times 2 \times 2 Gaussian quadrature rule.

Solve the following problems using a computer program.

11.13 Determine the deflections at the four corners of the free end of the structural steel cantilever beam shown in Figure P11–13. Also determine the maximum principal stress.

text_image

z 800 lb 0.8 in. x 20 in. 2.5 in. y E = 30 × 10⁶ psi v = 0.3

Figure P11–13

11.14 A portion of a structural steel brake pedal in a vehicle is modeled as shown in Figure P11–14. Determine the maximum deflection at the pedal under a line load of 5 lb/in. as shown.

text_image

0.5 in. 2 in. 10 in. 5 lb/in. 1 in. 5 in. 0.25 in.

Figure P11–14

11.15 For the compressor flap valve shown in Figure P11–15, determine the maximum operating pressure such that the material yield stress is not exceeded with a factor of safety of two. The valve is made of hardened 1020 steel with a modulus of elasticity of 30 million psi and a yield strength of 62,000 psi. The valve thickness is a uniform 0.018 in. The value clip ears support the valve at opposite diameters. The pressure load is applied uniformly around the annular region.

text_image

Ø2.250 in. Pressure R0.040 in. Ø1.750 in. 2.750 in. 0.600 in. R0.250 in. R0.100 in. 0.200 in. 0.080 in.

△ denotes fixed boundary.

Figure P11–15

11.16 An S-shaped block used in force measurement as shown in Figure P11.16 is to be designed for a pressure of 1000 psi applied uniformly to the top surface. Determine the