基于PETSc的区域分解法数值模拟圆柱势流绕流问题

第九届全国水动力学学术会议暨第二十二届全国水动力学研讨会文集基于ＰＥＴＳｃ的区域分解圆柱势流绕流问题

王吉飞万德成（．１：海交通人学船舶海洋‘ｊ建筑Ｔ程学院海洋Ｔ程国家蕈点实验室２０００３０）摘要：本研究基于有限元法和并行计算工具箱ＰＥＴｓｃ（Ｐｏｎａｂｌｅ，ＥＸｔｅｎｓｉｂｌｅＴ００ｌｋｉｔｆｏｒ

ｓｃｉｅｎｔｉｆｉｃｃｏｍｐｕｔａｔｉｏｎ），开发和编制区域分解法程序，对理想不可压流体中的圆柱绕流问题进

行了并行数值计算，验证了网格的分辨率、计算精度、并行加速比，以及不同区域划分形式对数值计算效果的影响，计算得到的流场速度分布和物体表面压力分布结果与解析解比较，两者吻合很好。数值计算结果表明，基于ＰＥＴｓｃ的区域分解法和有限元法结合可以高效率实现理想不可压流体圆柱势流绕流问题的并行数值计算，为今后开展复杂黏性流体流动问题的并行计算奠定了基础．关键词：势流；有限元法；钝体绕流；ＰＥＴｓｃ；区域分解；并行计算

ｌ引言

圆柱绕流在实际ｊ【程应用中非常普遍。在船舶与海洋－』：程中，洋流对柱形桥墩的冲击，海洋平台立柱、立管绕流等均为圆柱绕流问题，因此研究圆柱绕流具有重要的实用意义。当前计算能力不足是ＣＦＤ发展所面临的一火挑战。为解决：ｆ程应用方面的大规模设计计算问题，并行化是大幅度提高计算效率的最有效手段。区域分解法是应用较广泛的一种粗粒度并行计算方法。在流体力学有限元计算中，区域划分可在网格生成之前进行，也可在网格生成之后进行，后者的优点是可以方便地对串行程序作并行改造而无需对网格生成形式进行修改，也适合于网格生成器与流场求解器分离的并行计算。并行计算１：具箱ＰＥＴＳｃ（ＰｏｒｔａｂＩｅ，ＥｘｔｅｎｓｉｂｌｅＴｏｏｌｋｉｔｆｏｒＳｃｉｅｎｔｊｆｉｃＣｏｍｐｕｔａｔｉｏｎ）为线性和非线性代数方程组的并行计算提供了可移植、可扩展的并行库，非常适合偏微分方程离散后人型稀疏代数方稗组的并行数值求解。本研究基丁．有限元法和并行计算：Ｌ具箱ＰＥＴｓｃ，开发和编制区域分解法程序，数值求解理想不可压流体中的圆柱绕流问题。首先对整个计算流域的网格进行区域划分，各子区域由重替过渡区域连接，然后对各子区域网格用有限单元离散，单元形式为任意人面体单元，单元形函数为Ｌａｇｒａｎｇｅ插值函数，有限元积分离散后的线性代数方程组，应用ＰＥＴＳｃ库函数（ＧＭＲＥＳ方法加ｌＬＵ预估条件）进行求解，并使用ＭＰｌ（ＭｅｓｓａｇｅＰ嬲ｓｉｎｇ

Ｉｎｔｅ而ｃｅ）库函数

对各子区域间并行计算数据进行通信。本研究结构安排如下：首先描述圆柱势流绕流的区域分解数值算法，然后给出数值模拟结果，并进行并行性能分析，最后给出结论。

．２７０．第九届全国水动力学学术会议暨第二十二届全国水动力学研讨会文集２数值算法２．１有限元空间离散对于三维理想不可压流体的流动，可求解以速度势缈为未知数的Ｌａｐｌａｃｅ方程边值问题：ｆ鲁＋窖＋鲁：ｏ（训㈡∈力

Ｉ良２。勿２‘瑟２。ｐ∥产７‘“

｛＾（１）

ｈ项矧。＝ｇ

其中Ｑ为计算流体区域，Ｓ为本质边界，＆为自然边界。方程（１）的弱解积分表达式为：蜷挈＋等掣＋警掣瑚＝协啦㈤

加。苏锄却却瑟瑟。坞一７～…

