MultiPhysicsVault/.raw/FiniteElementProcedures/FiniteElementProcedures_046.md

The transverse shear strain interpolation matrix is obtained from the shear interpolation given in Fig. 5.27 (and the tying procedure) written as

\left[ \begin{array}{c c c c c c c c c c} 1 & r & s & r s & s ^ {2} & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & r & s & r s & r ^ {2} \end{array} \right] \boldsymbol {\alpha} \tag {e}

where \alpha^T = [a_1\quad b_1\quad c_1\quad d_1\quad e_1\quad |\quad a_2\quad b_2\quad c_2\quad d_2\quad e_2]

The values in the vector \alpha are expressed in terms of the nodal point displacements and rotations using the tying relations. For example, since point A is at x = \frac{3}{2} [1 + 1 / \sqrt{3}], y = 2 , we have

\begin{array}{l} \gamma_ {x z} \mid_ {A} = a _ {1} + b _ {1} \left(\frac {3}{2}\right) \left(1 + \frac {1}{\sqrt {3}}\right) + c _ {1} (2) + d _ {2} (3) \left(1 + \frac {1}{\sqrt {3}}\right) + e _ {1} (4) \\ = \left. \left(\frac {\partial w}{\partial x} - \beta_ {x}\right) \right| _ {\text {at} r = 1 / \sqrt {3}, s = 1} \tag {f} \\ \end{array}

Of course, \partial w / \partial x is given by (a) and the section rotation \beta_{x} is given by (5.99) with the h_i corresponding to nine nodes.

Using all 10 tying relations in Fig. 5.27 as in (f), we can solve for the entries in (e) in terms of the nodal point displacements and rotations.

The numerical performance of the MITCn elements has been published by K. J. Bathe, M. L. Bucalem, and F. Brezzi [A]. However, let us briefly note that

The element matrices are all evaluated using full numerical Gauss integration (see Fig. 5.27).

The elements do not contain any spurious zero energy modes.

The elements pass the pure bending patch test (see Fig. 4.18).

To illustrate the performance of the elements and introduce a valuable test problem consider Figs. 5.28 to 5.32. In Fig. 5.28 the test problem is stated. We note that the transverse displacement and the section rotations are prescribed along the complete boundary of the square plate and that in this problem there are no boundary layers (as encountered in practical analyses; see B. Häggblad and K. J. Bathe [A]). Therefore, the numerically calculated orders of convergence should be close to the analytically predicted values. Figure 5.29 shows results obtained using uniform meshes, and these results compare well with the analytically predicted behavior (these predictions assume uniform meshes). Figures 5.30 and 5.31 show results obtained using a sequence of quasi-uniform ^{6} meshes, and we observe that the orders of convergence are not drastically affected by the element distortions. Finally, the convergence of the transverse shear strains, as predicted numerically, is shown in Fig. 5.32. In these specific finite element solutions the shear strains are predicted with surprisingly high orders of convergence (which in general of course cannot be expected).

text_image

(-1, 1) (1, 1) (-1, -1) (1, -1) y x

(a) Square plate considered in ad hoc plate bending problem; transverse loading p = 0, nonzero boundary conditions. A typical 4-node element is shown. The dashed line indicates the subdivision used for the triangular element meshes; h = 2/N, where N = number of elements per side.

(b) Exact transverse displacement and rotations: w = \sin kx e^{ky} + \sin k e^{-k} ; \theta_x = k \sin kx e^{ky} ; \theta_y = -k \cos kx e^{ky}

(c) Test problem: Prescribe the functional values of w , \theta_x , and \theta_y on the complete boundary and p = 0 , calculate interior values; k is a chosen constant; we use k = 5

Figure 5.28 Ad-hoc test problem for plate bending elements

line

Group	MITC7 (Log h)	MITC12 (Log h)	MITC4 (Log h)	MITC9 (Log h)	MITC16 (Log h)
MITC7	1.8	3.2	2.3	2.8	3.1
MITC12	0.5	2.9	2.5	2.7	2.9
MITC4	1.0	2.4	1.8	2.2	2.5
MITC9	0.0	1.4	0.8	1.2	1.6
MITC16	-1.0	0.0	-0.2	-0.4	-0.2

