第三章-势流理论教学内容

第六章：势流理论一．内容总结：二元流动包括平面流动和轴对称流动。

对于不可压缩流体的平面定常势流可以引入流函数和速度势函数。

而不可压缩平面势流速度势函数和流函数均满足拉普拉斯方程。

速度势函数的等值线与流函数等值线正交，流函数的等值线与流线重合。

本章研究物体在静止理想流体中平面运动时,流体对物体的作用力。

求解势流问题的思路为：当物体在流体中运动，即物体与流体之间产生相对运动时,物体受到流体的作用力。

对于理想流体的运动不存在切应力，理想流体中运动的物体表面上只受到法向的压力作用。

因此要解决在流场中物体所受的作用力,只要把物体表面上合压力求出即可。

由伯努利方程可知，若物面上（理想流体中无分离绕流时物面与流线重合）的速度分布已知可求出物面上压力分布，再沿物面积分便可求出物体受到的合压力。

因此，问题归结为求出流场的速度分布，对于不可压缩平面流动，求速度分布的问题又可归结为求速度势函数和流函数问题。

1． 势流问题求解的思路 基本方程 ： 20ϕ∇= 无旋流动20ψ∇=二维不可压缩流动V grad φ=G即得到三个速度分量u v 伯努立方程压力,,w →→P 再由边界条件→ 积分 spds ∫便求得了合力，因此只要确定V ϕ→→p G就可积分求合力了。

对于二维不可压缩无旋流动，整个问题的关键在于找到满足边界条件的ϕ或ψ。

求速度势ϕ的方法：因为方程是线性方程, 几个解的线性之和仍满足拉普拉斯方程。

20ϕ∇=根据已知知识确定应选的势流. 简单平面势流的表示式 1) 等速直线运动等速V 平行x 轴的平行流动速度势和流函数为： 0V x ϕ= 0V y ψ=2) 源和汇源心在坐标原点时速度势和流函数在平面极坐标下为： ln 2Q r ϕπ= 2Q ψθπ= 式中为源 为汇0Q >0Q <3) 旋涡速度势和流函数在平面极坐标下为： 2ϕθπΓ= ln 2r ψπΓ=−4)偶极子速度势和流函数为：222M x z x y ϕπ=+ 222M yx yψπ=−+ 221214sin p p p c V θρ∞∞−==− 在位置上,指向与X 轴成β角. 0z M ：称偶极矩,由汇指向源。

《流体力学》教学大纲一、课程基本信息二、课程概述中文：本课程是工程力学专业的学类核心课程，以高等数学、理论力学、材料力学为前导课程，着重培养学生分析解决实际工程中流体力学问题的能力。

本课程主要包括流体的平衡、流体力学的基本方程、不可压缩无粘流动、涡旋运动、平面势流等，强调应用这些基本概念及定律分析与流体力学相关的工程问题，学生需了解流体力学的发展现状和趋势，理解流体力学中的基本概念、基本理论及基本定律，掌握流体力学的实验、分析与数值计算的基本技能与基本方法，并能灵活运用这些基本概念及定律分析与流体力学相关的工程问题。

通过学习本课程，让学生学会流体力学基本理论，获得解决流体工程问题的基本技能，锻炼和提升对复杂的流体工程问题进行简化，从而建立数学模型并进行求解的能力。

英文：This is a bas ic course for majors of engineering mechanics, aiming at students’ physical concepts and basic principles commonly used to analyze engineering problems related to fluid mechanics, thus laying a solid foundation for their research and design in aerospace, mechanical, civil, chemical, environmental and ocean. Theapplications of the dimensional and order analysis method in engineering are emphasized in this course. The study of this course develops the students’ ability to simplify the complex problems, prese nt and solve the mathematic model of related engineering problems. The main contents of this course are the basic equations of fluid mechanics, incompressible in-viscid flow, the motion of vortex, dimensional analysis, incompressible viscid flow. Prerequisites: Advanced Mathematics, Mathematics Physics Equation, Field Theory，Theoretical Mechanics，Mechanics of Materials.三、课程内容（一）课程教学目标设置本课程是为了让工程力学专业的学生对工程力学专业知识体系的重要组成板块之一的流体力学进行较为系统的学习，并深度掌握与理解，具备应用流体力学的基本知识和基本理论分析解决生产实际工程问题的能力。

高等流体力学第一章 流体力学的基本概念连续介质：流体是由一个紧挨着一个的连续的质点所组成的，没有任何空隙的连续体，即所 谓的连续介质。

