MultiPhysicsVault/.raw/AFirstCourseInTheFiniteElementMethod/AFirstCourseInTheFiniteElementMethod_054.md

uniform thickness of the block needed such that the sensor is compressed no more than 0.05 in. Also make sure that the maximum stress from the maximum distortion energy failure theory is less than the yield strength of the material. Use a factor of safety of 1.5 on the stress only. The overall size of the block must fit in a 1.5-in.-high, 1-in.-wide, 1-in.-deep volume. The block should be made of steel.

text_image

1000 psi 1.5 in. 1 in. 1.5 in. t= (uniform) 1 in. Fixed base

Figure P11–16 S-shaped block

11.17 A device is to be hydraulically loaded to resist an upward force P ¼ 6000 lb as shown in Figure P11–17. Determine the thickness of the device such that the maximum deflection is 0.1 in. vertically and the maximum stress is less than the yield strength using a factor of safety of 2 (only on the stress). The device must fit in a space 7 in. high, 3 in. wide, and 2.3 in. deep. The top flange is bent vertically as shown, and the device is clamped to the floor. Use steel for the material.

text_image

1.8 P A A

text_image

6.2 0.80 2.3

text_image

3.0 0.80 0.80 Section A-A

Figure P11–17 Hydraulically loaded device

Introduction

In this chapter, we will begin by describing elementary concepts of plate bending behavior and theory. The plate element is one of the more important structural elements and is used to model and analyze such structures as pressure vessels, chimney stacks (Figure 1–5), and automobile parts. This description of plate bending is followed by a discussion of some commonly used plate finite elements. A large number of plate bending element formulations exist that would require a lengthy chapter to cover. Our purpose in this chapter is to present the derivation of the sti¤ness matrix for one of the most common plate bending finite elements and then to compare solutions to some classical problems from a variety of bending elements in the literature.

We finish the chapter with a solution to a plate bending problem using a computer program.

12.1 Basic Concepts of Plate Bending

A plate can be considered the two-dimensional extension of a beam in simple bending. Both beams and plates support loads transverse or perpendicular to their plane and through bending action. A plate is flat (if it were curved, it would become a shell). A beam has a single bending moment resistance, while a plate resists bending about two axes and has a twisting moment.

We will consider the classical thin-plate theory or Kirchho¤ plate theory [1]. Many of the assumptions of this theory are analogous to the classical beam theory or Euler–Bernoulli beam theory described in Chapter 4 and in Reference [2].

Basic Behavior of Geometry and Deformation

We begin the derivation of the basic thin-plate equations by considering the thin plate in the x-y plane and of thickness t measured in the z direction shown in Figure 12–1.

text_image

z 2c y q x 2b t

Figure 12–1 Basic thin plate showing transverse loading and dimensions

The plate surfaces are at z = \pm t / 2 , and its midsurface is \mathrm { a t } \ z = 0 . The assumed basic geometry of the plate is as follows: (1) The plate thickness is much smaller than its inplane dimensions b and c (that is, t \ll b \ \mathrm { o r } \ c ) . (If t is more than about one-tenth the span of the plate, then transverse shear deformation must be accounted for and the plate is then said to be thick.) (2) The deflection w is much less than the thickness t (that is, w / t \ll 1 ) .

Kirchhoff Assumptions

Consider a di¤erential slice cut from the plate by planes perpendicular to the x axis as shown in Figure 1 2 \mathrm { - } 2 ( \mathrm { a } ) . Loading q causes the plate to deform laterally or upward in the z direction, and the deflection w of point P is assumed to be a function of x and y only; that is, w = w ( x , y ) and the plate does not stretch in the z direction. A line a { - } b drawn perpendicular to the plate surfaces before loading remains perpendicular to the surfaces after loading [Figure 12–2(b)]. This is consistent with the Kirchho¤ assumptions as follows:

