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

%%%cell_s為每個(gè)cell的面積，為N維列向量

%%%%f為每個(gè)cell對(duì)應(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個(gè)數(shù)

N = length(f);

u = [10,0,0];%u描述三維來(lái)流速度

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

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

Rad = 0.45;

%%創(chuàng)建點(diǎn)源坐標(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));%%%點(diǎn)源位置嚴(yán)重影響準(zhǔn)確性

end

%% 右側(cè)列向量

rhs = zeros(N,1);

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

?for j=1:N

??sub = cell_cen(i,:,:)-dy(j,:,:);%第j個(gè)點(diǎn)源到第i個(gè)節(jié)點(diǎn)的相對(duì)坐標(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為一個(gè)單位向量。

??rhs(i) = -dot(u,f(i,:));%%rhs(i）向左移項(xiàng)，右邊0即圓上Vr

?end

end

m = svd_solve_zjj(Aij,rhs);

disp('點(diǎn)源強(qiáng)度求解結(jié)束!');

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

s = zeros(N,1);%%

for i = 1:N%%%%%第i個(gè)cell

??for j = 1:N

????sub = cell_cen(i,:)-dy(j,:);%第j個(gè)點(diǎn)源到第i個(gè)節(jié)點(diǎn)的相對(duì)坐標(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ù)值計(jì)算附加質(zhì)量為',num2str(M)])

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

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