基于温度异重流模型的连续分层流数值模拟方法研究

牛明昌，丁勇，马卫状，韩盼盼

（1.哈尔滨工程大学船舶工程学院，哈尔滨150001；2.上海交通大学海洋工程国家重点实验室，上海200240）

基于温度异重流模型的连续分层流数值模拟方法研究

牛明昌1，丁勇1，马卫状1，韩盼盼2

（1.哈尔滨工程大学船舶工程学院，哈尔滨150001；2.上海交通大学海洋工程国家重点实验室，上海200240）

文章基于温度异重流模型，通过指定流场中温度的分布及密度与温度的函数关系实现密度的连续分层，并假设流场中的热传导在数值过程中可以忽略不计以确保密度分层的稳定。分别采用RANS方法和LES方法对低傅汝德数和高傅汝德数情况下的圆柱绕流进行数值模拟，其尾流场中的lee波、葫芦状流线以及高傅汝德数情况下的波包与实验值相符，低傅汝德数情况下圆柱上下游速度剖面曲线的最小值与实验值仅存在2%的误差，且其分布规律相同。通过改变流体的导热系数，探讨了Peclet数对该方法收敛稳定性的影响，得到了收敛的数值条件。

分层流；温度异重流；尾流场；大涡模拟

0 引言

由于太阳的热辐射及大洋盐热环流等因素，海洋通常是分层的，其密度、温度和盐度沿垂直方向有明显的变化。当潜体在密度分层流体中运动时，流体质点会受到体积效应下的浮力作用，产生异于均匀环境的水动力学现象，如可以明显观察到流场中源生内波的存在。连续分层作为海洋中典型的密度分层形式受到学者们格外的关注。Chomaz,Bonneton和Hopfinger[1]通过对连续分层流中的球体进行拖曳实验，研究了近场尾流在不同傅汝德数下的形态；王进[2]在连续分层流中对三种不同长径比拖曳潜体激发的内波进行了实验，得到了内波转捩傅汝德数经验公式；Long[3]研究了连续分层流中不同尺度障碍物对定常绕流控制方程求解的影响；Miles[4]在Long模型的基础上用渐近分析的方法研究了不同傅汝德数条件下半圆柱绕流形成的lee波；Stevenson[5-6]通过对连续分层流中的圆柱进行拖曳实验，先后研究了尾流场中的定常内波和流体粘性对lee波的影响，其结果与Lighthill[7]的弥散波模型相吻合；Boyer[8-9]在连续分层流水槽中系统地研究了圆柱绕流的流动现象；丁勇[10]基于多相流混合模型研究了线性分层流中的圆柱绕流现象。目前，国内外对连续分层流的研究主要集中在理论和实验方面，尤其缺少相关的数值模拟研究。

本文通过指定流场中温度的分布及密度与温度的函数关系实现密度的连续分层，并假设数值过程中通过热传导改变的温度值很小，建立了密度连续分层流的数值模拟方法。采用RANS方法和LES方法对其中的圆柱绕流现象进行数值模拟，探讨研究了该方法获得收敛结果的数值条件、收敛结果的稳定性以及尾流场时间历程的演变，将流场特征与Boyer[8]的实验值进行对比，证实了此种方法在数值模拟密度连续分层流体时的可行性。

1 数值模型

1.1 温度场假设

为了与流体动力学控制方程中的对流、扩散项保持一致，本文中的热对流是指由于对流引起的温度物理量输运，与热力学中热对流的意义相异，热传导是指由于扩散引起的的温度物理量输运。

流体的密度主要由温度所决定，为了保证密度垂向的稳定分层，需保证温度在垂向的分布是稳定的，而根据热力学第二定律，经过足够长时间之后，温度在垂向的分布必然趋于均匀。因此本文假设通过热传导改变的温度值在潜体尾迹演化的时间跨度内足够小，即热传导的影响可忽略不计，温度在垂向上的分布短时间内是稳定不变的，如此便可实现稳定的密度分层形式。此时，温度输运方程（能量方程）为对流占优扩散方程，方程如下所示。

