bendtorsionelementstiffness.m

薄壁结构有限元计算

% 弯扭耦合分析单刚形成

% 输入剪切面积、剪切静矩、绕扭心的极惯性矩、布雷特扭转惯性矩
Axx=input('Please input the x-direction of shear-area of the section: ');
Ayy=input('Please input the y-direction of shear-area of the section: ');
Axy=input('Please input the mixed shear-area of the section: ');
Ssx=input('Please input the x-direction of shear-moment of the section: ');
Ssy=input('Please input the y-direction of shear-moment of the section: ');
Ip=input('Please input the polar inertia-moment of the section: ');
JB=input('Please input the Bredt torsion-constant of the section: ');

% 计算r矩阵
r=[Axx Axy Ssx;Axy Ayy Ssy;Ssx Ssy (Ip-JB)]\eye(3);

% 输入弹性模量、剪切模量、绕x轴的惯性矩、绕y轴的惯性矩、扇性惯性矩、圣维南扭转惯性矩
E=input('Please input the elastic modulus of the material: ');
G=input('Please input the shear modulus of the material: ');
Ix=input('Please input the x-direction of inertia-moment of the section: ');
Iy=input('Please input the y-direction of inertia-moment of the section: ');
Iw=input('Please input the flabellate inertia-moment of the section: ');
JS=input('Please input the StVenant torsion-constant of the section: ');

% 计算扭转惯性矩、mju常数、k常数
J=JB+JS;
mju=1/(1+r(3,3)*J);
mju1=(Ip-JB)/(Ip+JS);
k=sqrt(mju1*G*J/(E*Iw));

% 计算lambda矩阵
lambda=zeros(3,3);
lambda(:,1)=(E*Iy/G)*r(:,1);
lambda(:,2)=(E*Ix/G)*r(:,2);
lambda(:,3)=(E*Iw/(mju*G))*r(:,3);

% 计算dx,dy,g1,g2,g3常数
dx=-r(2,3)*J;
dy=-r(1,3)*J;
g1=dx+lambda(2,3)*k*k;
g2=dy+lambda(1,3)*k*k;
g3=1+lambda(3,3)*k*k;

% 输入单元长度
l=input('Please input the length of the element: ');

% 计算D矩阵
D=[1 0 0 0 0 0 0 0 0 0 0 0 0 0;
   0 0 0 0 1 0 0 0 0 0 0 0 0 0;
   0 0 0 0 0 0 0 0 1 0 0 0 0 0;
   0 0 0 -6*lambda(2,1) 0 -1 0 -6*lambda(2,2) 0 0 0 -dx -g1*k 0;
   0 1 0 6*lambda(1,1) 0 0 0 6*lambda(1,2) 0 0 0 dy g2*k 0;
   0 0 0 0 0 0 0 0 0 0 1 0 0 1;
   0 0 0 6*lambda(3,1) 0 0 0 6*lambda(3,2) 0 0 0 1 g3*k 0;
   1 l l*l l*l*l 0 0 0 0 0 0 0 0 0 0;
   0 0 0 0 1 l l*l l*l*l 0 0 0 0 0 0;
   0 0 0 0 0 0 0 0 1 l 0 0 0 0;
   0 0 0 -6*lambda(2,1) 0 -1 -2*l -3*l*l-6*lambda(2,2) 0 0 0 -dx -g1*k*cosh(k*l) -g1*k*sinh(k*l);
   0 1 2*l 3*l*l+6*lambda(1,1) 0 0 0 6*lambda(1,2) 0 0 0 dy g2*k*cosh(k*l) g2*k*sinh(k*l);
   0 0 0 0 0 0 0 0 0 0 1 l sinh(k*l) cosh(k*l);
   0 0 0 6*lambda(3,1) 0 0 0 6*lambda(3,2) 0 0 0 1 g3*k*cosh(k*l) g3*k*sinh(k*l)];

% 计算Ipx常数
Ipx=Ip-Ssx*Ssx/(Axx-Axy)-Ssy*Ssy/(Ayy-Axy);

% 计算zeta矩阵
zeta=zeros(3,3);
zeta(1,:)=lambda(1,:)+lambda(3,:)*Ssx/(Axx-Axy);
zeta(2,:)=lambda(2,:)+lambda(3,:)*Ssy/(Ayy-Axy);
zeta(3,:)=lambda(1,:)+lambda(2,:);

% 输入截面面积
A=input('Please input the area of the section: ');

% 计算单元弯扭耦合刚度矩阵K
k10=zeros(14,14);k10(10,10)=l;k1=E*A*(inv(D))'*k10*inv(D);