Figure 5.29 (a) Convergence of section rotations in analysis of ad-hoc problem using uniform meshes. The error measure is E = \|\beta - \beta_{h}\|_{1} . (b) Convergence of gradient of vertical displacement in analysis of ad-hoc problem using uniform meshes. The error measure is E = \|\nabla w - \nabla w_{h}\|_{0} .

natural_image

Pure geometric grid pattern with diagonal lines and a dashed diagonal line (no text or symbols)

▲MITC7
oMITC12

natural_image

Geometric pattern of intersecting diagonal lines forming a grid (no text or symbols)

♦MITC4
×MITC9
Figure 5.30 Two typical distorted meshes used in analysis of ad-hoc problem.----, indicates the subdivision used for the triangular element meshes.
▲MITC7
oMITC12
◇MITC4
×MITC9
+MITC16

+MITC16

line

Log h	Log E (Circle)	Log E (Triangle)
-1.0	1.2	2.3
-0.5	1.8	2.7
0.0	2.5	3.0

(a)

line

Log h	Series 1	Series 2	Series 3
-1.0	0.0	0.0	0.0
-0.5	0.0	0.0	0.0
0.0	0.0	0.0	0.0

line

Log h	Log E (Triangle)	Log E (Circle)
-1.0	1.5	0.5
-0.5	2.0	1.0
0.0	2.5	2.0

(b)

line

Log h	Series 1	Series 2	Series 3
-1.0	0.0	0.0	0.0
-0.5	0.5	0.75	0.6
0.0	1.0	1.25	1.0

Figure 5.31 (a) Convergence of section rotations in analysis of ad-hoc problem using distorted meshes. The error measure is E = \left\| \beta - \beta_{h} \right\|_{1} . (b) Convergence of gradient of vertical displacement in analysis of ad-hoc problem using distorted meshes. The error measure is E = \left\| \nabla w - \nabla w_{h} \right\|_{0} .

line

Condition	MITC7	MITC12	MITC4	MITC9	MITC16
Log h	-0.5	-1.8	0.0	-0.5	-1.8
Log h	1.5	0.5	1.5	0.5	1.5
Log h	2.5	1.8	2.5	1.8	2.5
Log h	3.0	2.0	3.0	2.0	3.0

(a)
(b)
Figure 5.32 Convergence of transverse shear strains in analysis of ad-hoc problem. The error measure is E = \left\| \gamma - \gamma_{h} \right\|_{0} . (a) Uniform meshes. (b) Distorted meshes.

General Shell Elements

Let us consider next the formulation of general shell elements that can be used to analyze very complex shell geometries and stress distributions. For this objective we need to generalize the preceding plate element formulation approach, much in the same way as we generalized the isoparametric beam element formulation from straight two-dimensional to curved three-dimensional beams. As in the case of the formulation of beam elements (see Section 5.4.1), we consider the displacement interpolation which leads to a pure displacement-based element (see S. Ahmad, B. M. Irons, and O. C. Zienkiewicz [A]), and we then modify the formulation so as to avoid shear and membrane locking.

The displacement interpolation is obtained by considering the geometry interpolation. Consider a general shell element with a variable number of nodes, q. Figure 5.33 shows a nine-node element for which q = 9. Using the natural coordinates r, s, and t, the Cartesian coordinates of a point in the element with q nodal points are, before and after deformations,

text_image

Shell midsurface (z, w) θz ex ey y, v x, u (a/2) 0Vn | at Gauss integration point = Σ ak/2 hk | at Gauss integration point 0Vn k 0Vk wk ak vk βk 0Vk uK αk r 0Vk1 t

Figure 5.33 Nine-node shell element; also, definition of orthogonal \overline{r} , \overline{s} , t axes for constitutive relations

text_image

0Vn t Top surface s s-coordinate line (r, t are constant) Gauss point r-coordinate line (s, t are constant) a/2 r Midsurface a/2 Bottom surface r, s, t = vectors tangent to r, s, t coordinate lines e7 = s × t / ||s × t||2; e5 = t × e7 / ||t × e7||2; e1 = t / ||t||2

Figure 5.33 (continued)

{ } ^ { \ell } x ( r , s , t ) = \sum _ { k = 1 } ^ { q } h _ { k } { } ^ { \ell } x _ { k } + \frac { t } { 2 } \sum _ { k = 1 } ^ { q } a _ { k } h _ { k } { } ^ { \ell } V _ { n x } ^ { k }

