Constant_DQM.zip_ABAQUS 复合层_MATLAB ABAQUS_abaqus write_matlab a

function [freq]=dqmvirabtionssss(a,b,n,E1,E2,v12,G12,T0,T1,t,p) %振动算例2 文献Non-linear vibrations of variable stiffness composite laminated plates对照 简支 %n是节点数，t是板厚，p是密度(kg/m^3) %节点选取 x0(1)=r/2; x0(2)=0.0001+r/2; x0(n-1)=r-0.0001; x0(n)=r; for i=3:n-2 x0(i)=r/2*[1+(i-1)/(n-1)]; end y0(1)=0; y0(2)=0.0001; y0(n-1)=theta-0.0001; y0(n)=theta; for j=3:n-2 y0(j)=[(j-1)/(n-1)]*theta; end %原始权系数矩阵 for i=1:n for j=1:n xp1=x0([1:i-1 i+1:n]); xp2=x0([1:j-1 j+1:n]); if i==j A1(i,j)=sum(1./(x0(i)-xp1')); else A1(i,j)=prod(x0(i)-xp1')/((x0(i)-x0(j))*prod(x0(j)-xp2')); end end end for i=1:n for j=1:n yp1=y0([1:i-1 i+1:n]); yp2=y0([1:j-1 j+1:n]); if i==j B1(i,j)=sum(1./(y0(i)-yp1')); else B1(i,j)=prod(y0(i)-yp1')/((y0(i)-y0(j))*prod(y0(j)-yp2')); end end end A2=A1^2; A3=A1^3; A4=A1^4; B2=B1^2; B3=B1^3; B4=B1^4; %内点 Aa1n=A1(2:n-1,2:n-1); Aa2n=A2(2:n-1,2:n-1); Aa3n=A3(2:n-1,2:n-1); Aa4n=A4(2:n-1,2:n-1); Ba1n=B1(2:n-1,2:n-1); Ba2n=B2(2:n-1,2:n-1); Ba3n=B3(2:n-1,2:n-1); Ba4n=B4(2:n-1,2:n-1); %方程中的权系数矩阵 for i=1:n-2 for j=1:n-2 for k=1:n-2 Ab1n((k-1)*(n-2)+i,(k-1)*(n-2)+j)=Aa1n(i,j); Ab2n((k-1)*(n-2)+i,(k-1)*(n-2)+j)=Aa2n(i,j); Ab3n((k-1)*(n-2)+i,(k-1)*(n-2)+j)=Aa3n(i,j); Ab4n((k-1)*(n-2)+i,(k-1)*(n-2)+j)=Aa4n(i,j); Bb1n(k-1+(i-1)*(n-2)+1,k-1+(j-1)*(n-2)+1)=Ba1n(i,j); Bb2n(k-1+(i-1)*(n-2)+1,k-1+(j-1)*(n-2)+1)=Ba2n(i,j); Bb3n(k-1+(i-1)*(n-2)+1,k-1+(j-1)*(n-2)+1)=Ba3n(i,j); Bb4n(k-1+(i-1)*(n-2)+1,k-1+(j-1)*(n-2)+1)=Ba4n(i,j); end end end %变系数求解 v21=v12*E2/E1; Q11=E1/(1-v12*v21); Q12=v12*E2/(1-v12*v21); Q22=E2/(1-v12*v21); Q66=G12; %系数求导 z=-4*t:t:4*t; syms x y theta1a=(2*(T1-T0)/r)*(x-r/2)+T0;%纤维变化方式 thetaa=[-theta1a,theta1a,-theta1a,theta1a,theta1a,-theta1a,theta1a,-theta1a];%铺层方式 %任意方向上的应力应变关系(关于x(y)的函数） for k=1:8 %层数 Q11h1(k)=Q11*cos(thetaa(k))^4+2*(Q12+2*Q66)*sin(thetaa(k))^2*cos(thetaa(k))^2+Q22*sin(thetaa(k))^4; Q12h1(k)=(Q11+Q22-4*Q66)*sin(thetaa(k))^2*cos(thetaa(k))^2+Q12*(sin(thetaa(k))^4+cos(thetaa(k))^4); Q16h1(k)=(Q11-Q12-2*Q66)*sin(thetaa(k))*cos(thetaa(k))^3+(Q12-Q22+2*Q66)*sin(thetaa(k))^3*cos(thetaa(k)); Q22h1(k)=Q11*sin(thetaa(k))^4+2*(Q12+2*Q66)*sin(thetaa(k))^2*cos(thetaa(k))^2+Q22*cos(thetaa(k))^4; Q26h1(k)=(Q11-Q12-2*Q66)*sin(thetaa(k))^3*cos(thetaa(k))+(Q12-Q22+2*Q66)*sin(thetaa(k))*cos(thetaa(k))^3; Q66h1(k)=(Q11+Q22-2*Q12-2*Q66)*sin(thetaa(k))^2*cos(thetaa(k))^2+Q66*(sin(thetaa(k))^4+cos(thetaa(k))^4); zz(k)=z(k+1)^3-z(k)^3; end %与x,y有关的系数的符号表达式 D11a1=1/3*Q11h1*zz'; D12a1=1/3*Q12h1*zz'; D16a1=1/3*Q16h1*zz'; D22a1=1/3*Q22h1*zz'; D26a1=1/3*Q26h1*zz'; D66a1=1/3*Q66h1*zz'; %求节点处的系数 for i=1:n for j=1:n x=x0(i); y=y0(j); D11(i,j)=eval(D11a1);%%%将字符串以表达式形式输出%%%%%%% D12(i,j)=eval(D12a1); D16(i,j)=eval(D16a1); D22(i,j)=eval(D22a1); D26(i,j)=eval(D26a1); D66(i,j)=eval(D66a1); end end %系数求解 for i=2:n-1 for j=2:n-1 a001(i-1,j-1)=D11(i,j); a002(i-1,j-1)=4*D16(i,j)*(1/x); a003(i-1,j-1)=2*(1/x^2)*(D12(i,j)+2*D66(i,j)); a004(i-1,j-1)=4*(1/x^2)*D26(i,j); a005(i-1,j-1)=(1/x^4)*D22(i,j); a006(i-1,j-1)=2*(1/x)*D11(i,j); a007(i-1,j-1)=-2*(1/x^3)*(D12(i,j)+2*D66(i,j)); a008(i-1,j-1)=-4*(1/x^4)*D26(i,j); a009(i-1,j-1)=-(1/x^2)*D22(i,j); a0010(i-1,j-1)=4*(1/x^3)*(D16(i,j)+D26(i,j)); a0011(i-1,j-1)=2*(1/x^4)*(D12(i,j)+D22(i,j)+2*D66(i,j)); a0012(i-1,j-1)=(1/x^3)*D22(i,j); a0013(i-1,j-1)=-4*(1/x^4)*(D16(i,j)+D26(i,j)); end end for i=1:n-2 for j=1:n-2 a01((j-1)*(n-2)+i)=a001(i,j);%将矩阵转换为数组%%%% a02((j-1)*(n-2)+i)=a002(i,j); a03((j-1)*(n-2)+i)=a003(i,j); a04((j-1)*(n-2)+i)=a004(i,j); a05((j-1)*(n-2)+i)=a005(i,j); a06((j-1)*(n-2)+i)=a006(i,j); a07((j-1)*(n-2)+i)=a007(i,j); a08((j-1)*(n-2)+i)=a008(i,j); a09((j-1)*(n-2)+i)=a009(i,j); a010((j-1)*(n-2)+i)=a0010(i,j); a011((j-1)*(n-2)+i)=a0011(i,j); a012((j-1)*(n-2)+i)=a0012(i,j); a013((j-1)*(n-2)+i)=a0013(i,j); end end a1=diag(a01);%%%%构建对角矩阵%%%%% a2=diag(a02); a3=diag(a03); a4=diag(a04); a5=diag(a05); a6=diag(a06); a7=diag(a07); a8=diag(a08); a9=diag(a09); a10=diag(a010); a11=diag(a011); a12=diag(a012); a13=diag(a013); m0=8*((3*theta*r^2)/8)*t*p/((3*theta*r^2)/8);%单位面积的质量 %振动方程和边界方程 ka1=a1*Ab4n+a6*Ab3n+a9*Ab2n+a12*Ab1n; kb1=Ab1n; %方程组化成矩阵形式，求kdd,kdb,kbb,kbd %%%%%%%%%%%止%%%%% %内点ka2 for k=1:n-4 ka2((k-1)*(n-4)+1:(k-1)*(n-4)+n-4,:)=ka1(k*(n-2)+2:k*(n-2)+n-3,:); end %内点矩阵kdd for k=1:n-4 kdd(:,(k-1)*(n-4)+1:(k-1)*(n-4)+n-4)=ka2(:,k*(n-2)+2:k*(n-2)+n-3); end %内点边界矩阵kdb kdb(:,1:n-2)=ka2(:,1:n-2); kdb(:,(n-2)+(n-4)*2+1:(n-2)+(n-4)*2+n-2)=ka2(:,(n-3)*(n-2)+1:(n-3)*(n-2)+n-2); for k=1:n-4 kdb(:,n-2+(k-1)*2+1)=ka2(:,k*(n-2)+1); kdb(:,n-2+(k-1)*2+2)=ka2(:,k*(n-2)+n-2); end %边界kb for k=1:n-4 kb(n-2+(k-1)*2+1,:)=kb1(k*(n-2)+1,:); kb(n-2+(k-1)*2+2,:)=kb1(k*(n-2)+n-2,:); end %边界矩阵kbb kbb(:,1:n-2)=kb(:,1:n-2); kbb(:,n-2+(n-4)*2+1:(n-2)+(n-4)*2+n-2)=kb(:,(n-3)*(n-2)+1:(n-3)*(n-2)+n-2); for k=1:n-4 kbb(:,n-2+(k-1)*2+1)=kb(:,k*(n-2)+1); kbb(:,n-2+(k-1)*2+2)=kb(:,k*(n-2)+n-2); end %边界内点矩阵kbd for k=1:n-4 kbd(:,(k-1)*(n-4)+1:(k-1)*(n-4)+n-4)=kb(:,k*(n-2)+2:k*(n-2)+n-3); end %求一阶频率freq k=kdd-kdb*inv(kbb)*kbd; m=m0*eye((n-4)*(n-4)); fre0=eig(inv(k)*m); fre1=max(fre0); freq=sqrt(1/fre1); %振型 [eigvector,eigvalue]=eig(inv(k)*m); wd=eigvector(:,1)/sum(eigvector(:,1)); wb=-inv(kbb)*kbd*wd; %化成标准形式w22,w32,w42,...,w23,w33,w43,... w1(1:n-2,:)=wb(1:n-2,:); w1((n-3)*(n-2)+1:(n-3)*(n-2)+n-2,:)=wb((n-2)+(n-4)*2+1:(n-2)+(n-4)*2+n-2,:); for k=1:n-4 w1(k*(n-2)+1,:)=wb(n-2+(k-1)*2+1,:); w1(k*(n-2)+n-2,:)=wb(n-2+(k-1)*2+2,:); end for k=1:n-4 w1(k*(n-2)+2:k*(n-2)+n-3,:)=wd((k-1)*(n-4)+1:(k-1)*(n-4)+n-4,:); end %画模态 for i=1:n-2 for j=1:n-2 w(i,j)=w1((j-1)*(n-2)+i); end end x1=x0(2:n-1); y1=y0(2:n-1); [X,Y]=meshgrid(x1,y1); f=mesh(X,Y,w); set(f,'FaceColor','white','EdgeColor','black'); colormap gray