理想流体力学课程设计(Hess Smith方法求附加质量)

一、物理背景无论是船舶还是海洋平台在海洋开发中都起着关键的作用，而开发海洋首先需要对海洋结构物进行深入地研究。

这其中，水动力学中的附加质量是研究的重要方面，掌握物体附加质量的计算无疑具有重要的意义。

附加惯性力的存在使物体在理想流体中的变速运动相当于物体自身质量上增加了一个附加质量而在真空中运动，换句话说，理想流体增大了物体的惯性，使物体很难加速也难减速。

计算机是求解附加质量的重要工具，本课程设计主要依据分布源模型的面元法等知识来对圆球、椭球、圆柱、双椭球的附加质量进行数值模拟计算，并进行相关讨论。

二、理论依据用s 表示无界流中的物体表面，来流为均匀流，其未扰动速度或无穷远处的速度为,1x y z V V i V j V k V V ∞∞∞∞∞∞=++== (2.1.1)用()Φx,y,z 表示定常速度势，它在物体外部空间域中适合拉普拉斯方程，在物面上适合不可进入条件，在无穷远处，应该与均匀来流的速度势吻合，即20∇Φ=（物体外） (2.1.2)0=∂∂nφ（物面s 上） (2.1.3) x y z xV yV zV ϕ∞∞∞Φ=+++（无穷远处）其中,单位法线向量n 指向物体内部。

在速度势Φ中分出已知的均匀来流项，记x y z xV yV zV ϕ∞∞∞Φ=+++ (2.1.4)这里的ϕ是扰动速度势，ϕ应适合以下定解条件：20ϕ∇=（物体外） (2.1.5)V n nϕ∞∂=-⋅∂（物面s 上） (2.1.6) 0ϕ→（无穷远处） (2.1.7)易知过物面s 的通量为零，即⎰⎰=∂∂sds n 0ϕ所以远方条件(2.1.7)可进一步具体化为⎪⎭⎫⎝⎛=21r O ϕ（r =→∞） (2.1.8) 用pq r 表示点p 和q 之间的距离，对函数()q ϕ和1/pq r 在物面s 外部和远方控制面c 的内部之空间域内用格林公式，当点p 在上述空间域内时()()()1114q pq qqpqs c p q q ds r n n r ϕϕϕπ+⎧⎫⎛⎫∂∂⎪⎪=-⎡⎤ ⎪⎨⎬⎣⎦ ⎪∂∂⎪⎪⎝⎭⎩⎭⎰⎰ (2.1.9) 从ϕ的远方条件(2.1.8)可知，c 上积分趋于零，式(2.1.9)成为()()()1114q pq qqpqs p q q ds r n n r ϕϕϕπ⎧⎫⎛⎫∂∂⎪⎪=-⎡⎤ ⎪⎨⎬⎣⎦ ⎪∂∂⎪⎪⎝⎭⎩⎭⎰⎰ (2.1.10) 其中, p 是物面s 外的任意一点。

在物体的内部域中构造一个合适的内部解i ϕ，它在s 内部适合拉普拉斯方程，在物面s 上适合某种物面条件，其具体形式将在下面给出。