{ } ^ { \ell } y ( r , s , t ) = \sum _ { k = 1 } ^ { q } h _ { k } { } ^ { \ell } y _ { k } + \frac { t } { 2 } \sum _ { k = 1 } ^ { q } a _ { k } h _ { k } { } ^ { \ell } V _ { n y } ^ { k } \tag {5.107}

{ } ^ { \ell } z ( r , s , t ) = \sum _ { k = 1 } ^ { q } h _ { k } { } ^ { \ell } z _ { k } + \frac { t } { 2 } \sum _ { k = 1 } ^ { q } a _ { k } h _ { k } { } ^ { \ell } V _ { n z } ^ { k }

where the h_{k}(r, s) are the interpolation functions summarized in Fig. 5.4 and

^{e}x, ^{e}y, ^{e}z = \text{Cartesian coordinates of any point in the element}

^{\ell}x_{k}, ^{\ell}y_{k}, ^{\ell}z_{k}= Cartesian coordinates of nodal point k

a_{k} = thickness of shell in t direction at nodal point k

^{t}V_{nx}^{k}, ^{t}V_{ny}^{k}, ^{t}V_{nz}^{k} = \text{components of unit vector } ^{t}\mathbf{V}_{n}^{k} \text{ “normal” to the shell midsurface in direction } t \text{ at nodal point } k; \text{ we call } ^{t}\mathbf{V}_{n}^{k} \text{ the normal vector}^{7} \text{ or, more appropriately, the director vector, at nodal point } k

and the left superscript \ell denotes, as in the general beam formulation, the configuration of the element; i.e., \ell = 0 and 1 denote the original and final configurations of the shell element. Hence, using (5.107), the displacement components are

u (r, s, t) = \sum_ {k = 1} ^ {q} h _ {k} u _ {k} + \frac {t}{2} \sum_ {k = 1} ^ {q} a _ {k} h _ {k} V _ {n x} ^ {k}

v (r, s, t) = \sum_ {k = 1} ^ {q} h _ {k} v _ {k} + \frac {t}{2} \sum_ {k = 1} ^ {q} a _ {k} h _ {k} V _ {n y} ^ {k} \tag {5.108}

w (r, s, t) = \sum_ {k = 1} ^ {q} h _ {k} w _ {k} + \frac {t}{2} \sum_ {k = 1} ^ {q} a _ {k} h _ {k} V _ {n z} ^ {k}

^{7} We call ^{6}V_{n}^{k} the normal vector although it may not be exactly normal to the midsurface of the shell in the original configuration (see Example 5.32), and in the final configuration (e.g., because of shear deformations).

where V_{n}^{k} stores the increments in the direction cosines of ^{0}V_{n}^{k} ,

\mathbf {V} _ {n} ^ {k} = ^ {1} \mathbf {V} _ {n} ^ {k} - ^ {0} \mathbf {V} _ {n} ^ {k} \tag {5.109}

The components of V_{n}^{k} can be expressed in terms of rotations at the nodal point k; however, there is no unique way of proceeding. An efficient way is to define two unit vectors ^{0}V_{1}^{k} and ^{0}V_{2}^{k} that are orthogonal to ^{0}V_{n}^{k} :

{ } ^ { 0 } \mathbf { V } _ { 1 } ^ { k } = \frac { \mathbf { e } _ { y } \times { } ^ { 0 } \mathbf { V } _ { n } ^ { k } } { \| \mathbf { e } _ { y } \times { } ^ { 0 } \mathbf { V } _ { n } ^ { k } \| _ { 2 } } \tag {5.110a}

where e_{y} is a unit vector in the direction of the y-axis. (For the special case ^{0}V_{n}^{k} parallel to e_{y} , we may simply use ^{0}V_{1}^{k} equal to e_{z} .) We can now obtain ^{0}V_{2}^{k} ,

{ } ^ { 0 } \mathbf { V } _ { 2 } ^ { k } = { } ^ { 0 } \mathbf { V } _ { n } ^ { k } \times { } ^ { 0 } \mathbf { V } _ { 1 } ^ { k } \tag {5.110b}

