Matlab解复合材料层合板的问题.是matlab源代码

N = input('N=') ; EL=180; AL=-0.025E-6; AT=27E-6; ET=10; RL=0.3; RT=0.03; GLT=5; Q11=EL/(1-RL*RT); Q22=ET/(1-RL*RT); Q12=ET*RL/(1-RT*RL); Q66=GLT; Q=[Q11 Q12 0;Q12 Q22 0;0 0 Q66] Q11L=zeros(1,N+1); Q22L=zeros(1,N+1); Q12L=zeros(1,N+1); Q66L=zeros(1,N+1); Q16L=zeros(1,N+1); Q26L=zeros(1,N+1); j=3;k=3; T2=zeros(j,k,N+1); T1=zeros(j,k,N+1); for i=1:N+1; c(i)=cos((i-1)*pi/180); s(i)=sin((i-1)*pi/180); Q11L(i)=Q11*(c(i)^4)+2*(Q12+2*Q66)*s(i)^2*c(i)^2+Q22*s(i)^4; Q22L(i)=Q11*(s(i)^4)+2*(Q12+2*Q66)*s(i)^2*c(i)^2+Q22*c(i)^4; Q12L(i)=Q12*(c(i)^4+s(i)^4)+(Q11+Q22-4*Q66)*(s(i)^2)*(c(i)^2); Q66L(i)=Q66*(c(i)^4+s(i)^4)+(Q11+Q22-2*Q12-2*Q66)*s(i)^2*c(i)^2; Q16L(i)=(Q11-Q12-2*Q66)*c(i)^3*s(i)+(-Q22+Q12+2*Q66)*c(i)*s(i)^3; Q26L(i)=(Q11-Q12-2*Q66)*c(i)*s(i)^3+(-Q22+Q12+2*Q66)*c(i)^3*s(i); Q21L(i)=Q12L(i); Q61L(i)=Q16L(i); Q62L(i)=Q26L(i); T1(:,:,i)=[c(i)^2 s(i)^2 c(i)*s(i);s(i)^2 c(i)^2 -s(i)*c(i);-2*s(i)*c(i) 2*s(i)*c(i) c(i)^2-s(i)^2]; end for i=1:N+1; c(i)=cos((i-1)*pi/180); s(i)=sin((i-1)*pi/180); Q11L1(i)=Q11*(c(i)^4)+2*(Q12+2*Q66)*s(i)^2*c(i)^2+Q22*s(i)^4; Q22L1(i)=Q11*(s(i)^4)+2*(Q12+2*Q66)*s(i)^2*c(i)^2+Q22*c(i)^4; Q12L1(i)=Q12*(c(i)^4+s(i)^4)+(Q11+Q22-4*Q66)*(s(i)^2)*(c(i)^2); Q66L1(i)=Q66*(c(i)^4+s(i)^4)+(Q11+Q22-2*Q12-2*Q66)*s(i)^2*c(i)^2; Q16L1(i)=(Q11-Q12-2*Q66)*c(i)^3*s(i)+(-Q22+Q12+2*Q66)*c(i)*s(i)^3; Q26L1(i)=(Q11-Q12-2*Q66)*c(i)*s(i)^3+(-Q22+Q12+2*Q66)*c(i)^3*s(i); Q21L1(i)=Q12L1(i); Q61L1(i)=Q16L1(i); Q62L1(i)=Q26L1(i); T2(:,:,i)=[c(i)^2 s(i)^2 c(i)*s(i);s(i)^2 c(i)^2 -s(i)*c(i);-2*s(i)*c(i) 2*s(i)*c(i) c(i)^2-s(i)^2]; end j=3;k=3; n=N+1; A1=zeros(j,k,n); A2=zeros(j,k,n); for i=1:n; A1(:,:,i)=[ Q11L(i) Q12L(i) Q16L(i);Q21L(i) Q22L(i) Q26L(i);Q61L(i) Q62L(i) Q66L(i)]; A2(:,:,i)=[ Q11L1(i) Q12L1(i) Q16L1(i);Q21L1(i) Q22L1(i) Q26L1(i);Q61L1(i) Q62L1(i) Q66L1(i)]; T=T1+T2; A=A1+A2; a(:,:,i)=inv(A(:,:,i)); u(i)=-a(1,2,i)/a(1,1,i); end j=3;k=1;ALF=zeros(j,k,n);W=zeros(j,k,1); for i=1:n; ALF(1,1,i)=AL*c(i)^2+AT*s(i)^2; ALF(2,1,i)=AL*s(i)^2+AT*c(i)^2; ALF(3,1,i)=2*s(i)*c(i)*(AL-AT); alf(:,:,i)=-a(:,:,i)*T(:,:,i)*Q*ALF(:,:,i); a1(i)=alf(1,1,i); a2(i)=alf(2,1,i); a3(i)=alf(3,1,i); end t=0:90; subplot(2,2,1) plot(t,u) title('泊松比','FontWeight','bold') xlabel('铺层角度') ylabel('泊松比') subplot(2,2,2) plot(t,a1,'r'); title('x方向热膨胀系数','FontWeight','bold') xlabel('铺层角度') ylabel('膨胀系数') subplot(2,2,3) plot(t,a2,'g'); title('y方向热膨胀系数','FontWeight','bold') xlabel('铺层角度') ylabel('膨胀系数') subplot(2,2,4) plot(t,a3,'b'); title('剪切方向热膨胀系数','FontWeight','bold') xlabel('铺层角度') ylabel('膨胀系数')