1. Normals remain normal. This implies that transverse shear strains \gamma _ { y z } = 0 and similarly \gamma _ { x z } = 0 . However, \gamma _ { x y } does not equal 0 ; right angles in the plane of the plate may not remain right angles after loading. The plate may twist in the plane.
2. Thickness changes can be neglected and normals undergo no extension. This means normal strain, \varepsilon _ { z } = 0 .
3. Normal stress \sigma _ { z } has no e¤ect on in-plane strains \varepsilon _ { x } and \varepsilon _ { y } in the stress-strain equations and is considered negligible.
4. Membrane or in-plane forces are neglected here, and the plane stress resistance can be superimposed later (that is, the constant-strain triangle behavior of Chapter 6 can be superimposed with the basic plate bending element resistance). That is, the in-plane deformations in the x and y directions at the midsurface are assumed to be zero; u ( x , y , 0 ) = 0 and v ( x , y , 0 ) = 0 .

text_image

z, w dx a P O t/2 t/2 b x, u z

(a)

text_image

w u = -zα a P z Midsurface α = ∂w/∂x O w b x, u

(b)
Figure 12–2 Differential slice of plate of thickness t (a) before loading and (b) displacements of point P after loading, based on Kirchhoff theory. Transverse shear deformation is neglected, and so right angles in the cross section remain right angles. Displacements in the y-z plane are similar

Based on the Kirchho¤ assumptions, any point P in Figure 12–2 has displacement in the x direction due to a small rotation a of

u = - z \alpha = - z \left(\frac {\partial w}{\partial x}\right) \tag {12.1.1}

and similarly the same point has displacement in the y direction of

v = - z \left(\frac {\partial w}{\partial y}\right) \tag {12.1.2}

The curvatures of the plate are then given as the rate of change of the angular displacements of the normals and are defined as

\kappa_ {x} = - \frac {\partial^ {2} w}{\partial x ^ {2}} \quad \kappa_ {y} = - \frac {\partial^ {2} w}{\partial y ^ {2}} \quad \kappa_ {x y} = - \frac {2 \partial^ {2} w}{\partial x \partial y} \tag {12.1.3}

The first of Eqs. (12.1.3) is used in beam theory [Eq. (4.1.1e)].

Using the definitions for the in-plane strains from Eq. (6.1.4), along with Eq. (12.1.3), the in-plane strain/displacement equations become

\varepsilon_ {x} = - z \frac {\partial^ {2} w}{\partial x ^ {2}} \quad \varepsilon_ {y} = - z \frac {\partial^ {2} w}{\partial y ^ {2}} \quad \gamma_ {x y} = - 2 z \frac {\partial^ {2} w}{\partial x \partial y} \tag {12.1.4a}

or using Eq. (12.1.3) in Eq.(12.1.4a), we have

\varepsilon_ {x} = - z \kappa_ {x} \quad \varepsilon_ {y} = - z \kappa_ {y} \quad \gamma_ {x y} = - z \kappa_ {x y} \tag {12.1.4b}

The first of Eqs. (12.1.4a) is used in beam theory [see Eq. (4.1.10)]. The others are new to plate theory.

text_image

z dy q dx τxz σx τxy x τyz t σy τyx (a)

text_image

Qx Mxy Mxy Mxy Mxy Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyy Mxyy Mxyy Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyy Mxyy Mxyy Mxyxy Mxyy Mxyy Mxyx Mxyx Mxyx Mxyxy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyx Mxyx Mxyx Mxyxy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyxy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyy Mxyy Mxyy Mxyy Mxyx Mxyx Mxyx Mxyy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyy Mxyy Mxyy Mxyy Mxyy Mxyxy Mxyx Mxyx Mxyx Mxyx Mxyy Mxy y Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxy y Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxy y Mxyy Mxy y Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyx Mxyx Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyx Mxyx Mxy y Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyx Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyx Mxy y Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyx Mxyy Mxy y Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyx Mxyy Mxyy Mxy y Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyx Mxy y Mxyy Mxy y Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyx Mxyx Mxyy Mxy y Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyx Mxyy Mxyy Mxyy Mxy y Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyx Mxyy Mxy y Mxyy Mxy y Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyx Mxyy Mxy y Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mxyy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy Mzy

Figure 12–3 Differential element of a plate with (a) stresses shown on the edges of the plate and (b) differential moments and forces

Stress= Strain Relations

Based on the third assumption above, the plane stress equations can be used to relate the in-plane stresses to the in-plane strains for an isotropic material as

\sigma_ {x} = \frac {E}{1 - \nu^ {2}} (\varepsilon_ {x} + \nu \varepsilon_ {y})

