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

理想流体力学课程设计(Hess Smith方法求附加质量) 理想流体力学课程设计(Hess Smith方法求附加质量) 一、物理背景 无论是船舶还是海洋平台在海洋开发中都起着关键的作用，而开发海洋首先需要对海洋结构物进行深入地研究。这其中，水动力学中的附加质量是研究的重要方面，掌握物体附加质量的计算无疑具有重要的意义。 附加惯性力的存在使物体在理想流体中的变速运动相当于物体自身质量上增加了一个附加质量而在真空中运动，换句话说，理想流体增大了物体的惯性，使物体很难加速也难减速。 计算机是求解附加质量的重要工具，本课程设计主要依据分布源模型的面元法等知识来对圆球、椭球、圆柱、双椭球的附加质量进行数值模拟计算，并进行相关讨论。 二、理论依据 用s表示无界流中的物体表面，来流为均匀流，其未扰动速度或无穷远处的速度为 V??V?xi?V?yj?V?zk, V????1 (2.1.1) 用??x,y,z?表示定常速度势，它在物体外部空间域中适合拉普拉斯方程，在物面上适合不可进入条件，在无穷远处，应该与均匀来流的速度势吻合，即 ?2??0(物体外) (2.1.2) ——————————————————————————————————————————————— ???0(物面s上) (2.1.3) ?n ??xV?x?yV?y?zV?z??(无穷远处) 其中,单位法线向量n指向物体内部。 在速度势?中分出已知的均匀来流项，记 ??xV?x?yV?y?zV?z?? (2.1.4) 这里的?是扰动速度势，?应适合以下定解条件: ?2??0(物体外) (2.1.5) ????V??n(物面s上) (2.1.6) ?n ??0(无穷远处) (2.1.7) 易知过物面s的通量为零，即 ?? ds?0 ???ns 所以远方条件(2.1.7)可进一步具体化为 ??O? ?1?r??) (2.1.8) (?2 ?r? 用rpq表示点p和q之间的距离，对函数??q?和1/rpq在物面s 外部和远方控制面c的内部之空间域内用格林公式，当点p在上述空 间域内时 1 ??p?? ——————————————————————————————————————————————— 4? ???1???1?? ?q??q??????????dsq (2.1.9) ???????nq?rpq??s?c??rpq?nq? 从?的远方条件(2.1.8)可知，c上积分趋于零，式(2.1.9)成为 1 ??p?? 4? ???1?1? ??q?????q?????????rr?n?npqqqs??pq? ??? ds (2.1.10) ???q??? 其中, p是物面s外的任意一点。 在物体的内部域中构造一个合适的内部解?i，它在s内部适合拉普拉斯方程，在物面s上适合某种物面条件，其具体形式将在下面给出。对于上述物体外部的点p函数1/rpq在物体内部域中没有奇点，在内部域中对函数?i(q)和1/rpq用格林公式，得到 ???1?1? 0???????i?q?????i?q??n??r?nq?rpqs??pqq ??? ds (2.1.11) ???q??? 式(2.1.10)和(2.1.11)中的p是物体外部同一个点，把两式相减，——————————————————————————————————————————————— 得到 ?1 4???p????? s??rpq ?????i ????n ?q?nq ? ??????i????nq? ?1??r?pq ?? ??dsq (2.1.12) ???? 在物面s上取?i适合下述两种物面条件，得到两种?i的定解条件，一种是: (s内部)?? ? ?i??(s上)?? (2.1.13) ?2?i?0 定解问题(2.1.13)是拉普拉斯方程的第一类边值问题，它的解是存在且唯一的。取式(2.1.12)中的内部解?i为式(2.1.13)所决定的函数，则式(2.1.12)成为 ??p???? ——————————————————————————————————————————————— s ??q? rpq q (2.1.14) 其中 ??q?? 14? ?????i? ??? (2.1.15) ??n?n? 式(2.1.14)表示扰动速度势?可以用物面s上的分布源表示，其中分布源密度 ?是未知函数，将由扰动势的物面条件(2.1.6)来决定。 物体的附加质量mij，表示物体沿i方向运动引起的j方向的附加质量，公 式如下: mji?????injds?????i sb sb ??j?n ds(i,j?1,2,6) (2.2.1) 利用式(2.1.14)，再结合物面条件，得到 ? 2???p??????q? ——————————————————————————————————————————————— ?np s?? ?1??r?pq ? d??V??np (2.3.1) ??sq? 这就是分布源密度?所适合的线性积分方程。 把积分方程(2.3.1)转换成线性代数方程组，即用离散量代替连续变量。把物面s分成N小块，记 s???sj (2.3.2) j?1N 用平面四边形或三角形来近似代替小曲面?sj。具体做法如下，取第j小块的四个顶点坐标之算术平均值，得到中心点pj的坐标。计算对角线向量的向量积(指向与曲面法线指向相符合)，用nj表示该方向上的单位向量，形成以nj为法线且通过中心点pj的平面，再把四个顶点向该平面作投影，以四个投影点为顶点组成平面四边形?Qj，用?Qj代替原来的小曲面?sj，称?Qj为单元。通常把小范围内的分布源密度?作为常数，因此只要分割不太粗，可以认为?在单元 ?Qj上为常数，记作?j，从而 ? ???q???np?sj ?1??r?pq????dsq??j?1 ????np??Qj??rpq? ——————————————————————————————————————————————— ?dsq (2.3.3) ?? 因此物面s上的积分可以用N个平面四边形(三角形)上积分之和来近似，即 ????q???nps ?1??r?pqN????dsq???j?1 ???j?1?np??Qj??rpq? ?dsq (2.3.4) ?? 上式左端的未知量?(q)是连续型变量，而上式右端的未知量是N个离散量 ?j(j?1N)。为了求解这N个未知数，须要N个方程。取积分方程(2.3.1)中的 动点p为N个单元?Qj的中心点pj(j?1N)，称之为控制点，即控制物面条件使之成立的点。用近似式(2.3.4)代替积分方程(2.3.1)的左端，便可以写出?j的N阶线性代数方程组: j?1N ?aij?j?bi(i?1,2, ,N) (2.3.5) 其中 ? aij??? ?npi?Qj ?1??rpq?i ——————————————————————————————————————————————— ? ?dsq(j?i) ?? aij?2?,bi??V??npi 称aij为影响系数，即第j个单元上的分布源在第i个控制点上的影响。 求解线性代数方程组(2.3.5)得到?j的值以后，便可以得到速度势??p?在控制点pi处的值，即 ??pi???cij?j (2.3.6) j?1 N cij? 另外，物面s上的诱导速度为 ?Qj ?? 1 q (2.3.7) rpiq ???p??2??ini????jvj?pi? j?1 N ?p??si? (2.3.8) 其中??表示求和是不计j?i这一项。，这里的曲面法线n指向物体内部。 三、数值模型 ——————————————————————————————————————————————— 将物体表面划分成四边形面面元，物面为S，每一个四边形面面元为?Sj。为了简化计算，将面网格投影到各自对应的平面上，使曲面网格?Sj变为平面网格?Qj。投影的方法为: 取?Sj四个顶点坐标之平均值，作为中心点pj的坐标。计算?Sj对角线连线 向量的向量积并使得积的方向与流域法向相同。用nj表示该方向上的单位向量。 设: nj?(njx,njy,njz) pj的坐标为: (pjx,pjy,pjz) 则取投影面为过pj并以nj为法向量的平面: njx(x?pjx)?njy(y?pjy)?njz(z?pjz)?0 设在该平面上的投影点为:(x,y,z) 而曲面四边形某个顶点为:(x0,y0,z0) 则有:(x?x0,y?y0,z?z0)//(njx,njy,njz) 因此得到由?Sj顶点坐标求解投影点(?Qj顶点)坐标的线性方程组: ?njx??njy?n?jznjy?njx0njz??x??njxpjx?njypjy?njzpjz??????0?y??njyx 0?njxy0??????njx??z???njzx0?njxz0??? 由此线性方程组可解出投影点坐标。 假设速度势和分布源在?Qj上是不变的，其值为该单元中点(控制点)处的 ——————————————————————————————————————————————— 速度势或分布源。 由于所求速度势和速度等物理量均为物面上的物理量，因此要令p点落在物面之上。式右端分布源的法向导数极限由两部分组成，一部分是P点附近小曲面 另一部分是屋面其余部分s??贡献。当p所趋近于的物面上的点p'作?的贡献， 为控制点的单元，积分时需要考虑奇异性;其余部分为s??。设其中一单元为单元i，其余模型为单元j。对于每一个控制点i，令j?1~N循环一次求得前述方程的积分项(包括奇异积分)。再由i?1~N可以得到N组方程，进而形成求解各个控制点处物理量的矩阵。 ，对于每一个控制点的物理量，通过在?和s??上p点表示控制点(编号i) 积分得到。 即:对于任意一i,,,,(((,，有 ????q???nps N j?1?1??r?pqN????dsq???j?1???j?1?np??Qj??rpq??dsq ???aij?j?bi(i?1,2,, N) ??dsq(j?i) aii?2?,bi??V??npi ?? a将模型控制点数据导入到程序中计ij和bi，可以得到方程组: ?其中， aij????npi?Qj?1??rpq?i ——————————————————————————————————————————————— ?a11???a?N1a1N???1??b1???????????? aNN?????N????bN?? 求解上面线性代数方程组得到?j的值以后，便可以得到速度势??p?在控制点pi处的值，即 ??pi???cij?j 其中，cij? j?1N ?Qj ?? 1 dsq rpiq 每一个i点(控制点)处诱导速度为(是向量): ??1(pi)???jvj(pi) j?1N (pi??si) 其中， 1? ?(?pi??rdsq)(i?j) ?Qjpq? vj(pi)?? 1?lim?(dsq)(i?j)???p?pipi?rQjpq? 在常分布单元假设条件下: p?pi lim?pi(?? 1 ——————————————————————————————————————————————— dsq)?2?nir?Qjpq N 可得， ??1(pi)?2??jni??'?jvj(pi) j?1 (pi??si) 于是便可求解出速度势和物面速度。 四、几何模型 4.1椭球、圆球、有限长圆柱、平行椭球 利用计算机编程来完成以上四类几何模型的建立并划分网格如图4.1、图4.2、图4.3、图4.4所示。其中，椭球长短轴之比为5:1，有限长圆柱柱体长和截面直径之比为5，平行椭球的两个椭球相同且长轴平行，间距为短半轴的3、5、7倍。 图4.1 椭球面网格划分 图4.2 圆球面网格划分 图4.3 有限长圆柱面网格划分 图4.4 平行椭球面网格划分 五、计算参数及结果讨论 5.1椭球 5.1.1椭球附加质量系数M11 经查阅资料可知，椭球附加质量系数M11的理论值为0.059。不同面元数时，M11的计算结果见表5.1。 ——————————————————————————————————————————————— 表5.1椭球的附加质量系数M11随面元数的变化 M11随面元数变化曲线见图5.1。 图5.1 M11随面元数变化曲线 5.1.2椭球附加质量系数M33 椭球附加质量系数M33的理论值为0.894。不同面元数时，M33的计算结果见表5.2。 表5.2椭球的附加质量系数M33随面元数的变化 M33随面元数变化曲线见图5.2。 图5.2 M33随面元数变化曲线 5.1.3 结果讨论 由计算结果可以看出，椭球在附加质量系数M11较M33小很多，这点由椭球几何形状可以容易看出，与实际有较好的符合。随着面元数量的增加附加质量系数M11和M33均更接近于各自理论值，误差都有减小。 在计算过程中发现，面元的划分形式存在质量好坏的区别，对于同样多数量的面元，面元形式的不同会造成不一样的结果，比如同是数量为600的单元，长轴方向划分为20份，周向30份得出的M33为0.964419，而长轴方向划分为30份，周向为20份，得出的附加质量M33为0.9514185，更接近与理论值，误差更小，这是由于对于本细长的椭球模型而言，周方向划分20份已经较为紧密，而细长的长轴方向需要划分更多的数量以使面元更接近于真实物面，所以在划分面元时，需要考虑物体的实际情况，以得到更高质量的面元。 ——————————————————————————————————————————————— 5.2圆球 5.2.1圆球附加质量系数M11 经查资料可知，圆球附加质量系数M11的理论值为0.5。不同面元数时，M11的计算结果见表5.3。 表5.3球的附加质量系数M11随面元数的变化 M11随面元数变化曲线见图5.3。 图5.3 M11随面元数变化曲线 5.2.2圆球附加质量系数M33 由于对称性，圆球附加质量系数M33的理论值也为0.5。不同面元数时，M33的计算结果见表5.4。 表5.4球的附加质量系数M33随面元数的变化 M33随面元数变化曲线见图5.4。 图5.4 M33随面元数变化曲线 5.2.3 结果讨论 由计算结果可以看出，随着圆球面元数的增加，计算而得的附加质量系数逐渐接近解析解，即理论值，误差逐渐减小，由面元法知识可知，这是由于面元数量增加，面元密度增加，单元也就更接近真实物面，从而结果更准确。M11与M33计算结果基本一致，差别较小。 5.3 有限长圆柱 5.3.1圆柱附加质量系数M11 不同面元数时，M11的计算结果见下表。 表5.5有限长圆柱的附加质量系数M11随面元数的变化 ——————————————————————————————————————————————— 有限长圆柱M11随面元数变化曲线图5.5 图5.5 M11随面元数变化曲线 5.3.2圆柱附加质量系数M33 不同面元数时，圆柱M33的计算结果见下表。 表5.6有限长圆柱的附加质量系数M11随面元数的变化 有限长圆柱M33随面元数变化曲线图5.6 图5.6 M33随面元数变化曲线 5.3.3 结果讨论 由计算结果可以看出，对于柱体长与直径长之比为5:1的有限长圆柱体而言，其沿柱体方向附加质量系数M11要小于垂直柱体方向附加质量系数M33，并且随着面元数量增加而增大并且有收敛趋势，而垂直于柱体长轴方向的附加质量系数M33随着面元数量增加而减小并且也呈现出收敛趋势。 5.4平行椭球 5.4.1两椭球平行时附加质量系数M11 每个椭球面元划分数量为400，并且划分形式与单个椭球时一样，以便比较二者的区别。单个椭球在以上计算中得到的附加质量系数M11的值为0.06149，以该值作为计算两椭球并行时相互影响的参考值。 在不同轴间距时时，M11的计算结果与对比如下。 表5.7两椭球并行时M11随长轴间距的变化 M11随长轴间距变化曲线见下图。 ——————————————————————————————————————————————— 图5.7两椭球并行时M11随长轴间距的变化曲线 5.4.2两椭球平行时附加质量系数M33 单个椭球在以上计算中得到的附加质量系数M33的值为0.9675557，以该值作为计算两椭球并行时相互影响的参考值。不同轴间距下的算结果见下表。 表5.8两椭球平行时M33随长轴间距变化 双椭球并行时M33随长轴间距变化时的曲线如下， 图5.8双椭球并行时M33随长轴间距变化曲线 5.4.3 结果讨论 由计算结果可以看出，两个椭球并行时，二者之间会产生干扰，使单个椭球附加质量系数增加，然而其影响会随着椭球之间轴间距增大而减小，以使每个椭球附加质量系数逐渐接近于单个椭球运动时的值，可以推测当其间距足够大时，每个椭球可以看成单独运动的效果。 参考文献 1. 戴遗山,段文洋. 船舶在波浪中运动的势流理论[M]. 国防工业出版社，2008. 2. 陈科. FORTRAN完全自学手册[M]. 机械工业出版社. 2009. 3. 孙寒冰,邹劲,庄家园,吴恭兴.深V型无人艇附加质量计算及其影响[D].2011年 中国智能自动化学术会议论文集(第一分册).2011. 4. 李刚,段文洋,郭志彬.复杂构型潜器附加质量的研究[D].哈尔滨工业大学学 ——————————————————————————————————————————————— 报.2010. 5. 林超友.潜艇近海底航行附加质量数值计算[D].SHIP ENGINEERING,2003 理想流体力学课程设计 应用Hess—Smith方法求附加质量 ———————————————————————————————————————————————