势流理论分析长短轴不同的椭球附加质量系数

最新推荐文章于 2023-04-01 12:05:40 发布

hser-chen

最新推荐文章于 2023-04-01 12:05:40 发布

阅读量4.9k

点赞数 14

文章标签：水动力分析势流理论面元法

本文链接：https://blog.csdn.net/weixin_41788456/article/details/84452654

版权

一、物理背景

当今世界，无论是船舶还是海洋平台在海洋开发中都起着关键的作用，而开发海洋首先需要对海洋结构物进行深入地研究。这其中，水动力学中的附加质量是研究的重要方面，掌握物体附加质量的计算无疑具有重要的意义。

附加惯性力的存在使物体在理想流体中的变速运动相当于在物体自身质量上增加了一个附加质量而在真空中运动，换句话说，理想流体增大了物体的惯性，使物体很难加速也难减速。附加质量有物体上的单位运动绝对速度势来决定，与流体密度，物体形状和运动方向有关，而与时间无关。

计算机是求解附加质量的重要工具，本次大作业主要是为了研究无界流中运动物体所受到的流体动力，在假设的不可压缩理想流体的无旋运动环境下，利用分布源模型的面元法并代入物面边界条件和远方条件，求流体速度势在控制体内的分布，进而可以得到所求物体的附加质量。

二、理论依据

根据格林公式、分布奇点方法、面元法、分布源模型理论求解定解问题。介绍如下：

2.1赫斯-史密斯(Hess-Smith)方法

用S表示无界流中的物体表面，来流为均匀流，其未扰动速度或无穷远处的速度为

         $\overline{V}_{\infty }=V_{\infty x}\overline{i} +V_{\infty y}\overline{j}+V_{\infty z}\overline{k}$

                  $V_{\infty }=\left |\overline{V}_{\infty }\right |=1$                                                (2.1.1)

用 $\phi (x,y,z)$ 表示定常速度势，它在物体外部空间域中适合拉普拉斯方程，在物面上适合不可进入条件，在无穷远处，应该与均匀来流的速度势吻合。即

$\bigtriangledown ^{2}\phi =0$     （物体外）                   (2.1.2)

$\frac{\partial \phi }{\partial n}=0$ （物面s上）                   (2.1.3)

其中,单位法线向量 $\overline{n}$ 指向物体内部。

在速度势 $\phi$ 中分出已知的均匀来流项，记

$\phi =xV_{\infty x} +yV_{\infty y}+zV_{\infty z}+\varphi$                 (2.1.4)

这里的 $\varphi$ 是扰动速度势， $\varphi$ 应适合以下定解条件：

$\bigtriangledown ^{2}\varphi =0$ （物体外）                      (2.1.5)

$\frac{\partial \varphi \phi }{\partial n}=0$ （物面s上）                    (2.1.6)

$\varphi \rightarrow 0$ （无穷远处）                      (2.1.7)

易知过物面s的通量为零，即

$\iint_{s}\frac{\partial n }{\partial x}ds=0$

所以远方条件(2.1.7)可进一步具体化为

$\varphi =0(\frac{1}{r^{2 }})$ （ $r=\sqrt{x^{2}+y^{2}+z^{2}}\rightarrow \infty$ ）           (2.1.8)

用 $r_{pq}$ 表示点p和q之间的距离，对函数 $\varphi (q)$ 和在1/ $r_{pq}$ 物面s外部和远方控制面c的内部之空间域内用格林公式，当点p在上述空间域内时

$\varphi (p)=\frac{1}{4\pi }\iint_{s+c}\left \{ \frac{1}{r_{pq}}\frac{\partial [\varphi (q)]}{\partial n}-\varphi (q)\frac{\partial }{\partial n_{q}}(\frac{1}{r_{pq}}) \right \}ds_{q}$                                (2.1.9)    从 $\varphi$ 的远方条件(2.1.8)可知，c上积分趋于零，式(2.1.9)成为

        $\varphi (p)=\frac{1}{4\pi }\iint_{s}\left \{ \frac{1}{r_{pq}}\frac{\partial [\varphi (q)]}{\partial n}-\varphi (q)\frac{\partial }{\partial n_{q}}(\frac{1}{r_{pq}}) \right \}ds_{q}$                            (2.1.10)