\sigma_ {y} = \frac {E}{1 - \nu^ {2}} (\varepsilon_ {y} + \nu \varepsilon_ {x}) \tag {12.1.5}

\tau_ {x y} = G \gamma_ {x y}

The in-plane normal stresses and shear stress are shown acting on the edges of the plate in Figure 12–3(a). Similar to the stress variation in a beam, these stresses vary linearly in the z direction from the midsurface of the plate. The transverse shear stresses \tau _ { y z } and \tau _ { x z } are also present, even though transverse shear deformation is neglected. As in beam theory, these transverse stresses vary quadratically through the plate thickness. The stresses of Eq. (12.1.5) can be related to the bending moments M _ { x } and M _ { y } and to the twisting moment M _ { x y } acting along the edges of the plate as shown in Figure 12–3(b).

The moments are actually functions of x and y and are computed per unit length in the plane of the plate. Therefore, the moments are

M _ {x} = \int_ {- t / 2} ^ {t / 2} z \sigma_ {x} d z \quad M _ {y} = \int_ {- t / 2} ^ {t / 2} z \sigma_ {y} d z \quad M _ {x y} = \int_ {- t / 2} ^ {t / 2} z \tau_ {x y} d z \tag {12.1.6}

The moments can be related to the curvatures by substituting Eqs. (12.1.4b) into Eqs. (12.1.5) and then using those stresses in Eq. (12.1.6) to obtain

M _ {x} = D (\kappa_ {x} + \nu \kappa_ {y}) \qquad M _ {y} = D (\kappa_ {y} + \nu \kappa_ {x}) \qquad M _ {x y} = \frac {D (1 - \nu)}{2} \kappa_ {x y} \tag {12.1.7}

where D = E t ^ { 3 } / [ 1 2 ( 1 - \nu ^ { 2 } ) ] is called the bending rigidity of the plate.

The maximum magnitudes of the normal stresses on each edge of the plate are located at the top or bottom at z = t / 2 . For instance, it can be shown that

\sigma_ {x} = \frac {6 M _ {x}}{t ^ {2}} \tag {12.1.8}

This formula is similar to the flexure formula \sigma _ { x } = M _ { x } c / I when applied to a unit width of plate and when c = t / 2 .

The governing equilibrium di¤erential equation of plate bending is important in selecting the element displacement fields. The basis for this relationship is the equilibrium di¤erential equations derived by the equilibrium of forces with respect to the z direction and by the equilibrium of moments about the x and y axes, respectively. These equilibrium equations result in the following di¤erential equations:

\frac {\partial Q _ {x}}{\partial x} + \frac {\partial Q _ {y}}{\partial y} + q = 0

\frac {\partial M _ {x}}{\partial x} + \frac {\partial M _ {x y}}{\partial y} - Q _ {x} = 0 \tag {12.1.9}

\frac {\partial M _ {y}}{\partial y} + \frac {\partial M _ {x y}}{\partial x} - Q _ {y} = 0

where q is the transverse distributed loading and Q _ { x } and Q _ { y } are the transverse shear line loads shown in Figure 12–3(b).

Now substituting the moment/curvature relations from Eq. (12.1.7) into the second and third of Eqs. (12.1.9), then solving those equations for Q _ { x } and Q _ { y } , and finally substituting the resulting expressions into the first of Eqs. (12.1.9), we obtain the governing partial di¤erential equation for an isotropic, thin-plate bending behavior as

D \left(\frac {\partial^ {4} w}{\partial x ^ {4}} + \frac {2 \partial^ {4} w}{\partial x ^ {2} \partial y ^ {2}} + \frac {\partial^ {4} w}{\partial y ^ {4}}\right) = q \tag {12.1.10}

From Eq. (12.1.10), we observe that the solution of thin-plate bending using a displacement point of view depends on selection of the single-displacement component w, the transverse displacement.

If we neglect the di¤erentiation with respect to the y coordinate, Eq. (12.1.10) simplifies to Eq. (4.1.1g) for a beam (where the flexural rigidity D of the plate reduces to EI of the beam when the Poisson e¤ect is set to zero and the plate width becomes unity).

Potential Energy of a Plate

The total potential energy of a plate is given by

