（环境工程专业论文）基于POM模型的三维潮流及物质输运数值模拟研究.pdf

摘要 本文在对近年来国内外三维水流数值模拟的有关文献进行了较系统地分析研究的基 础上，以p o m 模型为基础，对三维水流及物质输运的数值模拟进行了一定的研究。论文 首先评述了天然水流三维数学模型的国内外研究现状；其次，采用定边界点的方法，通过 求解拉普拉斯方程组建立了定边界点下正交曲线网格的生成方法，并分别应用于边界规则 的非对称渠槽截面及多连通区域的正交曲线网格的生成；然后，从笛卡尔坐标系下的三维 数学模型基本方程组出发，通过坐标变换，导出了水平正交曲线坐标、垂向盯坐标系下的 三维数学模型基本方程组：给出p o m 模型的离散网格及内外模式方程的差分格式，采用 过程分裂法求解三维数学模型，并对紊动粘性系数的确定以及不同数值大小、不同结构型 式下的垂向紊动粘性系数对数值计算结果的影响进行了讨论和分析：基于p o m 模型，对 矩形水槽风生环流、明渠恒定流两个理论算例进行了计算，得到了与解析解及经验公式较 一致的结果，并将模型应用于具有实际地形的潮汐水域长江太仓段的水流及物质输运的数 值模拟，取得了较满意的结果。该研究成果对潮汐水域的水污染控制具有一定的理论价值 和实际意义。 关键词： 潮汐水域p o m 模型三维水流物质输运数值模拟 a b s t r a c t b a s e do nt h ed e t a i l e da n da b r o a di n v e s t i g a t i o no ft h ei n l a n da n do v e r s e a sl i t e r a t u r e so f t h r e e d i m e n s i o n a lt u r b u l e n c ec o m p u t a t i o ni nr e c e n ty e a r s ，t h en u m e r i c a ls i m u l a t i o no f3 dt i d a l f l o wa n dm a s st r a n s p o r t a t i o nb a s e do np o mi sb u i l t f i r s t ，t h em e t h o d so ft i d a lc u r r e n tn u m e r i c a l m o d e l i n ga r ec o m m e n t e d ，a n dt h ea n a l y s i sf o rt h ed e v e l o p m e n ts i t u a t i o no ft i d a lc u r r e n t n u m e r i c a lm o d e l i n gi sm a d e s e c o n d ，t h en u m e r i c a lp r o c e d u r ef o rt h eg e n e r a t i o no fo r t h o g o n a l b o d y f i t t e dc o o r d i n a t es y s t e m sw i t hd i r e c td e t e r m i n a t i o no fg r i dp o i n t so nt h eb o u n d a r yi s p r e s e n t e d m o r e o v e r t h e t u r b u l e n c ec o n t r o le q u n i o n si nb o d y - f i t t e dc o o r d i n a t ea n d 口 c o o r d i n a t ew e r ec o n s t r u c t e df r o me q u a t i o n si nc a r t e s i a nc o o r d i n a t e ，a n dt h eb o u n d a r yc o n d i t i o n s a n di n i t i a lc o n d i t i o n sa r ea n a l y z e d t h e nt h ed i s c r e t i z a t i o n 鲥da n dd i f f e r e n c es c h e m e sa r e p r e s e n t e d t h em o d e s p l i t t i n gs o l u t i o nt e c h n i q u ei sa p p l i e di n t h i sm o d e l e d d yv i s c o s i t y c o e f f i c i e n t sa r ea n a l y z e da n dt h ec a l c u l a t i o no f t w ot h e o r e t i c a le x a m p l e sa r et a k e nf o ra n a l y s i so f t h ei n f l u e n c eo ft h ev e r t i c a le d d yv i s c o s i t yc o e f f i c i e n t ( k w ) o n3 - dn u m e r i c a ls i m u l a t i o n b a s e d o np o mm o d e lt h ec o m p u t a t i o n so fw i n di n d u c e dc u r r e n ti nac l o s e dr e c t a n g u l a rb a s i n ，c u r r e n t i no p e nc h a n n e la n dt i d a lc u r r e n to ft a i c a n gp a r to fy a n g z ir i v e ra r et a k e n ，t h ec o m p u t e dr e s u l t s a r es a t i s f i e d t h es t u d yo f t h i sp a p e ri ss i g n i f i c a n tt os o l v et h ep o l l u t i o ni nt i d a lr i v e r k e y w o r d s ：t i d a lr i v e lp o m ，t h r e e d i m e n s i o n a l ，f l o w , m a s st r a n s p o r t a t i o n ，n u m e r i c a l s i m u l a t i o n 学位论文独创性声明： 本人所呈交的学位论文是我个人在导师指导下进行的研究工作及取得的 研究成果。尽我所知，除了文中特别加以标注和致谢的地方外，论文中不包 含其他人已经发表或撰写过的研究成果。与我一同工作的同事对本研究所做 的任何贡献均已在论文中作了明确的说明并表示了谢意。如不实，本人负全 部责任。 论文作者( 签名) ： 学位论文使用授权说明 年月日 河海大学、中国科学技术信息研究所、国家图书馆、中国学术期刊( 光 盘版) 电子杂志社有权保留本人所送交学位论文的复印件或电子文档，可以 采用影印、缩印或其他复制手段保存论文。本人电子文档的内容和纸质论文 的内容相一致。除在保密期内的保密论文外，允许论文被查阅和借阅。论文 全部或部分内容的公布( 包括刊登】授权河海大学研究生院办理。 论文作者( 签名) ： 年月日 基于p o m 模型的三维涮流及物质输运数值模拟研究 1 绪论 1 1 引言 目前，随着我国人口增长、工业迅速发展和城市化进程，各种工业废水及生活污水大 量直接排入河流、湖泊、水库和海湾等水体，超过了天然水体的自净能力，使水体环境质 量日益恶化，尤其对于大江河口及湖泊、近海地区，水环境恶化程度不断加剧，生态系统 遭到严重破坏，这不仅影响着居民的生活质量，而且制约着区域的经济发展：水污染加剧、 水质下降和水生态恶化已经成为当前制约社会经济持续发展的重要因素之一1 , 2 1 。水环境 治理己成为亟待解决的问题，并引起人们越来越多的关注。我国对水资源的开发及环境保 护问题已有了很高的重视，为了准确预报河流、湖泊和近海等天然水体的水环境质量，合 理地进行水环境容量及污染物允许排放量计算；科学制定污染物总量控制规划，这就需要 深入进行天然水体水动力状况的研究以及污染物在水环境中的输运迁移规律的研究，而数 值模拟作为解决上述问题的一个重要手段，已经得到越来越多的应用。 数值模拟至今虽然只有几十年的历史，但随着计算机和计算方法的进步，其得到了高 速的发展与应用，已经成为环境、水利、海洋工程问题的重要研究手段1 。数值模拟实质 上就是采用数值近似的计算方法来模拟研究对象的自然现象和动力过程。水流的数值模拟 通过对描述水流运动的控制方程用数值离散求近似解的手段来模拟水流运动，而对大体积 天然水体水流现象的数值模拟，历来是大家十分感兴趣的课题。 由于计算机发展及数值求解技术的限制，以往对大体积水体的研究重点主要是在流动 的二维模拟上。然而，天然水体中的绝大多数流动具有三维特性，尤其是在边界变化较为 剧烈的区域，水体流动的三维性尤为突出，二维数学模型就难以保证相应的准确度。近年 来，工程实践对数值模拟提出了越来越高的要求，进行水环境的三维模拟势在必行，而计 算机及计算技术的发展也为进行三维数值模拟提供了可能性。在环境及水利工程中，为合 理解决污染物输运扩散特性、河口盐水入侵及河床演变等问题，都有必要预测水流的三维 结构。 9 7 3 年，l e e n d e r t s e 等人首先进行三维水流数值模拟的研究工作“1 ，提出数值模拟三 维水流的分层方法，此后也出现了一些简化的三维模型。但是，由于问题本身的复杂性和 计算机运算能力的限制，直到上世纪7 0 年代末期，三维模型才真正得以发展，近些年发 河海大学硕士学位论文 展尤其迅速。迄今为止，三维问题的数值模拟已经成功的用于工程实际计算，在环境保护、 水利工程、海洋工程等诸多领域中得到了广泛的应用。 1 2 三维模型基本方程及常用假定 1 2 1 常用假定 对天然水体直接求解三维n - - s 方程目前尚不可能，为此，人们提出了一些近似假定， 使实际天然水体三维问题具有求解的可能性，常用的假定有： 1 ) b o s s i n e s q 近似由于盐度、温度和物质浓度变化较小，除在动量方程的重力项中 必须考虑密度变化的影响外，其余各项一般可忽略不计，即在动量方程中，除重力项外， 视密度为常数。 2 ) b o s s i n e s q 假定即b o s s i n e s q 在1 8 7 7 年提出的关于水流紊动应力可类比于层流粘 性应力的假定。 3 ) 静水压强假定 天然河流、浅海及湖泊水域具有宽浅特征，垂向加速度远远小于 重力加速度，因而在垂向动量方程中忽略垂向加速度而近似假定为静水压强关系。目前水 流运动的三维数值模拟大多都采用该假定。 在上述假定的基础上，建立并发展了各种不同的三维水流数学模型，并因各种数值方 法的应用和改进，使三维水流计算发展成各种不同类型的数学模型。 1 2 2 基本方程 通过以上假定，可以写出描述天然水域三维水流运动的基本控制方程组( 常见或常见 书写形式) ，可见式( 1 - 1 ) ，而所谓的数值模拟的方法就是运用适当的初始条件和边界条件， 采用恰当的数值离散格式和求解方法，求解模型的控制方程组，得到一定程度上满足工程 实际要求的数值解。 a ua 矿a w _ + _ + _ = 0 o x e y 0 z o u + u 型+ y 型+ 丝一：一上望+ 监+ 堡+ 篮 8 te x a yezpa x a x a y 0 ： ( 、1 ) 里+ u 里+ ，里+ 里+ 彤：一上至+ 堡+ 堡+ 堡 a t a x a ye：pa y o x a y o z 上娑= 一g 口o z 2 一 堆于p o m 模型的三维潮流及物质输运数值模拟研究 其中k ：2 厶掣 d x f 。：2 a m _ = o v c r 铲厶f 警+ 豢1 出砂 f 。吨o uf 户k 。鼍 式中：t y ，z 构n - 或x o y 面置于未扰动静止水面、z 轴铅直向上的笛卡尔坐标系：u ，v ，矿分别 为流速沿，y ，= 轴方向的分量；p 为水密度：p 为水压强；a 。为水平紊动粘性系数；k 。为 垂向紊动粘性系数。 1 3 三维数学模型研究概况 国内外学者近年来对三维水动力学数学模型的研究赋予了极大的热情，目前，三维数 学模型的研究已经得到了较大的发展，本文对三维数学模型的研究进展综述如下： 1 3 1 分层二维法 这种方法就是将水体沿水深方向( 垂向) 分成若干层，将三维方程化为各水层的二维 方程来处理。l e e n d e r t s e 等。1 于1 9 7 3 年最早推出了这种方法，他假定各层中的变量在垂向 上变化不大( 可以忽略) ，可简单的用积分平均来表示变量的一阶近似，从而将三维问题 转化为多个二维问题，并在垂直方向采用固定分层法，即将计算水域划分为固定的多层， 在每层中沿水深积分使之成为二维问题，并用a d i 格式进行数值离散；赵士清( 5 1 于克俊 等：1 也应用分层二维法分别对长江口南槽和口外海域的三维潮流场和渤海潮波运动进行了 数值模拟：卢启苗。1 采用分层二维的方法及平面集中质量有限元建立了海岸河口的三维潮 流数学模型；白玉川“1 在变换的垂向坐标下，采用分层二维法和算予分裂法，将三维控制 方程分解为易于求解的物理意义单一的方程，用特征差分近似和g a l e r k i n 集中质量有限元 相结合的方法对廉州湾潮流进行了三维模拟计算：张涤明等。”则将分层二维中的每层中的 变量用三次样条插值表示，在每层内对三维控制方程积分便得到一组二维方程，与上面的 积分平均分层不同，样条分层的每层水深是随时间变化的。 分层二维法简单实用，但该法在水深较大的水域计算比较理想，在浅水部分的垂向分 辨率很低，并且引入了人为交界面，交界面上必须给定摩阻系数，摩阻系数的选用往往带 有很大的经验性。 河海大学硕士学位论文 1 3 2 整体三维模型 1 3 2 1 原始坐标系( x ，y ，z ) 法 该方法就是在未经变换的原始坐标系( 置y ，? ) 中求解三维控制方程( i - 1 ) 而获得三 维流场。逢勇等“”采用( 丘j ，z ) 坐标系下的三维水动力模型，考虑太湖局地风场的影响， 模拟了太湖风生流、风涌增水的变化情况：李狮等“1 采用混合破开算子法，并引入透水 介质法处理动边界、自由水面边界和床底边界，对大尺度水体流动进行了成功的模拟：王 晓建等“由三维浅水方程出发，利用有限体积法的积分模式并引入g r e e n - g u a s s 公式， 将控制体单元的体积分转化为面积分，采用一种任意六面体单元划分计算域，较好的适应 了复杂边界，对辽东湾海区流场进行了数值模拟，并获得理想结果。 但该类模型应用有时不是十分方便，例如在水深变化较大的浅水区，其层数可能较少， 层面不同，模型的侧边界条件位置不同，难以较好的反映潮流垂直结构，使程序的适用性 受到一定限制；当计算区域边界比较复杂时，采用矩形计算网格，则会使模型无法较好的 模拟边界处水流，出现虚假的曲折水流。为了解决以上问题坐标变换的方法引入了三维 水动力模型的研究中。 1 3 2 2 ( 工，y ，盯) 坐标系法 该法“”是在z 向引入口= ( z f ) ( h + ) 或口= ( 一z ) ( h + f ) ( 或f 的其他形式) 变 换，将整个计算域的水体厚度处处变为单位l ，同时将( ty ，：，f ) 坐标系中的三维控制方 程变成( x ，y ，口，f ) 坐标系中的方程：在口坐标系中求解方程并将计算结果转换到( x ，y ，= ，f ) 坐标系中从而完成三维控制方程的求解，a 坐标变换法将自由水面和不规则水底变成了a 坐标系中的表层和底层坐标平面，“水深”为1 ，这不仅使得整个计算水域垂向具有相同的 网格数且可以随意分层，从而保证浅水部分有了更高的垂向分辨率，而且从数值方法上讲， a 坐标系中方程的离散求解要容易得多，因此，a 坐标变换法在三维水流数值计算中获得 了广泛的应用。d a r a n c e 和h u g h e s “4 1 在1 9 8 3 年将仃坐标变换方案应用于欧洲北海；窦振 兴“”将仃坐标变换引入三维潮流计算中，并基于过程分裂法( 也称模态分裂法) “卜1 ， 即将三维流动中的快过程一表面重力长波与慢过程一缓行的内重力波分开，并对每部分采 用适合其物理特性和数值行为的计算格式耦联求解得到了三维流场：朱耀华“1 将a 坐标 綦十p o m 模型的= 维潮漉及物质输运数值模拟研究 变换应用于我国渤海、黄海、东海的海洋环流的研究：韩国其“”1 在该坐标变换下，采用 算子分裂法对长江口南支进行了成功的模拟，并使用k s 双方程紊流闭合模型确定紊动粘 性系数；胡维平”1 应用a 坐标变换，采用半隐差分格式对太湖风生流进行了模拟，反映 了太湖水位及湖流随风场的动态变化，揭示了湖流的垂直结构；李孟国”则在水平方向 使用正方形网格和不规则三角形网格的差分离散格式，垂向采用隐式中心差求解海岸工程 海区的三维流场：石磊等。在三维浅海流体动力学模型中也引入盯坐标，在时间上采用 分步法，同时模型采用三角形网格，模型有较好的灵活性和通用性：董文军”“采用平面 逆风有限元法与垂向隐式差分相结合的方法，并引入过程分裂法，对天津新港附近渤海海 域的三维潮流运动进行了数值模拟。 1 3 2 3 ( 毒，仉：) 坐标系法 为了克服了矩形网格差分法不能精确拟合天然水域不规则水陆边界的困难，研究人员 在水平方向使用边界拟合的正交曲线坐标，垂向未进行坐标变换，在( 善，q ，：) 坐标系求 解三维模型。华祖林”“采用该坐标变换模式，对长江三峡工程弯曲河段进行三维数值模 拟，在模式引入紊流模型以反映局部小范围水流运动状况，并应用动水压力校正法1 2 5 1 将 三维潮流控制方程中的压力p 分解为动水压力与静水压力之和，由控制方程导出三维控制 方程的离散格式，采用p a t a n k e r 和s p a l d i n g 提出的压力校f 法求解动水压力场，然后求解 垂向积分型连续方程获得水位，用最新求出的水位替代压力校正法中的水位，如此循环迭 代直至收敛。 1 3 2 4 ( 毒，叩，盯) 坐标系法 该法在垂向上使用a 坐标变换，具有a 坐标变换的优点，水平方向使用边界拟合的正 交曲线坐标，可以精确拟合天然水域的复杂边界。该方法以p r i n c e t o n 大学的海洋动力环 境数值模拟小组的p o m ( p r i n c e t o no c e a nm o d e l ) 模型”“。1 为代表，该模型在垂向采用口 坐标系，在水平方向采用边界拟合的正交曲线坐标系，采用m e i i o r - - y a m a d a 的双方程紊 流闭合模型确定垂向紊动粘性系数，并采用过程分裂法求解，该模型目前在海洋河口潮流 的数值模拟中有着较为成功的应用：t a lg z e r 。”等首次将p o m 模型应用于大西洋的整体 流场及温j ；：场的数值模拟：万振文”1 等采用p o m 模型较成功的模拟了渤、黄、东海的潮 波运动：任华堂等应用p o m 模型模拟了季风在台湾海峡引起的风海流以及三维流场结 河海大学硕l 学位论文 构，在冬季和夏季两种条件下，得出了季风引起的水位变化和三维流场结构：管卫兵”等 在p o m 模型基础上添加示踪粒子三维l a g r a n g e 运动轨迹计算及水质模块，研究了榆林湾 水交换能力及污染物输运过程；刘子龙等“在( # ，_ ，叮) 坐标系下，采用算予分裂，对长 江e l 南支部分水域的三维流动进行了较成功的模拟：史峰岩等。“1 在该坐标系下，对e c o m 模型进行改进，建立了一个以杭州湾和长江口为整体的三维联合模型，用于潮流和物质输 运的计算研究；华祖林“等采用该坐标法除对长江口南支河段的水流模拟外，还对物质 输运进行了模拟。本文将以p o m 模型为基础，在平面采用正交曲线变换，垂向采用仃坐 标变换，开展潮汐水域水流及物质输运数值模拟的研究。 1 3 2 5 ( x ，y ，z ) 坐标伸缩变换法 该法”1 除了在z 向使用a 坐标变换外，在x 、y 两个方向上办使用伸缩坐标变换，使 整个计算区域在x 和y 向的尺度变为单位1 ，从而使整个计算域变成单位正方体。该法具 备a 坐标变换的优点，在水平方向上能拟合水陆边界，但变换后的方程相对要复杂些。匡 翠萍。将整个计算水域通过简单的坐标变换转换成单位正方体，三个坐标方向的空间网 格均可采用非均匀网格，对变换后的基本方程进行差分离散，水平采用中心差分，垂向采 用隐式中心差分，建立了河口区赫水入侵汁算的三维数学模型，并将该模型用于长江口水 域的盐水入侵计算。 上述几种坐标变化法中，x ，y ，z 均采用伸缩坐标变换法应用并不多见；而垂向盯坐标 变换较简单，应用晟为广泛，其中以( x ，y ，盯) 坐标系下的数学模型应用较多：但由于在 平面采用正交曲线坐标变换能够精确拟合计算边界，所以( 善，r l ，= ) 和( 孝，叩，盯) 坐标系下 的模型也得到了较多的应用，其中( 善，叩，盯) 坐标系下的三维数学模型又以其在水平方向 和垂向的良好的地形拟合性得到了越来越多的应用。 1 _ 3 3 解析方法 该法1 1 的基本思想是先用水深积分的二维模型计算出整个计算域的垂向平均流速， 将u 和v 的垂向分布用垂向平均流速和表示流速垂向变化的函数之积表示，即用一维垂向 模型来计算各节点的速度剖面，这种方法计算简单，速度快，但计算结果的合理性与表示 流速垂向变化的函数密切相关。 基于p o m 模型的三维潮流发物质输运数值模拟研究 1 3 4 谱方法 该法最初是由h e a p s 。渭提出的，他将两个水平流速分量u 和v 在垂直方向的变化以本 征函数为基函数作为一维谱展开式( 级数形式) ，以本征函数为权函数使运动方程自水底 至水面的权余量为0 ，获得谱展_ 丌式中系数的二维方程，与垂直积分型连续方程一起构成 一组联立方程，从而将三维问题转化为一组二维问题来求解。d a v i e s ”1 在谱展开式中不用 本征函数，而是采用各种广义傅氏级数甚至样条函数作基函数( 这样就不必解出本征问题， 甚至也不必顾及基函数是否具有完备正交性) 。谱方法能将三维问题转化为二维问题来处 理，并能获得u 和v 沿垂向的连续分布，其缺点是在谱展开式中要展开很多项才能达到精 度要求。 1 3 5 流速分解法 该法“”1 的基本思想是引入复型流速向量q = ( u + i v l ，将三维运动方程化为q 的定解问 题，并进一步处理后得到q 的流速分解定解问题，将q 直接分解为几项具有一定物理意义 的流动速度，每一项中均含有一个垂直剖面函数，根据剖面函数的定解问题可以获得其解； 将q 沿水深积分得到垂向平均复型流速向量，将有关各量代入垂直积分型连续方程则得到 单纯水位的二维问题和求解有限个垂宜剖面函数得一维问题，流速分解法属于“二维半” 方法，它亦能获得u 和v 沿垂向的连续分布，但该法必须将垂向紊动粘性系数写成线性的 变量分离形式，并假定紊动粘性系数的垂向剖面形式。 1 4 本文主要工作 本文在分析前人的三维水流数值模拟研究的基础上，基于p o m 模型，在垂向采用 o 一- - - ( z 一# ) ( h 十；j 的坐标变换法，使其更好的拟合实际地形：在水平方向采用边界拟合的 正交曲线坐标较好的拟合实际的计算域边界，以提高计算的精度。由此将模型应用于潮汐 河段的水流、物质输运的三维数值模拟中。 本文第一章对三维水流数学模型的进展作出了分析和总结，并简要说明了本文的主要 研究内容。第二章对贴体的正交曲线网格的生成进行了讨论，在给定边界网格点的情况下， 通过求解拉普拉斯方程建立了定边界点下萨交曲线网格的生成方法。论文的第三章通过数 学推导，建立了新坐标系下的三维数学模型，导出了新坐标系中的紊流控制方程组，给出 了其推导过程及相应的边界条件及初始条件。第四章采用有限差分法对新坐标系中模型控 河海大学顺l 学位论文 制方程进行了数值离散，给出了内外模式方程组的离散过程及差分格式，并对稳定性条件 进行了简要分析，同时对不同数值、不同结构型式的垂向紊动粘性系数对数值模拟的影响 进行了较详细的讨论。第五章分别对假想的矩形水槽风生环流、理想水槽恒定流和潮汐河 段长江太仓段的潮流及物质输运进行了计算。第六章对本文的工作进行了简要的总结，并 指出了今后有待进一步深入工作。 基于p o m 模型的三维潮流及物质输运数值模拟研究 2 定边界网格点下正交曲线网格的生成 在实际天然河流、湖泊、海湾地区水流输运问题的求解中，人们最早以矩形网格来划 分计算区域，在直角坐标系中由于矩形网格节点间关系容易确定，也有利于用差分逼近导 数，因而也便于组织数据结构和各种算法的实现。但是天然水域岸边界通常都为不规则， 使用矩形网格会将岸边界概化成锯齿形边界，陆地附近出现虚假的曲折水流，有可能出现 计算奇点，边界附近数值解误差较大：若要提高边界计算精度，只能通过加密网格来实现， 这又必然增加网格数量，导致计算量的增大。针对这种情况，边界拟合坐标系“1 ( 或称 贴体坐标系b o u n d a r y 占d d y f i t t e dc u r v i l i n e a rc o o r d i n a t es y s t e m s ) 技术应运而生，所 谓的边界拟合坐标系，实际上是一类数值网格生成方法，这类方法能够较容易、精确的描 述复杂几何形状物体的边界，其大致可分有结构与无结构网格两类。有结构网格的特点是 每一节点与其邻点之间的联结关系固定不变且隐含在所生成的网格中，因而不需要专门设 置数据去确认节点与邻点之间的这种联系，便于以后进行水流及物质输运的计算。对于无 结构网格，单元与节点的编号无固定规则可遵循，因而除了每一单元及其节点的几何信息 必须存储外，与该单元相邻的那些单元的编号等也必须作为联接关系的信息存储起来，使 无结构化的存储信息量较大。正交曲线网格属于有结构网格的一种，除了具有一般有结构 网格存储量小的优点外，还有下列几个特点：可以相对简化曲线坐标下的控制方程，还 可以避免由于网格的非正交性而引起的某些离散误差：虽然实际流体的流动方向并不与 网格线完全一致，但两者间的偏差可望小于按一般方法所生成的非正交网格中的偏差，有 利于减少假扩散误差；有利于数值解的收敛等。因而，近几十年来，正交曲线网格以其 良好的贴体性及较高的计算精度而在水流数值模拟中得到了广泛的应用。 2 ，1 曲线网格生成方法概述 曲线网格生成方法主要有三种：保角变换、代数生成法、微分方程生成法三大类。其 中保角变换法的主要优点是，可以保持物理平面上的生成的网格的正交性，但由于所应用 的变换是复变函数理论中的一些成熟的关系式，所以仅适用于某些特定的问题；代数法生 成网格是利用各种插值公式来建立起计算平面与物理平面之间节点的对应关系，其优点是 方法简便、节省机时，同时也能在一定范围内控制网格的形状与密度，但与椭圆型偏微分 方程方法相比，虽然线性的无限插值方法简单；但是边界上的斜率的不连续性会传递到计 河海大学顾士学位论文 算区内部去，而在采用h e r m i t e 插值来控制边界上的f 交性时又会引起计算区域内部网格 的重叠及不光顺性。 微分方程生成法是将贴体坐标的网格生成看成一个边值问题，通过求解偏微分方程生 成网格，利用偏微分方程的一些性质使生成的网格更加完善、合理。该方法最早由w i n s l o w 在1 9 6 7 年提出1 4 2 1 而后由t h o m p s o n 等人于1 9 7 4 年比较全面而系统的提出 4 3 1 。 t h o m p s o n 等通过求解泊松方程组来生成曲线坐标系。其方程如下： 姿+粪：p(翻ox ! o y 2 、7 。 宴+ 窖：q ( 翱o瓠： y ! 、 ( 2 1 ) 其中：p 、q 是具有“合并”或“聚集”控制函数( p a c k i n gf u n c t i o n ) ，t h o m p s o n 建议p 、 q 函数的表达式为： 尸( 善，1 ) = 一qs g n ( f 一眚) e x p ( 一t l 一善i ) 毒，s 驴( 5 - 乞) e x p 卜肛每) ：+ ( 嗍h q ( # ，叩) = 一a is g n ( 叩一研) e x p ( 一q h - r l ) 一喜qs s n ( - 一仉) e x p 一d ， ( 一) 2 + ( ”一”，) 2 只 ( 2 2 ) p i l t h o m p s o n 提供的曲线网格生成方法并不能确保网格的正交性，而p 、q 函数中可调参数太 多，当网格生成不理想时，很难确定该调节哪一个参数才能改善网格；并且网格疏密程度的 调节是以一定数量的等善，叩线为中心，对大面积水域，局部加密不够灵活；虽然可以通过构 造不同的p 、q 函数“；。来调节网格线的位置变化，并且积累选择p 、q 函数中调节因子 数值的一定经验后，可获得较为满意的结果，但也有局限性：物理坐标系中网格点的位置不 能预先给定，只能通过调节p 、q 函数中的系数来调节网格线的疏密；如果需要在计算域某 边界指定位置布置网格，或布置按某种规律变换的网格，用该方法就难以实现。例如，在水 流数值模拟时，计算者往往希望在水工建筑物的边界变化处，或河流口门处布置特定的网格， 使用t h o m p s o n 方法就会遇到一定的困难。 在实际工程中，有时需要首先确定边界处的网格，本文在n a a t h a n a s s i a d i s ，余利仁“6 4 1 等的基础上利用物理坐标系在变换平面中求解拉普拉斯( l a p l a c e ) 方程组( 2 - 2 ) 式进行 基于p o m 模型的三维潮流；5 乏物质输运数值模拟研究 了定边界点下正交曲线网格的生成。 些+ 氅：o 缸。a v ： a 二口a qn a x z + 萨刈 ( 2 2 ) 该法优点在于可在两邻边上根据实际需要选择、布置网格点，对水利工程、环境工程的 数值计算有突出作用。因为在这一类问题中，往往一侧边界较为复杂，有入流或建筑物，计 算时需要布置特定的网格点。本文应用该方法，分别进行了非对称渠槽截面及多连域的网格 生成，得到了较满意的效果。 2 2 定边界网格点下正交曲线网格的生成 2 2 1 拉普拉斯( l a p l a c e ) 方程性质 在如图2 1 a 所示的物理区域内，考虑调和函数f = ( x ，y ) ，即v 2 善= 0 ，并设p 为域中任 意一点。 为 bc d ad 图2 - - 1 a ，b 物理平面及映射平面示意图 现将 ( x ，j ，) 在p 点邻域内按t a y l o r 级数展开： 4 ( x , y ) = 鼻+ ( 纠，x + ( 专，卜j ，+ 圭 ( 。) ，x 2 + 2 ( 岛) ，砂十( 厶) ，y ! + 特别的在p 点附近邻域内取作2 h 2 h 的小正方形领域( 如图2 - l b ) ，那么该邻域的均值 改写为 j = 去lh lh d 】c 毋气p + 专蝗。+ ，、p 河海大学硕l 擘位论文 ( 。+ 专，) ，= 素( # 一缶) ( 2 3 ) 由此，式( 2 3 ) 体现了l a p l a c e 算子的物理含义，其表征所在点微小邻域内函数均值与该点 函数值之间的差值。由式( 2 3 ) 可以推演满足l a p l a c e 方程的调和函数 值存在极大一一极 小值法则，即“如果函数 在区域q 内满足v2 掌= 0 ，则f 的极大( 小) 值不可能发生在域内 而只可能发生在边界上”。无疑，满足拉氏方程的调和函数f 的变化梯度是单向性的( 递增 或递减) ，是有顺序的。这样在整个定义域内就不同点( t y ) 不可能出现相同的f 值，这保 证了以拉氏方程作为坐标变换方程时，绝对不会出现多值对应的情形。 2 2 2 网格生成原理 利用l a p l a c e 方程的性质，本文生成f 交曲线网格。设( x ，y ) 为物理平面上的二维笛卡 尔坐标系，( 害，叩) 为变换( 或称映射) 平面上的直角坐标系，通过数值变换，希望将物理平 面上的正交曲线网格变换为映射平面上单位正方形域中的规则正交网格。根据共轭调和函数 的性质，当孝，叩是一对共轭函数时，孝和叩在( x ，y ) 平面上的等值线簇就是一对正交曲线， 在映射平面( 孝，玎) 上的相应值构成所需要的规则正交网格。 下面通过共轭函数满足柯西一黎曼条件和映射与逆映射的存在关系，推导出映射平面规 则正交网格中z ，y 所应满足的方程，首先，设一对共轭函数： f = f ( x ，y ) ( 2 - 4 一a ) 【_ = 口( 丘y ) ( 2 - 4 - b ) 这时，孝和7 7 在( j ，y ) 平面上的等值线簇为： 4(x,y)：=co(2-5r(x ) l ，) = e 而这就是一对彼此正交的曲线，它们在映射平面( 善，叩) 上的相应值构成了我们所需要的规则 正交网格。这一对共轭函数应满足拉普拉斯方程组： 当给定了计算域四周的边界条件之后， 这对共轭函数办满足柯西一黎曼条件： 爹专 任。， 窘= 。 一 基于p o m 模型的三维潮流及物质输运数值模拟研究 根据逆映射的存在，设 f z = z ( 善，刁) l y = y ( ，1 ) 对( 2 - 7 一a ) 式两边求孝的导数，( 2 - 8 口) 式两边求x 的导数，得 将( 2 - 9 ) 中两式相减得 又根据柯西黎曼条件( 2 - 7 一b ) 1 ：丝鱼+ 巩葩 1 ：鱼堕+ 葩教 硝砂 砂鸳 瓠雒 a f 出 篓皇：里塑 砂骘a 出 箜立：鱼f 一堕1 毋鹫t 却l 砂 由于映射( 孝，叼) 的逆映射的存在，显然有些0 ，故： 咖 同理可求出 巩卸 a r 西 ( 2 7 一口) f 2 7 - b ) ( 2 8 一口) ( 2 8 6 ) ( 2 9 ) ( 2 1 0 ) ( 2 1 1 ) ( 2 一1 2 ) 兰：兰 ( 2 1 3 ) a fd 、 即映射g ，y ) 满足柯西一黎曼条件，因此，它们在皓，叩) 平面上也是一对共轭调和函数，满足 拉普拉斯方程。 为使毒，叩无量纲化，使其在o 1 中取值，在上式中引入一尺度因子h ，方程变为： 根据a n t o n o p o u l o s “”的计算，驳 = ，这早 是善等值线的总长，是叩等值线的总长。 “ 却一砂忉一出 | i 一 鸳一缸西一妙 舢 却 枷 焉啬 投一妒西一 1 l 海人学顺l 学位论文 2 2 3 定边界点方法 如图2 一l 所示，采用上述方法，可将物理平面的区域a b c d 转换成映射平面的单位面 积域a b c d ，由于物理平面区域的边界线也是网格线，因此边界线需满足正交条件，即a 、b 、 c 和d 两侧局部边界应相互垂直。水利工程和环境工程中的绝大部分流体数值计算问题都 可以满足这一条件。我们将一对对边( 如a b ，c d 边) 取作进口和出口断面。 所谓定边界点，即根据实际工程需要，先在物理平面的一对邻边上布置网格点，这些 点在计算过程直至f 交曲线网格生成完成后位置都不发生变化，通过另一对邻边上网格点的 滑动使得网格线f 交。确定好一对邻边的网格点后，再从这些点出发分别向对边作网格线， 得出初始网格的坐标值，同时给出初始的亭，7 7 值。由于一对邻边的网格点固定，而初始给 定的孝，叩值并不一定跟网格点真实的掌，叩，照相等，因此，我们的目的是求出各网格点上 对应 ，玎的真实值，再通过求解拉普拉斯方程组( 2 - 1 4 ) 得出内部网格点。 2 2 4 数值求解过程 2 2 4 1 方程离散 在以尸( 善，叩) 点为中心的控制体积中，将( 2 - 1 4 ) 式离散为： a p x p2 以+ 句+ 爿s x , 7 + a n x , v( 2 1 5 ) 【爿，y p = a e y + 山y w + a s y s + a y “ 式中系数分别为 以2 瓦乏素习 以2 瓦乏钢 a , v2 瓦i 顶石i 耵 l _ 一儿叩一叩pj n 一 4 2 瓦i 丽丽r sl 叩一，7 s 儿玎p 一，一 其中下标e 、矽、v 、s 分别代表p 点的四个相邻点。 ( 2 1 6 ) 基十p o m 模型的二三维潮流发物质输运数值模拟研究 2 2 4 2 求解过程 真实亭，叩值的求解是一个逐步校f 逼近的过程，即：在初始网格中，网格与边界一般 不正交，假定欲求解的掌为常数的线是从南部边界a i d 上任选点出发，如图2 - 2 所示： 舌 曼 岛 4 - l ，2 图2 - 2 角度修正不意图 其中实线为初始网格线：虚线是j = l 线上通过边界点的垂线，并与j = 2 线相交：空心点为初 始网格点；实心点为新点；垂线与初始网格线有一夹角，善值的校正就是通过此夹角来进行： 鸳亍c 。勋 比例因子c 一般可取0 , 0 0 0 0 5 “，根据初始网格线与边界垂线的关系，以图2 2 为例，修正 后的掌+ 为 舌= 夤一苗，最+ = 参一鸳：，茧= 岛+ 由上，假设固定左边界和下边界网格点，二维正交曲线坐标系的数值生成步骤可归结 为以下几点： 1 给出初始的网格点坐标和f 、7 7 值： 2 从下边界点出发，作边界垂线与j = 2 线相交，得出j = 2 线上的新的网格点坐标，根据 勋求出劈，进而求出新的手值； 3 上边界进行张力样条函数插值，并且过j = n 1 线上的点，对插值后的上边界作垂线， 垂足即为上边界的新点。在此过程中，上边界的点是可滑动的，这一步保证了上边界的正交； 4 同理，从左边界出发，求i t i 新的叩值和声边界的新点： 5 保证了边界的正交和得出新的、7 7 值之后，求解拉普拉斯方程组，得出新的网格坐 标点，判断是否收敛，如果没有，回到步骤2 。 2 3 实例 为说明定边界点二维正交曲线网格生成方法的特点，下面举了两个例子，第一例是边 1 5 河海大学硕l 学位论义 界规则的非对称渠槽截面的网格生成，根据需要在其左边界和下边界定网格点，生成如下网 格，由图2 3 、2 _ 4 可见，该方法可以通过确定边界点的密度而人为控制网格密度，如图2 3 中对a 处网格加密。 aa 图2 - 3 物理区域上的正交曲线网格( x ，y )图2 4 计算区域上的网格( 善，玎) 第二例是分叉区域或多连域的网格生成。利用公共边将此区域分成两个或多个单连域， 并在公共边上布置网格点，然后在每个单连域上分别求解生成网格。由于公共边网格点位置 不变，因此每个单连域正交曲线网格生成后将其合并，公共边的网格点重合，从而构成整体 网格。如图2 5 、2 6 所示，其中a 、b 、c 、d 与a 、b 、c 1 、d 1 分别对应为同一点。 一一 霹凌筌麟龌篱 2 4 小结 图2 - 5 分开生成的网格图2 - 6 合并的完整网格 数值网格生成技术是近年来十分活跃的研究领域，它大大的开拓了有限差分等数值求 解技术的应用范畴，正交曲线网格作为有结构网格的一种，由于表述的基本控制方程简单， 能使用于许多已经发展成熟的标准算法和求解技术，因此在工程上有广阔的应用前景。 本文提出的定边界点下正交曲线网格的数值生成方法用物理坐标在变换平面上求解拉 普拉斯方程生成正交曲线坐标系，物理概念明确，方法简单，易于掌握，适用于在工程上经 常遇到的大量的不规则区域中f 交曲线网格的生成，特别是在区域内某一边网格需要特别给 定的情形。在f 交曲线网格生成过程中，一对邻边边界网格点固定，因此这对边界点与实际 边界点完全重合，其对边采用张力样条函数插值，从而保证了边界的光滑性，并且网格边界 线与实际边界线重合；该方法可随意控制网格的疏密度。利用这一特性可将分汉区域或多连 域分割成多个单连域分别求解，使得对几何图形复杂的计算区域网格生成简单化。 1 6 基于p o m 模型的三维潮流及物质输运数值模拟研究 3 三维数学模型基本方程的导出 三维水流数学模型在维度上与实际水流一致，能尽可能真实的表现水流情况，近年来 在河口、海岸水动力及物质输运研究中得到了广泛的重视和研究。其中m e l l o r 等人的海洋 紊流模型p o m ( p r i n c e t o no c e a nm o d e l ) 在海洋潮流模拟中有着较多的应用。该海洋紊流 模型具有以下特点“”1 ：采用坐标变换技术，较好的拟合实际地形，网格划分在水平方向采 用正交曲线网格，在垂直方向采用t o 坐标网格，平面上使用曲线网格使不规则区域转化为 规则区域，在垂向上使用仃坐标后使整个水域具有相同的垂向分层数，给数值计算带来便 利；模型求解采用过程分裂法，将求解过程分为内模式和外模式，即先求解二维外模式获 得自由表面及垂向平均流速，接着求解三维内模式获得三维流场。目前，该模型及其改进 形式已经在国内外相关研究中得到了一定的应用。“。本文基于该模型建立了适用于潮汐 河段的三维潮流物质输运数学模型，首先导出模型的基本方程组。 3 1 笛卡尔坐标系下三维数学模型基本方程 连续万程： o u o + 盟+ o w o ：o( 3 1 ) o x o砂。瑟。 动量方程： 等+ 饥等+ 券+ 等一形一吉毒+ 鼍+ 玺+ 玺 c s 圳 些+ 警+ k 警+ 比誓+ f u o ：一土要+ 誓+ 堡+ 挈 ( 3 - 3 ) o t o吼出。p 饥吼o y 。 上罢：一g ( 3 4 ) 物质输运方程： 鲁+ 乩鼍+ 喏+ 鼍= 杀 。坠& oj + 去( q 象 + 毒 b