其中,p是物面s外的任意一点。

在物体的内部域中构造一个合适的内部解 $\varphi _{i}$ ，它在s内部适合拉普拉斯方程，在物面s上适合某种物面条件，其具体形式将在下面给出。对于上述物体外部的点p函数1/ $r_{pq}$ 在物体内部域中没有奇点，在内部域中对函数 $\varphi _{i}(q)$ 和1/ $r_{pq}$ 用格林公式，得到

        $0=\frac{1}{4\pi }\iint_{s}\left \{ \frac{1}{r_{pq}}\frac{\partial [\varphi_{i} (q)]}{\partial n}-\varphi _{i}(q)\frac{\partial }{\partial n_{q}}(\frac{1}{r_{pq}}) \right \}ds_{q}$ (2.1.11)

式(2.1.10)和(2.1.11)中的是物体外部同一个点，把两式相减，得到

     $4\pi \varphi (p)=\iint_{s}\left \{ \frac{1}{r_{pq}}\left ( \frac{\partial\varphi }{\partial n_{q}} -\frac{\partial \varphi _{i}}{\partial n_{q}}\right )-\left ( \varphi -\varphi _{i} \right )\frac{\partial }{\partial n_{q}}(\frac{1}{r_{pq}}) \right \}ds_{q}$ (2.1.12)

在物面s上取 $\varphi _{i}$ 适合下述两种物面条件，得到两种 $\varphi _{i}$ 的定解条件，一种是：

                 $\bigtriangledown ^{2}\varphi _{i}=0$ （s内部）

                    $\varphi _{i}=\varphi$ （s上）    (2.1.13)

定解问题(2.1.13)是拉普拉斯方程的第一类边值问题，它的解是存在且唯一的。取式(2.1.12)中的内部解为式(2.1.13)所决定的函数，则式(2.1.12)成为

                $\varphi (p)=\iint_{s}\frac{\delta (q)}{r_{pq}}ds_{q}$        (2.1.14)

其中

　       $\sigma (q)=\frac{1}{4\pi }\left ( \frac{\partial \varphi }{\partial n}-\frac{\partial \varphi _{i}}{\partial n} \right )$       　(2.1.15)

式(2.1.14)表示扰动速度势可以用物面上的分布源表示，其中分布源密度是未知函数，将由扰动势的物面条件（2.1.6）来决定。

2.2附加质量

物体的附加质量 $m_{ji}$ ，表示物体沿i方向运动引起的j方向的附加质量，公式如下:

       $m_{ji}=\rho \iint_{sb}\phi _{i}n_{j}ds=\rho \iint_{sb}\varphi _{i}\frac{\partial\phi _{j}}{\partial n}ds$    (2.2.1)

2.3计算模型

利用式(2.1.14)，再结合物面条件，得到

            $2\pi \delta (p)+\iint_{s-\varepsilon }\sigma (q)\frac{\partial }{\partial n_{q}}\left ( \frac{1}{r_{pq}} \right )ds_{q}=-\vec{V_{\infty }}\vec{n}_{p}$    $\left ( i,j=1,2,...,6 \right )$ (2.3.1)

这就是分布源密度 $\sigma$ 所适合的线性积分方程。

把积分方程(2.3.1)转换成线性代数方程组，即用离散量代替连续变量。把物面s分成N小块，记

                           $s=\sum_{j=1}^{N}\bigtriangleup s_{j}$ (2.3.2)