Let \alpha_{k} and \beta_{k} be the rotations of the director vector ^0\mathbf{V}_n^k about the vectors ^0\mathbf{V}_1^k and ^0\mathbf{V}_2^k . We then have, because \alpha_{k} and \beta_{k} are small angles,

\mathbf {V} _ {n} ^ {k} = - ^ {0} \mathbf {V} _ {2} ^ {k} \alpha_ {k} + ^ {0} \mathbf {V} _ {1} ^ {k} \beta_ {k} \tag {5.111}

This relationship can readily be proven when ^{0}V_{1} = e_{x} , ^{0}V_{2} = e_{y} and ^{0}V_{n} = e_{z} , but since these vectors are tensors, the relationship must also hold in general (see Section 2.4). Substituting from (5.111) into (5.108), we thus obtain

u (r, s, t) = \sum_ {k = 1} ^ {q} h _ {k} u _ {k} + \frac {t}{2} \sum_ {k = 1} ^ {q} a _ {k} h _ {k} \left(- ^ {0} V _ {2 x} ^ {k} \alpha_ {k} + ^ {0} V _ {1 x} ^ {k} \beta_ {k}\right)

v (r, s, t) = \sum_ {k = 1} ^ {q} h _ {k} v _ {k} + \frac {t}{2} \sum_ {k = 1} ^ {q} a _ {k} h _ {k} \left(- ^ {0} V _ {2 y} ^ {k} \alpha_ {k} + ^ {0} V _ {1 y} ^ {k} \beta_ {k}\right) \tag {5.112}

w (r, s, t) = \sum_ {k = 1} ^ {q} h _ {k} w _ {k} + \frac {t}{2} \sum_ {k = 1} ^ {q} a _ {k} h _ {k} \left(- ^ {0} V _ {2 z} ^ {k} \alpha_ {k} + ^ {0} V _ {1 z} ^ {k} \beta_ {k}\right)

With the element displacements and coordinates defined in (5.112) and (5.107) we can now proceed as usual to evaluate the element matrices of a pure displacement-based element. The entries in the displacement interpolation matrix H of the shell element are given in (5.112), and the entries in the strain-displacement interpolation matrix can be calculated using the procedures already described in the formulation of the beam element (see Section 5.4.1).

To evaluate the strain-displacement matrix, we obtain from (5.112),

\left[ \begin{array}{l} \frac {\partial u}{\partial r} \\ \frac {\partial u}{\partial s} \\ \frac {\partial u}{\partial t} \end{array} \right] = \sum_ {k = 1} ^ {q} \left[ \begin{array}{l l l} \frac {\partial h _ {k}}{\partial r} \left[ 1 \quad t g _ {1 x} ^ {k} & t g _ {2 x} ^ {k} \right] \\ \frac {\partial h _ {k}}{\partial s} \left[ 1 \quad t g _ {1 x} ^ {k} & t g _ {2 x} ^ {k} \right] \\ h _ {k} \left[ 0 \quad g _ {1 x} ^ {k} \quad g _ {2 x} ^ {k} \right] \end{array} \right] \left[ \begin{array}{l} u _ {k} \\ \alpha_ {k} \\ \beta_ {k} \end{array} \right] \tag {5.113}

and the derivatives of v and w are given by simply substituting for u and x the variables v, y and w, z , respectively. In (5.113) we use the notation

\mathbf {g} _ {1} ^ {k} = - \frac {1}{2} a _ {k} ^ {0} \mathbf {V} _ {2} ^ {k}; \quad \mathbf {g} _ {2} ^ {k} = \frac {1}{2} a _ {k} ^ {0} \mathbf {V} _ {1} ^ {k} \tag {5.114}

To obtain the displacement derivatives corresponding to the Cartesian coordinates x, y, z , we use the standard transformation

\frac {\partial}{\partial \mathbf {x}} = \mathbf {J} ^ {- 1} \frac {\partial}{\partial \mathbf {r}} \tag {5.115}

where the Jacobian matrix J contains the derivatives of the coordinates x, y, z with respect to the natural coordinates r, s, t. Substituting from (5.113) into (5.115), we obtain

