Code covered by the BSD License

# ABD matrix

### Vijay Kumar Badagi (view profile)

to calculate ABD matrix in composite materials

ABD_mat_revised.m
```%%Program for Calculation of A, B and D matrices for a Uniform laminate.
clc;
close all;
clear all;
digits(32);
format short e;
disp('1  NCT 301')
disp('2  New Material')
comm=input('Enter the number correspond to your choice>> ');
n=input('Enter Number of Layers (Total Number)>> ');
disp('----------------------------------------------------------')
disp('please enter layer stacking, one example is provided below');
disp('If the Layer Stacking(LS) is [+-30 0]s enter " [30 -30 0 0 -30 30] " Do Not Forget to Use Brackets [] !!!')
disp('---------------------------------------------------------------------------------------------------------')
l=input('Enter LS >> ');
clc

%Mechanical Properties of NCT 301
if comm==1
E1=155E+9;
E2=12.1E+9;
E3=E2;
n21=0.248;
n31=n21;
n23=0.458;
G12=4.4E+09;
G13=G12;
G23=3.20E+09;

else
disp('Enter the values for extensional and shear modulus in GPa')
mm=1e+9;
E1=input('E1 = ');
E1=E1*mm;
E2= input('E2 = ');
E2=E2*mm;
E3= input('E3 = ');
E3=E3*mm;
n21= input('n21 = ');
n31= input('n31 = ');
n23= input('n23 = ');
G12= input('G12 = ');
G12=G12*mm;
G13= input('G13 = ');
G13=G13*mm;
G23= input('G23 = ');
G23=G23*mm;
end

t   = 125e-6;

all_ang = l;
k = {};

N=length(all_ang);
h=zeros(N+1);
h(N+1)=N*t/2;
lstf=zeros(3*N,3);

for k = 1 : N

%%Compliance Matrix elements
S = [1/E1     -n21/E1    -n31/E1    0         0       0;
-n21/E1   1/E2       -n23/E2    0         0       0;
-n31/E1   -n23/E2     1/E3      0         0       0;
0           0           0       1/G23     0       0;
0           0           0        0       1/G13    0;
0            0          0        0       0       1/G12];

%%Stiffness Matrix elements
C = inv(S);

%Reduced Stiffness
Q = zeros(3,3);
Q(1,1) = C(1,1) - (C(1,3))^2 / C(3,3);
Q(2,2) = C(2,2) - (C(2,3))^2 / C(3,3);
Q(1,2) = C(1,2) - C(1,3)*C(2,3) / C(3,3);
Q(2,1) = Q(1,2);
Q(3,3) = C(6,6);

%Transformatioin Matrix

m = cosd (all_ang(k));
n = sind (all_ang(k));

T = zeros(3,3);
T(1,1)=m^2;
T(1,2)=n^2;
T(1,3)=2*m*n;
T(2,1)=n^2;
T(2,2)=m^2;
T(2,3)=-2*m*n;
T(3,1)=-m*n;
T(3,2)=m*n;
T(3,3)=m^2-n^2;

%Transformed Stiffness
q=zeros(3,3);
q([1:2],[1:2])=Q([1:2],[1:2]);
q(3,3)=2*Q(3,3);

q=inv(T)*q*T;
for i = 1:3
q(i,3)=q(i,3)/2;
end

lstf([3*k-2:3*k],[1:3])=q([1:3],[1:3]);

h(k)=(k-N/2-1)*t;
end

% To Calculate A B and D Matrices for Uniform Laminate

A=zeros(3,3);
B=zeros(3,3);
D=zeros(3,3);

for i=1:3
for j=1:3
q([1:3],[1:3])=lstf([1:3],[1:3]);
A(i,j) = q(i,j) * (h(2) - h(1));
B(i,j) = 1/2*(q(i,j) * (h(2)^2 - h(1)^2));
D(i,j) = 1/3*(q(i,j) * (h(2)^3 - h(1)^3));

for k = 2 : N
q([1:3],[1:3]) = lstf( [3*k-2:3*k] , [1:3] );
A(i,j) = q(i,j) * (h(k+1) - h(k)) + A(i,j);
B(i,j) = 1/2*(q(i,j) * (h(k+1)^2 - h(k)^2)) + B(i,j);
D(i,j) = 1/3*(q(i,j) * (h(k+1)^3 - h(k)^3)) + D(i,j);
end
end
end

LamnStf=zeros(6,6);
a=zeros(3,3);
b=zeros(3,3);
d=zeros(3,3);
LamnStf([1:3],[1:3])=A([1:3],[1:3]);
LamnStf([4:6],[4:6])=D([1:3],[1:3]);
LamnStf([1:3],[4:6])=B([1:3],[1:3]);
LamnStf([4:6],[1:3])=B([1:3],[1:3]);
LamnCmp=inv(LamnStf);
a([1:3],[1:3])=LamnCmp([1:3],[1:3]);
b([1:3],[1:3])=LamnCmp([1:3],[4:6]);
d([1:3],[1:3])=LamnCmp([4:6],[4:6]);

if rem(length(all_ang),2)==0
A= A.*[1 1 0 ; 1 1 0 ; 0 0 1];
B= B.*[0 0 0 ; 0 0 0 ; 0 0 0];
D= D.*[1 1 1 ; 1 1 1 ; 1 1 1];

a= a.*[1 1 0 ; 1 1 0 ; 0 0 1];
b= b.*[0 0 0 ; 0 0 0 ; 0 0 0];
d= d.*[1 1 1 ; 1 1 1 ; 1 1 1];

end;

A
B
D

a
b
d

% To Calculate Force and Moment Resultants or Reference Strains or
% reference curvatures.

Nx= A(

```