其中：T为温度，k为流体的导热系数，c为流体的比热容，Γ=k/c为相对于T的广义扩散系数。方程左边为温度场的局部变化率，右边两项分别为热对流及热扩散项。热传导可以忽略的程度在方程中体现为对流项及扩散项对局部场变化率影响程度的相对大小。

对流与扩散的强度之比可以用Pe（Peclet数）来衡量，其表达式如下所示：

其中：u为特征速度，L为特征长度，α=Γ/ρ为特征扩散系数。Pe数越大，对流占优程度越高，那么必然存在一个临界Peclet数Pecr，当Pe＞Pecr时，假设成立，流场中可以形成稳定的密度分层。

1.2 模拟方法

整个流场的温度、密度分布及计算模型示意图如图1所示。规定该图的原点在圆柱圆心位置，x为流向，z为垂向，H为计算域的高度，U为远方来流速度，g为重力加速度，ρ0和T0分别为流体域中的平均密度和平均温度。为了与实验保持一致，计算域的高度设置为0.2 m，圆柱的直径设置为0.024 m，入口及出口边界足够远。

文中用的无量纲参数：内傅汝德数：Fr=U/Nd；雷诺数：Re=Ud/ν；无量纲时间：t′= tU/d，t在数值模拟条件下为计算时间，无量纲时间用tn′表示，在实验条件下为绕流时间，无量纲时间用te′表示；N为浮频率，N=（g△ρ/ρ0H）1/2；ν为流体的运动粘性系数，其它参数的含义如图1所示。

图1 计算模型示意图Fig.1 Sketch of the calculation model

温度T的垂向分布形式如（3）式所示。当T=300 K时，ρ=1 018.58 kg/m3；当T= 300.2 K时，ρ=998 kg/m3。由此可得整个流域温度的垂向分布函数如（4）式所示。

此时，浮频率的值为1，与实验相符。值得注意的是，流体的运动粘性系数ν虽为流体的固有属性，但是受温度的影响很大。根据第10届ITTC关于淡水和盐水粘性系数的规定，在T=300 K附近，如果温度变化5 K，粘性系数相差约0.07×10-6m2·s-1，本文流场的最大温差为0.2 K，粘性系数在该范围内的变化非常小，因此可以将粘性系数视作常数。

Boyer的实验表明，当Fr=0.4时，lee波已经不是尾流场的主要特征，这时尾流场中出现明显的涡结构，所以规定当Fr＞0.4时为高傅汝德数，Fr＜0.4时为低傅汝德数。对低傅汝德数下的情形，可通过RANS方法予以模拟实现，对高傅汝德数下的情形，则应通过LES方法实现。在高傅汝德数条件下，尾流场中湍流脉动成份的水动力作用会增强，而RANS方法抹去了瞬时脉动成份，LES方法则没有这种弊端。理论研究[11-13]认为，当横向距离大于πd时，利用大涡模拟进行计算，三维的影响就变得非常小，因此本文横向设置为0.084 m（3.5d）。

采用O型网格对流体域进行空间离散，二维及三维情况下的网格数量分别为30万和150万。对于RANS方法，压力离散格式采用“PRESTO！”，对流项离散格式采用“Second Order Upwind”，扩散项离散格式采用中心差分格式，时间离散格式采用“First Order Upwind”。对于LES方法，压力离散格式采用“Standard”，对流项离散格式采用“Bounded Central Differencing”，扩散项采用中心差分格式，时间离散格式采用“Second Order Implicit”。同时两种方法都是采用PISO算法进行求解，时间步长为0.005 s。

2 数值结果与讨论

2.1 低傅汝德数条件下的数值结果

在低傅汝德数条件下，计算了三组工况，具体的参数如表1所示。

工况1条件下的流线图如图2所示，图中（a）为数值结果，（b）为实验结果。该图表明低速条件下虽然在圆柱的下游形成一对与实验值较为相符的lee波，但是下游的流线形状并不能与实验值完全相符。说明在低速条件下对流项的影响较弱，此时热扩散并不能完全忽略不计。

