边界元法在固体与结构数值分析领域是有限元法的一种重要补充^[1-3]。边界元法在计算中利用了基本解，因而比有限元等区域解法有着降低维数的优势，且便于模拟复杂的几何模型。在实施过程中，边界单元的类型和离散方式对计算结果的精度有重要影响。例如，为了减小离散带来的误差，得到收敛的高精度边界元解，在局部应力集中或梯度变化剧烈的区域往往要划分更密的边界网格来准确反映该处物理场的变化^[4]。为此，很多学者研究和提出了不同的网格划分和细化方式，并根据相关误差准则来提高计算精度的方法^[5-7]。此类方法的目的在于用一定的误差准则来指导单元的细化，采用较少单元获得更高精度的解，提高计算效率，也被称为高精度边界元法^{[4, 8-9]}。然而，在复杂的模型中，采用某一准则对单元进行划分时，往往会忽略另外一个问题：部分单元积分的源点离单元过近，导致单元积分的精度不足，并进一步影响代数方程组求解的精度。这类问题被称为近奇异问题或边界层效应(boundary layer effect)^[10-13]。通常来说，近奇异积分主要存在于以下几种情况^[14]：

1) 薄壁结构、涂层结构^{[12, 15-18]}和裂纹问题^[19-20]。此类问题通常出现在不同边界相距很近的情况。若不对该处的单元进行细化和相关处理，则很容易出现由于近奇异性而导致的精度下降，但细化会导致计算量增加，采用坐标变换方式处理近奇异问题是一个有效手段。

2) 近边界点物理量的计算^[21]。力学问题通常比较关注边界上的物理量，因为最大应力通常发生在边界上。当内部节点非常接近边界时，内点应力计算则涉及近奇异问题。

3) 网格尺寸变化较大处^[22]。当采用一种离散准则去划分边界单元时，往往容易忽略由于单元长度相差较大而导致的近奇异问题。

有效地处理近奇异积分是保证边界元法获得良好数值结果的前提之一，也是国内外学者研究的一个重要问题。目前处理的方法主要有：单元细分法^[22-25]、坐标变换法^[26-32]、解析半解析算法^{[21, 33-35]}等。单元细分法主要是通过细化单元的方式来减小单元积分的近奇异性，例如Mustoe^[25]提出的自适应Gauss积分方法，根据误差上界公式细化单元，保证结果满足指定精度。坐标变换法则是通过变量代换来降低单元积分的近奇异性，主要有三次多项式变换^[26-27]、Sigmoidal变换^[28-30]、sinh坐标变换^[31-32]等。解析半解析算法是将近奇异积分转化为单元首末节点的解析函数，而弱近奇异积分采用非线性变换后用常规Gauss求积方法计算^[34]，强近奇异积分和超近奇异积分则采用复平面正则化方法处理，其中规则化部分采用常规Gauss数值积分计算，形式简单的奇异积分部分则进行解析求解^[35]。

对于近奇异积分，通常只关注上述第1种情况，即两个边界相近的情况，而忽视了复杂模型中其他可能存在的源点距被积单元很近而导致的近奇异问题，也就是第3种情况。针对这个问题，本文以位势问题为例，提出基于几何指标的方法对边界法中所有单元积分的精度进行监测，并对由于源点离单元过近而导致单元积分精度不足的情况利用修正近奇异积分精度的方法对其进行修正，从而保证单元积分和最终求解结果的精度。

1 位势问题边界元法的基本理论

以二维位势问题为例，其边界积分方程表达式为

$ \begin{array}{*{20}{c}} {c\left( p \right)\phi \left( p \right) = }\\ {\int_\mathit{\Gamma } {\left[ {{\phi ^{\rm{S}}}\left( {p;q} \right)\frac{{\partial \phi \left( q \right)}}{{\partial n\left( q \right)}} - \frac{{\partial {\phi ^{\rm{S}}}\left( {p;q} \right)}}{{\partial n\left( q \right)}}\phi \left( q \right)} \right]{\rm{d}}\mathit{\Gamma }\left( q \right)} }. \end{array} $

(1)