\left[ \begin{array}{l} \frac {\partial u}{\partial x} \\ \frac {\partial u}{\partial y} \\ \frac {\partial u}{\partial z} \end{array} \right] = \sum_ {k = 1} ^ {q} \left[ \begin{array}{l l l} \frac {\partial h _ {k}}{\partial x} & g _ {1 x} ^ {k} G _ {x} ^ {k} & g _ {2 x} ^ {k} G _ {x} ^ {k} \\ \frac {\partial h _ {k}}{\partial y} & g _ {1 x} ^ {k} G _ {y} ^ {k} & g _ {2 x} ^ {k} G _ {y} ^ {k} \\ \frac {\partial h _ {k}}{\partial z} & g _ {1 x} ^ {k} G _ {z} ^ {k} & g _ {2 x} ^ {k} G _ {z} ^ {k} \end{array} \right] \left[ \begin{array}{l} u _ {k} \\ \alpha_ {k} \\ \beta_ {k} \end{array} \right] \tag {5.116}

and the derivatives of v and w are obtained in an analogous manner. In (5.116) we have

\frac {\partial h _ {k}}{\partial x} = J _ {1 1} ^ {- 1} \frac {\partial h _ {k}}{\partial r} + J _ {1 2} ^ {- 1} \frac {\partial h _ {k}}{\partial s} \tag {5.117}

G _ {x} ^ {k} = t \left(J _ {1 1} ^ {- 1} \frac {\partial h _ {k}}{\partial r} + J _ {1 2} ^ {- 1} \frac {\partial h _ {k}}{\partial s}\right) + J _ {1 3} ^ {- 1} h _ {k}

where J_{ij}^{-1} is element (i,j) of \mathbf{J}^{-1} , and so on.

With the displacement derivatives defined in (5.116) we now directly assemble the strain-displacement matrix B of a shell element. Assuming that the rows in this matrix correspond to all six global Cartesian strain components, \epsilon_{xx} , \epsilon_{yy} , \ldots , \gamma_{zx} , the entries in B are constructed in the usual way (see Section 5.3), but then the stress-strain law must contain the shell assumption that the stress normal to the shell surface is zero. We impose therefore that the stress in that direction is zero. Thus, if \tau and \epsilon store the Cartesian stress and strain components, we use

\boldsymbol {\tau} = \mathbf {C} _ {\mathrm{sh}} \boldsymbol {\epsilon} \tag {5.118}

where

\boldsymbol {\tau} ^ {T} = \left[ \begin{array}{l l l l l l} \tau_ {x x} & \tau_ {y y} & \tau_ {z z} & \tau_ {x y} & \tau_ {y z} & \tau_ {z x} \end{array} \right]

\boldsymbol {\epsilon} ^ {T} = \left[ \begin{array}{l l l l l l} \epsilon_ {x x} & \epsilon_ {y y} & \epsilon_ {z z} & \gamma_ {x y} & \gamma_ {y z} & \gamma_ {z x} \end{array} \right]

\mathbf {C} _ {\mathrm{sh}} = \mathbf {Q} _ {\mathrm{sh}} ^ {T} \left(\frac {E}{1 - \nu^ {2}} \left[ \begin{array}{c c c c c c} 1 & \nu & 0 & 0 & 0 & 0 \\ & 1 & 0 & 0 & 0 & 0 \\ & & 0 & 0 & 0 & 0 \\ & & \frac {1 - \nu}{2} & & 0 & 0 \\ \text { Symmetric } & & & k \frac {1 - \nu}{2} & & 0 \\ & & & & & k \frac {1 - \nu}{2} \end{array} \right]\right) \mathbf {Q} _ {\mathrm{sh}} \tag {5.119}

and Q_{sh} represents a matrix that transforms the stress-strain law from an \overline{r} , \overline{s} , t Cartesian shell-aligned coordinate system to the global Cartesian coordinate system. The elements of

the matrix Q_{sh} are obtained from the direction cosines of the \overline{r} , \overline{s} , t coordinate axes measured in the x, y, z coordinate directions,