本研究使用任意六面体单元离散计算流域，采用八结点的Ｌａ可绷ｇｅ线性插值函数进行单元分析。进行单元分析时，为了使不同形状和尺寸的任意六面体单元均能使用规则的插值函数进行单元分析，首先要进行坐标变换，将笛卡尔坐标（ｘ，少，ｚ）变换为无量纲的局部坐标（孝，７，ｆ），即将任意六面体单元变换为边长为２的正方体单元（一１≤孝≤１，一１≤，７≤ｌ，一１≤ｆ≤１）。取插值函数为：

瞬＝亡（１＋参考）（１＋仍７７）（１＋缶ｆ）（ｆ＝１，…，８）（３）

其中（参，仍，毒）是ｆ号结点的坐标值，如图１所示。

幽ｌ边长为２的正方体标准单元单元内任意一点的势函数值可表示为：矽＝∑哆哆，令却＝嘭，则在任一单元内弱解积分表达式可转化为：ｋｃ警等＋詈等＋警等，出妙出卜＝如田衄

㈤

应用不同坐标系的微分量变换关系，单元分析方程中系数矩阵的生成可采用Ｇａｕｓｓ求积公式进行近似的数值积分。Ｇａｕｓｓ求积公式的表达式为：

．２７１．第九屑全国水动力学学术会议暨第二十二届全国水动力学研讨会文集ｆ。ｆ。上，厂（ｆ靠ｆ）ｄ笋碉ｆ＝喜喜喜ｗ＿Ｍ厂（当，矾，幺）（５）

其中（专，协，缶）是积分点的坐标值，Ｗ，ｗ，，Ⅵ是相应积分点三个坐标的权系数，Ｚ是孝方向的积分点个数，朋是刁方向的积分点个数，门是ｆ方向的积分点个数。对丁二自然边界面上的积分，即单元分析方程的矗端项，可通过坐标变换将任意四边形变换为边长为２的正方形，然后采ＨｊＧａｕｓｓ求积公式进行近似的数值积分。当单元分析完成后即可进行总体合成，最后对总体有限元方稗进行本质边界条件修正。

３．２区域分解法本研究使用重叠型区域分解算法进行并行计算。考虑在区域Ｑ上二阶椭圆型方程的边值问题：Ｉ三“＝／（“∈力）ＩＩ一

抛ｌ（６）

…ｓ刮’瓦ｌ。邓

Ｌ’’ｚ将区域刃分解为Ⅳ个相互重叠的子区域｛岛，…，ｑ，…，ｑ），Ⅳ个子区域中非重叠的

区域依次为｛岛，…，岛，…，ｑ），重叠的区域则依次为｛Ｑ＋１＇…，Ｑ＂…，钇＋川）。对于子区域乌（１≤，≤Ⅳ），ｑ＝ｑＵｑ＋，．。｜『＞１ＵＱ＋，Ｉ阳（，＝ｆ），其本质边界为Ｓ，，Ｓ，∈Ｓ，其自然边界为是，，是，∈是，其入．１：内边界分别为‘一（歹＞１）和厂，＋（歹＜Ⅳ），如图２所示。

ｂｖｅｒｌａｐｐｉｒＩＩ）ｖｅｒｌａｐｐｉ啦

呻∞’ｌａｐｐｉ珈

ｒｆ，｜Ｉ，．ｒ｛，．１Ｉ－Ｌ．ｒｌＭ）．ｒ‘＾ｌ，．

ｎ＂ｌｎ．Ｑ－．ＪＱ…Ｑ＿．ｆ．Ｉ

’’、＼ｌ，—一

，，，’、、、＼ｌ／

＿＿＿一一

Ｑ．ｏ，．Ｉ

图２区域分解计算模型

记“；为第七次迭代后乌上的近似解，甜ｊＩｒ，一（，＞１）和“；ｌｒ。（，＜Ⅳ）为第七次迭代后日

内边界上的边界值，“‘为第七次迭代后整个求解域Ｑ的近似解，下面给出区域分解的并行算法：首先，选择初始近似内边界值“？ｈ（，＞１）和”？Ｉｒ，＋（，＜Ⅳ）；其次，对七＝１，２，３…，在第，号进程上求解区域乌上的边值问题：

