% Script for Assignment 1, MEK4540 Autumn 2009
% Brian Hayman 2009
%
% The material properties and laminate lay-ups must be replaced with the correct data.
%
clear all;    % Sets all variables to zero

% Material data
EL=40000;
ET=7700;
vLT=0.25;
GLT=2800;
t=0.5;

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


% 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


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

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

    % 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

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]