（环境科学专业论文）波浪水槽和平面泊肃叶流中扰动空间演化的数值模拟.pdf

a b s t r a c t f o rt h eg l o b a lb i o g e o c h e m i c a lc y c l e s ，g a se x c h a n g ei so n ei m p o r t a n tl i n k b e t w e e nt h ea t m o s p h e r ea n dt h eo c e a n t h u s ，i nr e c e n t l y ，t h eg a se x c h a n g eb e t w e e n t h ea t m o s p h e r ea n dt h eo c e a nh 酏r e e e t v e dm u c ha t t e n t i o n a i r - w a t e rg a st r a n s f e ri s a l s ov e r yi m p o r t a n ti no t h e ra r e a s ，s u c ha se n v i r o n m e n t a le n g i n e e r i n ga n dc h e m i c a l e n g i n e e r i n g a st h ef i r s ts t e po fs t u d y i n gt h eg a st r a n s f e ro nt h ei n t e r f a c i a ls u r f a c e ，a n u m e r i c a lw a v et a n ki sb u i i ti nt h i sd i s s e r t a t i o n i nt h i sn u m e r i c a lw a v et a n k ，t h et w o p h a s ef l o wt h e o r yi sf i r s t l yi n t r o d u c e dt o m o d e lt h eb e h a v i o ro fw a t e ra n da i r t h ei n t e r f a c ei st r a c k e db yv o fm e t h o d t h e w a v ei si n t r o d u c e db ya no s c i l l a t i n gp l a n ea n db e f o r et h eo u t l e tb o u n d a r y ，aa r t i f i c i a l d a m p i n gl a y e r i s p r e s e n t e d t oa b s o r bt h er e f l e c t i o nw a v e t h e c o n t r o l v o l u m e - b a s e dt e c h n i q u ei su s e dt oc o n v e r tt h eg o v e r n i n ge q u a t i o n st o a l g e b r a i ce q u a t i o n st h a tc a nb es o l v e dn u m e r i c a l l y t h en o n - u n i f o r mr e c t a n g l eg r i d s a r ea r r a n g e da c c o r d i n gt or e q u i r e m e n t s s u c h 船f i n e rg n d sa ta i r - w a t e ri n t e r f a c e ， t h eo u t l e tp r e s s u r ec o n d i t i o ni su s e da tt h et o pb o u n d a r y b yc o m p a r i n gt h et h e o r y w i t hw a v ep r o f i l ea n dt h et i m e - d e p e n d i n gp r e s s u r ea n dv e l o c i t ya tac e r t a i np o i n t ， t h ec o r r e c t n e s so fo u rn u m e r i c a lw a v et a n ki sp r o v e d n u m e r i c a lr e s u l t ss h o wt h e a b s o r bi n t e n s i t yi sm a i n l yd e p e n do nt h ea m p l i t u d ea n dl e n g t ho f w a v ea n dt h ed e p t h o fw a t e r t h ee f f e c t so ft h ea i rp r e s s u r ea n da i ra n dw a t e rv i s c o u so nt h ew a t e rw a v e a r es t u d i e d i ti sf o u n dt h a tt h ea b o v ef a c t o r sh a v el i t t l ei n f l u e n c eo nw a t e rw a v e t h es p a t i a ld e v e l o p m e n to fd i s t t t r b a n c e si np l a n ep o i s e u i l l ef l o wi ss i m u l a t e d n u m e r i c a l l y t h ec o n t r o l l i n ge q u a t i o n sa r en - se q u a t i o n s t h em e a nv e l o c i t yp r o f i l e i sa tt h ei n l e tb o u n d a r ya n dt h eo u t f l o wi su s e da tt h eo u t l e tb o u n d a r y a t o l l m i e n - s c h l i c h t i n g ( t s ) w a v ew i t hs m a l la m p l i t u d ek n o w nf r o mo r r - s o m m e r f e l d e q u a t i o ni ss u p e r i m p o s e do nt h em o a nf l o wf i e l dt oi n t r o d u c eas m a l ld i s t u r b a n c e 。 n u m e r i c a lr e s u l t ss h o wt h a tt h es p a t i a ld e v e l o p m e n to fd i s t u r b a n c ei si nag o o d a g r e e m e n tw i t ht h a to fo r r s o m m e r f e l de q u a t i o n k e yw o r d s ：a i r - w a t e ri n t e r f a c e ，c f ds i m u l a t i o n ，n u m e r i c a l w a t e rt a n k ， p o i s e u i l l ef l o w , s p a t i a ld e v e l o p m e n to fd i s t u r b a n c e ，f l o ws t a b i l i t y 独创性声明 本人声明所呈交的学位论文是本人在导师指导下进行的研究工作和取得 的研究成果，除了文中特别加以标注和致谢之处外，论文中不包含其他人已经 发表或撰写过的研究成果，也不包含为获得墨生盘鲎或其他教育机构的学 位或证书而使用过的材料。与我一同工作的同志对本研究所做的任何贡献均已 在论文中作了明确的说明并表示了谢意。 学位论文作者签名：飞务母日 签字日期： 。厂年胡、w 学位论文版权使用授权书 本学位论文作者完全了解盘鲞盘茔有关保留、使用学位论文的规定。 特授权苤鲞盘茎可以将学位论文的全部或部分内容编入有关数据库进行检 索，并采用影印、缩印或扫描等复制手段保存、汇编以供查阅和借阅。同意学 校向国家有关部门或机构送交论文的复印件和磁盘。 ( 保密的学位论文在解密后适用本授权说明) 学位论文作者签名： 知t 卜h 呤 导师签名： 签字日期： 。睁调椭 签字日期：b r 年1 月1 日 第一章绪论 1 1研究背景及意义 第一章绪论 气液界面处气体传输是广泛存在于自然界与工程技术实践中的问题。为了 提高气体吸附、蒸发、冷凝设备的效率，把握气泡柱、搅拌器中气液相间传质 过程与反应过程之间的关系，人们非常关心影响气液界面间气体传输过程的因 素，化学工程师们进行了大量的有关实验和理论研究，可以说，7 0 年代以前， 气液界面处气体传输研究主要由化学工程界所推动。8 0 年代以来，另外一个社 会需求给气液界面处气体传输研究带来新的推动力，从全球环境角度看，海洋 是温室效应气体最大的汇，通过海气界面，海洋以何等速率吸收大气中的二氧 化碳等影响地球气候的气体? 这种速率与风速、海面水文状况关系如何? 这是 人们十分关心的问题。国际社会对环境保护及全球气候变化的严重关注，使气 液界面气体传输研究出现了一个新高潮。相继召开了数次关于水气界面处气体 传输国际学术会议，不同领域的研究者汇聚一堂，充分展示了当前人们对这个 问题的浓厚研究兴趣及其研究成果。 大气和海洋之间的气体交换是生物地球化学循环中的一个重要链条。另外 气液界面气体传输问题对于环境工程领域也非常重要，例如溶解氧是表征水体 污染程度的重要指标，溶解氧对于枵染物的迁移、转化和降解有着重要的作用， 而大气复氧是水体中溶解氧的重要来源，因此人们对于大气复氧速率的影响因 素非常关心。 总之，界面传质问题非常重要，值得我们去深入的研究。然而至今还不能 建立普适气体传输模型，各种现存传质模型还不能准确反映界面传质的机理， 更深入的机理研究和新模型建立是十分必要和紧迫的。主要原因就是问题本身 的复杂性和相关研究的滞后，尤其是实验技术的滞后限制了界面传质的研究， 因为气液界面处气体传输受控于界面处约几十微米厚度的物质边界层，而要直 接测量物质边界层附近的流场和浓度场是比较困难的。本文想通过数值模拟的 手段对此问题进行研究，从而避开了物理实验中存在的种种弊端。自然界中由 于有风和其它因素的存在，水面总是波动的，所以首要的任务是实现波浪的数 值模拟。 第一章绪论 1 2 波浪数值模型概述 计算机技术和计算方法的迅速发展，为采用数值方法求解波浪运动提供了 可能，采用数值模型研究波浪问题越来越受到人们的重视。与物理模型相比， 数值模型具有较强的灵活性，非常经济( 只需要一台计算机和必需的软件即可 展开工作) ，可以很方便的获得详细的流场信息，减小比尺效应，消除了物理 模型中由于传感器及模型变形等因素的影响。 近年来，国内外学者在波浪数值模型的研究和应用上得到了很大的发展， 通过建立合理的数值模型，求解波浪在不同地形上的传播、变形及破碎“3 ，波 浪与海工结构物的相互作用0 3 ，以及船行波问题等等。根据模型所依据的基本 方程不同，可分为四种：缓坡方程模型，b o u s s i n e s q 浅水短波模型，l a p l a c e 方 程模型和n s 方程模型。n s 方程。1 是流体力学中很著名的数学方程，近年来 发展的数值计算方法，为其求解提供了可能。各种模型所采用的离散化方法主 要有有限差分法，有限元法和有限体积法“，。 有限差分法( f i n i t ed i f f e r e n c em e t h o d ) 是数值解法中最经典的方法。它是 将求解域划分为差分网格，用有限个网格节点代替连续的求解域，然后将偏微 分方程( 控制方程) 的导数用差商代替，推导出含有离散点上有限个未知数的 差分方程组。求差分方程组( 代数方程组) 的解，就是微分方程定解问题的数 值近似解，这是一种直接将微分问题变为代数问题的近似数值解法。 有限元法( f i n i t ee l e m e n tm e t h o d ) 也是广泛应用的流体动力学数值计算方 法，有限元法是将一个连续的求解域任意划分成适当形状的许多微小单元，并 于各小单元分片构造插值函数，然后根据极值原理，将问题的控制方程转化为 所有单元上的有限元方程，把总体的极值作为各单元极值之和，即将局部单元 总体合成，形成嵌入了指定边界条件的代数方程组，求解该方程组就得到了各 节点上待求的函数值。 有限体积法( f i n i t ev o l u m em e t h o d ) 是近年来发展非常迅速的一种离散化 方法，其特点是计算效率高。其基本思路是：将计算区域划分为网格，并使每 个网格点周围有一个互不重复的控制体积，将待解微分方程( 控制方程) 对每 个控制体积积分，从而得到一组离散方程。其中的未知数是网格点上的因变量 中。为了求出控制体积的积分，必须假定中值在网格点之间的变化规律。 有限体积法的基本思想易于理解，并能得到直接的物理解释。离散方程的 物理意义，就是因变量垂在有限大小的控制体积中的守恒原理，如同微分方程 第一章绪论 表示因变量在无限小的控制体积中的守恒原理一样。 有限体积法得到的离散方程，要求因变量的积分守恒对任意一组控制体积 都得到满足，对整个计算区域，自然也得到满足。这是有限体积法吸引人的优 点。有一些离散方法，例如有限差分法，仅当网格极其细密时，离散方程才能 满足积分守恒，而有限体积法即使在粗网格情况下，也显示出准确的积分守恒。 本文中数值模型所依据的基本方程为n s 方程模型，离散化方法采用有限 体积法。 1 3 f l u e n t 软件简介 本文的计算工具为f l l r e n t 软件，它是f l u e n t 公司于1 9 8 3 年推出的 c f d 软件。f l u e n t 软件是目前功能最全面、适用性最广、国内使用最广泛的 c f d 软件之。f l u e n t 软件提供了非常灵活的网格特性“1 ，让用户可以使用 非结构网格来解决具有复杂外形的流动，甚至可以使用混合型非结构网格。它 允许用户根据解的具体情况对网格进行修改( 细化或粗化) 。f l u e n t 软件可 用于二维平面、二维轴对称和三维流动分析，可完成多种参考系下流场模拟、 定常与非定常流动分析、可压缩与不可压缩流动计算、层流和湍流模拟、传热 和热混合分析、化学组分混合和反应分析、多相流分析、多孔介质分析等，( 详 见f l u e n t 软件中相关说明文件”3 ) 。 f l u e n t 软件可让用户定义多种边界条件，如流动入口边界条件、壁面边 界条件等，所有边界条件可随时间和空间变化，包括轴对称和周期变化等。 f l u e n t 软件提供的用户自定义子程序功能”，可让用户自行设定连续方程、 动量方程、能量方程或其它标量方程中的体积源项，自定义边界条件、初始条 件、流体的物性和添加新的标量方程等。 f l u e n t 软件是用c 语言编写的，可实现动态内存分配及高效数据结构， 具有很大的灵活性和很强的处理能力。f 1 1 礓n t 软件可以在w i n d o w s 2 0 0 0 x p 、 l i n u x u n i x 操作系统下运行，支持并行处理。 f l u e n t 程序包包括f l u e n t ( 解算器) 、g a m b i t ( 画物理模型及产生网 格的预处理程序) 、p r e p d f ( 模拟非预混合燃烧的预处理程序) 、t g r i d ( 产生体 网格的预处理程序) 等。f l u e n t 软件中将不同领域的计算软件组合起来，成为 c f d 计算机软件群，软件之间可以方便地进行数值交换，并采用统的前、后 处理工具，这就省却了科研工作者在计算方法、编程、前后处理等方面投入的 重复、低效的劳动，而可以将主要精力集中于物理问题本身规律的探索上。 第一章绪论 用f l u e n t 软件解决问题的基本步骤： ( 1 ) 画出物理模型的图形及划分网格。 ( 2 ) 选择解算器( 二维或三维、单精度或双精度) 。 ( 3 ) 导入网格。 ( 4 ) 检查网格。 ( 5 ) 选择解算器公式( 分离解算器或耦合解箅器) 。 ( 6 ) 选择所需要解决的基本方程( 控制方程) 。 ( 7 ) 说明物质性质( 密度、粘度等) 。 ( 8 ) 说明边界条件。 ( 9 ) 调整解的控制参数。 0 0 ) 起始化流场。 a d 计算出结果。 0 d 检查结果。 o $ 保存结果。 如果需要，细化网格或修改数值模型或物理模型重新进行计算。 1 4本文的研究工作 本文的具体工作如下： ( 1 ) 概括了气液界面气体传输的主要影响因素，总结了自耿膜模型以来出 现的几个重要的界面传质模型，特别是最近的研究成果，并指出了其存在的 优缺点，从而为今后用数值方法研究界面传质问题打下基础。 ( 2 ) 选择了适当的控制方程、边界条件和数值方法，建立了二维数值波浪 水槽。 ( 3 ) 验证了数值波浪水槽，对造波系统和人工衰减层的有效性进行了验 证，对数值耗散和粘性耗散等问题进行了研究，本文还研究了气相压力变化 对波浪的影响。从而为今后研究波动对气液界面气体传输的影响走出了第一 步。 ( 4 ) 用求解n s 方程的方法对平面波肃叶流中小扰动的空间演化过程进行 了数值模拟。并将计算结果与o s 方程的计算结果进行了比较。 ( 5 ) 总结本文的主要结论并对本课题的进一步研究提出建议。 ( 5 ) 总结本文的主要结论并对本课题的进一步研究提出建议。 4 第二章气液界面处气体传输研究概况 第二章气液界面处气体传输研究概况 本章中对气液界面气体传输的主要影响因素进行了总结，然后回顾了自 双膜模型以来特别是近几年来出现的几个重要的界面传质模型，并指出了其 存在的优缺点，最后指出了气液界面气体传输研究的发展方向。 2 1 影响气液界面气体传输的主要因素 影响气液界面气体传输的主要因素有风和雨的影响、污染物的影响( 包 括表面活性剂、浮油花等) 、生物和温度。下面分别讨论这些因素对界面传质 的影响。 2 1 1 风的影响 风对气液界面的影响是非常复杂的，主要表现为风的剪切作用。它不仅 推动水的流动和生成表面波，还会导致风生波的破碎，造成水滴的飞溅和气 泡的卷入，这些现象都对气液界面处气体传输产生重要的影响。 2 1 1 1 风生流 风的剪切作用可以产生近界面剪切流，使界面附近产生强烈的湍动，从而 增大气体的界面传输速度。根据c h i a 等叩3 的实验研究，气体传输速度与风的 剪切速度成正比。 2 1 1 2 风生波 风的剪切作用可产生表面波，d o w n i n g 等1 9 在风波槽实验中观察到：在 风速达到3 m s 对，水面开始出现微波和波纹；与此同时，通过气液界面传入 水中的氧气的传输速度也突然增大。s a y l o r 等o o 也研究了波动对气液界面气 体传输速度的影响，他们使水槽做垂直震动而在水表面处产生毛细波纹，发 现：当波的坡度从0 增至o 2 m m ，气体传入液体的速度增大近两个数量级。 由于实验的独特兴波方式，在水槽中只有毛细波纹的产生，从而排除了风波 槽实验中影响传质的湍动和浮质效应。因此可以得出毛细波纹可大大增加气 体传输速度，证明了表面波在气体传输中的重要性。这可能是由于表面波增 加了气液接触面积，还可能是表面波加强了界面附近的湍流。 箜三童气液界面处气体传输研究概况 2 ，1 1 3 破碎波 波的破碎是界面附近湍流的一个主要能量来源。w a l l a c e 等1 1 1 对大西洋海 湾进行了长达4 个月的实地监测，发现风暴可打破溶解氧的时间连续性。这 是由于在风暴期间波浪破碎，大量的气泡夹杂在海水中，当气泡随着水流运 动到几米处的水深时，在水压的作用下，气泡中的氧气溶于水中，溶解氧达 到过饱和状态。在风暴过后，溶解氧会经过海气界面缓慢的逃逸海洋。 2 1 1 4 气泡 目前气泡对于气液界面气体传输作用还没有得出可以信赖的模型。 k e e l i n g1 ”】认为通过气泡实现的气体传输量关键取决于直径大于t m m 的气泡 量的产率；同时，他也暗示以往忽略大气泡的做法可能大大低估了气泡在气 体传输中的作用。w o o l f ”1 证实了k e e l i n g 的结论。但是根据w a r m i n k h o f 1 的研究，对于二氧化碳这种溶解度相对较大的气体由于气泡作用而增加的气 体传输量甚微，也就是说气泡过程只能大大提高微溶气体的传输速度。 2 1 2 雨的影响 d a v i d 等的研究显示降雨能大大增强气液界面间的气体交换，并且传质 速度与降雨强度和雨滴尺寸有关。导致界面传质增强可能是由于以下原因： 降雨增强了界面附近的水流湍动；雨滴的作用而使水中夹杂着大量气泡； 雨水与海水的温度、盐度不同而产生密度流与分层现象；降水加强了作 用于水面上的风应力；降水分散了表面活性物质微层。 2 1 3 污染物的影响 污染物的存在对于气液界面的气体传输有很大影响。z h o u 等 16 】对去离 子水、自来水和纸浆废水三种水质中的臭氧传输进行了研究，结果表明当使 用纸浆废水时，由于臭氧与废水中污染物的化学反应消耗掉了气液膜中的溶 解臭氧，从而使气液界面传质增强。氧气的氧化能力虽不及臭氧，但可以推 测气液界面间氧气传输会受到溶解氧参与的化学反应和生化反应的影响。 v a z q u e z i l n a 等【1 7 实验研究了液相粘度对传质的影响。他们认为增加液相 粘度有正反两个方面的作用：方面可使实验装置中激发点的振幅增大，波 形稳定性提高，从而使产生大振幅波成为可能，增强了传质速度；但另一方 面也增加了阻尼系数，使得振幅随着波动远离激发点而迅速减小，这将不利 第二章气液界面处气体传输研究概况 于传质。实验结果显示：随着粘度的增加，传质速度下降。因而，由于污染 物可以改变海水的粘度，而影响气液界面间的氧气传输。 污染物还可以通过其它途径影响气液界面的传质，例如：石油在水中可 以以浮油形式存在，浮油漂浮于水面，易扩散成油膜。当油膜厚度大于1 微 米时，可以隔断气液界面间的气体交换，并可大大降低界面附近的湍动，同 时，石油自身的氧化要耗氧，这些都影响气液界面间的传质。 2 1 4 生物的影晌 庄万金“”根据厦门西海域溶解氧和浮游植物数量周年变化的实测资料得 出：溶解氧含量跟浮游植物数量成正相关性，即溶解氧含量随着浮游植物细胞 数的增殖而上升，随着浮游植物的消亡而下降。这种正相关性在赤潮发生期间 表现的更加明显。由此可见，海洋生物通过影响溶解氧浓度改变界面两侧的浓 度梯度，从而影响界面间的气体通量。 生物还可以通过其它途径影响界面传质。f r e w 等“”发现在海洋表面富集 的表面活性物质微层能降低海气界面的气体交换速率。而z u t i c 等。”研究认为 海洋中表面活性物质主要来源于浮游植物渗出液及它们的降解物质。 2 1 5 温度影响 界面传质系数与温度一般符合b e n n e t t 和r a t h b u n 。”提出的关系式： k r = k 2 0 0 ( ”) ( 2 1 1 ) 其中，k 0 为温度t 时的传质系数，k 2 0 为2 0 时的传质系数，o 为常数，一 般取1 0 2 4 1 。 除此之外，温度还可通过影响污染物的生化降解速度和生物数量等对气液 界面间的氧气交换产生影响。 2 2 早期传质模型研究 w h i t e m a n 2 2 】提出了目前广为使用的双膜模型。模型假定在界面两侧各有 一个薄膜，气体在两膜层中的传质方式为分子扩散，界面间的传质过程阻力 就是气体通过双膜层的阻力。由于氧气在液相中的传输阻力远大于其在气相 中的传输阻力，故称为液相控制的气体传输。研究结果显示：当分子扩散方 式为组分a 通过停滞组分b 的扩散时， 第二章气液界面处气体传输研究概况 世。：兰殳 ( 2 2 1 ) z l cb 其中k l 为液膜内的传质系数，d 分子扩散系数，c 。混合物的总平均物质 的量浓度，z l 为液膜的厚度，c b m 为停滞组分b 的对数平均浓度。用上述模 型描述具有固定相界面及流动速度不高的两流体间的传质过程与实际情况相 符较好，但对于湍动相当厉害时的情况，此模型预测结果不好。此外，此模 型未能涉及膜层外湍流与膜层厚度的关系，所以不能直接应用于复杂传质过 程。 为了弥补双膜模型的不足，h i g b i e 2 3 1 提出了溶质渗透模型。他认为在湍 流作用下，内部流体单元运动到爨面，经过暴露时间( 0 。) ，又被排挤回到流 体内部。当流体单元扩散到界面时，气体不稳态扩散到液膜中。经过推导， 得到传质系数k 的表达式为： 厂元一 世= 2 f 导 ( 2 2 2 ) 、f 礁 一。 其中，0 。为暴露时间。溶质渗透模型虽考虑了非稳态传质过程，但假定流体 微元暴露时间为常数，忽略流体微元年龄的可变性。 后来，d a n c k w e r t s 。”提出了表面更新模型。他假定液体表面是由不同暴 露时间的液面单元构成，不管它们的年龄如何，被置换的概率相等。在假定 流体微元年龄分布函数中( 日) 满足m ( 口) = s e “8 分布后，得到传质系数k 的 表达式为： k ：丽( 2 2 3 ) 其中s 为表面更新率，e 为年龄。 溶质渗透模型和表面更新模型都正确导出了k o c d ”2 的结论，也合理地把 湍流活动的影响模式化为表面更新作用，但是在应用这两个模型时，很难确 定暴露时间0 。或表面更新率s ，使得模型难于实际应用。 6 0 年代末至7 0 年代初，为了合理计及湍流的作用，f o r t e s c u e 和p e a r s o n 2 5 提出了大涡模型，随后，l a m o n t 和s c o t t1 2 6 1 提出了小涡模型。大涡模型认 为：湍动流体是由波长不一的旋涡构成，它们的平均传质特性可以由一系列 规则的旋子( r o l lc e l l ) 来模拟，这些旋子的长度尺度( a ) 可取界面下湍动流场的 积分长度尺度，也就是含能涡的长度尺度，即假定大尺度涡控制界面传质过 程。从而导出传质系数k 为： k = 1 4 6 ( d 、 u2 a 1 1 72 ( 2 2 4 ) 第二章气液界面处气体传输研究概况 其中，u 为湍动速度。与其它模型和实验值的比较发现，大涡模型理论预测 值比实验值高1 1 - 1 8 ，比溶质渗透模型预测值高3 0 ，结果说明大涡模型提 供了一个预测的上限值。 与大涡模型不同，小涡模型假定界面附近的大尺度涡上熏叠着许多小尺 度涡，这些小涡控制着界面附近的局部混合，并且进而假定这种局部混合控 制着整个传质过程。这些涡是粘性耗散的k o l m o g o r o v 涡。导出的气体传质系 数为： k = 0 4 0 ( s y ) 1 7 4 s c 一1 7 2 ( 2 2 5 ) 其中，s c 为s c h m i d t 数，湍流的能量耗散率，u 粘性系数。 以上两个涡模型在概念上使表面更新与湍流概念建立了联系，但是如何 确定参数以及到底是什么尺度的涡控制着传质过程等问题还没有得到解决。 t h e o f a n o u s 等认为大涡模型和小涡模型分别适用于不同的湍动雷诺数 ( r l ) 范围，并估计出大涡模型的适用范围：r l 亭时 “= 0 ( 3 2 4 ) v = 0 ( 3 2 ，5 ) 其中孝= a c o s k ( x - c t ) ，知振幅，d 为水深，c 为波速，k 为波数，。为圆频 率，e 为波面。但是如果观察式( 3 2 2 ) 和( 3 2 3 ) 可以发现，造波板应具有沿水 深变化的速度，这在实际水波实验中很难实现，是一种理想中的造波方法。 第三章二维数值波浪水槽的建立和求解 f u r s e l l 和r g d e a n 【39 】提出了利用水平振动平板产生水波的推板造波理 论，但此方法在处理计算域内存在障碍物时遇到了困难。b r o r s e n 等“”提出了 源造波方法，即在计算域内设置一个造波源，其源项等于生成波相应的水平 速度，在源两侧同时产生方向相反的两列波，源项处可透过波浪遇障碍物形 成的反射波，解决了上述困难。高学平等“”利用同样的造波方法建立了不规 则数值波浪水槽。王永学“2 1 基于线性造波理论应用v o f 法给出了可吸收造波 机数值边界条件，即造波板的运动除了产生行进波外，同时还产生一个抵消 反射波的局部运动。p e t e r 等“”通过监测计算域内反射波信息，随时调整入射 波的大小，而建立了具有消除二次反射功能的造波机数值边界。王永学和p e t e r 等的方法也可以有效处理波浪遇障碍物形成的反射波。 笔者曾试用过王永学的可吸收式数值造波机作为入流边界条件，根据文 献 2 中的计算结果和笔者经验可知，采用这种入流边界的数值水槽需要较长 的计算区域( 因为水槽上游需要数个波长距离作为过渡区域) ，这样会造成计 算区域增大，网格数目多，从而使计算代价增大。另一方面，可吸收式数值 造波机的主要优点是消除波浪遇到计算域中构筑物的反射在造波机边界造成 的二次反射，而本文计算域中没有构筑物，故而没有采用此法。为了减小计 算时间，本文中采用了f u r s e l l 和r g d e a n 3 9 提出的推板造波理论，其主要 原理为： 当造波板在某点处作周期性运动，其水平位移为： x ( t ) = s o 2 xs i n o ) t 其中s o 为造波板的水平最大位移，u 为造波板的运动频率 论，在水深为d 的水池中，波面方程为： ( 3 2 ，6 ) 根据线性造波理 f = s 2 0 、2 埘4 s “i n h l l h 2 ( k ( d 2 而) c 。s ( 缸c o t ) ，+ 4 s i n 2 ( “d ) 万- u 。s i n ( c o t ) ( 3 其中堙t a n h ( 删) 一国2 = 0 ( 3 2 8 ) u 。gt a n ( b 。d ) + 0 9 2 = 0 ( 3 2 9 ) 其中式( 3 2 7 ) 中第二项在离开造波板约3 倍水深后其值很快衰减，可以忽 略。在距离造波板较远处，波面方程简化为： f = s 2 0 、2 矧4 s + i n s h i 血( 2 ( k d 2 蕊) c 。s ( 勉一鲥) = 鲁c 。s ( 缸一研) ( 3 2 1 0 ) 因此，造波板水平位移的最大值同距离造波板3 倍水深以外的小振幅波波高 第三章二维数值波浪水槽的建立和求解 之比为： i h = 瓦鬲4 s i n 丽h 2 k d ( 3 2 1 1 ) 2 柳+ s i i 】1 l ( 2 榭) 本文中为了保证计算的稳定性，避免计算开始产生大波，在造波机边界 条件中( 即入流速度和相应的波面值中) 加入了一个双曲正切函数，使入射 波在前三个周期内缓慢增大。 3 2 2 波浪开边界的无反射处理 为了使波浪能顺利穿过开边界而不在开边界处产生反射，许多学者对开 边界间题进行了深入的研究。对于小振幅的线性波浪问题，可以采用 s o m m e r f e l d 辐射边界条件处理，即； 掣+ c 掣：0 ( 3 2 1 2 ) o t田 其中q 为变量( 例如速度，波高等) ，c 为波速。但在实际计算中，由于波浪 的非线性。通常波速不易确定，而使得该条件不能完全顺利放出波浪，从而 在开边界处产生反射，对计算结果有一定的影响。 同时需要指出的是单独使用s o m m e r f e l d 辐射边界条件很多时候不能完全 起到无反射辐射波浪的作用。即使对于小幅值不规则波浪，由于其波速不好 确定，s o m r n e r f e l d 辐射边界条件也不好应用。一种有效的办法是在使用 s o m m e r f e l d 辐射边界条件前，设置人为衰减层，消减大部分波能后，再利用 s o m m e r f e l d 边界条件放波或直接消除所有波能。人工衰减层又称海绵层。 r o m a t e l 4 4 】对表面波的能量吸收方法进行了较为全面的总结。刘应中等【4 5 ) 通过在出口边界前设置阻尼耗散区的方法消除波能。其基本思想为让垂向速 度在耗散区内受到阻尼，该阻尼沿波浪传播方向从零开始按二次曲线关系增 大，而在垂直方向从海底的零向上线性增加。最后，为了消除残余波，在出 口边界使用s o m m e r f e l d 辐射条件。刘长根 2 】采用了海绵层吸收处理，即在海 绵层内对速度和波面进行消减，而对压强没有进行消减。具体消减方法为： u = “( x ) ，v = v r ( x ) ，叩= 1 7 , u ( x ) ，k = k “( x ) ，g = 占( z ) ，p = p ( 3 2 1 3 ) 其中( x ) 消波系数，k 和e 分别为湍动动能和湍动耗散率。 本文中开边界采用无滑移墙壁，在墙壁前采用消除波能的方法，目的是防 止波浪的反射。采用文献 4 1 和【4 6 中提出的消波滤波器方法，即在出流边界前 第三章二维数值波浪水槽的建立和求解 设一倍波长的海绵层消波段，在该段内动量方程中加进源项一uv ，其中u 为消 波系数。消波系数沿波浪传播方向取线性分布，即 f l = a + 芝羔( 3 2 1 4 ) 山 其中，x 。海绵层起始处的水平坐标，x 海绵层水平坐标，l 海绵层厚度，a 为常 数，可以理解为吸收强度。吸收强度的取值会对波浪的反射程度有影响，它 的确定详见4 1 2 。 笔者曾尝试用人为加大粘性系数的方法消除波能，即在人工衰减层中以 线性或其它函数形式增大粘性系数，从而达到耗散波能，然而很难在很短距 离内消除波能而又不产生波反射，因为只有大粘性值才能显著地消除波能量， 这可以从图4 2 4 中看出。实际操作中，人工衰减层中必须设置很大的粘性系 数( 通常是有效计算区域内粘性系数的数十万倍) ，而f l u e n t 又不允许有 大梯度的粘性系数出现，导致粘性增大到其值的数十万倍势必需有较长的距 离，造成计算区域很大，故这种方法没有成功。 3 2 3 其它边界条件和初始条件 数值水槽底部边界设置为无滑移边界条件，在此边界应满足 “= o ，v = 0 ( 3 2 1 5 ) 上边界设为等压面，这一点与以往的单相流数值波浪水槽不同。自由表 面不再是边界条件，无需在自由表面处设置边界条件。初始时刻计算区域内 的速度和压力均设为0 ，水面为静水水面。 3 3自由表面的跟踪方法 为了追踪自由界面的运动，本文采用v o f 法。在计算域内的每一点上定 义一个函数f ，f 是时间和空间的函数，当该点被流体质点占据时，此点处f = - i ， 反之f = - 0 。在计算域离散为计算网格后，在二维情况下，体积分数f 为一个单 元内的流体面积与该单元体面积之比，即 1， f ( z ，y ，f ) 2 焘j j , ( q ，r ，t ) d c d , 7 ( 3 3 1 ) 如果f = i ，说明该单元充满流体；如果f = 0 说明该单元不舍该流体，即充 满气体，为外部单元；如果f 的值在0 与1 之间，说明该单元是含有自由面 单元。通过体积分数的输运方程可跟踪界面，即 第三章二维数值波浪水槽的建立和求解 睾v 吁o a ( 3 3 2 ) 其中f 2 为第二相的体积分数。通过求解f 的值，就可以确定自由面在那些单元 内。另外，利用f 值沿自由面法线方向变化最快，用f 的梯度最大值可以确 定自由面的方向。v o f 法的基本思想是通过构造一个流体体积函数来跟踪每 个单元网格内的流体量，并根据这个函数值和其导数值构造自由面形状。在 采用了流体体积函数方法后，计算运动自由表面问题转化为计算整个计算域 内的f 值问题。 f 函数是一个阶梯函数，不是连续函数，因此方程不能用通常的差分格 式进行离散求解，需要采用自由面重构的方法来处理自由面。本文中自由表 面重构采用分段线性体积跟踪方法“1 。它采用分段线性逼近自由面，即在每一 个控制体内，用一个斜线段表示自由面，它的具体做法是先根据一个控制体 内的流体体积函数f 的梯度值确定自由面的方向，然后根据此控制体积内的 流体体积函数值，确定自由面在此控制体内的位置。对于在一个时间步长内 的流体通过控制体某一边的体积流量是采用类似于l a g r a n g i a n 方法确定的， 而不是采用直接差分方法求解f 的控制方程确定的，根据在一个时间步长内 流过四周边界到相邻网格中的流体体积量，修改本网格和相邻网格的流体体 积函数值，从而确定下一时刻的自由面形状。 本文中采用k o t h e 等。1 的方法计算流体体积函数f 的梯度。如图3 2 孙+ 一 i + 1 。j + 1 缓j ! | _ t i 黝孙 i + 1 。j 缓 戮小， i 1 ，j - 1 - 施 2 匮3 率 图3 2 自由表面和自由表面在自由面单元内的形状 对于含有自由面的网格( i j ) ，其自由面的法线方向矢量分量为： n ，= d 1 ( + ，一+ - ) + 2 ( 一) + ( f + i , j _ 1 - - 一，) ( 3 3 3 ) ，= 石1 阮，- + 2 鼻+ - ) 一( 一+ 2 e ，t 一一，) 】 ( 3 川 从而对于网格( 列) 内的自由面与水平方向的夹角可以写为 第三章二维数值波浪东槽的建立和求解 卢= a r c t a n ( 一竺苎) ，( 一霈 p 万) n y 自由面在一个单元内的形状有1 6 种情况 向的夹角为a ，且 n = a r c t a n ( 芸- yt a n p ) ，( o 0 时，则f l u x ( r ) = 以a t ( f , 缈一寺u , a t t a n f l ) 当u b 0 时，若u b a t s e 缈，则f l u x ( b ) = u 6 a t l x x 1 若a t s c a y ，则f l u x ( b ) = u b m a x 一妄( 缸一巧妙) 2 c t g f l 否则f l u x ( 6 ) = f j 每 1 当 0 时，则f l u x ( o = u t 出( 耳缈+ 寺r s , a t m f l ) 二 用同样方法可以求出其它三种情况下通过网格边上的流量，对于时间层 ( n + 1 ) 上，计算网格的流体体积函数值为时间层( n ) 上的流体体积函数值 减去通过网格的四条边上的净流量。进而可以用上文中的方法重构出自由面 的位置和方向。 第三章二维数值波浪水槽的建立和求解 3 4 控制方程的离散和数学模型的求解 f l u e n t 软件中数值方法分为三个步骤： ( 1 ) 将计算区域分成离散的控制单元； ( 2 ) 将控制方程在每个控制单元上积分，产生离散因变量( 例如速度、压 力) 的代数方程： ( 3 ) 将这些代数方程线性化，求解线性方程组得到因变量的值。 控制方程的离散采取了有限体积法，在此以二维瞬态对流扩散问题为例， 说明这个离散过程。控制方程的通用形式： 挈+ 掣+ 掣：昙( r 罢) + j a ( 、f a 。o l ，+ s ( 3 4 1 ) a t 敏 卸a x ：瓠。a v 、a v 这里，中是广义变量，厂是相对于由的广义扩散系数，s 是与巾对应的广义源 项。对于动量方程，我们把压力梯度项暂且放在源项s 中去。 我们使用图中所示的计算网格来划分整个计算域，网格中实线的交点是 计算节点，由虚线所围成的小方格是控制体积。将控制体积的界面放置在两 个节点中间的位置，这样每个节点由一个控制体积所包围。 h | n w 一。钐孑笏一 e ) 、f 殇荔渗缓：l 1 s l 圭上一 f 耳司 图3 4 计算网格及控制体积 我们用p 来标识一个广义的节点，其东西两侧的相邻节点分别用e 和w 标识，南北两侧的相邻节点用s 和n 标识，与各节点对应的控制体积也用相 应字符标识。图中阴影部分表示节点p 处的控制体积p ，控制体积的东西南北 四个界面分别用e 、w 、s 和n 标识。控制体积在x 和y 方向的宽度分别用x 和y 标识，控制体积的体积值a v = 血缈。针对图中的计算网格，在控制 箜三童三丝塑焦鎏遗三! 三焦笪塞皇塑塞堡 体积p 及时间段a t 上积分控制方程，有 ri ，挚附+ rl v d i v ( p u ( b ) d v d t = 厂l ，v ( f g r a d 庐) d v d t d i ( f g r a d 庐) d v d t + r l v d t 2 jl ，v厂 ( 3 4 2 ) 方程( 3 4 2 ) 中左边第一项和第二项分别称为瞬态项和对流项，方程右 边第一项和第二项分别称为扩散项和源项。在处理瞬态项时，假定物理量中 在整个控制体积p 上均具有节点处的中，同时假定密度p 在时间段at 上的 变化量极小，则( 3 4 2 ) 中的瞬态项变为： r f ，挚m = 叫p d t d y = 群( 砟一卵) 矿 ( 3 4 3 ) 上式中上标0 表示物理量在时刻t 的值，而在什t 时刻的物理量没有用上标 标记，下标p 表示物理量在控制体积p 的节点p 处取值。源项变为： f +