工况2条件下结果稳定时的流线图及相同实验条件下的流线图如图3所示，图中（a）为数值结果，（b）为实验结果。（a）图中可以观察到在圆柱的后方有两对较明显的lee波峰线，其特征与（b）图的实验结果相吻合。图（a）、（b）在靠近x轴线处lee波的形状与大小基本一致，不同的是数值结果中x轴线附近的流线上下两侧间距较大，实验条件下x轴线附近上下两侧的流线在紧邻圆柱及圆柱后方第一个波包后端几近接触。

表1 工况表（低傅汝德数）Tab.1 Conditions of calculation (low Froude number)

图2 流线图对比（Fr=0.018，Re=12）Fig.2 Comparison of streamline(Fr=0.018，Re=12)

图3 流线图对比（Fr=0.08，Re=60）Fig.3 Comparison of streamline(Fr=0.08，Re=60)

工况3条件下的流线图及实验条件下的流线图如图4所示，图中（a）为数值结果，（b）为实验结果。图（a）表明紧邻圆柱的尾流场流线呈葫芦状，与图（b）中Boyer实验得到的流线图非常吻合，这是该阶段尾流场的主要特征。另外，图（a）中还存在斜向传播的lee波，这种斜向传播的lee波在实验结果中也可以观察到，只是lee波的位置略有差异，数值模拟过程中没有观察到转子的出现。在该工况下，数值结果与实验结果符合较好。

图4 流线图对比（Fr=0.17，Re=98.7）Fig.4 Comparison of streamline(Fr=0.17，Re=98.7)

2.2 高傅汝德数条件下的数值结果

在高傅汝德数条件下，计算了两组工况，具体的参数如表2所示。

图5为Fr=0.88，Re=480条件下流线图的实验结果及数值结果，图（a）为实验结果，图（b）为数值结果。该图表明数值和实验条件下的尾流场流线形状大致相同，在紧邻圆柱的后方都有一对涡出现，在距离圆柱6倍直径和16倍直径的尾流场中存在着两个波包，第一个波包下包络着一对涡，图（b）中第二个波包所在的位置比图（a）中第二个波包出现的位置距离圆柱近，这是由于波包是一个不稳定的尾流场特征，它在水平及垂向位置上不断震荡，震荡幅度可以达到5倍的圆柱直径。该数值方法基本模拟出了实验工况下的结果。

表2 工况表（高傅汝德数）Tab.2 Conditions of calculation (high Froude number)

图5 流线图对比（Fr=0.88，Re=480）Fig.5 Comparison of streamline(Fr=0.88，Re=480)

图6 涡量图对比（Fr=1.77，Re=960）Fig.6 Comparison of vortex(Fr=1.77，Re=960)

表3 不同导热系数下的工况表Tab.3 Conditions of calculation at different heat conductivity

图7为Fr=0.018，Re=12时不同Pe数条件下圆柱绕流流线图，图（a）为Pe=7 212的情形，图（b）为Pe=72 120的情形，两幅图中的流线明显不同于Pe=72.12时的情形（图2（a））。当Pe=7 212及Pe=72 120时，尾流场流线形状与实验结果更加相符（图2（b））。图8为Fr=0.08，Re=60时不同Pe数条件下圆柱绕流流线图，图（a）为Pe=32 054的情形，图（b）为Pe=320 540的情形，两幅图也不同于Pe= 320.54的情形（图3（a）），当Pe＞32 054时，尾流场流线形状与实验结果相符（图3（b））。图9为Fr= 0.17，Re=98.7时不同Pe数条件下圆柱绕流流线图，图（a）为Pe=68 114的情形，图（b）为Pe=681 140的情形，结果与图7、图8的规律相符。但是当傅汝德数增大到一定程度时，Pe数的值已经很大了，改变流体导热系数对流场的影响变得非常有限。由于连续分层流中的圆柱绕流是一个准稳态过程，流线随时间变化，图7～9各图中（a）与（b）数值时刻均相同。

图7 流线图（Fr=0.018，Re=12）Fig.7 Diagram of streamline(Fr=0.018，Re=12)

图8 流线图（Fr=0.08，Re=60）Fig.8 Diagram of streamline(Fr=0.08，Re=60)

