求解势流的正则化快速多极子边界元法
翟杰, 祝宝山 , 曹树良    
清华大学 热能工程系, 水沙科学与水利水电工程国家重点实验室, 北京 100084
摘要:该文将快速多极子算法和处理强奇异积分的正则化算法应用于传统边界元法中, 开发了正则化快速多极子边界元法。该方法既可以解决传统边界元法计算量和存储量会随着单元数量的增加而快速增加的问题, 也可以处理边界元法求解势流速度和速度梯度时产生的强奇异性积分问题。将所开发的方法应用于绕球势流的数值计算中, 计算结果证明了方法的可靠性和高效性; 对相关计算参数影响的分析为复杂边界流动问题的计算提供了参考依据。
关键词势流问题    边界元法    快速多极子法    正则化算法    
Regularization fast multipole boundary element method for solving potential flow problems
ZHAI Jie, ZHU Baoshan , CAO Shuliang    
State Key Laboratory of Hydroscience and Engineering, Department of Thermal Engineering, Tsinghua University, Beijing 100084, China
Abstract: The fast multipole method and the regularization algorithm are combined to process the strong singular integral in the conventional boundary element method for potential flow problems. This method reduces the number of calculations and the storage which increases sharply with the number of elements in the conventional boundary element method. This method can also handle strongly singular integrals for calculating the velocities and the velocity gradient for potential flow by directly differentiating the boundary integral equation. The method is applied to simulate potential flow over a sphere. The results show that this method is accurate and efficient. This model is used to analyze the influence of the calculation parameters for other complicated boundary condition problems.
potential problem    boundary element method (BEM)    fast multipole method (FMM)    regularization algorithm    

根据Helmholtz分解,流体运动的速度u可以分解为势流速度up和涡流速度uω。涡流速度uω可以利用Biot-Savart公式计算求解,而势流速度up则可通过边界元积分来计算。流场的积分计算方法具有良好的自适应性且能够精确处理无穷远边界条件。若对流场中涡量采用Lagrange涡方法计算处理[1],则可以避免流场中网格划分引起的数值粘性,从而能精确求解大Reynolds数的流动[2, 3]

边界元法(boundary element method,BEM)在求解复杂几何形状的势流问题方面具有一定的优势[4, 5, 6],但是边界元在大尺度模型的计算中,当边界单元数增加时,存储空间急剧增加,计算效率降低[7, 8, 9, 10]。另外,采用边界元法求解势流速度以及势流速度梯度时,会产生强积分的奇异性问题[11, 12, 13, 14]

本文提出了一种正则化快速多极子边界元法(regularization fast multipole boundary element method,rFMBEM)计算势流问题,采用N体问题中常用的快速多极子法合并边界元法[15],同时引进可以将强奇异积分转化为半解析形式积分的正则化算法[13]。将所开发的方法应用到绕球势流的计算中,验证了方法的可靠性和高效性;同时,分析了快速多极子参数对计算结果的影响,为复杂问题计算的参数选择提供了参考依据。

1 数学模型 1.1 边界元法

3维势流问题的基本方程可由如下速度势Φ的Laplace方程描述:

$\Delta \Phi = 0.$ (1)

边界条件为:

$\begin{gathered} \Phi = \bar \Phi ,\;\;\;\;\;\;\;\;\;\;\;\;Dirichlet\;边界条件; \hfill \\ q = \frac{{\partial \bar \Phi }}{{\partial n}} = \bar q,\;\;\;\;\;Neumann\;边界条件. \hfill \\ \end{gathered} $

其中:n为边界面的法向单位向量,q为速度在边界面上的法向分量。这里采用Neumann边界条件。

方程(1)的积分解为

$\begin{gathered} \alpha \left( x \right)\Phi \left( x \right) = \int_\Gamma {G\left( {x,y} \right)\frac{{\partial \Phi \left( y \right)}}{{\partial {n_y}}}d{s_y} - } \hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\int_\Gamma \Phi \left( y \right)F\left( {x,y} \right)d{s_y}. \hfill \\ \end{gathered} $ (2)