用平面四边形或三角形来近似代替小曲面。具体做法如下 $\bigtriangleup s_{j}$ ，取第j小块的四个顶点坐标之算术平均值，得到中心点 $p_{j}$ 的坐标。计算对角线向量的向量积（指向与曲面法线指向相符合），用 $\vec{n}_{j}$ 表示该方向上的单位向量，形成以 $\vec{n}_{j}$ 为法线且通过中心点的平面，再把四个顶点向该平面作投影，以四个投影点为顶点组成平面四边形 $\bigtriangleup Q_{j}$ ，用 $\bigtriangleup Q_{j}$ 代替原来的小曲面 $\bigtriangleup s_{j}$ ，称 $\bigtriangleup Q_{j}$ 为单元。通常把小范围内的分布源密度 $\sigma$ 作为常数，因此只要分割不太粗，可以认为 $\sigma$ 在单元 $\bigtriangleup Q_{j}$ 上为常数，记作，从而

        $\iint_{\bigtriangleup s_{j}}\sigma (q)\frac{\partial }{\partial n_{q}}\left ( \frac{1}{r_{pq}} \right )ds_{q}\approx \sigma (j)\iint_{\bigtriangleup Q_{j}}\frac{\partial }{\partial n_{q}}\left ( \frac{1}{r_{pq}} \right )ds_{q}$       (2.3.3)

因此物面s上的积分可以用N个平面四边形（三角形）上积分之和来近似，即

          $\iint_s\sigma (q)\frac{\partial }{\partial n_{q}}\left ( \frac{1}{r_{pq}} \right )ds_{q}\approx\sum_{j=1}^{N} \sigma (j)\iint_{\bigtriangleup Q_{j}}\frac{\partial }{\partial n_{q}}\left ( \frac{1}{r_{pq}} \right )ds_{q}$ (2.3.4)

上式左端的未知量 $\sigma (q)$ 是连续型变量，而上式右端的未知量是N个离散量 $\sigma _{j}(j=1\sim N)$ 。为了求解这N个未知数，须要N个方程。取积分方程(2.3.1)中的动点p为N个单元 $\bigtriangleup Q_{j}$ 的中心点 $p_{j}(j=1\sim N)$ ，称之为控制点，即控制物面条件使之成立的点。用近似式(2.3.4)代替积分方程(2.3.1)的左端，便可以写出 $\sigma _{j}$ 的N阶线性代数方程组：

              $\sum_{j=1}^{N}a_{ij}\sigma _{j}=b_{i}(i=1\sim N)$      (2.3.5)

其中

$a_{ij}=\iint_{\bigtriangleup Q_{j}}\frac{\partial }{\partial n_{q}}\left ( \frac{1}{r_{pq}} \right )ds_{q}$ $\left ( j\neq i \right )$

$a_{ij}=2\pi ,b_{i}=-\vec{V}_{\infty }\cdot \vec{n}_{p_{i}}$

称 $a_{ij}$ 为影响系数，即第j个单元上的分布源在第i个控制点上的影响。

求解线性代数方程组(2.3.5)得到 $\sigma _{i}$ 的值以后，便可以得到速度势 $\phi (p)$ 在控制点 $p_{i}$ 处的值，即

                   $\varphi (p_{i})\approx \sum_{j=1}^{N}c_{ij}\sigma _{j}$         (2.3.6)

                 $c_{ij}=\iint_{\bigtriangleup Q_{j}}\frac{1}{r_{p_{i}q}}ds_{q}$         (2.3.7)

