（环境科学专业论文）有限体积法在平面二维水流—水质数值模拟计算中的应用.pdf

中山大学硕士学位论文 a p p l i c a t i o n o ff v mi n d e p t h - a v e r a g e d 2 d f l o w - p o l l u t a n t sn u m e r i c a l m o d e l i n g m a j o r ： e n v i r o n m e n t a ls c i e n c e n a m e ：z h a n gj i a f u s u p e r v i s o r ：p r o f h u a n gp i n g a b s t r a c t w k ht h er a p i dd e v e l o p m e n to f a g r i c u l t u r ea n di n d u s t r ya n dt h er a p i di n c r e a s eo f c i v i l i z e dd e g l e e ，a g r e a td e a lo f w a s t ew a t e rf l o wi n t or i v e r sa n dl a k e s a s ar e s u k ， w a t e rb o d yi sp o l l u t e da n dt h eu s e so fw a t e rr e s o u r c e sa r ea f f e c t e d i no r d e rt o e l i m i n a t ep o l l u t i o n sa n dp r o t e c tw a t e rr e s o u r c e s ，w ea r eb o u n dt oe n c o u n t e rp r o b l e m s o f a s s e s s i n gw a t e rq u a l i t i e sq u a n t i f i c a t i o n a l l y w a t e rq u a l i t yn u m e r i c a lm o d e l sa r e w i l d l y u s e di nq u a n t i f i c a i o n a l l ya s s e s s m e n t so f w a t e r q u a l i t i e sf o rt h e r e a s o nt h a tt h e y a r es i m p l e ，1 e s st i m ea n dc h a r g en e e d e d ，h i g h l yr e p e a t a b l ea n dh a v es t r o n gf o r e c a s t a b i l i t i e s t h e r e f o r e ，t oe s t a b l i s has i m p l e ，p r a c t i c a ia n de f f e c t i v ew a t e rq u a l i t y n u m e r i c a lm o d e l ，a n ds o i v e dw i t h p r e c i s e ，s i m p l e b u tp r a c t i c a lm e t h o d ，i st h ek e yf o r w a t e r q u a l i t yf o r e c a s t i n g ，p l a n n i n g ，w a t e rp o l l u t a n t se m i s s i o ns t a n d a r de s t a b l i s h i n g a n de n v i r o n m e n t a l e n g i n e e r i n gd e s i g n i n g t h i st h e s i sn o to n l yd i s c u s s e st h ep r i m a r yn u n l e r i cm e t h o d sf o rf l o w - p o l l u t a n t s c o u p l e dm o d e l i n g b u ta l s oc o m p a r e st h e mw i t he a c ho t h e ra n de x p l a i nt h es t r o n g p o i n t a n ds h o r t c o m i n go fe a c ho n e b a s e do nt h ec o m p a r i s o n ，t h i st h e s i se m p h a s i z eo n d e d u c i n gt h e f v ma n di t sf v st h e m e ，a n dt h e nd e v e l o pad e p t h - a v e r a g e d2 d f l o w - p o l l u t a n t sc o u p l e dm o d e lt os i m u l a t ef l o wa n dc o n c e n t r a t i o nf i e l d ，t a k i n gt h e y a m i a nh a r b o rf o re x a m p l e o u rr e s e a r c hs h o w st h a tf v mw h i c hi sb a s e do nu n s t r u c t u r e d 鲋c a np e r f e c t l y a d a p t t o c o m p l e xb o u n d a r y a n d t o p o g r a p h y ，a n ds t r i c t l y m e e tt h e p h y s i c a l c o n v e r s a t i o nl a w s t h o u g hf v mi se a s yt ob ec o n s t r u c t e d ，i th a sh i g hn u m e r i c a l p r e c i s i o n t h e r e f o r e ，f v m i sah i g h q u a l i t yf l o w - p o l l u t a n t sn u m e r i c m e t h o d k e yw o r d s ：f l o w - p o l l u t a n t s n u m e r i c a l m o d e l i n g ，u n s t r u c t u r e dg r i d ，f v m ，f v s i i 第l 章引言 第1 章引言 1 1 问题的提出及研究意义 水是地球上分布最广泛的一种物质，也是一种非常重要的自然资源。随着工 农业生产的迅速发展与城市化程度的快速提高，大量的废水流入江河湖海，造成 了水体污染，影响了人们对水资源的使用。为了消除污染，保护水资源，必然要 碰到定量评价水环境质量的问题1 。比如为了改善一条有机污染严重，溶解氧 浓度低于给定标准的河流的水质，可能要建设一座污水处理厂。在规划设计这个 污水处理厂之前应该明确回答处理厂应当建在什么地方，规模多大及处理率多高 的问题；座规划中的发电厂准备把它的冷却水排入附近河流，为了评价电厂 建成后是否会对河流构成热污染，环保部门必须根据它的设计排水量和水温确定 所能造成的河水温升，以便迸一步确定这种温升是否会对水生生态系统，特别是 对鱼类产生明显影响。过量的营养物质( 如含磷、氮的物质) 流入湖泊或水库 会造成水体富营养化。为防止某一湖泊或水库发生富营养化，应该把流入的营养 物质限制在什么水平上也是实际工作中常常需要解决的问题。要解决这些问题可 以采用物理模型或数学模型进行模拟预测，物理模型由于受经费、时间和观察精 度等的限制，应用的较少；而数学模型则由于简单快捷，需要的时间和费用较少， 可重复性高且预测性能强而被广泛采用，数学模型经特殊简化可以得到解析解， 但更多依赖于数值解，数值方法是求解水环境数学模型最有效的方法t 2 1 。 因此，建立简便、实用、真实而有效的水质模型，并运用精确、守恒、简单 实用的数值解法进行求解，是进行水质预测、规划，水污染排放标准制定以及环 境工程设计的关键所在。本文重点研究有限体积法在平面二维水流一水质数值模 拟计算中的应用，为水环境数学模型的精确求解提供强有力的支持和保障，此也 为本论文的研究意义所在。 中山大学硕士学位论文 1 2 模型和方法的选择以及本论文的主要内容 从水质模型的空间维数来说，虽然所有物体在空间上都是三维的，但水质模 型可以简化为零维、一维、二维和三维的数学模型。水质模型的维数，主要取决 于所研究的范围及其水体中污染物的混和情况。如果对区域性水质进行粗略的规 划估算，零维模型就可以了。如果对一个较长的河段或河流进行水质预测或规划， 一维模型已能得出较好的结果。如果要研究河流局部河段范围内的水质情况、排 污口附近的污染物分布、功能区、混和区及污染带宽等这类问题时，就要考虑采 用二维、甚至三维的模型。通常情况下，对于宽而浅的河流或河口来说平面二维 模型就已经可以获得足够精度的流场和浓度场了，如果不切合实际地使用三维模 型，不但工作量大，具体紊动参数难以确定，而且也没有此必要n 川 5 1 0 所以 本文从实际应用出发，主要讨论平面二维的水流一水质数值模拟。 从水质模型计算网格的剖分方法来说，水质模型的计算方法可以分为基于矩 形网格的计算方法和基于拟合曲线网格的计算方法。天然水域岸边界通常都是不 规则的，以传统的矩形网格来剖分天然水域，岸边界将被概化为折线型，若网 格过粗，在近岸水流的模拟上会产生较大的误差；若采用加密网格的办法，必将 使网格数量增加，从而带来计算量的增加。拟合贴体技术的引入是解决这一问题 的较好途径之一，拟合曲线网格可以分为有结构和无结构两类，有结构网格通过 某种数学变换( 又分为代数变换和微分变换等) 将不规则平面转换成规则平面， 其虽有较好的拟合精度，但曲线网格实现起来的工作量相对较大；无结构网格则 根据实际地形任意布设，减轻了网格生成的工作量哺1 7 1 1 8 1 因此本文采用了国 内外先进的基于无结构网格的有限体积法进行二维水流一水质的数值模拟计算。 本文主要包括的内容有： ( 1 )介绍水流一水质数学模型的发展及其优缺点，简单介绍数值模拟 计算的步骤。 ( 2 )比较和分析数学模型常用的几种数值方法，并选择最合适本文所 研究问题的有效的数值离散方法。 ( 3 )根据控制方程和数值方法，选用性能优越、守恒型的有限体积法 第1 章引言 ( 4 ) 离散控制方程，并简单介绍模拟计算的技术要点。 利用选定的数学模型程序，结合实例，对所关注的问题进行数值 计算，根据计算结果分析和评价。 中山大学硕士学位论文 第2 章水流一水质数学模型及其计算方法 2 1 水流水质数学模型简介 水流一水质数学模型是将水动力问题和污染物在水体中的迁移转化问题用 数学方程进行描述，并在一定的定解条件下求解这些方程，从而模拟出某个理论 和工程实际问题的水流、水质状况。 最早的水质模型由s t r e e t e r h w 和p h e l p s e b 在1 9 2 5 年在俄亥俄 河水污染及自然扩散的研究中提出，至今已有接近8 0 年的历史。1 9 2 5 年1 9 6 5 年期间，根据s t r e e t e r h w 和p h e l p s e b 提出的理论，水环境科学工作者 开发了较为简单的b o b d o 模型，对河流、河口采用了一维计算。1 9 6 5 年1 9 7 0 年，随着电子计算机的应用及人们对生化耗氧认识的深入，除继续发展b o d d o 模型的多参数估值外，水质模型发展为6 个线性系统，计算方法由一维发展到二 维，并开始研究湖泊及海湾问题。1 9 7 0 年1 9 7 5 年，研究发展了相互作用的非 线性系统，涉及到营养物的循环过程，浮游动植物系统以及生物生长率同营养物 质、阳光、温度的关系。计算方法一般用数值解法。1 9 7 5 年以后，除了进一步 研究食物链问题外，还发展了多种相互作用系统，空间维数发展到三维。这些已 有模型能较好地反映客观实际情况；一般有机物的水质模型已基本成熟坤川们。 实际上，由于水质数学模拟的基础是水流数值模拟，因此水质模型的发展也 紧跟着水流模型的发展，在一、二、三维不稳定流数值模拟的基础上，水质模型 也从一维稳态模型逐渐发展到一、二、三维的动态模型。应用范围也由单一河道 扩展到网河、水库、湖泊和海湾川。水质数学模型是与水流数学模型紧密结合 的，如果没有水流模型计算提供的流场，水质数学模型的数值计算就无法进行 下去。水流模型和水质模型的结合通常分为耦合和非耦合两种方式，水质数学模 型和水质数学模型的数值解法基本上是一致的，本文采用水流一水质耦合的数学 模型。 4 第2 章水流一水质数学模型及其计算方法 2 1 1 数学模型的优点 水流一水质数学模型可以说是在计算机上进行水流一水质试验的技术，但 同物理模型试验相比有其独特的优势： ( 1 ) 可以完全自由地改变或控制流体以及污染物的性质： ( 2 ) 可以进行严格的一维、二维和三维流动的“试验” ( 3 ) 可以完全自由地选择流动、扩教参数，如糙率、扩散系数等 ( 4 ) 可以进行全尺度“试验”，不存在比尺效应 ( 5 ) 可以进行物理现象各种理论近似处理的有效性和敏感性试验 ( 6 ) 在计算机上模拟，无需其他仪器和设备，方便经济； 2 1 2 数学模型的不足与缺陷 利用数学模型求解控制方程组，必须借助各种计算方法离散方程和计算区 域，因而数值计算存在一下一些不足： ( 1 ) 不同的离散方法求解同一个问题，结果可能不同 计算区域的方程的离散，不仅会引起量上的误差，而且处理不当可能 会改变方程的性质； 描述无粘性流体流动的方程离散后，会引入无物理意义的数值粘性， 这在一定程度上改变了流体的性质和流场的形态； 许多研究问题目前从数学上还不能严格证明解的存在性，唯一性，以 及解的收敛性，因而其结果具有不确定性； 目前多种紊流模式被引用到实际计算，但均基于一系列假定，须引入 较多的试验常数，实验常数的引入和处理与个人的经验和技巧有关。 ) ) ) ) 他 中山大学硕士学位论文 2 1 3 数值模拟的基本步骤 数值计算的离散方法很多，但主要的计算步骤基本一致 问题的定义：即进行数学建模工作，确定控制方程和定解条件，界定 计算区域，预选参数，如糙率、粘滞系数、扩散系数等； 离散化：首先选定离散方法，将方程离散和区域几何离散( 划分网格 或单元) ，离散后的方程转化为以节点未知量表示的代数方程组； 解代数方程组：选定一种解法，编制计算程序，求出节点数值解； 结果整理：列表、绘图、显示计算结果，必要时进行第二级计算。 2 2 水流一水质数学模型计算方法 描述水流运动和物质输运的控制方程为非线性方程组，通常没有解析解，只 能借助于计算方法求得数值解。数值解的相容性、收敛性、稳定性和实际水流、 物质输运的符合程度，很大程度上取决于数学方程的离散形式，即采用的数值 解法。因而在数值计算中，计算方法选择得是否合理是数值计算成败的关键所在。 本章工作是比较现有各种主要的数学模型计算方法，并根据它们的特点提出适合 本文所研究问题的方法。 数学模型的发展历程，很大一部分取决于各种离散方法的发展过程。由于离 散基本原理和方法不同，离散方法大致分为： 有限差分法 有限单元法 有限分析法 有限体积法 边界单元法 6 ) ) ) ) n 心 h ，f、，l 型 型 域 界 区 边 厂iiij。、ll，i一， 法方散离 第z 章水流一水质数学模型及其计算方法 区域型方法是将所讨论问题的区域分成许多单元( 子区域) ，在子区域上将 控制微分方程用全部或部分满足边界条件的函数来近似，其特点是在边界上满足 边界条件，域内近似满足微分方程。以下对主要的数值计算方法做一些简单的介 绍和比较。 2 2 。1 有限差分法( f i n i t e - - d i f f e r e n c em e t h o d ，简写f d m ) f d m 是最早提出的方程离散方法。早在1 9 2 8 年c o u r a n t 等人就提出了有限差 分理论，但真正应用于流体力学的计算是在六十年代后，最初多用正规网格和松 弛解法，1 9 6 8 年引入交替方向隐式差分法，七十年代末提出上游加权有限差分 法，至今已形成很多成熟的算法格式和软件。目前的研究前沿是非线性问题的有 效算法及其理论分析。 f d m 是将微分方程中的各微分项离散成在微小网格上各邻接格点的差商形 式，得到一个以各节点上函数值为未知变量的代数方程，称差分方程。并根据 原问题的初值和边界值合理地给出离散化代数方程的初、边值条件，从而求出控 制方程的数值解。在潮汐河道数值计算中，当前最广泛的有以下几种方法：就一 维问题，有显、隐式差分，特征线法等；雨对于多维闯题，有a d i 法，破开算子 法等。 2 。2 。1 。l 直接差分法 直接差分法是将微分方程中的微商直接用差商代替，从而得到差分方程。其 优点是离散原理简单，数学推演和编写程序的工作量小。缺点是嬲格布置不灵活， 网格无法与不规则边界完全符合。 直接差分法有显式差分和隐式差分两种方法。显式差分是最最应用的一种方 法，这种方法简单，易理解，是其他各种计算方法的基础。然而，显式差分收敛 性和稳定性较差，时间步长受到限制，且往往得不到较好的结果。 中山大学硕士学位论文 隐式差分是在显式差分的基础上发展起来的，它不能直接求解，而必须利用 联立方程组才能求出函数值。隐式差分克服了显式差分稳定性差、时间步长受限 制、精度低等的缺点。现在广泛使用的三种格式有l a x - - w e n d r o f f 格式、a b o t t 格式和p r e i s s m a n n 格式。实践表明，隐式差分格式在一维水流一水质计算中是 一种行之有效的方法1 。对于一般的河流，作为初步计算，均可获得较为满意 的结果。但在多维问题中处理比较繁琐，较少使用。 2 。2 1 2 特征线法 2 0 世纪5 0 年代初，林秉南首先提出一维水流计算的特征线法。它在数学上 是通过坐标的转换寻求可以将偏微分方程转化为全微分的一组特征线，然后沿特 征线积分求解。由于特征线法的计算网格不规则，同时要存储节点在( x ，t ) 平 面上的位置坐标和节点上的水力参数，要求存储量大。所以，在实际计算中一般 均采用特征线差分法，它综合了特征线法和差分法两种方法的优点，网格的节点 位置仍按照一般差分法布置，但格点之间的关系按照特征线法原理建立。 特征线法具有计算精度高，不受稳定性限制的优点。它最大的优点是能应用 于涌潮地区，这是其他方法所不能比拟的。因为涌潮地区，潮头破裂，沿x 方向 发生了不连续现象，这时用一般差分方法不能求解，但它的特征线仍然存在，因 此可用特征线法求解2 ”。然而，和隐式差分一样，这种方法在多维问题 中显得很繁琐。 2 2 1 3a d i 法 a d i 法也称交替隐式，是由p e a c e m a n - - r a c h f o r d 和d o u g l a s 于1 9 5 5 年提出 的一种用于专解二维问题的特殊分步法。这种方法是将一个时间步长分为两个半 步，用显、隐交替的方法求解，这样做可以把原来较大的系数矩阵转化为两个三 对角系数矩阵，然后用追赶法求解。可以证明a d i 法的隐式格式具有( a t ，x 2 ) 阶精度，而且是无条件稳定的。由于在差分格式中引入中间步长，在每半个步长 中只需解一个三对角矩阵的代数方程，所以计算量比较小。但是a d i 法由于坐 第2 章水流一水质数学模型及其计算方法 标方向的限制，水域边界不得不概化成锯齿形，加上网格划分必须是等间距的， 因而在计算中，由于边界网格化所引起的失实比较严重，而且重点计算水域不能 加密网格，因此a d i 法的应用也有一定的局限性“4 “”1 。 2 2 1 4 破开算子法 破开算子法，也称分步法，是由苏联学者y a n a n k o 等人，于五十年代末到六 十年代初提出的，并在理论上逐渐加以完善。破开算子法是数值计算中的一种分 裂解法。在数学上，它将一个微分算子分裂成几个简单算子的组合，在各分步长 上，分别求出各简单算子的解，从而得出整个步长的函数值。破开算子法，不像 a d i 法，比较灵活，可以使用不同组的方程各自按其最佳方式求解，甚至可以将 欧拉空间与拉格朗同空间结合起来，以改善稳定性。而且破开算子法可以用于多 位问题的数值计算。但由于方程组的分割，也形成了一些新的矛盾，反而影响计 算精度6 “。 2 2 1 5 边界拟合法 边界拟合法是在计算区域上求解椭圆方程，得出物理区域上网格的位置，然 后在该计算区域上导出求解方程。它能够拟合不规则边界，提高了差分法的使用 优势。边界拟合法大概可以分为有结构和无结构网格两类。 有结构网格是通过某种数学变换( 一般为代数或微分变换) ，在物理计算域 内生产曲线网格，使得求解的边界与网格曲线相重合，这样再将贴体坐标转换 为直角坐标，实现将不规则平面转化成规则平面，再次基础上对基本方程进行离 散求解。自1 9 6 7 年w i n s l o w 首次提出以微分方程法建立贴体坐标以来，不少学 者对此展开深入研究，在不同方面完善了发展了曲线贴体网格技术。但曲线网格 推导、离散过程繁琐复杂，计算工作量大，所以应用不是很广泛1 1 8 3 t 1 9 1 。 近年来有学者提出在无结构网格上对基本方程进行离散。如1 9 8 3 年赵士清 提出的三角形网格下差分格式，用三角形网格将计算区域分解为许多单元，通过 9 中山大学硕士学位论文 函数在单元内线性化假设对基本方程进行离散，实现无结构网格上的边界拟合技 术2 0 “2 1 物。 总的来说，有限差分法离散原理简单，数学概念明确，编程简单，计算速度 快，因而得到广泛应用。但有的差分格式是有条件收敛的，尤其对于不规则边界 很难得到较好的处理，从而影响计算精度。 2 。2 。2 有限元法( f i n i t e - - e l e m e n tm e t h o d ，简写f e m ) f e m 是1 9 5 0 年由一些飞机结构师首先提出的一种数值方法，六十年代有限元 法开始应用于流体力学，最初是联解一些线性问题，到七十年代已推广应用于非 线性问题的研究。 f e m 是在积分最小值原理的基础上建立起来的，是将求解问题的连续区域任 意划分为适当形状( 三角形或四边形) 的许多微小单元( 称为子域) ：各小单元 分片构造插值函数，然后根据极值原理( 变分原理或加权余量法) ，将微分方程 组化为控制所以孤立单元的有限元方程。然后，将所有局部单元上的方程汇集成 总体的微分方程组或代数方程组，再换上应有的边界条件和初始条件，形成一个 完备的代数方程组，求解该方程组就得到微分方程在整个计算区域上的数值解。 根据建立有限元方程的依据不同，f e m 可以分为g a l e r k i n 法，变分法、最小二 乘法、配置法等。 f e m 的优点是采用不规则的无结构网格，易于较准确地逼近水体边界和水下 地形，后者正是决定水体在重力下如何流动的主要因素。但对于非恒定流，每一 时间步要求解一个大型的线性方程组，消耗的机时多。同时，传统的f e m 在数学 上适于求解椭圆型方程的边值问题，不适于求解以对流为主的输运问题，因此后 来又发展了逆风、最小二乘、特征耗散等改良f e m 方法，但总的来说算法过于复 杂，消耗机时也很多。国外水环境工作者在这方面的研究比较多，如美国的 r a 。w a l t e r s ，英国的z i e n k i e w i c s 等，他们在计算时计算网格动辄数十万个以 上，多用于海湾、河口等宽大水域的计算心卜枷。 l o 第2 章水流一水质数学模型及其计算方法 为了解决f e m 在计算非恒定流时的繁琐和计算时间长等问题，1 9 8 4 年吴江航 提出一种分步杂交法，主要特点是在不规则的有限元三角形网格上，对于对流算 子和扩散算子，先用破开算子法把他们分开，再分步采用特征线法和集中质量有 限元法求解，用于渤海湾水流一水温数值模拟计算，效果不错“6 “3 1 1 。 2 2 3 有限分析法( f i n i t ea n a l y t i cm e t h o d ，简写f a m ) f a m 在1 9 8 0 年由美国i o w a 大学的陈景仁教授等人首先提出，是解析方法与 数值方法结合的一种近似方法，其实质是把数学物理方程的解析方法结合至偏微 分方程的数值解中去，因此有限分析法又叫半解析法。它克服了有限差分法在求 解不可压缩粘性流体时在大雷诺数下的各种困难，避免了差分近似中的各种数值 效应，时求解大雷诺数的各种流体力学问题行之有效的方法。目前，f a m 在求解 n s 方程、地下水运动方面有较为广泛的应用。 f a m 的基本思想是：首先将待求问题的总体区域分为许多小的子区域，在这 些子区域中将方程局部线性化后求解局部的解析解；然后从局部解析解导出一个 代数方程，使子域上的内节点值与相邻的节点值联系起来，利用在求解区域中以 及离散了的子域之问的搭接覆盖，使子域与子域之间平滑过渡，接着把所有的局 部解析解汇集在一起，就得到所求问题的有限分析数值解眦“3 3 1 。 f a m 的优点是：有自动迎风的数值特性，能准确地模拟对流项，不存在数值 扩散现象，计算稳定性好，收敛快；但对于高雷诺数流动往往难以得到精确而稳 定的解。 2 2 4 有限体积法( f i n i t ev o l u m em e t h o d ，简写f v m ) f v m 早在7 0 年代初开始用于平面不可压数值模拟，形成矩形网格上的s i m p l e 类隐式算法。此时，计算气体动力学却采用f e m 网格，沿每个控制体分别建立质 量、动量和能量平衡，得到f v m 方程组，用于求解计算时段末控制体平均的数值 解。f v m 和f e m 一样将计算域划分成若干规则或不规则形状的单元或控制体，可 中山大学硕士学位论文 以很好地拟合边界和水下地形。在计算出通过每个控制体边界沿法向输入( 出) 的流量和动量通量后，对每令控制体分别进行水量和动量平衡计算。便得到计算 时段末各控制体平均水深和流速。因此，f v m 正是对于推导原始微分方程所用控 制体途径的回归，与f d m 和f e m 的数值逼近相比其物理意义更直接明晰哳1 。 设计f v m 格式的关键在于如何计算跨控制体界面的通量。如果采取相邻控制 体形心处通量的平均，便相当于中心格式。此时，格子中心式( c c ) 的f v m 在矩 形网格上，相当于二阶中心f d m 格式，而格子顶点式( c v ) 的f v m 和线性三角形 及双线性四边形单元的g a l e r k i nf e m 等价。但如果采用特征逆风格式计算通量， f v m 便适于处理对流占优的输运问题，且在矩形网格上相当于守恒逆风f d m 格式。 因此，f v m 能像f e m 一样适用于任意的不规则网格。且着眼于控制体上的逼近， 且具有守恒性，又能象特征线格式一样具有以特征为基础的逆风性。并且，在具 有上述优良性能的同时，处理效率与f d m 相近，远高于f e m 。在f v i d 中，根据空 间离散化的形式又可分为中心格式、经典逆风格式、通量向量分裂格式( f v s 格 式) 、利用黎曼解的格式、t v d 格式等。”。 谭维炎、胡四一( 南京水利科学研究院) 采用f v m 成功地进行了长江中游洞 庭湖和钱塘江口涌潮的数值模拟b ；陈阳宇( 珠江水资源保护科学研究所) 也 用f v m 对广卅市南部潮汐河网区进行了水质模拟，不仅精度较高，求解速度也较 快3 6 1 ，蒋艳、杨钰等用有限体积法的o s h e r 格式对长江南京八卦洲江段浓度场 进行模拟计算，验证了模型的合理性和模拟能力n ”。 2 2 5 边界单元法( b o u n d a r ye l e m e n tm e t h o d ，简写b e m ) b e m 是七十年代英国南安普敦大学土木工程系首创的，是一种很有发展潜力 的方法。该方法不需对整个计算区域进行网格剖分，只需区域边界。b e m 的基本 思想是：基于g r e e n 公式和定解问题的g r e e n 函数，将求解的微分方程问题化为 边界积分问题，以边界积分方程为基础，用有限元思想将边界积分方程离散化， 在每个边界单元上将待定函数用其节点值表示，导出近似解的代数方程组。 第2 章水流一水质数学模型及其计算方法 b e m 的优点是：将全区域的计算化为区域边界上的计算，能使问题的维数降 低一维，减小代数方程系数矩阵的阶次，减少计算工作量。b e m 近似范围仅在区 域边界上，而有限元法在全区域单元内采用局部逼近函数，因而边界单元法的精 度一般高于f e m 。实际计算中单元与节点数较f e m 少得多，因而可以合理地加密 网格以提高精度。 b e m 的缺点是：方程组的系数矩阵不对称并为满阵，有时是近似的奇异阵， 求解这类方程组的方法受到限制，计算机内存和计算时间不能缩短，常常还需数 值计算许多积分，对于复杂几何形状问题，本方法不经济“4 “4 ”。 2 2 6 小结 以上对几种主要的水模型计算方法作了简要的介绍和比较。从对各种方法的 对比分析可知，各种方法有其优点，也有其缺点。有限差分法和有限元法提出的 时间早，发展成熟；有限分析法、边界元法和有限体积法还有待完善，发展潜力 大。相比而言，有限差分法对于一维问题具有简便高效的优势；而有限体积法采 用无结构网格的边界拟合技术，离散守恒型的控制方程，在边界拟合、计算精度、 计算速度方面有一定的优势，特别适合求解几何、物理条件复杂的多维问题。根 据以上分析，本文采用有限体积法对平面二维水流水质模型进行离散、求解。 中山大学硕士学位论文 第3 章f v m 在二维水流一水质数值模拟中的应用 本章重点介绍平面二维水流一水质数学模型，并介绍f v m 的基本原理及其求 解过程。通过f v m 的通量向量分裂格式( f v s 格式) 离散控制方程，把二维问题 转化为一系列的一维问题进行求解，不仅提高了数值模拟的精度，而且能模拟包 括恒定、非恒定或急流、缓流等水流状态。 3 1 控制方程 3 1 1 直角坐标系中以原始变量表示的方程 连续方程 丝+ 塑+ 塑：0 a缸砂 宴+ 。_ o u + ，拿+ g 孕： 否州西w 瓦+ g 夏劬z 安+ “却+ ，景+ g 磊o h - g ：钆t 西州+ v 面+ g 面2 q 对流扩散方程 笙+ “笙+ v o c ：b a缸 o y 其中： = 一吉警一g 警+ 字+ ( 3 2 ) ( 3 3 ) ( 3 4 ) 第3 章f v m 在平面二维水流一水质数值模拟中的应用 v = 一上p 挈o y g 等+ 竿+ 7 砂 。 = 昙( e 面0 c ) + 丽o 【髟面o c ) + s戚啦洲鲫 式中：x 、y 、t 分别为空间及时间坐标 u 、v 为x 、y 方向沿水深平均流速分量： p 为水体密度； c 为水体污染物浓度 s 为源汇项；h 为水深 h ，e y 为沿x 、y 方向混和系数 p 为水面大气压力 z b 为河底高程； t a x 、7 a y 为水面风应力的x 分量、y 分量； l 、嘞为水底摩阻应力的x 分量、y 分量 名、为作用在单位质量水体上的体积力( 如科氏力) 的x 分量、 y 分量； 方程组中有关参数的计算公式为 ( 1 ) 风应力和水底摩擦应力 l 的常用公式为：t o = p c 。j 眈l 叱 中山大学硕士学位论文 式中：p 为空气密度；n 0 为水面以上l o m 处风速；c o 为风阻力系数 “在水面无风时可用曼宁公式的二维推广，表达式如下： 铲觥止：呼：巫h 匝i 3 ：硝痧：1 p g v 拓u 孕2 - i - v 2 ：学v 当水面有风时，风应力引起附加水底摩i j h p r 。，系数由经验公式确定。 ( 2 ) 地转力 地转力的公式为( 北半球) f b 。= 务f h = 一如 式中：，为c o r i o l i s 系数，f = 2 m s i n c p ，珊= 7 2 9 x 1 0 r a d s ，即地球自 转角速度，妒为纬度。 ( 3 ) 涡粘项 有时，运动方程还增加一个水平涡粘项，其表达式可以写为： x 方向：鼻v 2 ，y 方向：掌，v2 。，其中量和孝，为x 、y 方向的涡粘系数。 方程( 3 1 ) ( 3 4 ) 不能写成守恒形式，故常用于计算连续流，且此时 动量平衡存在误差。 3 1 2 直角坐标系中守恒形式的方程 1 6 第3 章f v m 在平面二维水流一水质数值模拟中的应用 计算浅水力学中应用较广的一种守恒形式，取 q = ( 厅，h u ，h v ，h c ) 7 = ( h ，q ，q ，q 。) 7 在浅水流动假设下，通过水量和动量平衡导出 连续方程 盟+ 拿( 坠2 + 丝) + 昙单) ： 6 ，：坟 a t缸、h27 却、h 4 等+ 景( 盟h ) + 若萍h + 丝2 v 咄8 ta ) c 、 。 卸、 4 对流扩散方程 等+ 昙c 吼鲁，+ 杀c a ，= v = 如 式( 3 - - 5 ) ( 3 8 ) 可写成向量形式： 塑+ 笪业+ o g ( q ) ：6 a to x加 ( 3 5 ) ( 3 6 ) ( 3 7 ) ( 3 8 ) ( 3 9 ) 以下主要采用这一形式。由此建立的守恒型离散格式，在计算连续流时保证 无守恒误差，在计算间断点时是唯一可用的格式。 其中：q = ( ，h u ，h v ，h c ) 7 = ( ，q ，q ，q 。) 7 胞m 譬+ 譬，半，挚7 = c h u , h u 2 + 譬加加 = 堕砂 望良 槐 研 捏 中山大学硕士学位论文 她h 竿，譬+ 譬，半r ：c h v , h u v , h v 2 + 譬m b = ( o ，g h ( s o 。一s i x ) ，g h ( s o ，一s 办) ，v ( e v ( 向c ) ) 一k c ) 7 s o x - - a 斑z b ，为x 方向的河底底坡比降 & ，：一孥 。 秒，为y 方向的河底底坡比降 f ( q ) 为x 方向通量；g ( q ) 为y 方向通量 = 譬p u 届u 24 - v 2 = 学m 降 耻譬= 1 , o n 2 v q r u 孚2 + v 2m 枷愀降 3 1 。3 定解条件 ( 1 ) 初始条件： l i ( x ，y ，0 ) = v ( x ，y ，0 ) = o c ( x ，y ，o ) = c 。 ( 2 ) 入流边界 入流边界一般取在河道顺直段。 流场边界：z ( x ，y ，t ) ir = z ( t ) ，z 为水位 或：u ( x ，y ，t ) jp = u ( t ) ；v ( x ，y ，t ) jr = v ( t ) ； 第3 章f v i v i 在平面二维水流一水质数值摸拟中的应用 浓度场边界：c ( x ，y ，t ) lr = c ( t ) ( 3 ) 出流边界 出流边界一般也取在比较顺直的河道，要求满足： 且罢= 罢= 祟：0 ( s 为流线方向) a sa s0 s j 。j ( 4 ) 固壁边界 采用不可入条件( 可滑条件) 和固壁物质通量为零条件。 3 2f v m 的原理和基本方程 3 2 。if v m 网格 f v m 采用任意三角形或四边形构成不规则无结构网格，目前倾向于用凸四边 形( 因为节点数相同时三角形网格的格子和边的数目为四边形的两倍或更多，因 而计算量大；拉长的三角形上一阶方法的精度和稳定性很差；用于显格式时，四 边形网格的计算时间步长可以比三角形的大) 。 f v m 二维网格主要有三种布置方式： ( 1 ) 格子中心式( c c ) ：只在每个格子形心处布设一个节点，控制体内的 各变量( 水深、流速等) 为一常数，并定位于形心。 ( 2 ) 格子顶点式( c v ) ：所有流动变量定义在所有格子顶点。 ( 3 ) 混和式：在格子的中心定义水深，在格子顶点定义流速分量。 本文采用格子中心式定义控制体的变量，亦即控制体内各变量为一常数，并 定位于格子形心，见图3 一】。 1 9 中山大学硕士学位论文 口 图3 1f v m 的控制点布置方式 无结构网格的优点是 ( 1 ) 与边界以及水下地形拟合较好，利于边界条件的实现； ( 2 ) 便于控制网格密度，易于修改和适应性调整： ( 3 ) 建网比曲线网格容易，大型三角网、矩形网及混和网可用程序自动生 成。 但无结构网格也有一些缺点： ( 4 )格子排列不规则，需建立适当的数据结构来检索格子问的邻接关系， 占用计算机内存多； ( 5 ) 隐格式求解效率低，一般采用跌代法，数值解后处理工作较大。 3 。2 2f v m 原理及基本方程 f v m 的水流一水质基本方程是对控制体写出的积分形式的物理守恒律。浅水 方程若与对流一扩散方程结合使用，可以用来计算浅水中污染物的输运和扩散。 通常水中介质可视为“被动介质”处理假设其浓度或含量较低，不影响水的密 度和粘性等物理性质，从而也不影响水的流动。此时，水流计算与对流扩散解耦 合，前者所得流场作为后者的输入。 二维浅水方程和对流扩散方程的守恒形式可统一表示为式( 3 9 ) ，采用无 结构网格，将整个计算区域分为许多单元，控制体为x y 平面上的m 边凸多边 形。 定义矩阵，( g ) = ( g ) ，g ( g ) r ，利用散度定理肛v a d c o = 屯a 。d s ，在任意形 n 状的单元n 上对式( 3 9 ) 进行积分并利用散度定理可得有限体积法的基本公式： 2 0 第3 章f v m 在平面二维水流一水质数值模拟中的应用 j j l 挚一掣十掣+ i a b d c o = - l 附舭肛圳， 其中：n = ( n 。f l y ) 1 为单元边界q 的外法线单位向量； f ( g ) n = ，( g ) + g ( g ) h ，为n 方向的通量，以下写为f n ( q ) d w 及d s 为面积分及线积分的微元 q 为控制体平面域( 面积a ) ，勰为其周边( 逆时针方向) ； 在每个单元中q 以常数近似( 即假设单元内的q 为定值) ，因此利用散度定 理，上式的左项和右边第二项可以写成a q t 及a b 。离散后的( 3 - - 1 0 ) 可以 写成： 一堕：一兰彤( g ) l i + a b d t 智1 ( 3 一1 1 ) 其中：l j 为单元第j 边的边长。对于m 边形而言，等式右边第一项表示m 项 之和，其中每一项等于单元各边单位长度上的法向通量与该边长的乘积。 则对控制体i ，可写出显式f v m 方程为： 4 ，( q ? “一g ? ) = , t e e 露( g ) + 4 眈) ( 3 1 2 ) 式( 3 1 2 ) 左边表示控制体内守恒变量在f 内的变化，右边第一项表示沿 各边法向通量之和，第二项表示控制体内源汇项( 入流及外力) h a t 内的作用。 f v m 方程反映了守恒物理量的守恒原理；守恒物理量g = ( 矗，h u ，加，h c ) 7 在控制 体内随时间变化等于各边法向数值通量的和再加上源汇项。不含源汇项时，可以 用图( 3 2 ) 来表示有限体积法中物理量的守恒。 2 1 中山大学硕士学位论文 f 图3 2 有限体积法的守恒原理 将单元边法向通量碟( 口) 简记为f n ( q ) ，为f ( q ) 及g ( q ) 投影到法向的 通量。可利用欧拉方程的旋转不变性，使f v m 方程的计算转化为一维问题的求解。 设n 与x 轴的夹角为舻( 逆时针方向) ，则有： e ( g ) _ m ) ”鼬) b = m ) 尝叫g ) 芸= m ) c o s a + g ( g ) s i n 妒( 3 1 3 ) 上式f n ( q ) 为f ( q ) 和g ( q ) 投影到法向的通量。s p e r k r e i j s e 证明f ( q ) 及g ( q ) 具有旋转不变性，即满足： 丁( 妒) e ( g ) = 厂 r ( 妒) 们= 厂函) 也即f 。( g ) = r ( p ) 一1 厂( 面 ( 3 1 4 ) 式中：丁( 妒) = 10 0c o s 够 0 一s i n l p 0o r ( 妒) = 1o 0c o s 口 0s i n 口 0o oo s i n 妒0 c o s 口0 o1 0 一s m c o s 口 o 第3 章f v m 在平面二维水流一水质数值模拟中的应用 g = r ( 伊) q = ( ，矗“。，h v 。， c ) ，是由q = ( 矗，h u ，h v ，加) 通过坐标旋转变换得到 的，“，和v 。为局部f 一刁坐标系( f 轴为外法向，逆时针旋转9 0 6 得到叩轴) 中 的法向及切向流速分量。因此法向通量f 。( q ) ，亦即f ( q ) 及g ( q ) 在法线的 投影，可转变为先投影q 到q ( 法向) ，然后代入f ( q ) 得到f ( q ) ，再经过一个 坐标逆变换就可以求得f n ( q ) ，由于此二者等效，因此g ( q ) 得以消去，也j 下 是因为这个特性，使原来的二维问题转化为一维问题来处理，即只需计算f ( q ) ， 大大简化了计算并提高效率。 f v m 水流一水质基本方程的最终形式为： 磅= 喜m 。1 而肌舶 ( 3 1 5 ) 本模型的时间离散采用显格式。由上述方程可见，问题归结为如何确定法向 通量，( ；) ，而，函) 可通过求解局部一维黎曼问