%%%cell_cen為cell中心的坐標(biāo)，N*3的矩陣

%%%cell_s為每個cell的面積，為N維列向量

%%%%f為每個cell對應(yīng)的法向量，為N*3的矩陣

%%%svd_solve_zjj為另外定義的SVD奇異值分解法函數(shù)，解病態(tài)方程組，避免奇異工作精度

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個數(shù)

N = length(f);

u = [10,0,0];%u描述三維來流速度

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

density = 1*10^3;%牛頓流體介質(zhì)的密度

Rad = 0.45;

%%創(chuàng)建點源坐標(biāo)矩陣

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));%%%點源位置嚴(yán)重影響準(zhǔn)確性

end

%% 右側(cè)列向量

rhs = zeros(N,1);

for i=1:N%先固定一個i，即考察所有點源對圓上一個節(jié)點的影響

?for j=1:N

??sub = cell_cen(i,:,:)-dy(j,:,:);%第j個點源到第i個節(jié)點的相對坐標(biāo)向量

??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('點源強度求解結(jié)束!');

disp('開始計算輸出信息!');

s = zeros(N,1);%%

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

??for j = 1:N

????sub = cell_cen(i,:)-dy(j,:);%第j個點源到第i個節(jié)點的相對坐標(biāo)向量

??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(['數(shù)值計算附加質(zhì)量為',num2str(M)])

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

disp(['解析附加質(zhì)量為',num2str(M0)])