圆球势流绕流的有限元方法

黄 璐，陈 立，邱辽原，许 辉

(中国舰船研究设计中心，湖北 武汉430064)

0 引 言

随着CFD和计算机的飞速发展及广泛应用，许多流体力学的复杂问题都通过计算机数值模拟进行解决，有限元法就是众多数值计算方法的一种［1］。

在工程实践中，钝体绕流是一个十分典型的问题。船舶与海洋工程中就有很多例子可以简化为钝体绕流的模型进行研究，如推进器，海洋平台的立管和柱体绕流等。圆球绕流是钝体绕流的典型例子。

圆球绕流的问题已有许多研究，吴望一等［2］利用有限元法和传统的分析方法相结合，研究了无界域上理想流体的圆球绕流问题。王吉飞和万德成［3］利用有限元法和并行计算相结合，研究了理想不可压流体的二维和三维钝体势流绕流问题，给出了二维无限域圆柱绕流和三维无限域圆球绕流的数值模拟结果，并做了并行效率的分析。岳蕾、张志国等［4］利用大涡模拟等湍流模型，分析了圆球绕流的尾部流场及圆球表面的压力变化。Hanazaki［5］通过差分方法求解三维非稳态N－S 方程，得到一些圆球绕流结果(Re ＜200)。Sungsu Lee［6］用有限元方法模拟Re 数在100 ～500 之间的圆球绕流。

本文利用Fortran 语言自编程序，采用有限元法求解圆管内圆球在理想不可压缩条件下的流动问题，将求解区域划分成离散的网格，采用1，0.5，0.1，0.05 等不同的网格密度，比较不同的网格密度对有限元法求解圆球绕流问题的影响，为有限元求解圆球绕流问题时采用合适的网格密度提供参考。

1 问题的提出

1.1 有限元法简介

有限元法利用泛函变分原理或方程余量与权函数正交化原理。利用“分块逼近”进行离散求解［7］，将流场的求解区域划分成互不重叠的单元;在每个单元内，选择若干点作为求解函数的插值点。单元中的函数将由线性化的基函数取代，然后通过单元体的积分就可得到单元有限元方程，经过累加获得总体有限元方程，通过求解总体有限元方程得到节点上的函数值，从而求解得到需要的物理量。

1.2 问题的物理模型

考虑位于圆管内无限长的圆球体的绕流问题，几何尺寸如图1所示。无穷远处来流为Vz=1，Vy=0。由于流场具有上下左右的轴对称性，取圆球中心处的纵向剖面，且只考虑左上角1/4 的计算区域a－b－c－d－e，把它作为有限元的求解区域Ω。

图1 求解问题的物理模型Fig.1 Physical model of the problem

坐标系按照图中所示建立，无穷远处来流速度和压强分别为V∞=1，P∞=0，分别求解网格密度为1，0.5，0.1，0.05 下的流场流函数和压强分布。

1.3 物理问题满足的微分方程式

由于圆管内流动具有上下、前后对称特点，计算区域Ω 取流场的1/4 区域;Γ1是Ω 的本质边界，由进口边界a－e，固壁边界e－d，圆球面固壁边界a－b－c 组成;Γ2是自然边界，它是前后对称轴线cd。

1.4 边界条件

区域Ω和相关边界条件值如下:

2 有限元数值求解

2.1 积分表达式

对流函数ψ 满足的Stokes 方程采用Galerkin 方法，可得到:

应用Green－Gauss 公式可将上述积分改写为:

其中:Ω 为(r，z)平面上的求解区域;Γ2为自然边界。

2.2 单元剖分

根据物理问题的特点以及区域的形状，把计算区域分成许多几何形状规则但大小可以不同的单元，确定单元节点的数目和位置，建立表示网格的数据结构。采用的单元形状和节点的分布，以及插值函数的选取还应考虑到计算精度和可微性的要求。本题中根据给定的网格密度划分相应的网格，采用三角形结构化网格。为使形成的总体有限元方程组系数矩阵具有较小的带宽，区域节点序号沿着区域短边方向，自下而上逐列编号。每个单元节点序号采用逆时针的编号方法，最后列出总体序号与单元序号之间的表格。

除此之外，单元剖分还要建立本质边界条件和自然边界条件的节点表。自然边界条件对解此题并无贡献，由网格表特征可知，本题的本质边界条件为a－b－c 段上，ψ 为0;d－e 段上，ψ 为2;e－a段上，ψ 为r2/2。

2.3 确定单元基函数

单元及结点数确定后，单元基函数也确定了。本算例中采用三角形三结点的单元基函数，单元基函数的线性插值函数为:

