m=input('y方向涡格数,板上涡格的数量');
n=input('x方向涡格数,板上涡格的数量');
km=input('y方向板长度');
kn=input('x方向板长度');
knn=input('x方向总长度');
nn=input('x方向涡格数,总涡格的数量');
alpha=input('尾迹修正系数');
x=sym(zeros(m,n));
x1=sym(zeros(m,n));
x2=sym(zeros(m,n));
y=sym(zeros(m,n));
y1=sym(zeros(m,n));
y2=sym(zeros(m,n));
a=sym(zeros(m,n,m,n));
b=sym(zeros(m,n,m,n));
c=sym(zeros(m,n,m,n));
d=sym(zeros(m,n,m,n));
e=sym(zeros(m,n,m,n));
t1=sym(zeros(m,n,m,n));
tt1=sym(zeros(m*n,m));
tt2=sym(zeros(m*n,(nn-n-3)*m));
for i=1:1:m; %y坐标
for j=1:1:n; %x坐标
for t=1:1:m; %y坐标
for s=1:1:nn; %x坐标
x(i,j)=(4*j-1)*kn/4/n;
y(i,j)=(2*i-1)*km/2/m;
x1(t,s)=(4*s-3)*knn/4/nn;
y1(t,s)=(t-1)*km/m;
x2(t,s)=(4*s-3)*knn/4/nn;
y2(t,s)=t*km/m;
a(t,s,i,j)=1/((x(i,j)-x1(t,s))*(y(i,j)-y2(t,s))-(x(i,j)-x2(t,s))*(y(i,j)-y1(t,s)));
b(t,s,i,j)=((x2(t,s)-x1(t,s))*(x(i,j)-x1(t,s))+(y2(t,s)-y1(t,s))*(y(i,j)-y1(t,s)))/sqrt(((x(i,j)-x1(t,s))^2)+((y(i,j)-y1(t,s))^2));
c(t,s,i,j)=((x2(t,s)-x1(t,s))*(x(i,j)-x2(t,s))+(y2(t,s)-y1(t,s))*(y(i,j)-y2(t,s)))/sqrt(((x(i,j)-x2(t,s))^2)+((y(i,j)-y2(t,s))^2));
d(t,s,i,j)=(y1(t,s)-y(i,j));
e(t,s,i,j)=(1+((x(i,j)-x2(t,s))/sqrt(((x(i,j)-x2(t,s))^2)+((y(i,j)-y2(t,s))^2))))/(y2(t,s)-y(i,j));
t1(t,s,i,j)=(a(t,s,i,j).*(b(t,s,i,j)-c(t,s,i,j))+d(t,s,i,j)-e(t,s,i,j))/4/pi;
end
end
end
end
t2=t1(:,1:n,:,:); '板的系数矩阵';
t3=reshape(t2,[],m*n); '该命令使得矩阵成为列向量,包含所有m*n的元素';
t4=t3.'; '矩阵拉直,此矩阵即为板上涡格系数矩阵';
t5=t1(:,n+1,:,:); '尾迹第一格的系数矩阵';
t6=reshape(t5,[],m*n);
t7=t6.';
t8=repmat(t7,1,n);
t9=t4-t8; 't9为板上t+1时刻矩阵,t8为板上t时刻矩阵';
t10=t1(:,nn,:,:);
t11=reshape(t10,[],m*n);
t12=t11.'; '尾迹最后一列系数矩阵,需要乘前一列t时刻涡强';
t13=alpha*t12; '乘最后一列t时刻涡强';
t14=t1(:,n+2:nn-1,:,:);
t15=reshape(t14,[],m*n);
t16=t15.'; '尾迹中间列系数矩阵,需要乘对应列t+1时刻涡强';
aa=[t9 tt1 t16 tt1]; 't+1涡强系数矩阵';
bb=[t8 tt1 tt2 t12 t13]; 't涡强系数矩阵';
评论3