U = \frac {1}{2} \int \left(\sigma_ {x} \varepsilon_ {x} + \sigma_ {y} \varepsilon_ {y} + \tau_ {x y} \gamma_ {x y}\right) d V \tag {12.1.11}

The potential energy can be expressed in terms of the moments and curvatures by substituting Eqs. (12.1.4b) and (12.1.6) in Eq. (12.1.11) as

U = \frac {1}{2} \int \left(M _ {x} \kappa_ {x} + M _ {y} \kappa_ {y} + M _ {x y} \kappa_ {x y}\right) d A \tag {12.1.12}

12.2 Derivation of a Plate Bending Element Stiffness Matrix and Equations

Numerous finite elements for plate bending have been developed over the years, and Reference [3] cites 88 di¤erent elements. In this section we will introduce only one element formulation, the basic 12-degrees-of-freedom rectangular element shown in Figure 12–4. For more details of this formulation and of various other formulations including triangular elements, see References [4–18].

text_image

n m z w_i, f_wi y θ_yi, f_θyi i x θ_xi, f_θxi j

Figure 12–4 Basic rectangular plate element with nodal degrees of freedom

The formulation will be developed consistently with the sti¤ness matrix and equations for the bar, beam, plane stress/strain, axisymmetric, and solid elements of previous chapters.

Step 1 Select Element Type

We will consider the 12-degrees-of-freedom flat-plate bending element shown in Figure 12–4. Each node has 3 degrees of freedom—a transverse displacement w in the z direction, a rotation \theta _ { x } about the x axis, and a rotation \theta _ { y } about the y axis.

The nodal displacement matrix at node i is given by

\{d _ {i} \} = \left\{ \begin{array}{l} w _ {i} \\ \theta_ {x i} \\ \theta_ {y i} \end{array} \right\} \tag {12.2.1}

where the rotations are related to the transverse displacement by

\theta_ {x} = + \frac {\partial w}{\partial y} \quad \theta_ {y} = - \frac {\partial w}{\partial x} \tag {12.2.2}

The negative sign on \theta _ { y } is due to the fact that a negative displacement w is required to produce a positive rotation about the y axis.

The total element displacement matrix is now given by

\{d \} = \{\underline {{d}} _ {i} \quad \underline {{d}} _ {j} \quad \underline {{d}} _ {m} \quad \underline {{d}} _ {n} \} ^ {T} \tag {12.2.3}

Step 2 Select the Displacement Function

Because there are 12 total degrees of freedom for the element, we select a 12-term polynomial in x and y as follows:

\begin{array}{l} w = a _ {1} + a _ {2} x + a _ {3} y + a _ {4} x ^ {2} + a _ {5} x y + a _ {6} y ^ {2} + a _ {7} x ^ {3} + a _ {8} x ^ {2} y \\ + a _ {9} x y ^ {2} + a _ {1 0} y ^ {3} + a _ {1 1} x ^ {3} y + a _ {1 2} x y ^ {3} \tag {12.2.4} \\ \end{array}

Equation (12.2.4) is an incomplete quartic in the context of the Pascal triangle (Figure 8–2). The function is complete up to the third order (ten terms), and a choice of two more terms from the remaining five terms of a complete quartic must be made. The best choice is the x ^ { 3 } y and x y ^ { 3 } terms as they ensure that we will have continuity in displacement among the interelement boundaries. (The x ^ { 4 } and y ^ { 4 } terms would yield discontinuities of displacement along interelement boundaries and so must be rejected. The x ^ { 2 } y ^ { 2 } term is alone and cannot be paired with any other terms and so is also rejected.) The function [Eq. (12.2.4)] also satisfies the basic di¤erential equation [Eq. (12.1.10)] over the unloaded part of the plate, although not a requirement in a minimum potential energy approximation.

Furthermore, the function allows for rigid-body motion and constant strain, as terms are present to account for these phenomena in a structure. However, interelement slope discontinuities along common boundaries of elements are not ensured.

To observe this discontinuity in slope, we evaluate the polynomial and its slopes along a side or edge (say, along side i \ – j , the x axis of Figure 12–4). We then obtain

w = a _ {1} + a _ {2} x + a _ {4} x ^ {2} + a _ {7} x ^ {3}

