（环境工程专业论文）复杂边界河道流场及污染物输移数值模拟研究.pdf

摘要 环境污染与环境保护是当今世界所面临的严重问题之一。由于工业的发展，大量有害于 人类和其它生物的工业废物( 废水、废气、废渣) 和生活污水、农业废水排入河流、湖泊、 水库与海湾，使天然水体受到严重的污染。为了保护环境，防止污染的危害，必须了解这些 污染物质的扩散输移规律。就是要探求当污染物投放于水体后，由于扩散、输移所造成的污 染物浓度随空间和时间的变化规律。本文研究了一种正交曲线网格生成技术导出了正交曲 线坐标系下水流及污染物浓度输移的控制方程。建立复杂边界河道污染物浓度输移的数学模 型，劳对美国c o l o r a d o 洲f a l lr i v e r 流场及汉江仙桃河段污染情况进行了数值模拟，计算结 果与实铡资料比较表明，本文提出的方法效果较好。 本文的主要研究成果如下： 1 ) 针对t h o m p s o n 的曲线网格生成方法中调节因子存在的问题，提出了一种新的调节因 子的表达式，并给出了正交曲线网格生成实例。实例检验表明，该调节因子能够对复杂边界 的单连通域或多连通域的水域生成理想的正交曲线网格，内部网格分布能够适应物理量场的 变化情形。 2 ) 利用坐标变换将笛卡尔坐标下的二维水流控制方程转换到正交曲线坐标系f ，并在正 交曲线坐标f 采用a d i 法离散、求解水流的控制方程，充分利用了贴体网格能较好拟合不 规则河岸线的特点。该网格生成技术能够有效克服河道的复杂边界给差分法带来的困难。 3 ) 采用a d i 法求解止交曲线坐标系下水流及污染物浓度输移的耦合方程时计及了水 流扩散项的影响。并用本文建立的数学模型对美国c o l o r a d o 洲f a l l r i v e r 流场及汉江仙桃河 段污染情况进行了数值模拟，计算结果与实测资料比较证明了该模型的可靠性。 总之，本文提出的模型具有形式简单、编程方便、计算收敛速度快等优点，且能基 本满足工程要求，对水利，环境工程领域中的工程实践有一定的指导作用，是一种较好的：i 。 程模型。 关键词 数值模拟，调节因子，正交曲线网格，污染物浓度输移，a d i 法，复杂边界河道 s t u d yo nn u m e i u c a ls 江u l a t i o n0 f2 dr i v e rf l o w a n ow a t e rq u a r l i t yw i t hc o m p l i c t e db o u n d a r y a b s t r a c t e n v i r o m e n t a lp o l l u t i o na n dp r o t e c t i o nj so n eo ft h eb i gp r o b l e m si nt h ew o r l d w i t ht h e d e v e l o p m e n to fi n d u s t r y , al a r g ea m o u n to fi n d u s t r yw a s t ea n dl i f ep o l l u t e dw a t e rh a r m f u l lt o p e o p l ea n do t h e rl i v i n gt h i n g si sd i s c h a r g e di n t or i v e r , l a k e r e s e r v o i ra n do c e 如a n dm a k et h e n a t u r a iw a t e rm u c hm o r ep o l l u t e d t op r o t e c te n v i r o m e n t , p e o p l em u s tl e a n la n ds t u d yt h e t r a n s p o r t a t i o no ft h ep o l l u t e dm a t t e ri nw a t e r , t oo b t a i nt h ed i s t r i b u t i o no ft h ep o l l u t e dm a t t e r c o n c e n t r a t i o n h e r e am a t h e m a t i c a lm o d e lf o rn u m e r i c a ls i m u l a t i o no f t w o d i m e n s i o n a if l o wa n d p o l l u t e dm a t t e rt r a n s p o r t a t i o nw i t hc o m p l e xb o u n d a r yu n d e ro r t h o g o n a lc u r v i l i n e a rc o o r d i n a t e s y s t e m ac o m p a r i s o nb e t w e e nc o m p u t e dr e s u l t sa n dm e a s u r e dd a t af o rp r a c t i c a le n g i n e e r i n g s h o w sas a t i s f a c t o r ya g r e e m e n t n er e s e a r c hd e t a i l so f t h i sp a p e ra r e f o i l o w s ： 1 s o m ep r o b l e m sc o n c e r n i n gt h ea d j u s t i n gf a c t o r si nt h ec o m m o n l y - u s e dt h o m p s o nm e t h o d i nt h eg e n e r a t i o no fc u r v i l i n e a rg a da r ep o i n t e do u t an e we x p r e s s i o no fa d j u s t i n gf a c t o r si s p r o p o s e d , a n ds o m ec u r v i l i n e a rg r i dg e n e r a t i o ne x a m p l e ss h o wt h a tb yu s i n gt h ep r o p o s e d a d j u s t i n gf a c t o r si d e a lg r i d sc a nb eg e n e r a t e df o rs i m p l e - - c o n n e c t e dr e g i o n sa n dc o n n e c t e d - r e g i o n s w i t hc o m p l e xb o u n d a r y t h ec u r v i l i n e a rg a di so r t h o g u n a la n dt h ei n t e d o r 西d sc a nb ea d a p t i v et o t h ev a i l a t i o no f p h y s i c a lp a r a m e t e r so f f l u i d 2 t h eg o v e r n i n ge q u a t i o n sa r et r a n s f o r m e df r o mp h y s i c a lp t a n e ( 】( yc o o r d i n a t es y s t e m ) t o c o m p u t a t i o n a lp l a n e ( e ，qo r t h o g o n a lc u r v i l i n e a rc o o r d i n a t es y s t e m ) 。a n dt h eg o v e r n i n g e q u a t i o n su n d e ro n h o g o n a lc u r v i l i n e a rc o o r d i n a t e s a r ed e e r e t i z e da n ds o l v e d t h ec o m p u t e d e x a m p l e ss h o wt h a to r t h o g o n a lc u r v i l i n e a rg a dg e n e r a t i o np r o p o s e dh e r ec a nb eu s e dt oo v e r c o m e t h ed i f f i c u l t i e sr e s u l t i n gf r o mn u m e r i c a ls i m u l a t i o no f 2 df l o ww i t hc o m p l e xb o u n d a r y 3 a d im e t h o di su s e dt os o l v et h et r a n s f o r m e de q u a t i o n su n d e ro r t h o g o n a lc u r v i l i n e a r c o o r d i n a t e s ，a n dt h ep o l l u t e dm a t t e rd i f f u s i o nt e r mi sc o n s i d e r e d a n dt h em a t h e m a t i c a lm o d e li s u s e dt on u m e r i c a ls i m u l a t i o no f a2 dp o l l u t i o np r o b l e ma n dt w o2 df l o w s t h er e s u l t ss h o wt h a t t h em o d e li sae f f e c t i v e ，c o i t e c t ，g o o dm o d e l k e y w o r d s ： n u m e r i c a ls i m u l a t i o n ，a d j u s t i n gf a c t o r s ，o r t h o g o n a lc u r v i l i n e a rg a d ，a d im e t h o d ， w a t e rq u a l i t y , c o m p l e xb o u n d a r y 第一章引言 第一章引言 1 1研究意义 环境污染与环境保护是当今世界所面临的严重问题之一。所谓环境，一般是指大气环境、 水环境和生态环境，而这三者之间又保持着密切联系。由于工业的发展，人类居住的相对集中， 大量有害于人类和其它生物的工业废物( 废水、废气、废渣) 和生活污水、农业废水排入河流、 湖泊、水库与海湾，使天然水体受到严重的污染。污染物质通过在水体中溶解扩散以及随水流 的输移而扩大影响范围。为了保护环境，防止污染的危害，必须了解这些污染物质的扩散输移 规律。即：要了解当污染物投放于水体后，由于扩散、输移所造成的污染物浓度随空间和时间 的变化规律。过去一个时期，人们在研究水资源的利用和开发时，多着眼于对水的数量做出预 测和规划。可是在当今的历史条件下，如果仍然仅仅注意对水量做出预报和规划，那就显然不 可能达到预期的目的。因为要把水作为资源，不仅有数量的要求，同时也有质量的要求。该项 研究正是为准确预报复杂边界河道水流的水质污染而提出的。 本文首先研究了一种正交曲线网格生成技术，能够有效克服河道的复杂边界给差分法带来 的困难；再导出了正交曲线坐标系下水流及污染物浓度输移的控制方程，采用a d i 法离散、求 解水流的控制方程。建立了复杂边界河道污染物浓度输移的数学模型。目的是能够准确预报复 杂边界河道水流的水质污染程度。即：预测污染物浓度随空间和时间的变化规律。因此，该项 研究成果将会为水质评价、水质预报、水质规划与管理以及水资源保护等提供决策性的依据， 对水环境的防治工程有一定的指导意义。 1 2河道流场与浓度场数值模拟的研究方法及现状 对于河道流场与浓度场的模拟计算，河道边界的复杂性及河道水流的非恒定性，给模拟河道 水流、水质的变化过程及规律带来困难。对于浅水型( 水面宽度远大于水深) 河道，由于流场是 独立于浓度场的( 除非污染物浓度太高以至于影响水的密度和粘滞系数) ，这样就可以先进行流 场的计算，然后以此为基础，计算浓度场。受计算机目前的发展水平和在很多情况下缺乏实测验 证资料的限制，现阶段对于河道流场与浓度场直接进行三维的数值计算在理论上虽然已解决，但 所需计算量和存贮空间较大，很难在工程实际中推广。目前人们处理天然河道污染物沿水体垂向 混合比较均匀的问题时，大都采用垂向平均的二维方程来描述水体运动和污染物的浓度分布。数 值计算方法主要是有限差分法( f d m ) ，包括特征线法、显式差分法、交替隐式法( a d i ) 、破 开算子法等：有限单元法( f e m ) ，已有多种单元形式及单元插值方法；有限分析法( f a m ) 、 有限体法( c v m ) 和边界单元法( b e m ) 。 有限差分法是一种传统的数值离散方法。其基本思想是”：在矩形网格上采用有限差分近似 代替微分方程中的各阶微分项进行数值离散要求所得的代数方程组在网格节点上得到满足。它 适用丁各种类型的微分方程，数学概念清晰、简单，便于编制程序，计算精度随差分格式的不同 而不同，误差估计、收敛性和稳定性理论趋于成熟和完善，是应用晟多和最成功的一种方法。但 是，在股的矩形网格差分计算中，对计算区域复杂的几何边界通常采用阶梯化近似。这种处理 方法使模拟边界难以精确地拟台实际边界，从而给计算精度带来很大影响。为了克服有限差分法 的局限性。许多学者致力于不规则边界处理问题的研究，诸如坐标变换法”“。9 0 - 9 “、l 谢利 曼所提出的任意网格有限差分法”“1 等文献( 1 8 1 9 】把不规则网格有限差分计算引入河道 流场及浓度场计算中，不规则网格差分法r 如三角形、四边形网格差分1 可以方便、灵活地处理复 西安理工大学硕士学位论文 杂计算区域的边界但是存在着差分格式比较繁琐，节点编号复杂，程序不易通用化等缺点。因 此，不规则网格差分法在河道流场及浓度场的数值模拟中有待于进一步发展和成熟。坐标变换法 ”1 是近年来常用的不规则边界的处理方法，其基本思想是：把物理平面上的任意求解区域， 通过坐标变换转化成计算平面上的矩形区域，然后在变换平面的矩形区域上，离散求解变换后的 水流控制方程。目前生成坐标变换的方法大致可以分为三类1 9 4 复变函数法；代数变换法： 解微分方程法。美国学者thomps o1 1 等人所提出的边界拟合坐标法属解微分方程法。第 一、二类变换法适用于比较简单的情况：而第三类方法适用于较复杂的情况，例如处理复杂的计 算域内含有物体、岛屿的情况。解微分方程法是处理复杂计算区域的有力工具，但其代价是求解 附加的椭圆方程来实现坐标变换，而且在计算平面( 变换平面) 上，待解的微分方程变得相当复杂。 但其最大的优势是获得了高精度的数值解。在thomps on 边界拟台坐标法中，p o s s i o n 方 程型坐标变换是一种较为常用的方法。该坐标变换理论与方法较为完善，其主要优点是数值计算 非常简单可靠，变换网格光滑趋势较好。通过源汇项的适当选取，可以方便地调节和选取网格的 集中和稀疏。通常一般的差分离散格式与数值迭代方法，及必要的边界条件处理便构成了变换的全 部计算内容。在边界拟合坐标变换中正交型坐标变换在河道流场及浓度场计算中具有特别的意 义，也正越来越受到重视，虽主要的原因在于变换网格的正交性使流体运动的变换控制方程简化 以及数值计算误差减小。尽管对于某些复杂的物理域而言，完全正交的坐标变换网格难以获得， 而正交度误差控制在一定范围内的正交型变换网格可以实现，并常能为流体数值计算所接受。 d e p o t t e r o 采用求解一定的初值问题的方法利用插值处理初始非正交网格，获得接近正交的 变换网格。h j h a u s s l i n g 4 6 1 以迭代数值求解方法类似于p o t t e r 的处理方式形成正交变换网格。 h j h a u s s l i n g 4 6 直接以变换网格的正交度构造数值形式的微分方程定解问题，将变换网格的正 交性实现于数值收敛解中，构成近似正交坐标变换。a a r o n g w a l l a 4 7 1 以保角变换的 c a u e h y - r i e m a n 方程出发，考虑变换坐标的变形因素构成椭圆型方程的坐标变换，其原理实际上 为从数值上拟台保角变换方法，g r y s k i n 4 8 在正交曲线坐标系统下以l a p l a c e 算子并同样考虑变 换坐标的变形因素构造坐标变换控制方程，显示出较强的实用性 3 1 1 。由此可见，thom ps on 边界拟合坐标法特别是正交曲线坐标变换法是改善有限差分法性能的有效方法。这些专门 技术与有限差分法配套使用，差分法就如虎添翼了。因此，有限差分法仍是流体数值计算中的一 种主要方法。 因有限差分法的具体实施过程不同可分为特征线法、显式差分法、交替隐式法( a d i ) 、 破开算子法等。特征线法的物理意义比较清楚，使用方便，计算精度高，而且不受稳定性条件的 限制，它的缺点是计算的网格不规则，要求的存贮量大，计算较繁琐，由于使用交叉网格点，也 可能有节点扰动问题发生。显式差分法有两种一为直接的显式差分，另一为交替显式法，简称 为a d e 法。虽然显式差分格式运算简便，但其稳定性较差，柯朗数要求严格，因此。该法的应 用也遇到一定的障碍。交替隐式法( a d i ) 是一个比较成熟且应用广泛的方法 2 0 】。稳定性能好。 当柯朗数较大时，会出现不应有的数值衰减和发散，文献【”1 也论证了当柯朗数等于9 5 时，流场 已高度失真。在实际计算中，一般都取柯朗数小于1 0 t “1 。算子分裂法将复杂数学物理问题分裂 简化为一系列较简单的易于求解的子问韪，其分裂方式可视所研究问题灵活选择，其出发点是要 克服a d i 法柯朗数的限制，并有可能对各子问题选用最适当的计算方法求解：但其缺陷是一维 化现象可能造成流场失真。 另一个比较常用的计算方法是有限单元法。由于一般工程中要解决的河道流场问题都具有较 复杂的边界和地形，而笛氏坐标差分法存在着网格不能弯曲，适应性差和边界条件难以处理等缺 点故人们便开始将有限元法应用到河道计算中。有限单元法的概念 8 5 - 8 6 是采用局部近拟的低 阶多项式作为试函数，构成包含因变量结点值的代数方程。有限元法能采用不同形式的不均匀网 格对于不规则边界和地形变化复杂的河道有较强的适应性。特别是对于半淹没的浅滩、河道的 分汉及沙岛所形成的弯曲岸线运用该方法都可以方便处理。但是应用有限元法求解浅水方程时， 2 第一章引言 它在空间上用有限元法，而在时间上是用有限差分法来获得的，因而需求解大型矩阵方程，它对 计算机的存贮、计算时间和费用的要求均要比有限差分法高这就限制了使用中小型计算机计算 节点较多的大范围的实际问题。另外，它在误差估计、收敛性和稳定性等方面的理论研究远不及 有限差分法严密和完善。因此，有限单元法目前仍不能满意地解决一些复杂的水流问题，特别是 用于紊流的数值计算中尚不成熟，处于探索阶段。 有限体积法是7 0 年代由s p a l d i n g 和p a t a n k e r 等人。“所提出和发展起来的一种离散方法。其 基本思想是：将计算区域划分成一系列连续但不重迭的控制体积，并使每个控制体积包围一个网 格点，将待解的微分方程对一个控制体积积分，得出一组离散方程其中的未知数是网格点上的 因变量m 的数值。为了求出控制体积的积分，必须假定中值在网格点之间的变化规律。即设定巾 值分段的分布剖面。有限体积法的基本思想易于理解，并能得出直接的物理解释：离散方程的物 理意义，就是因变量由在有限大小的控制体积中的守恒原理，如同微分方程表示因变量在无限小 的控制体积中的守恒原理一样。有限体积法得出的离散方程要求因变量的积分守恒对任意一组控 制体积都得到满足，对整个计算区域自然也得到满足。这是有限体积法的显著优点。正因如此， 有限体积法在的数值计算中也得到了广泛的应用。其缺点是在处理复杂的几何边界时，遇到的困 难与有限差分法一样。 有限分析法是美籍华人陈景仁于1 9 8 0 年提出的”，其基本思想是将古典解析法纳入偏微分 方程的数值解中。首先将待解问题的总体区域划分成许多小的子区域，在这些子区域上求解析 解，然后从局部解析解导出一个代数方程，把子区域上的内结点值与相邻的结点值联系起来汇 集成一组代数方程，再加上边界条件可解出区域内各点因变量。有限分析法具有明显的自动迎风 性质，克服了在高rey1 3 _ olds 数下有限差分数值解容易振荡或发散的缺点，计算稳定性 好，收敛速度快。但是对于双曲方程，由于其规整域上的解的表达式不易通过边界点的值表达域 内的值，所以存在一定的误差【2 5 3 2 6 】，而且有限分析法适用非规则域的性能较差，目前的处理方 法是采用贴体坐标变换，边界值外推内插等。尽管如此，有限分析法受到国内外学者的高度重视， 到1 9 8 5 年，李炜和吴江航对有限分析法的收敛性和稳定性进行了分析和证明，使其理论日益完善。 边界单元法”是与有限单元法几乎同时发展起来的一类离散方法。早先称为边界积分方程 法。在流体力学中称为有限奇点法或有限基本解法。其实质是解数学物理方程的格林函数法，据 此给出解的积分表达式，再利用定解条件建立边界积分方程。所得的边界积分方程通常不可能解 析求解，因此利用有限元离散，将其转化为由边界结点上未知量表示的代数方程组，解该代数方 程组，便求得边界结点所有的未知值。对于利用直接边界单元法并且仅对边界物理量变化的问题， 求解到此为止；对于利用间接边界单元法或对计算域内物理量变化的问题，则应将已求得的边 界结点上的未知值同定解条件一起代入解的积分表达式。计算相应的积分，得出问题的数值解。 由此可见，边界单元法可以看作格林函数法( 或边界积分方程法) 与有限元法相结合的产物。边界 单元法可以使求解问题的空间维数降低一维，从而使计算工作量与所需内存明显减少，这是边界 单元法的晟大特点。边界单元法是计算椭圆性问题的有效方法，但由于需要控制方程的基本解， 所以对于复杂的问题，如解完整的n s 方程尚未得到广泛的应用。 综上所述，采用任何一种数值方法模拟复杂边界河道流场及浓度场时。都必须解决好如下三 个问题：( 1 ) 如何处理复杂的河道几何形状。使模拟边界精确拟合实际边界? ( 2 ) 如何方便、 准确地处理边界条件? ( 3 ) 如何离散、求解水流的控制方程? 本文正是为解决好这三个问题而 提出的。近年来数值网格生成技术得到了迅速发展特别是j e t h o m p s o n 等人发展了边界拟合 坐标系的网格生成方法以来i ”】i 卅，这个方法更是被广泛地应用于实际工程问题的计算中。该方 法主要是利用数值计算方法在计算域内自动生成与边界形状一致的曲线型网格( 物理坐标系) ， 该网格与个矩形的规则网格( 变换坐标系) 互为映射，数值计算在变换坐标系统中进行。本文 基于t h o m p s o n 的边界拟合坐标变换法i 蚓，研究了一种正交曲线网格生成技术与传统的有限 3 蕊安理工夫学硕士学垃旋文 差分法结合，离散、求解水流的控制方程；与非j e 受网格技术相比较，不仅能简化控制方程，而 且可以使边界条件处理方便，还可以减少误差来源。采用a d i 法离散、求解正燮曲线坐标系下的 水流控制方程，稳定性能好，计算精度高；继承了传统有限差分法的优点。克撇了其缺骼，悬模 接复杂迷雾i 霉遂漉场教污染耪浓度塌戆一静较理慧瓣数学模型； 1 0 本文的数学模型及研究方法 零文对于天然辩道漉场的数簸模拟，采用垂觏警鹄的二维浅农方程。对于浓度糖的数值摸拨， 采焉警嚣二维运区搂麓。 为了适应天然河道地形及水力祭件的复杂性的特点，克服河道边界复杂所g i 起的数值计算的 困难- 基于边界拟坐标变换技术，研究一种正交曲线网格生成技术，将天然河道复杂、不规则 的区域变为规则区域使流场的数缎计算在坐标变换系统中进行。此外，由于撼述河道非恒定流 起蒸零方程是一疆糍线性双蠡方程筑，受格式豹稳定毪窝赣凄麓彩畹较大，零文采臻a d i 技术 解包禽扩散项影响静控耙水流的浅水方程。实际箨例验证了本模式的可靠性。 基于计算流场。求解污染的输移方程，可以得到污染物随空间的分布和随时间的演化。求 解方粳时- 将污染物输移方程转化成正交曲线坐标系下的形式，对流项用逆风格式离散，扩散项 用中心差分撂式离激，最嚣整成三对角薄方程，袋髑t d m a 技术求释。本文藏耀上逮方法，对 美国c o l o r e d o 潮f a l l r i v e r 流场及汉江伯挑弼段浮染情况进行了数值模攒。数德诗算的结果与嶷 测赉料的比较表明本文的方法可靠性好、精度高，实用意义较大。对水环境的防治工程有一寇的 指导意义。 4 第二章正交曲线哥悟生成技术研究 第二章正交曲线网格生成技术研究 实际工程问题往往具有复杂的几何边界要求计算工作者解决复杂边界条件下的数值模 拟和预测。对于非规则边界进行处理主要有如下的一些方法：近点取值法、边界插值法、邻 点加权法、假定物理参数法，边界拟合坐标变换法等。其中边界拟合坐标变换法在处理边界 时对编制统一形式的计算机软件具有较大的优越性。在应用坐标变换法求解具体问题时， 可选定某一变换形式，使物理区域边界线与坐标线重合，这样无需对边界进行特殊处理，可 以使计算在规则的变换平面内进行。但需指出，曲线网格生成的质量将直接影响问题求解的 精度，甚至影响数值计算的收敛性。理想的曲线网格应具有下列四种特性：网格线与边界的 正交性、光滑性( 连续性) 、网格疏密程度的可控性及广泛的通用性。笔者在应用t h o m p s o n 边界拟合坐标变换法进行曲线网格生成时，发现该方法中的调节因子p ，q 函数存在一些缺 陷”。，比如：可调整的参数太多、权映射现象等；不能保证生成网格的质量。鉴于此，本文 导出了一种新的调节因子的表达式，克服了上述缺陷。曲线网格生成实例表明，该调节因子 能够对复杂边界的单连通域或多连通域的水域生成理想的正交曲线网格，内部网格分布能够 适应物理量场的变化情形。该方法是种成功的方法。 2 1 边界拟合坐标变换基本理论 2 , 1 1 边界拟合坐标的基本概念 据前所述，在作物理问题的数值计算时，晟理想的坐标系是各坐标轴与所计算物体的边 界相符合的坐标系称为边界拟合坐标。当没有现成的坐标系可以利用时就希望通过计算的方 法构造这样的坐标系。试看图2 - 1 a ，设在物理平面的x v 坐标系中有一不规则区域，为了构造 个与该区域相适应的坐标系，把该区域相交的两个边界作为曲线坐标系的两个轴，记为e 及n 。在该区域的四个边上，规定不同地点的e 、n 值i 在作这种规定时要注意：在一 条边上只能一个坐标单值地发生变化，而另一个坐标则保持为常数；在两条对应边上 同一曲线坐标的最大值与最小值应当对应相等，以便在计算平面上能得出矩形区域f 图 2 1 b 、。 ( a ) 物理域 图2 - 1 坐标变换 ( b ) 计算域 西安理工大学硕士学位论文 其余的问题是：在物理平面的求解区域1 8 与任一点( x 、y ) 的相应( e 、n ) 值是多少? 一旦建 立了计算平面与物理平面求解区域内部点与点之间的对应关系，则在计算平面上所得的求解 结果就很容易转换到物理平面上去。如果e 、n 看成是物理平面上的两个未知函数，则上述 确定e 、r l 的问题就是物理平面上的一个边值问题。因此，从物理平面上来说，所谓要生成 一个边界拟台坐标，实际相当于要求解物理平面上的一个边值问题。也可以反过来叙述这一 问题：首先把物理平面上的a b c d 区域按己规定的边界上的l 、n 值，画成为计算平面的 e 、n 直角坐标系中对应的矩形f 图2 - 1 b ) ，然后以均匀的网格把矩形离散化。于是问题就变 为：已知在计算区域边界上各节点( e 、q ) 相应的( x 、y ) ，问在计算区域内部与任一点( e 、 q ) 相应的( x 、y ) 值是多少? 这样如把( x 、v ) 看成是计算平面上的未知函数，则生成边界拟 台坐标的问题也就是计算平面中的一个边值问题。从数值计算的观点，对生成的边界拟合坐 标有以f 几个要求：物理平面上的节点应与计算平面上的节点一一对应，同一簇中的曲 线不能相交，不同簇中的两曲线仅能相交一次。；在边界拟合坐标系中的每个节点麻当 是一系列曲线坐标轴的交点而不是一群三角形元素的顶点或个无序的点群，以便设计有 效、经济的算法及程序”“，要做到这一点，只要在计算平面中采f ；矩形阿格即可；物理 平面求解区域内部的网格疏密程度要易于控制：在边界拟合坐标的边界上罔格线最好与 边界线正交或接近于正交，以便于边界条件的离散化 0 0 1o 实现、两点要求的方法将在下 章中讨论。 目前生成边界拟合坐标的方法大致可以分为三类：复变函数法、代数变换法和解微分方 程法。分述如下： 1 1 复变函数法。利用复变函数的理论可以把二维不规则区域变换成矩形区域，而且可以 得出解析的或部分解析的变换关系式。例如要把物理平面上的环形区域变成计算平面上的矩 形区域，可将物理平面上的极坐标，一口f 图2 - 2 a ) 作为计算平面中的直角坐标f 图2 - 2 b ) ，于是环 形域即转化成为计算平面上的矩形域。r 一目与x y 间的关系为 x = r c o s o y = r s i n 臼 或 r = ；7口= t a n 一1 ( y x ) 实际上，这就是复变函数中的变换 w = p 2 = p 。+ 9 这里肜及z 各为物理平面及计算平面上的复数。 ( a ) r r 2 r l 图2 - 2变换举例 6 ( b ) 第二章正交曲线同格生成技术研究 用复变函数理论来构造适体坐标的方法仅限于二维问题。 2 ) 代数变换法。这是利用一些代数关系式来进行区域变换的方法。 3 、解微分方程的方法。通过求解边值问题的微分方程来建立物理平面与计算平面上各点 间的对应关系。至于这边值问题控制方程的类型物理问题本身并无任何限定。这就给我 们一定的自由度，以便按照对所生成网格的要求来选择控制方程。这类方法将是我们讨论的 重点。 2 1 2 生成网格的微分方程法 前己指出，边界拟合坐标的网格生成问题，实际上是一个边值问题。边值问题的求解是 偏微分方程领域中的一个经典课题。应用这种方法时，我们可以利用微分方程的一些性质使 所生成的网格更完善、合理。这个方法最早是由w i n s l o w 在1 9 6 7 年提出的”j ，以后不少研 究者都对此法的发展做出过贡献。但比鞍全面而系统地研究这一方法的当推t h o m p s o n ( 汤 姆逊1 、t h a m e s ( 泰姆斯) 及m a r t i n ( 马丁) 1 9 7 4 年的论文( 9 3 3o 此后，在流体力学与传热学的 数值计算研究中就逐渐形成了一个分支领域网格生成技术。文献中所谓的t r m 方法就 是指通过求解微分方程生成网格的方法( t f m 系t h o m p s o n 、t h a m e s 及m a r t i n 三人姓的第 一字母1 ，在文献 9 1 e 9 4 3 中有详细的综述，有关的数学理论可见文献 9 5 。 二阶的椭圆型方程是描写边值问题的最简单的椭圆型方程就是l a p l a c e 方程。根据 l a p l a c e 方程解的唯一性及极值原理，可以把孝，矿作为物理平面上l a p l a c e 方程的解，即 v | = h + 。= 0 ( 2 1 a ) v 2 r = + = 0 ( 2 1 b ) 同时在物理平面的求解区域边界上规定e ( x ，y ) 、n ( x ，y ) 的取值方法，于是就形成了物 理平面上的第一类边界条件的l a p l a c e 问题。f 也可以在部分边界上规定、n 对x 、y 的导 数) 。在具体实施时，取定物理边界上有限个节点上的fe 、n ) 值，然后用数值方法确定内 部节点上的值。关于l a p l a c e 方程的数值计算问题已研究得很成熟，但由于物理平面上是个 不规则区域，于是在物理平面上解这一问题又碰到了不规则边界的困难。如果从计算平面上 的边值问题出发来考虑，则情况就大为改观，因为在计算平面上可以永远取成一个规则区域。 所谓计算平面上的边值问题就是指在计算平面的矩形边界上规定x ( ，n ) 、y ( ，n ) 的 取值方法，然后通过求解微分方程来确定计算区域内部各点的( x 、y ) 值，即找出与计算平面 内各点相应的物理平面上的坐标。另外，在计算平面上来求解时如果在部分边界上规定了x ， y 对e 、n 的导数而不是x ，y 之值，这就可能在曲线坐标系上使网格线与这部分边界正交， 但此时这部分边界上网格线与边界的交点就由计算确定。实际上大多数生成网格的情形都是 求解计算平面上的边值问题。为此，需把物理平面上的l a p l a c e 方程转换到计算平面上以e 、 n 为自变量的方程。我们以甲2 掌= 0 为例来进行转换。 为了保证两个平面上点的一一对应关系，要求雅可比行列式 旦业o( 2 2 ) o ( x ，y ) 因为正变换雅可比行列式成立，则其逆变换存在，且其雅可比行列式 嚣安理工丈学醺士学位论文 j ：垫业啦o a ( 善，叩) 正遂交接熬耩笼嚣列式毒 _ 矧关系 塑! 垫艘：i a ( x ，) a ( 考，叩) 列霹得 丝一l o y 丝。一三鱼 院j8 研 卸 je q 同理可得到它们的二= 阶导数的表达式 ( 2 3 、 ( 2 4 ) 票三j 嚣塑o y = 妻t ：固教8 j8 e 、。 警= 昙髻，= 妄e 孚，= 亨 y ，b 时( z y 。一r 。y f ) 一x 封y ；一，孵x ；y 。+ z 时y 一，+ y 错x 。y 。 一 岁# y r t r t ( x 夕口一聋口y p 一；，：秘一y w x c y 口+ x r t y e y # + y 秘x 夕q 】 等= 号h k 瓢矿训飞即。氇。x q + x ；, t y ；x q + y 并】+ k ，( 吩y 。一x , y g ) 一芏。y 。麓。一x c x 。y ，+ x ，强x ，+ x ：j 锄j 将此两式代入v 2 善= 0 ，并经整理得 善，( j ，；x ：) 一2 y 孵( y c y ，+ x f x 。) y 。( + y ；) j x y 。k ( y ：+ x ：) 一2 x 蝣( _ ) ， y ，+ x ；x q ) + x 。( + _ y ；) 】 ( a ) 类叛建v 2 翠= 0 可纯巍 七k ( y ；+ x ；) 一2 y r ( y y 。+ 赡) + 。( x ；+ y ；) j ； 坎k ( 2 + 2 ) 一2 x 。e ( y ，妊+ ) + 南。( + 露) 】 浚慧掰式( a ) 与( b ) 的等号表边方括号内瓣项烧稽等的，等号右边方括号内鹣部分也是相 等的，把等号两侧这两个方括号内的部分分别记为a 及b ，然艏以。( 口) 一x 。( 6 ) ，得 瓠；y 一b x f y 4 = 0b ( x i x ；、= b j = 0 夜转换区中j 0 ，因而只能是b = o 于是有 x 错( x ；+ y ；) 一2 x 时( 鼍。口+ y y _ ) + x 目( x ；+ y h = 0 酆 第二耄难交韭线曩接擞或技术研究 类似地可以证明 测茹一2 皿卸+ 口口= 0 蝌饕一2 国轴+ t y m ，= 0 熟中a = x ；十y ；，p = x x 口+ y y _ ，y = x ：+ y ；，j = 工；y 口一x 目y 。 式f 2 ，囝虽然嫩2 1 】复杂，毽计算匿域却麓肇。蠲鼗建方法强诗舞乎嚣豹矩形城是枣簿式 f 2 6 ) ，裁可获褥与计算平露上箨节点f ；，q ) 褶对应斡穆理平蕊上节点缒坐标( x ，n ，腻磊可 嗣有限差分法进一步计算所需的一系列几何参数。 最后要指出，虽然计算平丽上网格是等分的，但并不意味糟物理平面上生成的网格也是 等分的。实际上由于太多数变换躲非线性特点，物理平面上憋嬲格多为菲均分的。至于如 露控镶穆蓬平嚣上嬲据薛蘸密程艘，将在下一节潜论。 2 1 - 3 网格分布的控制 在进行数值计算时，常常希望在物理平面上的阔格划分能避成区域中物理鬣场的变化情 形，帮在交往剧烈艇稠密一些，焱交纯平缓姓秘蘸一些；另静，飙逑雾条咎离散健豹角瘦， 希麓两格线与镪瑷区域的边器线藏变，班稠子较准确缝计算逑器上豹物理量。掰有这些要求 部属于网格分布的控制问题。网格分布对于获得一个较好的数德解有很大的影响，因而它是 网格生成理论中的个重要研究内容。最简单的调节物理平面上网格疏密程度的方法是在计 撵平面豹边界上恰当地设置x ，y 值。铡如，在物理平瑶上要求嘲格线比较稠密的遗区，在 谤髯乎垂熬稳寂缝焦上相罄嚣瓣接线翔x 蕞y 黛拣赫交托裁取褥枣一些。毽是逸器节点蓬懿 设置仅在靠近边界的区域内有影响，而在离开边界较远的内部则主要由方程本身的特性所决 定。由于l a p l a c e 方程具有使被求解函数均匀与光滑化的特点，在远离边界内部由l a p l a c e 方程所生成的物理平面上的网格将趋于均匀分稚。因而想用l a p l a c e 方程，通过使计算平面 迭赛土x 。y 不均麓分蠢静方法采获褥熟蜀2 3 8 豹潮捂努毒是不萄爱熬f 为蕊硬起冕，蛰理平 露上静求解区域 疑赫成了矩形，褥更能产生如豳2 3 b 所示斡傣形。由l a p l a c e 方程生成的 网格还有一个特点，即在曲面边界附近，随曲磷凹向的不同，网格线会发生渐稀戚渐密的变 化，前者发生于曲面向计算医域内凸出时( 图2 - 4 a ) ，后者发生于曲面凹向计簿区域时( 图 2 4 蜕。由l a p l a c e 方程所生成嬲格的这种分布特性，对提高数德计算结果的搬确性常是不 羁的。 ，函- - 7 函 ( s t ) 图2 - 3l a p l a c e 方程嫩成嗣格的特性 ( b ) 西安理工大学硕士学位论文 ( a ) ( b ) 图2 - 4 曲面附近阿线琉密程度的变化 我们如果把物理平面上的等r l 或等e 线看成是稳态导热问题的等温线，则等温线密集的 程度就反映了局部地区热流密度的不同。如果某一局部地区有热源存在，则该地区的等温线 必定更密集。这一联想启发我们应用带“源项”的l a p l a c e 方程即p o i s s o n 方程来生成网格。 即在物理平面上采用下列方程 孝+ f y = p ( 亭，玎) r 。+ r ) r = q ( 六口) ( 2 7 a ) ( 2 7 b ) 这里源函数( 或分布函数) p 、q 是用来调节区域内部网格分布的。同样从数值计算的角度 在计算平面上解方程比较方便。采用类似于导出式( 2 6 ) 的方法可导得在计算平面上关丁x ， y 的偏微分方程 “嚣一2 k , f 口+ 肄q 目+ d 2 ( + o x 目) = 0( 2 8 a ) c o , 鬈一2 p 妇+ 目”+ j 2 ( a + q 譬口) = 0 其中o f = x ：+ y _ 2 ，声= 工 x q + y y 口，y = x ；+ y ；，j = x y _ 一x 口y 。 ( 2 8 b ) 这样，函数p 、q 的构造就成t 控制网格分布的主要因素。 关于控制函数p 和q 的选择，t h o m p s o n 等建议3 2 1 尸( 硎= 一叩咖( 一身) e x “咱睁劬一 ? ( 2 9 a ) 如j i g n ( 一珈x “一以( f 一卯+ ( 口哪) 2 ) o ( 4 ，1 7 ) = 一 鲫( 口一r ，) e x p ( 1 卜聃1 ) 一 ? ( 2 9 b ) 6 ，嘲( 口一, 1 , ) e x p ( - - d ，( f 一2 + ( 口一m ) 2 ) j - 1 式中ai ，ci ，b i ，d i 是调整网格线的疏密程度的参数。a i ，b 7 - 为密集强度；c ，d 广 衰减因子：s i g n ( x ) - - 符号函数；n 一网格线数日；m 一网格数。所有的参数都是为调整网格 0 第二章正交曲线罔格生成技术研究 的质量。 以上所述，是二维曲线网格生成的一般表示法，其e 、n 线的方向可以任意布设，不能保 证正交。正交曲线坐标系有许多显著的优点：比如：( 1 1 可以简化计算平面上的控制方程， 还可以避免由于网格的非正交性而引起的某些离散误差；( 2 1 虽然实际流体的流动方向并不 与网格线完全一致，但两者间的偏差可望小丁二按一般方法所生成的非正交网格中的偏差，有 利于减少假扩散误差；( 3 ) 有利于数值解的收敛。( 4 ) 利于方便准确处理边界条件等。因此， 如何构造一种调节因子e o 函数，生成理想的正交曲线网格便是本文的一项重要研究内容。 关于正交曲线网格生成技术的研究将在下节介绍。 2 2 正交曲线网格生成研究 笔者在实际的网格生成中，发现t h o m p s o n 调节因子可能存在f 列几个问题：可调整 的参数过多，当网格生成质量不理想时，很难确定调整哪一个参数才能改善网格的质量； 网格生成迭代过程中不易收敛；参数选取不当易引起双映射现象；不能保证生成的网格 在边界处正交。若克服了这些缺陷，t h o m p s o n 的曲线网格坐