Meshless 1无网格方法介绍

第一章整体弱形式、加权余量、有限元、边界元和局部弱形式1.1引言微分方程（一维物理域的常微分方程和二维或三维的偏微分方程）和边界条件是工程中物理过程的扼要的、通常也是高度近似的特征。通常只有对于定义在昀简单的几何域上的线性问题（即微分方程和边界条件、初始条件中只出现未知变量及其导数的线性组合）才有可能求得精确的数学解。但是，大多数工程过程发生在复杂的几何域上，而且通常包含未知变量及其导数的非线性。为了求解这样的实际问题，来辅助工程分析、设计和综合，必须采用基于数字计算机建模和模拟技术。我们将此称为基于计算机建模和模拟的工程。其目的是利用功能强大的现代数字计算机、高速科学计算、可视化技术、虚拟现实和伴随的信息技术来达到对于复杂工程现象的接近实时的模拟，来辅助真实寿命工程系统的产品和过程完整性设计。这些系统包括飞机、飞船、火箭、汽车、船舶、电站、桥梁、生物系统、微电子装置、计算机芯片、微机电系统、包含传热、传质的系统、压缩机、汽轮机、包含在闭域和开域流体流动的系统、大气污染、可持续生态技术、工程材料的纳米、微米、细观和宏观力学等。既然通过适当定义每种过程的未知变量，在这些过程所满足的微分方程中有大量的共性，人们能否想象一种通用的工具箱，使能用于各种各样工程过程的计算机建模和模拟。在过去三十年中工程计算取得了令人吃惊的进展。其中昀成功的方法是有限元法和边界元法。但是，与其可靠的理论基础相比，有限元法还有诸如闭锁（剪切、薄膜、不可压缩等）和在裂纹扩展、剪切带形成、大变形等问题中费时的网格划分和重划分等缺点。由于这些原因，近年在求解边值问题中无网格法作为克服其中某些缺点的新途径而得到了相当的注意。无网格法，作为克服有限元法和边界元法的熟知缺点的可供选用的数值方法，由于它的适应性、更重要的是由于它能免除在域中构造几何网格的费时费力的过程，在昀近十年倍受注意。无网格法的主要目标是摆脱，至少是减轻在整个结构划分和重分网格的困难，代之以只要在整个结构中增加或消除节点。无网格法还能减轻某些与有限元法有关的其它问题，例如闭锁、单元畸变等等。无网格法的原始思想可追溯到用于模拟天体物理现象的平稳粒子流体力学（SPH）法(Gingold,Monaghan,1977).仅当Nayroles,Touzot,Villon(1992)发表了扩散单元法之后，无网格法的研究才变得非常活跃。其后又发表了许多所谓无网格法：Belytschko,Lu,Gu(1994)的无单元伽辽金(EFG)法，Liu,Chen,Uras,Chang(1996)的再生核粒子法(RKPM)，Babuska,Melenk(1997)的单位分解有限元法(PUFEM)，Duarte,Oden(1996)的hp-云法，Suku-mar,Moran,Belytschko(1998)的自然单元法(NEM)，Wendland(1999)的利用径向基函数的无网格伽辽金法(RBF)。在本专著的书末提供了对于无网格法详细的文献目录。这些无网格法的主要差别只是来源于试探函数插值采用的技术。虽然在这些方法中对于求解变量的试探函数和检验函数的插值都不需要网格，但在这些方法中对于弱形式或能量的积分采用影子单元是不可避免的。因此，这些方法不是真正的无网格。昀近，对于求解线性和非线性边值问题，Zhu,Zhang,Atluri(1998a,b),Atluri,Zhu(1998a,b)andAtluri,Kim,Cho(1999)发展了两种真正的无网格法，即无网格局部边界积分方程(LBIE)法和无网格局部彼得洛夫－伽辽金(MLPG)法。在这两种方法中无论为了对求解变量的试探函数和检验函数的插值，或者为了对弱形式的积分，在区域或边界都不需要网格，因此这两种方法是真正的无网格法。所有相关的积分很容易在互相覆盖的、形状规则的区域（在三维问题中通常是球）上求积。实际上，LBIE法可以简单地看作MLPG法(Atluri,Kim,Cho,1999)的一种特例。Zhang,Yao(2001)综合LBIE和边界点法的优点，提出了一种新的规则杂交边界点法。MLPG法令人注目的成功在一些问题中已经得到了报道，其中包括：求解对流－扩散问题[Lin&Atluri(2000)]，断裂力学问题[Kim,Atluri(2000),Ching,Batra(2001)]，Navier-Stokes流动问题[Lin,Atluri(2001)]，剪切变形梁[Cho,Atluri(2001)]以及板弯曲问题[Gu&Liu(2001),Long&Atluri(2002)]等。总之，MLPG是一种真正的无网格法，它不仅包含对于试探函数的一种无网格插值（诸如移动昀小二乘MLS、单位分解PU、Shepard函数或RBF），而且包含一种弱形式的无网格积分（即所有积分总是在规则形状的子域，例如三维问题中的球、平行六面体或椭球上求积）。本专著的目的是展示MLPG法，包括它的一般原理和各种实施方案。可以看到，MLPG法的某些方案不仅完全消除了区域划分网格的人工成本，而且其计算成本也比有限元法和边界元法低得多。1.2整体弱形式和加权余量法(WRM)为了搭建展示MLPG法的舞台，我们通过一个例题讨论利用整体弱形式的各种建模策略来开始。可以看到，支配工程中许多线性边值问题的数学描述可以写成如下形式：其中mL是一个m阶的线性微分算子。mL可以是一个一维问题的常微分算子，或者一个二维、三维物理域的偏微分算子。在一般的三维域中，设想一条穿过该域、连结边界两点的直线。沿着这条直线，可以设想一个m阶的常微分方程。然后，在边界每点有m/2个边界条件，每个边界条件包含一个不高于m−1阶的微分算子。图1.1一个实际的边值问题例如，我们考虑一个用笛卡尔坐标ix定义的二维空间内的线性Poisson方程。考虑域Ω如图1.1所示。域的边界记作ΩΓ(或者有时记作∂Ω)，由曲线A-B-C-D-E-F-A,以及曲线GH定义。沿着AB段，给定边界条件φφ=。类似地，沿着其它各段的边界条件为：沿着BC段nqφ∂∂=；沿着CD段一部分为φφ=，另一部分为nqφ∂∂=；沿着切口（或裂纹）D-E-F0nφ∂∂=；沿着FA段一部分为φφ=，另一部分为nqφ∂∂=；沿着GH孔边0nφ∂∂=。在图1.1中定义的边值问题是相当实际的工程问题。因此，AB属于本质边界uΓ，而EF,ED和GH属于自然边界。Poisson方程可以写成：qΓ)iiR其中(x=xe是域内点的位置矢量，Ωp是一个给定的源函数，而域被带有如下边界条件的边界所包围：ΩuΓ=Γ∪Γq条件(1.3a)常称为Dirichlet边界条件；而(1.3b)称为Nermann边界条件。许多工程领域中出现Poisson方程问题，例如：任意截面固体直杆受扭的的小扭转，均匀介质中的传热，渗流，流体膜润滑，静电场，静磁场等等。令为问题(1.2)的一个近似解，下面称之为试探函数。既然u可能不是精确解，将代入微分方程和边界条件(1.2-1.3)，我们得到如下误差余量的公式：uu域内误差：边界误差：其中IR是域内误差，而(1,2BiRi)=是边界条件的误差。“昀好”近似就是以某种方式使误差IR和BiR为零。通常，试探函数必须满足若干昀低要求来保证消除误差的过程有意义。满足这种昀低要求的函数称为容许函数。例如，在方程(1.4)中，u必须二阶可微，而必须连续（或u是C连续函数）；如果u不是连续的，在u不连续处u将为无限大。在哪种情况下，在这些点,iuii11C,i,()Iu是无限的，()IRu取昀小值就不再有任何意义。考虑个容许函数N()1,2,,iuiN=…的一个线性组合，使其中()1,2,,iaiN=…为待定系数。我们采用了对于重复指标的求和约定。这里，在任意形状的整体域Ω上选择容许函数()u[在方程(1.4)中x()ux为连续]是一个重要的命题，本专著的第2章主要讨论这个问题。但是，如果域是简单几何形状的，诸如圆、正方形等，1C()ux可在三角函数或其它函数的形式中选择。图1.2:配点法示意图我们现在讨论几种使误差IR和BiR尽可能小、在极限情况使误差在整个域上均匀地趋于零的方法。Ω1.2.1配点法我们假设是方程(1.4)中在整个域uΩ上C连续的。我们考虑使误差在域1Ω内的一些点,1,2,,jjM=x…处为零，由此可导出如下对于待定系数的方程：和假设我们利用两种形式的边界余量为零[方程(1.8)]来建立个方程，那么，利用域内的NM−ΩM个配点，我们可以建立对于个未知系数的个方程。虽然这样一个方程组的解同时使误差NiaNIR和BiR在有限个点处等于零，但是余量IR和BiR在这些点之间可能振荡，而且可能真的相当大。利用比M更多的配点将建立多于未知量数的方程。在这种情况下，通常不可能求得一组系数使在配点处的误差为零。ia*ja假设方程(1.7)的配点取任意大的配点数，同时强制(1.8)式在任意点数处满足，通常会导致超定方程组：其中，且。我们寻求的一个近似解，使1,2,,;1,2,,iKj==……NkNja而使误差的平方*jaiiεε昀小。于是或或其中，而1,2,,;1,2,,iKjK==……1,2,,qN=…。当jqB正定时，上述方法是求解超定方程组的一种代数昀小二乘法。但是，对于不定方程组，只能得到部分未知量用其它未知量表示的关系。可以指出，方程(1.11)也等于使逐点相加的误差平方和昀小。1.2.2误差平方的加权积分如同前面讨论的那样，在配点法中，我们可以利用多于试探函数中未知系数数的配点，采用误差昀小二乘法。作为极限，当配点数趋于无限时，我们可以考虑误差平方的积分：其中σβγ,,为加权常数，也常称为罚参数。将上式ε对取极小，就可得到欲求的个代数方程的方程组来确定。iaNia1.2.3子域误差积分/平均法作为另一种选择，我们可以考虑使误差在一组每个子域上的平均值为零。例如，我们可令图1.3子域示意图子域数的选择使我们得到对于,1,2,,iaiN=…N的一个可解的方程组（即选个子域或更多子域，以得到个方程或多于个方程）。图1.3是示意图。注意到子域NN()1,2,,iiΩ=…k是邻接而不覆盖的，并且布满整个域。其中大部分是规则形状的（比如三角形），其它的也接近于三角形。ΩiΩ1.2.4有限体积法此法和子域积分法相似，并且在非定常空气动力学和热传导方面的文献中已经很普及。有限体积法的原理是在一个控制体积（即子域）kΩ上积分方程(1.2)，得到其中是子域k∂ΩkΩ的边界，Ω∈，其次，流量要在控制体积的边界上近似。因此，我们将近似网格每边上的的积分。再次注意到，子域是彼此kΩk∂ΩkΩ,iuni邻接不覆盖的，如图1.3所示。有限体积法的另一种方案可以将方程改写成混合形式来导出：其中φ和*iφ看作独立的变量。如果u和u分别为*iφ和*iφ的试探函数，我们可以写出类似(1.16)式的公式：以及代替整体试探函数()kkuau=x和()**ikkiuau=x，可以假设和在每个单元上的一个线性插值，如图1.3所示。再次注意，图1.3中的所有子域是邻接而不覆盖的，这种方法无论为了插值，还是为了积分，都包含着要生成网格，如方程(1.18),(1.19)所示。u*iukΩ1.2.5加权余量法和整体弱形式如果函数()Fx对于域上的所有Ω()gx均满足条件那么，在Ω上()0F≡x。因此，我们可以反过来求试探函数()iiu=ua，通过对域上的所有可能函数Ω()vx强制满足弱形式方程，使()0IRu≡（其中IR在方程(1.7)中定义）。我们称()vx为检验函数。通常，试探函数u选为一组基函数()x()1,2,,iNui=…的带有待定系数的线性组合，即ia(),iiuua=∈Ωxx。同样地，检验函数()vx选为，，其中为基函数，b为待定系数。其“昀好”近似通过u的组合满足对于选定检验函数的弱形式方程(1.21)来确定。对于试探函数和检验函数的基函数的不同选择导致不同的近似方法。例如，当检验函数采用Dirac(jjvbv=xi)1,2,,jM=…jvj∆函数()()jvδ=−xxx时，这就归结为配点法。当如下分