应用srgk模型对六直叶涡轮搅拌桨流动场的数值模拟

应用srgk模型对六直叶涡轮搅拌桨流动场的数值模拟

搅拌槽内流动场的模拟搅拌槽广泛应用于各种工艺操作中，如液体混合、气体分散、固体悬浮和反应处理。近年来,利用计算流体力学(CFD)方法研究搅拌槽内的流动场已经成为搅拌反应器研究的新热点。对搅拌槽内流动场的数值研究多数采用标准k-ε模型。由于该模型基于各向同性假设,对搅拌槽内流动场,尤其是搅拌桨附近流动场的模拟并不合适。实验表明,该区域的流动非常复杂,呈现强烈各向异性。为了改善对此区域的预报,许多学者进行了各方面的努力,如采用高精度差分格式,或是对k-ε模型本身的改进和修正。其中比较有代表性的是由Yokhot和Orszag等人应用重整化群理论发展并改进的RNGk-ε模型,其基本思想是把湍流视为受随机力驱动的输运过程,通过频谱分析的方法消去其中的小尺度涡并将其影响归并到涡粘性中,以得到所需尺度上的输运方程。在高雷诺数时RNGk-ε模型与标准的k-ε模型具有相同的形式,只不过在方程中出现一个附加生成项,这一项随流动的畸变程度而变化,从而会在一定程度上补偿标准k-ε模型的某些不足。本文应用该模型对六直叶涡轮搅拌桨的三维流动场进行了数值模拟,通过对搅拌桨附近流动场的研究,分析了该模型计算结果的一些特点。1湍流模型及应力模型对变量ϕ的守恒方程可以用(1)式表示∂∂t(ρ⋅ϕ)+∂∂xi(ρ⋅ui⋅ϕ)=∂∂xi(Γϕ⋅∂ϕ∂xi)+SϕⅠⅡⅢⅣ(1)∂∂t(ρ⋅ϕ)+∂∂xi(ρ⋅ui⋅ϕ)=∂∂xi(Γϕ⋅∂ϕ∂xi)+SϕⅠⅡⅢⅣ(1)式中(Ⅰ)为时间项,(Ⅱ)为对流项,(Ⅲ)为扩散项,(Ⅳ)为源项。上式对应于质量,动量,湍流动能和耗散率的守恒方程时,ϕ分别为1,uj,k和ε。其展开形式如表1所示。有关模型常数的数值如表2所示。同标准k-ε模型一样,RNGk-ε模型也采用基于Boussinesq假设的雷诺应力关联式。在高雷诺数时有μΤ=Cμρk2εμT=Cμρk2ε。因此,在形式上它与标准k-ε模型几乎完全相同,区别仅在于:标准k-ε模型中的系数均由实验确定,而RNGk-ε模型中的系数则由理论推导产生;在标准k-ε模型中,系数C*1是一个简单的常数,而在RNGk-ε模型中,C*1则是应变率的函数,由此间接改进了对耗散率方程的模拟,在一定程度上考虑了湍流的各向异性效应,从而能够改善对复杂湍流的预报结果。2计算策略2.1搅拌槽直径计算计算所采用的搅拌槽结构与文献的实验装置一致。槽体为圆柱形,均布4块挡板。搅拌桨为六直叶涡轮。搅拌槽直径T=0.27m,液位高H=T,挡板宽为T/10。搅拌桨直径D=T/3,桨叶离底距离C=T/3。工作介质为水。计算中搅拌转速为200r/min。在此条件下,叶端线速度为Utip=0.94m/s,搅拌雷诺数Re=2.7×104。根据流动的对称性,计算域选取了槽体的一半。2.2算过程中的收敛性计算中采用的网格是结构化的六面体网格,这种网格的划分比较复杂,但在计算过程中的收敛性比较好。计算中采用的网格是:66×36×77(r×θ×z),共有172090个网格。为了能够更准确地预报桨叶附近的流动场,对此区域内的网格进行了加密,在叶片表面采用了2×18的网格分布。2.3外部静电网格的时均计算计算使用的软件是CFX4.3。计算方法采用滑移网格法。这种方法将计算域分为两部分,如图1所示。一部分包含了运动的桨叶(a),另一部分包含静止的挡板与槽体(b)。两部分网格之间要求彼此独立。在CFX4.3中是通过定义非匹配边界条件来实现。外部的网格静止,内部的网格随搅拌桨一起转动,如图2所示。两部分网格之间通过滑移界面进行插值处理。由于滑移网格法计算的是瞬时值,为了能将计算结果与稳态的实验数据作比较,需要对计算结果进行处理。对外部静止网格,可以作简单的时均处理。时间平均所选取的时间间隔是桨叶转动一圈所需时间的1/6,若它对应Nθ个时间步,则变量Φ可用(2)式进行平均ΦL=1ΝθΝ0+Νθ∑n=Ν0+1ΦnL(2)式中,ΦnL是变量Φ在第n个时间步和切向位置L时的值;N0是一个足够大的时间步,以确保初始的瞬态效应消失。对旋转部分的处理有许多方法,文献有详细介绍。求解方程采用控制容积法。利用Rhie-Chow算法来避免发生振荡,压力-速度的耦合求解采用SIMPLEC算法,对流项的离散使用混合-上风差分格式。由于滑移网格法计算量非常大,为节省时间,在计算开始时选取一个较大的时间步进行计算,以消除初始效应。在最后一轮计算时,选用小的时间步进行计算,同时利用自己开发的用户接口程序对计算结果进行时均转换。3数值模拟结果与分析3.1叶轮区流体流动特征图3是两挡板中间平面上的宏观流动场。六直叶涡轮属径向流搅拌器,旋转的叶轮产生高速径向射流,撞击到槽壁后,流体出现分离,一部分沿槽壁向上流动,一部分向下流向槽底。到达液面或槽底后,流体发生转向,沿轴向又流回叶轮区。从图中流体的整体流动形式可以看出,计算结果与许多实验研究的报道一致,成功再现了涡轮搅拌桨流场的典型流动特征。图4是不同高度水平截面的宏观流动场。由于靠近液面与槽底处的速度与桨叶中心所在平面的速度相差较大,为了能够清楚看到上下平面的速度分布,在此两平面上的矢量图使用了更大的倍数来放大。从图中可以看出,在桨叶中心所在平面流体是向槽壁处流动的,其径向与切向速度均较大。在靠近液面与槽底处的两个平面上,流体从外部向槽中心方向流动。在液面处靠近挡板的地方可以观察到旋涡的存在,与在实验中观察到的现象一致。在靠近槽底平面的中心位置,有一个较大的旋涡存在,流体的流动比较复杂。3.2测量结果及比较图5、6为不同湍流模型的计算结果与实验数据的比较。为了进一步考察RNGk-ε模型的特点,并便于比较,将标准k-ε模型的计算结果也加了进来。图5为径向速度分布的比较,对于径向速度分布,计算结果比起同位置的实验测量结果略小。对于最大径向速度,由于测试方法的差别,不同研究者的研究结果有一定的差距。Wu与Dyster的实验结果一致,大约是0.75Utip。而Kemoun的结果则与Mahouast的结果比较接近,大约是0.53Utip。侯栓弟对激光多普勒测量结果进行了统计得到,最大径向速度的平均值为0.7Utip。由图5可见,标准k-ε模型与RNGk-ε模型对径向速度的预报比较接近,差别主要是在最大径向速度上。采用标准k-ε模型时,最大径向速度0.63Utip;而RNGk-ε模型的计算结果是0.68Utip,比标准k-ε模型更加接近实验结果。图6为切向速度分布的比较。对切向速度分布的预报与实验数据的比较要比对径向速度分布的预报更加接近。由图6可见,计算结果与实验测量结果吻合较好。对于最大切向速度的测量,不同研究者之间的差别也相对较小。Kemoun的测量结果为0.71Utip,Mahouast与Wu的测量结果接近,约为0.68Utip。标准k-ε模型与RNGk-ε模型对切向速度的预报基本一致。采用标准k-ε模型时,最大切向速度为0.66Utip;而RNGk-ε模型的计算结果为0.68Utip,比标准k-ε模型略有改善。由此可以看出,RNGk-ε模型在耗散率方程中通过系数C*1引入描述流场畸变效应的附加源项后,确实改善了对搅拌桨附近复杂流动场速度分布的预报,表明RNGk-ε模型在速度分布预报方面要优于标准k-ε模型。3.3实验结果和讨论搅拌槽内湍流动能的分布对许多混合过程来说是很重要的。但是,利用CFD方法所获得的湍流动能分布却远远没有对速度分布的预报那样准确,尤其是对桨叶附近湍流动能较高的区域,所得到的结果均较实验值偏小。Schafer对涡轮桨的测量结果是k/U2tip=0.158。Ng用滑移网格法计算了直径0.1m搅拌槽内涡轮搅拌桨排出流区的湍流动能,计算值比实验值小。增加网格密度后,预报精度有所提高。但是,即使采用了很高的网格密度(约40万),计算值也大约仅有实验值的50%。本文计算结果如图7所示,由图7中可见,RNGk-ε模型对湍流动能的计算值比实验结果要小得多,而且比标准k-ε模型的计算结果也要小。这一方面是由于k-ε模型本身的缺陷造成了对湍流动能的预报偏低,另一方面可能是由于系数C*1的引入抑制了湍流动能的产生,造成了湍流动能的进一步减小。其他许多研究者也得到了类似的结论。对搅拌槽内计算结果的改善需要采用新的计算方法才能从根本上提高计算结果的准确性。随着计算机计算能力的提高,大涡模拟和直接数值模拟方法用来研究搅拌槽内的流动场已逐步开展起来,并将逐步取代现行的计算方法,计算结果准确性的提高将会实现。4预报与实验的局限性RNGk-ε模型在耗散率方程中通过系数C*1引入描述流场畸变效应的附加源项后,在一定程度上改善了对桨叶附近各向异性湍流的预报能力,对速度场的预报比k-ε模型更加准确。但是由于RNGk-ε模型只不过是标准k-ε模型的一种修正形式,其并未放弃各向同性的假定,采用的仍然是标量形式的有效涡粘性系数概念,致使该模型仍具有较大的局限性,尤其是对桨叶附近湍流动能的预报与实验仍然差距很大。因此,要想改善对桨叶附近流动场的预报结果,尤其是对湍流动能的预报,有必要完全放弃k-ε模型,转而采用基于各向异性湍流的应力模型如代数应力模型或雷诺应力模型,或采用其他方法如大涡模拟,直接数值模拟等。结构参数的测定C——桨叶离底距离,mC1,C2,Cμ——模型常数D——搅拌桨直径,mG——湍流动能发生率,Pa/sH——槽内液位高度,mk——湍流动能,m2/sP——压力,Par——径向坐标,mRe——雷诺