流体质点：是指微小体积内所有流体分子的总和。

欧拉法质点加速度：时变加速度与位变加速度和zuu y u u x u u t u dt du a x z x y x x x x x ∂∂+∂∂+∂∂+∂∂==质点的随体导数：质点携带的物理量随时间的变化率称为质点的随体导数，用dtd表示。

在欧拉法描述中的任意物理量Q 的质点随体导数表述如下：x kk Qu t Q dt dQ ∂∂+∂∂= 式中Q 可以是标量、矢量、张量。

质点的随体导数公式对任意物理量都成立，故将质点的随体导数的运算符号表示如下：x kk u t dt d ∂∂+∂∂= 其中t∂∂称为局部随体导数，x k k u ∂∂称为对流随体导数，即在欧拉法描述的流动中，物理量的质点随体导数等于局部随体导数与对流随体导数之和。

体积分的随体导数：质点携带的物理量随时间的变化率称为质点的随体导数。

则在由流体质点组成的流动体积V 中标量函数Φ（x, t ）随时间的变化率就是体积分的随导函数。

由两部分组成①函数Φ 对时间的偏导数沿体积V 的积分，是由标量场的非恒定性引起的。

②函数Φ通过表面S 的通量。

由体积V 的改变引起的。

()dV divv dt d dV v div t dS u dV t dV dt d v v n s v v ⎥⎦⎤⎢⎣⎡Φ+Φ=⎥⎦⎤⎢⎣⎡Φ+∂Φ∂=Φ+∂Φ∂=Φ⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰()dV adivv dt da dV av div t a dS au dV t a adV dt d v v n s v v ⎥⎦⎤⎢⎣⎡+=⎥⎦⎤⎢⎣⎡+∂∂=+∂∂=⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰⎰ 变形率张量： 11ε12ε13εD ij = 21ε 22ε 23ε 31ε 32ε 33ε其中ii ε表示所在方向的线性变形率，其余ij ε（j i ≠）为角变形率。

第六章势流理论本章内容：1．势流问题求解的思路2．库塔----儒可夫斯基条件3. 势流的迭加法绕圆柱的无环绕流，绕圆柱的有环绕流4．布拉休斯公式5．库塔----儒可夫斯基定理学习这部分内容的目的有二：其一，获得解决势流问题的入门知识，即关键问题是求解速度势。

求出速度势之后，可按一定的步骤解出速度分布、压力分布，以及流体和固体之间的作用力。

其二，明确两点重要结论：１）园柱体在理想流体中作等速直线运动时，阻力为零（达朗贝尔疑题）；升力也为零。

２）园柱本身转动同时作等速直线运动时，则受到升力作用（麦格鲁斯效应）。

本章重点：1、平面势流问题求解的基本思想。

2、势流迭加法3、物面条件，无穷远处条件4、绕圆柱有环流，无环流流动的结论，即速度分布，压力分布，压力系数分布，驻点位置，流线图谱，升力，阻力，环流方向等。

5、四个简单势流的速度势函数，流函数及其流线图谱。

6、麦马格鲁斯效应的概念7、计算任意形状柱体受流体作用力的卜拉修斯定理8、附加惯性力，附加质量的概念本章难点：1．绕圆柱有环流，无环流流动的结论，即速度分布，压力分布，压力系数分布，驻点位置，流线图谱，升力，阻力，环流方向等。

2．任意形状柱体受流体作用力的卜拉修斯定理3．附加惯性力，附加质量的概念§６-1 几种简单的平面势流平面流动：平面上任何一点的速度、加速度都平行于所在平面，无垂直于该平面的分量；与该平面相平行的所有其它平面上的流动情况完全一样。

例如：1）绕一个无穷长机翼的流动，2）船舶在水面上的垂直振荡问题，由于船长比宽度及吃水大得多，且船型纵向变化比较缓慢，可以近似认为流体只在垂直于船长方向的平面内流动，如图６-２所示。

如果我们在船长方向将船分割成许多薄片，并且假定绕各薄片的流动互不影响的话，则这一问题就可以按平面问题处理。

这一近似方法在船舶流体力学领域内称为切片理论。

