MATLAB Examples

# Vibration of a Square Plate

This example shows how to calculate the vibration modes and frequencies of a 3-D simply supported, square, elastic plate. The dimensions and material properties of the plate are taken from a standard finite element benchmark problem published by NAFEMS, FV52 ([1]).

## Create a PDE Model

The first step in solving any 3-D PDE problem is to create a PDE Model. This is a container that holds the number of equations, geometry, mesh, PDE coefficients, and boundary conditions for your PDE. The equations of linear elasticity have three components, so the number of equations in this model is three.

N = 3; model = createpde(N); 

## Construct the Geometry

Import an STL file of a simple plate model using the importGeometry function. This function reconstructs the faces, edges and vertices of the model. It can merge some faces and edges, so the numbers can differ from those of the parent CAD model.

importGeometry(model,'Plate10x10x1.stl'); 

Plot the geometry and turn on face labels. You will need the face labels to define the boundary conditions.

figure hc = pdegplot(model,'FaceLabels','on'); hc(1).FaceAlpha = 0.5; title('Plate with Face Labels') 

## Examine the Elasticity Equations

As explained in 3-D Linear Elasticity Equations in Toolbox Form, you can express the elasticity equations for the deflection of a linear isotropic solid as a three-component system

where when there are no body forces, and is the symmetric matrix

The symbol means the entry is symmetric. Here is the elastic modulus, is Poisson's ratio, and

You can create the c coefficient for a linear isotropic solid using the elasticityC3D function, which is included in your software (see 3-D Linear Elasticity Equations in Toolbox Form).

## Define the Coefficients in Toolbox Syntax

Define the elastic modulus of steel, Poisson's ratio, and the material density.

E = 200e9; % Modulus of elasticity in Pascals nu = .3; % Poisson's ratio m = 8000; % Material density in kg/m^3 

Incorporate these coefficients in toolbox syntax.

c = elasticityC3D(E,nu); a = 0; 

## Specify PDE Coefficients

Include the PDE COefficients in model.

specifyCoefficients(model,'m',m,'d',0,'c',c,'a',a,'f',0); 

## Define the Boundary Conditions

In this example, the only boundary conditions are that the four edge faces have zero displacement in the direction. The edge faces have labels 1 through 4. All other boundary conditions, by default, are free Neumann boundaries.

Set the boundary conditions.

applyBoundaryCondition(model,'mixed','Face',1:4,'u',0,'EquationIndex',3); 

## Create a Mesh

Create a mesh that uses 10-node tetrahedral elements with quadratic interpolation functions. This element type is significantly more accurate than the linear interpolation (four-node) elements, particularly in elasticity analyses that involve bending.

hmax = 1.2; % Maximum element length for a moderately fine mesh generateMesh(model,'Hmax',hmax,'GeometricOrder','quadratic'); figure pdeplot3D(model); title('Mesh with Quadratic Tetrahedral Elements'); 

## Calculate the Solution

For comparison with published values, load the reference frequencies in Hz from [1].

refFreqHz = [0 0 0 45.897 109.44 109.44 167.89 193.59 206.19 206.19]; 

Use solvepdeeig to calculate all eigensolutions in a specified range. Define the upper value of the range slightly as slightly larger than the highest frequency from [1].

maxLam = (1.1*refFreqHz(end)*2*pi)^2; r = [-.1 maxLam]; result = solvepdeeig(model,r); eVec = result.Eigenvectors; eVal = result.Eigenvalues; numEig = size(eVal,1); 
 Basis= 10, Time= 0.28, New conv eig= 0 Basis= 12, Time= 0.29, New conv eig= 1 Basis= 14, Time= 0.30, New conv eig= 1 Basis= 16, Time= 0.30, New conv eig= 1 Basis= 18, Time= 0.32, New conv eig= 1 Basis= 20, Time= 0.33, New conv eig= 1 Basis= 22, Time= 0.33, New conv eig= 1 Basis= 24, Time= 0.34, New conv eig= 2 Basis= 26, Time= 0.35, New conv eig= 6 Basis= 28, Time= 0.36, New conv eig= 4 Basis= 30, Time= 0.36, New conv eig= 10 Basis= 32, Time= 0.37, New conv eig= 10 Basis= 34, Time= 0.38, New conv eig= 10 Basis= 36, Time= 0.42, New conv eig= 10 Basis= 38, Time= 0.43, New conv eig= 15 End of sweep: Basis= 38, Time= 0.43, New conv eig= 15 Basis= 25, Time= 0.47, New conv eig= 0 Basis= 27, Time= 0.47, New conv eig= 0 Basis= 29, Time= 0.48, New conv eig= 1 End of sweep: Basis= 29, Time= 0.48, New conv eig= 1 Basis= 26, Time= 0.53, New conv eig= 0 End of sweep: Basis= 26, Time= 0.53, New conv eig= 0 

Calculate the frequencies in Hz from the eigenvalues.

freqHz = sqrt(eVal(1:numEig))/(2*pi); 

## Compare the Solution to Reference

Compare the reference and computed frequencies for the lowest 10 modes. The lowest three mode shapes correspond to rigid-body motion of the plate. Their frequencies are close to zero.

fprintf('Frequency(Reference 1),Hz Frequency(solvepdeeig),Hz\n'); numToPrint = min(length(freqHz),length(refFreqHz)); for i=1:numToPrint fprintf('%10.5g %10.5g\n', ... refFreqHz(i),freqHz(i)); end 
Frequency(Reference 1),Hz Frequency(solvepdeeig),Hz 0 0 0 0 0 2.6043e-05 45.897 44.794 109.44 109.5 109.44 109.53 167.89 168 193.59 193.72 206.19 207.41 206.19 207.45 

You see good agreement between the computed and published frequencies.

## Plot the Mode Shapes

For plotting and analysis, create a PDEResults object from the solution. Plot the third () component of the solution for the seven lowest nonzero-frequency modes.

h = figure; h.Position = [100,100,900,600]; for i = 4:numToPrint subplot(4,2,i-3); pdeplot3D(model,'ColorMapData',eVec(:,3,i)); axis equal title(sprintf(['Mode=%d, z-displacement\n', ... 'Frequency(Hz): Ref=%g FEM=%g'], ... i,refFreqHz(i),freqHz(i))); end 

## Reference

[1] National Agency for Finite Element Methods and Standards. The Standard NAFEMS Benchmarks. United Kingdom: NAFEMS, October 1990.