MATLAB三维圆球均匀绕流附加质量计算（代码自取）

%%%cell_cen为cell中心的坐标，N*3的矩阵

%%%cell_s为每个cell的面积，为N维列向量

%%%%f为每个cell对应的法向量，为N*3的矩阵

%%%svd_solve_zjj为另外定义的SVD奇异值分解法函数，解病态方程组，避免奇异工作精度

clear all

close all

cell_s = xlsread("C:\Users\86182\Desktop\R0.45.csv",'A2:A345');

cell_cen = xlsread("C:\Users\86182\Desktop\R0.45.csv",'B2:D345');

f = xlsread("C:\Users\86182\Desktop\R0.45.csv",'E2:G345');

%%%%%% N为三维物面划分的cell个数

N = length(f);

u = [10,0,0];%u描述三维来流速度

U = sqrt(u(1)^2+u(2)^2+u(3)^2);

density = 1*10^3;%牛顿流体介质的密度

Rad = 0.45;

%%创建点源坐标矩阵

dy = zeros(N,3);

for i = 1:N

  dy(i,1:3) = cell_cen(i,1:3)-f(i,1:3)*sqrt(cell_s(i));%%%点源位置严重影响准确性

end

%% 右侧列向量

rhs = zeros(N,1);

for i=1:N%先固定一个i，即考察所有点源对圆上一个节点的影响

 for j=1:N

  sub = cell_cen(i,:,:)-dy(j,:,:);%第j个点源到第i个节点的相对坐标向量

  r = sqrt(sub(1)^2+sub(2)^2+sub(3)^2);

  Aij(i,j) = dot(f(i,:),sub/r^3/4/pi);%注意sub/r为一个单位向量。

  rhs(i) = -dot(u,f(i,:));%%rhs(i）向左移项，右边0即圆上Vr

 end

end

m = svd_solve_zjj(Aij,rhs);

disp('点源强度求解结束!');

disp('开始计算输出信息!');

s = zeros(N,1);%%

for i = 1:N%%%%%第i个cell

  for j = 1:N

    sub = cell_cen(i,:)-dy(j,:);%第j个点源到第i个节点的相对坐标向量

  r = sqrt(sub(1)^2+sub(2)^2+sub(3)^2);

 s(i) = s(i)+m(j)/4/pi/r*dot(u,f(i,:))*cell_s(i);

  end

end

p = 0;

for i = 1:N

  p = p+s(i);

end

M = -p/U^2*density;

disp(['数值计算附加质量为',num2str(M)])

M0 = density*2/3*pi*Rad^3;

disp(['解析附加质量为',num2str(M0)])