圆柱绕流 附加质量 matlab计算程序

N=input（‘输入N：‘）; 
M=input（‘输入M：‘）;  
K=input（‘输入K：‘）;  
R=1;    %单元划分数（2N+M）*K圆柱底面半径为1；
%划分结构化网格
L=10*R;
rr=linspace（0RN+1）;
ll=linspace（0LM+1）;
thet=linspace（02*piK+1）;
x=zeros（2*N+M+1K+1）;
y=zeros（2*N+M+1K+1）;
z=zeros（2*N+M+1K+1）;
P=cell（2*N+MK）;
A=cell（2*N+M+1K+1）;
delta_S=zeros（2*N+MK）;
delta_r=R/N;
delta_l=L/M;
delta_thet=2*pi/K;

for i=1:N;
    x（i:）=0;
    for j=1:K+1;
        z（ij）=rr（i）*sin（thet（j））;
        y（ij）=rr（i）*cos（thet（j））;
        A{ij}=[x（ij）y（ij）z（ij）];
    end
end
for i=N+M+1:2*N+M+1;
    x（i:）=L;
    for j=1:K+1;
        z（ij）=rr（2*N+M+2-i）*sin（thet（j））;
        y（ij）=rr（2*N+M+2-i）*cos（thet（j））;
        A{ij}=[x（ij）y（ij）z（ij）];
    end
end
for i=N+1:N+M;
    x（i:）=ll（i-N）;
    for j=1:K+1;
        z（ij）=R*sin（thet（j））;
        y（ij）=R*cos（thet（j））;
        A{ij}=[x（ij）y（ij）z（ij）];
    end
end

for j=1:K;
     P{1j}=1/3*（A{1j}+A{2j}+A{2j+1}）;
     P{2*N+Mj}=1/3*（A{2*N+Mj}+A{2*N+Mj+1}+A{2*N+M+1j}）;
end

for i=1:N;
    delta_S（i:）=sin（delta_thet）*（rr（i+1）^2-rr（i）^2）/2;
    delta_S（2*N+M-i+1:）=delta_S（i:）;
end
sw=sin（delta_thet/2）*R*2*delta_l;
for i=N+1:N+M;
    delta_S（i:）=sw;
end
for i=2:2*N+M-1;
    for j=1:K;
        P{ij}=0.25*（A{ij}+A{ij+1}+A{i+1j}+A{i+1j+1}）;
    end
end 
surf（xyzz）
xlabel（‘x‘）
ylabel（‘y‘）
zlabel（‘z‘）
axis image

PP=2*N+M;
%开始进行单元计算，解方程组
%设远方速度为1，x正方向
a=zeros（PP*KPP*K）;
b=zeros（PP*K1）;
c=zeros（PP*KPP*K）;
sigma=zeros（PP*K1）;
PHI=zeros（PP*K1）;
r=zeros（PP*KPP*K）;

%求圆柱m_11
n_p=cell（PP*K1）;
PHI_1=zeros（PP*K1）;

for k=