Vibration of 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. See Reference.

First, create an femodel object for modal structural analysis and include the geometry of a plate.

model = femodel(AnalysisType="structuralModal", ...
                Geometry="Plate10x10x1.stl");

Plot the geometry with the face labels.

pdegplot(model.Geometry,FaceLabels="on", ...
                        FaceAlpha=0.5);

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

model.MaterialProperties = ...
    materialProperties(YoungsModulus=200e9, ...
                       PoissonsRatio=0.3, ...
                       MassDensity=8000);

In this example, the only boundary condition is the zero $$z$$-displacement on the four edge faces. These edge faces have labels 1 through 4.

model.FaceBC(1:4) = faceBC(ZDisplacement=0);

Create and plot a mesh. Specify the target minimum edge length so that there is one row of elements per plate thickness.

model = generateMesh(model,Hmin=1.3);
figure 
pdemesh(model);
title("Mesh with Quadratic Tetrahedral Elements");

For comparison with the published values, load the reference frequencies in Hz.

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

Solve the problem for the specified frequency range. Define the upper limit as slightly larger than the highest reference frequency and the lower limit as slightly smaller than the lowest reference frequency.

maxFreq = 1.1*refFreqHz(end)*2*pi;
result = solve(model,FrequencyRange=[-0.1 maxFreq]);

Calculate frequencies in Hz.

freqHz = result.NaturalFrequencies/(2*pi);

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

tfreqHz = table(refFreqHz.',freqHz(1:10));
tfreqHz.Properties.VariableNames = {'Reference','Computed'};
disp(tfreqHz);

    Reference     Computed 
    _________    __________

          0      1.6471e-05
          0      6.7345e-05
          0      7.4219e-05
     45.897          44.887
     109.44          109.78
     109.44          109.81
     167.89          168.64
     193.59          193.74
     206.19          207.46
     206.19          207.49

You see good agreement between the computed and published frequencies.

Plot the third component ($$z$$-component) of the solution for the seven lowest nonzero-frequency modes.

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