势流理论笔记:03 Hess-Smith积分方法

书接上回势流理论笔记:02 直接法与间接法

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} \]

  1. \(Q\)面元建立局部坐标系,将待积点\(P\)转换至该局部坐标系。
  2. 在局部坐标系下,利用圆柱坐标系推导积分公式(四个三角积分之和=面元积分)
  3. 可以拓展得到该目标函数对于局部坐标系三个方向导数的积分(三个坐标轴的投影)

具体步骤

建立局部坐标系

按照顺序输入大地坐标系\(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
posted @ 2021-09-27 13:07  陈橙橙  阅读(1167)  评论(0编辑  收藏  举报