\frac {\partial w}{\partial x} = a _ {2} + 2 a _ {4} x + 3 a _ {7} x ^ {2} \tag {12.2.5}

\frac {\partial w}{\partial y} = a _ {3} + a _ {5} x + a _ {8} x ^ {2} + a _ {1 2} x ^ {3}

The displacement w is a cubic as used for the beam element, while the slope { \hat { o } } w / { \hat { o } } x is the same as in beam bending. Based on the beam element, we recall that the four constants a _ { 1 } , a _ { 2 } , a _ { 4 } , and a _ { 7 } can be defined by invoking the endpoint conditions of ( w _ { i } , w _ { j } , \theta _ { y i } , \theta _ { y j } ) . Therefore, w and { \partial w } / { \partial x } are completely defined along this edge. The normal slope { \hat { o } } w / { \hat { o } } y is a cubic in x. However, only two degrees of freedom remain for definition of this slope, while four constants ( a _ { 3 } , a _ { 5 } , a _ { 8 } . , and a _ { 1 2 } ) exist. This slope is then not uniquely defined, and a slope discontinuity occurs. Thus, the function for w is said to be nonconforming. The solution obtained from the finite element analysis using this element will not be a minimum potential energy solution. However, this element has proven to give acceptable results, and proofs of its convergence have been shown [8].

The constants a _ { 1 } through a _ { 1 2 } can be determined by expressing the 12 simultaneous equations linking the values of w and its slopes at the nodes when the coordinates take up their appropriate values. First, we write

\left\{ \begin{array}{c} w \\ + \frac {\partial w}{\partial y} \\ - \frac {\partial w}{\partial x} \end{array} \right\} = \left[ \begin{array}{c c c c c c c c c c c c} 1 & x & y & x ^ {2} & x y & y ^ {2} & x ^ {3} & x ^ {2} y & x y ^ {2} & y ^ {3} & x ^ {3} y & x y ^ {3} \\ 0 & 0 & + 1 & 0 & + x & + 2 y & 0 & + x ^ {2} & + 2 x y & + 3 y ^ {2} & + x ^ {3} & + 3 x y ^ {2} \\ 0 & - 1 & 0 & - 2 x & - y & 0 & - 3 x ^ {2} & - 2 x y & - y ^ {2} & 0 & - 3 x ^ {2} y & - y ^ {3} \end{array} \right]

\times \left\{ \begin{array}{l} a _ {1} \\ a _ {2} \\ a _ {3} \\ \vdots \\ a _ {1 2} \end{array} \right\} \tag {12.2.6}

or in simple matrix form the degrees of freedom matrix is

\{\psi \} = [ P ] \{a \} \tag {12.2.7}

where [P] is the 3 \times 12 first matrix on the right side of Eq. (12.2.6).

Next, we evaluate Eq. (12.2.6) at each node point as follows

\{d \} = \left\{ \begin{array}{c} w _ {i} \\ \theta_ {x i} \\ \theta_ {y i} \\ w _ {j} \\ \vdots \end{array} \right\} = \left[ \begin{array}{c c c c c c c c c c c c} 1 & x _ {i} & y _ {i} & x _ {i} ^ {2} & x _ {i} y _ {i} & y _ {i} ^ {2} & x _ {i} ^ {3} & x _ {i} ^ {2} y _ {i} & x _ {i} y _ {i} ^ {3} & y _ {i} ^ {3} & x _ {i} ^ {3} y _ {i} & x _ {i} y _ {i} ^ {3} \\ 0 & 0 & + 1 & 0 & + x _ {i} & + 2 y _ {i} & 0 & + x _ {i} ^ {2} & + 2 x _ {i} y _ {i} & + 3 y _ {i} ^ {2} & + x _ {i} ^ {3} & + 3 x _ {i} y _ {i} ^ {2} \\ \vdots & & & & & & & & & & \\ \vdots & & & & & & & & & & \\ \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots & \dots \end{array} \right]

\times \left\{ \begin{array}{l} a _ {1} \\ a _ {2} \\ \vdots \\ a _ {1 2} \end{array} \right\} \tag {12.2.8}

In compact matrix form, we express Eq. (12.2.8) as

\{d \} = [ C ] \{a \} \tag {12.2.9}

