（固体力学专业论文）渗流自由面问题的无网格方法.pdf

摘要 本文主要研究无单元伽辽金法( e f g ) 和无单元彼得罗夫一伽辽会法( m l p g ) 在渗流自由面问题中的应用，并针对具体实施过程中出现的问题提出了相应的解 决方法，同时研究了不同参数对计算结果的影响。 首先，从加权残量法的角度出发，通过加权残量法和变分原理推导出e f g 和m l p g 方法的基本形式，同时给出了用移动最小二乘法建立近似场函数的方 法，并采用多项式正交基解决导致求解过程失败的问题。 其次，基于无网格方法的基本思想，导出了渗流自由面问题的无单元伽辽会 法和无单元彼得罗夫一伽辽金法的控制方程。无单元彼得罗夫一伽辽金法是首次 应用于渗流自由面问题的计算。对于计算中涉及的自由面拟合和自由面边界积 分，提出了数值处理方法。计算结果同试验数据吻合得很好。同其他几种数值方 法比较，无网格方法所得到的计算结果误差小，精度高。这表明，对于渗流自由 面问题，无网格方法是行之有效的计算方法。通过数值算例来分析比较节点分布、 权函数的选取、正交多项式基的次数以及初始自由面选取等因素对于计算精度和 效率的影响，并为以后的渗流无网格方法研究提供经验性总结。 无网格方法具有计算精度高、前处理过程简单、无需划分单元等优点，适用 于动边界的渗流问题计算。但是计算量较大、计算效率较低等问题都是无网格方 法有待进一步研究改进的地方。 【关键词】无网格法，无单元伽辽金法，无单元彼得罗夫一伽辽金法 渗流，自由面，移动最小二乘法 【中图分类号l 0 3 0 2 a b s t r a c t t w om e s h f r e em e t h o d s e f gm e t h o da n dm l p gm e t h o d f o r s e e p a g ea n a l y s i sw i t ha f r e es u r f a c ea r ep r e s e n t e di nt h i sp a p e r s o m e n u m e r i c a la p p r o a c h e sa r ec a r r i e do u tt og e ts a r i s f a c t o r yr e s u l t s s e v e r a l n u m e r i c a le x a m p l e sw i l lb ep r o v i d e dt os h o wt h ee 虢c to ft h eo r d e ro f b a s i s ，w e i g h tf u n c t i o n sa n de t c o nt h ec o m p u t a t i o n a lp r e c i s i o na n dt h e c o m p u t a t i o n a le f f i c i e n c y t h i sp a p e rp r o v i d e sab r i e fi n t r o d u c t i o nt ot h em e s hf r e em e t h o d s i n c l u d i n gb a c k g r o u n da n dm o t i v a t i o n st h a th a v el e a dt ot h ed e v e l o p m e n t o f 匝r e em e t h o d s g a l e r k i n w e i g h t e dr e s i d u a lm e t h o d c o n s t r a i n e d w e a kf o r ma n dv a r i a t i o n a lp r i n c i p l ea r eu s e df o rc r e a t i n gd i s c r e t i z e d s y s t e me q u a t i o na n dm f r e ef o r n a u l a t i o n ad e t a i l e dd e s c r i p t i o no fe f g a n dm l p gi sg i v e n t h em e s h l e s sl o c a lp e t r o g a l e r k i n ( m l p g lm e t h o d i si n t r o d u c e dt ot h ep r o b l e mo fs e e p a g ef o rt h ef i r s tt i m e i ti s d e m o n s t r a t e dt h a t 几p gi sau s e f u lt o o if o rs t u d yo fs e e p a g ew i t hf r e e s u r f a c e a n dan u m b e ro ft e c h n i q u e s ，s oa su s i n gm o v i n gl e a s t - s q u a r e a p p r o x i m a t i o nt oc o n s t r u c ts h a p ef u n c t i o n s ，a r ea l s od i s c u s s e di nd e t a i l s r e em e t h o d sf o rs e e p a g ea n a l y s i sw i t haf l e es u r f a c ep r e s e n t e dc a l l a v o i dt r o u b l e s o m em o d i f i c a t i o n so ft h em e s ha si nt h ef i n i t ee l e m e n t m e t h o d s t e a d ys e e p a g ei na u n i f o r me a r t hd a mi sa n a l y z e ds e p a r a t e l yb y e f gm e t h o da n dm l p gm e t h o d t h r o u g hn u m e r i c a le x p e r i m e n t s ，t h e e f f e c t i v e n e s so ft h e s et w om e t h o d sf o rs e e p a g ea n a l y s i si sv e r i f i e d t h e n u m e r i c a lr e s u l t sa r ei ng o o da g r e e m e n tw i t ht h ee x p e r i m e n t a lo n e s k e y w o r d s ：m e s h f r e em e t h o d ，e f g ，m l p g ，s e e p a g e ，f l e e s u r f a c e ， m o v i n g - l e a s ts q u a r ea p p r o x i m a t i o n 第一章绪论 第一章绪论 1 1 无网格方法概述 1 1 1 提出无网格方法的背景 在工程计算中，计算机模拟技术通常被用来建立物理模型，通过数值方法求 解相应的微分方程以研究、解释和预测各种物理现象。在数值计算中，常用的求 解方法包括有限元法、有限差分法和有限体积法等。以上这些方法都是基于网格 的数值方法，在分析涉及特大变形( 如加工成型、高速碰撞、流固耦合) 、奇异 性、腐蚀、裂纹动态扩展和可动边界等问题时，遇到了很多困难。在用l a g r a n g e 方法求解金属冲压成形、高速碰撞、流固耦合和局部化等涉及特大变形的问题时， 网格可能会产生严重变形甚至扭曲，从而大大降低了求解的精度。为避免这“睛 况，求解过程需要不断地重构网格。这不仅浪费大量的机时，而且对于三维问题， 技术上还存在一定的障碍；对于高速冲击等动态问题，显式时间积分的步长取决 于网格的最小尺寸，网格的扭曲将使时间积分步长过小，大幅度地增加计算工作 量；对于裂纹的动态扩展闯题，出于裂纹的扩展方向不能事先确定，因此在计算 过程中需要不断地重新划分网格以模拟裂纹的动态扩展过程。这些阻网格为基础 的方法要保证边界和裂缝不能穿过单元网格，对于边界变化和裂缝生长这类问 题，必须在计算的每个时间步内重新划分网格。复杂的三维结构的网格生成也是 基于网格的计算方法的一大问题。 不同于有限元法，无网格法的近似函数建立在一系列离散点上，用全域场变 量函数进行插值近似，不需要借助于网格( 至少插值函数不需要网格) 。这一特 点使得无网格法克服了对于网格的依赖，摆脱不连续性对问题的束缚( 如网络的 重构等) ，不仅保证了求解的精度，而且大大减少前处理的工作时间，省去后处 理中对不连续结果的光顺化处理工作。 1 1 2 无网格方法的研究进展 对无网格法的研究可以追溯到2 0 世纪7 0 年代对非规则网格有限差分法的研 耕1 t 舶，1 9 7 7 年l u e y t 和g i n g o l d 等分别提出了光滑粒子流体动力学法( s m o o t h e d p a r t i c l eh y d r o d y n a m i c s ，简称s p h ) 3 一，并成功地应用于天体物理领域中。但在 很长一段时期内，这一无网格方法的思想直没有引起人们的重视，这主要是因 为有限元、有限差分等方法在工程计算中取得了极大的成功。直到十多年前，随 第一章绪论 着网格类方法中问题的逐渐暴露，人们开始寻求新的替代方法。1 9 9 2 年n a y r o l e s 等将移动最小二乘近似1 5j ( m o v i n gl e a s ts q u a r e ，简称m l s ) 引入g a l e r k i n 法， 提出了散射元法( d i f f u s ee l e m e n tm e t h o d ，简称d e m ) 忡j 。b e l y t s c h k o 等对d e m 进行了改进，并利用拉格朗日乘子法引入本质边界条件，提出了无单元g a l e r k i n 法( t h ee l e m e n t f r e eg l e r k i nm e t h o d ，简称e f g ) 1 7 1 。较s p h 方法，这类方法计 算成本高，但具有较好的稳定性。e f g 方法的计算精度和收敛速度都高于有限 元，但是计算量较大，且需要借助于背景网格进行数值积分。1 9 9 8 年，a t l u r i 等 提出了无网格局部p e t r o v g a l e r k i n 法( m e s h l e s sl o c a lp e t r o v g a l e r k i nm e t h o d ， m l p g ) 8 - t l 】。m l p g 方法在局部子域内建立无网格方程，积分时不需要背景网 格，是比较理想的无网格方法。 表1 - 1 几种无网格方法的比较 无网格方法名称场函数近似方法离散方法背景网格 无单元g a l e r k i n 法移动最小二乘法g a l e r k i n 法有 局部p e t r o g a l e r k i n 无网格法移动最小二乘法p e t r o v - g a l e r k i n 法无 散射单元法移动最小二乘法 g a e r k i n 法 有 有限点法移动最小二乘法配点法 无 无网格配点法1 73 重构核近似配点法无 无网格法，作为一种独立的数值分析工具，发展至今已经形成了很多种方法。 除了上述几种，还有有限点法 1 2 - 1 4 j 、再生核粒子法 嘲等1 0 余种等。这些方法的 共同点就是采用定义在离散节点而非网格上的一组权函数和基函数来构造近似 场函数。这个共同点使得无网格法具有前处理简单、自适应性强和计算结果光滑 连续等优点。而这些无网格方法之间的区别在于所使用的试函数( 如移动最小二 乘法等) 和微分方程的等效形式( 如g a l e r k i n 法、p r e t r o g a l e r k i n 法等) 。表1 1 比较了几种常见的无网格方法的特点。 第一章绪论 1 1 3 无网格方法的应用 无网格方法相对于有限元等传统数值方法最突出的优点就是摆脱了对网格 的依赖性，因此在一些有限元方法求解困难的领域，如大变形【4 2 1 、裂纹扩展【4 3 45 1 、 冲击碰撞等，都取得了较大的进展。又因为计算精度高，前后处理方便等优点， 无网格方法也被应用于固体力学h 6 4 踟、流体力学4 川等力学领域甚至传热过程、 电磁场1 5 1 - 5 3 】，材料5 帅等非力学领域的数值分析中。 1 1 4 无网格方法的主要特点 建立近似场函数不需要借助网格，基函数逼近近似而非插值是无网格法与有 限元法的主要区别。采用定义在离散节点上具有紧支特性的函数来构建近似函 数，而不用定义在全域上的级数展开形式是无网格方法和经典加权残数法的主要 区别。 无网格方法具有以下优点： 1 ) 无网格法的近似函数没有网格依赖性，减少了因网格畸变而引起的困难， 适用于处理高速碰撞、断裂、塑性流动、流固耦合等涉及大变形和需要动态调熬 节点位置的各类应用问题； 2 ) 采用紧支函数的无网格方法和有限元法一样具有带状稀疏系数矩阵的特 点，适用于求解大型科学与工程问题； 3 ) 无网格方法的前处理只要节点位簧信息，不需要网格信息，容易构建复 杂三维结构； 4 ) 无网格方法的基函数可以包含能够反映待求问题特性的函数作为基函 数，适用于分析各类具有高梯度、奇异性等特殊性质的应用问题； 5 ) 无网格方法的计算结果是光滑连续的，不必再进行应力光顺化等后处理。 作为一个还在发展中的数值方法，无网格法在严格的数学论证、计算效率、 边界条件处理和大量应用实例等方面都还不成熟，有待进一步研究完善。不同于 有限元方法有确定的基函数形式，无网格方法的形函数是在计算过程中，在背景 网格或节点子域中使用高斯积分得到的，为保证计算精度，还会用到高阶高斯积 分，因此无网格方法的计算量一般大于有限元等方法。另外，用移动最小二乘法 建立无网格方法的近似函数时，涉及到系数矩阵的求逆，计算量较大，有时还会 出现矩阵不可逆的情况。如何提高无网格的计算效率是近年来的研究热点。虽然 无网格方法还不够成熟，但由于它能克服有限元分析中所面临的困难，因此在超 高速碰撞、爆炸、金属加工成形、流体动力学等领域中具有广阔的发展前景。 第一章结论 1 2 渗流问题 渗流是指多孔介质内的流体运动，它广泛存在于自然界和土木、水利工程中。 对渗流问题的研究有着重要的理论和应用价值，它不仅成为了土坝、堤防等工程 领域的一个重要课题，还为石油开采、核废料储存和生命科学技术的研究提供了 重要的理论基础。早在1 8 5 6 年，法国工程师达西就通过试验建立了多孔介质渗 流的基本定律达西线性渗流定律。2 0 世纪2 0 年代，建立了单相气体在均质介 质中的渗流方程，其渗流方程具有拉普拉斯方程的形式。2 0 世纪8 0 年代，郭尚 平、陈钟祥等深入研究了双重、三重及多重介质的渗流问题u 8 - 2 0 】。2 0 世纪9 0 年 代以来，物理化学渗流、非等温渗流、非达西非线性渗流等问题都取得了可喜的 进展。随着非线性科学理论和实验方法及计算机技术的应用，渗流理论及其在工 程中的应用研究将会得到进一步的发展。 由于问题的复杂性，工程实际主要运用数值计算方法来计算模拟。数值计算 能较好地解决油、气、水开发工程及地质工程和环境工程等领域中提出的复杂力 学问题，这与应用和研究新的数值方法是分不开的。2 0 世纪5 0 年代初，gh b r u c e 、d w p e a c e m a n 用有限差分法在计算机上求得了理想气体非线性非定常 渗流的数值解；1 9 5 9 年，d o u g l a s 研究了两相多维渗流非线性微分方程组的数值 计算方法；1 9 6 6 年，z i e n k i e w i c z 把有限元方法用于二维定常渗流计算；1 9 6 8 年， j a v e n d a l 、w i t h e r s p o o n 进一步用有限元法求解了非定常渗流问题1 。在各种数值 方法中，有限元法以其对复杂边界、备向异性、非均匀性和三维渗流等问题的强 大适应能力而得到广泛的重视 2 1 - 2 3 1 。 有限元法确定渗流自由面的方法主要分为移动网格法和固定网格法这两大 类方法。虽然移动网格法已取得了成功的经验，但是仍然存在以下问题：( 1 ) 当 初始自由面与最终自由面相差较大时，网格变形过大，会引起单元畸形，误差太 大；( 2 ) 网格变形改变了渗流域不同介质的边界；( 3 ) 渗流域内存在结构物时， 网格变形也会改变结构物的边界；( 4 ) 无法计算自由面以上区域其他物理量，因 此无法进行耦合分析。为了解决以上问题，有限元法又提出了多种固定网格法1 4 ， 如剩余流量法【37 1 ，单元渗透矩阵调整法，初流量法1 3 8 3 明和虚单元法等。剩余流量 法通过计算自由面单元内过自由表面的流量来修改各节点的势，直至过节点的流 量小于某允许值为止。其全部调整均依赖于第一次有限元计算结果，因此计算精 度较低；单元渗透矩阵调整法对跨自由面单元按复合材料单元处理。复合材料单 元渗透系数在复合面突变，单元渗透矩阵不能代表这一特性，且矩阵主对角元常 不占优，因此计算精度和计算稳定性均受到影响；初流量法在计算跨自由面单元 的节点初流量时，自由面以下高斯点不予计算，实际上把跨自由面单元当作复合 4 第一章绪论 材料单元处理，计算精度也受到影响；虚单元法在形成渗透矩阵时，不考虑移动 跨自由面单元的某些节点和自由面以上的单元，这样在迭代过程中，自由面向上 移动时，迭代将无法继续。由于渗流自由面的位置预先未知或动态变化，必须通 过迭代方法求解，用有限元法等基于网格的计算方法求解，都容易造成网格扭曲， 从而影响计算精度。 1 9 7 6 年，n i w a ，k o b a y a s h i 和f u k u i 2 4 】提出了采用边界元法确定渗流自由面 的基本迭代方法。由于边界元只在渗流边界上剖分单元，所以在确定渗流自由丽 方面较有限元具有明显优势1 2 5 。 2 0 0 4 年，薛馨、金吾根和李珏等p 副提出了用t r e f f t z 方法求解渗流自由面和 流场。该方法计算量小、收敛快，计算结果接近实验值，具有实际应用价值。 随着无网格方法的发展，人们开始试图采用无网格方法来求解带自由面的渗 流问题。目前主要用无单元g a l e r k i n 方法求解稳态渗流问题1 27 4 “。 1 3 本文的主要工作 1 ，对于渗流自由面问题的e f g 解法，提出自由面的多项式拟合方法以优化计 算结果，并通过数值算例验证渗流自由面问题e f g 解法的稳定性。 2 首次将无单元彼得罗夫一伽辽金法应用于渗流自由面的计算中，分析推导 m l p g 解法的渗流控制方程，并提出边界积分处理方法以提高计算精度。 3 通过数值算例分析比较节点分布、权函数的选取、正交多项式基的次数以及 初始自由面选取等因素对于计算精度和速度的影响，并为以后的渗流无网格 方法研究提供经验性总结。 第二章尤咐格方_ ；去的柯造 第二章无网格方法的构造 基于微分方程等效积分形式的加权残量法是求解线性和非线性微分方程近 似解的一种有效方法。无网格方法是一种独立的数值求解方法，它可以用加权残 量法来建立求解方程。紧支试函数加权残量法，即加权残量法中的试函数采用紧 支函数，可以作为无网格法的基础。只要试函数是利用离散点来建立的，那么由 紧支试函数加权残量法导出的所有方法都属于无网格法的范畴。本章首先从加权 残量法的角度出发，导出本文中用于求解渗流自由面问题的e f g 和m l p g 两种 无网格方法的基本形式。 2 1 加权残量法 加权残量法是求解微分方程的一种有效方法，许多问题可以归结为一组微分 方程及已知的边界条件。即未知函数u ( x ) 满足微分方程 a u ( x ) 】 a 1 【u ( x ) 】1 爿：【u ( x ) 】 ：o ( 在q 内)( 2 _ 1 ) 4 。，【u ( x ) 】j 和边界条件 b 【u ( x ) 】= q 【u ( x ) 】1 岛咿1l = 0 ( 在f 上) ( 2 2 ) b 。： u ( x ) 】j 其中f 是域q 的边界，x = i x ， z p 表示空间点。待求函数u ( x ) 可以是标量 场( 例如温度) ，也可以是几个变量组成的向量场( 例如位移、应变、应力等) 。 微分方程可以是单个方程，也可以是联立方程组，因此上面采用了矩阵符号。爿， 和b 是关于独立变量( 如空间坐标、时间坐标) 的微分算子。 由于要求方程( 2 1 ) 在q 域内任意点都满足，式( 2 2 ) 在边界1 1 上任意点都满 足，因此对任意函数v 和v 都有： 【v 7 a 【u ( x ) 】擒+ 【；1b u ( x ) d r = 0 ( 2 - 3 ) 反之，如果对任意函数v 和v ，式( 2 3 ) 都成立，则式( 2 1 ) 必定在域q 内任意 点满足，且式( 2 2 ) 必定在边界r 上任意点满足。函数v 和v 称为捡殓酉数，它们 分别是m ，阶和鸭阶的函数列阵。式( 2 3 ) 称为微分方程( 2 1 ) 和边界条件( 2 2 ) 的笠 煎积盆显式。 第二章光网格方法的构造 对于复杂问题，式( 2 一i ) 和( 2 2 ) 无法精确求解，只能近似求解。设u 6 ( x ) 为方 程( 2 一1 ) 和边界条件( 2 - 2 ) 的一个近似解，它可以表示为一组线性无关的已知试探函 数( 或称为形函数) n ，( x ) 的线性组合，即 u ( x ) “u h ( x ) = n 。( x ) u ，= n ( x ) u( 2 - 4 ) i = l 其中u = u l u ；，u ：r ，n ( x ) = i n 。( x ) ，n ：( x ) ，n 。( x ) 】，u 。可以是节点上的 待定函数值，也可以是待定参数。试函数的项数n 越多，近似解的精度就越高。 理论上，当项数n 趋于无穷大时，近似解将收敛于精确解。 般，近似解u 6 ( x ) 不能精确满足微分方程( 2 1 ) 和边界条件( 2 2 ) ，它们将产 生残量r 和r ： r ( x ) = a u “( x ) 】， r ( x ) = b u 6 ( x ) 】( 2 - 5 ) 为得到未知场函数u ( x ) 的最佳近似解，应以某种方式使残量r 和r 为零。 由式( 2 3 ) 可知，如果对任意检验函数v 和v ， l v 7 a i u “( x ) d q + i 。b u “( x ) 订= 0 ( 2 6 ) 都成立，则残量r 在域q 中任意点必定为零，且残量r 在边界r 上任意点 必定为零。 实际计算时，不可能对式( 2 - 6 ) 取无穷多个检验函数，通常是将检验函数v 和 v 取为一组线性无关基函数的线性组合，即 v ：杰6 ，w ，；：窆6 ，可 ，2n ( 2 7 ) ；lj = l 将式( 2 7 ) 代入式( 2 - 6 ) 中，并考虑到系数b ，的任意性，有 【w ，t a u 6 ( x ) d f 2 + 【可7b 【u ( x ) 】订= 0 ( _ ，= 1 , 2 ，) ( 2 - 8 ) 上式的意义是通过选择合适的待定参数u ，强迫残量r 和r 在加权平均意 义上等于零。在极限情况下，残量r 和r 在整个求解域内及其边界上趋于零。 式( 2 8 ) 给出了，个求解方程，用以求解近似解中的珂个待定系数u 。，从而得到原 问题的近似解。当，= 时，方程有唯一解；当， n 时，方程( 2 8 ) 是超定的，需 要用最小二乘法求解。 对式( 2 - 8 ) 进行分部积分，可以得到等效积分的弱形式，其一般表达式为 】q w 7 脚“( x ) 脚+ f e w 7 帅“( x ) 弦+ f 可b i n h ( x ) 协= 0 ( ，= 1 ，2 ，) ( 2 9 ) 其中微分算子c 、d 中包含的导数阶数要低于式( 2 1 ) 中的a ，e 、f 中包含 第二章无网格方法的构造 的导数阶数要低于式( 2 2 ) 中的b 。式( 2 9 ) 称为微分方程( 2 一1 ) 和边界条件( 2 2 ) 的弱 形式。从形式上看，弱形式对函数u “( x ) 的连续性要求降低了，但对实际问题往 往能给出较原微分方程更好的近似解，因为原始方程对解提出了过分“平滑”的 要求。 在经典的加权残量法中，试函数是定义在整个求解域上的，因此所得的系数 矩阵为满阵，存储和计算量都较大，同时这也给定义在形状复杂的区域上的微分 方程的求解带来困难。如果采用紧支函数作为加权残量法的试函数，则可以得到 紧支试函数加权残量方程，其系数矩阵一般为稀疏阵，可大大节约存储空间，减 少求解方程的时间。在紧支试函数加权残量法中选择不同的检验函数和试函数， 就得到不同的近似求解方法。无网格方法可以理解为是一种紧支的加权残量法。 下面介绍本文所采用的两种无网格方法。 2 2 无单元伽辽金( e f g ) 法的一般方程 2 2 1 由加权残量法导出无单元伽辽金方程 伽辽金法是俄i 虱工程师伽辽金提出的，它将式( 2 - 8 ) 中的试函数和检验函数取 为同一函数空间，r r = n ，即w j = n ，在边界上w ，= 一w = 一n ，。 于是式( 2 - 8 ) 可写成 l n f f a 耋n i u ，卜n 7 b 扣咖，卜屯z ，州z 邶， 对式f 2 1 0 ) 进行分部积分，得 c ( n t 瞎m x ) u f 扪一f e ( n ，瞎m x ) u ，d f = 。 ( 2 - 1 1 ) 由式( 2 4 ) ，可以定义近似解u 6 ( x ) 的变分丽6 ( x ) 为 钿“( x ) = n i 加i + n 2 钿2 + n 。抽。 ( 2 - 1 2 ) 痂，是任意的，因此式( 2 1 0 ) 又可以写成 翮胪a u 6 a a f 砸胪b u 6 a t = 0 ( 2 1 3 ) 而式( 2 - 11 ) 可以写成 c 7 ( 踟6 ) d u 6 a n f e 7 ( 钿“) f 订一v u 舻b u 6 刃= 0 ( 2 - 1 4 ) 式( 2 1 4 ) 是伽辽金提法的等效积分的弱形式。在许多力学问题中，微分算子 a 是自伴随的，因此由伽辽金法得到的求解方程的系数矩阵是对称的，这给数值 分析带来了很大的方便。当伽辽金法中的试探函数使用移动最小二乘法近似，那 第二章无网格方法的构造 么由伽辽金法就能得到无单元伽辽金( e f g ) 方程。当微分方程存在相应的泛函时， 伽辽金法与变分法往往得到相同的求解方程。下一节介绍通过变分原理导出无单 元伽辽金方程，并在下一章中将其应用于渗流自由面问题的求解。 2 2 2 由变分原理导出无单元伽辽金方程 无单元伽辽金方程的另种构造方法是，构造与控制微分方程相对应的泛函 兀，使得微分方程的求解等价于求该泛函的驻值问题。 为此，将控制微分方程写成如下形式， n u ( x ) - l ( u ) + f = 0 ( l 是线性、自伴随微分算子)( 2 - 1 5 ) 假如泛函存在并写成如下的一般形式， 兀= f ( u ，竺，) 加+ e ( u ，竺，) d r ( 2 1 6 ) i 强 ； 强 其中f 和e 是特定的算子。 将式( 2 4 ) 的近似函数代入( 2 1 6 ) 式，就得到用试探函数和待定参数表示的泛 函。令泛函取驻值，得到 m ：罢加1 + 要钿2 + + 要钿。：0 ( 2 - 1 7 ) a i i 1 加， 加。 “ 由于踟的任意性，必然有 饥 咖 a 兀 钿 铘 钿2 踟 a u 。 = 0 ( 2 1 8 ) 由这组方程可求解出参数u ，从而得出近似函数的解。 微分方程a 扣( x ) 】= o 和边界条件b 【u ( x ) 】= o ( 在r 上) 的等效伽辽金提法可以 表示为 羽= 脚7 a ( u ) d f 2 一协7 b ( u ) d f = 0 ( 2 1 9 ) 6i 其中b u ( x ) 】；l 【u ( x ) 】+ g = o ( 在r 上) l 。是微分算子 ( 2 - 2 0 ) 将( 2 1 5 ) 式和( 2 - 2 0 ) 式代入( 2 - 1 9 ) 式，可得到 9 第二章无网格方法的构造 斑7 ( l ( u ) + f 1 ( 舰一曲7 l ( u ) + g 】( 扩= o ( 2 2 1 ) 利用算予是线性、自伴随的，可以导出 c 如7 l ( u ) c 哲2 = 且丢t m 7 l ( u ) + 1 2 6 h r l ( u ) 班1 noo = 儿告砸7 l ( u ) + 吾一l ( 加) d q + 撕( u ，加) = j 弓加r l ( u ) + 丢u r a l , ( u ) d q + n ( u ，丽) n n = 艿仁u ，l ( u ) 西2 + 6 r ( u ，钿) 代入式( 2 2 1 ) ，得到 a n ( u ) = 0 兀( u ) ：【i n t l ( u ) + u 7 f d f 2 + 6 ( u ，g ) ( 2 - 2 2 ) 其中，b t ( u ，g ) 是由幻( u ，两) 和( 2 2 1 ) 式中的边界积分项两部分组成。 未知函数u ( x ) 应满足的边界条件分为强制边界条件b ，f u ( x ) 】= o ( 在r l 上) 和 自然边晃条件b ： u ( x ) 】一o ( 在r 2 上) 。对于2 m 阶的微分方程，含o m 1 阶导数 的边界条件称为强制条件，近似函数应该事先满足。含m 2 m 一1 阶导数的边界 条件称为自然边界条件，近似函数不必事先满足。 当近似函数事先满足了强制边界条件的情况下，6 ，( u ，g ) 就是在自然边界条 件下的积分项。于是( 2 2 2 ) 式中的n ( u ) 便是在近似函数事先满足强制边界条件时 的泛函。 当近似函数事先并不满足强制边界条件的时，引入l a g r a n g e 乘予旯用于施 加强制边界条件b 。【u ( x ) j - o ( 在r l 上) ，于是泛函为 兀( u ) = ( u ) + b l u 】d r ( 2 - 2 3 ) 一 从以上讨论可知，问题的等效积分的伽辽金提法等效于它的变分原理，即泛 函取驻值。满足问题的微分方程和边界条件的解也可以通过泛函取驻值得到。在 下一章中将采用这种方法构造渗流自由面问题的e f g 方程。 因此，对于一个带边界条件的微分问题，我们可以建立它的伽辽金提法的等 效积分弱形式，也可以对相应的泛函取驻值。这两种方法都能得出问题的解。当 这两个方法的试探函数使用移动最小二乘近似时，我们就能得到无单元伽辽金法 ( e f g 方法) 的形式。 1 0 第二章无网格方法的构造 2 3 局部彼得罗夫一伽辽金( p e t r o v - g a l e r k i n ) 法 由上面的介绍我们可以知道，虽然e f g 方法的试函数是紧支的，但是它所 建立的离散方程是基于微分方程的全局等效，所以需要额外的背景积分网格，因 此无单元伽辽金法并不是理想的无网格方法。 局部彼得罗夫一伽辽金法的思想是，要求控制方程的残量在r 个子域q ，内 及其边界r ，上消除，且选取不同的检验函数和试函数，即 r 广，1 l ，w y r a l 高n f ( x ) u 一棚+ f ，w ，7 b i 喜w x ) u f | d f = o ( _ ，= 1 ，2 ，2 一( 2 - 2 4 ) 式中子域q 是求解域与第j 个节点局部积分域的相交区域，而 r ：a n n r ，是位于求解域q 边界1 1 上的那部分子域边界，如图2 1 所示。对 于一个完全位于求解域内部的局部积分域，子域即为局部积分域。 甜曼一。1 r，目t 一一i - b ，旷 。、钆 ，07 。鹅。q 、：、 o 7 i ， ，q j 。- 。r j j 每n 。j j 。 c ) 。：。i ? 。、r ： q3 ：j 、。b誉。 、jf i 主二_ 一， 、碱濂一。一少 r - t 一j 图2 - 1 求解域和局部积分域 与伽辽金法相i z ，局部彼得罗夫一伽辽金法形式的积分是在每个子域q ，中 进行的，不需要引入背景网格。对( 2 - 2 4 ) 式进行分部积分，得 q w 。睡n ，蛳，卜+ e 畔7 坷卦卜+ 雨卦c 咖，卜= 。 ( _ ，= 1 , 2 ，) ( 2 - 2 5 ) 式( 2 。2 4 ) 或( 2 2 5 ) 给出了r 个求解方程，如果r = ”时，方程有唯一解，可以直 接求出待定参数u ；当， 以时，则需要用最小二乘法求解。 如果选取的w ，、丽，在q ，边界上为零，则对于一个完全位于求解域内部的 子域，式( 2 2 5 ) 中的第二项为零，成为 第二章无网格方法的构造 ，c ( w 7 ) d i n ，( x ) u 一艘= 0 ( 2 _ 2 6 ) l l = lj 如图2 - 1 所示，对于不同的边界条件，当q 与计算区域边界相交，式( 2 2 5 ) 中的第二项的积分边界由三部分组成，0 。= 铀，nr l 、1 1 ，：= 铀，n f 2 和内部边 界f 。 由于彼得罗夫一伽辽金法的检验函数和试函数取自不同的函数空间，这使得 最终得到的求解方程的系数矩阵不对称，这是它的不足之处。但同时这给施加强 制边界条件提供了方便，可以将给定的强制边界条件直接写入最终的方程而不必 顾忌破坏系数矩阵的对称性。 求解渗流问题的无单元彼得罗夫一伽辽金法的具体构造详见第三章。 2 4 移动最d - - 乘近似 移动最小二乘近似( m o v i n g l e a s t - s q u a r e a p p r o x i m a t i o n ，m l s ) 是目前无网 格法中比较常用的构造近似函数的方法之一。 设待求函数“( x ) 在计算点x 领域q ，内的近似函数为 州 一 “6 ( 薰，x ) = p ，( i ) a ，( x ) = p 。( i ) a ( x )( 2 - 2 7 ) - 1 其中x = 【x ，弘= 】7 是计算点x 的领域q 。内各点的坐标。 a ( x ) = b ( x ) ，吒( x ) ，( x ) r ，a i ( x ) 是待定系数，与普通的插值方法不同，这些系 数随x 的变化而变化，b ( x ) 是基函数，p 7 ( x ) = 吼( x ) ，p - z ( x ) ，；( x ) 】，m 是基函数的 个数。通常我们使用单项式作为基函数。 对于二维空间，单项式基函数为 线性基：p ( i ) = 1 x y r m = 3 二次基：p ( x ) = 【l xyx 2 x yy 2 埘= 6 系数口，( x ) 的选取必须使得近似函数“( i ，x ) 是待求函数u ( x ) 在某种最小二 乘意义下的最佳近似。设计算域q 内包含个节点x ，( ，= 1 , 2 ，) ，在每个 节点处定义一个非负的权函数w ，( x ) ，该权函数只在节点x ，周围的一个有限区域 q ，中大于零，而在该区域外为零。区域q ，称为权函数w a x ) 的支持域。另外，x 越接近于x ，权函数w ，( x ) 的值越大，在x ，处达到最大，它是x 与x ，之间距离 的函数，因此可写成w ，( x ) = w ( x x ，) 的形式。 设计算点x 的领域q 。包括野个节点，n 个节点上的函数值为“，。定义近似 函数“6 ( x ) 在这些节点x = x ，上的误差的加权平方和 第二章光删格方法的构造 以x ) ，( x ) = w ，( x x m “6 ( x ，x ) 一“，】2 写成矩阵形式， p l ( 1 i ) p l ( 1 2 ) p 2 ( x 1 ) p 2 ( x 2 ) p m ( x 1 ) p m ( x 2 ) p l ( x 。) p 2 ( x 。) p 。( x 。) 记：u = p = p 1 ( x 1 ) p 1 ( x 2 ) p 2 ( x 1 ) p 2 ( x 2 ) nm = z w ，( x x m p 小加，( x ) - u 1 = 1 j = 1 嘶( x ) 盯2 ( x ) a m ( x ) p 。( x 1 ) p 。( x 2 ) p l ( x 。) p 2 ( x 。) p 。( x 。) i l ，w ( x ) ： w 1 ( x ) o 0 w 2 ( x ) jl 0 p l ( x 1 ) p 2 ( x 1 ) p 。( x 1 ) p 1 ( x 2 ) p 2 ( x 2 ) p 。( x 2 ) p l ( x 。) p 2 ( x 。) p 。( x 。) 贝i ，以x ) = ( p - a ( x ) 一u ) 7 w ( x ) ( p a ( x ) 一u ) “l “2 ： “” 0 ： 0 0 w a x ) = a t ( x ) p 7 w ( x ) p a ( x ) 一2 a 7 ( x ) p 7 w ( x ) u + u t w ( x ) u 记，a ( x ) = p x w ( x ) p b ( x ) = p t w ( x ) 则，( x ) = a t ( x ) a ( x ) a ( x ) 一2 a r ( x ) b ( x ) u + u r w ( x ) u r 2 2 8 ) w 。( x ) f 2 2 9 ) ( 2 - 3 0 ) ( 2 3 1 ) 系数a ( x ) 的选取要使误差的加权平方和，( x ) 取极小，即差等= 。，于是得 a ( x ) = a ( x ) 。1 b ( x ) ur 2 3 2 ) n 蜥屹 rll ； q 吒 第二章无网格方法的构造 代回到( 2 2 7 ) 式，得到近似函数的表达式， “( x ) = p ( x ) a ( x ) _ 1 b ( x ) u = n u = m ( x ) “， ( 2 3 3 ) i f f i l 其中n 为形函数矩阵， n = p r ( x ) a ( x ) 一1b ( x )( 2 3 4 ) 在实际应用中，需用到形函数关于坐标的偏导数， n 。= p t d ( x ) a ( x ) _ 1 b ( x ) + p t ( x ) a ( x ) _ 1 ，b ( x ) + p t ( x ) a ( x ) - 1 b ( x ) ， ( 2 3 5 ) 其中a ( x ) 一，= 一a ( x ) 山1a ( x ) ，a ( x ) 。( 2 3 6 ) 形函数的构造涉及到对a 阵的求逆，因此要求a 是非奇异阵。权函数的选 取、影响域的大小、节点的分布情况等都会对a 矩阵的性质产生影响，从而影 响到形函数的形式和计算结果。 2 5 紧支近似函数 紧支近似函数是定义在局部区域上的函数，它只在其支持域中有定义，在其 支持域外为零。当我们用移动最小二乘法来构造形函数时，权函数的紧支性使得 形函数也具有紧支性。在二维问题中，支持域一般取为圆形域或矩形域，它与求 解域q 相比，是一个很小的区域，如图2 2 所示。这样，函数u ( x ) 在x 点的值就 可以用该点支持域内的一组离散节点x ，( = 1 , 2 ， i x ) 上的紧支函数，( x ) 的线 性组合来近似，即 h ， u ( x ) * u h ( x ) = n ，( x ) “，= n u ( 2 - 3 7 ) ，一i 式中n x 是x 点支持域内的节点数，“，是函数u ( x ) 在节点x ，处的值，即 u ( x ，) = “，。 ( a ) 圆形支持域 第二章无网格方法的构造 ( b ) 矩形支持域 图2 - 2 节点的支持域 不同于有限元方法的插值形函数，无网格方法的由移动最小二乘法构成的形 函数，( x ) 是一种加权近似，所以一般不具有k r o n e c k e r - 8 性质，即“6 f x ，) “，。 而且对于某个节点x ，它的场函数“( x ，) 不是由这个节点上的参数“，确定，而 是同节点x ，支持域内的所有节点上的参数有关，这就给施加本征边界条件带来 一些困难。另外，各节点对应的紧支函数的支持域可以是相互重叠的，如图2 2 所示。支持域的大小是影响无网格方法计算结果的重要因素之一，但还没有一个 统一的选取规则，因此成为目前无网格方法研究的一个热点。 2 6 权函数的选取 权函数的选取在无网格方法中起到很重要的作用。文献7 1 分析表明无网格方 法的计算精度与节点权函数及其影响域大小有关。总体上，权函数( x ) = w ( x x ，) 应该遵循这些原则： 非负性，即w ( x x ，) 0 ； 单调性，即w ，在节点x ，处的值最大，当x 离x ，越远，w 。的值就越小。 紧支性，即当j x - - x ，l = d m ，时，w ( x x ，) = 0 。为紧支域的半径。 定义r = i x x i i j d m l ，则权函数可写为w ( r ) 的形式。 在无网格计算中，节点权函数影响域的大小通常根据节点分布密度