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