k20=zeros(14,14);k20(3,3)=4*l;k20(4,4)=12*l*l*l;k20(3,4)=6*l*l;k20(4,3)=6*l*l;k2=E*Iy*(inv(D))'*k20*inv(D);

k30=zeros(14,14);k30(7,7)=4*l;k30(8,8)=12*l*l*l;k30(7,8)=6*l*l;k30(8,7)=6*l*l;k3=E*Ix*(inv(D))'*k30*inv(D);

k40=zeros(14,14);k40(13,13)=k*k*k*sinh(2*k*l)/4-k*k*k*k*l/2;k40(14,14)=k*k*k*sinh(2*k*l)/4+k*k*k*k*l/2;
    k40(13,14)=k*k*k*sinh(k*l)*sinh(k*l)/2;k40(14,13)=k*k*k*sinh(k*l)*sinh(k*l)/2;k4=E*Iw*(inv(D))'*k40*inv(D);
    
k50=zeros(14,14);k50(12,12)=l;k50(13,13)=k*sinh(2*k*l)/4+k*k*l/2;k50(14,14)=k*sinh(2*k*l)/4-k*k*l/2;k50(12,13)=sinh(k*l);k50(13,12)=sinh(k*l);
    k50(12,14)=cosh(k*l)-1;k50(14,12)=cosh(k*l)-1;k50(13,14)=k*sinh(k*l)*sinh(k*l)/2;k50(14,13)=k*sinh(k*l)*sinh(k*l)/2;k5=G*J*(inv(D))'*k50*inv(D);
    
k60=zeros(14,14);k60(4,4)=36*lambda(3,1)*lambda(3,1)*l;k60(8,8)=36*lambda(3,2)*lambda(3,2)*l;k60(4,8)=36*lambda(3,1)*lambda(3,2)*l;k60(8,4)=36*lambda(3,1)*lambda(3,2)*l;
    k60(13,13)=lambda(3,3)*lambda(3,3)*(k*k*k*k*k*sinh(2*k*l)/4+k*k*k*k*k*k*l/2);k60(14,14)=lambda(3,3)*lambda(3,3)*(k*k*k*k*k*sinh(2*k*l)/4-k*k*k*k*k*k*l/2);
    k60(4,13)=6*lambda(3,1)*lambda(3,3)*k*k*sinh(k*l);k60(13,4)=6*lambda(3,1)*lambda(3,3)*k*k*sinh(k*l);k60(8,13)=6*lambda(3,2)*lambda(3,3)*k*k*sinh(k*l);
    k60(13,8)=6*lambda(3,2)*lambda(3,3)*k*k*sinh(k*l);k60(4,14)=6*lambda(3,1)*lambda(3,3)*k*k*(cosh(k*l)-1);k60(14,4)=6*lambda(3,1)*lambda(3,3)*k*k*(cosh(k*l)-1);
    k60(8,14)=6*lambda(3,2)*lambda(3,3)*k*k*(cosh(k*l)-1);k60(14,8)=6*lambda(3,2)*lambda(3,3)*k*k*(cosh(k*l)-1);
    k60(13,14)=lambda(3,3)*lambda(3,3)*k*k*k*k*k*sinh(k*l)*sinh(k*l)/2;k60(14,13)=lambda(3,3)*lambda(3,3)*k*k*k*k*k*sinh(k*l)*sinh(k*l)/2;k6=G*(Ipx-JB)*(inv(D))'*k60*inv(D);
    
k70=zeros(14,14);k70(4,4)=36*zeta(1,1)*zeta(1,1)*l;k70(8,8)=36*zeta(1,2)*zeta(1,2)*l;k70(13,13)=zeta(1,3)*zeta(1,3)*(k*k*k*k*k*sinh(2*k*l)/4+k*k*k*k*k*k*l/2);
    k70(14,14)=zeta(1,3)*zeta(1,3)*(k*k*k*k*k*sinh(2*k*l)/4-k*k*k*k*k*k*l/2);k70(4,8)=36*zeta(1,1)*zeta(1,2)*l;k70(8,4)=36*zeta(1,1)*zeta(1,2)*l;
    k70(4,13)=6*zeta(1,1)*zeta(1,3)*k*k*sinh(k*l);k70(13,4)=6*zeta(1,1)*zeta(1,3)*k*k*sinh(k*l);k70(8,13)=6*zeta(1,2)*zeta(1,3)*k*k*sinh(k*l);
    k70(13,8)=6*zeta(1,2)*zeta(1,3)*k*k*sinh(k*l);k70(4,14)=6*zeta(1,1)*zeta(1,3)*k*k*(cosh(k*l)-1);k70(14,4)=6*zeta(1,1)*zeta(1,3)*k*k*(cosh(k*l)-1);
    k70(8,14)=6*zeta(1,2)*zeta(1,3)*k*k*(cosh(k*l)-1);k70(14,8)=6*zeta(1,2)*zeta(1,3)*k*k*(cosh(k*l)-1);
    k70(13,14)=zeta(1,3)*zeta(1,3)*k*k*k*k*k*sinh(k*l)*sinh(k*l)/2;k70(14,13)=zeta(1,3)*zeta(1,3)*k*k*k*k*k*sinh(k*l)*sinh(k*l)/2;k7=G*(Axx-Axy)*(inv(D))'*k70*inv(D);
    