其中：ϕ^S(p; q)和$\frac{{\partial {\phi ^{\rm{S}}}(p;q)}}{{\partial n(q)}}$分别是二维位势问题的位势及其梯度的基本解，p为场点，q为源点，Γ为积分边界，c(p)为边界的与p点处几何特征有关的常数，n为节点所在单元的外法线方向，ϕ表示实际物理场的位势解。

将边界Γ离散为N_e个边界单元，引入形函数N_l^(m)(ξ)，则离散后的边界积分方程表示为

$ \begin{array}{*{20}{c}} {c\left( {{p_i}} \right)\phi \left( {{p_i}} \right) = }\\ {\sum\limits_{j = 1}^{{N_{\rm{e}}}} {\sum\limits_{l = 1}^m {\frac{{\partial \phi }}{{\partial n}}} } \left( {{q_{jl}}} \right)\int_{{\mathit{\Gamma }_j}} {N_l^{\left( m \right)}\left( \xi \right){\phi ^{\rm{S}}}\left( {{p_i};q\left( \xi \right)} \right){\rm{d}}\mathit{\Gamma }\left( \xi \right)} - }\\ {\sum\limits_{j = 1}^{{N_{\rm{e}}}} {\sum\limits_{l = 1}^m \phi } \left( {{q_{jl}}} \right)\int_{{\mathit{\Gamma }_j}} {N_l^{\left( m \right)}\left( \xi \right)\frac{{\partial {\phi ^{\rm{S}}}\left( {{p_i};q\left( \xi \right)} \right)}}{{\partial n\left( q \right)}}{\rm{d}}\mathit{\Gamma }\left( \xi \right)} .} \end{array} $

(2)

式中：m是单元节点数，ξ表示局部坐标。对单元Γ_j的积分可以表示为：

$ \begin{array}{*{20}{c}} {h_{{p_i}{q_j}}^l = \int_{{\mathit{\Gamma }_j}} {N_l^{\left( m \right)}} \left( \xi \right){\phi ^{\rm{S}}}\left( {{p_i};q\left( \xi \right)} \right){\rm{d}}\mathit{\Gamma }\left( \xi \right),}\\ {g_{{p_i}{q_j}}^l = \int_{{\mathit{\Gamma }_j}} {N_l^{\left( m \right)}\left( \xi \right)\frac{{\partial {\phi ^{\rm{S}}}\left( {{p_i};q\left( \xi \right)} \right)}}{{\partial n\left( q \right)}}{\rm{d}}\mathit{\Gamma }\left( \xi \right)} .} \end{array} $

(3)

单元积分中的奇异性和近奇异性是由ϕ^S(p_i; q(ξ))和$\frac{{\partial {\phi ^{\rm{S}}}\left( {{p_i};q(\xi )} \right)}}{{\partial n(q)}}$这两项引起的。

2 单元与源点几何关系的指标λ

本文采用几何指标λ反映单元积分精度。λ为源点到被积单元最短距离与单元长度的比值，

$ \lambda = \frac{s}{l}. $

(4)

其中：s为源点到单元的最短距离，l为单元长度。

单元长度l的表达为

$ l = \int {\rm{d}} l = \int_{ - 1}^1 {\sqrt {{{\left( {\frac{{{\rm{d}}{x_1}}}{{{\rm{d}}\xi }}} \right)}^2} + {{\left( {\frac{{{\rm{d}}{x_2}}}{{{\rm{d}}\xi }}} \right)}^2}} {\rm{d}}\xi } . $

(5)

对于源点到单元的最小距离s，不同情况要进行分别讨论，如图 1所示。

若源点在单元上投影的位置不在单元内(包括投影点在单元的延长线上)(图 1a、图 1b、图 1c)，s取源点到单元两端距离的较小者；若源点的投影在单元上(图 1d)，则s的表达式为

$ s\left( \xi \right) = \sqrt {{{\left( {\sum\limits_{i = 1}^k {{x_i}} {\varphi _i}\left( \xi \right) - {x_p}} \right)}^2} + {{\left( {\sum\limits_{i = 1}^k {{y_i}} {\varphi _i}\left( \xi \right) - {y_p}} \right)}^2}} . $

(6)

其中：φ_i(ξ)是单元的插值形函数，k=1，2，3分别表示常值单元、线性单元、二次单元。当且仅当ξ=ξ_p时，s(ξ)取到最小值，即源点到单元的最小距离。(x_i，y_i)和(x_p，y_p)分别为单元节点和源点的整体坐标。

