比例边界计算机有限元算法理论及研究发展分析

余芝轩

(黄冈职业技术学院，湖北 黄冈 438002)

1 引言

随着社会的发展与进步，新技术、新材料、新方法不断出现，人类在社会生产活动中所面临的工程日益复杂，越来越多的数值计算方法正在被人们所提出、接受以及研究。对于每一种数值计算方法，计算精度、计算效率以及适用范围往往都是研究的重点。

计算机结合比例边界有限元法是由John.P.Wolf 和Chongmin Song提出的一种结合了边界元和有限元诸多优点的数值计算方法[1]。目前已有学者对此种方法进行了较为深入的研究，在工程领域进行了应用，并且和传统的数值计算方法进行了对比。结果表明，对于某些工程实际问题，相比较于边界元和有限元，此种方法可以显著地减少计算时间，提高计算精度[2-3]。

2 比例边界有限元法的理论基础

对于二维弹性静力学来说，应力{σ(x,y)}和体积力{p(x,y)}在域内任一点应满足如下关系：

[L]T{σ(x,y)}={p(x,y)}

(1)

这里[L]是线性微分算子。

应变{ε(x,y)}和位移{u(x,y)}的表达式如下：

{ε(x,y)}=[L]{u(x,y)}

(2)

而应力与应变的关系表达式如下：

{σ(x,y)}=[D]{ε(x,y)}

(3)

其中[D]是弹性矩阵。

图1 比例边界有限元坐标系统

在推导比例边界有限元方程时，需在比例中心(x0,y0)建立相应的坐标系。比例中心可以在对所有边界可见的前提下放在域内任意位置[4]。比例边界元的坐标系统包括径向坐标ξ和环向坐标s。径向坐标规定，在比例中心处ξ=0，在边界上ξ=1，域内其他地方取0<ξ<1；环向坐标指沿着边界逆时针方向的距离。

笛卡尔坐标下的点(x,y)可以由ξ和s表示为：

x=x0+ξxs(s)=x0+ξ[N(s)]{x}y=y0+ξys(s)=y0+ξ[N(s)]{y}

(4)

运用雅克比矩阵变化，可得：

(5)

其中，雅克比矩阵

(6)

在ξ=1的边界上，式(6)的[J(ξ,s)]将变成[J(s)]，且[J(s)]只与边界上的节点坐标有关，矩阵形式为：

(7)

即

[J(s)]=x(s)y(s)s-y(s)x(s)s

(8)

结合式(5)和(6)，可得：

(9)

其中，

|J(ξ,s)|=x(s)ξy(s)s-y(s)ξx(s)s=ξ[x(s)y(s)s-y(s)x(s)s]=ξ|J(s)|

(10)

将式(10)带入式(9)，可得：

(11)

对于二维的线弹性问题，

(12)

将式(11)代入(12)，可得：

(13)

其中，令

(14)

(15)

将式(14)和(15)代入(13)，可得：

(16)

同时，令

[B1(s)]=[b1(s)][N(s)],[B2(s)]=[b2(s)][N(s)]s

(17)

采用如下的虚功原理表达式进行相应的变化和推导：

(18)

对式(17)进行相应的积分和运算，可得到比例边界有限元下关于{uh(ξ)}ξξ、{uh(ξ)}ξ、{uh(ξ)}以及系数[E0]、[E1]、[E2]的Euler-Cauchy位移平衡微分方程，并可得到相应解为

{uh(ξ)}=c1ξ-λ1{φ1}+c2ξ-λ2{φ2}+…+ciξ-λi{φi}

(19)

其中，ci代表的是每一种独立的位移模式对于方程解的加权系数，-λi则表示了对于径向坐标的比例因子。

通过对方程自由度数的翻倍，可构建一个关于不同模态下的位移以及节点力的特征方程并进行求解。

对于式(19)中的ci，一般通过在边界上(ξ=1)的方程{c}=[Φ1]-1{uh}来求得。并最终可以得到位移和应力的表达式分别为：

(20)

(21)

3 比例边界有限元法研究进展

