AN ADAPTIVE ELEMENT SUBDIVISION METHOD FOR SINGULAR INTEGRALS BASED ON AFFINE TRANSFORMATIONS
-
摘要:
边界元法已广泛应用于解决工程实际问题中，精确高效地计算奇异积分对于求解边界积分方程至关重要。为此，提出了一种基于仿射变换的奇异积分自适应单元细分法，用于解决边界积分方程中的奇异性问题。其基本思想是通过仿射变换对单元进行特征分区，并结合自适应二叉树细分技术生成高质量的积分单元。利用仿射变换和特征分区技术将初始单元划分为细分区域和投影区域，不同区域的单元细分算法分别执行、计算效率高。源点附近的细分子块采用Serendipity积分单元构建，其余积分单元则采用传统积分方法进行计算。与传统单元细分方法相比，该方法具有自适应细分、积分精度高、易于实现等优点。数值算例验证了该文方法具有良好的准确性、鲁棒性和可行性。
Abstract:Boundary element method (BEM) has been widely used for solving the practical engineering problems. Accurate and efficient evaluation of singular integral is of crucial importance for solving the boundary integral equations. In the BEM implementation, an adaptive element subdivision method for singular integrals based on affine transformations is presented for evaluating the singular integrals. The basic idea consists of partitioning the input surface element via affine transformation and then generating a set of high-quality patches by adaptive binary-tree subdivision. By using the domain partitioning technique, the surface element can be divided into several element projection and subdivision regions under affine transformations. It is far more efficient to separately perform the element subdivision for different regions where the desirable patches are required. The ultimate sub-elements in the vicinity of the singular point are constructed by the serendipity patches, while the remaining patches are evaluated accurately by the conventional quadrature techniques. The proposed method has some advantages over the conventional element subdivision methods, such as the adaptive element subdivision, the improved accuracy and the straight-forward implementation. Numerical results are provided to validate the accuracy, robustness and availability of the proposed method.
-
边界元法[1-2]是一种将区域中微分方程边值问题归化为边界上的积分方程[3-4]并在边界上数值求解该积分方程的数值计算方法，以高精度、降低问题维度等优势在数值计算中广泛地应用。在边界积分方程求解计算中，单元的积分方式和单元的几何形状对数值计算精度有重要影响，例如，在物理变量变化梯度较大处需要进行单元加密[5]；薄型结构等问题[6]不同的边界之间距离很近。此类情况均会带来源点距离单元边界过近导致产生的积分块质量较差的问题(如邻边角较大或狭长的积分块)，即奇异问题。奇异积分[7-8]是指源点在单元内部或边界上的积分，计算奇异积分时，随着源点到场点的距离变小，被积函数的值变换越剧烈，此时常规的积分方法无法保证积分精度[9-10]。当求解问题规模较大、几何边界形状复杂时，无法准确计算质量较差的计算单元[11]，从而影响最终计算结果的精度。因此，精确高效地计算奇异积分对边界元法在实际工程应用中具有重要意义。
为了提高边界元法的计算性能，近年来人们提出并发展了各种积分方法，以消除边界积分方程中积分的奇异性。主要包括以下几类：① 解析和半解析方法。牛忠荣等[12-14]提出的解析和半解析方法，该方法成功的耦合到平面单元上并取得了很好的数值计算结果。② 非线性变换法。TELLES[15]提出的二次多项式和三次多项式进行变换可以消除积分的奇异性。③ 距离变换方法。马杭等 [16]提出的距离变换方法可以应用到各阶次的奇异积分，应用范围较广。
针对奇异积分的计算方法，提出了各种单元细分技术，如传统的细分方法、球面细分方法[17]、二叉树/四叉树细分方法[18]等。以上方法的核心思想是通过将单元细分成一组子单元，再对含奇异点的子单元进行数值计算。传统的细分方法对质量不好的单元进行细分时，细分后得到的积分子单元质量较差。二叉树/四叉树细分方法的通用性强稳定性好，但在源点周围存在冗杂的子单元使计算效率降低。球面细分法与结构网格生成的前沿推进法有些类似，效率好精度高，但该方法引入网格生成的模板法，在自适应单元细分过程中采用的模板数量较多，在处理三维体单元的奇异积分计算时情况较为复杂。
结合上述方法的特点，本文提出了一种基于仿射变换的奇异积分自适应单元细分法，用于计算边界积分，与现有的细分方法相比，该方法无需重构投影腔面或使用模板对单元进行细分。该细分方法可对任意类型单元进行仿射变换，将单元划分为细分区域和投影区域，依据细分准则在细分区域内对子单元递归细分，采用投影算法填充源点和投影区域边界的空隙，最终得到高质量的积分子单元。
1 奇异积分
考虑位势问题，区域内源点的边界积分方程可表示为：
c(P)u(P)=∫Γu∗(P,Q)q(Q)dΓ(Q)−∫Γq∗(P,Q)u(Q)dΓ(Q) (1) 式中：P为源点；Q为场点，Q∈Γ；c为系数，其值与源点P所在边界Γ的几何形状有关；u*、q*为基本解，各基本解的表达式如下：
u∗=14πr (2) {q^ * } = \frac{{\partial {u^ * }}}{{\partial {{\boldsymbol{n}}}}} (3) 式中： r = \left| {P - Q} \right| 为源点P到场点Q的距离；n为单位法向量。
式(1)中只涉及边界积分，在边界离散后，主要计算以下两种形式的积分：
G = \int_\varGamma ^{} {{u^ * }} (P{\text{,}}Q)N(Q){\text{d}}\varGamma {\text{(}}Q{\text{)}} (4) H = \int_\varGamma ^{} {{q^ * }} (P{\text{,}}Q)N(Q){\text{d}}\varGamma {\text{(}}Q{\text{)}} (5) 式中，N为单元形函数。
当源点P位于积分单元内部或边界上时，式(4)和式(5)中的积分为奇异积分，其基本解具有奇异性。源点距离场点越近，被积函数的数值变化越剧烈，使用常规的高斯积分准则会失效。本文以三角形单元和四边形单元为例，研究基于仿射变换的奇异积分有效算法。
2 基于仿射变换的单元细分
2.1 仿射变换
仿射变换是一种在空间坐标系中的线性转换，可以保持单元的平行性、平直性和共线性，主要涉及旋转、缩放、剪切、平移等几何变换，如图1所示。一个矩阵表示单元进行一次仿射变换，一组矩阵表示单元进行多种仿射变换，本文采用仿射变换中的缩放、平移等几何变换方法实现对三角形、四边形单元的自适应细分。仿射变换定义为：
{{{\boldsymbol{V}}}^*} = {{\boldsymbol{V}}} \times {{\boldsymbol{D}}} \times {{\boldsymbol{S}}} {{{\boldsymbol{V}}}^*} = \left[ {\begin{array}{*{20}{c}} x \\ y \\ {\textit{z}} \\ 1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} 1&0&0&{{T_x}} \\ 0&1&0&{{T_y}} \\ 0&0&1&{{T_{\textit{z}}}} \\ 0&0&0&1 \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{S_x}}&0&0&0 \\ 0&{{S_y}}&0&0 \\ 0&0&{{S_{\textit{z}}}}&1 \\ 0&0&0&1 \end{array}} \right] (6) 式中：D和S为平移矩阵和缩放矩阵；Sx、Sy、Sz为缩放比例因子；Tx、Ty、Tz为平移向量；V为未进行仿射变换的单元齐次坐标；V*为仿射变换后得到的齐次坐标。
2.2 基于仿射变换的自适应单元细分
本文提出的基于仿射变换的自适应单元细分方法，不同于传统的单元细分，是通过仿射变换对初始单元进行区域分割。其基本思想为通过仿射变换对面单元的边界进行缩放、平移得到投影区域和细分区域，再分别对两个区域进行投影及自适应细分，最终得到高质量的子单元，提高积分精度和计算效率。基于仿射变换的自适应单元细分法的主要步骤如下：
Step 1：单元预处理。由于源点相对于单元的位置不同，直接对单元进行分区较为困难，将源点与初始单元的顶点相连构建分区后，通过确定单元中心点和源点的相对位置关系来对单元的分区进行质量评判。
Step 2：构造单元投影区域。采用缩放变换按照一定缩放比例构建包围盒，再采用平移变换将包围盒沿指定方向移动到包含源点的位置，这有利于提高积分的计算效率。
Step 3：构建单元细分区域。构建完成投影区域后，通过将投影区域的顶点与单元顶点相连接来构建单元细分区域。依据单元细分准则以及自适应细分技术对细分区域进行细分。
Step 4：采用投影算法构建源点附近的积分子单元。为了消除积分的奇异性，要求单元投影区域的最终子单元需通过奇异点，采用径向投影算法来填补单元投影区域边界与源点之间的空腔。
由于仿射变换的缩放比率定义为α=0.3~0.6，对单元边界沿径向缩放，所以在离源点越近的地方细分子单元分布越密集，离源点越远的地方子单元分布越稀疏。此外，采用α-β变换子单元进行精确计算来消除积分的奇异性，其余单元采用常规的高斯积分进行计算，这样的单元细分既保证了积分精度又保证了积分效率[19-20]。
图2显示了一个线性四边形单元细分的流程，图3为基于仿射变换的自适应单元细分流程图。
![]()
图 2 线性四边形单元细分流程Figure 2. The element subdivision process of a linear quadrilateral element![]()
图 3 基于仿射变换的自适应单元细分流程图Figure 3. Flow diagram of adaptive element subdivision based on affine transformation2.3 初始单元的预处理
由于源点相对于单元的位置不同，若直接对单元进行区域划分，导致最终生成的子单元质量差，子单元分布不均匀。为了保证能够对单元进行合理的分区，需要在单元细分之前确定源点在单元内的位置，并对分区的质量进行评判。如图4所示，将源点与单元顶点相连后通过下列条件对分区进行质量评判。
![]()
图 4 初始单元预处理及分区质量判断Figure 4. The initial element preprocessing and domain partition quality assessment1) 纵横比。纵横比是指源点到单元边界的垂直距离与单元边界两顶点的线段长度的比值。一个合格的区域需满足条件Fref <F，其中F为纵横比，Fref为参考比率，通常取值为0.1。如图4所示，F=H/LAB，H为源点到分区边界的垂直距离，LAB为分区边界的两顶点之间的距离。
2) 扁平率。扁平率是指源点到单元边界的垂直距离与中心点到单元边界垂直距离之比。如图4所示，一个合格的区域需满足ηref <η，其中η=H/h。ηref为扁平率的参考值，通常取值为0.3。
3) 邻边角。邻边角是指分区中以源点O为中心，以源点和顶点相连的两线段为边组成的夹角。图4中∠AOB为邻边角，对于质量较好的细分区域，通常邻边角约为60°，当分区中的邻边角为钝角时，需对其进一步分割为锐角，保证积分精度。
4) 边长比。边长比是指分区中源点与顶点相连的短边与长边之比。合格的分区应满足条件Uref <U。其中，Uref为参考比率，通常取值为0.3；U为边长比，如图4所示，U=LOA/LOB，LOB为源点到顶点B的距离，LOA为源点到顶点A的距离。
本文通过定义缩放比对初始单元进行缩放，缩放比是初始单元的边界长度与仿射变换后包围盒对应边界的长度之比。缩放比α的值与源点到单元顶点之间的距离有关，即α=H/d，参数d表示源点到单元顶点的平均距离，缩放比通常取0.3~0.6。
2.4 考虑源点及单元形状的特征分区
本文根据源点的初始位置，采用仿射变换中的缩放、平移对初始单元进行自适应细分，构建单元细分区域和单元投影区域。如图5所示，以线性三角形单元为例，根据源点与单元的相对位置关系，通常有3种基本类型的分区方案：第一种类型是源点位于单元顶点附近，这种类型通常将单元划分为1个投影区域和1个细分区域；第二种类型是源点邻近三角形单元的边界，这种类型通常将三角形单元划分为1个投影区域和2个细分区域；第三种类型是源点位于单元内部，这种类型通常将单元划分为1个投影区域和3个细分区域。
![]()
图 5 不同源点位置单元细分方法示意图Figure 5. Schematic diagram of element partition for different locations of source points如图5所示，以线性四边形单元为例，第一种类型是源点位于单元顶点附近，这种类型通常将单元划分为1个投影区域和2个细分区域；第二种类型是源点邻近四边形单元的边界，这种类型通常将四边形单元划分为1个投影区域和3个细分区域；第三种类型源点位于单元内部，这种类型通常将单元划分为1个投影区域和4个细分区域。
基于单元形状和源点位置的分区方案不是唯一的，可以根据单元细分的具体要求来进行分区，具体步骤在后文单元细分区域的构建中详细说明。
2.5 单元投影区域的构建
采用仿射变换在单元的源点处对初始单元进行缩放，再根据源点与初始单元的相对位置关系，计算出包围盒内参考点的坐标，计算参考点到源点之间的平移向量，使用平移变换将包围盒沿着指定方向移动到源点所在的位置，得到包含源点的投影区域。
构建完成投影区域后，通过源点和投影区域中心点的坐标位置，计算出投影区域内两点的扁平率，判断构建的投影区域质量是否满足投影填充的条件，如果满足则采用投影算法填充投影区域和源点外球面之间的空隙，如不满足则继续对其进行仿射变换或二叉树细分处理。如图6所示，源点靠近投影区域的边界或顶点时，可对投影区域进行二叉树细分，将源点在投影区域的范围进一步缩小，使重构的细分子单元更好的贴近球面。该细分有利于积分点在单元内部更合理的分布，从而提高单元积分的计算效率。
图7描述了单元投影区域构建的示意图。由于仿射变换的共线性，使得单元投影区域与初始单元的相对边相互平行，这有利于提高积分子单元生成的质量。
![]()
图 7 单元投影区域构建示意图Figure 7. Schematic diagram of the element projection region construction2.6 单元细分区域的构建
通过对初始单元进行预处理，确定每个分区均为规则的细分区域后，通过仿射变换来构造细分区域。单元细分区域构建算法的主要步骤如下：
Step 1：确定源点和初始单元的相对位置关系，通过仿射变换构造投影区域。
Step 2：对于源点位置规则的初始单元，将单元的顶点与投影区域的顶点相连接，按照基本的分区方案构造初始单元的细分区域。
Step 3：通过计算每个细分区域的单元细分参数来评判每个细分区域的质量，若分区质量不好，则进行标记并做退化处理。
在构建细分区域过程中，单元的类型以及源点在单元内的位置关系影响细分区域构建的质量，对于源点位置不规则的初始单元，直接连接投影区域与单元顶点划分的细分区域会存在低质量的区域，如图8所示。源点附近的最小线框表示投影区域，单元中的其他空白部分表示细分区域，其中阴影部分为低质量细分区域。以线性四边形单元和线性三角形单元为例，导致构建的细分区域质量较低的情况有以下几种：
1) 对于线性四边形单元，当源点位于单元的内部且靠近单元顶点时。如图8(a)所示，直接连接投影区域与单元顶点时，细分区域存在两处低质量区域，需要将此类型退化为源点靠近单元顶点的类型。
2) 对于线性四边形单元，当源点位于单元的内部且靠近单元边界时。如图8(b)所示，直接连接投影区域与单元顶点时，细分区域存在一处低质量区域，需要将此类型退化为源点靠近单元边界的类型。
3) 对于线性四边形单元，当源点位于单元的内部且靠近单元顶点以及单元边界时。如图8(c)所示，直接连接投影区域与单元顶点时，细分区域存在一处低质量区域，需要将此类型退化为源点靠近单元顶点的类型。
4) 对于线性三角形单元，当源点位于单元内部且靠近单元边界以及单元顶点时。如图8(d)所示，直接连接投影区域与单元顶点时，细分区域存在两处低质量区域，需要将投影区域的边界延长至线性三角形单元的边界，将其划分为一个高质量细分区域。
5) 对于线性三角形单元，当源点位于单元内部同时靠近单元边界时。如图8(e)所示，直接连接投影区域与单元顶点时，细分区域存在一处低质量区域，需要将投影区域边界延长至线性三角形单元的边界，将其划分为2个高质量细分区域。
![]()
图 8 基于仿射变换细分方法对低质量细分区域的处理Figure 8. The low quality subdivision region processing based on affine transformation subdivision method3 仿射变换的单元细分方案
3.1 仿射变换的单元细分准则
结合区域划分技术，通过仿射变换可以将初始单元划分为单元投影区域和单元细分区域，这使单元内不同的区域进行积分更加高效。初始单元在细分后，应该保证积分子单元在单元内部分布合理，使积分子单元在源点处呈现近密远疏的分布状态，本节将详细说明在细分区域内用于单元细分的细分准则。
在几何计算方面，与其他的空间分解法相比较，利用仿射变换进行的细分方法具有简单、鲁棒性好的优点，可以更高效地实现单元的自适应细分。细分准则示意图如图9所示，为了更高效地计算奇异积分，单元细分准则可以表示为d<Rmax/2i，其中几何参数d为源点到中心点的距离，d1为源点到单元边界的最大垂直距离，d2为源点到单元顶点的最大距离，因此Rmax表示源点到单元边界的最大距离。依照细分准则检验每一个子单元是否满足细分要求，如果子单元不满足细分准则，则继续对子单元进行递归细分，直到满足细分准则停止细分。
3.2 基于仿射变换的单元细分
由于非结构网格的单元形状较为复杂[21]，采用单一的细分技术直接对其细分比较困难。为了可以准确地获取单元细分区域的几何特征，采用自适应生成技术来构建初始细分结构。这样可以保证单元细分区域的积分子单元分布合理。实现自适应单元细分算法的主要步骤：
1) 计算源点到单元边界的最大距离和最小距离；
2) 依据细分准则，对所有的单元细分区域进行细分；
3) 为了得到高质量的积分子单元，在细分时优先沿径向细分；
4) 确定细分的子单元是否满足细分要求，并将形状较小的子单元与相邻的子单元进行合并，形成质量较好的子单元。
4 源点附近的积分子单元构建
利用仿射变换对单元细分，每个单元投影区域都可以得到所需的细分结构。为了消除积分的奇异性，单元投影区域中生成的积分子单元需要通过源点。本节通过在投影区域构建投影多边形以及采用径向投影算法填充了单元投影区域边界与源点之间的空白区域。
4.1 构建投影多边形
由于单元类型以及源点位置不同，使初始单元划分的细分区域和投影区域类型也不同。通常情况下，若源点与投影区域的边界之间的距离较大，直接将源点与投影区域的顶点相连会使最终产生的积分子单元的质量较低，所以为保证积分的准确性，需要单元投影区域的边界尽可能的靠近源点。
当源点位于投影区域的中心，直接对其进行投影便可到得高质量的积分子单元。当源点位于单元投影区域的顶点或边界附近，需要先构建几何投影多边形进行优化处理，再对其进行投影。利用仿射变换在源点周围将投影区域划分为若干个子单元，投影多边形是由与源点相邻的子单元的边界构成，为满足奇异积分精度的要求，通过将投影多边形的顶点移动至单元边界的优化技术来提高积分子单元生成的质量。
4.2 径向投影算法
为了更好的计算奇异积分，采用径向投影算法重新填充源点和投影区域之间的细分子域。径向投影算法的核心思想是将投影区域沿球的半径方向逐层插入三角形、四边形的单元来填充单元投影区域边界和球面之间的空隙[22]。球心和球径的参数由源点和参考半径决定，其中在空隙中插入的三角形和四边形单元又称为缓冲层，可以提高积分子单元生成的质量。为了保证源点周围的子单元和球面相连，采用三角形单元和四边形单元作为最终的子单元来进行空隙填充。
5 数值算例
本节将给出三角形单元和四边形单元奇异积分的数值算例，并将本文方法与传统积分方法进行比较，验证了本文方法的准确性和收敛性。为了验证本文方法的准确性，定义以下相对误差的计算公式：
{I_{\rm error}}{\text{ = }}\left| {\frac{{{I_{\rm n}} - {I_{\rm e}}}}{{{I_{\rm e}}}}} \right| (7) 式中，In和Ie分别为奇异积分的近似解和参考解。
参考解是通过对积分单元进行较大数量的子单元细分，然后在子单元上使用大量的高斯积分点计算得到，本文用于奇异积分计算的表达式如下：
I = \int_{\varGamma }^{} {\frac{1}{{4\pi r}}N{\text{d}}{\varGamma }} (8) 式中：N为给定单元的形函数；Γ为积分域。
为了准确计算奇异积分，对源点周围的三角形单元和四边形单元进行数值积分，其余规则的区域采用标准的高斯积分以及Hammer积分来计算，积分点的数目m由式(9)计算：
m = - \frac{1}{{10}}\ln \left( {\frac{e}{2}} \right)\sqrt {\frac{2}{3} + \frac{2}{5}} \left[ {{{\left( {\frac{{8L}}{{3R}}} \right)}^{\frac{3}{4}}} + 1} \right] (9) 式中：m为积分点数；e为积分相对误差；L为单元积分方向最大长度；R为源点到几何边界的距离。采用本文方法的计算结果与传统方法的数值结果进行对比，分别给出各自的积分点数目以及相对误差大小。
5.1 线性三角形单元
此算例展示了线性三角形单元的计算结果，三角形单元的节点坐标分别为(0.0, 1.0, 0.0)、(8.0, 0.0, 0.0)、(8.0, 2.0, 0.0)。源点位置基本上沿图10所示的直线均匀分布，源点编号在全局坐标系中对应的坐标从左至右依次为(3.1, 0.9, 0.0)、(4.1, 0.8, 0.0)、(5.1, 0.7, 0.0)、(6.1, 0.6, 0.0)。
![]()
图 10 线性三角形单元节点及源点坐标Figure 10. Coordinates of source points and nodes of the linear triangular element图11为在不同源点位置，线性三角形单元在基于仿射变换下的自适应单元细分法细分后得到的子单元。从图11中可以看出，用该方法得到的子单元形状和尺寸良好，可直接用于奇异积分的计算。
表1给出了本文方法与传统方法在不同源点下积分点数目以及相应误差的数值结果。由表1可得，在积分点数目相当的情况下，本文方法的精度高于传统方法，这说明基于仿射变换的自适应单元细分法能够准确地计算线性三角形单元的奇异积分。
表 1 线性三角形单元的积分计算结果Table 1. Numerical results of the linear triangular element单元类型 源点位置 积分点数量 相对误差 传统
方法本文
方法传统
方法本文
方法线性
三角形单元No.1 289 243 1.70×10−2 4.60×10−5 No.2 300 266 1.75×10−2 7.55×10−5 No.3 403 363 1.99×10−3 1.75×10−5 No.4 300 243 1.79×10−2 4.61×10−5 5.2 线性四边形单元
此算例展示了在不同源点位置的线性四边形单元的自适应细分结果。细长四边形单元的节点坐标分别是(0.0, 0.0, 0.0)、(6.0, 0.0, 0.0)、(6.0, 3.0, 0.0)、(0.0, 3.0, 0.0)。源点位置基本上沿图12中的直线均匀分布，源点编号在全局坐标系中对应的坐标从左至右依次为(2.3, 0.4, 0.0)、(3.3, 0.4, 0.0)、(4.3, 0.4, 0.0)、(5.3, 0.4, 0.0)。
![]()
图 12 线性四边形单元节点及源点坐标Figure 12. Coordinates of source points and nodes of the linear quadrilateral element图13显示了在不同的源点位置，线性四边形单元在仿射变换后得到的子单元。从图13中可以看出，该方法得到的积分子单元形状、尺寸良好。
表2中将本文方法与传统方法的数值结果进行对比，分别列出各自在不同源点处所需的积分点数目以及相应的误差大小。由表2可得，在相同源点位置，积分点数目相当的情况下，本文方法的精度高于传统的方法，这说明基于仿射变换的自适应单元细分法能准确地计算线性四边形单元的奇异积分。
表 2 线性四边形单元的积分计算结果Table 2. Numerical results of the linear quadrilateral element单元类型 源点位置 积分点数量 相对误差 传统
方法本文
方法传统
方法本文
方法线性四边形单元 No.1 576 432 2.32×10−3 5.63×10−5 No.2 324 283 2.72×10−2 4.26×10−5 No.3 400 300 3.00×10−2 3.76×10−5 No.4 324 288 1.64×10−2 3.08×10−5 5.3 二次三角形单元
此算例展示了不同源点位置的二次三角形单元的自适应细分结果，其节点坐标分别为(0.0, 0.0, 0.0)、(8.0, 0.0, 0.0)、(4.0, 3.0, 0.0)、(4.0, 0.0, 0.0)、(6.0, 2.0, 0.0)、(2.0, 2.0, 0.0)。如图14所示，在全局坐标系中，基本上沿直线均匀分布的4个源点坐标从左到右依次为(2.3, 0.9, 0.0)、(3.6, 1.0, 0.0)、(4.8, 1.2, 0.0)、(6.0, 1.3, 0.0)。
![]()
图 14 二次三角形单元节点及源点坐标Figure 14. Coordinates of source points and nodes of the quadratic triangular element图15显示了在不同源点位置，二次三角形单元在基于仿射变换下的自适应单元细分法细分后得到的子单元。从图15中可以看出，用该方法得到的子单元形状和尺寸良好。
数值结果如表3所示，通过将本文方法与传统方法所用积分点数量及相对误差进行对比，证明了在积分点数量相当的情况下，本文方法的精度高于传统方法，是处理二次三角形单元奇异积分的有效方法。
表 3 二次三角形单元的积分计算结果Table 3. Numerical results of the quadratic triangular element单元类型 源点位置 积分点数量 相对误差 传统
方法本文
方法传统
方法本文
方法二次三角形单元 No.1 255 243 1.35×10−3 7.03×10−5 No.2 347 308 5.68×10−3 1.91×10−5 No.3 326 300 4.45×10−2 2.32×10−5 No.4 335 323 1.94×10−3 8.41×10−5 5.4 二次四边形单元
此算例展示了二次四边形单元在不同源点处的自适应细分结果。二次四边形单元的节点坐标分别为(0.0, 0.0, 0.0)、(20.0, 0.0, 0.0)、(20.0, 2.0, 0.0)、(0.0, 2.0, 0.0)、(10.0, 0.3, 0.0)、(20.1, 1.0, 0.0)、(10.0, 2.3, 0.0)、(0.1, 1.0, 0.0)。源点的位置基本上沿图16所示的直线均匀分布，源点编号在全局坐标系中对应的坐标从左到右依次为(3.0, 0.9, 0.0)、(6.0, 0.9, 0.0)、(9.0, 0.9, 0.0)、(12.0, 0.9, 0.0)。
图17展示了在不同源点位置，二次四边形单元在细分后得到的子单元。从细分结果可以看出，通过仿射变换细分得到的子单元形状和尺寸比较好。
表4列出了本文方法和传统方法所用积分点数量以及相对误差。通过将两种方法进行对比得出，在积分点数量相当的情况下，本文方法的精度高于传统方法，是处理二次四边形单元奇异积分的有效方法。
![]()
图 16 二次四边形单元节点及源点坐标Figure 16. Coordinates of source points and nodes of the quadratic quadrilateral element![]()
图 17 二次四边形单元细分结果Figure 17. Element subdivision for the quadratic quadrilateral element表 4 二次四边形单元的积分计算结果Table 4. Numerical results of the quadratic quadrilateral element单元类型 源点位置 积分点数量 相对误差 传统
方法本文
方法传统
方法本文
方法二次四边形单元 No.1 400 328 1.44×10−2 8.65×10−5 No.2 400 313 2.80×10−2 7.70×10−5 No.3 402 350 1.65×10−2 2.60×10−5 No.4 400 337 7.88×10−2 5.04×10−5 6 结论
针对边界积分方程中存在的奇异积分难题，本文提出了基于仿射变换的奇异积分自适应单元细分法，通过与传统单元细分方法进行对比分析，得到以下结论：
(1)该方法将单元划分为投影区域和细分区域，通过仿射变换对细分区域进行划分，采用径向投影算法重构源点附近的积分区域，细分子单元形状和尺寸良好。
(2)该方法适用于不同类型的单元，能够实现任意单元形状和源点位置的高效高精度奇异积分计算，且数据结构简单，易于实现，鲁棒性好。
(3)该方法在任意情况下均能保证单元细分的收敛性和高质量的单元细分，积分点在单元内部分布更为合理，具有良好的计算精度和稳定性。
-
![]()
图 2 线性四边形单元细分流程
Figure 2. The element subdivision process of a linear quadrilateral element
![]()
图 3 基于仿射变换的自适应单元细分流程图
Figure 3. Flow diagram of adaptive element subdivision based on affine transformation
![]()
图 4 初始单元预处理及分区质量判断
Figure 4. The initial element preprocessing and domain partition quality assessment
![]()
图 5 不同源点位置单元细分方法示意图
Figure 5. Schematic diagram of element partition for different locations of source points
![]()
图 7 单元投影区域构建示意图
Figure 7. Schematic diagram of the element projection region construction
![]()
图 8 基于仿射变换细分方法对低质量细分区域的处理
Figure 8. The low quality subdivision region processing based on affine transformation subdivision method
![]()
图 10 线性三角形单元节点及源点坐标
Figure 10. Coordinates of source points and nodes of the linear triangular element
![]()
图 12 线性四边形单元节点及源点坐标
Figure 12. Coordinates of source points and nodes of the linear quadrilateral element
![]()
图 14 二次三角形单元节点及源点坐标
Figure 14. Coordinates of source points and nodes of the quadratic triangular element
![]()
图 16 二次四边形单元节点及源点坐标
Figure 16. Coordinates of source points and nodes of the quadratic quadrilateral element
![]()
图 17 二次四边形单元细分结果
Figure 17. Element subdivision for the quadratic quadrilateral element
表 1 线性三角形单元的积分计算结果
Table 1 Numerical results of the linear triangular element
单元类型 源点位置 积分点数量 相对误差 传统
方法本文
方法传统
方法本文
方法线性
三角形单元No.1 289 243 1.70×10−2 4.60×10−5 No.2 300 266 1.75×10−2 7.55×10−5 No.3 403 363 1.99×10−3 1.75×10−5 No.4 300 243 1.79×10−2 4.61×10−5 
下载: 导出CSV
表 2 线性四边形单元的积分计算结果
Table 2 Numerical results of the linear quadrilateral element
单元类型 源点位置 积分点数量 相对误差 传统
方法本文
方法传统
方法本文
方法线性四边形单元 No.1 576 432 2.32×10−3 5.63×10−5 No.2 324 283 2.72×10−2 4.26×10−5 No.3 400 300 3.00×10−2 3.76×10−5 No.4 324 288 1.64×10−2 3.08×10−5
下载: 导出CSV
表 3 二次三角形单元的积分计算结果
Table 3 Numerical results of the quadratic triangular element
单元类型 源点位置 积分点数量 相对误差 传统
方法本文
方法传统
方法本文
方法二次三角形单元 No.1 255 243 1.35×10−3 7.03×10−5 No.2 347 308 5.68×10−3 1.91×10−5 No.3 326 300 4.45×10−2 2.32×10−5 No.4 335 323 1.94×10−3 8.41×10−5
下载: 导出CSV
表 4 二次四边形单元的积分计算结果
Table 4 Numerical results of the quadratic quadrilateral element
单元类型 源点位置 积分点数量 相对误差 传统
方法本文
方法传统
方法本文
方法二次四边形单元 No.1 400 328 1.44×10−2 8.65×10−5 No.2 400 313 2.80×10−2 7.70×10−5 No.3 402 350 1.65×10−2 2.60×10−5 No.4 400 337 7.88×10−2 5.04×10−5
下载: 导出CSV
-
[1] GAO X W. An effective method for numerical evaluation of general 2D and 3D high order singular boundary integrals [J]. Computer Methods in Applied Mechanics and Engineering, 2010, 199(45/46/47/48): 2856 − 2864.
[2] LI C, HU B, NIU Z R. Effectiveness of the stress solutions of 3-D V-notched/cracked structures by using extended boundary element method [J]. Journal of Engineering Mathematics, 2022, 135(1): 5. doi: 10.1007/s10665-022-10228-5
[3] YANG K, LI H Y, PENG H F, et al. New interface integration BEM for solving multi-medium nonlinear heat transfer problems [J]. Engineering Analysis with Boundary Elements, 2020, 117: 66 − 75. doi: 10.1016/j.enganabound.2020.03.015
[4] 左冲, 姚鸿骁, 姚伟岸. 时域径向积分边界元法在平面单相凝固问题中的应用[J]. 工程力学, 2019, 36(3): 33 − 39. doi: 10.6052/j.issn.1000-4750.2018.01.0003 ZUO Chong, YAO Hongxiao, YAO Weian. The application of the time domain radial integral boundary element method in 2D one-phase solidification [J]. Engineering Mechanics, 2019, 36(3): 33 − 39. (in Chinese) doi: 10.6052/j.issn.1000-4750.2018.01.0003
[5] 王富顺, 池宝涛, 贾志超, 等. 一种适用于邻近特征的非结构体网格生成方法[J]. 机械工程学报, 2023, 59(24): 118 − 130. doi: 10.3901/JME.2023.24.118 WANG Fushun, CHI Baotao, JIA Zhichao, et al. An unstructured mesh generation algorithm for arbitrary geometry with proximity features [J]. Journal of Mechanical Engineering, 2023, 59(24): 118 − 130. (in Chinese) doi: 10.3901/JME.2023.24.118
[6] 谢伟, 贺旭东, 吴建国, 等. 二维光滑边域有限元法在弹性力学中的应用研究[J]. 西北工业大学学报, 2017, 35(1): 7 − 12. doi: 10.3969/j.issn.1000-2758.2017.01.002 XIE Wei, HE Xudong, WU Jianguo, et al. An edge-based smoothed finite element method for 2D mechanics problems [J]. Journal of Northwestern Polytechnical University, 2017, 35(1): 7 − 12. (in Chinese) doi: 10.3969/j.issn.1000-2758.2017.01.002
[7] ZHANG J M, CHI B T, SINGH K M, et al. A binary-tree element subdivision method for evaluation of singular domain integrals with continuous or discontinuous kernel [J]. Engineering Analysis with Boundary Elements, 2020, 116: 14 − 30. doi: 10.1016/j.enganabound.2020.03.023
[8] 高效伟, 冯伟哲, 杨恺. 边界元中计算任意高阶奇异线积分的直接法[J]. 力学学报, 2014, 46(3): 428 − 435. GAO Xiaowei, FENG Weizhe, YANG Kai. A direct method for evaluating line integrals with arbitrary high order of singularities [J]. Chinese Journal of Theoretical and Applied Mechanics, 2014, 46(3): 428 − 435. (in Chinese)
[9] 邓琴, 李春光, 王水林, 等. 边界元中近奇异积分的一种解析方法[J]. 工程力学, 2010, 27(9): 49 − 54. DENG Qin, LI Chunguang, WANG Shuilin, et al. An analytical method for nearly singular integrals in boundary element method [J]. Engineering Mechanics, 2010, 27(9): 49 − 54. (in Chinese)
[10] CHI B T, GUO Q J, ZHANG L G, et al. An adaptive binary-tree element subdivision method for evaluation of volume integrals with continuous or discontinuous kernels [J]. Engineering Analysis with Boundary Elements, 2022, 134: 298 − 314. doi: 10.1016/j.enganabound.2021.10.010
[11] 池宝涛, 张见明, 鞠传明. 基于T-Spline的全自动几何拓扑修复方法[J]. 自动化学报, 2019, 45(8): 1511 − 1526. CHI Baotao, ZHANG Jianming, JU Chuanming. An automatic topology recovery method using T-Spline [J]. Acta Automatica Sinica, 2019, 45(8): 1511 − 1526. (in Chinese)
[12] 牛忠荣, 胡宗军, 葛仁余, 等. 二维边界元法高阶元几乎奇异积分半解析算法[J]. 力学学报, 2013, 45(6): 897 − 907. doi: 10.6052/0459-1879-13-215 NIU Zhongrong, HU Zongjun, GE Renyu, et al. A new semi-analytic algorithm of nearly singular integrals in high order boundary element analysis of 2D elasticity [J]. Chinese Journal of Theoretical and Applied Mechanics, 2013, 45(6): 897 − 907. (in Chinese) doi: 10.6052/0459-1879-13-215
[13] 李佳龙, 李钢, 余丁浩. 多边形比例边界有限元非线性高效分析方法[J]. 工程力学, 2020, 37(9): 8 − 17. doi: 10.6052/j.issn.1000-4750.2019.10.0634 LI Jialong, LI Gang, YU Dinghao. Polygon scaled boundary finite element nonlinear efficient analysis method [J]. Engineering Mechanics, 2020, 37(9): 8 − 17. (in Chinese) doi: 10.6052/j.issn.1000-4750.2019.10.0634
[14] 牛忠荣, 王左辉, 胡宗军, 等. 二维边界元法中几乎奇异积分的解析法[J]. 工程力学, 2004, 21(6): 113 − 117. NIU Zhongrong, WANG Zuohui, HU Zongjun, et al. Analytic algorithm for nearly singular integrals in two-dimensional boundary element analysis [J]. Engineering Mechanics, 2004, 21(6): 113 − 117. (in Chinese)
[15] TELLES J C F. A self-adaptive co-ordinate transformation for efficient numerical evaluation of general boundary element integrals [J]. International Journal for Numerical Methods in Engineering, 1987, 24(5): 959 − 973. doi: 10.1002/nme.1620240509
[16] MA H, KAMIYA N. A general algorithm for the numerical evaluation of nearly singular boundary integrals of various orders for two-and three-dimensional elasticity [J]. Computational Mechanics, 2002, 29(4): 277 − 288.
[17] 李光耀, 何建平, 董云桥, 等. 基于球面细分法的高效高精度近奇异积分计算[J]. 固体力学学报, 2018, 39(5): 453 − 461. LI Guangyao, HE Jianping, DONG Yunqiao, et al. Evaluating nearly singular integrals accurately and efficiently based on sphere subdivision method [J]. Chinese Journal of Solid Mechanics, 2018, 39(5): 453 − 461. (in Chinese)
[18] ZHANG J M, CHI B T, SINGH K M, et al. A binary-tree element subdivision method for evaluation of nearly singular domain integrals with continuous or discontinuous kernel [J]. Journal of Computational and Applied Mathematics, 2019, 362: 22 − 40. doi: 10.1016/j.cam.2019.04.027
[19] SLADEK V, SLADEK J, TANAKA M. Optimal transformations of the integration variables in computation of singular integrals in BEM [J]. International Journal for Numerical Methods in Engineering, 2000, 47(7): 1263 − 1283. doi: 10.1002/(SICI)1097-0207(20000310)47:7<1263::AID-NME811>3.0.CO;2-I
[20] XIE G Z, ZHOU F L, ZHONG Y D, et al. Bi-directional sinh transformations based on the generalized Duffy space for nearly singular integrals [J]. Journal of Computational and Applied Mathematics, 2020, 380: 112981. doi: 10.1016/j.cam.2020.112981
[21] ZHANG J M, CHI B T, LIN W C, et al. A dual interpolation boundary face method for three-dimensional potential problems [J]. International Journal of Heat and Mass Transfer, 2019, 140: 862 − 876. doi: 10.1016/j.ijheatmasstransfer.2019.06.011
[22] 池宝涛, 张见明, 鞠传明. 基于仿射算术和区间运算的直线与NURBS曲线/曲面求交[J]. 中国机械工程, 2019, 30(9): 1026 − 1033. doi: 10.3969/j.issn.1004-132X.2019.09.003 CHI Baotao, ZHANG Jianming, JU Chuanming. An improved fast algorithm for solving intersections of lines and NURBS curves and surfaces using interval analysis with affine transform [J]. China Mechanical Engineering, 2019, 30(9): 1026 − 1033. (in Chinese) doi: 10.3969/j.issn.1004-132X.2019.09.003
计量
- 文章访问数: 285
- HTML全文浏览量: 37
- PDF下载量: 44