k80=zeros(14,14);k80(4,4)=36*zeta(2,1)*zeta(2,1)*l;k80(8,8)=36*zeta(2,2)*zeta(2,2)*l;k80(13,13)=zeta(2,3)*zeta(2,3)*(k*k*k*k*k*sinh(2*k*l)/4+k*k*k*k*k*k*l/2);
    k80(14,14)=zeta(2,3)*zeta(2,3)*(k*k*k*k*k*sinh(2*k*l)/4-k*k*k*k*k*k*l/2);k80(4,8)=36*zeta(2,1)*zeta(2,2)*l;k80(8,4)=36*zeta(2,1)*zeta(2,2)*l;
    k80(4,13)=6*zeta(2,1)*zeta(2,3)*k*k*sinh(k*l);k80(13,4)=6*zeta(2,1)*zeta(2,3)*k*k*sinh(k*l);k80(8,13)=6*zeta(2,2)*zeta(2,3)*k*k*sinh(k*l);
    k80(13,8)=6*zeta(2,2)*zeta(2,3)*k*k*sinh(k*l);k80(4,14)=6*zeta(2,1)*zeta(2,3)*k*k*(cosh(k*l)-1);k80(14,4)=6*zeta(2,1)*zeta(2,3)*k*k*(cosh(k*l)-1);
    k80(8,14)=6*zeta(2,2)*zeta(2,3)*k*k*(cosh(k*l)-1);k80(14,8)=6*zeta(2,2)*zeta(2,3)*k*k*(cosh(k*l)-1);
    k80(13,14)=zeta(2,3)*zeta(2,3)*k*k*k*k*k*sinh(k*l)*sinh(k*l)/2;k80(14,13)=zeta(2,3)*zeta(2,3)*k*k*k*k*k*sinh(k*l)*sinh(k*l)/2;k8=G*(Ayy-Axy)*(inv(D))'*k80*inv(D);
    
k90=zeros(14,14);k90(4,4)=36*zeta(3,1)*zeta(3,1)*l;k90(8,8)=36*zeta(3,2)*zeta(3,2)*l;k90(13,13)=zeta(3,3)*zeta(3,3)*(k*k*k*k*k*sinh(2*k*l)/4+k*k*k*k*k*k*l/2);
    k90(14,14)=zeta(3,3)*zeta(3,3)*(k*k*k*k*k*sinh(2*k*l)/4-k*k*k*k*k*k*l/2);k90(4,8)=36*zeta(3,1)*zeta(3,2)*l;k90(8,4)=36*zeta(3,1)*zeta(3,2)*l;
    k90(4,13)=6*zeta(3,1)*zeta(3,3)*k*k*sinh(k*l);k90(13,4)=6*zeta(3,1)*zeta(3,3)*k*k*sinh(k*l);k90(8,13)=6*zeta(3,2)*zeta(3,3)*k*k*sinh(k*l);
    k90(13,8)=6*zeta(3,2)*zeta(3,3)*k*k*sinh(k*l);k90(4,14)=6*zeta(3,1)*zeta(3,3)*k*k*(cosh(k*l)-1);k90(14,4)=6*zeta(3,1)*zeta(3,3)*k*k*(cosh(k*l)-1);
    k90(8,14)=6*zeta(3,2)*zeta(3,3)*k*k*(cosh(k*l)-1);k90(14,8)=6*zeta(3,2)*zeta(3,3)*k*k*(cosh(k*l)-1);
    k90(13,14)=zeta(3,3)*zeta(3,3)*k*k*k*k*k*sinh(k*l)*sinh(k*l)/2;k90(14,13)=zeta(3,3)*zeta(3,3)*k*k*k*k*k*sinh(k*l)*sinh(k*l)/2;k9=G*Axy*(inv(D))'*k90*inv(D);
    
Ke=k1+k2+k3+k4+k5+k6+k7+k8+k9;