\mathbf {Q} _ {\mathrm{sh}} = \left[ \begin{array}{c c c c c c} l _ {1} ^ {2} & m _ {1} ^ {2} & n _ {1} ^ {2} & l _ {1} m _ {1} & m _ {1} n _ {1} & n _ {1} l _ {1} \\ l _ {2} ^ {2} & m _ {2} ^ {2} & n _ {2} ^ {2} & l _ {2} m _ {2} & m _ {2} n _ {2} & n _ {2} l _ {2} \\ l _ {3} ^ {2} & m _ {3} ^ {2} & n _ {3} ^ {2} & l _ {3} m _ {3} & m _ {3} n _ {3} & n _ {3} l _ {3} \\ 2 l _ {1} l _ {2} & 2 m _ {1} m _ {2} & 2 n _ {1} n _ {2} & l _ {1} m _ {2} + l _ {2} m _ {1} & m _ {1} n _ {2} + m _ {2} n _ {1} & n _ {1} l _ {2} + n _ {2} l _ {1} \\ 2 l _ {2} l _ {3} & 2 m _ {2} m _ {3} & 2 n _ {2} n _ {3} & l _ {2} m _ {3} + l _ {3} m _ {2} & m _ {2} n _ {3} + m _ {3} n _ {2} & n _ {2} l _ {3} + n _ {3} l _ {2} \\ 2 l _ {3} l _ {1} & 2 m _ {3} m _ {1} & 2 n _ {3} n _ {1} & l _ {3} m _ {1} + l _ {1} m _ {3} & m _ {3} n _ {1} + m _ {1} n _ {3} & n _ {3} l _ {1} + n _ {1} l _ {3} \end{array} \right] \tag {5.120}

where

\begin{array}{l} l _ {1} = \cos \left(\mathbf {e} _ {x}, \mathbf {e} _ {\overline {{r}}}\right); \quad m _ {1} = \cos \left(\mathbf {e} _ {y}, \mathbf {e} _ {\overline {{r}}}\right); \quad n _ {1} = \cos \left(\mathbf {e} _ {z}, \mathbf {e} _ {\overline {{r}}}\right) \\ l _ {2} = \cos \left(\mathbf {e} _ {x}, \mathbf {e} _ {s}\right); \quad m _ {2} = \cos \left(\mathbf {e} _ {y}, \mathbf {e} _ {s}\right); \quad n _ {2} = \cos \left(\mathbf {e} _ {z}, \mathbf {e} _ {s}\right) \tag {5.121} \\ l _ {3} = \cos \left(\mathbf {e} _ {x}, \mathbf {e} _ {t}\right); \quad m _ {3} = \cos \left(\mathbf {e} _ {y}, \mathbf {e} _ {t}\right); \quad n _ {3} = \cos \left(\mathbf {e} _ {z}, \mathbf {e} _ {t}\right) \\ \end{array}

and the relation in (5.119) corresponds to a fourth-order tensor transformation as described in Section 2.4.

It follows that in the analysis of a general shell the matrix Q_{sh} may have to be evaluated anew at each integration point that is employed in the numerical integration of the stiffness matrix (see Section 5.5). However, when special shells are considered and, in particular, when a plate is analyzed, the transformation matrix and the stress-strain matrix C_{sh} need only be evaluated at specific points and can then be employed repetitively. For example, in the analysis of an assemblage of flat plates, the stress-strain matrix C_{sh} needs to be calculated only once for each flat structural part.

In the above formulation the strain-displacement matrix is formulated corresponding to the Cartesian strain components, which can be directly established using the derivatives in (5.116) . Alternatively, we could calculate the strain components corresponding to coordinate axes aligned with the shell element midsurface and establish a strain-displacement matrix for these strain components, as we did in the formulation of the general beam element in Section 5.4.1. The relative computational efficiency of these two approaches depends on whether it is more effective to transform the strain components (which always differ at the integration points) or to transform the stress-strain law.

It is instructive to compare this shell element formulation with a formulation in which flat elements with a superimposed plate bending and membrane stress behavior are employed (see Section 4.2.3). To identify the differences, assume that the general shell element is used as a flat element in the modeling of a shell; then the stiffness matrix of this element could also be obtained by superimposing the plate bending stiffness matrix derived in (5.94) to (5.99) (see Example 5.29) and the plane stress stiffness matrix discussed in Section 5.3.1. Thus, in this case, the general shell element reduces to a plate bending element plus a plane stress element, but a computational difference lies in the fact that these element matrices are calculated by integrating numerically only in the r-s element midplanes, whereas in the shell element stiffness calculation numerical integration is also performed in the t-direction (unless the general formulation is modified for this special case).

We illustrate some of the above relations in the following example.

