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 ().
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);
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.
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')
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 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;
Include the PDE COefficients in
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.
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');
For comparison with published values, load the reference frequencies in Hz from .
refFreqHz = [0 0 0 45.897 109.44 109.44 167.89 193.59 206.19 206.19];
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 .
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= 1.61, New conv eig= 0 Basis= 12, Time= 1.66, New conv eig= 1 Basis= 14, Time= 1.73, New conv eig= 1 Basis= 16, Time= 1.83, New conv eig= 1 Basis= 18, Time= 1.89, New conv eig= 1 Basis= 20, Time= 2.02, New conv eig= 1 Basis= 22, Time= 2.18, New conv eig= 1 Basis= 24, Time= 2.28, New conv eig= 2 Basis= 26, Time= 2.40, New conv eig= 4 Basis= 28, Time= 2.48, New conv eig= 5 Basis= 30, Time= 2.56, New conv eig= 11 Basis= 32, Time= 2.62, New conv eig= 11 End of sweep: Basis= 32, Time= 2.62, New conv eig= 11 Basis= 21, Time= 3.07, New conv eig= 0 Basis= 23, Time= 3.21, New conv eig= 0 Basis= 25, Time= 3.30, New conv eig= 0 Basis= 27, Time= 3.35, New conv eig= 1 End of sweep: Basis= 27, Time= 3.36, New conv eig= 1 Basis= 22, Time= 3.84, New conv eig= 0 Basis= 24, Time= 3.93, New conv eig= 0 Basis= 26, Time= 4.06, New conv eig= 0 Basis= 28, Time= 4.17, New conv eig= 1 End of sweep: Basis= 28, Time= 4.17, New conv eig= 1 Basis= 23, Time= 4.86, New conv eig= 0 Basis= 25, Time= 5.01, New conv eig= 0 End of sweep: Basis= 25, Time= 5.01, New conv eig= 0
Calculate the frequencies in Hz from the eigenvalues.
freqHz = sqrt(eVal(1:numEig))/(2*pi);
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 0 45.897 44.794 109.44 109.48 109.44 109.55 167.89 167.9 193.59 193.72 206.19 207.37 206.19 207.43
You see good agreement between the computed and published frequencies.
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
 National Agency for Finite Element Methods and Standards. The Standard NAFEMS Benchmarks. United Kingdom: NAFEMS, October 1990.