김경종
20 KiB
5.19. The problem in Example 5.13 is modeled by two five-node plane stress elements and one three-node bar element:
text_image
P
Same data as in Example 5.13
(a) Establish in detail all matrices used in the formulation of the governing equilibrium equations but do not perform any integrations.
(b) Assume that you have evaluated the unknown nodal point displacements. Present graphically in an oblique/aerial view all displacements and stresses in the elements.
5.20. The 20-node brick element shown is loaded by a concentrated load at the location indicated. Calculate the consistent nodal point loads.
text_image
12 12 8 P 5 1 7.5 15 2 6 18 8 17 15 4 13 3 7 10 10 19 14 20 16 9 11 15 12 All dimensions in centimeters x y z
5.21. The element in Exercise 5.20 is to be used in dynamic analysis. Construct a reasonable lumped mass matrix of the element; use \rho = 7.8 \times 10^{-3} kg/cm ^{3} .
5.22. The 12-node three-dimensional element shown is loaded with the pressure loading indicated. Calculate the nodal point consistent load vector for nodes 1, 2, 7, and 8.
text_image
6 cm 1 3 cm 7 6 2 5 8 4 3 11 12 9 10 p0 z y
Pressure loading varies linearly with y and is constant in x
5.23. Evaluate the Jacobian matrix J of the following element as a function of r, s and plot det J "over" the element (in our oblique/aerial view).
text_image
3/2 1/2 2 5 1 1 6 1 y x 3 7 4 1 1 8 1/2 3/2
text_image
det J 2 6 y x 3 7 4 5 1 8 1/2 3/2
5.24. Plot, for node 9, the displacement interpolation functions and their x-derivatives for the nine-node element and the assemblage of two six-node triangles (formed using the interpolation functions in Fig. 5.11).
text_image
4 9 y 4 x
text_image
4 9 4
5.25. Prove that the interpolation functions in (5.36), with \Delta h defined in (5.37), define the same displacement assumptions as the functions in Fig. 5.11 (note that the origins of the coordinates used in the two formulations are different).
5.26. Collapse a 20-node brick element into a spatially isotropic 10-node tetrahedron (use the collapsing of sides in Fig. E5.16). Determine the correction that must be applied to the interpolation function h_{16} of the brick element in order to obtain the displacement assumption h_{6} of the tetrahedron (given in Fig. 5.13).
5.27. Consider the six-node isoparametric plane strain finite element shown.
line
| Point | x | y |
|---|---|---|
| 1 | 2 | 3 |
| 2 | 1 | 2 |
| 3 | 1 | 1 |
| 4 | 2 | 1 |
| 5 | 1.5 | 2.5 |
| 6 | 1.5 | 1 |
(a) Construct the interpolation functions, h_i(r, s), i = 1, \ldots, 6 , of the element.
(b) Prove in detail that this finite element does (or does not) satisfy all convergence requirements when used in a compatible finite element assemblage.
5.28. Consider a general isoparametric four-node element used in an assemblage of elements as shown.
text_image
R₂ R₁ m
Either plane strain or plane stress condition
(a) Prove that the nodal point forces defined as
\mathbf {F} ^ {(m)} = \int_ {V ^ {(m)}} \mathbf {B} ^ {(m) ^ {T}} \boldsymbol {\tau} ^ {(m)} d V ^ {(m)}
are in equilibrium for element m , where \pmb{\tau}^{(m)} = \mathbf{CB}^{(m)}\mathbf{U} has been calculated.
(b) Show that the sum of the nodal point forces at each node is in equilibrium with the applied external loads R_{i} (including the reactions). (Hint: Refer to Section 4.2.1, Fig. 4.2.)
5.29. Consider Table 5.1 and the case of angular distortion. Prove that the terms listed for the 12- and 16-node elements are indeed correct.
5.30. Consider Table 5.1 and the case of curved edge distortion. Prove that the terms listed in the column for the 8-, 9-, 12-, and 16-node elements are correct.
5.31. Consider the 4/1 isoparametric u/p element shown. Construct all matrices for the evaluation of the stiffness matrix of the element but do not perform any integrations.
Plane strain condition
text_image
s 2 1 3 r 5 3 4
Bulk modulus \kappa Shear modulus G
5.4 FORMULATION OF STRUCTURAL ELEMENTS
The concepts of geometry and displacement interpolations that have been employed in the formulation of two- and three-dimensional continuum elements can also be employed in the evaluation of beam, plate, and shell structural element matrices. However, whereas in the formulation of the continuum elements the displacements u, v, w (whichever are applicable) are interpolated in terms of nodal point displacements of the same kind, in the formulation of structural elements, the displacements u, v, and w are interpolated in terms of midsurface displacements and rotations. We will show that this procedure corresponds in essence to a continuum isoparametric element formulation with displacement constraints. In addition, there is of course the major assumption that the stress normal to the midsurface is zero. The structural elements are for these reasons appropriately called degenerate isoparametric elements, but frequently we still refer to them simply as isoparametric elements.
Considering the formulation of structural elements, we have already discussed briefly in Section 4.2.3 how beam, plate, and shell elements can be formulated using the Bernoulli beam and Kirchhoff plate theory, in which shear deformations are neglected. Using the Kirchhoff theory it is difficult to satisfy interelement continuity on displacements and edge rotations because the plate (or shell) rotations are calculated from the transverse displacements. Furthermore, using an assemblage of flat elements to represent a shell structure, a relatively large number of elements may be required in order to represent the shell geometry to sufficient accuracy.
Our objective in this section is to discuss an alternative approach to formulating beam, plate, and shell elements. The basis of this method is a theory that includes the effects of shear deformations. In this theory the displacements and rotations of the midsurface normals are independent variables, and the interelement continuity conditions on these quantities can be satisfied directly, as in the analysis of continua. In addition, if the concepts of isoparametric interpolation are employed, the geometry of curved shell surfaces is interpolated and can be represented to a high degree of accuracy. In the following sections we discuss first the formulation of beam and axisymmetric shell elements, where we can demonstrate in detail the basic principles used, and we then present the formulation of general plate and shell elements.
text_image
dw/dx Neutral axis w Beam section x Deformation of cross section
text_image
Element 2 Element 1 x
Boundary conditions between beam elements
w \left| _ {x ^ {- 0}} = w \right| _ {x ^ {+ 0}}; \left. \frac {d w}{d x} \right| _ {x ^ {- 0}} = \left. \frac {d w}{d x} \right| _ {x ^ {+ 0}}
(a) Beam deformations excluding shear effect
text_image
dw/dx - γ = β dw/dx γ Neutral axis w Beam section x Deformation of cross section
text_image
Element 1 Element 2 x
Boundary conditions between beam elements
w \left| _ {x ^ {- 0}} = w \right| _ {x ^ {+ 0}}
(b) Beam deformations including shear effect
\beta \left| _ {x ^ {- 0}} = \beta \right| _ {x ^ {+ 0}}
Figure 5.18 Beam deformation assumptions
5.4.1 Beam and Axisymmetric Shell Elements
Let us discuss first some basic assumptions pertaining to the formulation of beam elements. The basic assumption in beam bending analysis excluding shear deformations is that a normal to the midsurface (neutral axis) of the beam remains straight during deformation and that its angular rotation is equal to the slope of the beam midsurface. This kinematic assumption, illustrated in Fig. 5.18(a), corresponds to the Bernoulli beam theory and leads to the well-known beam-bending governing differential equation in which the transverse displacement w is the only variable (see Example 3.20). Therefore, using beam elements formulated with this theory, displacement countinuity between elements requires that w and dw/dx be continuous.
Considering now beam bending analysis with the effect of shear deformations, we retain the assumption that a plane section originally normal to the neutral axis remains plane, but because of shear deformations this section does not remain normal to the neutral axis. As illustrated in Fig. 5.18(b), the total rotation of the plane originally normal to the neutral axis of the beam is given by the rotation of the tangent to the neutral axis and the shear deformation,
\beta = \frac {d w}{d x} - \gamma \tag {5.56}
where \gamma is a constant shearing strain across the section. This kinematic assumption corresponds to Timoshenko beam theory (see S. H. Crandall, N. C. Dahl, and T. J. Lardner [A]). Since the actual shearing stress and strain vary over the section, the shearing strain \gamma in (5.56) is an equivalent constant strain on a corresponding shear area A_{s} ,
\tau = \frac {V}{A _ {s}}; \quad \gamma = \frac {\tau}{G}; \quad k = \frac {A _ {s}}{A} \tag {5.57}
where V is the shear force at the section being considered. Different assumptions may be used to evaluate a reasonable factor k (see S. Timoshenko and J. N. Goodier [A] and K. Washizu [B]). One simple procedure is to evaluate the shear correction factor using the condition that when acting on A_{s} , the constant shear stress in (5.57) must yield the same shear strain energy as the actual shearing stress (evaluated from beam theory) acting on the actual cross-sectional area A of the beam. Consider the following example.
EXAMPLE 5.23: Evaluate the shear correction factor k for a beam of rectangular cross section, width b, and depth h.
The shear strain energy \mathcal{U} of the beam per unit length is
\mathcal {U} = \int_ {A} \frac {1}{2 G} \tau_ {a} ^ {2} d A \tag {a}
where \tau_{a} is the actual shear stress, G is the shear modulus, and A is the cross-sectional area, A = bh.
In our finite element model, by assumption, the shear strain is constant over the cross-sectional area of the beam [see (5.56)]. Since in reality the shear strain varies over the beam cross section, we want to find an equivalent beam cross-sectional area A_{s} for our finite element model. This equivalency will be based on equating shear strain energies.
Hence using \mathcal{U} , given in (a), with the actual shear stress distribution, we can calculate A , from
\int_ {A} \frac {1}{2 G} \tau_ {a} ^ {2} d A = \int_ {A _ {s}} \frac {1}{2 G} \left(\frac {V}{A _ {s}}\right) ^ {2} d A _ {s} \tag {b}
where V is the total shearing force at the section,
V = \int_ {A} \tau_ {a} d A \tag {c}
If we use k = A_{s} / A , we obtain from (b),
k = \frac {V ^ {2}}{A \int_ {A} \tau_ {a} ^ {2} d A}
We now use (b) and (c) for the rectangular cross-section beam. Elementary beam theory gives (see S. H. Crandall, N. C. Dahl, and T. J. Lardner [A])
\tau_ {a} = \frac {3}{2} \frac {V}{A} \left[ \frac {(h / 2) ^ {2} - y ^ {2}}{(h / 2) ^ {2}} \right]
which gives k = \frac{5}{6} .
The finite element formulation of a beam element with the assumption in (5.56) is obtained using the basic virtual work expressions in (4.19) to (4.22). In the following we consider first, for illustrative purposes, the specific formulation of the beam element matrices corresponding to the simple element in Fig. 5.19, and we discuss afterward the formulation of more general three-dimensional beam elements, and the formulation of axisymmetric shell elements.
Two-Dimensional Straight Beam Elements
Figure 5.19 shows the two-dimensional rectangular cross-section beam considered. Using the general expression of the principle of virtual work with the assumptions discussed above we have (see Exercise 5.32)
E I \int_ {0} ^ {L} \left(\frac {d \beta}{d x}\right) \left(\frac {d \overline {{{\beta}}}}{d x}\right) d x + G A k \int_ {0} ^ {L} \left(\frac {d w}{d x} - \beta\right) \left(\frac {d \overline {{{w}}}}{d x} - \overline {{{\beta}}}\right) d x = \int_ {0} ^ {L} p \overline {{{w}}} d x + \int_ {0} ^ {L} m \overline {{{\beta}}} d x \tag {5.58}
where p and m are the transverse and moment loadings per unit length. Using now the interpolations
w = \sum_ {i = 1} ^ {q} h _ {i} w _ {i}; \quad \beta = \sum_ {i = 1} ^ {q} h _ {i} \theta_ {i} \tag {5.59}
where q is equal to the number of nodes used and the h_{i} are the one-dimensional interpolation functions listed in Fig. 5.3, we can directly employ the concepts of the isoparametric formulations discussed in Section 5.3 to establish all relevant element matrices. Let
w = \mathbf {H} _ {w} \hat {\mathbf {u}}; \quad \beta = \mathbf {H} _ {\beta} \hat {\mathbf {u}} \tag {5.60}
\frac {\partial w}{\partial x} = \mathbf {B} _ {w} \hat {\mathbf {u}}; \quad \frac {\partial \beta}{\partial x} = \mathbf {B} _ {\beta} \hat {\mathbf {u}}
E = Young's modulus; G = shear modulus;
k = \frac{5}{6};A = hb;I = \frac{bh^3}{12}
(b) 2-, 3-, and 4-node models; \theta_{i} = \beta_{i}, i = 1, \ldots, q (Interpolation functions are given in Fig. 5.3)
Figure 5.19 Formulation of two-dimensional beam element
where \hat{\mathbf{u}}^T = [w_1\ldots w_q\theta_1\ldots \theta_q]
\mathbf {H} _ {w} = [ h _ {1} \dots h _ {q} \quad 0 \dots 0 ] \tag {5.61}
\mathbf {H} _ {\beta} = [ 0 \dots 0 \quad h _ {1} \dots h _ {q} ]
and \mathbf{B}_w = J^{-1}\left[\frac{\partial h_1}{\partial r}\quad \dots \quad \frac{\partial h_q}{\partial r}\quad 0\quad \dots \quad 0\right] (5.62)
\mathbf {B} _ {\beta} = J ^ {- 1} \left[ \begin{array}{l l l l l l l} 0 & \dots & 0 & \frac {\partial h _ {1}}{\partial r} & \dots & \frac {\partial h _ {q}}{\partial r} \end{array} \right]
where J = \partial x / \partial r ; then we have for a single element,
\mathbf {K} = E I \int_ {- 1} ^ {1} \mathbf {B} _ {\beta} ^ {T} \mathbf {B} _ {\beta} \det J d r + G A k \int_ {- 1} ^ {1} \left(\mathbf {B} _ {w} - \mathbf {H} _ {\beta}\right) ^ {T} \left(\mathbf {B} _ {w} - \mathbf {H} _ {\beta}\right) \det J d r \tag {5.63}
\mathbf {R} = \int_ {- 1} ^ {1} \mathbf {H} _ {w} ^ {T} p \det J d r + \int_ {1} ^ {- 1} \mathbf {H} _ {\beta} ^ {T} m \det J d r
Also, in dynamic analysis the mass matrix can be calculated using the d'Alembert principle [see (4.23)]; hence,
\mathbf {M} = \int_ {- 1} ^ {1} \left[ \begin{array}{l} \mathbf {H} _ {w} \\ \mathbf {H} _ {\beta} \end{array} \right] ^ {T} \left[ \begin{array}{c c} \rho b h & 0 \\ 0 & \frac {\rho b h ^ {3}}{1 2} \end{array} \right] \left[ \begin{array}{l} \mathbf {H} _ {w} \\ \mathbf {H} _ {\beta} \end{array} \right] \det J d r \tag {5.64}
In these evaluations we are using the natural coordinate system of the beam because this is effective in the formulation of more general beam, plate, and shell elements. However, when considering a straight beam of constant cross section, the integrals can also be evaluated efficiently without using the natural coordinate system, as demonstrated in the following example.
EXAMPLE 5.24: Evaluate the detailed expressions for calculation of the stiffness matrix and the load vector of the three-node beam element shown in Fig. E5.24.
Figure E5.24 Three-node beam element
The interpolation functions to be used are listed in Fig. 5.3. These functions are given in terms of r and yield
x = \sum_ {i = 1} ^ {3} h _ {i} x _ {i}
Using x_{1} = 0, x_{2} = L, x_{3} = L / 2 , we obtain
x = \frac {L}{2} (1 + r)
Hence, the interpolation functions in terms of x are
h _ {1} = \frac {2 x ^ {2}}{L ^ {2}} - \frac {3 x}{L} + 1
h _ {2} = \frac {2 x ^ {2}}{L ^ {2}} - \frac {x}{L}
h _ {3} = \frac {4 x}{L} - \frac {4 x ^ {2}}{L ^ {2}}
Using the notation ( )' ≡ ∂/∂x it follows that
h _ {1} ^ {\prime} = \frac {4 x}{L ^ {2}} - \frac {3}{L}
h _ {2} ^ {\prime} = \frac {4 x}{L ^ {2}} - \frac {1}{L}
h _ {3} ^ {\prime} = \frac {4}{L} - \frac {8 x}{L ^ {2}}
Hence, with the degrees of freedom ordered as in (5.61), we have
\mathbf {K} = \frac {E h ^ {3}}{1 2} \int_ {0} ^ {L} \left[ \begin{array}{l} 0 \\ 0 \\ 0 \\ h _ {1} ^ {\prime} \\ h _ {2} ^ {\prime} \\ h _ {3} ^ {\prime} \end{array} \right] \left[ \begin{array}{l l l l l l} 0 & 0 & 0 & h _ {1} ^ {\prime} & h _ {2} ^ {\prime} & h _ {3} ^ {\prime} \end{array} \right] d x
+ \frac {5 G h}{6} \int_ {0} ^ {L} \left[ \begin{array}{c} h _ {1} ^ {\prime} \\ h _ {2} ^ {\prime} \\ h _ {3} ^ {\prime} \\ - h _ {1} \\ - h _ {2} \\ - h _ {3} \end{array} \right] \left[ \begin{array}{c c c c c c} h _ {1} ^ {\prime} & h _ {2} ^ {\prime} & h _ {3} ^ {\prime} & - h _ {1} & - h _ {2} & - h _ {3} \end{array} \right] d x
and \mathbf{R}^T = -P[-\frac{1}{9} \frac{2}{9} \frac{8}{9} 0 0 0]
The element in Fig. 5.19 is a pure displacement-based element (assuming exact integration of all integrals) and can be employed provided three or four nodes are used (and the interior nodes are located at the midpoint and third-points, respectively). However, if the two-node element is employed or the interior nodes of the three- and four-node elements are not located at the midpoint and third-points, respectively, the use of the element cannot be recommended because the shearing deformations are not represented to sufficient accuracy. This deficiency is particularly pronounced when the element is thin.
In order to obtain some insight into the behavior of these elements when the beam becomes thin, we recall that the principle of virtual work is equivalent to the stationarity of the total potential (see Example 4.4). For the beam formulation the total potential is given by
\Pi = \frac {E I}{2} \int_ {0} ^ {L} \left(\frac {d \beta}{d x}\right) ^ {2} d x + \frac {G A k}{2} \int_ {0} ^ {L} \left(\frac {d w}{d x} - \beta\right) ^ {2} d x - \int_ {0} ^ {L} p w d x - \int_ {0} ^ {L} m \beta d x \tag {5.65}
where the first two integrals represent, respectively, the bending and shearing strain energies and the last two integrals represent the potential of the loads.
Let us consider the total potential \Pi
\tilde {\Pi} = \int_ {0} ^ {L} \left(\frac {d \beta}{d x}\right) ^ {2} d x + \frac {G A k}{E I} \int_ {0} ^ {L} \left(\frac {d w}{d x} - \beta\right) ^ {2} d x \tag {5.66}