This example shows how to calculate the deflection of a structural plate acted on by a pressure loading.

The partial differential equation for a thin, isotropic plate with a pressure loading is

$${\nabla}^{2}(D{\nabla}^{2}w)=-p$$

where $$D$$ is the bending stiffness of the plate given by

$$D=\frac{E{h}^{3}}{12(1-{\nu}^{2})}$$

and $$E$$ is the modulus of elasticity, $$\nu $$ is Poisson's ratio, and $$h$$ is the plate thickness. The transverse deflection of the plate is $$w$$ and $$p$$ is the pressure load.

The boundary conditions for the clamped boundaries are $$w=0$$ and $${w}^{\prime}=0$$ where $${w}^{\prime}$$ is the derivative of $$w$$ in a direction normal to the boundary.

The Partial Differential Equation Toolbox™ cannot directly solve the fourth order plate equation shown above but this can be converted to the following two second order partial differential equations.

$${\nabla}^{2}w=v$$

$$D{\nabla}^{2}v=-p$$

where $$v$$ is a new dependent variable. However, it is not obvious how to specify boundary conditions for this second order system. We cannot directly specify boundary conditions for both $$w$$ and $${w}^{\prime}$$. Instead, we directly prescribe $${w}^{\prime}$$ to be zero and use the following technique to define $${v}^{\prime}$$ in such a way to insure that $$w$$ also equals zero on the boundary. Stiff "springs" that apply a transverse shear force to the plate edge are distributed along the boundary. The shear force along the boundary due to these springs can be written $$n\cdot D\nabla v=-kw$$ where $$n$$ is the normal to the boundary and $$k$$ is the stiffness of the springs. The value of $$k$$ must be large enough that $$w$$ is approximately zero at all points on the boundary but not so large that numerical errors result because the stiffness matrix is ill-conditioned. This expression is a generalized Neumann boundary condition supported by Partial Differential Equation Toolbox™

In the Partial Differential Equation Toolbox™ definition for an elliptic system, the $$w$$ and $$v$$ dependent variables are u(1) and u(2). The two second order partial differential equations can be rewritten as

$$-{\nabla}^{2}{u}_{1}+{u}_{2}=0$$

$$-D{\nabla}^{2}{u}_{2}=p$$

which is the form supported by the toolbox. The input corresponding to this formulation is shown in the sections below.

Create a pde model for a PDE with two dependent variables.

numberOfPDE = 2; model = createpde(numberOfPDE);

E = 1.0e6; % modulus of elasticity nu = .3; % Poisson's ratio thick = .1; % plate thickness len = 10.0; % side length for the square plate hmax = len/20; % mesh size parameter D = E*thick^3/(12*(1 - nu^2)); pres = 2; % external pressure

For a single square, the geometry and mesh are easily defined as shown below.

gdm = [3 4 0 len len 0 0 0 len len]'; g = decsg(gdm,'S1',('S1')');

Create a geometry entity.

geometryFromEdges(model,g);

Plot the geometry and display the edge labels for use in the boundary condition definition.

figure; pdegplot(model,'EdgeLabels','on'); ylim([-1,11]) axis equal title 'Geometry With Edge Labels Displayed';

The documentation on PDE coefficients shows the required formats for the a and c matrices. The most convenient form for c in this example is $${n}_{c}=3N$$ from the table where $$N$$ is the number of differential equations. In this example $$N=2$$. The $$c$$ tensor, in the form of an $$N\times N$$ matrix of $$2\times 2$$ submatrices is shown below.

$$\left[\begin{array}{cccc}c(1)& c(2)& \cdot & \cdot \\ \cdot & c(3)& \cdot & \cdot \\ \cdot & \cdot & c(4)& c(5)\\ \cdot & \cdot & \cdot & c(6)\end{array}\right]$$

The six-row by one-column c matrix is defined below. The entries in the full $$2\times 2$$ a matrix and the $$2\times 1$$ f vector follow directly from the definition of the two-equation system shown above.

c = [1 0 1 D 0 D]'; a = [0 0 1 0]'; f = [0 pres]'; specifyCoefficients(model,'m',0,'d',0,'c',c,'a',a,'f',f);

`k = 1e7; % spring stiffness`

Define distributed springs on all four edges.

bOuter = applyBoundaryCondition(model,'neumann','Edge',(1:4),... 'g',[0 0],'q',[0 0; k 0]);

`generateMesh(model, 'Hmax', hmax);`

The solution is calculated using the `solvepde`

function and the transverse deflection is plotted using the `pdeplot`

function. For comparison, the transverse deflection at the plate center is also calculated using an analytical solution to this problem.

res = solvepde(model); u = res.NodalSolution; numNodes = size(model.Mesh.Nodes,2); figure pdeplot(model,'XYData',u(1:numNodes),'Contour','on'); title 'Transverse Deflection'

numNodes = size(model.Mesh.Nodes,2); fprintf('Transverse deflection at plate center(PDE Toolbox) = %12.4e\n',... min(u(1:numNodes,1)));

Transverse deflection at plate center(PDE Toolbox) = -2.7632e-01

Compute analytical solution.

```
wMax = -.0138*pres*len^4/(E*thick^3);
fprintf('Transverse deflection at plate center(analytical) = %12.4e\n', wMax);
```

Transverse deflection at plate center(analytical) = -2.7600e-01