对于上述物体外部的点p 函数1/pq r 在物体内部域中没有奇点，在内部域中对函数()i q ϕ和1/pq r 用格林公式，得到()()110i iq pq q q pqs q q ds r n n r ϕϕ⎧⎫⎛⎫∂∂⎪⎪=-⎡⎤ ⎪⎨⎬⎣⎦ ⎪∂∂⎪⎪⎝⎭⎩⎭⎰⎰ (2.1.11) 式(2.1.10)和(2.1.11)中的p 是物体外部同一个点，把两式相减，得到()()q s pqq i qiq pqds r n n nr p ⎰⎰⎥⎥⎦⎤⎢⎢⎣⎡⎪⎪⎭⎫ ⎝⎛∂∂--⎪⎪⎭⎫⎝⎛∂∂-∂∂=114ϕϕϕϕπϕ (2.1.12) 在物面s 上取i ϕ适合下述两种物面条件，得到两种i ϕ的定解条件，一种是：20(i i s s ϕϕϕ⎫∇=⎪⎬=⎪⎭内部)(上） (2.1.13) 定解问题(2.1.13)是拉普拉斯方程的第一类边值问题，它的解是存在且唯一的。

取式(2.1.12)中的内部解i ϕ为式(2.1.13)所决定的函数，则式(2.1.12)成为()()q spqds r q p ⎰⎰=σϕ (2.1.14)其中()⎪⎭⎫⎝⎛∂∂-∂∂=n n q i ϕϕπσ41 (2.1.15) 式(2.1.14)表示扰动速度势ϕ可以用物面s 上的分布源表示，其中分布源密度σ是未知函数，将由扰动势的物面条件（2.1.6）来决定。

物体的附加质量ij m ，表示物体沿i 方向运动引起的j 方向的附加质量，公 式如下:(,1,2,6)j ji i j isbsbm n ds ds i j nρρ∂Φ=Φ=Φ=∂⎰⎰⎰⎰ (2.2.1)利用式(2.1.14)，再结合物面条件，得到()()12q s p ppq s p q d V n n r επσσ∞-⎛⎫∂+=-⋅ ⎪ ⎪∂⎝⎭⎰⎰ (2.3.1) 这就是分布源密度σ所适合的线性积分方程。

把积分方程(2.3.1)转换成线性代数方程组，即用离散量代替连续变量。

把物面s 分成N 小块，记1Nj j s s ==∑∆ (2.3.2)用平面四边形或三角形来近似代替小曲面j s ∆。

具体做法如下，取第j 小块的四个顶点坐标之算术平均值，得到中心点j p 的坐标。

计算对角线向量的向量积（指向与曲面法线指向相符合），用j n 表示该方向上的单位向量，形成以j n 为法线且通过中心点j p 的平面，再把四个顶点向该平面作投影，以四个投影点为顶点组成平面四边形j Q ∆，用j Q ∆代替原来的小曲面j s ∆，称j Q ∆为单元。

通常把小范围内的分布源密度σ作为常数，因此只要分割不太粗，可以认为σ在单元j Q ∆上为常数，记作j σ，从而()q pq Q p j q pq p s ds r n ds r n q j j⎪⎪⎭⎫⎝⎛∂∂≈⎪⎪⎭⎫ ⎝⎛∂∂⎰⎰⎰⎰∆∆11σσ (2.3.3) 因此物面s 上的积分可以用N 个平面四边形（三角形）上积分之和来近似，即()q pq Q p j N j q pq p sds r n ds r n q j ⎪⎪⎭⎫⎝⎛∂∂∑≈⎪⎪⎭⎫ ⎝⎛∂∂⎰⎰⎰⎰∆=111σσ (2.3.4) 上式左端的未知量()q σ是连续型变量，而上式右端的未知量是N 个离散量(1)j j N σ=。

为了求解这N 个未知数，须要N 个方程。

取积分方程(2.3.1)中的动点p 为N 个单元j Q ∆的中心点(1)j p j N =，称之为控制点，即控制物面条件使之成立的点。

用近似式(2.3.4)代替积分方程(2.3.1)的左端，便可以写出j σ的N 阶线性代数方程组：1(1,2,,)Nij j i j a b i N σ=∑== (2.3.5)其中1()ii jij q p p q Q a ds j i n r ∆⎛⎫∂=≠ ⎪ ⎪∂⎝⎭⎰⎰2,i ij i p a b V n π∞==-⋅称ij a 为影响系数，即第j 个单元上的分布源在第i 个控制点上的影响。

求解线性代数方程组(2.3.5)得到j σ的值以后，便可以得到速度势()p ϕ在控制点i p 处的值，即()j ij Nj i c p σϕ1=∑≈ (2.3.6)1i jij q p qQ c ds r ∆=⎰⎰(2.3.7) 另外，物面s 上的诱导速度为()()()12Ni i j j i i j p n v p p s ϕπσσ='∇≈+∑∈∆ (2.3.8)其中'∑表示求和是不计j i =这一项。

，这里的曲面法线n 指向物体内部。

三、数值模型将物体表面划分成四边形面面元，物面为S ，每一个四边形面面元为j S ∆。

为了简化计算，将面网格投影到各自对应的平面上，使曲面网格j S ∆变为平面网格j Q ∆。

投影的方法为：取j S ∆四个顶点坐标之平均值，作为中心点j p 的坐标。

计算j S ∆对角线连线向量的向量积并使得积的方向与流域法向相同。

用j n 表示该方向上的单位向量。

设：(,,)j jx jy jz n n n n =j p 的坐标为：(,,)jx jy jz p p p则取投影面为过j p 并以j n 为法向量的平面：()()()0jx jx jy jy jz jz n x p n y p n z p -+-+-=设在该平面上的投影点为：(,,)x y z 而曲面四边形某个顶点为：000(,,)x y z 则有：000(,,)//(,,)jx jy jz x x y y z z n n n ---因此得到由j S ∆顶点坐标求解投影点（j Q ∆顶点）坐标的线性方程组：000000jx jx jy jy jz jz jx jy jz jy jx jy jx jzjx jz jx n p n p n p n n n x n n y n x n y n n z n x n z ⎡⎤++⎛⎫⎡⎤⎢⎥ ⎪⎢⎥-=-⎢⎥ ⎪⎢⎥⎢⎥ ⎪-⎢⎥-⎣⎦⎝⎭⎣⎦由此线性方程组可解出投影点坐标。

假设速度势和分布源在j Q ∆上是不变的，其值为该单元中点（控制点）处的速度势或分布源。

由于所求速度势和速度等物理量均为物面上的物理量，因此要令p 点落在物面之上。

式右端分布源的法向导数极限由两部分组成，一部分是P 点附近小曲面ε的贡献，另一部分是屋面其余部分s ε-贡献。

当p 所趋近于的物面上的点'p 作为控制点的单元，积分时需要考虑奇异性；其余部分为s ε-。

设其中一单元为单元i ，其余模型为单元j 。

对于每一个控制点i ，令1~j N =循环一次求得前述方程的积分项（包括奇异积分）。

再由1~i N =可以得到N 组方程，进而形成求解各个控制点处物理量的矩阵。

p 点表示控制点（编号i ），对于每一个控制点的物理量，通过在ε和s ε-上积分得到。

即：对于任意一i ＝１,２．．．Ｎ，有()q pq Q p j N j q pq p sds r n ds r n q j ⎪⎪⎭⎫⎝⎛∂∂∑≈⎪⎪⎭⎫ ⎝⎛∂∂⎰⎰⎰⎰∆=111σσ 1(1,2,,)N ij j i j a b i N σ=∑==其中， 1()iij ij q p p q Q a ds j i n r ∆⎛⎫∂=≠ ⎪ ⎪∂⎝⎭⎰⎰ 2,i ii i p a b V n π∞==-⋅ 将模型控制点数据导入到程序中计ij a 和i b，可以得到方程组：111111N N NN N N b a a a a b σσ⎡⎤⎡⎤⎛⎫ ⎪⎢⎥⎢⎥= ⎪⎢⎥⎢⎥ ⎪⎢⎥⎢⎥⎝⎭⎣⎦⎣⎦求解上面线性代数方程组得到j σ的值以后，便可以得到速度势()p ϕ在控制点i p 处的值，即()j ij Nj i c p σϕ1=∑≈ 其中，1i jij q p qQ c ds r ∆=⎰⎰每一个i 点（控制点）处诱导速度为（是向量）：11()()()Ni j j i i i j p v p p s σ=∇Φ≈∈∆∑其中，1(d )()()1lim (d )()j i j pi q pq Q j i pi q p p pq Q s i j r v p s i j r ∆→∆⎧∇≠⎪⎪=⎨⎪∇=⎪⎩⎰⎰⎰⎰ 在常分布单元假设条件下：1lim (d )2ij pi q i p p pqQ s n r π→∆∇=⎰⎰可得，11()2'()()Ni j i j j i i i j p n v p p s πσσ=∇Φ≈+∈∆∑于是便可求解出速度势和物面速度。