一、均匀流流体质点沿ｘ轴平行的均匀速度Vo ，如图６-5所示，V ｘ＝V o ， V y ＝０平面流动速度势的全微分为dx V dy V dx V dy ydx x d y x 0=+=∂∂+∂∂=ϕϕϕ 积分：φ＝V ox （６-4） 如图６-3流函数的全微分为，dy V dy V dx V dy ydx x d o x y =+-=∂∂+∂∂=ψψψ 积分：ψ＝V o ｙ （６-5） 如图６-4由（６-4）和（６-5）可得： 流线：ｙ=const ，一组平行于ｘ轴的直线，如图６-3中的实线。

势流理论笔记：01势流理论基础前⾔：势流理论复习笔记，没想到⾃⼰⼜重新学了⼀遍势流理论。

所以记个笔记笔记内容基本摘抄⾃朱仁传⽼师的《船舶在波浪上的运动理论》，写得好哇基础理论均匀、不可压缩理想流体的流场中，连续性⽅程与欧拉⽅程可以描述为：∇v=0∂∂t+v⋅∇v=−∇pρ+gzv(x,y,z)与p(x,y,z)分别为速度⽮量与压⼒场，存在向量关系∇v22=∇v⋅v2=(v⋅∇)v+v×(∇×v)欧拉⽅程可以改写为以下形式，称为兰姆⽅程(Lamb′sEquation)∂v∂t+∇v22−v×(∇×v)=−∇pρ+gz以上共四个⽅程，四个未知数，⽅程封闭。

如果流体流动⽆旋有势，⽅程可以进⼀步简化v=∇ϕ(x,y,z,t)∇2ϕ(x,y,z,t)=∂2ϕ∂x2+∂2ϕ∂y2+∂2ϕ∂z2=0pρ+gz+v22+∂ϕ∂t=C(t)格林函数法船舶在波浪中的运动问题关键在于求解流畅中的速度势，即求在确定边界条件下的拉普拉斯⽅程。

格林函数法(Green′s function method)是⼀类成熟常⽤的求解⽅法。

格林函数法的基础势格林公式（散度定理）推导得到，对三维空间中有界区域τ，有以下关系式∭τ∇⋅A dτ=∬S n⋅A d S其中，S为空间域τ充分光滑的边界⾯；n为曲⾯S的单位外法向⽮量（从流体域内指向外部），⽮量A在封闭区域τ+S上连续。

