gistlib - calculate the fem modes and theanalytical modes for the rectangular plate simply supported on the four edges in matlab

To calculate the finite element method (FEM) modes and analytical modes for a rectangular plate simply supported on the four edges in MATLAB, you can follow these steps:

Define the plate geometry and material properties, such as length, width, thickness, Young's modulus, Poisson's ratio, and density.


main.m
% Plate dimensions
L = 1; % Length (m)
W = 0.5; % Width (m)
t = 0.01; % Thickness (m)

% Material properties
E = 2.07e11; % Young's modulus (Pa)
nu = 0.3; % Poisson's ratio
rho = 7850; % Density (kg/m^3)
204 chars
10 lines

Create a mesh of the plate using the meshgrid function and define the nodal coordinates and element connectivity.


main.m
% Mesh generation
x = linspace(0, L, 21); % X-coordinates
y = linspace(0, W, 11); % Y-coordinates
[X, Y] = meshgrid(x, y); % Meshgrid
nodes = [(1:numel(X))' X(:) Y(:)]; % Nodal coordinates
elements = delaunay(nodes(:, 2), nodes(:, 3)); % Element connectivity
259 chars
7 lines

Assemble the stiffness matrix (K) and mass matrix (M) using the FEM approach.


main.m
% Assemble stiffness and mass matrices
K = zeros(2*size(nodes, 1)); % Initialize stiffness matrix
M = zeros(2*size(nodes, 1)); % Initialize mass matrix
for e = 1:size(elements, 1) % Loop over elements
    % Element coordinates and dof
    xe = nodes(elements(e, :), 2); % X-coordinates
    ye = nodes(elements(e, :), 3); % Y-coordinates
    dof = [2*elements(e, 1)-1 2*elements(e, 1) ...
           2*elements(e, 2)-1 2*elements(e, 2) ...
           2*elements(e, 3)-1 2*elements(e, 3)];
    
    % Element stiffness matrix
    Ke = plateStiffness(E, nu, t, xe, ye);
    
    % Element mass matrix
    Me = plateMass(rho, t, xe, ye);
    
    % Assemble to global matrices
    K(dof, dof) = K(dof, dof) + Ke;
    M(dof, dof) = M(dof, dof) + Me;
end
749 chars
22 lines

The plateStiffness and plateMass functions are user-defined functions that calculate the element stiffness and mass matrices, respectively.

Solve the generalized eigenvalue problem K*U = M*U*D, where U is the matrix of eigenvectors and D is the diagonal matrix of eigenvalues (squared natural frequencies).


main.m
% Solve eigenvalue problem
[U, D] = eig(K, M);
[D, idx] = sort(diag(D)); % Sort eigenvalues in ascending order
U = U(:, idx); % Sort eigenvectors accordingly
158 chars
5 lines

Calculate the analytical modes by solving the governing equation of a rectangular plate, which is a fourth-order partial differential equation in terms of deflection (w). The analytical solutions are given as a product of sinusoidal functions of x and y and polynomials of a and b, where a and b are constants that depend on the length, width, and thickness of the plate (see Timoshenko and Woinowsky-Krieger, Theory of Plates and Shells, 1959). Here, we will calculate the first six modes (mode shapes) analytically.


main.m
% Analytical mode shapes (first six)
a = L; % Length
b = W; % Width
h = t; % Thickness
n = 1:6; % Mode numbers
m = repmat(n, 4, 1)';
p = repmat([1; 2; 1; 2], 1, 6)';
q = repmat([1; 1; 2; 2], 1, 6)';
A = h/2*sqrt((E/(1-nu^2))./(m.^2*pi.^4/a^2 + p.^2*pi.^2/b^2 + q.^2*pi.^2/a^2));
U_analytical = zeros(size(nodes, 1), numel(n));
for i = 1:numel(n)
    U_analytical(:, i) = A(:, i).*sin(m(:, i).*pi.*X(:)./a).*sin(p(:, i).*pi.*Y(:)./b);
end
438 chars
14 lines

Plot the first six FEM modes and analytical modes to compare their shapes and frequencies.


main.m
% Plot the first six modes
figure;
for i = 1:6
    % FEM mode
    subplot(2, 6, i);
    trisurf(elements, nodes(:, 2), nodes(:, 3), U(:, i));
    shading interp;
    colorbar;
    title(sprintf('FEM Mode %d\nFreq = %.2f Hz', i, sqrt(D(i))/(2*pi)));
    axis equal;
    axis([0 L 0 W]);
    
    % Analytical mode
    subplot(2, 6, i+6);
    trisurf(elements, nodes(:, 2), nodes(:, 3), U_analytical(:, i));
    shading interp;
    colorbar;
    title(sprintf('Analytical Mode %d\nFreq = %.2f Hz', i, sqrt((m(i,:).^2*pi.^4/a^4 + ...
          p(i,:).^2*pi.^2/b^4 + q(i,:).^2*pi.^2/a^4)*E*h^2/(12*rho*(1-nu^2)))));
    axis equal;
    axis([0 L 0 W]);
end
653 chars
23 lines

The resulting figure shows the first six FEM modes and analytical modes, with their frequencies listed in the titles. Note that the analytical modes are only valid for certain boundary conditions and aspect ratios, while the FEM modes are general and can handle arbitrary boundary conditions and geometries.

Modes

calculate the fem modes and theanalytical modes for the rectangular plate simply supported on the four edges in matlab

similar matlab code snippets

related categories