复合材料力学层合板计算铺层应力以及最先一层失效载荷程序

function Main
sj=dlmread('input.txt'); %读入数据
A=0;B=0;D=0; % A 面内刚度矩阵 B 耦合刚度矩阵 D 弯曲刚度矩阵
syms R; % R 强度比
zl=sj(1,1); % zl 材料种类
cs=sj(zl+2,1); % cs 层数
t=sj(cs+zl+3,cs+1)-sj(cs+zl+3,1); % t 总厚度
for m=1:zl % Q1(:,:,m) 第m层正轴刚度
v12=sj(1+m,2)*sj(1+m,3)/sj(1+m,1);
Q1(:,:,m)=Z_G(sj(1+m,1),sj(1+m,2),v12,sj(1+m,3),sj(1+m,4)); % 计算每层正轴刚度
end
for m=1:cs % 计算每层偏轴刚度
sj(zl+2+m,2)=sj(zl+2+m,2)*pi/180;
[Q(:,:,m),L1(:,m),L2(:,m)]=P_G(Q1(:,:,sj(zl+2+m,1)),sj(zl+2+m,2),sj(sj(zl+2+m,1)+1,5),...
sj(sj(zl+2+m,1)+1,6),sj(sj(zl+2+m,1)+1,7),sj(sj(zl+2+m,1)+1,8));
end
for m1=1:cs % 计算 A B D
A=A+Q(:,:,m1)*(sj(cs+zl+3,m1+1)-sj(cs+zl+3,m1))/t;
B=B+Q(:,:,m1)*(sj(cs+zl+3,m1+1)^2-sj(cs+zl+3,m1)^2)/2/t;
D=D+Q(:,:,m1)*(sj(cs+zl+3,m1+1)^3-sj(cs+zl+3,m1)^3)/3/t;
end
q=[A B % q 总刚度
B D];
s=q^(-1); % s 总柔度
Ew=s*[sj(zl+cs+5,1);sj(zl+cs+5,2);sj(zl+cs+5,3);sj(zl+cs+5,4);sj(zl+cs+5,5);sj(zl+cs+5,6)];
Ez=vpa(Ew*R); % Ez 外力引起的中面应变
for m=1:cs
C(:,m)=Q(:,:,m)*([Ez(1);Ez(2);Ez(3)]); %C(:,m) 第m层中面应力
end
for m=1:cs %带入Hill-蔡强度理论求屈服载荷
C(:,m)=Z_H(C(:,m),sj(zl+2+m,2));
Hx(:,m)=solve(C(1,m)^2/sj(cs+6,1)^2-C(1,m)*C(2,m)/sj(cs+6,1)^2+C(2,m)^2/sj(cs+6,2)^2+C(3,m)^2/sj(cs+6,3)^2-1);
end
Hx=double(Hx);
Nx=Hx(2,:); %计算结果取正值
if Hx(1,1)>Hx(2,1)
Nx=Hx(1,:);
end
[Nmin,zb]=min(Nx); %计算最小屈服载荷
fid=fopen('output.txt','wt'); %以下为输出结果到文本
fprintf(fid,'一.求单层刚度\n');
for m=1:cs
fprintf(fid,'\nQ%d:',m);
for m1=1:3
fprintf(fid,'\n');
for m2=1:3
fprintf(fid,'%10.5f ',Q(m1,m2,m));
end
end
end
fprintf(fid,'\n 求总刚Q和柔度S \n Q:');
for m1=1:6
fprintf(fid,'\n');
for m2=1:6
fprintf(fid,'%10.5f ',q(m1,m2));
end
end
fprintf(fid,'\n 柔度S \n S:');
for m1=1:6
fprintf(fid,'\n');
for m2=1:6
fprintf(fid,'%20.15f ',s(m1,m2));
end
end
fprintf(fid,'\n 求中面应变\n Ez: \n');
for m1=1:3
o=char(vpa(Ez(m1),6));

fprintf(fid,'%s\n',o);
end
fprintf(fid,'\n 求各层应力 C \n');
for m=1:cs
fprintf(fid,'\nC%d:\n'

,m);
for m1=1:3
o=char(vpa(C(m1,m),6));
fprintf(fid,'%s\n',o);
end
end
fprintf(fid,'\n 带入Hill-蔡强路理论求屈服载荷 \n');
for m1=1:cs
fprintf(fid,'(Nx/t)%d: ',m1);
fprintf(fid,'%10.5f\n',Nx(m1));
end
fprintf(fid,'\n 则最先一层失效载荷为 \n (Nx/t)%d: ',zb);
fprintf(fid,'\n %10.5f \n',Nmin);
fclose(fid);

end

function Q=Z_G(E1,E2,V12,V21,G12) %计算正轴刚度
Q(1,1)=E1/(1-V21*V12);
Q(2,2)=E2/(1-V12*V21);
Q(1,2)=E2*V21/(1-V12*V21);
Q(2,1)=Q(1,2);
Q(3,3)=G12;
end

function [q,c,d]=P_G(Q,a,e,f,g,h) %计算偏轴刚度、正轴转换为偏轴的转换矩阵
T=[cos(a)*cos(a) sin(a)*sin(a) 2*sin(a)*cos(a)
sin(a)*sin(a) cos(a)*cos(a) -2*sin(a)*cos(a)
-sin(a)*cos(a) sin(a)*cos(a) cos(a)*cos(a)-sin(a)*sin(a)
];
q=T^(-1)*Q*(T^(-1))';
c=T'*[e;f;0];
d=T'*[g;h;0];
end

function b=Z_H(e,c) %偏轴转换到正轴的转换矩阵
a=-c;
T=[cos(a)*cos(a) sin(a)*sin(a) 2*sin(a)*cos(a)
sin(a)*sin(a) cos(a)*cos(a) -2*sin(a)*cos(a)
-sin(a)*cos(a) sin(a)*cos(a) cos(a)*cos(a)-sin(a)*sin(a)
];
b=T'*e;
end