where [C] is the 12 \times 12 matrix on the right side of Eq. (12.2.8).

Therefore, the constants (a 's) can be solved for by

\{a \} = [ C ] ^ {- 1} \{d \} \tag {12.2.10}

Equation (12.2.7) can now be expressed as

\{\psi \} = [ P ] [ C ] ^ {- 1} \{d \} \tag {12.2.11}

or

\{\psi \} = [ N ] \{d \} \tag {12.2.12}

where [N]=[P][C]^{-1} is the shape function matrix. A specific form of the shape functions N_{i}, N_{j}, N_{m} , and N_{n} is given in Reference [9].

Step 3 Define the Strain (Curvature)/ Displacement and Stress (Moment)/Curvature Relationships

The curvature matrix, based on the curvatures of Eq.(12.1.3), is

\{\kappa \} = \left\{ \begin{array}{l} \kappa_ {x} \\ \kappa_ {y} \\ \kappa_ {x y} \end{array} \right\} = \left\{ \begin{array}{c} - 2 a _ {4} - 6 a _ {7} x - 2 a _ {8} y - 6 a _ {1 1} x y \\ - 2 a _ {6} - 2 a _ {9} x - 6 a _ {1 0} y - 6 a _ {1 2} x y \\ - 2 a _ {5} - 4 a _ {8} x - 4 a _ {9} y - 6 a _ {1 1} x ^ {2} - 6 a _ {1 2} y ^ {2} \end{array} \right\} \tag {12.2.13}

or expressing Eq. (12.2.13) in matrix form, we have

\{\kappa \} = [ Q ] \{a \} \tag {12.2.14}

where ½Q is the coe‰cient matrix multiplied by the \boldsymbol { a ^ { \prime } \mathrm { s } } in Eq. (12.2.13). Using Eq. (12.2.10) for \{ a \} , we express the curvature matrix as

\{\kappa \} = [ B ] \{d \} \tag {12.2.15}

where [ B ] = [ Q ] [ C ] ^ { - 1 } ð12:2:16Þ

is the gradient matrix.

The moment/curvature matrix for a plate is given by

\{M \} = \left\{ \begin{array}{l} M _ {x} \\ M _ {y} \\ M _ {x y} \end{array} \right\} = [ D ] \left\{ \begin{array}{l} \kappa_ {x} \\ \kappa_ {y} \\ \kappa_ {x y} \end{array} \right\} = [ D ] [ B ] \{d \} \tag {12.2.17}

where the ½D matrix is the constitutive matrix given for isotropic materials by

[ D ] = \frac {E t ^ {3}}{1 2 (1 - \nu^ {2})} \left[ \begin{array}{c c c} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & \frac {1 - \nu}{2} \end{array} \right] \tag {12.2.18}

and Eq. (12.2.15) has been used in the final expression for Eq. (12.2.17).

Step 4 Derive the Element Stiffness Matrix and Equations

The sti¤ness matrix is given by the usual form of the sti¤ness matrix as

[ k ] = \iint [ B ] ^ {T} [ D ] [ B ] d x d y \tag {12.2.19}

where [ B ] is defined by Eq. (12.2.16) and [ D ] is defined by Eq. (12.2.18). The sti¤ness matrix for the four-noded rectangular element is of order 1 2 \times 1 2 . A specific expression for ½k is given in References [4] and [5].

The surface force matrix due to distributed loading q acting per unit area in the z direction is obtained using the standard equation

\left\{F _ {s} \right\} = \iint \left[ N _ {s} \right] ^ {T} q d x d y \tag {12.2.20}

For a uniform load q acting over the surface of an element of dimensions 2 b \times 2 c Eq. (12.2.20) yields the forces and moments at node i as

\left\{ \begin{array}{l} f _ {w i} \\ f _ {\theta x i} \\ f _ {\theta y i} \end{array} \right\} = 4 q c b \left\{ \begin{array}{c} 1 / 4 \\ - c / 1 2 \\ b / 1 2 \end{array} \right\} \tag {12.2.21}

with similar expressions at nodes j , m , and n. We should note that a uniform load yields applied couples at the nodes as part of the work-equivalent load replacement, just as was the case for the beam element (Section 4.4).