Problem 5 Write a matlab code to solve for the laminate [0/45/-45/90/30/60/30/-6
ID: 2088940 • Letter: P
Question
Problem 5 Write a matlab code to solve for the laminate [0/45/-45/90/30/60/30/-60] Each lamina has a fiber volume fraction of 63% with properties as follows. Both materials are isotropic E-200 GPa, v-0.3 Em-7 Gpa Vn-0.2 Calculate the ABD matrix, assume cach layer is lmm thickness. Then for an applied load and moments of: N,-14x106 Nm Ny-3x106 Nm Ns=4 x 106 Nm M- 5x103 N M,-0 M-16x103 N. Find the resultant strains. Finally caleulate the stress field in each lamina. Sketch the stresses as a function of z for the laminate.Explanation / Answer
%% Code for ABD
clear; format; clc;
format short g;
Nplies = 6;
thetadt = [0 45 -45 90 30 60 30 -60]; % ply angles in degrees, from top
thetadb = fliplr(thetadt); % ply angles in degrees, from bottom
h_ply = 0.003; % SI units, meters
h = Nplies * h_ply ;
for i = 1:Nplies
zbar(i) = - (h + h_ply)/2 + i*h_ply;
end
%% Parameters
E1 = 200.e9 ; % Pa
nu12 = .3 ;
E2 = 7.e9 ; % Pa
G12 = 5.61e9 ; % Pa
nu21 = nu12 * E2 / E1 ;
a1 = -0.3e-6 ; % coefficients of thermal expansion
a2 = 28.1e-6 ;
deltaT = 1 ;
% Q matrix (material coordinates)
denom = 1 - nu12 * nu21 ;
Q11 = E1 / denom ;
Q12 = nu12 * E2 / denom ;
Q22 = E2 / denom ;
Q66 = G12 ;
Q = [ Q11 Q12 0; Q12 Q22 0; 0 0 Q66] ;
%Q=[20 .7 0; .7 2 0; 0 0 .7]
a = [a1 a2 0]' ;
S=inv(Q) ;
% Qbar matrices (laminate coordinates) and contributions to
% ABD matrices
A = zeros(3,3);
B = zeros(3,3);
D = zeros(3,3);
NT = zeros(3,1);
MT = zeros(3,1);
for i = 1:Nplies
theta = thetadb(i) * pi / 180; % ply i angle in radians, from bottom
m = cos(theta) ;
n = sin(theta) ;
T = [ m^2 n^2 2*m*n; n^2 m^2 -2*m*n; -m*n m*n (m^2 - n^2)];
Qbar = inv(T) * Q * (inv(T))' ;
Sbar = inv(Qbar);
abar = T' * a ;
A = A + Qbar * h_ply;
B = B + Qbar * h_ply * zbar(i);
D = D + Qbar * (h_ply * zbar(i)^2 + h_ply^3 / 12);
NT = NT + Qbar * abar * h_ply * deltaT ;
MT = MT + Qbar * abar * h_ply * zbar(i) * deltaT ;
end
%% Printing Results ABD
Qbar
Sbar
A
B
D
ABD = [A B; B D] ;
ABDinv = inv(ABD) ;
%%%%% Code for Stress %%%%%%%%%%%%%%%
function Laminate(In,Out)
clc
if nargin==1
Out='output.out';
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%% Reading Program Inputs from (Laminate.dat) %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[F1t,F2t,F1c,F2c,F6,Del_C,Del_T]=textread(In,'%f %f %f %f %f %f %f',1,'headerlines',2);
[Nx,Ny,Nxy,Mx,My,Mxy]=textread(In,'%f %f %f %f %f %f',1,'headerlines',10);
[th,t,E1,E2,G12,v12,Alpha_1,Alpha_2,Beta_1,Beta_2]=textread(In,'%f %f %f %f %f %f %f %f %f %f','headerlines',15);
N = [Nx ; Ny ; Nxy];
M = [Mx ; My ; Mxy];
A=0;B=0;D=0;NT=0;NH=0;MT=0;MH=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Setting the (h) Matrix %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h=zeros(size(t,1)+1,1);
h(1,1)=sum(t)/2;
for i=2:size(t,1)+1
h(i)=h(i-1)-t(i-1);
end
h
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Calculating (Qbar)'s, (A,B,& D) Matrices,NT, NH, MT, & MH %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k=1;
disp(['A, B, & D Matrices for [' num2str(th') '] Laminate:'])
for i=size(t,1):-1:1
Alpha = [Alpha_1(i);Alpha_2(i);0];
Beta = [Beta_1(i);Beta_2(i);0];
S12 =[1/E1(i,1) -v12(i,1)/E1(i,1) 0; ...
-v12(i,1)/E1(i,1) 1/E2(i,1) 0
0 0 1/G12(i,1)];
Q12 = inv(S12);
disp(['Lamina # ' int2str(k) ': 'num2str(th(i,1))])
eval(['Qbar' int2str(k) '= T_sig(-th(' int2str(i) '))*Q12*T_eps(th(' int2str(i) '))']);
Qbar=eval(['Qbar' int2str(k)]);
Sbar = inv(Qbar);
A = A - Qbar * (h(i+1,1) - h(i,1)); % Negative summation since i=size(t,1):-1:1
B = B - 1/2 * Qbar * (h(i+1,1)^2 - h(i,1)^2);
D = D - 1/3 * Qbar * (h(i+1,1)^3 - h(i,1)^3);
NTp(:,k) = Del_T * Qbar * T_eps(-th(i)) * Alpha * (h(i,1) - h(i+1,1));
NHp(:,k) = Del_C * Qbar * T_eps(-th(i)) * Beta * (h(i,1) - h(i+1,1));
MTp(:,k) = 1/2 * Del_T * Qbar * T_eps(-th(i)) * Alpha * (h(i,1)^2 - h(i+1,1)^2);
MHp(:,k) = 1/2 * Del_C * Qbar * T_eps(-th(i)) * Beta * (h(i,1)^2 - h(i+1,1)^2);
NT = NT + Del_T * Qbar * T_eps(-th(i)) * Alpha * (h(i,1) - h(i+1,1));
NH = NH + Del_C * Qbar * T_eps(-th(i)) * Beta * (h(i,1) - h(i+1,1));
MT = MT + 1/2 * Del_T * Qbar * T_eps(-th(i)) * Alpha * (h(i,1)^2 - h(i+1,1)^2);
MH = MH + 1/2 * Del_C * Qbar * T_eps(-th(i)) * Beta * (h(i,1)^2 - h(i+1,1)^2);
k=k+1;
end
Nbar = N + NT + NH;
Mbar = M + MT + MH;
ABD = [A B;B D];
abd = inv(ABD);
a = abd(1:3,1:3);
b = abd(1:3,4:6);
bT= abd(4:6,1:3);
d = abd(4:6,4:6);
eps0_k = abd * [Nbar ; Mbar];
eps0 = eps0_k(1:3);
k = eps0_k(4:6);
Alpha_bar = 1/Del_T * ( a * NT + b * MT);
Beta_bar = 1/Del_C * ( a * NH + b * MH);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% Setting the (z) Matrix %%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
z(1,1)=sum(t)/2;
for i=2:size(t,1)
z(i-1,1)=z(i-1);
z(i-1,2)=z(i-1)-t(i-1);
z(i,1) = z(i-1,2);
end
z(size(t,1),2) = -sum(t)/2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculating [sig]k & [eps]k for each ply %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
j=size(t,1);
for i=1:size(t,1)
Alpha_p = T_eps(-th(j)) * [Alpha_1(j);Alpha_2(j);0];
Beta_p = T_eps(-th(j)) * [Beta_1(j);Beta_2(j);0];
eval(['eps' int2str(i) '_upper = eps0 + z(' int2str(j) ',1)*k - Alpha_p * Del_T - Beta_p * Del_C;']);
eval(['eps' int2str(i) '_lower = eps0 + z(' int2str(j) ',2)*k - Alpha_p * Del_T - Beta_p * Del_C;']);
eval(['sigxy' int2str(i) '_upper = Qbar' int2str(i) ' * eps' int2str(i) '_upper;']);
eval(['sigxy' int2str(i) '_lower = Qbar' int2str(i) ' * eps' int2str(i) '_lower;']);
j=j-1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculating Tsai-Hill Number %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
j=size(t,1);
for i=1:size(t,1)
eval(['sig' int2str(i) '_upper = T_sig(th(' int2str(j) ')) * sigxy' int2str(i) '_upper;']);
eval(['sig = sig' int2str(i) '_upper']);
if sig(1) > 0
F1 = F1t;
else
F1 = F1c;
end
if sig(2) > 0
F2 = F2t;
else
F2 = F2c;
end
%%%%%%%%%%%%%%%%%
% Max Stress Test
%%%%%%%%%%%%%%%%%
if abs(sig(1))>F1 | abs(sig(2))>F2 | abs(sig(3))>F6
MSRes(i,:) = 'Failed';
else
MSRes(i,:) = 'Passed';
end
%%%%%%%%%%%%%%%%%
% Tsai Hill Test
%%%%%%%%%%%%%%%%%
eval(['TH' int2str(i) '_upper = (sig(1)/F1)^2 + (sig(2)/F2)^2 - (sig(1)*sig(2)/F1^2) + (sig(3)/F6)^2';]);
if eval(['TH' int2str(i) '_upper']) >= 1
THRes(i,:) = 'Failed';
else
THRes(i,:) = 'Passed';
end
eval(['sig' int2str(i) '_lower = T_sig(th(' int2str(j) ')) * sigxy' int2str(i) '_lower;']);
eval(['sig = sig' int2str(i) '_lower;']);
if sig(1) > 0
F1 = F1t;
else
F1 = F1c;
end
if sig(2) > 0
F2 = F2t;
else
F2 = F2c;
end
%%%%%%%%%%%%%%%%%
% Max Stress Test
%%%%%%%%%%%%%%%%%
if abs(sig(1))>F1 | abs(sig(2))>F2 | abs(sig(3))>F6
MSRes(i,:) = 'Failed';
else
MSRes(i,:) = 'Passed';
end
%%%%%%%%%%%%%%%%%
% Tsai Hill Test
%%%%%%%%%%%%%%%%%
eval(['TH' int2str(i) '_lower = (sig(1)/F1)^2 + (sig(2)/F2)^2 - (sig(1)*sig(2)/F1^2) + (sig(3)/F6)^2;']);
if eval(['TH' int2str(i) '_upper']) >= 1
THRes(i,:) = 'Failed';
else
THRes(i,:) = 'Passed';
end
j=j-1;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Outputing to File
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
OUT=fopen(Out,'w');
fprintf(OUT,'%s ',datestr(now));
fprintf(OUT,'Results for [ %s ] Laminate: ',num2str(th'));
MAT1 = [Nbar,N,NT,NH,Mbar,M,MT,MH];
fprintf(OUT,' Nbar N NT NH Mbar M MT MH ');
fprintf(OUT,' %+7.6e %+7.6e %+7.6e %+7.6e %+7.6e %+7.6e %+7.6e %+7.6e ',MAT1');
MAT2 = [eps0,k,Alpha_bar,Beta_bar];
fprintf(OUT,' eps0 k Alpha_bar Beta_bar ');
fprintf(OUT,' %+7.6e %+7.6e %+7.6e %+7.6e ',MAT2');
fprintf(OUT,' Ply epsxy_upper epsxy_lower sigxy_upper sigxy_lower sig12_upper sig12_lower NT NH MT MH TH_upper TH_lower ');
for i=1:size(t,1)
MAT3 = [i*ones(3,1),eval(['eps' int2str(i) '_upper']),eval(['eps' int2str(i) '_lower']),eval(['sigxy' int2str(i) '_upper']),eval(['sigxy' int2str(i) '_lower']),eval(['sig' int2str(i) '_upper']),eval(['sig' int2str(i) '_lower']),...
NTp(:,i),NHp(:,i),MTp(:,i),MHp(:,i),eval(['TH' int2str(i) '_upper'])*ones(3,1),eval(['TH' int2str(i) '_lower'])*ones(3,1)];
fprintf(OUT,'%d %+7.6e %+7.6e %+7.6e %+7.6e %+7.6e %+7.6e %+7.6e %+7.6e %+7.6e %+7.6e %+7.6e %+7.6e ',MAT3');
fprintf(OUT,' ');
MAT3g(3*i-2:3*i,:)=MAT3;
end
fprintf(OUT,' Ply Status(Tsai-Hill) Status(Max-Stress) ');
for i=1:size(th,1)
fprintf(OUT,'%d %s %s ',th(size(th,1)+1-i),THRes(i,:),MSRes(i,:));
end
fprintf(OUT,' Effective Moduli of Laminate: ');
Ex = 1/a(1,1)/sum(t);
Ey = 1/a(2,2)/sum(t);
Gxy = 1/a(3,3)/sum(t);
vxy = -a(2,1)/a(1,1);
fprintf(OUT,' Ex Ey Gxy vxy ');
fprintf(OUT,' %+7.6e %+7.6e %+7.6e %+7.6e ',Ex,Ey,Gxy,vxy);
fprintf(OUT,' ABD ');
fprintf(OUT,' %+7.6e %+7.6e %+7.6e %+7.6e %+7.6e %+7.6e ',ABD);
fprintf(OUT,' abd ');
fprintf(OUT,' %+7.6e %+7.6e %+7.6e %+7.6e %+7.6e %+7.6e ',abd);
for i=1:size(unique(th),1)
ang = unique(th);
x=find(th==ang(i));
x=sort(x);
fprintf(OUT,' Qbar%s: ',int2str(x(1)));
fprintf(OUT,' %+7.6e %+7.6e %+7.6e ',eval(['Qbar' int2str(x(1))]))
end
fclose(OUT);
clc;
% --------------------------------
% INTERNAL FUNCTIONS
% --------------------------------
function T_sig=T_sig(ang)
th=ang*pi/180;
m=cos(th);
n=sin(th);
T_sig=[ ...
m^2 n^2 2*m*n
n^2 m^2 -2*m*n
-m*n m*n m^2-n^2];
function T_eps=T_eps(ang)
th=ang*pi/180;
m=cos(th);
n=sin(th);
T_eps=[ ...
m^2 n^2 m*n
n^2 m^2 -m*n
-2*m*n 2*m*n m^2-n^2];