3 指标λ与单元积分精度的关系

式(3)所示的单元积分的计算通常采用Gauss求积公式，

$ I = \int_a^b f (x){\rm{d}}x = \frac{{b - a}}{2}\sum\limits_{k = 1}^n {{w_k}} f\left( {{x_k}} \right) + {E_n}. $

其中是E_n计算误差，其表达式为

$ {E_n} = \frac{{{L^{2n + 1}}{{(n!)}^4}}}{{\left( {2n + 1} \right){{[(2n)!]}^3}}}{f^{(2n)}}(x). $

(7)

为了进一步分析单元积分的相对误差与λ的关系，下面以1/r型的单元积分$\int_{ - 1}^1 {\frac{1}{{\sqrt {4{\lambda ^2} + \left( {\xi - {\xi _p}} \right)} }}} {\rm{d}}\xi $为例进行考察。Gauss积分点数取为20个，考察ξ_p=0.0，0.5，1.0这3种情况，计算单元积分的相对误差与λ的关系，结果如图 2所示。

由图 2可以看出，单元积分的源点距离单元过近时会导致单元积分精度的下降，并且ξ_p=0.0时的误差最大。为了使λ能准确地反映单元积分的精度，本文利用单元积分的误差上界公式，来表示单元积分相对误差与λ之间的对应关系。误差上界公式有以下几种准则：

3.1 Lachat与Watson的准则

通过对式(7)的推导和近似，Mustoe^[25]得到式(8)的Gauss求积误差上界公式，

$ {E_n} < 2{\left( {\frac{l}{{4s}}} \right)^{2n}}\frac{{(2n + p - 1)!}}{{(2n)!(p - 1)!}} < e. $

(8)

其中e表示目标积分精度。

Gao和Davies^[22]在此基础上，进一步提出关于Gauss积分点数的实用公式，

$ n = \frac{{{p^*}\ln \left( {e/2} \right)}}{{2\ln [l/(4s)]}}. $

(9)

其中: ${p^*} = \sqrt {\frac{2}{3}p + \frac{2}{5}} $，p=1，2，3表示奇异性阶次。

3.2 Davies与Bu的准则

Bu与Davies^[23]在1995年提出了一个基于曲面积分数值实验的准则。该准则在积分精度e=10^-4时，关于Gauss积分点数的表达式如下：

$ n = {\left( {\frac{{8{l_i}}}{{3s}}} \right)^{\frac{3}{4}}} + 1,\;\;\;\;\left( {\frac{1}{r}型} \right) $

(10)

$ n = {\left( {\frac{{4{l_i}}}{{3s}}} \right)^{\frac{3}{4}}} + 1.\;\;\;\;\left( {\frac{1}{{{r^2}}}型} \right) $

(11)

Gao和Davies^[22]在式(10)和(11)的基础上通过大量数值实验给出一个近似公式，

$ n = {p^*}\left[ { - \ln \left( {e/2} \right)10} \right]\left[ {{{\left( {\frac{{8l}}{{3s}}} \right)}^{\frac{3}{4}}} + 1} \right]. $

(12)

该公式是对Davies与Bu准则的拓展，适用性更广。

3.3 误差上界准则

本文基于近奇异型积分的Gauss求积最大相对误差E_n与源点到单元的相对距离λ在对数坐标下的近似线性关系，提出了如下误差上界表达式

$ {E_n} = {{\rm{e}}^{[(0.025p - 1.05)n + 2]4\lambda + 1.3p - 2}} $

(13)

及对应的Gauss积分点n的表达式，

$ n = \frac{{2 - \frac{{\ln {E_n} - 1.3p + 2}}{{4\lambda }}}}{{1.05 - 0.025p}}. $

(14)

3.4 3种上界公式的对比

为了能让λ更好地反映单元积分精度，本文进一步考察了3种误差上界公式估算结果与实际误差的关系。通过选取20个Gauss积分点，以核函数1/r型为例，比较了在不同λ下的结果，如图 3所示。