另外，物面s上的诱导速度为

          $\bigtriangledown \varphi (p)\approx 2\pi \sigma _{i}n_{i}+\sum_{j=1}^{N}^{'}\sigma _{j}v_{j}(p_{i})$      $\left ( p\in \bigtriangleup s_{i} \right )$ (2.3.8)

其中 $\sum_{j=1}^{N}^{'}$ 表示求和是不计j=i这一项。，这里的曲面法线 $\vec{n}$ 指向物体内部。

三、数值模型

将物体表面划分成四边形面面元，物面为S，每一个四边形面面元为 $\bigtriangleup S_{j}$ 。为了简化计算，将面网格投影到各自对应的平面上，使曲面网格 $\bigtriangleup S_{j}$ 变为平面网格 $\bigtriangleup Q_{j}$ 。投影的方法为：

取 $\bigtriangleup S_{j}$ 四个顶点坐标之平均值，作为中心点 $p_{j}$ 的坐标。计算 $\bigtriangleup S_{j}$ 对角线连线向量的向量积并使得积的方向与流域法向相同。用 $\vec{n}_{j}$ 表示该方向上的单位向量。

设：

$\vec{n}_{j}=({n}_{jx},{n}_{jy},{n}_{jz})$

$p_{j}$ 的坐标为：

$({p}_{jx},{p}_{jy},{p}_{jz})$

则取投影面为过 $p_{j}$ 并以 $\vec{n}_{j}$ 为法向量的平面：

                            $n_{jx}(x-p_{jx})+n_{jy}(y-p_{jy})+n_{jz}(z-p_{jz})=0$

设在该平面上的投影点为：

而曲面四边形某个顶点为： $(x_{0},y_{0},z_{0})$

则有： $\left ( x-x_{0} ,y-y_{0}, z-z_{0} \right )//\left ( n_{jx},n_{jy} ,n_{jz}\right )$

因此得到由顶点 $\bigtriangleup S_{j}$ 坐标求解投影点（顶点 $\bigtriangleup Q_{j}$ ）坐标的线性方程组：

                      $\begin{pmatrix} n_{jx} &n_{jy} & n_{jz} \\ n_{jy} &- n_{jx} & 0\\ n_{jz} & 0 &- n_{jx} \end{pmatrix}\begin{bmatrix} x\\ y \\ z \end{bmatrix}=\begin{bmatrix} n_{jx} p_{jx} +n_{jy} p_{jy} +n_{jz} p_{jz} \\n_{jy} x_{0}- n_{jx} y_{0} \\ n_{jz} x_{0}- n_{jx} z_{0}\end{bmatrix}$

     由此线性方程组可解出投影点坐标。

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

由于所求速度势和速度等物理量均为物面上的物理量，因此要令p点落在物面之上。式右端分布源的法向导数极限由两部分组成，一部分是P点附近小曲面的贡献，另一部分是其余部分 $s-\varepsilon$ 贡献。当p所趋近于的物面上的点 $p^{'}$ 作为控制点的单元，积分时需要考虑奇异性；其余部分为 $s-\varepsilon$ 。设其中一单元为单元i，其余模型为单元j。对于每一个控制点i，令j=1~N循环一次求得前述方程的积分项（包括奇异积分）。再由I=1~N可以得到N组方程，进而形成求解各个控制点处物理量的矩阵。

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

即：对于任意i＝１,２．．．Ｎ，有

$\iint_s\sigma (q)\frac{\partial }{\partial n_{q}}\left ( \frac{1}{r_{pq}} \right )ds_{q}\approx\sum_{j=1}^{N} \sigma (j)\iint_{\bigtriangleup Q_{j}}\frac{\partial }{\partial n_{q}}\left ( \frac{1}{r_{pq}} \right )ds_{q}$

$\sum_{j=1}^{N}a_{ij}\sigma _{j}=b_{i}(i=1\sim N)$

其中，         $a_{ij}=\iint_{\bigtriangleup Q_{j}}\frac{\partial }{\partial n_{q}}\left ( \frac{1}{r_{pq}} \right )ds_{q}$    $\left ( j\neq i \right )$            $a_{ij}=2\pi ,b_{i}=-\vec{V}_{\infty }\cdot \vec{n}_{p_{i}}$

将模型控制点数据导入到程序中计 $a_{ij}$ 和 $b_{i}$ ，可以得到方程组：

$\begin{bmatrix} a_{11} &\cdots & a_{1N}\\ \vdots &\ddots &\vdots \\ a_{N1}&\cdots & a_{NN} \end{bmatrix}\begin{bmatrix} \sigma _{1}\\ \vdots \\ \sigma _{N} \end{bmatrix}=\begin{bmatrix} b _{1}\\ \vdots \\ b _{N} \end{bmatrix}$

求解上面线性代数方程组得到 $\sigma _{j}$ 的值以后，便可以得到速度势 $\varphi _{p}$ 在控制点 $p_{i}$ 处的值，即

          $\varphi \left ( p_{i} \right )\approx \sum_{j=1}^{N}c_{ij}\sigma _{j}$    其中， $c_{ij}=\iint_{\bigtriangleup Q_{j}}\frac{1}{r_{p_{i}q}}ds_{q}$

每一个i点（控制点）处诱导速度为（是向量）：

$\bigtriangledown \phi _{1}(p_{i})\approx \sum_{j=1}^{N}\sigma _{j}\vec{v}_{j}(p_{j})$          $\left ( p_{i} \in \bigtriangleup s_{i}\right )$

其中， $\vec{v}_{j}(p\left ( i \right ))=\left\{\begin{matrix} \bigtriangledown _{p_{i}}(\iint_{\bigtriangleup Q_{I}}\frac{1}{r_{pq}}ds_{q})(i\neq j)\\\\ \lim_{p\rightarrow p_{i}}\bigtriangledown _{p_{i}}(\iint_{\bigtriangleup Q_{I}}\frac{1}{r_{pq}}ds_{q})(i=j) \end{matrix}\right.$

在常分布单元假设条件下：

$\lim_{p\rightarrow p_{i}}\bigtriangledown _{p_{i}}(\iint_{\bigtriangleup Q_{I}}\frac{1}{r_{pq}}ds_{q})=2\pi \vec{n}_{i}$

    可得，

$\bigtriangledown \phi _{1} (p_{i})\approx 2\pi \sigma _{i}\vec{n}_{i}+\sum_{j=1}^{N}^{'}\sigma _{j}\vec{v}_{j}(p_{i})$      $\left ( p_{i}\in \bigtriangleup s_{i} \right )$

于是便可求解出速度势和物面速度。以上给出的数值离散求解方法称为基于分布源模型的面元法，又称赫斯-史密斯方法。

四、几何模型

对于计算附加质量M11而言，来流方向是正对着长轴平面的X轴的，而对于M33而言，来流方向是正对着短轴平面的X轴。本次的程序是沿着两个方向划分网格的，沿着长轴a的方向将长轴等分成m份，在短轴b的方向，是将其轴向的2π角度平局分成n份，所以总的网格数是m×n。由于汉斯史密斯方法主要的研究对象是平面四边形，所以该程序在长轴的顶点处会切掉一小部分以避免三角形单元的出现，由于切掉的只是一小部分，所以对椭球体的计算结果并不会产生太大的影响。

利用计算机编程来完成几何模型的用origin建立并划分网格如图4.1、图4.2、图4.3、图4.4和图4.5所示。其中，椭球长短轴之比分别为1:1，2:1，3:1，4:1，5:1 。

图4.1 椭球长短轴之比为1:1

图4.2 椭球长短轴之比为1:2

图4.3 椭球长短轴之比为1:3

图4.4 椭球长短轴之比为1:4

图4.5 椭球长短轴之比为1:5

五、计算参数及结果讨论

对于旋转的椭球体而言，在长短轴相等的情况下，其相当于是圆球。在长短轴之比从1到5的变化过程中，分别对每种情况计算五种网格的划分形式。网格数依次为170，340，500，625，875。分别计算M11和M33的数值大小，得出多组数据以后分析结果，并在本报告的结尾给出网格划分程序的代码。

5.1计算结果

表5.1计算结果

m*n	面元个数	长轴	短轴	体积V	M11	M11理论值	M33	M33理论值
10*17	170	3	3	106.968	0.569422	0.5	0.569421	0.5
		6		213.936	0.242866	0.21	0.765914	0.704
		9		320.904	0.142085	0.122	0.874952	0.804
		12		427.872	9.52E-02	0.082	0.936038	0.86
		15		534.84	6.90E-02	0.059	0.973641	0.894
20*17	340	3	3	109.5079	0.535184	0.5	0.539671	0.5
		6		219.0159	0.226374	0.21	0.763893	0.704
		9		328.524	0.132113	0.122	0.87332	0.804
		12		438.0318	8.86E-02	0.082	0.934557	0.86
		15		547.5402	6.43E-02	0.059	0.97231	0.894
25*20	500	3	3	110.5527	0.528844	0.5	0.532691	0.5
		6		221.1053	0.22332	0.21	0.753829	0.704
		9		331.6579	0.130197	0.122	0.861764	0.804
		12		442.2106	8.73E-02	0.082	0.922165	0.86
		15		552.7632	6.34E-02	0.059	6.34E-02	0.894
25*25	625	3	3	111.2113	0.527909	0.5	0.526308	0.5
		6		199.7135	0.149147	0.21	0.743707	0.704
		9		333.6339	0.129974	0.122	0.849719	0.804
		12		444.8451	8.71E-02	0.082	0.909032	0.86
		15		556.0565	6.33E-02	0.059	0.945606	0.894
35*25	875	3	3	111.5317	0.521544	0.5	0.521544	0.5
		6		223.0631	0.219852	0.21	0.742732	0.704
		9		334.5953	0.128008	0.122	0.849041	0.804
		12		446.1268	8.58E-02	0.082	0.908541	0.86
		15		557.6587	6.23E-02	0.059	0.94524	0.894

利用程序总共计算了对应于长轴短轴比从5到1，网格数从170到875的50个数据，将这些数据描述与上面的的表格中。由表中数据可以看出，在各个长短轴比不同的情况下，附加质量M11和M33随网格数目变化都有着一样的趋势即下降趋势。因此，我在此处仅描绘出长短轴之比a：b=5的情况下的附加质量和网格点数的变化图，如下图3.4和3.5所示。

5.2数据处理与分析

5.2.1 数据处理

由于各个长短轴比的情况下附加质量随网格数目有着相同的变化趋势，因此图5.1、图5.2附加质量随网格数变化曲线中可以看出，随着网格数目增加，无论是M11还是M33，其都会越来越接近于理论值，因此在程序允许的最大网格数（N<900）的情况下，网格数越大，越能反映附加质量的真实值。

经查阅资料可知，椭球附加质量系数M11的理论值为0.059、M33的理论值为0.894。不同面元数时，M11、M33的计算结果见表5.2。

表5.2 M11、M33计算值与理论值误差统计

单元数	170	340	500	625	875
M11	6.90E-02	6.43E-02	6.34E-02	6.33E-02	6.23E-02
M1理论值	0.059	0.059	0.059	0.059	0.059
误差	17.02%	9.03%	7.41%	7.24%	5.51%
M33	0.973641	0.97231	0.94524	0.945606	0.94524
M3理论值	0.894	0.894	0.894	0.894	0.894
误差	8.91%	8.76%	5.73%	5.77%	5.73%

图 5.1 M11和M11理论值随网格数变化曲线

图5.2 M11和M11理论值随网格数变化曲线

图5.3 长短轴随面元数变化曲线

由图5.3可以看出，椭球越细长，其附加质量系数越小。随着长短轴比的增加，附加质量系数逐渐收敛。

5.3 结果分析

由计算结果可以看出，椭球在附加质量系数M11较M33小很多，这点由椭球几何形状可以容易看出，与实际有较好的符合。随着面元数量的增加附加质量系数M11和M33均更接近于各自理论值，误差都有减小。

在计算过程中发现，面元的划分形式存在质量好坏的区别，对于同样多数量的面元，面元形式的不同会造成不一样的结果，比如同是数量为875的单元，长短轴比为4时，长轴方向划分为35份，周向25份得出的M33为0.908541，而长轴方向划分为25份，周向为35份，得出的附加质量M33为0.8946489，更接近与理论值，误差更小，这是由于对于本细长的椭球模型而言，周方向划分25份已经较为紧密，而细长的长轴方向需要划分更多的数量以使面元更接近于真实物面，所以在划分面元时，需要考虑物体的实际情况，以得到更高质量的面元。

参考文献

戴遗山,段文洋. 船舶在波浪中运动的势流理论[M]. 国防工业出版社，2008
孙寒冰,邹劲,庄家园,吴恭兴.深V型无人艇附加质量计算及其影响[D].2011年中国智能自动化学术会议论文集（第一分册）.2011.
林超友.潜艇近海底航行附加质量数值计算[D].SHIP ENGINEERING,2003