现令A=ϕ∇ψ，于是有∇A=∇(ϕ∇ψ)=∇ϕ⋅∇ψ+ϕ∇2ψn⋅A=n⋅(ϕ∇ψ)=ϕ∂ψ∂n将A=ϕ∇ψ代⼊格林公式得到∬Sϕ∂ψ∂n d S=∭τ∇ϕ⋅∇ψdτ+∭τϕ⋅∇2ψdτ将A=ψ∇ϕ代⼊格林公式得到∬Sψ∂ϕ∂n d S=∭τ∇ϕ⋅∇ψdτ+∭τψ⋅∇2ϕdτ两式作差得到{()() ()()()()∬S ψ∂ϕ∂n −ϕ∂ψ∂n d S =∭τϕ⋅∇2ψ−ψ⋅∇2ϕd τ若ϕ,ψ在τ内处处调和，即: ∇2ϕ=0,∇2ψ=0，则有∬S ψ∂ϕ∂n −ϕ∂ψ∂n d S =0∬S ψ∂ϕ∂n d S =∬S ϕ∂ψ∂n d S称作格林第⼆公式。

第三章 流体运动学本章在连续介质假设下，讨论描述流体运动的方法，根据运动要素的特性对流动进行分类。

本章的讨论是纯运动学意义上的，不涉及流动的动力学因素。

连续方程是质量守恒定律对流体运动的一个具体约束，也在本章的讨论范围之中。

§3—1 描述流动的方法一. 拉格朗日法和欧拉法● 拉格朗日法是质点系法，它定义流体质点的位移矢量为：r r a b c t =(,,,)，其中(,,)(,,,)a b c r a b c t =0是拉格朗日变数，即t 0时刻质点的空间位置，用来对连续介质中无穷多个质点进行编号，作为质点标签。

● 欧拉法是流场法，它定义流体质点的速度矢量场为：u u x y z t =(,,,)，其中(,,)x y z 是空间点（场点）。

流体的其它物理特性和运动要素也都用对应于时间与空间域的场的形式描述。

二. 流体质点的加速度、质点导数● 在拉格朗日观点下，流体质点加速度的求法是比较简单的。

求速度和加速度只须将位移矢量直接对时间求一、二阶导数即可，求导时a,b,c 作为参数不变，意即跟定流体质点。

u r t rt a u t u t r t =====d d ,d d ∂∂∂∂∂∂22.● 欧拉法中流体质点加速度的表达必须特别注意，求加速度需要跟定流体质点，于是 x,y,z均随 t 变，而且),,(d ),,d(z y x u u u tz y x =，所以加速度 u u tz u u y u u x u u t u t z z u t y y u t x x u t u t u a z y x)(d d d d d d d d ∇⋅+=+++=+++==∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂∂. ● 建立 t 时刻和 t+dt 时刻的流场图，假设一流体质点在 t 时刻位于场点 M ，t + dt 时刻它到达场点M ’，在 t+dt 时刻的流场图上再标上与点M 处于同一位置的场点M 1，此时有另一个流体质点占据该场点。

势流理论笔记：03Hess-Smith积分⽅法书接上回Hess-Smith⽅法采⽤⾯向对象编程的思路，Matlab程序脚本，实现以下功能：输⼊⾯元（四边形⾯元顶点坐标）输出系数矩阵[H][M]以及[V x],[V y]Hess-Smith积分⽅法思路待积分的⽬标函数：φ=∬A d A rφ(P)=∬A1r(P,Q)d A(Q)1. 在Q⾯元建⽴局部坐标系，将待积点P转换⾄该局部坐标系。

2. 在局部坐标系下，利⽤圆柱坐标系推导积分公式（四个三⾓积分之和=⾯元积分）3. 可以拓展得到该⽬标函数对于局部坐标系三个⽅向导数的积分（三个坐标轴的投影）具体步骤建⽴局部坐标系按照顺序输⼊⼤地坐标系O−XYZ下⾯元四个顶点坐标p1，p2，p3，p4，以z轴⾓度看，四个点顺序呈顺时针排列，以四边形⾯元中⼼O为坐标原点，x⽅向平⾏于→p1p3，建⽴⾯源局部坐标系将四个点投影⾄局部坐标系ξO′η平⾯下的得到(ξ1,η1),(ξ2,η2),(ξ3,η3),(ξ4,η4),则r可以表⽰为r=√(x−ξ)2+(y−η)2+z2⽤于计算下式的积分（选择简单的1r作为积分的格林函数）φ=∬Ad Ar坐标系之间的转换矩阵坐标原点:O′=p1+p2+p3+p44x′→x′=p3−p1 |p3−p1|z′→z′=(p4−p2)×(p3−p1) |(p4−p2)×(p3−p1)|y′→y′=→z′×→x′坐标变换:X′=(X−O′)⋅AA=[x′;y′;z′]′最后的’表⽰矩阵转置建⽴局部圆柱坐标系以(X,Y,0)为坐标原点建⽴圆柱坐标系，则r可以表⽰为r=√R2+z2原积分式可以写作\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 矩阵tichess_smith=Hess_Smith(Panel); %创建Hess_Smith对象H=hess_smith.Vz; %系数矩阵HM=hess_smith.M; %系数矩阵MN=hess_smith.Num; %⾯元数⽬toc%% 圆球绕流算例C=0.5*eye(N,N); %系数矩阵CU=[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:NQ1=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); %⼊射势endA=4*pi*C+H;b=M*Phi_n;Phi_D=A\b; %绕射势Phi=Phi_D+Phi_in; %总速度势=绕射势+⼊射势%% 绘图figure(1111)plot(Phi,'r');hold onplot(Phi_t,'b');legend("Phi Cal","Phi Theory");Hess\_Smith.m\quad classclassdef Hess_Smith < handle%HESS_SMITH 求解积分系数矩阵% 输⼊：四边形⾯元坐标% 输出：系数矩阵[H]即Vz、[M]propertiesPanel; %⾯元坐标点Num; %⾯元数⽬M; %输出系数矩阵Vx;Vy,Vz; %其他⽅向导数矩阵%% ⾯元局部坐标系参数Center; Tran;A;center;endmethodsfunction 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();endfunction [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; %坐标变换endfunction [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、JP=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;endfunction [Phi,Vx,Vy,Vz]=Cal_Phi(self,i,j)%公式(4.4.18)&(4.4.19)% 输⼊:P⾯元索引i，Q⾯元索引jP=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;endPhi=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);endfunction Cal_H_M(self)%计算 H & V_z Vx Vy矩阵for j=1:1:self.Numfor 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);endself.Tran{j,1}=self.A;self.Center(j,:)=self.center;endfprintf("H-M矩阵建⽴完成!\n");end%%%post plotendend Loading [MathJax]/extensions/TeX/mathchoice.js。