从图 3可以看出，Lachat与Watson准则对计算结果精度的估算是非常消极的，在λ＜1的区间中，估算结果总是比实际结果大3、4个量级以上，且随着λ的变小，偏差更大，不利于对单元积分精度的估算。Bu与Davies准则并不能确保作为所有数据点的误差上界，因为该准则并未考虑到Gauss积分中精度最差的情况，导致在比较大的范围内，该准则会失效。因此，本文采用所提出的式(13)来作为几何指标与单元积分对应关系的表达式。

4 全局求解精度指标

本文进一步考查单元积分的误差对最终求解结果的影响。利用矩阵误差传递公式得到计算结果的最大可能相对误差，作为一个度量全局单元积分精度的指标E。

对于代数方程组Ax=b，在考虑误差的情况下，其实际计算的结果可以表示为

$ \left( {\mathit{\boldsymbol{A}} + \delta \mathit{\boldsymbol{A}}} \right)\mathit{\boldsymbol{x'}} = \mathit{\boldsymbol{b}} + {\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{b}}. $

(15)

其中：δA和δb分别表示相应的偏移量，x′表示方程组的数值解。

方程组发生扰动时解的误差的估计式为^[36]

$ \begin{array}{*{20}{c}} {E = \frac{{{{\left\| {\mathit{\boldsymbol{x}} - {\mathit{\boldsymbol{x}}^\prime }} \right\|}_\alpha }}}{{{{\left\| \mathit{\boldsymbol{x}} \right\|}_\alpha }}} \le }\\ {\frac{{ {\rm{cond}} (\mathit{\boldsymbol{A}})}}{{1 - {\rm{cond}} (\mathit{\boldsymbol{A}})\left( {\left\| {{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{A}}} \right\|/\left\| \mathit{\boldsymbol{A}} \right\|} \right)}}\left( {\frac{{\left\| {{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{A}}} \right\|}}{{\left\| \mathit{\boldsymbol{A}} \right\|}} + \frac{{{{\left\| {{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{b}}} \right\|}_\alpha }}}{{{{\left\| \mathit{\boldsymbol{b}} \right\|}_\alpha }}}} \right).} \end{array} $

(16)

其中：cond(A)=‖A‖_α‖A^-1‖_α表示矩阵A的条件数，‖·‖_α表示矢量的α范数。

在边界元法的实施过程中，δA和δb是指数值积分结果和精确解的差值，可用于估算方程组解的最大可能相对误差。在实施求解之前，通过式(13)和几何指标λ估算误差值δA和δb，并通过式(16)求得E，本文称为估算全局指标；求解后用实际误差值求得的式(16)结果称为实际全局指标，记为E′。据此，本文将E(或E′)作为一个全局的求解精度指标。

5 单元积分精度修正方法

对于局部的单元积分因精度不足而影响了最后结果的求解精度情况，本文利用sinh坐标变换的方法进行精度修正^[32]。

sinh坐标变换方法是利用坐标变换将原本具有近奇异的积分在变换后降低其奇异性，从而可以采用常规数值积分而确保积分精度。例如，对于原积分表达式

$ g = \int_{ - 1}^1 f (\xi )\ln r(\xi ){\rm{d}}\xi , $

(17)

现用一个新变量η去代替原来的变量ξ,

$ \xi = a + s\sinh \left( {{k_1}\eta + {k_2}} \right). $

(18)

其中：a是源点到单元投影点的局部坐标，s是源点到单元的最短距离。k₁和k₂的取值使得积分区域保持为[-1, 1], 便于直接使用常规Gauss求积公式：

$ {k_1} = \frac{1}{2}\left( {{\rm{arcsinh}}\left( {\frac{{1 + a}}{s}} \right) + {\rm{arcsinh}}\left( {\frac{{1 - a}}{s}} \right)} \right), $

(19)

$ {k_2} = \frac{1}{2}\left( {{\rm{arcsinh}}\left( {\frac{{1 + a}}{s}} \right) - {\rm{arcsinh}}\left( {\frac{{1 - a}}{s}} \right)} \right). $

(20)

于是可得采用新变量的表达式为

$ g = \int_{ - 1}^1 f \left( {\xi \left( \eta \right)} \right)\ln r\left( {\xi \left( \eta \right)} \right)J\left( \eta \right){\rm{d}}\eta . $

(21)

其中

$ J(\eta ) = \frac{{{\rm{d}}\xi }}{{{\rm{d}}\eta }} = s{k_1}\cosh \left( {\eta {k_1} + {k_2}} \right) $

(22)

是变换的Jacobi值。采用8个Gauss积分点，核函数取1/r、形函数取(1+ξ)/2的数值模拟结果如图 4所示。

由图 4可以看出，单元积分在坐标变换后的相对误差明显降低，在λ比较小时，提高精度的效果尤其明显，说明采用sinh坐标变换的方式可以提高λ较小时的单元积分精度，并且具有不改变既有计算流程和计算成本基本不增加的优点。

6 数值算例

为了检验本文提出的检验单元积分并进行修正方法的有效性，对以下两种可能出现近奇异性的问题进行计算、检验和修正。

算例1  薄壁问题。如图 5所示的矩形区域，长为L=100，宽为h，逐渐减小h值使从0.1~1变化。在矩形的上下边施加流场q=$ \mp $ 1，左右边施加势场u=x。L上划分10个单元，h上划分1个单元。计算结果如表 1所示。

由表 1可知，在薄壁问题中，当h取值越来越小时，其全局的精度指标E将会变得越来越差。当E过小时，直接采用sinh坐标变换对其进行修正，可以使计算结果的精度得到很好的提升，且该过程实施流程简单，额外计算成本可以忽略。

算例2  网格尺寸相差较大的问题。如图 6所示，势场为$u = \ln \;\frac{1}{r}$的矩形区域，边界条件为${\left. u \right|_{x = 0, \;L}} = \ln \;\frac{1}{r}, \;{\left. q \right|_{y = 0, \;L}} = - \frac{1}{r}\frac{{{\rm{d}}r}}{{{\rm{d}}y}}$。

采用以下两种方式(单元总数均为40个线性单元)对边界进行划分：

1) 对边界上所有单元进行等长度划分，如图 7a所示；

2) 根据该势函数梯度场的变化实施离散，如图 7b所示。