．２７２．第九届全国水动力学学术会议暨第二十二届全国水动力学研讨会文集￡“；＝厂（“；∈Ｑ，）巩氓乱～㈩

“ｊｈ＝“㈡ｒ，．

（，＞１）

“‰，＋＝“㈡ｒ，＋（，三Ⅳ）

该边值问题使用有限元方法进行离散，并应用ＰＥＴＳｃ库函数（ＧＭＲＥＳ方法加ＩＬＵ预估条件）进行求解。有限元离散后的系数矩阵和右端项向量均可重复使用，对于不同的尼值，仅需进行内边界的本质边界条件处理即可。求解后，通过进程间的数据传递，取得整个计算域Ｑ上的解，并更新各子区域的内边界值：Ｉｏ．＝甜；Ｉｑ。Ｑ，（１≤ｆ≤门，，＝ｆ）

Ｉｑ“＝ｏ・５书‘誓；，Ｉｑ“ｃｎ，＿“鼻・Ｉｎ。“ｃ。，＋．）（１≤‘＜刀，，＝７）（８）Ｉｒ，一＝ｏ．５母（“ｊ叫Ｉｒ，－＋＂卫ｌＩｒ，一ｎｑ．。）（，＞１）

ｈ＝ｏ．５幸（“㈡ｒ，＋＋“欠ｌｌｒ，＋ｎＱ，＋．）（Ｊ＜Ⅳ）

依次循环得到序列｛甜１，甜２，“３…），当七专∞，“‘一致收敛，即为方程式的解。３数值模拟结果３．１三维等截面圆柱势流绕流本研究首先模拟三维等截面圆柱势流绕流，其计算模型如图３所示。计算流域的取值范围是：ｘ∈【一ｌｏ，ｌｏ】，ｙ∈【一ｌｏ，ｌｏ】，ｚ∈卜２０，２０】。圆柱沿ｚ轴方向，半径为ｌ。均匀来流速度

为甜＝ｌ，１，＝Ｏ，ｗ＝０。边界条件的定义如下：当ｘ：１ｏ时，万：ｌｏ；当ｘ：一１０时，业：一ｌ；锄

当ｙ＝一１０，ｙ＝ｌｏ，ｚ＝一２０，ｚ＝２０或ｘ２＋ｙ２＝１时，娑：ｏ。计算网格由ｇａｍｂｉｔ软件ａ玎生成，结点数为１０８３２２，任意六面体单元数为１０２１６０。

图４为三维等截面圆柱绕流势函数值分布剖面。图５为圆柱绕流的速度大小分布剖面。图６为圆柱绕流的压力值分布剖面。

．２７３．

矿矿嘭畔第九届全国水动力学学术会议暨第＿－二十一二届全国水动力学研讨会文集删３圆柱绕流模型图４速度势分布剖面图５速度大小的分布剖面闰６慊力的分布剖面为了验证结果的准确性，将本研究的计算结果与解析解进行比较。图７为ｚ＝０平面内沿，轴的述度分布与解析解的比较．可见数值解与解析解吻合良好。幽８为ｚ＝Ｏ平面内圆柱表面压力分布与解析解的比较，两者也吻合良好。

誊ｊ－鞋

“ｌ

Ｘ删∞ｔｖｅｌ。ｃｌｔ㈣ｓｏ“

幽７速度的数值解２口解析解的比较凹８抹山的数值解’，解析解的比较３１三维变截面圆柱势流绕流本研究对变截面圆杆势流绕流进行了数值模拟．圆柱沿：轴方向，最小截面圆半径为Ｏ４，

最人截面圆半释为ｌ６，计算模删如图９所示。边界条什与等截面｜丑｜韩势流绕流相同。计算网格

由髓ｍｂｉｔ软１，｜＝生成，结点数为６８８８０，任意Ａ面体单元数＿：ｌ『６４０００。图】Ｏ所示为三维变截面【县｜杜绕流势函数值分布剖面。幽１ｌ所示为变截面圆柱绕流的速度人小分布剖面。幽１２所示为变截面圆柱绕流的压力值分布剖面。对ｒ变截面圆＿｜＝｛！势游ｉ绕流，流动早现明显的二维特性，图１３所示为变截面圆柱表面及并剖面的ｚ方向速度分量分布。图１４为变截面劂柱附近速度父量分布，同样可以看到变截厩恻柱绕流的二维特性。