比例边界有限元法(SBFEM)是一种基于计算机有限元方法，并且结合了有限元和边界元诸多优点的较为新型的数值方法。相比较于边界元，比例边界有限元同样只需要在边界上进行离散，但它针对边界元无法生成对称刚度矩阵的问题做出了相应改进。目前，已有研究和工程实践表明，比例边界有限元对于解决应力奇异性问题以及无限域问题，相对于传统的有限元方法来说具有很大的优势。对于应力奇异性问题，比例边界有限元法不需要进行额外的处理。在无限域问题上，由于SBFEM只需在边界上进行离散，从而可以在处理相关问题时通过更少的自由度数来提高计算效率，减少计算时间[5]。Hauke Gravenkamp、Carolin Birk和Chongmin Song在2015年对弹性导波在带有缺陷的无限域板中的传播进行了模拟，并在建模效率和计算精度这两方面将模拟结果与传统有限元在同等条件下的模拟结果进行了比较。结果表明，在精度相同的情况下，比例边界有限元的建模效率可以提高两到三倍。

图2 比例边界有限元关系图

除此之外，比例边界有限元在裂缝问题的处理上也有独到的优势，其可以通过将比例中心设置在裂尖上来解决应力奇异性的问题。而且裂缝位置的变化以及裂缝深度的改变也仅需要进行少量的网格重划分即可完成，相比较于有限元法来说大大减少了工作量。Hauke Gravenkamp、Jens Prager、Albert A.Saputra和Chongmin Song对Lamb波在带有裂缝板中的传播进行了模拟，并在自由度数(DOF)、误差(Deviation)和运行时间(CPU)这三个方面对SBFEM与传统有限元的模拟结果进行了综合比较，结果表明在误差率一致的情况下，比例边界有限元的自由度数仅有传统有限元法的5.88%，运行时间仅有传统有限元法的1.50%。从以上工程实例可以得知，在裂缝问题的模拟以及计算过程中，比例边界有限元的计算效率相比较于传统有限元法有了较为明显的提升。

正是由于比例边界有限元法(SBFEM)在各种工程实际问题中得到了不断的应用并且相比较于传统有限元方法在某些方面有了较为明显的进步，所以近年来，比例边界有限元法得到了越来越广泛的认可、研究、应用以及发展。

比例边界有限元最初是用来计算无边界介质的动力刚度矩阵，然后Chongmin Song和J.P.Wolf又用此方法对多项异性材料进行裂缝分析[6]。随后，在此前工作的基础之上，J.P.Wolf和A.J.Deeks对比例边界有限元的原理以及坐标变化做了系统的阐述，即推导了从笛卡尔坐标系统到径向、环向坐标系统的转变，且进一步扩宽了其适用范围。随后Chongmin Song对静力学问题又提出了一种新的矩阵函数计算方法，并表明此种计算方法具有更高的计算精度[5]。大连理工大学的林皋院士和他的博士生杜建国将比例边界有限元应用于计算坝面动水压力，通过对二维以及三维的动水压力进行计算，并且和有限元的计算结果进行对比发现，比例边界有限元可以在提高计算效率、减少工作量的同时减小计算误差[7]。Hauke Gravenkamp和Albert A.Saputra等人首次将比例边界有限元应用于模拟兰姆波在带裂缝板中的传播，结果表明对于裂缝问题，比例边界有限元在计算误差减小、计算时间大幅度缩减的同时，计算模型仍然具有比有限元少的自由度数[8]。之后，Hauke Gravenkamp和Carolin Birk等人继续模拟了弹性导波在任意边界的带缺陷条形板中的传播并推导出了动力刚度矩阵[9]。Ean Tat Ooi和Chongmin Song等人又提出一种基于比例边界的多边形网格划分方法，和以往的研究对象均是线弹性材料不同，此方法对弹塑性材料进行了模拟、分析，从而大大扩宽了比例边界有限元方法在未来的应用领域[10]。此后，陈灯红和戴上秋等用比例边界有限元做了土体-结构交互作用的动态断裂分析，并且得出了相对于传统方法更加准确和有效的计算结果[11]。

4 总结

本文针对目前一种较新型的数值计算方法——比例边界有限元法(SBFEM)，对其理论基础、坐标变化等方面做了推导，对比例中心的选取及适用范围做了详细的阐述。目前已有研究以及相应的工程算例在计算精度、计算效率等方面将SBFEM和传统的有限元方法进行了比较。计算结果表明由于SBFEM仅在边界上进行离散，具有降低维度、减少自由度的特点，故可以在提高计算效率的同时减小计算误差。

随着研究的深入，比例边界有限元法在解决应力奇异性、无限域以及裂缝问题上已不断地显示出良好的适应性、便利性和快捷性。相信随着研究的不断深入与在工程领域的不断应用，结合计算机比例边界有限元法会在今后得到越来越广泛的应用。