Help with the Deflection of a Plate in MATLAB

Dear all,

I am trying to solve for the deflection of plate and to do this I am using finite differences in MATLAB. The plate is 6m by 6m and has a uniform load of 2 kN/m^2.

The equation to solve the plate is:

u(i+1,j) + u(i-1,j) + u(i,j+1) + u(i,j-1) - 4*u(i,j) = q/D  

The input data is:

 E = 3.3E+07 % Modulus of elasticity
 D = E*0.2^3/(12*(1-0.2^2)) % Flexural rigidity
 q = 2 % Uniform load
 L = 6;  % plate dimensions
 dx = 1;  % Step size x direction
 dy = 1;  % Setp size y direction
 x = 0:dx:L; % x direction vector
 y = 0:dy:L; % y direction vector
 nx = length(x);
 ny = length(y);

 % Boundary Conditions at the edges - deflection is zero
 u(:,1) = 0
 u(1,:) = 0
 u(:,7) = 0
 u(7,:) = 0
 

So I set up a nested for loop as follows:

for i = 2:nx-1
     for j = 2:ny-1
        u(i,j) = (-q/D + u(i+1,j) + u(i,j+1) + u(i-1,j) + u(i,j-1))/4   
     end
end

U = u*1000

From this the answer I get is:

         0         0         0         0         0         0         0
         0   -0.0218   -0.0273   -0.0286   -0.0290   -0.0291         0
         0   -0.0273   -0.0355   -0.0378   -0.0385   -0.0387         0
         0   -0.0286   -0.0378   -0.0407   -0.0416   -0.0419         0
         0   -0.0290   -0.0385   -0.0416   -0.0426   -0.0430         0
         0   -0.0291   -0.0387   -0.0419   -0.0430   -0.0433         0
         0         0         0         0         0         0         0

This is not the solution I am aiming for as the correct solution will be a symmetrical matrix - with the magnitude increasing towards the centre - if you imagine the deflection of a plate under a uniform load fixed at the edges - the centre will deflect the most.

I am hoping some of you intelligent people out there could help me.

Many thanks

Scott

Answer 1

[Torsten]

You are trying to solve the heat conduction equation with a homogeneous heat sink -q/D and boundary temperature 0. Is this the same equation that has to be solved for the deflection of a plate ?

Further, you have to solve a linear system of equations to get u(i,j). Doing a fixed-point iteration in the loop

for i = 2:nx-1
    for j = 2:ny-1
        u(i,j) = (-q/D + u(i+1,j) + u(i,j+1) + u(i-1,j) + u(i,j-1))/4   
    end
end

will not solve this linear system "in one go".

Try this code instead:

E = 3.3E+07; % Modulus of elasticity
D = E*0.2^3/(12*(1-0.2^2)); % Flexural rigidity
q = 2; % Uniform load
L = 6;  % plate dimensions
dx = 1;  % Step size x direction
dy = 1;  % Setp size y direction
x = 0:dx:L; % x direction vector
y = 0:dy:L; % y direction vector
nx = length(x);
ny = length(y);
% Boundary Conditions at the edges - deflection is zero
u = zeros(nx,ny);
u(:,1) = 0;
u(1,:) = 0;
u(:,7) = 0;
u(7,:) = 0;
itermax = 50;
error = 1;
iter = 0;
while iter < itermax & error > 1e-8     
  iter = iter + 1;
  for i = 2:nx-1
    for j = 2:ny-1
       u(i,j) = (-q/D + u(i+1,j) + u(i,j+1) + u(i-1,j) + u(i,j-1))/4 ;  
    end
  end
  error = 0;
  for i = 2:nx-1
    for j = 2:ny-1
       error = error + (u(i,j) - (-q/D + u(i+1,j) + u(i,j+1) + u(i-1,j) + u(i,j-1))/4)^2 ;  
    end
  end
  error = sqrt(error);
end

error

error = 9.8363e-09

iter

iter = 33

U = u*1000

U = 7×7

0 0 0 0 0 0 0 0 -0.0831 -0.1225 -0.1343 -0.1225 -0.0831 0 0 -0.1225 -0.1854 -0.2047 -0.1854 -0.1225 0 0 -0.1343 -0.2047 -0.2266 -0.2047 -0.1343 0 0 -0.1225 -0.1854 -0.2047 -0.1854 -0.1225 0 0 -0.0831 -0.1225 -0.1343 -0.1225 -0.0831 0 0 0 0 0 0 0 0

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

[John D'Errico]

@Torsten

Actually, thin plate deflection would be governed by del^4, not del^2.

[Torsten]

@John D'Errico

Thank you. I couldn't imagine physics would be that easy.

[Scott Banks]

@John D'Errico Yes, that's correct, John. That it why I have to repeat the loop procedure again, but instead of the variable 'u', I define a new variable 'z' for the second order finite difference and make it being equal to the results from 'u'.

I have the correct answers now.

[Scott Banks]

@Torsten

Thank you, Torsten, everything is good now.

Could you talk me through what you did? Like what is the purpose of the errror variable you have defined?

[Torsten]

You have to solve

u(i+1,j) + u(i-1,j) + u(i,j+1) + u(i,j-1) - 4*u(i,j) - q/D = 0

for 2<=i<=6, 2<=j<=6, and "error" is the 2-norm of the deviation of this vector from 0 during the iteration for u.