其中:G=G(x,y)=1/(4π|x-y|)是Laplace方程的3维基本解,F(x,y)=∂G(x,y)/∂ny是基本解的方向导数。边界形状系数α(x)取值如下:

$\alpha \left( x \right) = \left\{ \begin{gathered} 1,\;\;\;\;\;\;x \in 流域内; \hfill \\ \frac{{\theta \left( x \right)}}{{4\pi }},\;\;\;\;x \in 边界; \hfill \\ 0,\;\;\;\;\;\;x \in 流域外. \hfill \\ \end{gathered} \right.$ (3)

这里,θ(x)是点x在边界上展开的立体角,如果点x处处光滑,则θ(x)=2π。

对于边界上势函数的求解采用配置点法,假定x在球表面Γ上,则方程(2)可写成如下的离散形式:

$\begin{gathered} \frac{1}{2}\Phi \left( {{x^i}} \right) + \sum\limits_{e = 1}^{{N_E}} {\int_{{\Gamma _e}} F } \left( {{x^i},y} \right)\sum\limits_{k = 1}^K {{N^k}\left( \xi \right)} {\Phi ^k}\left( {{y^k}} \right)d{\Gamma _e} = \hfill \\ \;\;\;\;\;\;\;\;\;\;\sum\limits_{e = 1}^{{N_E}} {\int_{{\Gamma _e}} {\bar q} G\left( {{x^i},y} \right)} d{\Gamma _e},i = 1,\cdots ,{N_P}. \hfill \\ \end{gathered} $ (4)

其中:K为插值节点数,NE为边界面元数,NP为边界上选择的配置点个数,Nk为插值函数。方程(4)中的两个积分方程分别为具有r-1r-2的Cauchy主值积分,采用子单元法[5]处理其积分奇异性问题。

对方程(4)进行化简整理,可得

$\left[H \right]\left\{ \Phi \right\} = B.$ (5)

求解线性方程组(5)可得边界上配置点的势函数值。对方程(2)求偏导数就可以得到势流场的速度解,

$\begin{gathered} \alpha \left( x \right){u_k}\left( x \right) = \int_\Gamma {\frac{{\partial G\left( {x,y} \right)}}{{\partial {x_k}}}{u_n}} \left( y \right)d\Gamma - \hfill \\ \;\;\;\;\int_\Gamma {\frac{{\partial F\left( {x,y} \right)}}{{\partial {x_k}}}\Phi } \left( y \right)d\Gamma ,k = 1,2,3. \hfill \\ \end{gathered} $ (6)

1.2 快速多极子算法

快速多极子算法的关键是对Green函数基本解的多极子展开[3, 4, 7]

$\frac{1}{{4\pi \left| {x - y} \right|}} = \frac{1}{{4\pi }}\sum\limits_{n = 0}^\infty {\sum\limits_{m = - n}^n {{{\bar S}_{n,m}}\left( {x - {y_c}} \right){R_{n,m}}\left( {y - {y_c}} \right)} } .$ (7)

式中,yc是场点y的展开中心,|yyc|<|xyc|。

两个双谐函数Rn,m(x)和Sn,m(x)可分别表示为:

$\begin{gathered} {R_{n,m}}\left( x \right) = \frac{1}{{\left( {n + m} \right)!}}P_n^{\left| m \right|}\left( {\cos \alpha } \right){e^{im\beta }}{\rho ^n},\hfill \\ {S_{n,m}}\left( x \right) = \left( {n - m} \right)!P_n^{\left| m \right|}\left( {\cos \alpha } \right){e^{im\beta }}{\rho ^{\frac{1}{{n + 1}}}}. \hfill \\ \end{gathered} $

这里,(ρ,α,β)是x在球坐标系下的坐标。Pnm为伴随Legendre多项式,它是n次Legendre多项式Pnm阶导数。

将方程(7)代入原积分项,可得:

$\begin{gathered} \;\;\;\;\;\;\;\;\;\int_{{\Gamma _e}} q \left( y \right)G\left( {{x^i},y} \right)d{\Gamma _e} = \hfill \\ \frac{1}{{4\pi }}\sum\limits_{n = 0}^\infty {\sum\limits_{m = - n}^n {{{\bar S}_{n,m}}} } \left( {x - {y_c}} \right){M_{n,m}}\left( {{y_c}} \right),\hfill \\ \end{gathered} $ (8)

$\begin{gathered} \;\;\;\;\;\;\;\;\;\int_{{\Gamma _e}} F \left( {{x^i},y} \right)\Phi \left( y \right)d{\Gamma _e} = \hfill \\ \frac{1}{{4\pi }}\sum\limits_{n = 0}^\infty {\sum\limits_{m = - n}^n {{{\bar S}_{n,m}}} } \left( {x - {y_c}} \right){{\tilde M}_{n,m}}\left( {{y_c}} \right). \hfill \\ \end{gathered} $ (9)

其中:

${M_{n{\text{,}}m}}\left( {{y_c}} \right) = \int_{{\Gamma _e}} q \left( y \right){R_{n{\text{,}}m}}\left( {y - {y_c}} \right)d{\Gamma _e}{\text{,}}$ (10)

${\tilde M_{n,m}}\left( {{y_c}} \right) = \int_{{\Gamma _e}} {\frac{{\partial {R_{n,m}}\left( {y - {y_c}} \right)}}{{\partial {n_y}}}} \Phi \left( y \right)d{\Gamma _e}.$ (11)

同样,可以对x的展开中心xL进行局部展开[7]。当展开中心变化时,可以采用中心转移的多级子展开公式,而不需要重新计算积分[7]

经过多极子展开,方程(5)的左边可表示为

$\begin{gathered} \;\;\;\;\;\sum\limits_{i = 1}^{{N_E}} {{H_{ji}}{\Phi _i} = \sum\limits_q {\sum\limits_{i \in {W_q}} {{H_{ji}}} } } {\Phi _i} + \hfill \\ \frac{1}{{4\pi }}\sum\limits_{n = 0}^\infty {\sum\limits_{m = - n}^n {{{\bar S}_{n,m}}\left( {x - {y_c}} \right){{\tilde M}_{n,m}}\left( {{y_c}} \right)} .} \hfill \\ \end{gathered} $ (12)

方程(5)的右边可表示为

${B_i} = \sum\limits_{i \in {W_q}} {{G_i}{q_i}} + \frac{1}{{4\pi }}\sum\limits_{n = 0}^\infty {\sum\limits_{m = - n}^n {{{\bar S}_{n,m}}\left( {x - {y_c}} \right){M_{n,m}}\left( {{y_c}} \right)} .} $ (13)

从式(12)和(13)可以看出,通过快速多极子算法,满秩矩阵已转化为带状稀疏矩阵,其计算效率降低为O(N)。数值计算过程中,采用自适应性八叉树结构对单元点进行结构划分[15]

1.3 奇异积分计算的正则化处理

在求解势流场的速度解时,方程(6)具有r-3r-5的奇异积分,当源点靠近球表面时就会产生强积分奇异性问题。采用子单元正则化算法来求解近边界点的速度。

图1所示,四边形等参边界元可近似分割成两个三角形等参元△123及△134。以三角形等参元△123为例,在其表面上建立参考正交坐标系ξoη。设源点x在三角形平面上的垂足为x0,并以x0为原点建立一个极坐标系ρθ,极轴初位置与平行,则有坐标转换关系:

$\left\{ \begin{gathered} \xi = {\xi _0} + \rho \cos \theta ,\hfill \\ \eta = {\eta _0} + \rho \sin \theta . \hfill \\ \end{gathered} \right.$ (14)

图1 三角形等参元示意图

三角形等参元的形函数和几何坐标系之间的关系如下[8, 9]

$\begin{gathered} \;\;\;\;\;{N_m}\left( {\xi ,\eta } \right) = \frac{1}{{aA}}\left( {{a_m} + {b_m}\xi + {c_m}\eta } \right) = \hfill \\ \frac{1}{{2A}}\rho \left( {{b_m}\cos \theta + {c_m}\sin \theta } \right) + {N_m}\left( {{\xi _0},{\eta _0}} \right),\hfill \\ \end{gathered} $ (15a)

$\begin{gathered} \;\;\;\;\;{y_i} = {N_m}\left( {\xi ,\eta } \right){y_{mi}} = \hfill \\ \frac{1}{{2A}}\rho \left( {{b_m}\cos \theta + {c_m}\sin \theta } \right){y_{mi}} + {x_{oi}},\hfill \\ \end{gathered} $ (15b)

其中:A是三角形等参元△123的面积;ambmcm是关于ξiηi的参数,i,m=1,2,3。

由式(15b)可得:

$\begin{gathered} {r_i} = {y_i} - {x_i} = \frac{1}{{2A}}\rho \left( {{b_m}\cos \theta + {c_m}\sin \theta } \right){y_{mi}} + {x_{oi}} - {x_i},\hfill \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{r^2} = {\rho ^2} + e_1^2. \hfill \\ \end{gathered} $

式中,e1=|xox|为源点至三角形等参元△123的最小距离。三角形等参元的面积微分为

$d\Gamma = d\xi d\eta = \rho d\rho d\theta .$

将上式代入方程(13),则各积分项可以归纳为

${I_n} = \int_{{\Gamma _e}} {\frac{{{Q_n}\left( {\rho ,\theta } \right)}}{{{r^n}}}\rho d\rho d\theta ,n = 1,3,5.} $ (16)

这里,Qn(ρ,θ)是关于ρ,θ的多项式函数[13]

e1较小时,Gauss数值积分失效,较小的e1是引起积分奇异性的主要原因。对方程(16)做关于变量ρ的积分,

${K_n}\left( {\rho ,\theta } \right) = \int {\frac{{{Q_n}\left( {\rho ,\theta } \right)}}{{{r^n}}}\rho d\rho ,n = 1,3,5.} $ (17)

对方程(17)反复运用分部积分,得到Kn(ρ,θ)的解析解。因此,方程(16)可写成

${I_n} = \int_{{\Gamma _e}} {{K_n}\left( {\rho ,\theta } \right)\left| \begin{gathered} \rho = \rho _2^{\left( \theta \right)} \hfill \\ \rho = \rho _1^{\left( \theta \right)} \hfill \\ \end{gathered} \right.d\theta .} $ (18)

其中:ρ1(θ)和ρ2(θ)是由三角形等参元3边的极坐标的表达式决定。此时,原强奇异积分已转化为关于变量θ的一系列线积分,且引起几乎奇异积分的因子e1也已被分离到线积分核外,可用Gauss积分计算这些线积分[13, 14]

2 数值结果与分析 2.1 单球势流

球的半径为r,均匀来流沿x轴方向速度为U,如图2所示,则球坐标下速度势的解析解为

$\phi = U\left( {r + \frac{1}{2}\frac{{{r^3}}}{{{R^2}}}} \right)\cos \theta .$ (19)

图2 数学模型示意图

对应的势流速度为:

$\begin{gathered} {u_r} = U\left( {1 - \frac{{{r^3}}}{{{R^3}}}} \right)\cos \theta ,\hfill \\ {u_\theta } = - U\left( {1 + \frac{1}{2}\frac{{{r^3}}}{{{R^2}}}} \right)\sin \theta ,\hfill \\ {u_\varphi } = 0. \hfill \\ \end{gathered} $ (20)

图3给出球表面的边界单元划分。从方程(19)可知,球绕流问题的速度势可分为均匀势和绕流势,为方便分析数值结果,这里只考虑绕流势。数值计算中aU都取量纲量为1。图4给出了φ=0平面内球表面绕流速度势函数分布,数值计算中的边界单元数为3 540。从图4中可以看出,BEM和rFMBEM的计算结果与理论解完全一致。

图3 球边界单元划分
图4 φ=0平面内球表面速度势函数分布

表1给出了正则化算法和直接Gauss积分法得到的势流速度的误差比较。定义接近度为到球边界距离d与球半径r的比值。从表1中可以看出,当d/r>0.1时,两种方法的精度相近。当d/r<0.1时,Gauss积分法的计算精度急剧恶化。当d/r<0.001 时,Gauss积分法已经无法计算出准确值,但正则化算法还可以得到很好的计算结果(O(10-2))。

表1 单球绕流问题计算精度比较
与球边界的
接近度(d/r)
计算误差/%
直接Gauss积分正则化算法
1.001.011.01
0.4641.011.01
0.2151.021.02
0.1001.871.87
0.046 436.522.14
0.021 5653.950.31
0.010 04 941 .231.87
0.002 1529 085 .511.28
0.001 0036 452 .071.83
0.000 4642.51
0.000 2153.66
0.000 1003.50
0.000 0105.21
0.000 0011.14

图5给出了BEM和rFMBEM在不同边界单元数下计算时间的比较。当边界单元数小于1 000 时,rFMBEM计算时间更长;当边界单元数在 1 000 至2 000之间时,两种数值方法计算效率相似;但随着边界单元数的增加,BEM的计算时间成指数性增长,而rFMBEM的计算时间成线性增长,当边界单元数大于3 000 时,rFMBEM比BEM的计算效率快得多。

图5 BEM和rFMBEM计算时间比较

图6为多极子展开系数p对计算效率的影响,实线表示在不同边界单元数下p为5~40的计算时间,虚线表示在不同边界单元数下p为5~40球表面势函数的均方根误差$\sigma {\text{ = }}\sqrt {{{\sum {\left( {{\varphi _i} - \bar \varphi } \right)} }^2}/N} $。可以看出,相同边界单元数下,计算时间随p值的增加而成指数增长,计算精度随p值的增加而提高。综合考虑两方面的影响,可以确定p的取值范围为15~30。

图6 多极子开展级数p的影响

图7为八叉树结构中每个单元内最大节点数s对计算效率的影响,实线表示在不同边界单元数下s为1~30的计算时间,虚线表示在不同边界单元数下s为1~30球表面势函数的均方根误差。可以看出,相同边界单元数下,计算时间随s值的增加而先降低再增加,s=20是计算时间的转捩点,而计算精度随s值的增加而降低。综合考虑两方面的影响,可以确定s的取值范围为10~20。

图7 最大节点数s的影响
2.2 多球势流

采用本文方法计算多球绕流的势流问题,这里选定p=20和s=15。图8给出φ=0平面内的双球绕流和4球绕流的速度分布图,其中水平方向和垂直方向的球心距离均为4a

图8 φ=0平面内多球绕流的速度分布图

图9给出BEM和rFMBEM在多球绕流问题中计算时间和势函数均方根误差的比较。可以看出,rFMBEM方法可以有效地加快计算速度,尤其是大边界单元数的情况,当边界单元数大于20 000时,BEM方法已超出本文所用计算机的计算能力,而rFMBEM方法仍然有效。从图9中给出的rFMBEM方法对于BEM方法得到球表面势函数结果的统计均方根误差可以看出,两种方法的结果相近(O(10-3))。

图9 BEM和rFMBEM计算多球绕流的数值结果比较
3 结 论

本文针对3维势流问题,提出了计算存储量小且计算速度快的正则化快速多极子边界元法。从理论上给出了数值方法的推导过程,并将方法应用到求解绕球势流问题中。从数值计算结果来看,本文方法具有以下优点:1) 加快了传统边界元法的计算速度;2) 解决了势流速度求解时产生的强奇异积分问题。在后续的工作中,将把本文开发的方法和已开发的快速3维涡方法相结合,研发无网格的快速积分法计算复杂边界的不可压缩3维流动问题。

