（水利水电工程专业论文）黄河宁夏沙坡头典型河段水流运动数值模拟.pdf

中文摘要 捅要 随着计算机技术的广泛应用和发展，对天然河道各方面的研究也更加深入。其 中，以河流动力学、流体力学、数值计算方法等学科知识为基础的水流数值模拟现 今已被纳入进一步研究的范围。对数学模型也已由一维问题逐步发展n - 维及三维 问题。本文主要利用雷诺应力模型，对黄河沙坡头典型河段水流运动进行数值模 拟，并与实测结果进行了比较分析。本文的主要工作有： ( 1 ) 对计算流体力学和水流运动数值模拟的研究现状和存在的问题进行了综述性回 顾，并提出了本文所要研究的主要问题。 ( 2 ) 采用守恒性能高的有限体积法对方程进行离散，并选用符合水流运动特性的混 合格式，利用交错网格下的s i m p l e 算法进行数值计算，提高了计算效率和精 度。 ( 3 ) 对不同工况下非对称突扩通道中的水流运动进行了平面二维数值模拟，得到了 一些结果，这些结果说明了本文数值方法的可靠性。 ( 4 ) 考虑河道的不规则特性，采用贴体坐标变换方法，较好地解决了天然河道不规 则边界的问题。 ( 5 ) 对2 0 0 8 年7 月沙坡头河段的实测资料进行了整理，并作出了分析。 ( 6 ) 对黄河宁夏沙坡头典型河段水流运动进行了三维数值模拟，并与实测结果进行 了对比，结果表明了计算结果的准确性和参数选取的合理性。 关键词：有限体积法，s i m p l e 算法，雷诺应力模型，黄河，沙坡头河段 英文摘要 a b s t r a c t w i t ht h ed e v e l o p m e n ta n de x t e n s i v ea p p l i c a t i o no fc o m p u t e rt e c h n o l o g y , e v e r yr e s p e c t o fn a t u r a lr i v e rr e s e a r c hh a sb e e nm o r ec o n c e r n e d n o w a d a y s ，t h en u m e r i c a ls i m u l a t i o no f r i v e r , w h i c hi sb a s e do nr i v e rf l o wd y n a m i c s ，f l u i dd y n a m i c s ，n u m e r i c a lm e t h o de t c ，h a sb e e n i nt h er a n g eo ff u r t h e rr e s e a r c h t h ed e v e l o p m e n to fm a t h e m a t i c a lm o d e lh a sb e e ne x p e n d e d f r o mo n e d i m e n s i o n a lp r o b l e m st ot w o d i m e n s i o n a lp r o b l e m s ，e v e nt h r e e d i m e n s i o n a lp r o b - l e m s s o m es i m u l a t i o nh a v eb e e nd o n ef o rt h ey e l l o wr i v e rf l o wi nn i n g x i a s h a p o t o ut y p i c a l r e a c hu s i n gr s m ，a n ds o m ec o m p a r i s o nt ot h em e a s u r e dr e s u l t s t h em a i nw o r ki sa sf o l l o w s 1 t h ec u r r e n ts t a t e o f - a r to fc f da n dw a t e r - s e d i m e n tf l o wn u m e r i c a ls i m u l a t i o ni sr e v i e w e da n dd i s c u s s e d s o m em a i np r o b l e m sa r ep r e s e n t e d 2 a d o p tf i n i t ev o l u m em e t h o dw i t hh i g hf l u xc o n s e r v a t i o nt od i s c r e t et h ep a r t i a ld i f - f e r e n t i a le q u a t i o n ，a n dc h o s et h eh y b r i ds c h e m eo fo n e o r d e l w h i c hc o n f o r mt ot h ef l u e n t m o v e m e n tp r o p e r t y t h es i m p l ea l g o r i t h m ，w h i c hb a s e do nb a c k w a r ds t a g g e r e dg r i d ，i s u s e df o rn u m e r i c a l c o m p u t a t i o n ，t og e th i g h e re f f i c i e n c ya n da c c u r a c y 3 t h eu n s y m m e t r i c a ls u d d e ne x p a n s i o n sc h a n n e lf l o wi ss i m u l a t e ds e v e r a ld i f f e r e n t s i t u a t i o n s t h es i m u l a t i o nr e s u l t ss h o wt h a tt h em o d e la n dn u m e r i c a lm e t h o d sa r er e l i a b l e 4 c o n s i d e r i n gt h ei r r e g u l a rp r o p e r t yo ff i v e r ，w eu s et h eb o d y f i t t e dc o o r d i n a t et r a n s - f o r m a t i o nm e t h o dt og e tb e t t e rs o l u t i o nt ot h en a t u r a lr i v e r i r r e g u l a rb o u n d a r i e s 5 w eh a v ed o n es o m er e a le x p e r i m e n t sa n dm e a s u r e m e n t ss i n c e2 0 0 8 6 s o m es i m u l a t i o nh a v eb e e nd o n ef o rt h ey e l l o wr i v e rf l o wi nn i n g x i as h a p o t o ut y p i c a lr e a c h ，a n ds o m ec o m p a r i s o nt ot h em e a s u r e dr e s u l t st oi l l u s t r a t et h ec o r r e c t n e s sa n dt h e r e a s o n a b i l i t yo ft h en u m e r i c a ls i m u l a t i o n k e yw o r d s ： f i n i t ev o l u m em e t h o d ，s i m p l ea l g o r i t h m ，r e y n o l d ss t r e s sm o d e l ，t h ey e l l o w r i v e r ，s h a p o t o u 一一 独创性声明 本人郑重声明：所呈交的学位论文，是本人在导师的指导下，独立进行研究工 作所取得的成果。除文中已经注明引用的内容外，本论文不含任何其他个人或集体 已经发表或撰写过的作品或成果。对本文的研究做出重要贡献的个人和集体，均已 在文中以明确方式标明。本声明的法律结果由本人承担。 论文作者签名：j 蹲至童l 日期：么甲年月j e t 关于学位论文使用授权的说明 本人完全了解宁夏大学学有关保留、使用学位论文的规定，即：学校有权保留 送交论文的复印件和磁盘，允许论文被查阅和借阅，可以采用影印、缩印或扫描等 复制手段保存、汇编学位论文。同意宁夏大学可以用不同方式在不同媒体上发表、 传播学位论文的全部或部分内容。 ( 保密的论文在解密后应遵守此规定) 学位论文作者签名：载歪盛 日期： 作者指导教师签名： 辜圭乞 日期： 训降6 月j e t b7 年月午日 一v 一 第一章绪论 1 1 选题的目的和意义 第一章绪论 水利枢纽的泄洪、排沙、发电和供水等设施以及进水口位置的设定与运行，都 要求对坝区水流运动有正确的了解。如果对坝区水流运动不能正确地预测就可能会 妨碍枢纽的防洪、灌溉、发电和航运等各种功能的发挥，有时甚至会威胁到人民生 命财产的安全。 对河道水流问题的研究有三种方法：用理论分析、物理模型试验和数值模拟计 算。理论分析是对大量原型观测资料进行定性和定量分析的基础上，根据所研究水 域的水文特征，利用水力学及河流动力学和相关学科理论，探讨河道流量水位的变 化规律，为水利工程的建设和运行方式提供科学依据。物理模型试验依据水流运动 学和动力学方程，建立一定的相似准则，将原型河流缩制成模型河流，在模型上进 行水流运动研究，预测拟建工程的作用与效果，以及工程对附近水流的影响，再依 据模型相似率推广到原型，为工程设计和运营提供科学依据。数值模拟计算是运用 一定的离散方法，数值求解水流运动方程，得出工程水域在设计水文情况下的水流 情况，预测拟建工程的作用与效果，为工程的前期规划和多方案比选提供依据。 理论分析方法优点在于所得结果具有普遍性，各种影响因素清晰可见，是指导 实验研究和验证新的数值计算方法的理论基础。但是它往往要求对计算对象进行抽 象和简化，才有可能得出理论解。 物理模型对一些复杂的问题，如局部急变水流，复杂边界条件等问题，可以提 供直观的观察。但物理模型其耗资和比尺带来的问题都往往使其在应用和扩展的范 围上受到限制，因此物理模型研究的对象往往是局部性的，对于人范围的工程规 划，就很难较准确地再现和重演物理现象的历史，适应性很差。 数值模拟计算的方法恰好克服了前两种方法的弱点，在计算机上实现一个特定 的计算，就好像在计算机上做了一次物理实验。数值模拟的优越性在于能开展长时 期、长河段河流水文情况的研究和预报，并具有周期短、投资少的优势f l l 。 1 2 河流数学模型的研究现状与发展方向 河流数学模型的前期工作始于5 0 年代，当时的计算工具运算速度和存储量都受 到很大的限制，在基本方程和计算方法上做了较多的假设和简化。现代河流数学模 型是建立在计算流体力学的基础上，伴随着电子计算机突飞猛进的发展，于7 0 年代 开始迅速发展起来f 2 1 。7 0 年代以来，在世界范围内，河流数学模型发展较迅速，其 中以美国、荷兰和法国发展最快，已研制出众多的一维模型。在我国，许协庆和朱 一l 一 宁夏大学硕十学位论文 鹏程、窦国仁、韩其为等人在6 0 年代和7 0 年代初期就开始了河流数学模型的前期工 作。目前，一、二维泥河流数学模型在理论上都得到较快的发展，在实际工程应用 中正越来越广泛得到应用。 一维模型模拟的是变量沿河长方向的平均值，基本上能够满足实际工程需要， 它是至今使用最广泛的一种模型。该模型引进一些假定后，沿断面积分三维n s 方程 便可得到一维水流运动方程和连续方程，即可进行一维计算。一维水流数学模型一 般采用差分法和特征线法进行数值计算，其基本思路是：将计算河段划分成若干小 河段，计算各断面的平均水力沿程变化及因时变化情况。一维水流计算中的关键之 一是糙率的确定，大多通过反求水面线的方法加以确定。 一维数模在理论及实践上都比较成熟，国内外使用已经很普遍，主要用于研究 长河段长时期的河道水流。一维数模以其经济灵活的优势见长，已经可以部分取代 物理模型实验，同时，经常一、二维耦合计算，利用一维结果为二维提供边界条 件。但是在实际应用中，如想了解局部河段的水力情况一维模型就显得无能为力： 必须用平面二维甚至三维数模才能解决。 二维水沙数学模型是从研究河口、海岸水流泥沙运动开始的，h a n s o n w 1 3 1 最早 进行了这一研究。二维模型克服了一维模型不能计算河流细部变化的缺点而得到迅 速发展，目前在工程中得到了较为广泛的应用，正逐步走向成熟，视工程的重要 性，经常与物理模型并用，并能部分替代物理模型试验。 三维水沙数学模型是近十几年才发展起来的，但发展很快。已有的三维泥河流 数学模型，大多数采用了静水压的假定。国内外也有一些学者在一定的简化情况 下，采用了近似的三维或准三维河流模型，研究实际工程水力问题。 在三维泥沙问题中，还存在一些计算上的难度。其一，某一时刻的流态取决于 初始条件、边界条件、时间相过程以及自由表面的波动等，这些因素的组合使得实 际流态表现为十分复杂的非恒定现象。其二，水流速度的三个分量随时间而改变。 其三，水流流速的梯度的空间不均匀性，且水流的水平尺度往往大于垂向尺度，使 得数值离散后网格的水平尺度远大于垂线尺度。在水流数值计算中，床面附近流速 计算是否准确，直接影响着泥沙的底床条件的准确与否。另外，由于水流速度在水 平方向和垂直方向上变化的数量级是不同的，因此数值计算过程在不同的方向上应 选用不同的网格尺寸。 河流数学模型涉及面颇广，目前国内外学者对于河流数学模型的研究虽己取得 了丰硕的成果，但是应该看到，仍有许多问题没有解决。其发展趋势就目前而言应 朝着以下几方面发展【4 1 。 ( 1 ) 高精度高效率的模型研究天然河流的水文条件具有极大的随机性，对于 这种复杂的水体变化过程要用一个精确并规则的模型去模拟需要进行很多理想化的 一2 一 第一章绪论 处理，这种近似并不能真实的反映水体变化的过程，所以提高模型的精度及效率是 极其重要的f 5 1 。 ( 2 ) 数值计算方法的改进各种各样的数学模型在应用和推广中往往受到数 值方法的限制，我们所要求的最基本的是数值方法的相容性、稳定性和收敛性。如 今，它应朝着守恒性、无震荡性、无弥散性、减少数值扩散等方向发展。因此寻找 一种集稳定性、相容性、精确性、健全性、通用性和经济性于一身的理想数值计算 格式是理论界和工程界所追求的目标。 ( 3 ) 模型的可视化计算结果的图像化动态显示已经成为流体计算结果后处理 的一个新的组成部分。它使计算结果显明、直观、生动。网格的自动生成、初边值 条件的自动给定以及与c a d ，g i s 等技术相结合是河流数学模型发展的一个重要趋 势。 ( 4 ) 高维模型的研究无论是一维还是二维模型，都在很大程度上对水流及泥 沙状态进行概化和近似，要想真实反映河流模型的优势，只有深入研究空间运动规 律，所以未来河流模型的发展趋势仍是三维f 6 1 。 1 3 数值计算方法 根据物理现象所建立的数学模型中，其基本方程大都是非线性的偏微分方程， 对于这些方程一般很难求出解析解，所以只能通过对控制方程的离散求解代数方程 组来得到近似的数值离散解。目前，常用的求解数学模型的方法主要有：特征线 法、有限差分法、有限元法、有限分析法有限体积法及等。 l 、特征线法2 0 世纪5 0 年代初该法由林秉南首先提出，是解双曲型偏微分方程 的最精确的数值方法。其优点在于它是一种更合理的迎风格式，具有优良的精度。 适用于解决短周期、变化急剧( 如涌潮) 的一维问题。其不足在于，特征方程常为非守 恒形式，用差分法离散特征方程时会带来守恒误差。而且也不能直接计算间断解。 2 、有限差分法它是计算机数值模拟采用的最早的方法，也是目前应用最为 广泛的一种计算方法。其优势在于处理问题效率高，方法成熟。缺点是对不规则边 界的适应性不强，有时不能严格保持守恒性质，而且常常也不能用来正确计算间断 解。在二维情形，由于使用t a y l o r 级数展开，有限差分法一般多用于矩形或正交曲线 网格，通常在计算域概化及数值解精度方面存在困难。 3 、有限元法有限元法是r c o u r a n t 于1 9 4 3 年首先提出的，7 0 年代后期开始应用 于水沙数学模型中。有限元法的优越性在于弥补了差分方法较难适应区域形状的任 意性的不足，计算稳定性好，精度高1 7 ，8 1 。但其存在计算格式复杂，计算量及存储量 大，尤其在多维计算中直接影响计算速度。 4 、有限分析法有限分析法是美籍华人陈景仁提出来的，该方法计算精度高， 一3 一 宁夏大学硕： ? 学位论文 具有自动迎风特性，计算稳定性好，收敛较快。但由于有限分析系数中含有交错级 数，这给实际计算与理论分析都代来了一些困难。虽然后来李炜等提出了混合有限 分析法，避免了无穷级数带来的不便，但仍存在着单元系数较复杂，计算速度比较 慢等缺点。 5 、有限体积法有限体积法1 9 l 又称有限控制容积积分法，它与有限元一样，将 计算域划分成若干规则或不规则形状的单元或控制体，而不像经典有限差分法那样 要求网格有结构的排列，因而处理过程具有较强的灵活性，能够满足复杂边界问题 的需要。目前有限体积法网格与有限元法网格( 主要指三角形网格) 一样，可由计算 机程序来辅助生成。其基本思想是把计算区域离散为若干点，以这些点为中心，把 整个计算区域划分为若干互相连接但不重叠的控制体。在计算出通过每个控制体边 界沿法向输入或输出的流量和通量后，对每个控制体分别进行水量和动量的平衡计 算，便得到计算时段末各控制体的平均水深和流速。因为跨控制体间界面运输的通 量对相邻控制体来说，大小相等，方向相反，故对整个计算区域沿所有内部边界的 通量相互抵消，对由一个或多个控制体组成的区域，以至整个计算区域，都严格满 足物理守恒定律，不存在守恒量的误差，并能正确计算间断点。 设计有限体积法格式的关键在于如何计算跨控制体界面的通量。若采取相邻控 制体形心处通量平均，在矩形网格计算中，便相当于二阶中心差分，而在三角形或 四边形网格中，若物理量定义在网格顶点，则又与线性三角形和双线性四边形单元 的g a l e r k i n 有限元法等价。若采用特征逆风格式计算通量，有限体积法适丁- 处理对流 占优的输运问题，且在矩形网格上相当于守恒逆风有限差分格式。由于有限体积法 从物理规律出发，每一离散方程都是有限大小体积上的某物理量的守恒表达式，在 推导过程中物理概念清晰，并可以保证离散方程的守恒特性，因此，它能象有限元 一样适用于任意的不规则网格，且着眼于控制体上的逼近，具有守恒性，又能象特 征法格式一样具有以特征为基础的逆风性。并且，处理效率与有限差分相近，而远 高于有限元法。可以说，有限体积发体现了有限元的灵活性、有限差分的效率和守 恒性以及特征发的精度，在数值模拟中有着丰富的发展潜力。 1 4 模型中存在的主要问题 1 、不规则边界的处理 天然河道的边界一般都是极其复杂的，从而使得一些数值格式的应用受到了很 大限制。对此，很多学者做了很多工作，应用比较多的是采用贴体坐标系计算网格 的生成方法，它使计算网格与计算域边界相重合，从而使边界条件能准确地引用到 计算网格节点上，并使得变化后的网格成为规则的矩形网格f l o l 。但数值网格生成 方法计算比较繁琐，变换后的基本方程也相当复杂，因而至今未得到广泛应用。此 后，喻国良于1 9 8 7 年提出了镜像法，周建军于1 9 8 8 年提出了边界对称点法。这两种 一4 一 第一章绪论 方法都是根据对称原理，简单易行，提高了不规则边界上运用差分方法的计算精 度，但由于域内外对称点一般不在节点上，仍需一定的插值过程，因此会带来一定 的误差。 2 、动边界的处理 由于河床冲淤变化和水位的波动，在河岸及洲滩等区域都会遇到动边界问题。 动边界数值模拟的困难主要在于沿动边界法向的流动不同于明渠均匀流，非恒定与 非均匀性强，常用的曼宁摩阻公式等在形式上难以套用，再者因水深很小，对离散格 式的要求很高，要求数值解不产生假振，保证水深总是大于零等，这给数值模拟带 来了一定的难度。 3 、阻力问题 阻力计算是阻力系数的确定问题。阻力系数是反映水流条件和河床形态的综合 系数，其合理与否直接影响到水力计算的精度，进而影响到含沙量及河床变形的计 算结果。对一维阻力问题，处理方法有两种，一是采用水力半径分割法或能坡分割 法划分阻力单元，分别计算各单元的阻力，然后再计算总阻力；另一种方法是根据 实测资料直接计算总阻力【1 。这两种方法只能用于一维水流的糙率计算，对二维水流 的研究，目前研究相对较少，其困难在于比降沿河宽分布资料很难取得。 4 、模型验证问题 一个数学模型能否有效地发挥作用，丰要在于能否给出具有足够精度的定量结 果，为此，模型的验证工作至关重要。一个模型如果没有天然河流实际资料的验 证，一般不能用来进行工程计算的。在选取验证资料时，要求所采用的资料与所研 究的问题属同一类型，并且有相似的地理水文条件，验证所需要覆盖的空间和时间 也应大体卜相互匹配。 5 、算法方面的问题 控制方程中均含有一些高度非线性的参数，在数值求解时容易发生不守恒、不 稳定和不收敛现象。 1 5 本课题所要解决的问题及拟采用的研究方法 1 5 1 选题基本内容概述 本课题研究的是二维和局部弯道处的三维水流数学模型，模拟的地段为黄河沙 坡头水利枢纽典型弯道段。在黄河宁夏段，这段河流具有代表性。建于该河段的黄 河宁夏沙坡头水利枢纽是一个以灌溉、发电为主的综合水利、水电工程。该枢纽距 上游拟建的大柳树水利枢纽1 2 1 k m ，距下游已建成的青铜峡水利枢纽1 2 2 k i n 。工程 于2 0 0 2 年1 2 月2 6 日开工，2 0 0 5 年5 月3 0 日最后一台机组发电，标志着枢纽主体工程建 一5 一 宁夏大学硕十学位论文 设任务的完成。沙坡头水利枢纽的修建，改原来的无坝引水为有坝引水，改善了沙 坡头灌区的供水水源，为灌区节水灌溉和分引黄河水创造了条件，因此该工程是一 个有着重要的社会价值、经济价值的综合工程。沙坡头水利枢纽投入运行以后，改 变了河床情况、来水、来沙条件等，对该河段和下游青铜峡库区的水流运动、泥沙 运动、及未来的河床演变带来深刻影响。 本课题以黄河沙坡头水利枢纽坝前桩号s 日1 0 至s 日4 段研究对象，用数值模拟的 方法探求该段水流运动规律，为该段水库水力预报和水库的合理运行起到一定的参 考作用。 1 5 2 研究的问题及方法 由于该领域存在着远不止上文所罗列的一些问题，而这些问题又需要通过深入 的探讨和研究去逐一的攻克和解决。因此，本课题针对其中的几个问题在展开深入 的研究，将要解决的问题呈列如下： l 、模型选择方面，本课题针对黄河沙坡头河段的特点，选用雷诺应力模型。 2 、针对区域的不规则性，本课题采用贴体坐标变换方法，将不规则的边界转化 成规则的计算区域，通过相应的变换控制方程的形式，在计算平面中对控制方程进 行离散求解，最后将求得的结果经逆变换变为物理平面上的解，于是从根本上解决 了处理不规则边界问题。 3 、本课题对于确定水动力参数问题，针对黄河的实际情况采用经验法，在模型 验证阶段对照实测资料对参数进行调试，使水力与实际情况相符。 4 、在计算方法上，采用物理守恒性能较好的有限体积法对控制方程进行离散， 并用压力耦合类s i m p l e 算法进行代数方程组的求解。 5 、在模型验证上借助北方民族大学数值计算与工程应用研究所的高科技设备实 地测取数据进行验证。 一6 一 第_ 章控制方祥 2 1 概述 第二章控制方程 河流运动属于湍流运动，遵守流体运动的基本方程。如果直接对其进行求解， 要求空间步长和时间步长特别小，目前无法用于真正意义上的工程计算。其实在水 利工程中一般对湍流的平均效应感兴趣，因此雷诺用时均法对湍流进行了研究，提 出了三维湍流时均方程。但是湍流时均方程不能封闭，必须引入新的湍流模型，本 文采用了雷诺应力方程模型。 2 2 河流运动的基本方程 2 2 1 质量守恒方程 河流运动必须满足质量守恒定律，根据这一定律可以得出质量守恒方程1 1 2 】： 警+ 掣+ 掣+ 掣= o 协1 ， 出如鼬舭 。 r 一7 式( 2 1 ) 口- - j p f l , 写成如下形式： 警+ d i v ( p u ) = o 对于不可压河流，密度p 为常数，式( 2 1 ) 变为： 警+ 等+ 警= 。 对于定常流动，密度p 为不随时f , - 变g a ，式( 2 1 ) 变为： 掣+ 掣+ 掣：o 如。 鼬 以 。 质量守恒方程又名连续性方程，习惯上使用连续性方程这个名字。 2 2 2 动量守恒方程 ( 2 2 ) ( 2 3 ) ( 2 4 ) 动量守恒定律也是河流运动必须满足的基本定律，根据这一定律可以得出z 、y 7 一 宁夏大学硕十学何论文 和z 三个方向的动量守恒方程i 1 2 1 ： 掣+ 警碧+ 警铲+ 百o ( p u w ) = 瓦0 ( p 瓦o u ) + 南( p 万o u ) + 未( p 瓦o u ) 一瓦o p + & ( 2 5 n ) 掣+ 掣+ 掣+ 掣= 瓦0c p 务瓦0c p 争瓣，一私5 6 ， 掣+ 警掣+ 掣+ 亟掣= 瓦0 ( p 筹) + 南( p 面o w ) + 南( p 警) 一瓦o p + & ( 2 5 c ) 式中，p 是动力粘度，瓯、& 和是动量守恒方程的广义源项，鼠= 足+ s $ ，鼠= 蜀+ 匈，& = 只+ ，足、日和e 是微元体上的体力，& 、黾和的表 达式如下： = 岳( p 瓦o u ) + 南宝) + 毫( p 瓦o w ) + 晏( a d i v 即) ( 2 6 n ) s = 丕o 口u ) + 瓦0 咧o v ) + 麦o 咧w ) + 南( m i v h ) ( 2 6 6 ) = 昙( p 瓦o u ) + 晏( p 塞) + 瓦0 ( p 瓦o w ) + 曼( a d i v 口) ( 2 6 c ) 2 瓦( p 瓦) + 瓦( p 瓦) + 瓦( p 瓦) + 瓦( a d i v 叫 ( 2 6 c ) 式中，a 是第二粘度，一般可取a = - 2 3 13 1 。、8 y 和s ：是小量，对于粘性为常 数的不可压水流，& = 卸= 8 ：= 0 。 方程( 2 6 ) n - - 以写成如下形式： 一8 一 百o ( p u ) + d i v ( 胛俄) = 出v ( p 酽a du ) 一赛+ & ( 2 7 口) 掣+ 州川= 酬p 口酬一塞+ & ( 2 7 6 ) 掣+ d i v ( 掣口) = 出v ( pg r a d 叫) 一宝+ ( 2 7 c ) 第_ 章控制方程 式( 2 6 ) 与( 2 7 ) 是动量守恒方程，简称动量方程，也称作运动方程，还称 为n a v i e r - s t o k e s 方程。 2 3 湍流时均方程 实际河道水流一般都属于湍流，具有随机的脉动性。雷诺采用时间平均法， 将各运动要素的瞬时值分解为时均值和脉动值两部分，o j a = 才+ ，并认为方 程( 2 1 ) 和( 2 6 ) 仍然适用于湍流的瞬时运动。在水利工程中，一般对湍流的平均效应 ( 如时均流速、时均压力等) 感兴趣，因此雷诺采用时均法对湍流进行了研究，推 出了湍流时均运动方程。 设某一运动要素的瞬时值为a ，时均值为万，脉动值为a 7 ，时均值的定义为 11 a = 击a d t ( 2 8 )工_ j、， 将以= 瓦+ 口7 ；u = 西+ u t ；口= 万+ 口7 ；w = 面+ 叫7 ；p = 矽+ p ，代入连续性方 程( 2 1 ) 与动量方程( 2 7 ) ，并对时间取平均得到湍流时均流动的控制方程如下： d i v 西= 0 ( 2 9 ) 篆删咖一去塞+ vd i v ( 倒- ) + _ 喾一塑o y 一挈】( 2 1 0 a ) 襄删咖一丢雾+ v d i v ( g r a dv ) + - 警一雩一挈】( 2 1 0 b ) 署抽( 而) _ _ 丢笔+ vd i v ( g r a d _ ) + 【- 警一孕一喾】( 2 1 0 c ) 方程( 2 6 ) 是时均形式的湍流连续方程，方程( 2 7 ) 是时均形式的湍流动量方程， 是雷诺1 8 9 5 年导出的，被叫做雷诺时均n a v i e r - s t o k e s 方程，常直接称为雷诺方程。 与n a v i e r - s t o k e s 方程相比，雷诺方程多出了雷诺应力项一面瓦，3 个正应力和3 个切应 力。 式( 2 9 ) 与( 2 1 0 ) 构成的方程组共有4 个方程，现在共有1 0 个未知量 ( u 、钐、伽、p 和6 个雷诺应力项) ，方程组不能封闭，因此必须引入新的湍流方 程( 模型) 。 一9 一 宁夏大学硕十学化论文 2 4 雷诺应力方程模型 目前湍流模型有两大类，雷诺应力模型和涡粘模型。雷诺应力模型又分为雷诺 应力方程模型和代数应力方程模型，本文在模拟天然河道中水流时选用的是前者。 雷诺应力模型由l1 个偏微分方程组成，其中雷诺时均连续性方程1 个，雷诺时均 动量方程3 个，雷诺应力输运方程5 个，湍动能耗散率输运方程1 个对充分发展的湍 流来说，稳态的连续性方程和动量方程如下： 竽：0( 2 1 7 ) ， 丐鬟= 一丢差+ 去( o - 瓦7 0 x j 一弼m = l 2 ，3 ( 2 m ) 这里x l ，z 2 和z 2 分别为纵向、横向和垂向流速，面为方向的时均流速，多为 时均压强，为水的运动粘滞系数，弼为雷诺应力，g 为重力加速度雷诺应力 ( = 玛) 可通过求解以下雷诺应力输运方程得到： w 功_ o p h 七j 一屯惫篆+ 蜀七象) + 嘞慨 ，2 3 ( 2 1 9 ) 这里d i j 为雷诺应力如由于扩散作用造成的输运项，m e l l o r 和h e r r i n g 提出以下 方程： 巩= g 去 譬( 甏+ 卺+ 等) 】 ( 2 2 。) 这里a 为常数( ：孚) ，k = r 2 i 是湍动能雷诺应力的耗散率如下式所示： = 云妨 ( 2 2 1 ) 这里e 是湍动能耗散率，如为k r o n e c k e r 算子符号g 可用以下方程求得： _ 两o e = 瓦0 ( q i k 蝴两o e ) + 丢( g t m g 2 ) ( 2 2 2 ) 这鼽为湍动能产生率，g = o 1 8 ，g l = 1 4 5 ，g 2 = 1 9 0 雷诺应力输运方程的晟后一项为压力应力项，其功能为在雷诺应力中重新分配 一1 0 一 第_ 章控制方稗 湍动能s p e z i a l ee ta 1 f 1 4 】提出以下公式： 巧= o t o e b i j + a l e ( b 诀+ b j k 一言k m 6 n m 如) + q 2 七岛+ o l 3 p k b i j u n + o q k ( b i k s j 七十b j k s 诀一丢6 斛鲰z 妨) + a s k ( b i k 七十b j k w i k ) u ( 2 2 3 ) 这里为反对称张量，为应力变化率张量，为旋转张量，咖= 一3 4 ，q 1 = 4 2 q 2 = 黑屯3 坠n 煳= - 1 8 , q 4 = 1 2 5 ，q 5 = o 5 幻= 鲁一吾如， = 丢( 筹+ 篆) ，= 互1l 面o 砺- i o 。豇z t j i ) “ 1 l 一 宁夏大学硕十学位论文 3 1 概述 第三章方程离散与求解 水流控制方程是微分方程，很难获得解析解，因此需要将方程在一定的网格系 统内离散，用网格节点处的场变量值近似描述微分方程中各项所表示的数学关系， 按一定的物理规律和数学原理构造和微分方程相关的离散代数方程组。引入边界条 件后求解离散代数方程组，得到各网格节点处的场变量分布，用这一离散的场变量 分布近似代替原微分方程的解析解。当前求解河流微分方程的计算方法有很多种， 如有限差分法、有限元法、特征线法和有限体积法等等。本文采用守恒性能高的有 限体积法。为了方便起见，本文介绍了求解二维问题的有限体积法，其思路可以容 易地推广到三维问题。 3 2 基于有限体积法的方程离散 3 2 1 有限体积法 有限体积法( f i n i t ev o l u m em e t h o d ，f v m ) 又称为控制体积法，最早是 由p a t a n k a r 在h a r l a w 和w e l c h 提出的“交错网格”技术的基础上进行研究的，进而 发展成为一套用于数值传热学和流体流动的计算方法【1 5 1 它已经成功解决了传热、传 质、燃烧以及流体流动过程的许多实际问题。 有限体积法( f v m ) 又称为控制体积法，其基本思路是：将计算区域划分成若干个 互不重叠的控制体，每个控制体包含一个计算网格点如图( 3 1 ) 所示，然后将控制方 程在每个控制体上进行积分，从而得到一组离散方程。有限体积法得出的离散方程 要求因变量在控制体积内积分守恒，因此得到的结果在整个计算区域内，所求物理 量局满足的积分守恒由于有限体积法是根据物质守恒的基本原理推导得出的，因 一1 2 一 。并 一 篁 冒 崔。勿绣髟笏el量 钐毵钐钐弗 = s 图3 1二维控制体示意图 第二章方稗离散与求解 此具有以下不可替代的优越性首先，具有良好的积分守恒性，离散方程的解表示的 变量如质量、动量等在任何数量的控制体积上的总体积分是守恒的，当然在整个计 算区域也是守恒的。这一点使得有限体积法不仅能在精细的网格上能获得很好的结 果，而且，即使在较粗的网格的情况下，也可以得到较好的结果。其次，有限体积 法具有计算过程稳定的特点，在计算网格不变的情况下，适当加大时间步长，只影 响精度，不影响计算过程的收敛性，便于程序的调试。 利用有限体积法建立离散方程时，必须遵循如下网项基本原则： ( 1 ) 控制体积界面上的连续性原则，当一个面为相邻两个控制体共有时，在这两个 控制体的离散方程中，通过该界面的通量表达式必须相同。 ( 2 ) 正系数原则，在大多数的流体动力学问题中，节点上的凶变量的数值通过对流 过程和扩散过程受到相邻点的影响，中心节点的系数a p 和相邻节点系数6 必须 恒为正值。 ( 3 ) 源项的负斜率线性化原则。 ( 4 ) 系数a ，等于相邻节点系数之和原则。 定理3 1 ：设空间二维单连通区域y 的边界曲面s 是光滑的，函数p ( x ，y ，z ) ，q ( x ，y ，z ) ， 兄( z ，y ，z ) 在y 及s 上具有关于x ，y ，z 的连续偏导数，则有： ( 筹+ 筹+ 蓑) d y = p 如+ q c o s ( n , y ) + 冗如粥 = l | p d y d z + q d x d z + r d x d y 其中：n 为曲面s 的外法线方向 证明见文献【1 6 1 上述定理也称为高斯散度定理，可表示成张量形式1 1 7 1 f vd i v ( 州y = zn a d s = z 岫d s c 3 m 对二维对流扩散问题控制方程可表示成通用形式： 塑咝+幽+警笋=旦(r00)ot o xo xo x + 南( r 等) + & 。 a 剪 、 a ”、a ”，。9 其张量形式： 掣+ 酬胛沪酬啊口d e ) + & ( 3 2 a ) ( 3 2 b ) 一1 3 一 宁夏大学硕十学位论文 根据如图( 3 1 ) 所示的计算网格，式( 3 2 a ) 在源项中分离出压力梯度项的形式在控制 体p 及时间段进行积分有： 掣+ 掣+ 等= 一塞+ 丢c r 鬈，+ 南c r 器，+ & c 3 渤 下面分别对式( 3 3 ) 各项进行积分计算 1 瞬态项 在处理瞬态项时，假定物理量咖在控制体p 上的值均为节点处的值，同时假定 密度p 在时间段t 内的变化量极小，因此有： 件越r n 掣蛐出叫沪霹) 龇可 4 ， 2 对流项 利用g u a s s 散度定理，将体积分化为面积分后得： 件越r n 掣捌础= 也七m 切 舭亡 汁& r n 掣蛐出= 九七m 。陋亡 ( 3 5 a ) ( 3 5 b ) 3 扩散项处理 件越r n 嘉( r 鬈) 如由出= ( r 鬈) 。一( r 鬈) 钮 却龀 c 3 6 a ， 绅r 挖丢( r o 则) d x 咖出= ( r 篙) n 一( r 器) 。 z & c 3 舳， 在计算流体力学中扩散项的处理通常采用中心差分格式，对于一般的流动i ；- 1 题这一 离散格式能够满足要求f 1 8 1 ，这样既保证了问题的精度，同时离散形式也比较简单，易 于计算。 4 源项 源项处理通常的做法是局部线性化处理，即假定在未知量微小的变动范围内， 源项s 可以表示为该未知量的线性函数，在控制体p 内表示为：s = & + 品咖p 其 中& 为常数部分，& 是s 随变化的曲线在尸点的斜率，一般要求品0 蚪础r 住s 蛐出娟+ ) 龇础 1 4 第三章方稗离散与求解 5 压力项的离散 一厂r n ( 塞) = 厂( m 刊舭 此处假定在控制容积的东、西界面上压力使各自均匀的，分别为仇，m 得到方程( 3 3 ) 的各项表达式后，只需要 ( 1 ) 在对流项中引入特定的离散格式将( 3 5 a ) ( 3 5 b ) 中界面物理量九，九，和九用节 点物理量来表示。 ( 2 ) 在对流项、扩散项和源项中引入全隐式的时问积分，如+ 出九d t = 九a t 可得到全隐式时间积分方案下的维瞬态对流扩散问题的离散方程： 口p p = o 缈w + n e e + 口 r + 口s s + b ( 3 7 ) 其中系数o w ，a f ，a n 取决于在对流项中引入的离散格式。 3 2 - 2 对流项离散格式 由上节分析，由于时间项采用的是全隐式时间积分方案，瞬态流动与稳态流动 的差别仅仅是对瞬态项的处理。因此为了计算简便先对稳态二维不可压缩对流扩散 方程进行分析 o ( 爰j + 百o ( p v o ) = 丕( r 鬈) + 瓦0 ( r o 口0 ) - - i l - 品 ( 3 8 ) 。十一= 一i l i 十一l i。 - a 可务z 、a z 各芏，、口 u 9 、， 对对流换热控制方程进行数值求解的关键就是对两个一阶导数项，非线性对流项及 动量方程中的压力梯度项的处理。不可压缩流场的数值求解中的关键问题都是由这 两个一阶导数项的离散所引起的。非线性对流项的处理涉及到对流项的离散格式问 题，而动量方程中的压力梯度项的处理则关系到压力与速度间耦合关系的问题。 总通量密度j 是指单位时间内、单位面积上由扩散及对流作用而引起的某一物理 量总的转移量。对于通用变量砂，总通量密度为： 五= 舢一r 鬈山= p v o r 箬 ( 3 9 ) 控制方程的通用表达式( 3 8 ) 可变形为： 篙+ 筹= 岛 ( 3 1 0 ) 如鼬 唧一7 将( 3 1 0 ) 对图( 3 2 ) 表示的控制体进行积分： 一1 5 一 宁夏大学硕十学伊论文 图3 2p 点控制体 五一l + 厶一以= ( s c + 品“) a z a y ( 3 1 1 ) 其中：莎p 为控制体中p 点处的多值，以，l ，以，厶为各自控制体上积分的总通量。同 样对连续性方程积分得： 尼一r + r 一只= 0( 3 1 2 ) 式中： e = ( 肚) 。a y r = ( ) n a x 咒= ( 膨) a y e = ( ) 。a x f 为对流质量流量，表示通过界面上单位面积的对流质量通量。 式( 3 1 2 ) c p 一( 3 1 0 ) 得： ( 以一尻蚱) 一( l 一凡九) + ( 厶一r p ) 一( 以一只如) = ( & + 品p ) z 秒( 3 1 3 ) 当对流项采用p a t a n k a r 的乘方格式时： a ( i p ) = m a z ( o ，( 1 0 1 i p l ) 5 ) fr 其中：p2 茜称为p e c l e t 数，d = 二5 x 为界面上的扩散导数。得到二维通用微分方程 的离散格式( 3 7 ) 其中： 一1 6 一 口e = d 。a ( i p d ) + m a x ( - f e ，0 ) a = 队a ( j r i ) + m a z ( - f n ，0 ) a p = a e - t - a w - t - a n + a s s p a x a y a w = 现a ( 1 r i ) + m a z ( f w ，0 ) n s = d , a ( i p , i ) - t - m a x ( f 。，0 ) b = & a x a y 第三章方稃离散与求解 见= 篙仇= 篙 对上式( 中的a ( i 尸i ) 取不同的值， 种格式： 队= 篙耻篙 可以得到不同的格式，表( 3 1 ) 列出最常用的五 表3 1 各种不同差分方案的函数a ( 1 尸i ) 格式 a ( i p i ) 中心差分 1 一o ，5 i p i 迎风格式l 混合格式 m a x ( o ，1 一o 5 1 p i ) 指数格式 i p l ( e x p ( i p i ) 一1 ) 乘方格式 m a x ( o ，( 1 一o 5 i p i ) 5 ) 3 3 基于s i m p l e 算法的流场计算 对控制方程进行数值模拟的关键是对两个一阶导数项的处理，上一章