哈尔滨工程大学理想流体力学大作业

理想流体力学大作业学生姓名：学号：2013年10月Hess—Smith方法计算物体附加质量作者：摘要：本文运用Hess-Smith方法计算了圆球、椭球和圆柱的附加质量系数以及椭球并行的干扰效应。同时，文章分析了网格变化对计算值的影响趋势。本文使用matlab语言对圆球、椭球与圆柱的模型进行了网格有限元的划分，得到各个单元的节点坐标，然后利用Hess-Smith方法对圆球、椭球及并行椭球的附加质量系数进行计算及分析。关键字：边界元；Hess-Smith;附加质量系数一、物理背景Hess-Smith方法是一种计算任意三维物体势流的方法，该方法由美国的Hess和Smith两人于20世纪60年代提出。Hess-Smith方法又称为分布奇点法，作为一种边界元方法，它用许多平面四边形或三角形表面单元来表示物体表面，并在每个单元上布置强度未知的源，然后在物体表面的某些考察点上满足法向速度为零的物面边界条件，得到求单元源密度的线性代数方程组。求解方程组得到源密度分布，进而可求流场内任意点的速度、压力等物理量。二、理论依据2.1分布源模型的建立s为无界流中的物体表面，来流为均匀流，在无穷远处流体的速度为：xyzVViVjVk(2.1.1)1VVx,y,z为定常速度势，并在物体外部空间域中满足拉普拉斯方程，在物面上适合不可穿透条件，在无穷远处，应该与均匀来流的速度势相同。即20（物体外）(2.1.2)0n（物面s上）(2.1.3)其中，单位法线向量n指向物体内部。在速度势中分出已知的均匀来流项，记xyzxVyVzV(2.1.4)这里的是扰动速度势，应适合以下定解条件：200(Vnnb（在物体外）（在物面s上)无穷远处）(2.1.5)用Rpq表示点p和点q之间的距离，根据格林第三公式，当p点位于物面s外部和远方控制面c的内部之空间域时，有如下公式：1()11()[()()]4qpqqpqSqpqdsnRnR(2.1.6)由远方边界条件可知，远方封闭控制面cs上的积分趋于零，从而上式化为：1()11()[()()]4bqpqqpqSqpqdsnRnR(2.1.7)又由式（2.1.5）可得：01111()()()()44bbqpqpqSSpqdsVndsnRR(2.1.8)得到混合分布模型，为了得到单一分布模型表示的扰动势，在物体内部域中构造一个合适的内部解i。于上述物体外部的点P，函数1/pqR在物体内部域中没有奇点，在物体内部域中对函数i和1/pqR用格林第三公式，得：()110()()()bbiiqpqqpqSSqqdsdsnRnR(2.1.9)将（2.1.7）与（2.1.11）相减，得：114(){[]()[]}BiiqpqqqpqSpdsnRnnR(2.1.10)取下式定解条件中的i：20i(在bs内部)(2.1.11)i(在bs)(2.1.12)则式(2.1.10)成为：()()BpqSqpdsR(2.1.13)其中，1()()4iqnn(2.1.14)2.2分布源密度的求解式(2.1.13)中右端分布源的法向导数极限由两部分组成，一部分是p点附近小曲面ε的贡献，另一部分是物面其余部分的贡献。法向指向取向物体内部，小曲面ε的贡献为2πσ（p）,则有如下关系式：()12qspppqsppqdnnr(2.2.1)再结合物面条件(2.1.5)，得到12qspppqspqdVnnr(2.2.2)这就是分布源密度所适合的线性积分方程。把积分方程(2.2.2)转换成线性代数方程组，即用离散量代替连续变量。把物面s分成N小块，记1Njjss(2.2.3)用平面四边形或三角形来近似代替小曲面js。具体做法如下，取第j小块的四个顶点坐标之算术平均值，得到中心点jp的坐标。计算对角线向量的向量积（指向与曲面法线指向相符合），用jn表示该方向上的单位向量，形成以jn为法线且通过中心点jp的平面，再把四个顶点向该平面作投影，以四个投影点为顶点组成平面四边形jQ，用jQ代替原来的小曲面js，称jQ为单元。通常把小范围内的分布源密度作为常数，因此只要分割不太粗，可以认为在单元jQ上为常数，记作j，从而qpqQpjqpqpsdsrndsrnqjj11(2.2.4)因此物面s上的积分可以用N个平面四边形（三角形）上积分之和来近似，即qpqQpjNjqpqpsdsrndsrnqj111(2.2.5)上式左端的未知量()q是连续型变量，而上式右端的未知量是N个离散量(1)jjN。为了求解这N个未知数，须要N个方程。取积分方程(2.2.2)中的动点p为N个单元jQ的中心点(1)jpjN，称之为控制点，即控制物面条件使之成立的点。用近似式(2.2.5)代替积分方程(2.2.2)的左端，便可以写出j的N阶线性代数方程组：1(1,2,,)NijjijabiN(2.2.6)当计算出影响系数ija后，即可解线性方程组得到分布源密度。2.3速度势与附加质量的求解根据速度势在控制点pi处的值，由公式：1()Niijjjpc(2.3.1)1qjijspiqQcdr(2.3.2)根据2.2得到的分布源密度j，求解线性方程组（2.2.7）可得速度势的值。物体的附加质量ijm，表示物体沿i方向运动引起的j方向的附加质量，公式如下:(2.3.3)12iimV(2.3.4)根据所求得的速度势的值可计算处附加质量的值。三、数值模型及参数计算3.1数值模型要求解流场中物体表面的速度势分布，需要先将物体的外表面进行网格划分。经过网格划分以后，原来的物体连续外表面被离散为N×M个相对独立的小平面，这些小平面构成了求解该问题的数值模型。Hess-Smith的基本思想是将连续曲面的积分离散为小单元来简化计算，其计算思路核心在于解该方程组：ijjiab，通过求解线性方程得到j,1,2,,ijNM。对于不同的计算目的，只需要改变控制面条件，即改变ib来实现。得到j后，进而由1NMiijjjPC求及附加质量，其中：jijiijssmdSndSn为求ijm，令pjbVn，则求得j。(,1,2,6)jjiijisbsbmndsdsijn故而1111111NMiiismndSPnPS，同理可得33m。3.2参数计算3.2.1计算ijn对于球面：,,ijijijopnPxyzop对于椭球面：设1111,,nxyz，2222,,nxyz，则易知12yyy，12zzz，21211xxxaa，故21,,ijijxyzanop。（a=1即为球）对于圆柱面：在圆柱顶面上1,0,0ijn，同理在其底面上有1,0,0ijn，侧面上则有0,cos(),sin()22ijn3.2.2计算ija对于上述几何模型中的四边形上均匀分布奇点的诱导速度公式计算，首先将四边形上的分布源密度取为1，四边形四个顶点逆时针方向排列，顶点pi的坐标为（ξi，ηi，0），其中坐标系oξηζ的原点为四边形形心，平面oξη即为四边形平面，则需计算如下三个积分式：111XQyQzQSddxrSddyrSddzr各积分式的计算公式为：,141,1111,1lniiiiiiiixiiiiirrlSlrrl,141,1111,1lniiiiiiiiyiiiiirrlSlrrl4,1,11111(arctanarctan)iiiiiiiiziiimchmchSzrzr在本几何模型中，需要将世界坐标系转换到局部坐标系中，查阅相关文献可知，当给定物体上的3个不共线的点,,1,2,3iiiiPxyzi，即可建立局部坐标系。其中局部坐标系的原点为点1P，X轴的正向为12PP的方向，Z轴的正向为1213PPPP的方向，Y轴的正向为121312PPPPPP的方向。设,,XYZ轴的单位矢量分别为,,,ijk它们在世界坐标系中的方向余弦分别用123,,1,2,3iiiuuui表示。对X轴，111213,,uuu可按下式给出：1211121312PPiuiujukPP，i，j，k为世界坐标系3根轴的单位矢量。31112112223211323331111121131211221231311321330001uuuuuuTuuuxuyuzuxuyuzuxuyuzu矩阵1T就是由世界坐标系到物坐标系的坐标变换矩阵。经过坐标转换可得每个网格单元的Sx、Sy、Sz，然后得到世界坐标系下的S，即ija由于编程实现上述过程非常繁复，现采用一种更为简便方法的求解ija当ij时，31ijiiiijjPPpijqjppqPPQrnadSSnrr当ij时，2ija，从而求得ija的系数矩阵。3.2.3计算ijC由1ijijqpqQCdSr可得当ij时，1ijijjppCSr当ij时，如图，设iP位于jQ上，取一小圆ε，转换为极坐标积分：212001122iRqpqIdSdrdrRSrr（广义积分收敛）故可近似122ijjCS，得到系数矩阵ijC。3.2.4计算m11、m33令1()ibVnP，则1111111NMiiiismndSPnPS同理求得33m。对于双椭球体，如下图所示，只需将第二个椭球的单元标记为（1NM，…2NM）即可。方程组的形式不变，在计算m11或m33时分别对两个椭球面的单元进行积分即可分别得到m11a、m11b、m33a、m33b。四、几何模型对于连续的积分方程，通常的数值处理是转换成线性代数方程组，即用离散量代替连续变量。在几何建模上，就是将物面s分成N小块，记为1ssNjj采用平面四边形或三角形近似代替小曲面sj。4.1球面网格划分4.1.1网格划分使用球坐标对球体进行网格划分，取球体半径R=1，以球心O为原点建立坐标系如下图所示取xoz平面上的圆，圆上任意一点与x轴的夹角范围为0-π，划分N份，每相邻两份夹角为δφ=π/N，在球体上以球心为圆心，划分N条纬线。纬线为在球体上的同心圆，取任意同心圆，圆半径为r，将圆划分M份，每相邻两份之间夹角为δθ=2π/M,N个同心圆都划分M份，形成M条经线。由纵横交错的线将球体划分形成网格。取球体任意网格交点A,设坐标为A(x，y，z),由球坐标可将A坐标表示为：cossincossinsinxRyRzR在matlab中根据以上网格划分原理划分网格，得到的球体几何模型为：4.1.2节点坐标取球体面上任意一个小网格，在小网格ABCD上，A标记为A(i，j)，即纵向第i个点，横向第j个点，则(,)(,)(,,)(cos,sincos,sinsin)AijfRRRyxz