%%%%%%%%ABBD MATRIX MATLAB CODE
%%%%%%DONE BY JAGATH
%%% WITH REFERENCE OF Aerospace Structural Composites
[Link]/haftka/
%University of Florida
close all;clc;clear all;
% Sir i have taken the problem of "ORTHOTROPIC (PLANE
STRESS)"condition
% If i take only orthotropic i need to initalize
% E1,E2,E3,G12,G31,G23,?21,?31,?32(9=indipendent constants) all
these values this will become a
% complex problem so for the simplification purpose i have taken
% "ORTHOTROPIC (PLANE STRESS )" case(4=indipendent constants)
% Initial Constants for the given materials
E1 = 38.6e9; % Youngs Modulus for material 1
E2 = 8.27e9; % Youngs Modulus for material 2
G12 = 4.14e9; % Shear Modulus for the given materials
nu12 = 0.26; % Poisson's Ratio for the given materials ,assume
that 'nu'is'?' Poisson's Ratio symbol
t = 0.001; % individual ply thickness (m)
theta = 30;% assume theta= mathamatical angle representation '?'
ang = [theta; theta+90; theta+90; theta; theta; theta+90;
theta+90; theta;%ang is nothing but angle
45; -45; -45; 45; 45; -45; -45; 45];% I have taken 16 layers of
composite laminate
ply_ang = ang.*pi./180; % converting to rad
n = length(ang); % Number of plies can be calculated by using
this formula
th = ply_ang;
h = n*t;
% Checking whether the number of plies match the number ply
angles
% But in general the below if not required because i have given
already ''ply_ang = ang.*pi./180;''
if(n~=length(ply_ang))
display('Insufficient number of ply angles provided')
exit;
end
nu21 = E2*nu12/E1; % This relation taken from jhonson book
% Q-matrix
Q11 = E1/(1-nu12*nu21);
Q22 = E2/(1-nu12*nu21);
Q12 = (nu12*E2)/(1-nu12*nu21);
Q66 = G12;
% Material invariants (Us)
U1 = (1/8)*(3*Q11 + 3*Q22 + 2*Q12 + 4*Q66);
U2 = (1/2)*(Q11 - Q22);
U3 = (1/8)*(Q11 + Q22 - 2*Q12 - 4*Q66);
U4 = (1/8)*(Q11 + Q22 + 6*Q12 - 4*Q66);
U5 = (1/8)*(Q11 + Q22 - 2*Q12 + 4*Q66);
% Finding the z's for each ply
layers = n/2;
for i=1:2*layers+1
z(i) = t*(layers-(i-1));
end
z = -z;
% Vs where the invariants are defined
% t is thickness of each individual layer and 'th' is ply angle
V0A = n*t;% n is number of layers
%[V0A to V4A invariants for A matrix V0B to V4B invariants for B
matrix; V0B to V4B invariants for B matrix again V0D to V4D
invariants for D matrix;
V1A = sum(t.*cos(2.*th));
V2A = sum(t.*sin(2.*th));
V3A = sum(t.*cos(4.*th));
V4A = sum(t.*sin(4.*th));
V0B = 0;
V1B = 0;
V2B = 0;
V3B = 0;
V4B = 0;
for i = 1:n
V1B = V1B + t*((z(i)+z(i+1))/2)*cos(2*th(i));
V2B = V2B + t*((z(i)+z(i+1))/2)*sin(2*th(i));
V3B = V3B + t*((z(i)+z(i+1))/2)*cos(4*th(i));
V4B = V4B + t*((z(i)+z(i+1))/2)*sin(4*th(i));
end
V0D = (h^3)/12;
V1D = 0;
V2D = 0;
V3D = 0;
V4D = 0;
for i = 1:n
V1D = V1D + t*(z(i)^2 + z(i)*z(i+1) + z(i+1)^2)/3*cos(2*th(i));
V2D = V2D + t*(z(i)^2 + z(i)*z(i+1) + z(i+1)^2)/3*sin(2*th(i));
V3D = V3D + t*(z(i)^2 + z(i)*z(i+1) + z(i+1)^2)/3*cos(4*th(i));
V4D = V4D + t*(z(i)^2 + z(i)*z(i+1) + z(i+1)^2)/3*sin(4*th(i));
end
% A-matrix
A11 = U1*V0A + U2*V1A + U3*V3A;
A22 = U1*V0A - U2*V1A + U3*V3A;
A12 = U4*V0A - U3*V3A;
A66 = U5*V0A - U3*V3A;
A16 = 0.5*U2*V2A + U3*V4A;
A26 = 0.5*U2*V2A - U3*V4A;
A = [A11 A12 A16; A12 A22 A26; A16 A26 A66];
% B-matrix
B11 = U1*V0B + U2*V1B + U3*V3B;
B22 = U1*V0B - U2*V1B + U3*V3B;
B12 = U4*V0B - U3*V3B;
B66 = U5*V0B - U3*V3B;
B16 = 0.5*U2*V2B + U3*V4B;
B26 = 0.5*U2*V2B - U3*V4B;
B = [B11 B12 B16; B12 B22 B26; B16 B26 B66];
% D-matrix
D11 = U1*V0D + U2*V1D + U3*V3D;
D22 = U1*V0D - U2*V1D + U3*V3D;
D12 = U4*V0D - U3*V3D;
D66 = U5*V0D - U3*V3D;
D16 = 0.5*U2*V2D + U3*V4D;
D26 = 0.5*U2*V2D - U3*V4D;
D = [D11 D12 D16; D12 D22 D26; D16 D26 D66];
% Sir is there any possibility to add non-linear terms for ABBD
matrix which you already
% explained in the class ,i have written and understood up to
this but i am unable to underestand
% how to add the extra terms for the non-linearity case
ABBD=[A B;B D] % i have written all the A B B D values in a
single matrix that is ABBD
A
B
B
D
%Results of ABD MATRIX
%
A =
% 1.0e+08 *
% 2.8725
1.2813
0.2307
% 1.2813
2.8725
-0.2307
% 0.2307
-0.2307
1.5946
%
B =
% 1.0e+04 *
% -5.3270
5.3270
-9.2267
%
5.3270
-5.3270
9.2267
% -9.2267
9.2267
5.3270
%
D =
%1.0e+03 *
%6.1895
2.7335
0.6069
%2.7335
6.0664
-0.3772
%0.6069
-0.3772
3.4019