三结点分别对应3 个基函数。记e 单元的三角形的3 个结点坐标为;则式(4)所满足的插值条件为:

将式(5)代入式(4)，即得到含有9 个待定系数ai，bi，ci(i=1，2，3)的9 个代数方程组:

其中下标i，j，k 按1，2，3 顺序循环取值。求解得到:

其中:

式中A(e)为e 单元三角形的面积。

2.4 单元分析

注意到有如下变形式:

其中(zc，rc)为单元三角形的形心坐标。将式(4)和式(10)代入式(9)，即可推导出有限元方程:

其中:

2.5 方程求解

将各个单元的方程按照顺序合成后形成总体方程，在形成的总体方程中利用消行修正法处理本质边界条件，得到修正方程。利用Gauss 列主元素消元法直接求解方程组。求出所有的待求量后，便得到了近似函数的表达式，并可以计算出相关的物理量。

各个单元速度的求解:

求出单元的速度后，进而可以得到各个单元节点的速度，通过Bernouli 方程计算出各节点的压力值，假设求解区域位于同一水平面内，介质密度ρ=1 来流压力P=0，那么结点压力:

3 计算结果及分析

3.1 流线图和压力云图

通过编制程序，得到点间隔1，0.5，0.1，0.05 的网格相关数据节点个数和单元个数，如表1所示。

表1 各个网格密度节点和单元数Tab.1 Node and element number of each mesh density

不同网格密度得到的流线图形如图2所示，不同网格密度得到的压力云图如图3所示。

图2 网格密度流线图Fig.2 Flow pattern when the mesh density equals

图3 网格密度压力图Fig.3 Pressure profile when the mesh density equals

3.2 不同网格密度对压力的影响

为了比较不同的网格密度对计算结果的影响，选取压力P 在不同网格密度下，分别在y=0，y=1，y=2 处计算结果的比较，比较结果如图4所示(图中t 代表网格密度)。

3.3 流线图和压力云图

从结果输出图中可以得到以下分析结果:

1)从流线图中可以看到，随着网格密度的增加，流线变得更加光滑，更加均匀，说明网格密度越密得到的流场流线图越接近真实的流场。从图2中可以看到，y=0 时是0 流线，随着y 值的增加，流线的值也在增加，在y=2 时达到最大。

2)从压力云图中可以看到，在圆球前端，压力最大。随着圆球的壁面向上，压力值慢慢变小，直达圆球顶端，压力达到最小。这与根据伯努利定理分析速度矢量得到的结果基本符合。

3)由图4 可以看出，在网格密度为1和0.5时，由于网格较稀疏，相关物理量的变化还比较大;在网格密度为0.05和0.1 时，相关物理量的变化比较小，说明随着网格的加密计算趋于稳定。

图4 y=0，1，2 时不同网格密度压力P 变化图Fig.4 Pressure profile of different mesh density when y=0，1，2

4 结 语

以上结果与文献［3］中计算结果比较，符合度较好，说明文中所用的有限元法适合于解决相关的流体力学问题，为不同密度网格的有限元法在解决流体力学问题研究时提供参考。以后可以进一步研究粘性流动下有限元解法，也可同CFD 模拟作进一步的比较。

［1］章本照．流体力学中的有限元方法［M］．北京:机械工业出版社，1986．

［2］吴望一，梁镭，温功碧，等．无界区域上理想绕流问题的混合元方程［J］．空气动力学学报，1989，7(4):449－454．

［3］王吉飞，万德成．钝体势流绕流的有限元并行计算［C］．第二十一届全国水动力学研讨会暨全国第八届水动力学学术会议，成都，2008:227－233．

［4］岳蕾，张志国，蒋奉兼．圆球绕流的大涡模拟分析研究［C］．第十一届全国水动力学学术会议暨第二十四届全国水动力学研讨会，无锡，2012:237－244．YUE lei，ZHANG Zhi-guo，JIANG Feng-jian．Investigation of flow cross sphere using large eddy simulation［C］．The 11st Session of the National Seminar on Hydrodynamics and National Eighth Hydrodynamics Academic Conference，Wuxi，2012:237－244．

［5］HANAZAKIH A．A numerical study of three-dimensional stratified flow past a sphere［J］．Journal of Fluid Mechanics，1998，192:393－419．

［6］SUNGSU L．A numerical study of the unsteady wake behind a sphere in uniform flow at moderate Reynolds numbers［J］．Computer and Fluids，2000，29:639－667．

［7］章本照，印建安，张宏基．流体力学数值方法［M］．北京:机械工业出版社，2003．