图9 流线图（Fr=0.17，Re=98.7）Fig.9 Diagram of streamline(Fr=0.17，Re=98.7)

分析表明，在同样导热系数条件下，随着来流速度的增加，Pe数越来越大，对流占优程度越高。虽然低傅汝德数工况下结果与实验值不完全相符，但在高傅汝德数工况下能够取得与实验一致的结果。因此如果在低速条件下的数值结果能够与实验值吻合，那么高速条件下该方法理论上也能够适用。

2.4 尾流场时间历程的研究

在模拟过程中，数值时间并不能与实验时间完全相对应，但在某数值时刻，若上游速度剖面曲线与某一实验时刻结果相对应，同时刻下游速度剖面数值结果也应与同实验时刻的结果相对应。图10～13为Fr=0.018，Re=12时数值及实验条件下上下游（x=±7.5d）速度剖面对比图，图10、图11为Pe= 7 212时的情形，图12、图13为Pe=72 120时的情形。速度曲线的极值代表该剖面处速度的最大值或最小值，在上游速度剖面曲线吻合的时刻，下游速度剖面基本吻合，四组结果中波谷值的相对误差最大为2%。另外，在Pe=7 212及Pe=72 120两种情形下，实验时刻te′=62对应的数值时刻均是实验时刻te′=21对应数值时刻的3倍左右（实验时刻对应的倍数也为3），在时间演化上，数值结果与实验结果相符。

图10 上下游速度剖面（Pe=7 212，tn′=7.88，te′=21)Fig.10 Velocity profiles of upstream and downstream(Pe=7 212,tn′=7.88,te′=21)

图11 上下游速度剖面（Pe=7 212，tn′=24.06，te′=62)Fig.11 Velocity profiles of upstream and downstream(Pe=7 212,tn′=24.06,te′=62)

图12 上下游速度剖面（Pe=72 120，tn′=6.64，te′=21)Fig.12 Velocity profiles of upstream and downstream(Pe=72 120,tn′=6.64,te′=21)

图13 上下游速度剖面（Pe=72 120，tn′=18.9，te′=62)Fig.13 Velocity profiles of upstream and downstream(Pe=72 120,tn′=18.9,te′=62)

3 结论

基于温度异重流模型，通过指定流场中温度的分布及密度与温度的函数关系实现了密度的连续分层，并假设流场中的热传导在数值过程中忽略不计以确保连续分层的稳定性。分别采用RANS方法和LES方法数值模拟了低傅汝德数和高傅汝德数工况下的圆柱绕流，其尾流场特征与实验相符合，证实了这种方法在数值模拟密度连续分层流时的可行性。具体结论如下：

（1）随着傅汝德数的增大，Peclet数逐渐变大，能量方程变为对流占优扩散方程，流场中的热传导可以忽略不计，密度分层能在较长时间内保持稳定，因此这种方法在数值模拟连续分层流时的稳定性和精度都比较高。

（2）采用RANS方法成功模拟出了低傅汝德数条件下尾流场中的lee波，紧邻圆柱的尾流场流线呈现葫芦状，其特征与实验值相符。圆柱上下游（x=±7.5d）速度剖面曲线出现了一个极小值和两个相等的极大值，其变化规律与实验定量一致，峰值的误差在2%以内，且在时间历程演变方面与实验相符合。

（3）采用LES方法对连续分层流中高傅汝德数条件下的圆柱绕流进行了数值模拟。当Fr=0.88时，距离圆柱6倍直径和16倍直径的尾流场中分别出现了一个波包，且其在水平及垂向位置上不断发生震荡。当Fr=1.77时，尾流场进入了完全湍流状态，圆柱后方出现了很强的涡结构，由于分层的抑制作用，涡街在垂向的宽度保持不变。

（4）RANS方法能够捕捉到流场中流线的变化，但对涡量的描述不够细致，LES方法则弥补了这一弊端，能够精确模拟出流场中的各种成分。因此，RANS方法适合于低傅汝德数下的流动，LES方法则同时适合低傅汝德数和高傅汝德数下的流动。

