大长高比腔体中强非线性状态混合流体渐进波对流图案的数值模拟

大长高比腔体中强非线性状态混合流体渐进波对流图案的数值模拟

1混合流体对流时空结构的研究在封闭的腔中，上下表面的温度固定，上下表面的加热，因此形成的温度差超过一定临界值，腔中的液体运动会导致腔体运动。对于这种流动现象,Benard等人早在1900年就进行了实验研究,Rayleigh最早(1916年)进行了理论分析,从而,也被称作Rayleigh-Benard对流。它是研究非平衡对流的稳定性、时空结构或图案(Pattern)、非线性动力学特性的典型模型之一。对这种流动现象的研究,同时可帮助人们对海洋,水库,湖泊中温度驱动的对流现象的理解。对于纯流体Rayleigh-Benard对流问题的研究,已经有了很长的历史,也获得了大量的成果。上世纪80年代以来,水和其他溶液混合形成的混合流体的热对流运动开始引起了研究者的重视。当表征流体非线性特性的参数——分离比(Separationratioψ)大于0时,对流系统出现了类似于纯流体对流(ψ=0)时的分叉特性,但其临界值小于纯流体对流时的临界值,流动图案全是定常的。对于ψ<0的混合流体对流,当温度梯度超过某个临界值之后,系统会出现对流滚动向某个方向传播的行进波对流图案,这极大的吸引了研究者的兴趣。引起了实验工作者的重视。利用观测一维对流时空结构的实验装置,研究人员在观测均匀行进波对流的时空结构的同时也获得了各种行进波的对流时空结构。可是,许多研究都是在中长高比腔体中进行的,对于大长高比腔体中行进波的对流时空结构研究的比较少。另一方面,利用流体力学方程,热传导方程及物质传输方程可以精确的描述混合流体Rayleigh-Benard对流运动。基于流体力学方程组的数值摸拟可以提供在实验中很难观测的流场,温度场及浓度场的结构特性。因此,流体力学方程组的数值摸拟将为理解对流的微细结构提供有力的工具。Yahata利用流体力学扰动方程组、比较了MAC法与迦辽金法的计算结果,并讨论了中长高比腔体中行进波与局部行进波的特性。Barten等基于流体力学基本方程组计算了周期边界条件下均匀行进波的场的结构,局部行进波的形成过程及相应的动力学特性。所有这些模拟,都提供了比实验更详实的资料和信息。通过分析一维对流时空结构(二维流动摸拟)及其动力特性,解释了一些非线性行进波对流现象。但摸拟研究才刚刚开始,其成果仍是初步的[5,6,7,8,9,10,11,12]。为了理解许多行进波对流现象及对流滚动的形成机理,通过流体力学方程组的数值摸拟系统的研究大长高比腔体中对流时空结构是很有必要的。2数学模型2.1soret效应的动力学方程在Boussinesq假设下,密度除了在浮力项中表示为温度和浓度变化的线性组合外,在流体力学方程的其它项都取为常量。对于温度T、浓度C离开参考值T0、C0一个小的偏差,质量密度的状态方程可表示为ρ=ρ0[1-α(Τ-Τ0)-β(C-C0)]ρ=ρ0[1−α(T−T0)−β(C−C0)]式中热引起的体积膨胀系数α和浓度引起的体积膨胀系数β分别为α=-1ρ0∂ρ(Τ,Ρ0‚C0)∂Τα=−1ρ0∂ρ(T,P0‚C0)∂Tβ=-1ρ0∂ρ(Τ0,Ρ0‚C)∂Cβ=−1ρ0∂ρ(T0,P0‚C)∂Cα的作用是单调的,因为不论液体或气体,温度升高时密度都会减小,α始终大于0。而β的作用却不一样。当β<0时,对流系统类似于纯流体的情况为定常流动。可是当β>0时,对流系统会出现行进波状态,局部行进波及定常流动。下标0表示传导状态的平均值,一般取为腔体高度二分之一处的数值作为计算值。对于二维腔体内的混合流体对流运动,设水平坐标轴x与腔体的下壁面重合,向右为正,垂直坐标轴z与腔体的左壁面重合,向上为正。如果所有的长度由d、时间由d2/v、速度场由v/d、温度场由kv/αgd3、浓度场由kvγ2/αgd3D、压强由v2/d3无因次化,则考虑了Soret效应的流体力学基本方程组可表示为:∇⋅U*=0(1)∂U*∂t+(U*⋅∇)U*=∇2U*-∇(pρ0)*+1Ρr[(Τ*-Τ*0)(1+Ψ)-Ψ(ζ*-ζ*0]ez(2)∂Τ*∂t+(U*⋅∇)Τ*=1Ρr∇2Τ*(3)∂ζ*∂t+(U*⋅∇)ζ*=LΡr∇2ζ*+1Ρr∇2Τ*(4)式中Ψ=-βα⋅γ2D为分离比,Pr=v/k为普朗特数,L=D/k为落依斯数,ζ*=C*+T*,γ2=STC0(1-C0)D;其中,U*,T*,ζ*,d,v,k,g,D,ST,ez分别为流速矢量场,温度,浓度流束,容器高度,运动黏性系数,热导率,重力加速度,浓度扩散系数,Soret系数,z方向的单位矢量。2.2根据不同的kv值的瑞利数计算表1对于无流动的传导状态,由方程组(1)—(4)可得传导解为:U*cond=0(5)Τ*cond=Τ*0+R(0.5-z)(6)C*cond=C*0+R(z-0.5)(7)ζ*cond=ζ*0(8)(pρ0)*cond=(pρ0)*0-12RΡr(1+Ψ)(12-z)2(9)式中R=(αgd3/kv)ΔT为瑞利数,ΔT为上下板间的温度差。2.3d的扰动方程进一步定义扰动量为:u=U*-U*cond,θ=Τ*-Τ*cond,ξ=C*-C*cond,η=ζ*-ζ*cond由方程组(1)—(4)和传导解(5)—(9),我们可得扰动方程为∇⋅u=0(10)∂u∂t+(u⋅∇)u=∇2u-∇(Ρρ0)+1Ρr[(1+Ψ)θ-Ψη]ez(11)∂θ∂t+(u⋅∇)θ-Rw=1Ρr∇2θ(12)∂η∂t+(u⋅∇)η=LΡr∇2η+1Ρr∇2θ(13)式中u(u,0,w),w分别为流速矢量场和垂直方向的流速。2.4流束的边界条件为了求解控制方程,需给出θ,η,u,w的边界条件,在z=0,1上,我们给出实际的无滑移、不可渗透的边界条件,在无滑移边界(z=0,1)上的浓度流束可表示为:∂η∂z=∂∂z(ζ-ζcond)=0在z=0,1上速度的边界条件为∂w∂z=w=u=0由于在上下边界上温度是固定的,温度的扰动量可表示为θ=0对于限定在矩形容器中的混合流体,侧向边界条件可表示为u=w=∂u∂x=0∂θ∂x=∂η∂x=02.5压力方程的求解计算域为1×46,计算中采用均匀网格Δx=Δz=116。为了求解偏微分方程,使用了MAC法,空间微分采用了中心差分法,时间微分采用向前差分,时间步长为Δt=0.0005(d2/v),由连续方程和动量方程推导的压力方程具有波松方程的形式,其边界条件由速度、温度及浓度的边界条件直接推导给出。利用ICCG法求解了压力方程。计算中采用的流体参数为,Pr=13.8,L=0.01,Ψ=-0.47。3数值匹配结果表明3.1进波场的非线性状态计算在相对瑞利数r=2.1(r=R/Rco,Rco=1708)处从很小的高斯分布初值开始,随时间的增加,我们成功的观测到了线性行进波,线性行进波向非线性行进波的过渡及非线性行进波状态。图1是在初始扰动作用下温度场沿腔体高度二分之一处分布的成长过程。扰动从中心向两边传播,在这一阶段行进波具有线性性质,行进波的最大振幅的成长率随时间线性增长,如图2所示。当浓度场沿腔体高度二分之一处的分布由调和函数或光滑曲线变成台形结构时,行进波场由线性阶段进入非线性结构。当达到对流振幅饱和的非线性行进波状态时,行进波速度减慢最终保持为常数,图3给出了计算区间36~46之间行进波场的等值线图。由图3中可以发现行进波的波数保持为K=3.548,波长为1.769。由上至下,图3(a)表示垂直速度分量等值线图,图3(b)表示流线的等值线图,可以看出即使在大振幅的强非线性状态,等值线图呈方圆状。图3(c)是扰动温度,其等值线图呈现三角形和倒三角形分布状。总的温度分布是不封闭的弯曲曲线如图3(d)。图3(e)的扰动浓度呈明显的台状分布。图3(f)的总的浓度场的分布,在方向相反的滚动内都具有浓度羽。最下边的图3(g)是shadowgraph强度分布,显示在每个滚动内都具有羽状分布,但相邻的滚动的浓度羽具有相反的方向。由于行进波处于稳定状态,因而,对流滚动的形状分布都是均匀的。3.2易发生区域的变化局部行进波是指在均匀的上下温度差作用下,计算域内部分区域存在对流行进波,部分区域是传导状态。这种现象称之为局部对流行进波。当我们从r=2.1缓慢的减小r,然后取r=1.319,在进行连续计算时,对流振幅经过一段过渡后,很快减小到0,说明对流已恢复到了传导状态。当从r=2.1减小到r=1.35时,均匀行进波的defect开始扩大,首先变成缺损区域,然后缺损区域不断扩大,最终变成稳定的缺损区间即传导区间,行进波区域保持着稳定的常长度。图4的斜线部分说明了行进波的前进方向和行进波所占空间长度,中间部分是传导状态,可以看出在t=100(d2/v)内局部行进波对流是稳定的。当从r=1.35的局部行进波状态,增加r到r=1.5时,经过长时间的计算,我们获得了稳定的局部行进波状态。图5是计算域高度二分之一处的对流行进波的分布。可以看出由上至下图5(a)的垂直速度,图5(b)的扰动温度与图5(d)的流线及图5(e)的局部Nusselt数均呈光滑的曲线状态,图5(c)的浓度场具有明显的台形分布,最下边图5(f)的shadowgraph强度是沿高度的平均值,其分布是与实验类似的。在两端对流行进波的存在域以外,中间部分是传导状态,行进波的对流区域长度较r=1.35时有所增加。当r进一步增加到r=1.57时,获得了类似于图5的分布曲线,如图6所示。这里行进波区域所控制的长度明显增加,传导区域明显减小。总之在Saddlenode点附近r=1.35~r=1.57之间,我们获得了稳定的局部行进波,并且行进波的对流区域长度随r的增加而增加。3.3局部并行波状态实验是在长高比为Γx=46,宽高比为Γx=0.6的腔体内进行的。实验采用的流体参数与计算时相等。利用阴影法观测了图案。对流实验中,当r=6.8时,对流发生。在对流临界点附近观察到的行进波图案如图7所示,可以看出CPW图案与计算结果是一致的。然后经过长时间的对流振幅成长和defect出现而伴随的波数调整,过渡到均匀行进波状态。随着r的减小,当r=5.5时,系统出现局部行进波状态,空腔的中间部分是传导状态,两端附近为对流状态。这种局部行进波状态在r=4.0~5.5的范围内是稳定的。左边的对流区域长度L1和右边的对流区域长度L2随着相对瑞利数r的变化如图8,9所示。可以看出,随着相对瑞利数r的增加,对流区域长度都在增加。这和我们的计算结果是一致的。说明对于相同流体参数的实验和计算结果其现象是一样的。需要指出的是,在实验中对流发生的临界瑞利数r,Saddlenode点附近存在着局部行进波的瑞利数范围r都比计算时的数值大。或者说实验获得的分叉曲线各控制点的数值都相对计算结果比较大。这是因为计算是二维的,实验中为了保证一维对流现象出现,腔体的宽高比选为Γx=0.6,因而,两侧壁面明显的影响着对流运动。另外,Barten的计算也说明二维的计算结果和Γx=1.3