以图 6的B点处的流场值的相对误差作为结果进行比较(Gauss积分点数取为4)，结果见表 2。

由表 2可以看出，采用相同数量的单元，等梯度划分的计算结果要好于等长度划分。等长度划分方式对应的全局指标非常小，可以认为其单元积分的精度没有对最终结果造成影响，目前的结果体现了40个单元情况下这一划分方式最好的求解精度；而等梯度划分方式尽管求解精度高于等长度划分，但它的指标E比等长度划分大，影响了计算结果精度，进一步分析可以发现这是由图 8所示的近奇异问题所导致的。

采用本文提出的近奇异修正方法对等梯度划分中精度过低的单元进行修正，结果见表 3。

可以发现，修正后全局指标得到了改善，提高了计算结果的精度，且该过程实施流程简单，相比于增加Gauss积分点的方式，能有效减少计算量。

再对等梯度网格划分进行加密，计算结果见表 4。

由表 4可以看出，加密后全局指标(E或E′)仍会改变，在本例中变大了，反映了单元积分精度变差。对单元积分进行修正后，全局指标和最终计算结果均得到有效改善。

可见，如果在边界元法的实施过程中，只加密单元而不关注本文提出的全局指标E，计算结果的精度可能无法得到有效提升，甚至可能会更差，而考虑全局指标，通过本文提出的精度修正方法恢复单元积分精度后，能有效解决这个问题。修正后的结果随着单元加密，精度不断提高，说明此时已不存在单元积分精度不足的问题。

7 结论

本文采用了几何指标λ反映单元积分精度，进一步建立了全局精度指标E，并根据E的大小判断是否需要对单元积分进行修正。本文提出采用sinh坐标变换法修正和恢复单元积分精度。本文所提出的基于几何指标和修正边界元法求解精度的方法，无需改变原计算流程，增加的计算成本可以忽略，实施简单方便。数值算例证明了该方法可行有效，能在既有网格离散的基础上保证边界元法整体刚度矩阵的准确性和最终求解精度。本文的结果表明：

1) 对于全局精度指标E过小的划分方式，采用sinh坐标变换法进行修正，能够有效提高计算结果的精度；

2) 网格加密后，全局指标和单元精度也会改变，必须在每次加密后都监测并在必要情况下进行单元精度的修正；

3) 在边界元法实施过程中，若不关注全局指标而盲目地加密单元，计算结果的精度可能无法得到有效提升，甚至可能会更差。