本文还分析比较了快速多极子算法中计算参数对计算时间和精度的影响,从数值结果来看:1) 多极子展开系数取值范围为15~30;2) 每个单元内最大节点数取值范围为10~20。这些结论可以为复杂形状的绕流问题以及流体机械内部流动问题提供参数选取的参考依据。

参考文献
[1] Zhu B S. Finite volume solution of the Navier-Stokes equations in velocity-vorticity formulation [J]. International Journal for Numerical Methods in Fluids, 2005, 48(6): 607-629.
[2] 祝宝山, 王旭鹤, 龟本乔司, 等. 流体机械非定常流动的涡方法数值模拟 [J]. 水力发电学报, 2011, 30(5): 178-185. ZHU Baoshan, WANG Xuhe, Kamemoto K, et al. Numerical simulation of unsteady flows of fluid machinery using discrete vortex method [J]. Journal of Hydroelectric Engineering, 2011, 30(5): 178-185. (in Chinese)
[3] Zhu B S, Wang H, Wang L B, et al. Three-dimensional vortex simulation of unsteady flow in hydraulic turbines [J]. International Journal for Numerical Methods in Fluids, 2012, 69(10): 1679-1700.
[4] Gharakhani A, Ghoniem A F. BEM solution of the 3D internal Neumann problem and a regularized formulation for the potential velocity gradients [J]. International Journal for Numerical Methods in Fluids, 1998, 24(1): 81-100.
[5] 蔡瑞瑛, 曾昭景, 黄文龙. 边界元法程序设计及工程应用[M]. 南京: 江苏科学技术出版社, 1996.CAI Ruiying, ZENG Zhaojing, HUANG Wenlong. Boundary Element Method Program and Engineering Application [M]. Nanjing: Jiangsu Science and Technology Press, 1996. (in Chinese)
[6] 姚振汉, 王海涛. 边界元法[M]. 北京: 高等教育出版社, 2010.YAO Zhenhan, WANG Haitao. Boundary Element Method [M]. Beijing: Higher Education Press, 2010. (in Chinese)
[7] Liu Y J, Nishimura N. The fast multipole boundary element method for potential problems: A tutorial [J]. Engineering Analysis with Boundary Elements, 2006, 30(5): 371-381.
[8] Liu Y. Fast Multipole Boundary Element Method: Theory and Applications in Engineering [M]. Cambridge, UK: Cambridge University Press, 2009.
[9] 宁德志, 滕斌, 勾莹. 快速多极子展开技术在高阶边界元方法中的实现[J]. 计算力学学报, 2006, 22(6): 700-704. NING Dezhi, TENG Bin, GOU Ying. Implementation of the fast multipole expansion technique in the higher order BEM [J]. Chinese Journal of Computational Mechanics, 2006, 22(6): 700-704. (in Chinese)
[10] 宁德志, 滕斌, 勾莹. 快速多极子边界元法在三维势流问题中应用[J]. 大连理工大学学报, 2005, 45(2): 243-247. NING Dezhi, TENG Bin, GOU Ying. Application of fast multipole boundary element method to 3-D potential flow problem [J]. Journal of Dalian University of Technology, 2005, 45(2): 243-247. (in Chinese)
[11] Liu Y, Rudolphi T. Some identities for fundamental solutions and their applications to weakly-singular boundary element formulations [J]. Engineering Analysis with Boundary Elements, 1991, 8(6): 301-311.
[12] Mukherjee S. CPV and HFP integrals and their applications in the boundary element method [J]. International Journal of Solids and Structures, 2000, 37(45): 6623-6634.
[13] 周焕林, 牛忠荣, 王秀喜. 三维位势问题边界元法中几乎奇异积分的正则化[J]. 计算物理, 2005, 22(6): 501-506. ZHOU Huanlin, NIU Zhongrong, WANG Xiuxi. Regularization of nearly singular integrals in the boundary element method for 3D potential problems [J]. Chinese Journal of Computational Physics, 2005, 22(6): 501-506. (in Chinese)
[14] 周焕林, 牛忠荣, 王秀喜. 位势问题边界元法中几乎奇异积分的正则化[J]. 应用数学和力学, 2003, 24(10): 1069-1074.ZHOU Huanlin, NIU Zhongrong, WANG Xiuxi. Regularization of nearly singular integrals in the boundary element method of potential problems [J]. Applied Mathematics and Mechanics, 2003, 24(10): 1069-1074. (in Chinese)
[15] Greengard L, Rokhlin V. A new version of the fast multipole method for the Laplace equation in three dimensions [J]. Acta Numerica, 1997, 6(1): 229-269.