（流体力学专业论文）床面上直立圆柱的三维湍流及局部冲刷数值模拟.pdf

摘要 当桩柱置于流场中时，它的存在将改变其附近的流动情况，馒其 附遮的流动速度鸯籍快，密诧弓 超一系列的变纯，造成桩柱附近的局部 冲刷，并由此丽危及诸如桥墩等建筑物的稳定。本论文就是依此必背 景，分橱了桩柱附近的流场遥动情况并计算了由此而造成的局部冲蒯 情况。 厂 7 在对溺题的处瑾过程中，我们将流场计算同底部泥沙的运动分开 ( 进行，由流场为床丽的泥沙运动提供相应的条件，在计算褥到床霆寒 程交仡以后再返回流场重新进行流动的计算。在处理的过程中考虑到 计算量的问题，本文取计算流场的时间步进不同予计算床垂变形的时 间步进，也就是说我们采用了半耦合的计算方法。y 绕桩柱三维流场的计算本文采用基于有限体积法( f i n i t ev o l u m e m e t h o d ) 的压力修正s i m p l e 算法。其中的湍流模式采用了 w i l c o x ( 1 9 9 3 ) k 一国模式。在求解离散的线性方程组的时候，本文采用 s t o n e ( 1 9 6 8 1s i p 法( s t r o n g l yi m p l i c i tp r o c e d u r e ) ，它是基于不完全l u 分解的一种迭代求解方程组的方法。本文分析了光滑床霆和粗糙床疆 两种不嗣的条件下的流动情况，并将计算得到的结果同前入的实验和 数值计算结果进行了比较。本文对反映流动信息的些主要的参数 ( 麴速度，压力，分离角，床谣切应力) 进行了详细的分析，得到了许 多有意义的结果。在计算过程中我们还对粗糙度对流场的影响进行了 一定豹分析。 - l 搏，j ，簪鼍j 嘴天 r 床面地形变化的数值模拟是最近几年新出现的一个研究方向，由 i 于泥沙运动的复杂性，要进行准确的数值模拟现在仍然有一定的难 度。本论文的数值模拟结果与前人的实验结果作了初步的比较。在泥 沙输运的计算过程中只考虑床面等粒径沙的推移质输运，忽略悬移质 的影响，采用了e n g e l u n d & f r e d s o e ( 1 9 7 6 ) 的两维输沙率计算公式。 床面地形的变化我们采用两维e x n e r 泥沙连续性方程。另外在计算过 程中考虑到泥沙水下休止角的影响，本文还进行了s a n d s l i d e 的计算。y 7 ” 关键词：桩柱，s i m p l e ，s i p ，k 一模式，局部冲刷，泥沙运动 i i t h r e e d i m e n s i o n a ln u m e r i c a lm o d e lf o rf l o wa n dl o c a l s c o u ra r o u n dab o t t o m - m o u n t e dp i l e w h e na p i l ei sp l a c e d i na s t e a d yc u r r e n t ，t h ep r e s e n c eo f t h ep i l ew i l lc h a n g et h e f l o w p a t t e r n i ni t s n e i g h b o r h o o d ，r e s u l t i n g i nt h e f o l l o w i n gp h e n o m e n a ：t h e c o n t r a c t i o no ft h ef l o w , t h ef o r m a t i o no fah o r s e s h o ev o r t e xi nf r o n to ft h ep i l ea n d l e e w a k ev o r t i c e sb e h i n di t t h e s ec h a n g e sn s l l a l l yc a u s ea ni n c r e a s ei nt h el o c a l s e d i m e n tt r a n s p o r tc a p a c i t y , a n dt h u sl e a dt ol o c a ls c o u r f o rt h es o l u t i o no f t h i sp r o b l e m ，t h ec a l c u l a t i o nf o rf l o wf i e l di sd i v i d e df r o mt h e m o r p h o l o g i cm o d e l c o n s i d e r i n gt h ec o m p u t i n gc o s t , w eu s es e m i - c o u p l e dm e t h o di n c a l c u l a t i o np r o c e d u r ew i t ht w ok i n d so f t i m es c a l e s t h et h r e e - d i m e n s i o n a lf l o wa r o u n dav e r t i c a lb o t t o m - m o u n t e dc i r c u l a r p i l e s t u d i e dn u m e r i c a l l yb ym e a n so fs i m p l em e t h o db a s e do nt h ef i n i t ev o l u m e m e t h o d t h ek 一t u r b u l e n c em o d e li su s e df o rc l o s u r e n l es i pi t e r a t i v em e t h o d i sa d o p t e df o rt h es o l u t i o no ft h el i n e a re q u a t i o ns y s t e m t h ec u l c u l a t i o nr e s u l t sa r e v a l i d a t e da g a i n s te x p e r i m e n t a ld a t aa n dn u m e r i c a lr e s u l t sf r o m p r e d e c e s s o r , w e d i s c u s st w oc a s e s s m o o t hb e da n dr o u g hb e d s o m ef l o wi n f o r m a t i o nl i k ev e l o c i t y , p m s s u r e s e p a r a t i o na n g l ef o rl e e - w a k ev o r t i c e sa i ea n a l y z e d i nt h et h e s i s a m o r p h o l o g i cm o d e li sd e v e l o p e dw i t hat w o - d i m e n s i o n a lb e d l o a ds e d i m e n t t r a n s p o r td e s c r i p t i o n w eo n l y c o n s i d e rt h eb e d l o a dm o v e m e n ta n d i g n o r e t h e s u s p e n d e ds e d i m e n t t h ee n g l u n df o r m u l a t i o ni sa p p l i e df o rt h eb e d i o a dd e s c r i p t i o n t h ec h a n g eo fb e de l e v a t i o ni so b t a i n e df r o mt h ee x n e re q u a t i o n t h ee f f e c to f s a n d s l i d ej sa l s ot a k e ni n t oa c c o u n t k e y w o a d s ：p i l e , s i m p l e ，s i p , k mm o d e l ，l o c a ls c o u r i i i 上海交通大学 学位论文原创性声明 本人郑重声明：所呈交的学位论文，是本人在导师的指导下 独立进行研究工作所取得的成果。除文中已经注明引用的内容外，本 论文不包含任何其他个人或集体已经发表或撰写过的作品成果。对本 文的研究做出重要贡献的个人和集体，均已在文中以明确方式标明。 本人完全意识到本声明的法律结果由本人承担。 学位论文作者签名： 日期：年月日 上海交通大学 学位论文版权使用授权书 塞学位论文作者完全了解学校有关保留、使用学位论文的规定， 同意学校保留并向国家有关部门或机构送交论文的复印件和电子 版，允许论文被查阅和借阅。本人授权上海交通大学可以将本学位 论文的全部或部分内容编入有关数据库进行检索，可以采用影印、 缩印或扫描等复制手段保存和汇编本学位论文。 保密口，在一年解密后适用本授权书。 本学位论文属于 不保密口。 ( 请在以上方框内打“4 ”) 学位论文作者签名：指导教师签名 日期：年月 日日期：年月 日 j ；嗣 交l 大k 习l 士- 擘文 物理量的说明 矢量 一维及两维数组 反映桩柱排列情况的因子 经验常数( = l o ) ，力矢量的矢径大小 桩柱的宽度，力矢量的矢径大小 拖拽力系数 压力系数 泥沙颗粒的直径( = d 。：泥沙级配中以重量计5 0 较之为小的粒径) 参考系数( = 1 0 5 ) 水流拖拽力 动摩擦阻力 上举力 尾涡脱落的频率 重力加速度 两种坐标系下的基矢量 冲淤深度，水深 c a r t e s i a n 坐标系下的三个单位基矢量 湍动能 斜面切线方向 n i k u r a d s e 平衡泥沙粗糙度 墙面单元上的平衡n i k u r a d s e 泥沙粗糙度 混合长度，进口段长度 泥沙质量 4 趔一 口。 q d 脚 工 g 榭。讲。f t k ， 坍 湍流度 泥沙孔隙率 斜面法线方向 r - y 压j j ( 等于流体压力减去静压力的分布：p = p + p 。一庸- f ) 粘性运动流体的压力函数 压力常数 泥沙起动概率 底沙的输沙率 圆柱r e y n o l d 数 矢径 冲刷深度 平衡冲刷深度( 最大冲刷深度) s t r o u h a l 数 泥沙相对密度( = 乡= 2 6 5 ) 桩柱的形状因子 冲刷时间尺度，变化流场时间尺度 无量纲时间尺度 时间( 底面地形变化尺度) 平衡冲刷时间 均匀来流速度 底面泥沙颗粒的运动速度 泥沙的起动流速 水流的摩阻流速 三个方向的水流速度 v i p p 风 p 玩 如 f s & s 。 曲 丁 r ， 一 玑 u u 匕海j ：，l 大掌囊j 呻e 文 u 。 u ， u w 删彻 x ，y ，z 口 a s e d 8 盯g “ p d ， 法向速度分量 相对速度 切向速度分量 相对重力( = 冬昭0 1 ) a 3 ) o 离散参数( = 2 或3 ) c a r t e s i a n 坐标系下的三个坐标 摩阻流速和斜面最倾斜方向之间的角度 柱侧面尾涡分离角 斜面的倾斜角 摩阻流速和泥沙运动方向的夹角 泥沙拣选系数 动力粘性系数 泥沙运动的动摩擦系数 泥沙运动的静摩擦系数 涡粘性系数 s h i e l d s 水流强度参数 斜面上临界s h i e l d s 水流强度参数 平面上临界s h i e l d s 水流强度参数 水的密度 泥沙密度 边界层厚度 泥沙运动的休止角( 3 2 0 ) 应力张量 蜥 臼 眭 p n 艿 0 j - 嗣 交丑大捌u 曩士埯文 底面切应力 无柱情况下的床面切应力的大小 特定湍流耗散率 湍流耗散率 插值函数 运动粘性系数 f 国 f 五 y 第一章前言 第一节论文背景 冲刷对河流和海流下的结构物的稳定会造成很大的威胁，有关水中结构物冲 刷的防护和维修的费用也是非常高的。 冲刷过程同流动情况和底床床面的性质有直接的联系。当桩柱( 如桥墩) 置 于水流中，它的存在将改变其附近的水流结构( 见图1 1 ) 。主要表现在：桩柱前 面的水面涌波，桩柱附近流动的压缩，桩柱前面的马蹄涡及后面的尾涡( 分为脱 落和不脱落两种情况) ，紊动的加强等。在桩柱的周围将出现尺度很大的旋涡体 系，旋涡的中心形成负压吸起床面的泥沙带入下游。这些变化将会引起当地水流 冲刷能力的增强，桩柱附近尤其明显( = 2 4 ，e t t e m a , ( 1 9 9 8 ) ；= 2 3 ， r a u d k i v i ，( 1 9 8 6 ) 。其中& 代表冲刷的最大深度，d 代表桩柱的直径) 。冲刷对结 构物的稳定性将造成巨大的威胁。 冲刷可以分为两种类型：清水冲刷和动床冲刷。在清水冲刷中，远离结构物 的泥沙由于口 以，在整个床面上的泥沙几乎都将 运动起来。冲刷还可以分为另外两种类型：整体冲刷和局部冲刷。整体冲刷反映 整个区域床面泥沙的运动情况，局部冲刷则着重体现由于结构物的存在而造成的 其附近的冲刷。我们所考虑的问题主要是桩柱附近的冲刷，所以我们将要讨论的 问题属于：动床局部冲刷。 水流作用下的桩柱附近的局部冲刷问题前人已经做了大量的工作( b r e u s e r s & r a u d k i v i ，1 9 9 1 ；w h i t e h o u s e1 9 9 8 ；m e l v i l l e1 9 9 9 ) 。他们主要是以实验的手段来 估计冲刷坑的形态及最大冲刷深度等问题。 水流作用下造成冲刷的最主要因素是马蹄涡的存在。由于逆压梯度的存在， 产生下降水流，边界层分离，并产生马蹄涡系。b a k e r ( 1 9 7 9 ) 利用可视技术对层 流马蹄涡的结构进行了细致的研究，分析了马蹄涡系中三个顺时针，三个逆时针 共六种涡运动变化的情况。之后，b a k e r ( 1 9 9 0 ) 又对马蹄涡系的振荡及由层流到 紊流的过渡进行了研究。d a r g a h i ( 1 9 8 9 ) 在实验水槽中对圆柱附近的流场进行了 详细的分析，讨论了马蹄涡和尾涡的一些性质，并分析了数同涡的数量及维 数之间的关系，同时通过分析拓朴结构，作者讨论了在床面附近马蹄涡对尾涡分 离角的影响。 由于整个流动是三维的，这样的流动将产生一个非常复杂的冲刷过程，直到 现在我们仍然很难用数值方法令人满意地完成对整个流动及冲刷过程的模拟。但 是由于没有处理好实验过程和实际情况中模型比尺的问题，由实验所得结果换算 成的冲刷坑的深度往往要比实际河流情况要深( e t t e m a e ta 1 1 9 9 8 ) ，造成实验结果 不可靠。为了更好地预测局部冲刷的地形情况，人们逐渐引入了三维流动的水动 力学数值模型。在这方面的工作首先是从对马蹄涡的数值模拟开始的。b r i l e y & m c d o n a l d ( 1 9 8 1 ) 用数值方法计算了平板同圆柱之间的层流马蹄涡。d e n g & p i q u e t ( 1 9 9 2 ) 用数值方法捕获到了马蹄涡的主要性质。 图1 1 圆柱附近的流动情况 f i g u r e l 1f l o w a r o u n dab o t t o m - m o u n t e dp i l e 上蔫 曩 ，l 大掌习l 士论文 对管线附近的冲刷闯题的研究早在八十年代就由w i n d 等人展开了。 l e e u w e n s t e i n & w i n d ( 1 9 8 4 ) 利用七一占两方程模式对管线近的冲刷进行了二维 的数值模拟，v a nb e e k & w i n d ( 1 9 9 0 ) 对同样的问题又进行了研究。而对圆柱 附近的局部冲刷问题的数值模拟工作首先是由o l s e n & m e l a a e n ( 1 9 9 3 ) 完成的。 他们利用有限体积法( f i n i t e - v o l u m em e t h o d ) 解决了三维定常清水冲刷问题。其 中的湍流应力项是用k 一占湍流模型来实现的，所得到的冲刷坑及相应的一些参 数同实验结果吻合的非常好。但最终的冲刷深度( 即平衡冲刷深度) 没有得到， 后来o l s e n & k j l l e s v i g ( 1 9 9 8 ) 利用七一占模型得到了动床冲刷的平衡深度。另外 r i c h a r d s o n & p a n c h a n g ( 1 9 9 8 ) 利用冲刷坑预测的方法从另一个途径来分析和求 解了此问题。在文中作者首先假设冲刷坑已经存在，采用f l o w 3 d 的计算流体 动力学模型，利用v o f 方法( v o l u m e o f - f l u i dm e t h o d ) 求解三维非定常 n a v i e r - s t o k e s 方程，其间利用了p r a n d t l 混合长度理论和双方程的k g 模型。 最近几年以来，丹麦技术大学的水力和水资源系( i s v a ) 在此方面做了大 量的工作。t j e r r y ( 1 9 9 8 ) 利用二维湍流模式d u n e 2 d 对冲积河流中的沙丘进行了 研究。a n d e r s e n ( 1 9 9 9 ) 将由w i l c o x ( 1 9 9 3 ) 年引入的k 一湍流模式用于d u n e 2 d ，改 变了原来模型中采用的k s 模式。r o u l u n d ( 2 0 0 0 ) 在其博士论文中也对圆柱附近 的三维流动及其造成的冲刷问题进行了详细的研究和分析。文中作者利用 e l l i p s y s 3 d 流动求解器及k 一0 3 ，s s t 湍流模型对整个流场进行了求解和分析， 并将自己所得到的结果同前人及自己的实验结果进行了分析比较，得到了满意的 结果。在文中他还对影响冲刷深度的各种因素进行了分析，得到了主要的影响因 素( & 数) ，最后作者的数值模拟结果得到了很好地反映了冲刷坑的形状，结果 令人振奋，也为我们今后的研究提供了许多宝贵的借鉴。 振荡流或波浪作用下的局部冲刷问题是一个新的研究方向。m a c h e m e h l ( 1 9 7 5 ) 在此方面采用实验的方法进行了旱期的研究，之后在此方向的研究日渐 增多。在这一方面s u m e r 做了大量的工作( s u m e r e ta 1 1 9 9 2 ；s u m e re ta l1 9 9 9 ) 。 人们首先是从对振荡流的分析入手的( 如k o z a k i e w i c ze t a l 1 9 9 4 ) ，得到了许多振 荡流条件下的结果，分析了波浪作用下局部冲刷的各种因素对冲刷过程的影响， 得到了主要的影响因素( k c 数) ，并得到了冲刷深度同k c 数之间的关系。在s u m e r ( 1 9 9 3 ) 一文中，他又分析了三种不同的桩柱形状( 圆形，9 0 0 正方形，4 5 0 正方形) 上海曩：噩大掌习【t 呻e 文 对波浪作用下的局部冲刷的影响，并分别得到了它们的冲刷深度同k c 数之间的 关系。在s u m e r ( 1 9 9 8 ) 一文中，他又分析了桩柱群在波浪作用下的冲刷问题，主 要分析了各桩柱不同的排列方式对冲刷的影响。另外小林智尚( 1 9 9 0 ) ( 见陈国平一 文) ，b i y k e r ( 1 9 8 8 ) 等人也做过一定的研究。在国内，陈国平等( 2 0 0 0 ) 分析了波浪 冲刷下的各种因素的影响，得到了最大冲刷深度随波浪的变化规律，并提出了最 大冲刷深度的计算方法。s u m e r ( 1 9 9 7 ) 对波浪作用下圆柱附近的马蹄涡和尾涡运 动做了详细的研究。在同一文章中作者还对造成冲刷的床面剪切力进行了分析。 对波浪作用下局部冲刷的三维数值模拟几乎还没有进行。k o b a y a s h i ( 1 9 9 2 ) 曾经 对振荡流问题进行了尝试，在他的研究中由于k c ( = 5 ) 数小于临界值( k c = 6 ) ， 根据s u m e r 的理论，此时没有马蹄涡和尾涡的发生。此后，i s v a 的j e n s e n 。m a d s e n & f r e d s o e ( 1 9 9 9 ) 和a n d e r s e n ( 1 9 9 9 ) 也分别对类似的问题进行过尝试。 但是直到现在还没有一个数学模型能对整个冲刷过程进行完整的描述，正如 s u m e r ( 1 9 9 9 ) 一文中指出只有对泥沙的输运过程有深入的了解，才能将强的不规 则湍流脉动引入计算中去。我们需要对造成冲刷的涡系有更加深刻的了解，计算 中对涡有更加精确的捕捉，才可能清楚地了解圆柱附近的流场，并通过流场得到 床面的变化情况。 第二秒研，i 方秦反意义 本人论文将主要分为两个部分：圆柱附近流场的计算和底面地形变化的计 算。 通过流场的计算，将对圆柱附近的流场进行具体的分析。在这里将主要考虑 定常情况下的平板圆柱绕流问题，其中底面平板分为光滑和粗糙两种情况。考虑 到问题的对称性，本文将只考虑一半流场的流体运动情况。 对于床面的变化情况，所要考虑的问题属于动床局部冲刷的三维问题。在此， 本文引入一个两维的推移质输沙公式，并对它进行适当的变形，然后用于计算底 面高程变化。 通过对圆柱附近局部冲刷的计算，我们希望能为工程实践提供一定的参考信 息。如果能提早知道结构物附近的冲刷情况，这对人们进行结构物的建设和如何 对底面进行适当的防护将有很大的参考价值。 上嗣；交曩大掌l 一哇n e 文 第二章三维湍流数掌模型 第一节擅铷方程 流体沿固定在床面上圆柱的运动是一个三维流动。它的控制方程包括连续性 方程和动量方程。其中连续性方程为： a p + 塑：o 0 t 0 x j ( 2 1 ) 动量方程( n a v i e r - s t o k e se q u a t i o n ，不可压缩均质并引入广义压力的概念) 为： 警+ 警一瓦0 1 十甜a u ，, 0 x0 x + 剀 口z ， 研 缸， 缸。 ，l 。l，缸，f 上面的方程中x ，( i = 1 , 2 ，3 ) 表示三个方向，u ，( i = 1 , 2 ，3 ) 为三个方向的速度，p 为流 体的密度，p 为广3 ( 玉- j j ( 流体动压力，其大小等于流体绝对压力减去静压力的 分布：p = p 。p o 一露彳) ，f 表示时间，为粘性系数。利用b o u s s i n e s q 假定， 引入涡粘性系数胁，可以得到r e y n o l d 平均的n s 方程为： 警+ 警= 等外，睁剀 s ， 第二节嘲| 嗣渊式 为了解决涡粘性系数，本文将适当的引入湍流模型。考虑不同的湍流模 式：混合长度理论对有流动分离的流场处理时会出现问题；k 一占两方程模式对 有逆压梯度的流场模拟误差也比较的大。w i l c o x ( 1 9 9 3 ) 弓l 入了k 一两方程模型， 在工程实际中发现它对以上两种问题的处理结果更接近于实际。考虑到现在所要 解决的问题，在圆柱前面床面边界层出现分离，在圆柱后也有尾涡的分离，而且 在圆柱表面处存在较强的逆压梯度。综合以上考虑，我们采用w i l c o x 的k 一湍 流模型来处理此问题。 k 一曲湍流模型表示如下： 湍动能k 的方程为： p 瓦o k + 纵i o k = 毛筹一+ 加+ 毒 c “考 c z 其中：湍动能t = 丢“，蚝= 丢( 万+ 矿+ 矿) ，它的单位为：( 长嗲爵间 2 。其中： 其中：湍动能t = 去“，蚝= 寺( “+ v “+ w “) ，它的单位为：i 农哆备问】。其中： 矿所降刳一弘 s ， 上式中一2 p k s , j 一项的出现是考虑到不可压缩流场下_ o u , ：0 ，而湍流切应力中 j蹦 存在“：一p 孑了：一2 p kj 页。( 可见w i l c o x ( 1 9 9 3 ) ) 。 特定湍流耗散率c o 的方程为： p 百o c o + 融虿3 c o2 口詈勺等一参舢2 + 毒l c + 掣r ，虿o c t ) j 。z 固p 百印u ，面钳i 勺蕾一声舢1 瓦i + 掣r 瓦i( 2 6 ) 其中：特定湍流耗散率2 詈，其单位为：间) 。湍流耗散率占= v 等等， ryhjlh j m ， r 其单鱿伥间) 3 ，其中y 为运动粘性系数( 室温下取1 e 柚s ) o 上面两个方程中的各参数的定义如下： = 0 0 7 5 ，= c 。= 0 0 9 ，口= 5 9 ，盯= 0 5 ，盯+ = 0 5 ， 由上面的方程可以得到辅助关系式： ：t 国，：k 2 么(27) 、7 其中f 为表示湍流尺度大小的混合长度。 由此可得涡粘性系数，其单位为：( 长譬篇间) ： 蜥：。生( 2 8 ) 以上给出了此问题的控制方程和计算中采用的湍流模式。通过使用k 一两 方程模式，希望尽量减少流场计算的误差，以便为最后的地形冲刷打下一个良好 上5 竞l 大尊u 曩j e 文 f e r z i g e r ( 1 9 9 6 ) 中指出：t u r b u l e n c e m o d e l ss h o u l db e r e g a r d e d 弱e n g i n e e r i n g a p p r o x i m a t i o n s r a t h e rt h a ns c i e n t i f i cl a w s 。要得到更加精确的流场，我们只能采用 更加精确的模式理论，比如：l a r g ee d d ys i m u l a t i o n ( l e s ) ，d i r e c tn u m e r i c a l s i m u l a t i o n ( d n s ) 。 第三节遗界| 伟 对于整个流场，本文将问题适当地简化，只提四种边界条件：入口边界条件， 出口边界条件，对称边界条件和墙面边界条件。在此，我们不提自由表面边界条 、- 一 件而是采用对称边界来替代。其中的墙面边界包括两种情况：底床的边界条件和 圆柱表面的边界条件。由于压力没有描述它的方程，在此不提它的边界条件，对 于它的处理可以参考后面的数值模拟部分。参考r o u l u n d ( 2 0 0 0 ) ，w i l c o x ( 1 9 9 3 ) ， f e r z i g e r ( 1 9 9 6 ) 我们提出以下的边界条件： 入口：由自己具体指定。r o u l u n d ( 2 0 0 0 ) 是通过计算段无柱的入口段来得到的。 考虑到本文是以均流为流动的最初条件，我们也可以采用类似的方法，首先计算 适当长度的进口段，使其在我们所需要的工作段中，边界层已经属于充分发展的 了，并以此作为入口的边界条件。 出口：r d u l u n d ( 2 0 0 0 ) 采用娑：0 ，其中表示u ，v ，矽，国等不同的物理量。在 f e r z i g e r ( 1 9 9 6 ) 指出，上面的式子是应用于定常流动的，对非定常流动是不满足 的。非定常流动下的边界条件为：娑一u 娑：0 。其中u 与出口的位置无关， a 它的选择使流量守恒，即进口流量= 出口流量。本文流场计算采用定常方式进行， 所取的出口边界条件同r o u l u n d ( 2 0 0 0 ) 的一致。 对称边界条件：在本文处理的问题中将存在三个对称面。采用此方法主要是替代 自由表面边界条件。在对称面的两边无物质的流进和流出，即有：u 。= 0 。考虑 对称性有：掣= 0 ，其中表示u 。，k , c o 等不同的物理量。 o n 墙面边界条件，分两种情况分别讨论： 底床的边界条件：考虑到无滑移条件则有，u i = 0 ，i = 1 , 2 ，3 。对于湍动能七采用 g _ k = 0 ，此种取法对粗糙床面更为精确，而对光滑床面如果采用壁面定律来处理 册 的话，也可以得到同样的结果( f e r z i g e r ( 1 9 9 6 ) ) 。 对于的边界条件w i l c o x ( 1 9 9 3 ) q p 指出：国：& c 。竖。其中： s 月= ( 5 州0 1 0 0 ，+ ， ( 肼 0 缸+ 2 5 2 5 t + 7 0 则表示床面粗糙，如果k s + 5 则表示床面是光滑的，如果i 。+ 处于上述的两个量 之间则说明床面处于过渡区。 本文采用r o l l l l l i l d ( 2 0 0 0 ) 建议的新方法，此两种方法大同小异，只是在对瓯 划分上面有所不同，采用此方法可以使计算网格适当的放大，减少计算的工作量， 但是求解的精度又不会下降。这种修改对粗糙床的流场计算尤其显示其优点。详 细见下式： 耻牌) 3 ，。 2 0 2 不过我们在处理床面磨阻流速的时候，没有采用r o u l u n d 所采用的计算u ，的方 法。在他的论文中作者是利用v a n d r i e s t ( 1 9 5 6 ) 的两维速度分布来求解的。 本文利用u ，的定义：u s = 止万，计算得到u ，的。其中f 为底面的剪切应力， 其大小为： r - ( 肿，埋 ( 2 1 1 ) 其中m 表示离底面最近的控制点的速度的大小，心为此控制点的高度。 圆柱表面边界条件( 考虑圆柱表面为光滑的情况，即颤+ 5 ) ：圆柱表面的无 滑移条件，u i ：o ，f - l ，2 ，3 。对湍动能j 采用婺：o 或j j ：o 均可以。而对国的边 界条件的提法本文参考w i l c o x ( 1 9 9 3 ) d p 所用的方法，有： ：旦 ( 2 1 2 ) 2 a z 2 【2 1 2 ) 其中心为控制点的高度。 以上我们建立了本问题流场部分的数学模型，并提出了相应的边界条件。 9 第置章置簟囊精藏值方镰 第一节一曩背曩 在分析得到了流场运动的控制方程及其相应的边界条件以后，我们来考虑如 何将上面的各个方程进行离散求解。首先对整个问题的背景进行分析，将上一章 所得到的控制方程进行一定的简化。 绕圆柱的三维流动是一个非定常的问题，它的非定常主要体现在圆柱后面的 尾涡的不断脱落。对于均匀来流，其绕圆柱流动的临界s t r o u h a l 数为s t = 0 2 。 对于本文所计算的问题，在对应流动条件下，可以知道此时尾涡脱落的频率约为： 工= 0 9 2 。但是由于水流作用下的冲刷主要同圆柱前端的马蹄涡有关，马蹄涡是 造成冲刷的主要因素。只有在波浪或振荡流作用时，由于来流的不定常，尾涡对 冲刷的贡献才是十分显著的，可以说此时尾涡是造成冲刷的主要因素。有关上面 的论述可以参考s u m e r ( 1 9 9 9 ) 或r o u l u n d ( 2 0 0 0 ) 。我们所考虑的是均匀来流下的 冲刷问题，此时需要提供的漉场信息主要是蹑柱前端的马蹄涡，尾涡的作用暂且 可以忽略，所以将问题定常化，也就是假设： 昙：o( 3 1 ) a 、7 本文计算采用的数值模拟方法是基于有限体积法( f i n i t ev o l u m em e t h o d ) 的 压力修正的s i m p l e 算法。s i m p l e ( s e m i i m p l i c i tm e t h o df o rp r e s s u r e l i n k e d e q u a t i o n s ) 算法是一种求解压力藕合方程组的半隐式方法。有限体积法的基本思 路是将待解的微分方程在每一个控制体上体积分，然后将体积分转化成控制点周 围控制面的面积分，将方程离散为可求解的线性方程组进行求解。其中控制体中 心点的选取方法分：控制点中心法和控制面中心法两种。两种方法各有优缺点， 不过第一种方法使用方便，对网格划分的要求比较的松，所以采用较多，本文即 采用此方法。为了得到变量的变化规律以便为方程的离散及求解打基础，我们需 要知道变量在控制体之间的变化规律。根据不同的精度要求可以得到各种不同计 算方法，比较常用的有l i d s ，c d s ，q u i c k 等方法，其中u d s 是将控制面上 的变量值取为其上游控制点的值，c d s 是用相邻的两个控制点值的线性插值来 得到此值，q u i c k 则需要用到三个控制点的值。采用有限体积法相比于其他的 方法，在较粗的网格下也能显示出较准确的积分守恒。 t 海交葺l 大掌【- b 玲文 为了保证解的物理合理性和满足总体的平衡要求，有限体积法的计算需要遵 循四条基本原则，用以指导插值公式和具体的离散分析。此四条基本原则是( 周 雪漪( 1 9 9 5 1 ) ： a 控制体积交界面的一致性，即当一个表面为相邻的两个控制体积所共有的时 候，在这个控制体积的离散方程中，通过该表面的通量的表达式必须相同。 b 对于离散方程4 中9 = 一中矿+ s 需要满足正系数原则，中心控制点的系 ， 数4 和相关控制点的系数如必须为正值。 c 源项的负坡线性化，既将源项线性化写为曩= s 。+ 邬耳时，系数s ，必须小 于或等于零。 d 相邻节点的系数之和应当等于控制点的系数，也就是说：a ，= a 盯。 ， 对于方程的离散可以用不同的网格形式，包括交错网格和非交错网格，其中 交错网格又称移动网格，是将计算速度分量的网格系统和计算压力分布的网格系 统错开。有关与此的论述可以参考周雪漪( 1 9 9 5 ) 。本文采用非交错网格，将要求 解的所有的未知量放置于同一套网格体系中进行计算，此方法可以较为简便的求 解复杂的几何求解区域。 第二节方j 置的，r ，鼻 上一章我们给出了有关的控制方程，从中可以看到k 方程和c o 方程同动量方 程的形式是相似的。这两个方程的离散也是采用与动量方程相同的方法。这里将 对动量方程的离散做具体的分析，方程和方程的离散类似。在下面的分析中， a 表示标量，j 表示矢量，4 表示数组，其中a 为方程中提到的变量，以后均采 用上面的描述方法。 取控制体如图3 1 ，首先对动量方程( 2 3 ) 定常化，并取u 。= o ，然后对一个 控制点p ，方程( 2 3 ) 两边体积分有( 其中i 代表单位方向矢量i ，7 ，正) ： p ( 露v ) 西d v - - - i v e i d y + 灿+ 脚) v - v q ) d v ( 3 2 ) w 图3 1 控制点p 和与其相关的点的分布 f i g u r e 3 1t h ec o n t r o lp o i n ta n d i t sr e l a t i o n a lp o i n t s e 利用g u a s s 定理及v ( 一) = v p i ，将体积分转化为面积分，可以得到： 辱砷厅厦= 一茁p i 面+ 辱( + 脚) r e 么 ( 3 3 ) 其中厦为同控制点相关的控制面的面积矢量。 首先考虑对流项重力打函，可以得到： 五胂厅d ；i 。莓( 历雨) 7 0 7 。军f 卅7 中7 ( 3 4 ) 其中7 = ( 历厢) 7 = p 厅7 厦7 为面f 上的质量流量。f 表示与控制点p 相关的 控制面，在三维问题中，如图3 1 ，就是东西南北上下六个面，分别用e ，w ，s ，1 1 ，t ，b 表示。0 7 为变量在面f 上的值，厅7 为面f 上的速度矢量，龋7 为控制面f 的面 积矢量。考虑上一节所说的原则四，将上式进行离散，对于其中的一个控制面f 分解为迎风部分和与其相对应的修正量，有： 7 m 7 = ( f 朋中) + ( 巴勋) 7 ( 3 5 ) 7 0 7 = 巴卜+ f 卅，m 矿+ 只，掷7 ( 3 6 ) 其中，+ = m a x ( 。，) ，fs - = m 叫。，瓦，) ，= 华。从上面的表 o i * i u u ；u；d 图3 2 定义修正量时各点间的相互关系 f i g u r e 3 2t h er e l a t i o n s h i pa m o n gp o i n t si nc o r r e c t i o n sd e f i n i t i o n 达式可以知道，对于l 司一个面f ，上面的两个僵之一一定为0 ，l 司时只能有一个 值存在( 此值为流动上游点前的系数，下游点前的系数为零，修正量前的系数就 是上游点前的系数) 。m 为中心控制点的变量值，中旷为同面f 相对应的相邻控 制点的变量值。式中修正量掷7 为： 6 m = 妒( ，) ( 巾“一。) ( 3 7 ) 其中o “中。分别表示上风点和下风点处变量西的值。取： d “中= m ”一o “ ( 3 8 ) d i = “一o ” 取r = 茄，定义系数r ) 为： ) = m a x o ，r a i n ( 2 r ，半+ t l - - k 啦) 】 ( 3 9 ) 其中参数k 2 菊函1 而五，其中碰4 n 仍为外部所给参数，可以取为2 或3 。妒( ，) 取值可以参考图3 3 。从图及上面的式子我们可以得到： f 0rs o 的) ： ! ：。o 。 r ? i + r ：， ( 3 l o ) 烈一。1 警+ 孚r擗 普 由此得到对流项的离散形式如下； 辱肿厅厦2 ；7 中75 莓 ，m ，+ m 峨，- m 矿+ ，勋7 ( 3 1 1 ) 在计算的时候，将修正量作显式处理移到方程的右端，迎风格式部分为隐式留在 方程的左端。 妒( ，) 2 1 址 2 0 图3 39 ( ，) 的取值范围 f i g u r e 3 3t h er a n g ef o r 妒( r ) 下面我们来看耗散项辱( 卢+ 肼) v 颤是如何离散的，将积分运算离散为 辱( + 所) v m 盔2 莩( + 所) 7 ( v m ) - 翻7 ( 3 1 2 ) 式中( v o ) 7 的处理方法可以参考图3 4 ，对于一个其中控制面f ，假设控制面f 上的三个基矢量分别为：羁( i = 1 , 2 ，3 ) ，可以得到： c w ) ，- 嚣焉+ 景碗+ 薏焉 d c ld c ，d c l ( 3 1 3 ) 首先分析上式右端的第一项： 詈确= 箐d 州加) 器d = 竺拈詈季1 ( 3 1 4 ) a 毒1元 1 、 7 翻 。 a g 。 、 图3 4 控制面f 上矢量的定义 f i g t t r e 3 4t h e d e f i n i t i o no f v e c t o r0 1 1f a c ef 1 4 若取季t5 d ，手。2 1 贝u 有： 季1 = 丽8 , 4 ，季1 蚕= l ( 3 1 5 ) 由此，可以得到： ( 砌) ，- 詈舟詈舟詈雪3 ( 3 ，