EXAMPLE 5.32: Consider the four-node shell element shown in Fig. E5.32.

(a) Develop the entries in the displacement interpolation matrix.
(b) Calculate the thickness at the midpoint of the element and give the direction in which this thickness is measured.

text_image

z 10 0Vₙ³ 1.2 3 45° 0Vₙ² 2 0.8√2 y 15 t s r 0Vₙ⁴ 0.8 4 0Vₙ¹ 45° 1 0.8√2 10 x

Figure E5.32 Four-node shell element

The shell element considered has varying thickness but in some respects can be compared with the plate element in Example 5.29.

The displacement interpolation matrix H is given by the relations in (5.112). The functions h_{k} are those of a four-node two-dimensional element (see Fig. 5.4 and Example 5.29). The director vectors ^{0}V_{n}^{k} are given by the geometry of the element:

{ } ^ { 0 } \mathbf { V } _ { n } ^ { 1 } = \left[ \begin{array} { c } 0 \\ - 1 / \sqrt { 2 } \\ 1 / \sqrt { 2 } \end{array} \right] ; \qquad { } ^ { 0 } \mathbf { V } _ { n } ^ { 2 } = \left[ \begin{array} { c } 0 \\ - 1 / \sqrt { 2 } \\ 1 / \sqrt { 2 } \end{array} \right] ; \qquad { } ^ { 0 } \mathbf { V } _ { n } ^ { 3 } = \left[ \begin{array} { c } 0 \\ 0 \\ 1 \end{array} \right] ; \qquad { } ^ { 0 } \mathbf { V } _ { n } ^ { 4 } = \left[ \begin{array} { c } 0 \\ 0 \\ 1 \end{array} \right]

Hence, ^{0}\mathbf{V}_{1}^{1} = ^{0}\mathbf{V}_{1}^{2} = ^{0}\mathbf{V}_{1}^{3} = ^{0}\mathbf{V}_{1}^{4} = \begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}

{ } ^ { 0 } \mathbf { V } _ { 2 } ^ { 1 } = { } ^ { 0 } \mathbf { V } _ { 2 } ^ { 2 } = \left[ \begin{array} { c } 0 \\ 1 / \sqrt { 2 } \\ 1 / \sqrt { 2 } \end{array} \right] ; \quad { } ^ { 0 } \mathbf { V } _ { 2 } ^ { 3 } = { } ^ { 0 } \mathbf { V } _ { 2 } ^ { 4 } = \left[ \begin{array} { c } 0 \\ 1 \\ 0 \end{array} \right]

Also, a_1 = a_2 = 0.8\sqrt{2};\quad a_3 = 1.2;\quad a_4 = 0.8

The above expressions give all entries in (5.112).

To evaluate the thickness at the element midpoint and the direction in which the thickness is measured, we use the relation

\left(\frac {a}{2}\right) ^ {0} \mathbf {V} _ {n} \Bigg | _ {\text { midpoint }} = \sum_ {k = 1} ^ {4} \frac {a _ {k}}{2} h _ {k} \Bigg | _ {r = s = 0} ^ {0} \mathbf {V} _ {n} ^ {k}

where a is the thickness and the director vector {}^0\mathbf{V}_n gives the direction sought. This expression gives

\frac {a}{2} ^ {0} \mathbf {V} _ {n} = \frac {0 . 8 \sqrt {2}}{4} \left[ \begin{array}{c} 0 \\ - 1 / \sqrt {2} \\ 1 / \sqrt {2} \end{array} \right] + \frac {1 . 2}{8} \left[ \begin{array}{l} 0 \\ 0 \\ 1 \end{array} \right] + \frac {0 . 8}{8} \left[ \begin{array}{l} 0 \\ 0 \\ 1 \end{array} \right] = \left[ \begin{array}{l} 0 \\ - 0. 2 \\ 0. 4 5 \end{array} \right]

which gives ^{0}V_{n}=\begin{bmatrix}0.0\\ -0.406\\ 0.914\end{bmatrix};\quad a=0.985

This shell element formulation clearly has an important attribute, namely, that any geometric shape of a shell can be directly represented. The generality is further increased if the formulation is extended to transition elements (similar to the extension for the isoparametric beam element discussed in Section 5.4.1). Figure 5.34 shows how shell transition

Figure 5.34 Use of shell transition elements