Hess-Smith方法
采用面向对象编程的思路，\(Matlab\)程序脚本，实现以下功能：
- 输入面元（四边形面元顶点坐标）
- 输出系数矩阵\([H][M]\)以及\([V_x],[V_y]\)
Hess-Smith积分方法思路
待积分的目标函数：
\[\begin{align} \varphi&=\iint_{A}\frac{\text d A}{r}\\ \varphi(P)&=\iint_A\frac{1}{r(P,Q)} \text dA(Q) \end{align} \]- 在\(Q\)面元建立局部坐标系，将待积点\(P\)转换至该局部坐标系。
- 在局部坐标系下，利用圆柱坐标系推导积分公式（四个三角积分之和=面元积分）
- 可以拓展得到该目标函数对于局部坐标系三个方向导数的积分（三个坐标轴的投影）
具体步骤
建立局部坐标系
按照顺序输入大地坐标系\(O-XYZ\)下面元四个顶点坐标\(p_1\)，\(p_2\)，\(p_3\)，\(p_4\)，以\(z\)轴角度看，四个点顺序呈顺时针排列，以四边形面元中心\(O\)为坐标原点，\(x\)方向平行于\(\overrightarrow{p_1p_3}\)，建立面源局部坐标系将四个点投影至局部坐标系\(\xi O'\eta\)平面下的得到\((\xi_1,\eta_1)\),\((\xi_2,\eta_2)\),\((\xi_3,\eta_3)\),\((\xi_4,\eta_4)\),则\(r\)可以表示为
\[r=\sqrt{(x-\xi)^2+(y-\eta)^2+z^2} \]用于计算下式的积分（选择简单的\(\frac{1}{r}\)作为积分的格林函数）
\[\varphi=\iint_{A}\frac{\text d A}{r} \]坐标系之间的转换矩阵
\[\begin{equation} \begin{aligned} \text{坐标原点:}\qquad O'&=\frac{p_1+p_2+p_3+p_4}{4}\\ x'\qquad \overrightarrow x'&=\frac{p_3-p_1}{|p_3-p_1|}\\ z'\qquad \overrightarrow z'&=\frac{(p_4-p_2)\times(p_3-p_1)}{|(p_4-p_2)\times(p_3-p_1)|}\\ y'\qquad \overrightarrow y'&=\overrightarrow z'\times\overrightarrow x'\\ \text{坐标变换:}\qquad X'&=(X-O')\cdot \boldsymbol {A}\\ \boldsymbol {A}&=[x';y';z']'\quad \text{最后的’表示矩阵转置} \end{aligned} \end{equation} \]建立局部圆柱坐标系
以\((X,Y,0)\)为坐标原点建立圆柱坐标系，则\(r\)可以表示为
\[r=\sqrt{R^2+z^2} \]原积分式可以写作
\[\begin{align} \varphi&=\oint\int_0^R\frac{R\text dR\text d\theta}{\sqrt{R^2+z^2}}\\ dr&=\frac{R\text dR}{\sqrt{R^2+z^2}}\\ \varphi&=\oint\left[ \int^{r}_{|z|} \text dr\right]\text d\theta\\ &=\oint(r-|z|)\text d \theta \end{align} \]Hess-Smith积分公式
四个三角形面积积分之和便是四边形的数值积分
\[\begin{align} \varphi=\iint_{A}\frac{\text d A}{r}=\oint r\text d\theta-|z|\Delta\theta \end{align} \]局部坐标系下的几何参数定义：
\[\begin{equation} \begin{aligned} d_{12}&=\sqrt{(\xi_2-\xi_1)^2+(\eta_2-\eta_1)^2}\\ C_{12}&=\frac{\xi_2-\xi_1}{d_{12}}\\ S_{12}&=\frac{\eta_2-\eta_1}{d_{12}} \end{aligned} \end{equation} \] \[\begin{equation} \begin{aligned} s_{12}^{(1)}&=(\xi_1-x,\eta_1-y)\cdot(C_{12},S_{12})\\ s_{12}^{(2)}&=(\xi_2-x,\eta_2-y)\cdot(C_{12},S_{12})\\ R_{12}&=(x-\xi_1,y-\eta_1)\cdot(S_{12},-C_{12}) \end{aligned} \end{equation} \] \[\begin{equation} \begin{aligned} r_1=\sqrt{(x-\xi_1)^2+(y-\eta_1)^2+z^2}\\ r_2=\sqrt{(x-\xi_2)^2+(y-\eta_2)^2+z^2}\\ \end{aligned} \end{equation} \]积分公式可以写作
\[\begin{align} \varphi_{12}=\oint r\text d\theta-|z|\Delta\theta=R_{12}Q_{12}+|z|J_{12} \end{align} \]式中：
第一项
\[\begin{equation} \begin{aligned} R_{12}\cdot Q_{12}=\oint r\text d\theta&=\oint \frac{R_{12}}{cos(\theta)}\text d\theta\\ &=\left .R_{12}\cdot\left( \text {ln}{\frac{1+sin\theta}{cos\theta}}\right)\right|^{\theta_{2}}_{\theta_{1}}\\ &=R_{12}\cdot\left(\text{ln}\frac{(1+sin\theta_2)\cdot cos(\theta_1)} {(1+sin\theta_1)\cdot cos(\theta_2)}\right)\\ &=R_{12}\cdot\left(\text{ln}\frac{r_2+s_{12}^{(2)}}{r_1+s_{12}^{(1)}} \right)\\ &=R_{12}\cdot\text{ln}\frac{r_1+r_2+d_{12}}{r_1+r_2-d_{12}} \end{aligned} \end{equation} \]式中：
\[\begin{align} sin\theta_1=\frac{s_{12}^{(1)}}{r_1},cos\theta_1=\frac{R_{12}}{r_1}\\ sin\theta_2=\frac{s_{12}^{(2)}}{r_2},cos\theta_2=\frac{R_{12}}{r_2} \end{align} \]第二项
\[\begin{align} |z|J_{12}=-|z|\Delta\theta&=\text{sgn}(R_{12}) \left[\text{tan}^{-1}\left( \left |\frac{z}{R_{12}}\right |\frac{s_{12}^{(2)}}{r_2}\right)- \text{tan}^{-1}\left( \left |\frac{z}{R_{12}}\right |\frac{s_{12}^{(1)}}{r_1}\right)\right]\\ &=\text{tan}^{-1}\left[\frac{R_{12}|z|(r_1s_{12}^{(2)}-r_2s_{12}^{(1)})} {r_1r_2R_{12}^{2}+z^2s_{12}^{(2)}s_{12}^{(1)}} \right] \end{align} \]综上，得到
\[\begin{equation} \begin{aligned} \varphi&=\iint_{A}\frac{\text d A}{r}\\ &=\varphi_{12}+\varphi_{12}+\varphi_{12}+\varphi_{12}-|z|\Delta\theta \end{aligned} \end{equation} \]\(\Delta\theta\)与局部坐标系下的\(P\)点与面元四顶点的相对位置相关（可以由四个三角积分之和与四边形面元积分之间的集合关系得出）
\[\begin{align} \Delta \theta=\begin{cases} 2\pi,& P\in Q\\ 0,& p\notin Q \end{cases} \quad or \quad \Delta \theta=\begin{cases} 2\pi,& R_{12}, R_{23}, R_{34}, R_{41}>0\\ 0,& else \end{cases} \end{align} \] \[\begin{equation} \begin{aligned} \frac{\partial \varphi}{\partial n}&=\iint_{A}\frac{\partial}{\partial n}\left(\frac{1}{r}\right)\text d A\\ &=V_z=-\frac{\partial \varphi}{\partial n}\\ &=\text{sgn}(z)[\Delta\theta-J_{12}-J_{23}-J_{34}-J_{41}] \end{aligned} \end{equation} \] \[\begin{align} V_x&=-\frac{\partial \varphi}{\partial x}=-S_{12}Q_{12}-S_{23}Q_{23}-S_{34}Q_{34}-S_{41}Q_{41}\\ V_y&=-\frac{\partial \varphi}{\partial y}=C_{12}Q_{12}+C_{23}Q_{23}+C_{34}Q_{34}+C_{41}Q_{41}\\ \end{align} \]参考文献
[1] Hess J L , Smith A . Calculation of potential flow about arbitrary bodies[J]. Progress in Aerospace Sciences, 1967, 8(none):1-138.
示例
\(demo.m\)
%% 面元
clc;clear all;
%导入面元
load('Panel.mat')
%% H-M 矩阵
tic
hess_smith=Hess_Smith(Panel); %创建Hess_Smith对象
H=hess_smith.Vz; %系数矩阵H
M=hess_smith.M; %系数矩阵M
N=hess_smith.Num; %面元数目
toc
%% 圆球绕流算例
C=0.5*eye(N,N); %系数矩阵C
U=[1,0,0]; %来流速度
a=1; %圆球半径
%速度势理论值
Phi_theory=@(x)U(1)*x(1)+0.5*U(1)*x(1)*(a^3)/(x(1)^2+x(2)^2+x(3)^2)^(1.5);
for i=1:1:N
Q1=Panel(4*i-3,:);Q2=Panel(4*i-2,:);Q3=Panel(4*i-1,:);Q4=Panel(4*i,:);
center=0.25*(Q1+Q2+Q3+Q4);
Phi_t(i,1)=Phi_theory(center); %速度势理论值
z_new=cross(Q4-Q2,Q3-Q1);
z_new=z_new/norm(z_new);
Phi_n(i,1)=-U*z_new'; %边界条件
Phi_in(i,1)=center(1)*U(1); %入射势
end
A=4*pi*C+H;
b=M*Phi_n;
Phi_D=A\b; %绕射势
Phi=Phi_D+Phi_in; %总速度势=绕射势+入射势
%% 绘图
figure(1111)
plot(Phi,'r');hold on
plot(Phi_t,'b');legend("Phi Cal","Phi Theory");
\(Hess\_Smith.m\quad class\)
classdef Hess_Smith < handle
%HESS_SMITH 求解积分系数矩阵
% 输入：四边形面元坐标
% 输出：系数矩阵[H]即Vz、[M]
properties
Panel; %面元坐标点
Num; %面元数目
M; %输出系数矩阵
Vx;Vy,Vz; %其他方向导数矩阵
%% 面元局部坐标系参数
Center; Tran;
A;center;
end
methods
function self = Hess_Smith(Panel)
%HESS_SMITH 构造函数
% 输入坐标点，四个点一组，尺寸:[4n,3]
self.Panel = Panel;
[self.Num,~]=size(self.Panel);
self.Num=floor(self.Num/4);
self.M=zeros(self.Num,self.Num);
self.Vz=zeros(self.Num,self.Num);
self.Vx=zeros(self.Num,self.Num);
self.Vy=zeros(self.Num,self.Num);
%面元局部坐标系参数
self.Center=zeros(self.Num,3);
self.Tran=cell(self.Num,1);
%计算系数矩阵
self.Cal_H_M();
end
function [P_c,Q_loc] = G2LOC(self,P,Q)
%Global to local 坐标转换
% 输入面元P，Q：以Q建立局部坐标系
% 输出：P_c(面元中心),Q_loc，局部坐标系下的坐标
Q_center=0.25*(Q(1,:)+Q(2,:)+Q(3,:)+Q(4,:));
self.center=Q_center;
P_center=0.25*(P(1,:)+P(2,:)+P(3,:)+P(4,:));
z_new=cross((Q(4,:)-Q(2,:)),Q(3,:)-Q(1,:));
z_new=z_new/norm(z_new);
x_new=(Q(3,:)-Q(1,:));
x_new=x_new/norm(x_new);
y_new=cross(z_new,x_new);
self.A=[x_new;y_new;z_new]'; %坐标转换矩阵
Q_loc=(Q-Q_center)*self.A; %坐标变换
P_c=(P_center-Q_center)*self.A; %坐标变换
end
function [Phi12,R12,Q12,J12,C12,S12]=Cal_Phi_Line(self,P_in,Q1_in,Q2_in)
%公式（4.4.17）
% 输入：P，Q为局部坐标系下的坐标，P_in=(1,3),Q1_in=(1,3),Q2_in=(1,3)
% 输出：三角区域积分值,中间参数R、Q、J
P=P_in(1:2);z=P_in(3);
Q1=Q1_in(1:2);Q2=Q2_in(1:2);
d12=norm(Q2-Q1);
direc=(Q2-Q1)/d12; %线段方向矢量
C12=direc(1);S12=direc(2);
s12_1=(Q1-P)*direc';
s12_2=(Q2-P)*direc';
R12=(P-Q1)*[S12,-C12]';
r1=sqrt(z^2+(P-Q1)*(P-Q1)');
r2=sqrt(z^2+(P-Q2)*(P-Q2)');
Q12=log((r1+r2+d12)/(r1+r2-d12));
J12=atan2(R12*abs(z)*(r1*s12_2-r2*s12_1),...
r1*r2*R12^2+z^2*s12_2*s12_1);
Phi12=R12*Q12+abs(z)*J12;
end
function [Phi,Vx,Vy,Vz]=Cal_Phi(self,i,j)
%公式(4.4.18)&(4.4.19)
% 输入:P面元索引i，Q面元索引j
P=self.Panel(4*i-3:4*i,:);
Q=self.Panel(4*j-3:4*j,:);
[P_c,Q_loc]=self.G2LOC(P,Q); %转换至局部坐标系
z=P_c(3);
%line 1-2
[Phi12,R12,Q12,J12,C12,S12]=self.Cal_Phi_Line(P_c,Q_loc(1,:),Q_loc(2,:));
%line 2-3
[Phi23,R23,Q23,J23,C23,S23]=self.Cal_Phi_Line(P_c,Q_loc(2,:),Q_loc(3,:));
%line 3-4
[Phi34,R34,Q34,J34,C34,S34]=self.Cal_Phi_Line(P_c,Q_loc(3,:),Q_loc(4,:));
%line 4-1
[Phi41,R41,Q41,J41,C41,S41]=self.Cal_Phi_Line(P_c,Q_loc(4,:),Q_loc(1,:));
theta=2*pi;
if (R12<0||R23<0||R34<0||R41<0)
theta=0;
end
Phi=Phi12+Phi23+Phi34+Phi41-abs(z)*theta;
Vx=-S12*Q12-S23*Q23-S34*Q34-S41*Q41;
Vy=C12*Q12+C23*Q23+C34*Q34+C41*Q41;
Vz=sign(z)*(theta-J12-J23-J34-J41);
end
function Cal_H_M(self)
%计算 H & V_z Vx Vy矩阵
for j=1:1:self.Num
for i=1:1:self.Num
[self.M(i,j),self.Vx(i,j),self.Vy(i,j),self.Vz(i,j)]=...
self.Cal_Phi(i,j);
end
self.Tran{j,1}=self.A;
self.Center(j,:)=self.center;
end
fprintf("H-M矩阵建立完成!\n");
end
%%
%post plot
end
end