Accelerating the pace of engineering and science

# Partial Differential Equation Toolbox

## Clamped, Square Isotropic Plate With a Uniform Pressure Load

This example shows how to calculate the deflection of a structural plate acted on by a pressure loading using the Partial Differential Equation Toolbox™.

PDE and Boundary Conditions For A Thin Plate

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

where is the bending stiffness of the plate given by

and is the modulus of elasticity, is Poisson's ratio, and is the plate thickness. The transverse deflection of the plate is and is the pressure load.

The boundary conditions for the clamped boundaries are and where is the derivative of 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.

where 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 and . Instead, we directly prescribe to be zero and use the following technique to define in such a way to insure that 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 where is the normal to the boundary and is the stiffness of the springs. The value of must be large enough that 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 and dependent variables are u(1) and u(2). The two second order partial differential equations can be rewritten as

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

Problem Parameters

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


Geometry and Mesh

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

gdmTrans = [3 4 0 len len 0 0 0 len len];
sf = 'S1';
nsmTrans = 'S1';
g = decsg(gdmTrans', sf, nsmTrans');
[p, e, t] = initmesh(g, 'Hmax', hmax);


Boundary Conditions

b = @boundaryFileClampedPlate;

type boundaryFileClampedPlate

function [ q, g, h, r ] = boundaryFileClampedPlate( p, e, u, time )
%BOUNDARYFILECLAMPEDPLATE Boundary conditions for heatTransferThinPlateExample
%   [ q, g, h, r ] = BOUNDARYFILECLAMPEDPLATE( p, e, u, time ) returns the
%   Neumann BC (q, g) and Dirichlet BC (h, r) matrices for the
%   clampedSquarePlateExample example.
%   p is the point matrix returned from INITMESH
%   e is the edge matrix returned from INITMESH
%   u is the solution vector (used only for nonlinear cases)
%   time (used only for parabolic and hyperbolic cases)
%
%   See also PDEBOUND, INITMESH

%   Copyright 2012 The MathWorks, Inc.

N = 2;
ne = size(e,2);
% Apply a shear force along the boundary due to the transverse
% deflection of stiff, distributed springs
k = 1e7; % spring stiffness
q = [0 k 0 0]'*ones(1,ne);
g = zeros(N, ne);
h = zeros(N^2, 2*ne);
r = zeros(N, 2*ne);

end



Coefficient Definition

The documentation for assempde shows the required formats for the a and c matrices in the section titled "PDE Coefficients for System Case". The most convenient form for c in this example is from the table where is the number of differential equations. In this example . The tensor, in the form of an matrix of submatrices is shown below.

The six-row by one-column c matrix is defined below. The entries in the full a matrix and the 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];


Finite Element and Analytical Solutions

The solution is calculated using the assempde 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.

u = assempde(b,p,e,t,c,a,f);
numNodes = size(p,2);
pdeplot(p, e, t, 'xydata', u(1:numNodes), 'contour', 'on');
title 'Transverse Deflection'

numNodes = size(p,2);
fprintf('Transverse deflection at plate center(PDE Toolbox)=%12.4e\n', min(u(1:numNodes,1)));
% 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(PDE Toolbox)= -2.7563e-01
Transverse deflection at plate center(analytical)= -2.7600e-01