clear
clc
syms E1 E2 Mu12 Mu21 t G12
E1=135000;
E2=13000;
Mu12=0.38;
t1=input('请输单层厚度:');
t=t1;
G12=4300;
Mu21=Mu12*E1/E2;
a=0;
l=-1;
Q=zeros(6,6);
Q11=E1/(1-Mu12*Mu21)
Q22=E2/(1-Mu12*Mu21);
Q12=Mu21*E2/(1-Mu12*Mu21);
Q21=Q12;
Q66=G12;
n=input('请输入总铺层数:');
z=zeros();
if(ceil(n/2)-n/2==0)
m=n;
for k=1:m+1
z(k)=-(m/2-k+1)*t;
%if k==m
% z(k+1)=-(m/2-k)*t;
% end
end
elseif(ceil(n/2)-n/2>0)
m=n;
for k=1:m+2
l=l+1;
z(k)=-m/2*t+l*t;
if(k==m)
z((k+3)/2)=0;
l=l-1;
end
end
end
A=zeros(6,6);
B=zeros(6,6);
D=zeros(6,6);
QQ=zeros();
for i=1:n
x=input('请一层一层输入铺层角度:');
m=x*pi/180;
QQ(1,1,i)=Q11*(cos(m))^4+2*(Q12+2*Q66)*(sin(m))^2*(cos(m))^2+Q22*(sin(m))^4
QQ(1,2,i)=(Q11+Q22-4*Q66)*(sin(m))^2*(cos(m))^2+Q12*((sin(m))^4+(cos(m))^4);
QQ(2,2,i)=Q11*(sin(m))^4+2*(Q12+2*Q66)*(sin(m))^2*(cos(m))^2+Q22*(cos(m))^4;
QQ(1,6,i)=(Q11-Q12-2*Q66)*sin(m)*(cos(m))^3+(Q12-Q22+2*Q66)*(sin(m))^3*cos(m);
QQ(2,6,i)=(Q11-Q12-2*Q66)*(sin(m))^3*cos(m)+(Q12-Q22+2*Q66)*sin(m)*(cos(m))^3;
QQ(6,6,i)=(Q11+Q22-2*Q12-2*Q66)*(sin(m))^2*(cos(m))^2+Q66*((sin(m))^4+(cos(m))^4);
QQ(2,1,i)=QQ(1,2,i);QQ(6,1,i)=QQ(1,6,i);QQ(6,2,i)=QQ(2,6,i);
end
for i=1:6
for j=1:6
for k=2:n+1
f=QQ(i,j,k-1)*(z(k)-z(k-1));
A(i,j)=f+A(i,j);
end
end
end
disp('拉伸刚度矩阵:');
A
for i=1:6
for j=1:6
for k=2:n+1
g=0.5*QQ(i,j,k-1)*((z(k))^2-(z(k-1))^2);
B(i,j)=g+B(i,j);
end
end
end
disp('耦合刚度矩阵:');
B
for i=1:6
for j=1:6
for k=2:n+1
h=1/3*QQ(i,j,k-1)*((z(k))^3-(z(k-1))^3);
D(i,j)=h+D(i,j);
end
end
end
disp('弯曲刚度矩阵:');
D