Boundary element method with parameterized elements for problems of potential flow velocity field
-
摘要:目的 边界元法在海洋工程水动力学中有着广阔的应用前景，为推广边界元法在海洋工程水动力学中的应用，方法 根据边界积分法建立积分方程，采用参数化单元边界元法对势流问题进行求解，得出流场速度势。对经典算例进行数值计算，与数学解析解比较，并进行误差分析。在二维问题下，分别采用非连续参数化单元和参数化单元边界元法求解势流速度场问题；在三维问题下，采用参数化单元边界元法求解势流速度场问题。结果 结果显示，在二维问题下，采用非连续参数化单元边界元法求解势流问题具有较高的精度和效率，可以在采用较少单元数的情况下得到较为理想的数值解；在三维问题下采用参数化单元边界元法虽然计算速度较快，并可以得到较好的平均相对精度，但有些点误差较大，需要改进算法或使用其他单元进行求解。结论 参数化单元边界元法在求解海洋工程势流问题时，数值计算实现过程更简洁，可发展成为求解船舶兴波等船舶水动力学问题的通用方法。Abstract:Objectives The Boundary Element Method (BEM) has broad application prospects in ocean engineering hydrodynamics. In order to promote the application of BEM in ocean engineering hydrodynamics,Methods the integral equation is established according to boundary integral method and a parameterized element BEM is adopted. This meta-method solves the potential flow problem and obtains the velocity potential of the flow field. The numerical calculations are performed on the basis of classic examples, the mathematical solutions are compared and error analysis is performed. Under the two-dimensional problem, the discontinuous parameterization element BEM and parameterized element BEM are used to solve the potential flow velocity field. Under the three-dimensional problem, a parameterized element BEM is used to solve the potential flow velocity field problem.Results The results show that under the two-dimensional problem, the discontinuous parameterized element BEM is used to solve the potential flow problem with high precision and efficiency, and the ideal numerical solution can be obtained with fewer elements. The parameterized element BEM under the three-dimensional problem is faster in calculation and can obtain better average relative accuracy, but some points have large margins of error and require the improvement of the algorithm or the introduction of another element to be solved.Conclusions When the parameterized element BEM is used to solve the potential problem in ocean engineering, the numerical calculation implementation process is more concise and can be developed into a general method for solving the hydrodynamic problems of ships, such as ship motion.
-
0. 引言
在势流理论下，海洋工程中的许多问题都可以归结为边界值问题，如波浪能发电装置、海洋石油钻井平台，以及船舶与波浪、水流的相互作用等。在现有的求解方法中：采用解析法求解精确、计算速度快，但对求解域的几何形状要求较为苛刻，只适用于一些几何形状较为简单的规则求解域，比如特征函数匹配展开法适用于截面为矩形的浮体与波浪的相互作用[1]，多极子法适用于截面为圆形的浮体与波浪的作用[2]；数值方法中的有限元法可以对复杂形状进行求解，但需要对整个求解域进行离散，计算量较大[3-4]。
相比于需要对整个求解域进行离散的方法，如有限元法、有限差分法等，边界元法具有显著的优势[5]。边界元法最显著的特点之一是在计算时更小的计算量和数据存储；此外，边界元法的数值精度通常要优于有限元法。采用边界元法仅需在求解域边界上进行离散，将三维问题转化为二维问题，二维问题转化为一维问题，即能很方便地处理无限域问题，这也是边界元法在水动力学问题中能得到广泛应用的原因之一[6]。
本文将分别采用二维下的非连续参数化单元和参数化单元边界元法解势流速度场问题，以便在较低的单元数下得到较为理想的求解精度，并采用一种三维下的参数化单元边界元法解势流速度场问题，以便快速得到可以接受的平均误差精度。
1. 数学模型及积分方程建立
1.1 数学模型
势流问题中的流体为无旋、无粘性、不可压缩的理想流体，流场域满足拉普拉斯方程：
∂2ϕ∂x2+∂2ϕ∂y2=0 (1) 流场域满足相应的边界条件如下：
\begin{array}{*{20}{c}} {\frac{{\partial \phi }}{{\partial \mathit{\boldsymbol{n}}}} = {f_1},}&{\phi = {f_2}} \end{array} (2) 式中：n为物面某点的法向向量，垂直于物面向外；ϕ为速度势；f为给定的函数，下标1，2表示不同的给定函数(以下同)。
式(2)中的法向导数由式(3)定义：
\frac{{\partial \phi }}{{\partial \mathit{\boldsymbol{n}}}} = {n_x}\frac{{\partial \phi }}{{\partial x}} + {n_y}\frac{{\partial \phi }}{{\partial y}} (3) 式中，nx，ny分别为物面单位法向向量在x，y轴上的分量，法向量方向垂直物面向外。需注意的是，单位法向量在不同的分量上是不同的，其是关于x和y的函数。求解域的控制方程(1)和边界条件已知，便构成了边界值问题，且存在特解。
1.2 积分方程的建立与离散
控制方程[7]存在基本解：
\phi \left( {x,y} \right) = \frac{1}{{4{\rm{ \mathsf{ π} }}}}\ln \left[ {{{\left( {x - \xi } \right)}^2} + {{\left( {y - \eta } \right)}^2}} \right] (4) 式中，ξ和η为选定点的坐标。根据格林公式，容易得边界积分方程
\begin{array}{l} \lambda \left( {\xi ,\eta } \right)\phi \left( {\xi ,\eta } \right) = \int\limits_C {\left[ {\phi \left( {x,y} \right)\frac{\partial }{{\partial \mathit{\boldsymbol{n}}}}\phi \left( {x,y;\xi ,\eta } \right) - } \right.} \\ \;\;\;\;\;\;\left. {\phi \left( {x,y;\xi ,\eta } \right)\frac{\partial }{{\partial \mathit{\boldsymbol{n}}}}\phi \left( {x,y} \right)} \right]{\rm{d}}S\left( {x,y} \right) \end{array} (5) 式中：S为曲线积分；C为积分路径。
在边界积分方程中，λ(ξ, η)的定义如下[8]：
\lambda \left( {\xi ,\eta } \right) = \left\{ {\begin{array}{*{20}{c}} \begin{array}{l} 0,\\ 1/2,\\ 1, \end{array}&\begin{array}{l} \left( {\xi ,\eta } \right)\;在求解域之外\\ \left( {\xi ,\eta } \right)\;在求解域边界上\\ \left( {\xi ,\eta } \right)\;在求解域内 \end{array} \end{array}} \right. (6) 将求解域的边界L使用非连续参数化单元近似为L1+L2+…+LN，各离散单元以逆时针顺序排列于求解域的边界上。离散单元L1，L2，…，LN由求解域边界上的点(x1, y1)，(x2, y2)，…，(xN, yN)分隔，且N为有限值。各离散单元的边界条件由式(2)，得
\begin{array}{*{20}{c}} {\frac{{\partial \phi }}{{\partial \mathit{\boldsymbol{n}}}} = f_1^k,}&{\phi = f_2^k} \end{array} (7) 每个单元的端点坐标分别是(xk, yk)和(xk+1, yk+1)，且需注意，此处的(xN+1, yN+1)=(x1, y1)。
贯穿整个非连续参数化单元，采用单元上点的f取值来近似fk。为了进行线性近似，需要在单元上的2个不同点取fk的近似值。这里选取2个点(ξk, ηk)和(ξN+k, ηN+k)，它们与点(xk, yk)，(xk+1, yk+1)之间的长度均为τlk。其中，lk为单元k的长度，τ为一给定系数(0＜τ＜1/2)。单元k的具体结构参见图 1。
点(ξk, ηk)，(ξN+k, ηN+k)处的ϕ值分别由ϕk和ϕN+k指代。根据Ang[9]的研究，为了用ϕk和ϕN+k近似表示出单元上ϕ值的线性变化，定义如下：
s\left( {x,y} \right) = \sqrt {{{\left( {x - {x^k}} \right)}^2} + {{\left( {y - {y^k}} \right)}^2}} (8) 其中，(x, y)属于Lk(Lk为第k个离散单元)。由式(8)可知，s(x, y)为Lk上的点(x, y)与端点(xk, yk)间的距离。注意，点(ξk, ηk)，(ξN+K, ηN+k)的s(x, y)值分别为τlk和(1-τ)lk。所需ϕ(x, y)值的线性近似可以以ϕk，ϕN+k的形式给出，见式(9)。类似地，有p= \frac{\partial }{{\partial \mathit{\boldsymbol{n}}}}ϕ(x, y)]，见式(10)。
\phi \left( {x,y} \right) \cong \frac{{\left[ {s\left( {x,y} \right) - \left( {1 - \tau } \right){l^k}} \right]{\phi ^k}}}{{\left( {2\tau - 1} \right){l^k}}} - \frac{{\left[ {s\left( {x,y} \right) - \tau {l^k}} \right]{\phi ^{N + k}}}}{{\left( {2\tau - 1} \right){l^k}}} (9) \begin{array}{*{20}{c}} {\frac{\partial }{{\partial \mathit{\boldsymbol{n}}}}\left[ {\phi \left( {x,y} \right)} \right] \cong \frac{{\left[ {s\left( {x,y} \right) - \left( {1 - \tau } \right){l^k}} \right]{p^k}}}{{\left( {2\tau - 1} \right){l^k}}} - }\\ {\frac{{\left[ {s\left( {x,y} \right) - \tau {l^k}} \right]{p^{N + k}}}}{{\left( {2\tau - 1} \right){l^k}}}} \end{array} (10) 当0＜τ＜1/2时，近似表达式(9)、式(10)定义为非连续单元。在对求解域边界L进行离散时，任意元素都有ϕk和pk，且其中有且仅有一个已知，如边界上ϕk已知，pk则未知。因此，式(9)、式(10)中有2N个未知项。将式(9)和式(10)代入式(5)，可得
\begin{array}{*{20}{c}} {\lambda \left( {\xi ,\eta } \right)\phi \left( {\xi ,\eta } \right) = \sum\limits_{k = 1}^N {\frac{1}{{\left( {2\tau - 1} \right){l^k}}}\left\{ {{\phi ^k}\left[ { - \left( {1 - \tau } \right){l^k}F_2^k\left( {\xi ,\eta } \right) + } \right.} \right.} }\\ {\left. {F_4^k\left( {\xi ,\eta } \right)} \right] + {\phi ^{N + k}}\left[ {\tau {l^k}F_2^k\left( {\xi ,\eta } \right) - F_4^k\left( {\xi ,\eta } \right)} \right] - }\\ {{p^k}\left[ { - \left( {1 - \tau } \right){l^k}F_1^k\left( {\xi ,\eta } \right) + F_3^k\left( {\xi ,\eta } \right)} \right] - }\\ {\left. {{p^{N + k}}\left[ {\tau {l^k}F_1^k\left( {\xi ,\eta } \right) - F_3^k\left( {\xi ,\eta } \right)} \right]} \right\}} \end{array} (11) 其中，
F_1^k\left( {\xi ,\eta } \right) = \int\limits_{{C^k}} {\phi \left( {x,y;\xi ,\eta } \right){\rm{d}}s\left( {x,y} \right)} (12) F_2^k\left( {\xi ,\eta } \right) = \int\limits_{{C^k}} {\frac{\partial }{{\partial \mathit{\boldsymbol{n}}}}\phi \left( {x,y;\xi ,\eta } \right){\rm{d}}s\left( {x,y} \right)} (13) F_3^k\left( {\xi ,\eta } \right) = \int\limits_{{C^k}} {s\left( {x,y} \right)\phi\left( {x,y;\xi ,\eta } \right){\rm{d}}s\left( {x,y} \right)} (14) F_4^k\left( {\xi ,\eta } \right) = \int\limits_{{C^k}} {s\left( {x,y} \right)\frac{\partial }{{\partial \mathit{\boldsymbol{n}}}}\phi\left( {x,y;\xi ,\eta } \right){\rm{d}}s\left( {x,y} \right)} (15) 式中，Ck为单元积分路径。
为了求解式(11)右侧的2N个未知项，建立了一个2N阶线性代数方程组，即可以依次取点(ξ, η)为(ξm, ηm)，其中m=1, 2, …, 2N，则又有式(16)。由图 1，可知点(ξk, ηk)和(ξN+k, ηN+k)为单元内的两点，则λ(ξm, ηm)=1/2，m=1, 2, …, 2N。式(16)又可写为式(17)，式中的zk，zN+k表示元素上的未知量。其中在边界单元上，当ϕk已知时，amk，am(N+k)，bmk分别由式(18)、式(19)和式(20)给出，δmk为Kronecker数。在边界单元上，当pk已知时，情况与上述类似，此处不再赘述。Fpk(ξm, ηm)(其中p=1, 2, 3, 4)系数的参数化计算方法参见文献[9]。
\begin{array}{*{20}{c}} {\frac{1}{2}{\phi ^m} = \sum\limits_{k = 1}^N {\frac{1}{{\left( {2\tau - 1} \right){l^k}}}\left\{ {{\phi ^k}\left[ { - \left( {1 - \tau } \right){l^k}F_2^k\left( {{\xi ^m},{\eta ^m}} \right) + } \right.} \right.} }\\ {\left. {F_4^k\left( {{\xi ^m},{\eta ^m}} \right)} \right] + {\phi ^{N + k}}\left[ {\tau {l^k}F_2^k\left( {{\xi ^m},{\eta ^m}} \right) - F_4^k\left( {{\xi ^m},{\eta ^m}} \right)} \right] - }\\ {{p^k}\left[ { - \left( {1 - \tau } \right){l^k}F_1^k\left( {{\xi ^m},{\eta ^m}} \right) + F_3^k\left( {{\xi ^m},{\eta ^m}} \right)} \right] - }\\ {\left. {{p^{N + k}}\left[ {\tau {l^k}F_1^k\left( {{\xi ^m},{\eta ^m}} \right) - F_3^k\left( {{\xi ^m},{\eta ^m}} \right)} \right]} \right\}}\\ {m = 1,2, \cdots ,2N} \end{array} (16) \sum\limits_{k = 1}^N {\left( {{a^{mk}}{z^k} + {a^{m\left( {N + k} \right)}}{z^{N + k}}} \right)} = \sum\limits_{k = 1}^N {{b^{mk}}} (17) 其中，
\begin{array}{*{20}{c}} {{a^{mk}} = {{\left( {2\tau - 1} \right)}^{ - 1}}\left[ { - \left( {1 - \tau } \right)F_2^k\left( {{\xi ^m},{\eta ^m}} \right) + } \right.}\\ {\left. {{{\left( {{l^k}} \right)}^{ - 1}}F_4^k\left( {{\xi ^m},{\eta ^m}} \right)} \right] - \frac{1}{2}{\delta ^{mk}}} \end{array} (18) \begin{array}{*{20}{c}} {{a^{m\left( {N + k} \right)}} = {{\left( {2\tau - 1} \right)}^{ - 1}}\left[ {\tau F_2^k\left( {{\xi ^m},{\eta ^m}} \right) - } \right.}\\ {\left. {{{\left( {{l^k}} \right)}^{ - 1}}F_4^k\left( {{\xi ^m},{\eta ^m}} \right)} \right] - \frac{1}{2}{\delta ^{\left( {m - N} \right)k}}} \end{array} (19) \begin{array}{*{20}{c}} {{b^{mk}} = - {p^k}{{\left( {2\tau - 1} \right)}^{ - 1}}\left[ {\left( {1 - \tau } \right)F_1^k\left( {{\xi ^m},{\eta ^m}} \right) - } \right.}\\ {\left. {{{\left( {{l^k}} \right)}^{ - 1}}F_3^k\left( {{\xi ^m},{\eta ^m}} \right)} \right] - {p^{N + k}}{{\left( {2\tau - 1} \right)}^{ - 1}}\left[ { - \tau F_1^k\left( {{\xi ^m},{\eta ^m}} \right) + } \right.}\\ {\left. {{{\left( {{l^k}} \right)}^{ - 1}}F_3^k\left( {{\xi ^m},{\eta ^m}} \right)} \right]} \end{array} (20) 2. 二维势流问题算例
2.1 矩形势流场算例
第1个算例，考虑矩形势流场问题，定解问题如下：
\left\{ \begin{array}{l} \frac{{{\partial ^2}\phi }}{{\partial {x^2}}} + \frac{{{\partial ^2}\phi }}{{\partial {y^2}}} = 0;x \in \left[ {0,5} \right],y \in \left[ {0,10} \right]\\ \phi \left( {x,0} \right) = 0\\ \phi \left( {0,y} \right) = 0\\ \phi \left( {x,10} \right) = 100\sin \left( {{\rm{ \mathsf{ π} }}x/100} \right)\\ \frac{\partial }{{\partial x}}\phi \left( {5,y} \right) = 0 \end{array} \right. (21) 由文献[10]可知，上述问题的解析解为
\phi \left( {x,y} \right) = \frac{{100\sin \left( {{\rm{ \mathsf{ π} }}x/100} \right){\rm{sh}}\left( {{\rm{ \mathsf{ π} }}y/100} \right)}}{{{\rm{sh}}\left( {\rm{ \mathsf{ π} }} \right)}} (22) 如图 2所示，在分析计算中，将矩形计算域边界离散为30个单元，将边界积分方程离散，可以得到一个线性方程组，解此方程组即可求得所需的速度势或其他未知量。
在划分30个单元、采用非连续参数化单元离散的情况下，速度势的最大相对误差为0.63%，平均相对误差为0.19%，具体计算结果如表 1所示；在划分单元数加倍、采用连续常数单元的情况下，速度势的最大相对误差为70.86%，平均相对误差为8.42%，计算结果如表 2所示。
表 1 非连续参数化单元速度势计算结果(案例1)Table 1. The calculation results of velocity potential by discontinuous parameterized elements(case 1)节点 x y 数值速度势计算结果 解析速度势计算结果 相对误差/% 1 0.1 0.1 0.008 537 0.008 546 0.11 2 0.1 2.5 0.236 313 0.236 265 0.02 3 0.1 5.0 0.628 773 0.625 917 0.46 4 0.1 7.5 1.422 434 1.421 926 0.04 5 0.1 9.9 3.062 776 3.043 561 0.63 6 4.9 0.1 0.271 407 0.271 940 0.20 7 4.9 2.5 7.534 049 7.518 070 0.21 8 4.9 5.0 19.896 945 19.917 008 0.10 9 4.9 7.5 45.318 735 45.246 430 0.16 10 4.9 9.9 96.962 292 96.847 673 0.12 11 0.2 0.1 0.017 073 0.017 084 0.06 12 2.0 0.1 0.159 467 0.159 921 0.28 13 2.5 0.1 0.192 659 0.192 385 0.14 14 4.0 0.1 0.258 137 0.258 758 0.24 15 4.9 0.1 0.271 407 0.271 940 0.20 16 0.2 9.9 6.096 510 6.084 118 0.20 17 2.0 9.9 57.112 596 56.953 737 0.28 18 2.5 9.9 68.502 724 68.515 455 0.02 19 4.0 9.9 92.406 620 92.153 083 0.28 20 4.9 9.9 96.962 292 96.847 673 0.12 平均误差 0.19 表 2 参数化单元速度势计算结果(案例1)Table 2. The calculation results of velocity potential by parameterized elements(case 1)节点 x y 数值速度势计算结果 解析速度势计算结果 相对误差/% 1 0.1 0.1 0.002 490 0.008 546 70.86 2 0.1 2.5 0.236 344 0.236 265 0.03 3 0.1 5.0 0.625 979 0.625 917 0.01 4 0.1 7.5 1.420 772 1.421 926 0.08 5 0.1 9.9 3.985 676 3.043 561 30.95 6 4.9 0.1 0.329 797 0.271 940 21.28 7 4.9 2.5 7.534 036 7.518 070 0.21 8 4.9 5.0 19.975 335 19.917 008 0.29 9 4.9 7.5 45.382 206 45.246 430 0.30 10 4.9 9.9 96.152 946 96.847 673 0.72 11 0.2 0.1 0.014 657 0.017 084 14.21 12 2.0 0.1 0.160 352 0.159 921 0.27 13 2.5 0.1 0.192 864 0.192 385 0.25 14 4.0 0.1 0.258 801 0.258 758 0.02 15 4.9 0.1 0.329 797 0.271 940 21.28 16 0.2 9.9 6.487 487 6.084 118 6.63 17 2.0 9.9 56.897 571 56.953 737 0.10 18 2.5 9.9 68.448 599 68.515 455 0.10 19 4.0 9.9 92.072 400 92.153 083 0.09 20 4.9 9.9 96.152 946 96.847 673 0.72 平均误差 8.42 2.2 圆形势流场算例
第2个算例，考虑圆形势流场问题，定解问题如下：
\frac{{{\partial ^2}\phi }}{{\partial {x^2}}} + \frac{{{\partial ^2}\phi }}{{\partial {y^2}}} = 0 (23) \phi \left| {_{r = {r_0}}} \right. = r_0^2\cos \theta \sin \theta (24) 式中:r，θ为极坐标中的变量；r0为圆形流场域的半径，此处r0=1。其解析解为：
\phi \left( {r,\theta } \right) = 0.5{r^2}\sin \left( {2\theta } \right) (25) 如图 3所示，在分析计算中，将圆形计算域边界离散为32个单元，将边界积分方程离散，亦可得到一个线性方程组，解此方程组即可求得所需的速度势或其他未知量。
在划分32个单元、采用非连续参数化单元离散的情况下，速度势的最大相对误差为0.02%，平均相对误差为0.01%，计算结果如表 3所示；在划分单元数加倍、采用连续常数单元的情况下，速度势的最大相对误差为0.97%，平均相对误差为0.08%，计算结果如表 4所示。
表 3 非连续参数化单元速度势计算结果(案例2)Table 3. The calculation results of velocity potential by discontinuous parameterized elements(case 2)节点 x y 数值速度势计算结果 解析速度势计算结果 相对误差/% 1 0.7 0.7 0.490 070 0.490 000 0.01 2 0.6 0.6 0.360 050 0.360 000 0.01 3 0.5 0.5 0.250 035 0.250 000 0.01 4 0.4 0.4 0.160 022 0.160 000 0.01 5 0.3 0.3 0.090 013 0.090 000 0.01 6 0.2 0.2 0.040 006 0.040 000 0.01 7 0.1 0.1 0.010 001 0.010 000 0.01 8 0 0 0.000 000 0.000 000 0.00 9 -0.1 -0.1 0.010 001 0.010 000 0.01 10 -0.2 -0.2 0.040 006 0.040 000 0.01 11 -0.3 -0.3 0.090 013 0.090 000 0.01 12 -0.4 -0.4 0.160 022 0.160 000 0.01 13 -0.5 -0.5 0.250 035 0.250 000 0.01 14 -0.6 -0.6 0.360 050 0.360 000 0.01 15 -0.7 -0.7 0.490 070 0.490 000 0.01 16 0.966 0.259 0.249 980 0.249 949 0.01 17 0.866 0.499 0.432 967 0.432 913 0.01 18 0.707 0.707 0.500 085 0.499 990 0.02 19 0.499 0.866 0.432 967 0.432 913 0.01 20 0.174 0.985 0.170 891 0.170 863 0.02 平均误差 0.01 表 4 参数化单元速度势计算结果(案例2)Table 4. The calculation results of velocity potential by parameterized elements(case 2)节点 x y 数值速度势计算结果 解析速度势计算结果 相对误差/% 1 0.7 0.7 0.490 009 0.490 000 0.00 2 0.6 0.6 0.360 007 0.360 000 0.00 3 0.5 0.5 0.250 005 0.250 000 0.00 4 0.4 0.4 0.160 003 0.160 000 0.00 5 0.3 0.3 0.090 002 0.090 000 0.00 6 0.2 0.2 0.040 001 0.040 000 0.00 7 0.1 0.1 0.010 000 0.010 000 0.00 8 0 0 0.000 000 0.000 000 0.00 9 -0.1 -0.1 0.010 000 0.010 000 0.00 10 -0.2 -0.2 0.040 001 0.040 000 0.00 11 -0.3 -0.3 0.090 002 0.090 000 0.00 12 -0.4 -0.4 0.160 003 0.160 000 0.00 13 -0.5 -0.5 0.250 005 0.250 000 0.00 14 -0.6 -0.6 0.360 007 0.360 000 0.00 15 -0.7 -0.7 0.490 009 0.490 000 0.00 16 0.966 0.259 0.250 809 0.249 949 0.34 17 0.866 0.499 0.432 442 0.432 913 0.11 18 0.707 0.707 0.499 981 0.499 990 0.00 19 0.499 0.866 0.432 442 0.432 913 0.11 20 0.174 0.985 0.172 523 0.170 863 0.97 平均误差 0.08 2.3 非连续参数化单元与参数化单元比较
在二维情况下，分别采用非连续参数化单元和参数化单元边界元法解势流速度场问题，其中非连续参数化单元可以在较低的单元数下得到较为理想的求解精度。
由于非连续参数化单元上的物理量并非与单元节点处近似，故相比于参数化单元边界元法，非连续参数化单元在求解多边界问题时不但具有常数单元交界点无需特殊处理的优点，而且还具有线性及高次单元在较少单元数下得到较高求解精度的特点。
当采用传统常数或线性单元的直接边界元法解势流速度场问题时，除了边界近似带来的误差，还有在求矩阵系数时由高斯积分法近似引起的误差。非连续参数化单元由于是将单元上的点进行参数化，求解系数时使用解析表达式求解，即柯西主值积分，故其只有边界近似误差以及系统矩阵求解时带来的误差。
3. 三维势流问题算例
3.1 数学模型
考虑无穷远边界的流场[11-12]，有沿x轴负向的、速度为1的均匀来流。流体无旋、无粘性、不可压缩，流场中有一个半径为1的圆球。流场速度势由下式表示[6]：
{\mathit{\Phi }_0} = - {V_\infty }x + \varphi (26) 式中：Φ0为总速度势；V∞为无穷远处的来流速度；φ为圆球的扰动速度势。将求解的直接变量选择为扰动速度势φ，扰动速度势φ满足定解方程
\left\{ \begin{array}{l} {\nabla ^2}\varphi = 0,\mathit{\Omega }\\ q = \frac{{\partial \varphi }}{{\partial {\mathit{\boldsymbol{n}}_Q}}} = - V \cdot {n_x}, \ \ Q \in S \end{array} \right. (27) 式中：Ω为流场区域；S为物面边界；nQ为Q点的单位法向量。
将基本解ϕ(x, y)=- \frac{1}{{4\pi r}}代入三维边界积分方程，离散并建立方程组，求解方程组后，便可得单元上的未知量。
3.2 离散与数值求解
使用三角形常数单元对求解域表面进行离散，并使用一种对单元上的点进行参数化的方法对问题进行求解。边界积分方程离散后，其求解过程与二维非连续参数化单元类似，此处不再赘述。图 4为圆球划分三角网格效果图。
在这里，使用式(28)对单元上的点进行参数化：
\begin{array}{*{20}{c}} {\left( {x,y,z} \right) = \left( {{X^k}\left( {u,v} \right),{Y^k}\left( {u,v} \right),{Z^k}\left( {u,v} \right)} \right)}\\ {0 < u < 1 - v,\;\;\;\;0 < v < 1} \end{array} (28) 式中：Xk，Yk，Zk为第k个单元上角点坐标的参数化表示；u，v为变换坐标中的量。
若|nzk|≥1/ \sqrt 3 ，则有
\left\{ \begin{array}{l} {X^k}\left( {u,v} \right) = \left( {x_2^k - x_1^k} \right)u + \left( {x_3^k - x_1^k} \right)v + x_1^k\\ {Y^k}\left( {u,v} \right) = \left( {y_2^k - y_1^k} \right)u + \left( {y_3^k - y_1^k} \right)v + y_1^k\\ {Z^k}\left( {u,v} \right) = - {\left( {n_z^k} \right)^{ - 1}}\left[ {n_x^k\left( {{X^k}\left( {u,v} \right) - x_1^k} \right) + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {n_y^k\left( {{Y^k}\left( {u,v} \right) - y_1^k} \right)} \right] + z_1^k \end{array} \right. (29) 若|nzk|＜1/ \sqrt 3 且|nyk|≥1/ \sqrt 3 ，则有
\left\{ \begin{array}{l} {X^k}\left( {u,v} \right) = \left( {x_2^k - x_1^k} \right)u + \left( {x_3^k - x_1^k} \right)v + x_1^k\\ {Z^k}\left( {u,v} \right) = \left( {z_2^k - z_1^k} \right)u + \left( {z_3^k - z_1^k} \right)v + z_1^k\\ {Y^k}\left( {u,v} \right) = - {\left( {n_y^k} \right)^{ - 1}}\left[ {n_x^k\left( {{X^k}\left( {u,v} \right) - x_1^k} \right) + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {n_z^k\left( {{Z^k}\left( {u,v} \right) - z_1^k} \right)} \right] + y_1^k \end{array} \right. (30) 若|nzk|＜1/ \sqrt 3 且|nyk|＜1/ \sqrt 3 ，则有
\left\{ \begin{array}{l} {Y^k}\left( {u,v} \right) = \left( {y_2^k - y_1^k} \right)u + \left( {y_3^k - y_1^k} \right)v + y_1^k\\ {Z^k}\left( {u,v} \right) = \left( {z_2^k - z_1^k} \right)u + \left( {z_3^k - z_1^k} \right)v + z_1^k\\ {X^k}\left( {u,v} \right) = - {\left( {n_x^k} \right)^{ - 1}}\left[ {n_y^k\left( {{Y^k}\left( {u,v} \right) - y_1^k} \right) + } \right.\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {n_z^k\left( {{Z^k}\left( {u,v} \right) - z_1^k} \right)} \right] + x_1^k \end{array} \right. (31) 式中，(x1, y1, z1)，(x2, y2, z2)，(x3, y3, z3)分别为三角形元素上角点1，2，3的坐标。
因此，可将单元积分转化为面上的积分D1k(ξ, η, ζ)和D2k(ξ, η, ζ)。
\begin{array}{*{20}{c}} {D_1^k\left( {\xi ,\eta ,\zeta } \right) = }\\ {\int\limits_0^1 {\int\limits_0^{1 - v} {{\mathit{\Phi }_{3{\rm{D}}}}\left( {{X^k}\left( {u,v} \right),{Y^k}\left( {u,v} \right),{Z^k}\left( {u,v} \right);\xi ,\eta ,\zeta } \right)} } }\cdot\\ {{J^k}{\rm{d}}u{\rm{d}}v} \end{array} (32) \begin{array}{*{20}{c}} {D_2^k\left( {\xi ,\eta ,\zeta } \right) = }\\ {\int\limits_0^1 {\int\limits_0^{1 - v} {\frac{\partial }{{\partial \mathit{\boldsymbol{n}}}}\left[ {{\mathit{\Phi }_{3{\rm{D}}}}\left( {x,y,z;\xi ,\eta ,\zeta } \right)} \right]\left| {_{\left( {x,y,z} \right) = \left( {{X^k}\left( {u,v} \right),{Y^k}\left( {u,v} \right),{Z^k}\left( {u,v} \right)} \right)}} \right.} } \cdot }\\ {{J^k}{\rm{d}}u{\rm{d}}v} \end{array} (33) 式中：Φ3D为三维拉普拉斯方程的基本解；Jk为雅可比常数，由下式给出：
{J^k} = 2\sqrt {{\sigma ^k}\left( {{\sigma ^k} - {\alpha ^k}} \right)\left( {{\sigma ^k} - {\beta ^k}} \right)\left( {{\sigma ^k} - {\gamma ^k}} \right)} 其中，
\begin{array}{l} {\sigma ^k} = \frac{{{\alpha ^k} + {\beta ^k} + {\gamma ^k}}}{2}\\ {\alpha ^k} = \sqrt {{{\left( {x_1^k - x_2^k} \right)}^2} + {{\left( {y_1^k - y_2^k} \right)}^2} + {{\left( {z_1^k - z_2^k} \right)}^2}} \\ {\beta ^k} = \sqrt {{{\left( {x_2^k - x_3^k} \right)}^2} + {{\left( {y_2^k - y_3^k} \right)}^2} + {{\left( {z_2^k - z_3^k} \right)}^2}} \\ {\gamma ^k} = \sqrt {{{\left( {x_3^k - x_1^k} \right)}^2} + {{\left( {y_3^k - y_1^k} \right)}^2} + {{\left( {z_3^k - z_1^k} \right)}^2}} \end{array} (34) 式(32)和式(33)又可以写为式(35)及式(36)，其中的t为前面所述的u。
\begin{array}{*{20}{c}} {D_1^k\left( {\xi ,\eta ,\zeta } \right) = }\\ {\int\limits_0^1 {\int\limits_0^1 {{\mathit{\Phi }_{3{\rm{D}}}}\left[ {{X^k}\left( {t\left( {1 - v} \right),v} \right),{Y^k}\left( {t\left( {1 - v} \right),v} \right),{Z^k}\left( {t\left( {1 - v} \right),v} \right);} \right.} } }\\ {\left. {\xi ,\eta ,\zeta } \right] \times \left( {1 - v} \right){J^k}{\rm{d}}t{\rm{d}}v} \end{array} (35) \begin{array}{*{20}{c}} {D_2^k\left( {\xi ,\eta ,\zeta } \right) = }\\ {\int\limits_0^1 {\int\limits_0^1 {\frac{\partial }{{\partial \mathit{\boldsymbol{n}}}}\left[ {{\mathit{\Phi }_{3{\rm{D}}}}\left( {x,y,z;\xi ,\eta ,\zeta } \right)} \right]\left| {_{\left( {x,y,z} \right) = \left( {{X^k}\left( {t\left( {1 - v} \right),v} \right),{Y^k}\left( {t\left( {1 - v} \right),v} \right),{Z^k}\left( {t\left( {1 - v} \right),v} \right)} \right)}} \right.} } \cdot }\\ {\left( {1 - v} \right){J^k}{\rm{d}}t{\rm{d}}v} \end{array} (36) 这样，就可以使用高斯积分公式(37)对上式进行数值求解[9]。
\int\limits_0^1 {\int\limits_0^1 {f\left( {t,v} \right){\rm{d}}t{\rm{d}}v} } = \frac{1}{{16}}\sum\limits_{i = 1}^{16} {f\left( {{t_i},{v_i}} \right)} (37) 3.3 速度求解
求得物面控制点势函数后，使用文献[6]所采用的方法计算物面的速度分布，详细步骤参见文献[13]。
本文选取不同的单元数，分别采用上述方法进行了数值实验，以参数化单元高斯积分法计算系数矩阵的结果如图 5~图 10所示。图中，相对误差即为相对解析解[11-12]的误差。
![]()
图 5 4 512个网格速度势计算结果Figure 5. The calculation results of velocity potential of 4 512 grids![]()
图 8 9 112个网格速度势计算结果Figure 8. The calculation results of velocity potential of 9 112 grids当单元数为4 512时，速度势的最大误差为19.76%，平均误差为1.72%；速度的最大误差为29.78%，平均误差为5.18%。
基于网格的划分方法，球的两个极点附近的单元数较少，即便增加总网格数，误差也较其他位置节点处的突出。由于所使用高斯积分公式的数值积分精度不高，因此，三维问题下的计算结果不如二维问题下的好，可以考虑结合采用改进网格划分和采用非连续单元计算，也可以采用其他数值积分方法来提高计算精度。
速度幅值计算结果的精度是可以接受的，但该计算结果也未能精确反映解析解给出的速度变化趋势。由图可以看到，速度势和速度的平均相对误差都可以接受，但某些节点的误差偏大。这主要是由引入的高斯积分公式精度不高、数值计算误差所致；另外，不仅速度势函数计算本身存在误差，物面速度计算也存在误差，故可以看到速度势与速度的误差特性有一定的相关性。
当单元数为9 112个时，速度势最大误差为17.62%，平均误差为1.23%；速度最大误差为27.98%，平均误差为4.92%。
数值实验表明，相对于文献[6]的7点高斯积分法，本文方法的优点是计算同数量网格所需时间更少，速度更快；缺点是求解精度存在问题，平均精度尚在可接受范围内，但某些节点的计算误差偏大。通过文中给出的2种情况可以看出，在网格加密的情况下，本文方法的精度提高并不大，数值实验表明，在约4 000网格数时就已达到相当于10 000网格数时的精度，也就是说，求解精度在4 000网格数左右达到最优结果；在计算机可计算范围内，随着网格的不断增加，计算结果中会出现误差非常大的异常点，笔者认为，这主要是由引入高斯积分公式及参数化单元计算时进行的坐标变换引起单元节点与实际节点不一致所引起。
4. 结语
边界元法的计算量主要集中在影响系数矩阵的计算上，误差也大多源于此。本文对多种参数化单元方法进行了数值实验，分析了各种方法对最终结果误差的影响。
二维下的非连续参数化单元和参数化单元边界元法，由于单元上的节点参数化，矩阵系数计算均为解析式计算，即柯西主值积分，速度快、精度高，没有引起数值计算误差。其中，非连续参数化单元在求解多边界问题时不但兼顾了常数单元交界点无需特殊处理，而且兼顾了二阶及以上单元在较少单元数下得到较高求解精度的优点。
三维下的参数化单元边界元法，由于单元上的节点参数化，计算速度相比于其他方法稍快，但矩阵系数计算引入了高斯积分公式，引进了数值计算误差，在计算系数积分时进行了坐标变换，使得变换后求得的单元节点与实际节点存在偏差，故又引入了新的误差。数值实验表明，计算平均误差可以接受，但存在着某些节点的误差偏大的问题，易导致求解不够精确。可以通过改进网格划分和采用非连续单元计算，或其他可通过坐标转换进行解析计算的单元来计算。
本文的二维非连续参数化单元可以用于对波浪与浮体，或者潜体的作用机理进行研究，单元数少，求解精度高；三维参数化单元虽然求解精度差，但求解速度快，若进一步改进，如非连续单元在三维时，可提高运算精度。
本文所涉及的参数化单元边界元法的编程与传统边界元法的编程程序相比具有通用性，可以解决一类在势流框架下的工程及科学问题。
-
![]()
图 5 4 512个网格速度势计算结果
Figure 5. The calculation results of velocity potential of 4 512 grids
![]()
图 8 9 112个网格速度势计算结果
Figure 8. The calculation results of velocity potential of 9 112 grids
表 1 非连续参数化单元速度势计算结果(案例1)
Table 1 The calculation results of velocity potential by discontinuous parameterized elements(case 1)
节点 x y 数值速度势计算结果 解析速度势计算结果 相对误差/% 1 0.1 0.1 0.008 537 0.008 546 0.11 2 0.1 2.5 0.236 313 0.236 265 0.02 3 0.1 5.0 0.628 773 0.625 917 0.46 4 0.1 7.5 1.422 434 1.421 926 0.04 5 0.1 9.9 3.062 776 3.043 561 0.63 6 4.9 0.1 0.271 407 0.271 940 0.20 7 4.9 2.5 7.534 049 7.518 070 0.21 8 4.9 5.0 19.896 945 19.917 008 0.10 9 4.9 7.5 45.318 735 45.246 430 0.16 10 4.9 9.9 96.962 292 96.847 673 0.12 11 0.2 0.1 0.017 073 0.017 084 0.06 12 2.0 0.1 0.159 467 0.159 921 0.28 13 2.5 0.1 0.192 659 0.192 385 0.14 14 4.0 0.1 0.258 137 0.258 758 0.24 15 4.9 0.1 0.271 407 0.271 940 0.20 16 0.2 9.9 6.096 510 6.084 118 0.20 17 2.0 9.9 57.112 596 56.953 737 0.28 18 2.5 9.9 68.502 724 68.515 455 0.02 19 4.0 9.9 92.406 620 92.153 083 0.28 20 4.9 9.9 96.962 292 96.847 673 0.12 平均误差 0.19 
下载: 导出CSV
表 2 参数化单元速度势计算结果(案例1)
Table 2 The calculation results of velocity potential by parameterized elements(case 1)
节点 x y 数值速度势计算结果 解析速度势计算结果 相对误差/% 1 0.1 0.1 0.002 490 0.008 546 70.86 2 0.1 2.5 0.236 344 0.236 265 0.03 3 0.1 5.0 0.625 979 0.625 917 0.01 4 0.1 7.5 1.420 772 1.421 926 0.08 5 0.1 9.9 3.985 676 3.043 561 30.95 6 4.9 0.1 0.329 797 0.271 940 21.28 7 4.9 2.5 7.534 036 7.518 070 0.21 8 4.9 5.0 19.975 335 19.917 008 0.29 9 4.9 7.5 45.382 206 45.246 430 0.30 10 4.9 9.9 96.152 946 96.847 673 0.72 11 0.2 0.1 0.014 657 0.017 084 14.21 12 2.0 0.1 0.160 352 0.159 921 0.27 13 2.5 0.1 0.192 864 0.192 385 0.25 14 4.0 0.1 0.258 801 0.258 758 0.02 15 4.9 0.1 0.329 797 0.271 940 21.28 16 0.2 9.9 6.487 487 6.084 118 6.63 17 2.0 9.9 56.897 571 56.953 737 0.10 18 2.5 9.9 68.448 599 68.515 455 0.10 19 4.0 9.9 92.072 400 92.153 083 0.09 20 4.9 9.9 96.152 946 96.847 673 0.72 平均误差 8.42
下载: 导出CSV
表 3 非连续参数化单元速度势计算结果(案例2)
Table 3 The calculation results of velocity potential by discontinuous parameterized elements(case 2)
节点 x y 数值速度势计算结果 解析速度势计算结果 相对误差/% 1 0.7 0.7 0.490 070 0.490 000 0.01 2 0.6 0.6 0.360 050 0.360 000 0.01 3 0.5 0.5 0.250 035 0.250 000 0.01 4 0.4 0.4 0.160 022 0.160 000 0.01 5 0.3 0.3 0.090 013 0.090 000 0.01 6 0.2 0.2 0.040 006 0.040 000 0.01 7 0.1 0.1 0.010 001 0.010 000 0.01 8 0 0 0.000 000 0.000 000 0.00 9 -0.1 -0.1 0.010 001 0.010 000 0.01 10 -0.2 -0.2 0.040 006 0.040 000 0.01 11 -0.3 -0.3 0.090 013 0.090 000 0.01 12 -0.4 -0.4 0.160 022 0.160 000 0.01 13 -0.5 -0.5 0.250 035 0.250 000 0.01 14 -0.6 -0.6 0.360 050 0.360 000 0.01 15 -0.7 -0.7 0.490 070 0.490 000 0.01 16 0.966 0.259 0.249 980 0.249 949 0.01 17 0.866 0.499 0.432 967 0.432 913 0.01 18 0.707 0.707 0.500 085 0.499 990 0.02 19 0.499 0.866 0.432 967 0.432 913 0.01 20 0.174 0.985 0.170 891 0.170 863 0.02 平均误差 0.01
下载: 导出CSV
表 4 参数化单元速度势计算结果(案例2)
Table 4 The calculation results of velocity potential by parameterized elements(case 2)
节点 x y 数值速度势计算结果 解析速度势计算结果 相对误差/% 1 0.7 0.7 0.490 009 0.490 000 0.00 2 0.6 0.6 0.360 007 0.360 000 0.00 3 0.5 0.5 0.250 005 0.250 000 0.00 4 0.4 0.4 0.160 003 0.160 000 0.00 5 0.3 0.3 0.090 002 0.090 000 0.00 6 0.2 0.2 0.040 001 0.040 000 0.00 7 0.1 0.1 0.010 000 0.010 000 0.00 8 0 0 0.000 000 0.000 000 0.00 9 -0.1 -0.1 0.010 000 0.010 000 0.00 10 -0.2 -0.2 0.040 001 0.040 000 0.00 11 -0.3 -0.3 0.090 002 0.090 000 0.00 12 -0.4 -0.4 0.160 003 0.160 000 0.00 13 -0.5 -0.5 0.250 005 0.250 000 0.00 14 -0.6 -0.6 0.360 007 0.360 000 0.00 15 -0.7 -0.7 0.490 009 0.490 000 0.00 16 0.966 0.259 0.250 809 0.249 949 0.34 17 0.866 0.499 0.432 442 0.432 913 0.11 18 0.707 0.707 0.499 981 0.499 990 0.00 19 0.499 0.866 0.432 442 0.432 913 0.11 20 0.174 0.985 0.172 523 0.170 863 0.97 平均误差 0.08
下载: 导出CSV
-
[1] ZHENG S M, ZHANG Y L. Wave diffraction and radiation by multiple rectangular floaters[J]. Journal of Hydraulic Research, 2016, 54(1):102-115. doi: 10.1080/00221686.2015.1090492
[2] LIU Y, LI H J. Analysis of wave interaction with submerged perforated semi-circular breakwaters through multipole method[J]. Applied Ocean Research, 2012, 34:164-172. doi: 10.1016/j.apor.2011.08.003
[3] SANNASIRAJ S A, SUNDARAVADIVELU R, SUNDAR V. Diffraction-radiation of multiple floating structures in directional waves[J]. Ocean Engineering, 2001, 28(2):201-234. doi: 10.1016/S0029-8018(99)00066-9
[4] 刘荣, 郑永红, 游亚戈, 等.波浪入射角及地形对浮体水动力学特性影响的有限元分析[J].中国工程科学, 2005, 7(3):42-48. doi: 10.3969/j.issn.1009-1742.2005.03.007 LIU R, ZHENG Y H, YOU Y G, et al. Finite element analysis for effects of incident wave angle and seabed on hydrodynamic characteristics of buoy[J]. Engineering Science, 2005, 7(3):42-48(in Chinese) doi: 10.3969/j.issn.1009-1742.2005.03.007
[5] BREBBIA C A. The boundary element method for engineers[M]. London:Pentech Press, 1980.
[6] 李井煜, 卢晓平, 赵鹏伟.直接边界元法解势流速度场问题[J].中国舰船研究, 2015, 10(1):68-75. doi: 10.3969/j.issn.1673-3185.2015.01.010 LI J Y, LU X P, ZHAO P W. Direct boundary element method for the problem of potential flow velocity field[J]. Chinese Journal of Ship Research, 2015, 10(1):68-75(in Chinese). doi: 10.3969/j.issn.1673-3185.2015.01.010
[7] KYTHE P K. An introduction to boundary element methods[M]. Boca Raton:CRC Press, 1995.
[8] ANG K C. Introducing the boundary element method with MATLAB[J]. International Journal of Mathematical Education in Science and Technology, 2008, 39(4):505-519. doi: 10.1080/00207390701722676
[9] ANG W T. A beginner's course in boundary element methods[M]. Boca Raton, USA:Universal Publishers, 2007.
[10] PENG M J, CHENG Y M. A boundary element-free method(BEFM)for two-dimensional potential problems[J]. Engineering Analysis with Boundary Elements, 2009, 33(1):77-82. doi: 10.1016/j.enganabound.2008.03.005
[11] 余灵, 李干洛, 羊少刚.球尾渔船表面流场的数值计算[J].华南理工学报(自然科学版), 1995, 23(10):155-163. http://www.cnki.com.cn/Article/CJFDTOTAL-HNLG510.023.htm YU L, LI G L, YANG S G. Numerical calculation of flow fields over ship hulls of bulbous stern fishing boats[J]. Journal of South China University of Technology(Natural Science), 1995, 23(10):155-163(in Chinese). http://www.cnki.com.cn/Article/CJFDTOTAL-HNLG510.023.htm
[12] 余灵, 李干洛, 羊少刚.船体表面流场的理论计算与数值分析[J].广东造船, 1994(3):1-7. http://www.wanfangdata.com.cn/details/detail.do?_type=perio&id=QK199400086437 [13] 王献孚.计算船舶流体力学[M].上海:上海交通大学, 1992:244-245. -
期刊类型引用(2)
1. 曾樱. 增设掺气坎对溢洪道水力特性的影响分析. 地下水. 2023(06): 315-317 . 
百度学术
2. 陈臣. 基于CFD模型的和粗糙度和尺度效应对溢洪道水流结构的影响研究. 水科学与工程技术. 2022(03): 36-41 . 百度学术
其他类型引用(0)
-
其他相关附件
-
本文双语pdf
2018-5-77_en 点击下载(2078KB)
-
计量
- 文章访问数: 634
- HTML全文浏览量: 225
- PDF下载量: 146
- 被引次数: 2
