（应用数学专业论文）地下水模拟中的无网格局部径向基点插值法.pdf

!f蠢 -繁爹甏鼍气l p t - 。 r 】 成果。论文中除特别加以标注和 ，其他同志的研究成果对本人的 局议 学位论文版权的使用授权书 本学位论文作者完全了解辽宁师范大学有关保留、使用学位论文的规定，及学校有 权保留并向国家有关部门或机构送交复印件或磁盘，允许论文被查阅和借阅。本文授权 辽宁师范大学，可以将学位论文的全部或部分内容编入有关数据库并进行检索，可以采 用影印、缩印或扫描等复制手段保存、汇编学位论文，并且本人电子文档的内容和纸质 论文的内容相一致。 保密的学位论文在解密后使用本授权书。 学位论文作者签名：白亟 指导教师签名：! 娄查塑 学位论文作者签名：囱么指导教师签名：f 堑塑型 签名日期： 沙【1 年莎月陟日 ，_ p j - 7 ： 辽宁师范大学硕士学位论文 摘要 目前无网格方法已经成为国内外研究的热点。无网格局部径向基点插值法( l r p i m ) 作为无网格方法中的一种，它采用局部子域上的加权残量形式，允许权函数取自不同空 间，由于积分在局部子域上进行，不需要额外的背景网格，因此是一种真正的无网格法。 无网格局部径向基点插值法是一种新兴的数值模拟方法，在水文地质领域的应用国内外 很少有报道。 本论文通过对求解微分方程的加权残量法和点插值形函数的分析，从局部 p e t r o v g a l e r k i n 法和径向基点插值形函数出发，构造了地下水二维稳定流和不稳定流的 无网格局部径向基点插值法，详细论述了二维稳定流和不稳定流无网格l r p i m 的算法 步骤，绘制了程序结构流程图。基于以上研究，将无网格局部径向基点插值法应用于地 下水二维稳定流和不稳定流实际问题的水头计算，利用m a a b 语言编程进行模拟得到 了满意的结果。结果表明，无网格局部径向基点插值法不但可以提高计算的精度，而且 具有很强的适应性，可以为地下水工程设计提供参考。 本论文对无网格局部径向基点插值法研究所得的成果，将有助于无网格法的进一步 研究。 关键词：加权残量法；点插值；局部径向基点插值法；地下水 ， 地下水模拟中的无网格局部径向基点插值法 m e s h l e s sl o c a lr a d i a lp o i n ti n t e r p o l a t i o nm e t h o di nt h e g r o u n d w a t e rs i m u l a t i o n a b s tr a c t a tp r e n s e n t ，m e s h l e s sm e t h o dh a sb e c o m er e s e a r c hh o t s p o ta th o m ea n da b r o a di n s c i e n t i f i cc o m p u t i n g a sak i n do fm e s h l e s sm e t h o d ，m e s h l e s sl o c a lr a d i a lp o i n ti n t e r p o l a t i o n m e t h o du s e sw e i g h t e dr e s i d u a lf o r mo nl o c a ls u b - d o m a i n a l l o w st h ew e i g h tf u n c t i o nf r o mt h e d i f f e r e n ts p a c e ，d o e sn o tn e e da d d i t i o n a lb a c k g r o u n dg r i ds i n c e t h ei n t e g r a lo nt h el o c a l s u b d a m a i n i ti sat r u em e s h l e s sm e t h o d m e s h l e s sl o c a lr a d i a lp o i n ti n t e r p o l a t i o nm e t h o di sa n e wn u m e r i c a ls i m u l a t i o nm e t h o d ，h a sb e e nr e p o r t e dr a r e l ya th o m ea n da b r o a di nt h ef i e l do f h y a r o g e o l o g y o nt h ea n a l y s i so ft h ew e i g h t e dr e s i d u a lm e t h o di ns o l v i n gd i f f e r e n t i a le q u a t i o n sa n d p o i n ti n t e r p o l a t i o nm e t h o d ，t h i sp a p e rc o n s t r u c t sm e s h l e s sl o c a lr a d i a lp o i n ti n t e r p o l a t i o n m e t h o do ft w o - d i m e n s i o n a ls t e a d ya n du n a b l eg r o u n d w a t e rf l o wf r o mt h el o c a l p e t r o v - g a l e r k i nm e t h o da n dt h er a d i a lp o i n ti n t e r p o l a t i o ns h a p ef u n c t i o ns t a r t i n g ，d e s i g n st h e p r o c e d u r ef l o wc h a r to ft h es t r u c t u r eo fm e s h l e s sl i 冲i mm e t h o da n d d r a w su pp r o c e d u r e so f c a l c u l a t i n gt w o d i m e n s i o n a lp r o b l e m b a s e do nt h e s es t u d i e s 。t l em e s h l e s sl o c a lr a d i a lp o i n t i n t e r p o l a t i o nm e t h o da r ea p p l i c a t e di nt h et w o d i m e n s i o n a ls t a b l ea n du n s t a b l eg r o u n d w a t e r p r o b l e m s ，g e ts a t i s f a c t o r yr e s u l t sb ys i m u l a t i o nu s i n gm a t l a b e x a m p l e ss h o wt h a tt h el o c a l r a d i a lp o i n ti n t e r p o l a t i o nm e t h o dn o to n l yc a ni m p r o v et h ec a l c u l a t i o na c c u r a c y , b u ta l s oh a s s t r o n ga d a p t a b i l i t y ，c a np r o v i d ear e f e r e n c ef o rt h ee n g i n e e r i n gd e s i g no fg r o u n d w a t e r i nt h i sp a p e r ，am e s h l e s sl o c a lr a d i a lp o i n ti n t e r p o l a t i o nm e t h o d sr e s u l t sw i l lp r o m o t et h e d e v e l o p m e n ta n dr e s e a r c ho ft h em e s h l e s sm e t h o d k e yw o r d s ：w e i g h t e d r e s i d u a lm e t h o d ；p o i n ti n t e r p o l a t i o n ；l o c a lr a d i a lp o i n ti n t e r p o l a t i o n m e t h o d ；g r o u n d w a t e r ，i f i 2 1 - 3 无网格l r p i m 的程序设计及流程图1 5 2 2 二维不稳定流的无网格l r p i m ：1 8 3 采用无网格l r p i m 的地下水问题2 0 3 1 地下水中的二维稳定流问题2 0 3 2 地下水中的二维不稳定流问题2 2 结论2 4 参考文献2 5 攻读硕士学位期间发表学术论文情况2 7 致谢2 8 i1 3 3 4 4 5 5 8 2 2 2 4 1 1 1 l _ - f 上 f 辽宁师范大学硕士学位论文 引言 2 0 世纪5 0 年代兴起的有限元法【1 1 ( f i n i t em e t h o d ，简称f e m ) 是数值方法领域的 重要里程碑之一。它是从加权残量法或变分原理出发，通过单元或网格剖分，把微分方 程的边值问题转化为一组等价的线性代数方程来进行求解。它的另一个重要特点是利用 在每一个单元内假设的近似函数分片地表示问题域上待求的未知函数。在对一个问题的 分析中，近似函数通常由未知函数及其导数在单元节点上的数值函数来表示，从而使一 个连续的无限未知量问题变成离散的有限未知量问题。求解出这些未知量，就可通过插 值函数计算出各个单元内场函数的近似，从而得到整个问题域上的近似解。目前在许多 实际工程问题中，有限元法虽然已经得到广泛的运用，但是f e m 也有其数值方法上所 固有的局限性：形成f e m 网格时的计算成本高；应力精度低：自适应分析困难；对某 些问题的分析有局限性，产生这些问题的原因是系统方程形成过程利用了单元或网格剖 分，为了克服这些困难，摆脱单元或网格的信息，无网格法的概念便慢慢形成了： 无网格法( m e s h l e s sm e t h o d ) ，也称为无单元法( e l e m e n t f r e em e t h o d ) ，它不需要 利用预先定义的网格或单元信息进行域的离散，而是利用一组散布在研究域中及研究域 边界上的节点来表示。它们不需要网格或单元的初始划分和重构，即不需要事先定义的 节点连接用于构造变量未知函数的近似表达式。无网格法产生于二十多年以前，但当时 发展很缓慢，而且应用的领域也非常有限。其中美国西北大学工程力学系的b e l y t s c h k 等人和西班牙数值分析中心o n a t e 等人在这方面的研究最为活跃，作了大量的工作。最 近几年，随着越来越多的学者投入到这方面的研究中，各式各样形式的无网格法不断涌 现出来，其中主要的无网格法有d y k a ，s w e g l e 等提出的光滑质点流体动力学方法f 2 1 ( s p h ) ：b e l y t s e h k o 等人在1 9 9 4 年提出的无网格伽辽金法f 3 】( e f g m ) ；a f l u r i 等提出 的无网格局部伽辽金法【4 j ( m l p g m ) 和局部边界积分方程无网格法f 5 1 ( m l b i e m ) ；移 动最d - - 乘积分核法【6 】( m l s r k m ) ；1 9 9 5 年美国计算数学学者b a b u s k a 和他的学生 m e l e n k 提出了单位分解法【7 】( p u m ) ：1 9 9 6 年西班牙数值分析中心的o n a t e 和i d e l s o h n 等人提出的有限点法【8 】( f p m ) ；1 9 9 8 年s u k u m a r 提出的自然单元法( 9 1 ( n e m ) ：有限 覆盖无单元法【l o l ( f c e f m ) ；2 0 0 1 年新加坡国立大学的g u 和l i u 提出的点插值法【l l 川 ( p i m ) ；径向基点插值法 1 3 , 1 4 i ( r p i m ) 。 目前无网格法作为一种新的有效数值计算方法已经取得良好的应用效果，并将成为 一种强有力的数值分析工具来解决各种各样实际复杂的工程问题及科学研究问题。 本文围绕无网格法展开，通过阅读大量有关无网格法的文献和资料，对无网格方法 地下水模拟中的无网格局部径向基点插值法 有了一定的认识。在此基础上，构造了地下水数值模拟的局部径向基点插值法，并且将 这种方法应用到求解二维地下水稳定和非稳定流的数值模拟问题中。主要内容概括如 下： 引言介绍了无网格方法的提出和发展。 第一章总结了求解微分方程的加权残量法，特别列出了两种特别的加权残量法， 同时详细叙述了多项式基点插值形函数和径向基点插值形函数的构造过程。 第二章针对径向基点插值形函数和局部p e t r o v g a l e r k i n 法，应用加权残量法推导 离散系统方程，建立了二维地下水稳定和不稳定流的无网格局部径向基点插值法 ( l i 强i m ) 方程，详细论述了无网格l r p i m 法的程序设计，绘制了程序结构流程图， 第三章基于以上研究，将基于径向基点插值法应用于二维地下水稳定和不稳定流 得具体数学模型中，给出了该类模型的求解算法，就实际问题利用m a t l a b 语言编制了 算法的实现程序，通过计算机上的计算对所得的求解算法进行了验证。 结论对全文进行了总结，并对无网格法的发展趋势进行展望。 辽宁师范大学硕士学位论文 1 预备知识 1 1 加权残量法 加权残量法是一种建立各类科学问题和工程问题有力的数学工具，是一种用于求解 微分方程近似解的有效方法，是一种将微分方程转化为积分方程的方法，大量数值方法 均基于加权残量法。 求解给定边界条件与初始条件的微分方程定解问题，考虑以下的微分方程 燃：置套霆裟内 n - ， 【b ( j 1 1 ) + g = o ，在边界r 上 、7 其中彳为微分算子，它定义了一种作用于未知函数h 的运算，b 为边界条件的微分 算子。 大多数问题是以常微分方程或偏微分方程表述的，人们总是希望得到这些方程的精 确解，不幸的是，对于大多数复杂问题常常无法求出其精确解，只能通过近似方法来求 解。具体步骤首先是将函数h ( x ) 近似表示为一组已知函数m ( x ) 的线性组合，即 h ( x ) h ( x ) = 口j n a x ) = n a ( 1 2 ) 其中m 是第棚i i 基函数或形函数，a ；是第i 项基函数或形函数的系数，玎是基函数 或形函数的项数。当项数以越大，近似解的精度就越高；当项数刀趋于无穷大时，近似 解h ( 曲将收敛于其精确解。 令在问题域q 上的残量函数和边界上的边界条件残量函数分别为 = a ( h ) + f ( 1 3 ) r r = b ( ) + g( 1 4 ) 许多实际问题，通常很难得到它们的精确解，所以r n ，r r 基本不为零。为了寻求其 最佳近似解，一般方法是使残量函数尽可能小，故我们采用令在平均意义上的残量加权 积分为零，即 w , p b d n + i v _ f r r d l = 0 ( 1 5 ) 一 f 式中i = 1 ，2 ，n ，而，k 为已知权函数，权函数的选取可以相等可以不等。 加权残量法中权函数的选取不同，其结果也会有所影响。通过选取不同的权函数将 会形成不同的近似方案，如：配点法，子域法，最小二乘法，力矩法，g a l e r k i n 法等。 下面介绍两种特殊类型的加权残量法。 1 1 1g ale r kin 法 g a l e r k i n 法可以被 数，即 式( 1 5 ) 可以变为 式中i = 1 ，2 ，n ， g a l e r k i n 法在各种数值方法中均得到广泛应用，但是它需要进行背景网格，这就增 加了计算难度，为避免这个问题，我们采用下面的加权残量法 1 1 2 局部p e t r o v g a i e r k i n 法 在局部p e t r o v - g a l e r k i n 法中，要求方程的残量在子域q ，内及其边界r ，上加权积分 为零，且试函数和权函数不同，即 l 形 4 ( n a ) + 厂f q + i k b ( n a ) + g 订= 0 ( 1 8 ) ( 式中j = 1 ，2 ，- ，上式给出了，个方程。如果，= n ，直接求解方程；如果r 拧， 需要用最d x - - 乘法求解。 注意，权函数的选取具有一定的随意性，只要满足实际条件即可，以下列出最常用 的几个权函数 三次样条函数( w 1 ) l2 3 - 4 r t 2 + 4 r t 3弓o 5 形( x ) = 4 3 - 4 r f + 4 ，；2 - 4 3 r f 3 0 5 1 l o p 1 它具有二阶连续性。 四次样条函数( w 2 ) 彬c x ，= 6 ，；2 + 8 巧3 3 ，；4i 三： 它具有三阶连续性。 五次样条函数 4 ( 1 9 ) ( 1 1 0 ) 辽宁师范大学硕士学位论文 形c x ，= 三一1 。2 + 2 。，；3 1 5 ，：4 + 4 5i 量： c t ， 它具有四阶连续性。 指数函数( w 3 ) 僻f 埘2焉 坳 式中口为形状参数，其中的 ，；：堕( 1 1 3 ) 式中喀= x - - x 。i 为节点x ，与计算点x 之间的距离，为权函数支持域的尺寸，这里采 用的式圆形支持域。如果是矩形支持域，现定义 彬( x ) = 虼( x ) ( x ) = 矽( 么) 矿( 名) ( 1 1 4 ) 其中= 掣，勺= 掣 而如，丸分别为沿x 和y 方向的支持域尺寸。 1 2 点插值方法 点插值方法【1 5 1 是由g r “u 和他的同事们提出的一种无网格插值技术，它是一种实 用的生成无网格形函数的方法。具体方法如下，问题域q 中的一未知函数 ( x ) 由一组节 点表示，在一计算点x 处未知函数j l z ( x ) 可以近似表示为 j i l ( x ) = p 。( x ) q ( 1 1 5 ) 式中的朋为基函数的数量，p ，( x ) 为基函数，a i 为基函数的系数。 利用多项式基函数和径向基函数现已形成了两种形函数，即多项式基点插值法形函 数和径向基点插值形函数。 1 2 1 多项式基点插值法( p i m ) 形函数 应用最早的插值方法之一：多项式插值基函数，已经被广泛运用于各种数值计算方 法中。具体插值过程如下，设在问题域q 中的一个未知函数 ( x ) 可由一组节点表示，在 计算点x 处未知函数 ( x ) 可以近似为 5 地下水模拟中的无网格局部径向基点插值法 ( x ) 石( x ) = 乃( x ) 巳( x ) = p r ( x ) a ( x ) ( 1 1 6 ) j = t 其中a ( x ) 为基函数的系数向量，可表示为a t ( x ) = 【口。( x ) a 2 ( x ) a 。( x ) 】，p r ( x ) 为多 项式基函数( p i m ) ，可以采用p a s c a l 三角形确定，并且尽量选用完备基函数。对于1 维( 1 d ) 和2 维( 2 d ) 空间，一次基函数分别为 p ( x ) = 【1 石rm = 2 ，p = 1 ( 1 d ) ( 1 1 7 ) p ( x ) = 【1 石y fm = 3 ，p = 1 ( 2 d ) ( 1 1 8 ) 二次基函数分别为 p ( x ) ：1 工工2 rm ：3 ，p ：2 ( 1 d ) p ( x ) ：ixy 工2x yy 2 rm = 6 ，p = 2 ( 2 d ) 完备的p 阶多项式可以表示为以下一般形式为 p ( x ) ：f 1xx 2 x p ix p r ( 1 d ) p ( x ) ：i 工yx 2x yy 2 石py p r( 2 d ) l ” i太弋， n ， 二 ， 三 r 曩力!，j 石3 砖x 4 撅x s 歹安久 图1 1 二维问题的p a s c a l 三角形 ( 1 1 9 ) ( 1 2 0 ) ( 1 2 1 ) ( 1 2 2 ) 为了求解基函数的系数a ；，需要在计算点x 处形成包含n 个节点的支持域。在传统 点插值法中，基函数个数m 和局部支持域中的节点数总是相等的，即拧= m 。式( 1 1 6 ) 中的基函数系数口，可通过令式( 1 1 6 ) q p 的h ( x ) 等于该，1 个节点上的函数值来求解，即 6 辽宁师范大学硕士学位论文 啊= 口f p ( x 1 ) = 口l + 口2 x l + 口3 y l + + 口用p 舶( x i ) i = 1 j l z 2 2 口j p ( x 2 ) = 口l + a 2 屯+ a 3 y 2 + 响朋p 所( x 2 ( 1 2 3 ) h 。= q p ( x 。) = 口l + 口2 + 口3 y 。+ + 口m p 埘( x 。) 上式用矩阵形式则可表示如下 u = l a 式中节点的函数值用向量可表示为 u r = j j l lj l l 2j i l 3 7 l 。】 多项式的矩阵形式为 壁w 五y l x 2y 2 by 3 x ny r x l y l p 。( x 1 ) x 2 y 2 p 。( x 2 ) 毛乃p 。( x 3 ) y 。p 。( x 。) ( 1 2 4 ) ( 1 2 5 ) ( 1 2 6 ) 基函数的系数用向量表示为 a 1 = a la 2a 3 口。】( 1 2 7 ) 由于在点插值法中，基函数个数和局部支持域中的节点数相等，即，l = m ，所以p 埘 是维数等于( 1 1 或m x m ) 阶的方阵。 求解式( 1 2 4 ) 可以得出 a = l 1 u ( 1 2 8 ) 只有当p 。- 1 存在，上面的方程才会有解；否则求不出基函数系数。但是这个条件并 不是总是成立的，这个问题以后可以讨论。 把式( 1 2 8 ) 代入式( 1 1 6 ) q b ，可以得到 五( x ) = p r ( x ) k 1 u = 。h f = 7 ( x ) u ( 1 2 9 ) f = i 其中m ( x ) 为形函数向量，表示为 m 1 ( x ) = p 1 ( x ) k 1 = i ( x ) 2 ( x ) 。( x ) 】( 1 3 0 ) 因为该点插值法形函数为多项式的形式，故上述形函数的导数很容易求得。 p i m 形函数的一阶导数可以得出为 7 地下水模拟中的无网格局部径向基点插值法 7 ( x ) = p i m 形函数的二阶导数可以得出为 西”( x ) = 。( x ) 2 ( x ) 。7 ( x ) 西。”( x ) 2 ( x ) 。”( x ) ：掣p = - 叙 ” ( 1 3 1 ) = 掣k 1 ( 1 3 2 ) 1 2 2 径向基点插值( r p i m ) 形函数 采用多项式基函数可能会引起奇异性问题，为了克服这种奇异性问题，我们采用径 向基函数( i 出f ) 代替多项式基函数。具体方法如下，添加多项式的i 冲蹦插值形式可 以表示为 | l ( x ) = 足( x ) q + p j ( x ) d j = r r ( x ) e + p r ( x ) d ( 1 3 3 ) l f f i l 1 = 1 其中r i ( x ) 为径向基函数，n 为径向基函数的个数，p ，( x ) 为多项式基函数，m 为多 项式基函数的个数。如果m = 0 ，则该式为单纯的r b f s ；否则为添加了m 个多项式基函 数的r b f 。q 和d ，分别为r b f 和p i m 的系数。 常用的径向基函数有 一 g a u s s 函数：r f ( x ) = e 8 2 m u l t i q u a d r i c 函数： l r j ( x ) = ( 2 + c 2 ) 2 l r j ( x ) = ( 2 + c 2 ) 2 t h i np l a t es p l i n e 函数： r ，( x ) = ，；2 l n r r f ( x ) = 2 n 1 其中万，c ，后是已知函数。 这里r 可取为 p i m 的矩阵可以表示为 p l = r l ( ) r 2 ( ，：1 ) r 。( ，：i ) 1 x t y l 而y l 1 x 2y 2x 2 y 2 1 x 3y 3x 3 y 3 1 毛y 。x y 。 p m ( x 1 ) p 。( x 2 ) p 。( x 3 ) p 。( x 。) r b f 的系数用向量可以表示为 二r = 【c lc 2 c n 】 p i m 的系数用向量可以表示为 d r = p 。d ：d 。】 其中r ，( 噍) 的气定义为 r k = ( 石t t ) 2 + ( j ，t y 1 ) 2 对于2 d 问题 咯= 4 ( x i 一工j ) 2 + ( j ，t - y f ) 2 + ( z 七一z f ) 2对于3 d 问题 然而式( 1 3 4 ) d p 有n + m 个未知数，而上面只有n 个方程，为了求解方程， 加r n 个方程，故可使用下面m 个约束条件来添加方程 p 砖加，= c = o ，= 1 , 2 ，m 9 ( 1 3 7 ) ( 1 3 8 ) ( 1 3 9 ) 需要再添 ( 1 4 0 ) 地下水模拟中的无网格局部径向基点插值法 到 u “= 苫 = 0 三 = s a 。 = il = i：10 i =， 其中 钌= i 苫l = 阪j 1 2 以 。】r s = 曙0 a 。：三 ：【q 乞c ；匾吐丸r ( 1 4 1 ) ( 1 4 2 ) ( 1 4 3 ) ( 1 4 4 ) 由于矩阵r 。是对称的，所以我们很容易得出矩阵s 也是对称的。求解( 1 4 1 ) p - 以得 a l = 阱旷哂 式( 1 3 4 ) 可以重写为 ( x ) = r r ( x ) c + p r ( x ) d = r r ( x ) p r ( x ) i 三i 把式( 1 4 5 ) 代入式( 1 4 6 ) n - - 以得出 五( x ) = r 7 ( x ) p r ( x ) s 一1 d = 西r ( x ) 寸 其中r p 。i m 形函数可表示为 面r ( x ) = ir r ( x ) p r ( x ) i s q = ，( x ) ：( x ) 。( x ) 川( x ) 。+ 。( x ) 】 相应的r p i m 形函数的导数可以求出为 塑塑：l o r r ( x ) o p t ( x ) i s - o x l o x 次j ：垫盟o a h ( x ) 望盟翌型盟塑幽盟i l o x叙 o xa x缸 j 最终对应n 个节点的r p i m 形函数m ( x ) 可表示为 m ( x ) = 【l ( x ) 2 ( x ) 。( x ) r r p i m 形函数m ( x ) 的一阶偏导数可求得为 1 0 ( 1 4 5 ) ( 1 4 6 ) ( 1 4 7 ) ( 1 4 8 ) ( 1 4 9 ) ( 1 5 0 ) 辽宁师范大学硕士学位论文 i 冲蹦形函数( x ) 的二阶偏导数可求得为 ( 1 5 3 ) ( 1 5 4 ) ( 1 5 5 ) 、，、， 1 1 二 i ， 5 _11 ，l，l 下lj甲lj 曲一 曲一 燮叙州一匆 曲一 曲一 嘲百嘲匆 燮叙 = 掣 妨一 盟匆 = 型匆 a 一 对一 对一 对一 (一2 (一2 (一y 立叙 竺可q 一向嗍可嗍可嗍百 对一 对一 对一 (一2 (一2 (一y 姒一帮 燮妒塑却 铲一 铲一 伊一 曲一 对一 垒叙 垒匆 铲一 铲一 。l。l = = ) 一 ) 一 伍 一叙 一匆 a a 一 盟如 a 一 = 熊埘 a 一 地下水模拟中的无网格局部径向基点插值法 2 无网格局部径向基点插值法基本原理 2 1 二维稳定流的无网格l r p i m 2 1 1 二维稳定流的无网格l r p l m 公式 考虑如下地下水二维稳定流问题 v ( r v h ) + 占= 0 ，( 石，y ) q( 2 1 ) hr i = h ( 2 2 ) 娑i = q ( 2 3 )_ l 2 u j ) o n l r , 其中r 为渗透系数，h 为水头函数，为越流补给强度，譬为补给量，q 为研究区 域，r 为本质边界，r ，为自然边界，而 v ( r v h ) ：i a 【z ，i o h ) + i a 【z ，i a h ) ( 2 4 ) u 磊u 玉 o yu y 在整个问题域及边界上分布甩个节点，对于每一节点i ( i = 1 ，2 ，疗) 采用局部加权残 量法，将得到该节点的局部弱式方程，具体形式如下 l 形( v ( 舟 ) + e ) d f 2 = 0 ( 2 5 ) 式中的矿为第i 个节点的权函数或检验函数。 为了计算方便，一般选择形状规则的子域，如对于二维问题可以选圆形、矩形等形 状规则的子域，而对于三维问题可以选择长方体、球体等形状规则的子域。 利用格林公式对式( 2 5 ) 左边第一项进行分部积分可得 l 彤( v ( 四l z ) ) d q = l 。矿( 昙( 丁瓦o h ) + 南( z 等) ) d q ( 2 6 ) = f 。( t 鬈彬n x + t u y ) a t - ，( 丁鬈彬，+ 丁筹彬。y ) d q 式中的n 刀，表示边界上单位外法线向量的坐标。x 将式( 2 6 ) 代入式( 2 5 ) 可得到如下的局部弱式方程 ( r 尝谚n x + t 等玩订一t 仃差彬，+ r 考彬，y 一矿占) a n = 。 ( 2 7 ) 1 2 辽宁师范大学硕士学位论文 局部积分域q 。的边界由三部分组成，即= u 。u l ，其中： l ，为积分域的内部边界：。为与积分域相交的自然边界； 1 q 1 8 为与积分域相交的本 质边界。 所以式( 2 7 ) 可改写为 工一( r 警谚以+ r 筹矿哆) 订+ l 仃警彬以+ 丁考谚b ) 订 + l 仃差谚以+ 丁筹形：) d r l 仃芸暧。+ r 考哦，二彬占) d q ：。 亿8 我们可以合理选取权函数或检验函数使其在r 上为零，这样就可以简化局部弱式形式。 为了得到系统方程，我们对水头函数采用r p i m 形函数近似表示为 h 6 = ，| j l ， ( 2 9 ) j_-0 j j 、 式中，为r p i m 形函数，聊为计算点的支持域内的节点数，它对于总结点数还有 总体编号。 将( 2 9 ) 式代入( 2 8 ) 式可得到如下离散系统方程 荟哆l 仃鲁彬攸+ ? 等矿r + 芸嘭l 仃等谚以+ r 詈彬咿r + 芸乃l 仃- 呜l - 彬吃+ r 等啊脚一芸乃( 丁i o c l ) ： 卞百o c l ) j 脚( 2 1 。) + l 彬枷= o 将( 2 1 0 ) 改写为矩阵形式如下 a h = b ( 2 11 ) 其中 以= l 仃- 抽l - 7 谚以+ r 等吼加+ l ( r 等呱+ r 百o f l i 啊 归 2 ) + l ( r 鲁呱+ r 鲁嘱加一l 仃等死+ r 等彬d q 占，= 一上。矿耐q ( 2 1 3 ) 我们利用高斯积分计算( 2 1 2 ) 和( 2 1 3 ) d ? 的积分。注意a 和b 的表达式中既有面积积 分又有曲线积分。 从而可解出h i 办：，h 。，代入( 2 9 ) 式即可得到问题的数值解。 地 2 1 2 施加本质边界条件 设本质边界上的第f 个节点的水头值为 鸟= 7 i( 2 1 4 ) 利用罚函数法，只是将矩阵a 中的对角元素以变为 。 4j 口。呜( 2 1 5 ) 其中口为惩罚因子，它的值较矩阵a 中的元素大许多，其它元素不变。矩阵a 则会 变为 a = 4 l 4 ( h ) 4 f - i ) j 4 班 4 f + 1 ) f 4 。 4 ， 4 l 口4 4 m ) f 4 ， 右端向量b 中的相应分量e 变为 卟始1 则右端向量变为 b = e a 4 , h , 4 ( f + 1 ) 4 h 肌1 ) 4 ( ) 4 x 件i ) 4 ( f + 1 ) 12j 1 j 卜第i 行 4 。 4 h ) 。 以 4 m ) 。 4 l 。 ( 2 1 6 ) ( 2 1 7 ) 罚函数法的优点是矩阵中只有两处发生变动，所以算法比较简单。但是它并不是精 确满足本质边界条件，而只能近似满足。这里惩罚因子的选取会大大影响到精确度的高 1 4 ) ) h d h 畎俨似；缸； b = 局 ： 二 红 ： e 卜第i 行 惩罚因子越大，施加的 的问题。所以在许多实 文献中取口= 2 0 1 0 3 。 ( 2 1 8 ) ( 2 1 9 ) 而得到新的矩阵a 为 ( 2 2 0 ) ( 2 2 1 ) 2 1 3 无网格l r p i m 的程序设计及流程图 根据以上理论推导，本文主要使用m a t l a b t l 6 】来编制程序，m a t l a b 是用来解决 工程与科学中的的实际问题的工程软件。 1 无网格l r p i m 的主要步骤如下： ( 1 ) 地下水流问题域及其边界上布置，1 个节点，并对其编号，同时获取它们的坐标； ( 2 ) 设定每个节点的支持域； ( 3 ) 采用二重循环组装系统矩阵： 外层循环针对地下水流中的所有节点，在该层首先构造局部积分域； 内层循环针对该积分域的所有高斯积分点，计算高斯积分点的径向基点插值形函数 及导数，并形成积分点子矩阵： 还是 的子 1 6 辽宁师范大学硕士学位论文 图2 1 程序流程图 地下水模拟中的无网格局部径向基点插值法 2 2 二维不稳定流的无网格l r p i m 考虑如下二维不稳定流问题 v ( t v h ) + 占= 瓦o h ( x ，y ) q ，f ( o ，t 】 h i r i = 万 锄i 石| r = g h ( x ，y ，0 ) = h o ( x ，y ) 舯v ( t v h ) ：i ( 9 ( t - - 凳- ) + i 0 【z - i o h ) ，= o o x 铆 贮水率 否孺蓊 为初始水头。 ( 2 2 2 ) 对于每一场节点i ( i = l ，2 ，n ) 采用局部加权残量法，将得到该节点的局部弱式方程， 具体如下形式 l 彬( v ( r v h ) + 占一瓦o h ) d q = o ( 2 2 3 ) 式中的哌为权函数或检验函数。 对式( 2 2 3 ) 利用格林公式进行分部积分得如下局部弱式 f ，( 丁差谚，+ 丁等谚，。) d i 一。( r - 锄石1 痧t , x + 丁筹彬，一矿s + 哦警) d q = 。( 2 2 4 ) 局部积分域q 。的边界0 由三部分组成，即l = l ，u l 。u l ，其中：0 为积分域 的内部边界；，为与积分域相交的自然边界；。为与积分域相交的本质边界。 f 。( 丁警谚，k + r 等彬勺) d r = f “( r 差彬他+ 丁嚣谚) 订+ f 弘仃差嘭吃+ z 嚣矿勺) d r + f 胛仃o 蚴绌。n ，+ r 考谚b ) d r ( 2 2 5 ) 对水头函数采用r p i m 形函数近似表示为 h “- - z 谚吩( f ) ( 2 2 6 ) 式中办为r p i m 形函数。 将( 2 2 6 ) 式代a ( 2 2 5 ) 式可得到常微分方程 辽宁师范大学硕士学位论文 喜鬼c f ，lc z 善哦，k + r 爹彬b ，d r 一喜忽。，l c r 差哦一+ r 等彬，d q 。2 2 7 ， 一窆t = 1l 嘞脚掣= 一羔i = il 谚谢q 一 喜曩龟，仃差彬? 。+ z 爹彬侈，d r 一喜c f ，l 仃差谚x + 丁爹彬d q 但2 8 ， 一窆。( 谚谚) d q 半= 一窆l ，哌耐q 。 喜噍7 q - 1 。，c 2 差，：，+ 2 1 爹彬咒，d r 一喜“1 。，lc 丁差谚，+ r 等哦，y ，d q 但2 9 ， 一窆i = 1l ( 彬谚) d q 掣= 一主i = 1l 谚耐q 。 1 9 地下水模拟中的无网格局部径向基点插值法 3 采用无网格l r p im 的地下水问题 3 1 地下水中的二维稳定流问题 考虑如下的二维混合边界稳定井流计算问题，研究区为矩形，范围为 2 0 0 0 mx1 0 0 0 m ，研究区左右边界为一类边界，且水头值为0 m ，上下边界为隔水边界， 以( 0 , 0 ) 为坐标原点，渗透系数为t = 1 0 0 m 2 d 。在坐标为( 1 0 0 0 ，5 0 0 ) 的节点上有一半径 p = 0 9 m 的抽水井( 见图3 1 ) ，流量为q = 2 0 0 0 m 2 d ，井的边界设为。则该问题可 表示为 i 0 【= fr i o h ) + i 0 【= ，r i o h ) ：0 ( 3 1 ) h ( x ，j ，) f 。o = 0 ，五( z ，y ) l ，2 o = 0 ( 3 2 ) 乱- 0 ，乱咖= o t o h l o n = q 孑 ( 3 3 ) 钿iy ；o细k 咖 d 、7 一一般节点o 井点节点 图3 1 井点的v o r o n o i 多边形 在该算例中有井点，所以我们应对井点做特殊的处理，这里我们采用控制面积法。 在井点处取v ( z v h ) = q t 6 ( x - 孑k ，y - y k ) ，使用无网格l r p i m 时，应该使每个井点都成 为场节点，我们用凳来近似等于吼万。一夏，y 一元) ，其中吼是各个井点的流量，( 夏，或) 是井点坐标，a 。是每个井点的控制面积，即根据井点附近的节点所构成的v o r o n o i 图形 b 8 所对应的多边形( 见图3 1 ) 的面积。要求此面积就需要进行网格剖分。 辽宁师范大学硕士学位论文 图3 2 节点布置 在计算区域上按图3 2 布置节点，采用的权函数为三次样条函数 i 0 5 0 5 1 1 采用的径向基函数为 三 r f ( x ) = ( ，；2 + c 2 ) 2 二维稳定井流计算问题的解析解见图3 3 ，计算结果见图3 4 ，从绝对误差( 见图3 5 ) 可以看出利用了网格剖分的井点附近误差有点大，而其它点处的模拟结果精度较高。 0 0 图3 3 解析解 2 l r