[1]Chomaz J M,Bonneton P,Hopfinger E J.The structure of the near wake of a sphere moving horizontally in a stratified fluid[J].Journal of Fluid Mechanics,1993,254:1-21.

[2]王进,尤云祥,胡天群,等.密度分层流体中不同长径比拖曳潜体激发内波特性试验[J].科学通报,2012,08(57): 606-617. Wang Jin,You Yunxiang,Hu Tianqun,et al.The characteristics of internal waves excited by towed bodies with different aspect ratios in a stratified fluid[J].Chin Sci Bull,2012,08(57):606-617.

[3]Long R R.Some aspects of the flow of stratified fluids[J].Tellus,1955,7(3):341-357.

[4]Miles J W H.Lee waves in a stratified flow.Part 2.Semi-circular obstacle[J].Journal of Fluid Mechanics,1968,33(4): 803-814.

[5]Stevenson T N.Some two-dimensional internal waves in a stratified fluid[J].Journal of Fluid Mechanics,1968,33(04): 715-720.

[6]Stevenson T N,Chang W L,Laws P.Viscous effects in lee waves[J].Geophysical&Astrophysical Fluid Dynamics,1979, 13(1):141-151.

[7]Lighthill M J.On waves generated in dispersive systems to travelling forcing effects,with applications to the dynamics of rotating fluids[M].Hyperbolic Equations and Waves,Springer,1970:124-152.

[8]Boyer D L,Davies P A,Fernando H,et al.Linearly stratified flow past a horizontal circular cylinder[J].Philosophical Transactions of the Royal Society of London A:Mathematical,Physical and Engineering Sciences,1989,328(1601):501-528.

[9]Xu Y,Fernando H J,Boyer D L.Turbulent wakes of stratified flow past a cylinder[J].Physics of Fluids(1994-present), 1995,7(9):2243-2255.

[10]丁勇,韩盼盼,段菲,马卫状.线性分层流中圆柱绕流数值模拟方法研究[J].哈尔滨工程大学学报,2016,9:1-5. Ding Yong,Han Panpan,Duan Fei,Ma Weizhuang.Numerical study of linearly stratified flow past a cylinder based on multiphase mixture model[J].Journal of Harbin Engineering University,2016,9:1-5.

[11]M B.Large eddy simulation of the subcritical flow past a circular cylinder:numerical and modeling aspects[J].International Journal for Numerical Methods in Fluids,1998,28(9):1281-1302.

[12]M B.Numerical and modeling influences on large eddy simulations for the flow past a circular cylinder[J].International Journal of Heat And Fluid Flow,1998,19(5):512-521.

[13]Kravchenko A G M P.Numerical studies of flow over a circular cylinder at ReD=3 900[J].Physics of Fluids,2000,12(2): 403-417.

Research on the numerical simulation methods of continuously stratified flows based on thermal density current model

NIU Ming-chang1,DING Yong1,Ma Wei-zhuang1,Han Pan-pan2
(1.College of Shipbuilding Engineering,Harbin Engineering University,Harbin 150001,China;2.State Key Laboratory of Ocean Engineering,Shanghai Jiao Tong University,Shanghai 200240,China)

Based on thermal density current model,the continuous stratification of density is achieved by specifying the distribution of temperature and regarding density as a function of temperature.It is assumed that the heat conduction is negligible to ensure the stability of density stratification.The RANS method and LES method are used respectively to simulate the flows past a cylinder at low Froude numbers and high Froude numbers.The lee waves,gourd-shaped streamlines and wave packets are consistent with the experimental results.For the minimum values of velocity profile curves,there are only 2%errors between the calculated and experimental values,and the curve shapes are the same.The influences of Peclet numbers on the stability of this method are discussed by changing the heat conductivity coefficient of fluid.Finally,the numerical conditions of convergence are obtained.

stratified fluid;thermal density current;wake field;large eddy simulation

U661.1

A

10.3969/j.issn.1007-7294.2017.08.002

1007-7294（2017）08-0941-09

2017-03-26

××减震降噪工程专项计划

牛明昌（1993-），男，硕士研究生；丁勇（1959-），男，博士，教授，通讯作者，E-mail:dingyong@hrbeu.edu.cn。