% MEK4540 Autumn 2009
%
% Matlab script to assist solution of Assignment 2 
%
% Brian Hayman 2009
%
% This script is based on the one provided for Assignment 1, with the 
% addition of a routine to estimate the lateral displacement w(x,y) at a 
% point with co-ordinates (x,y) on a simply supported, rectangular plate
% with edges x=0,a and y=0,b, subjected to a uniform lateral pressure p.
% Here we have set a=500 mm, b=1000 mm, p=10 kPa.  This illustrates the
% use of "nested for-loops" to perform a double series summation.
%
% The material properties, laminate lay-ups and plate dimensions must be
% replaced with the correct data.
%
% The last part must be modified to calculate the bending and twisting
% moments, as well as the deflection, at the points of interest.  You will
% also need to find a way of establishing where the respective moments
% have their maximum absolute values (remembering that these parameters
% may have positive or negative values). This may be done within the
% Matlab solution or by trial-and-error.
%
% Finally some extra lines must be added to find the critical value of the
% in-plane compressive load Nx for elastic buckling.  Here you must try
% appropriate values of m and n in the buckling mode shapes.

clear all;    % Sets all variables to zero

% Material data in N, mm units (including MPa = N/mm^2)
EL=40000;
ET=7700;
vLT=0.25;
GLT=2800;
t=0.5;

% Lay-up(s)
% You can use one Matlab file for all lay-ups, just "comment out" with "%"
% all except one lay-up each time you run it.

% ply angles in radians; alternative is R=[45 0 -45]*rad;
R=[45 0 -45]*3.14159/180;      
% z-coordinates for top and bottom surfaces of plies as defined in AB&C
% Fig. 6-5
h=[-1.5*t -0.5*t 0.5*t 1.5*t];  

% Compliance (S-matrix)
S(1,1)=1/EL;
S(2,2)=1/ET;
S(1,2)=-vLT/EL;
S(3,3)=1/GLT;
S(2,1)=S(1,2);

% Stiffness (Q-matrix)
Q=inv(S);

Q(3,3)=2*Q(3,3);        % convert from engineering to tensor shear strains

% Initialise A,B,D
A(1:3,1:3)=0;
B(1:3,1:3)=0;
D(1:3,1:3)=0;

% Calculate T and Q for each ply 'k'. Assemble A, B and D 
for i=1:length(R)                 

% Transformation matrix (T-matrix)
    T(1,1)=cos(R(i))^2;
    T(2,2)=cos(R(i))^2;
    T(1,2)=sin(R(i))^2;
    T(2,1)=sin(R(i))^2;
    T(3,1)=-sin(R(i))*cos(R(i));
    T(3,2)=sin(R(i))*cos(R(i));
    T(1,3)=2*sin(R(i))*cos(R(i));
    T(2,3)=-2*sin(R(i))*cos(R(i));
    T(3,3)=cos(R(i))^2-sin(R(i))^2;
    
    Qk(:,:,i)=inv(T)*Q*T;             % Stiffness matrix for tensor strains, following rotation
    Qk(:,3,i)=Qk(:,3,i)/2;            % Convert back to engineering shear strains
    A=A+Qk(:,:,i)*(h(i+1)-h(i));                % Calculate A-matrix
    B=B+.5*Qk(:,:,i)*(h(i+1)^2-h(i)^2);         % Calculate B-matrix
    D=D+1/3*Qk(:,:,i)*(h(i+1)^3-h(i)^3);        % Calculate D-matrix
end

%Print A,B,D
A
B
D

% Assemble total stiffness matrix for the laminate
% This gives the relation between stresses and strains with engineering shear strains
Qtot=[[A B];[B D]]	% Places contributions A,B,D with [B D] in rows below [A B]
    
%***********************************************
% Plate analysis
%***********************************************

% Plate dimensions (mm) and applied pressure (MPa)
a=500;			
b=1000;
p=.01;         % Here we have input 10 kPa = 0.01 MPa.

% Bending and twisting stiffnesses from D-matrix
D11=D(1,1);
D12=D(1,2);
D22=D(2,2);
D66=D(3,3);
DA=2*(D12+2*D66);

% Set maximum values of m and n (taken here as 11 for both)
M=11;
N=11;

% Initialise
q(1:M,1:N)=0;	% Initialise matrix q(m,n) representing coefficients qmn
w(1:M,1:N)=0;	% Initialise matrix w(m,n) representing contribution to W from (m,n) term

W=0;			% Initialise displacement W

% Specify co-ordinates of point at which displ. is calculated
x=250;		
y=500;		% Here we have chosen the centre of the plate.

X=x/a;
Y=y/b;

% Use nested for-loops to calculate and sum contributions from respective terms in series
for m=1:2:M
     for n=1:2:N
         coeffw=sin(m*pi*X)*sin(n*pi*Y);
         q(m,n)=16*p/(m*n*pi*pi);
         mp4=(m*pi/a)^4;
         np4=(n*pi/b)^4;
         mnp4=(m*n*pi*pi/(a*b))^2;
         w(m,n)=coeffw*q(m,n)/(D11*mp4+DA*mnp4+D22*np4);
         W=W+w(m,n);
     end
end

% Print results
p
x
y
W