超空化流动数值模拟

超空化流动数值模拟

0空化器特性仿真研究超空气泡法是一种高效、低可见水支持技术的方法。其原理是在水与航行体表面之间形成稳定的气层,达到减小水的粘性阻力,从而提高水下航行体的运行速度。为了能够生成稳定的、高效减阻的、形态可控的超空泡,位于航行体头部的空化器起着决定性作用。空化器形状直接影响超空泡形态及水下航行体减阻效果,因此优化空化器形状是获得稳定高效减阻超空泡的前提,其中圆盘空化器具有控制超空泡形态的特点。因此国内外的学者针对空化器作了大量实验研究即理论分析工作Savchenko利用控制空化器参数来获得稳定的超空泡,并总结出经大量实验验证过的经验公式。Stinebring等对头部为圆盘空化器进行了研究,并给出了空泡长度及厚度与空化器参数关系的公式。国内何友声等对系列头体的空化器作了实验研究,在速度比较低即空泡数比较大的情况下,给出了一系列的经验公式。贾力平、王海斌对空化器与阻力关系进行了研究。袁绪龙、张宇文应用商业软件FLUENT建立多相流CFD模型对典型空化器进行了数值模拟。以往对空化器进行的数值模拟都是针对二维轴对称体空化器进行研究。但是,当考虑重力影响以及带攻角空化器生成的超空泡问题时,二维仿真不能解决此问题,因此本文针对三维圆盘空化器超空化流动进行数值模拟。首先,对在不同空化数下无攻角圆盘空化器生成超空泡进行仿真,将仿真结果同Savchenko经验公式进行对比,总结其变化规律;其次,针对带攻角圆盘空化器进行数值模拟,将模拟结果与无攻角圆盘空化器所得结果进行对比,总结其变化规律。文中主要是针对带攻角圆盘空化器进行研究。为了解决三维数值模拟中庞大计算量问题和获得更好的仿真结果,本文采用构架于深腾1800高性能服务器上的商业软件CFX进行仿真。选择深腾1800高性能服务器是利用其多个高速CPU以及大容量内存的优点解决庞大计算量问题,而择商业软件CFX是因为其高效的并行计算方法,以及可以获得较精确仿真结果。1控制方程和算法运用N-S方程及k-ε湍流模型对均质流场求解,引入Rayleigh-Plesset方程模拟超空泡问题。1.1连续方程∂ρ∂t+∇⋅(ρU)=0,(1)式中:ρ表示密度;U表示速度;t表示时间。1.2流体动力sm广义源项{∂∂t(ρU)+∇⋅(ρU⊗U-μ(∇U+(∇U)Τ))=SΜ；ρ=ΝΡ∑a=1γaρa；μ=ΝΡ∑a=1γaμa‚(2)式中:a是相的序号且1≤a≤Np(Np是相数);μ为流体动力粘度;SM为广义源项;γa为a相体积分数。1.3t+u的选取{∂(ρk)∂t+∇(ρUk)=∇[(μ+μtσk)∇k]+pk-ρε；∂(ρε)∂t+∇(ρUε)=∇[(μ+μtσε)∇ε]+εk(Cε1pk-Cε2ρε).(3)式中:pk=μt∇U(∇U+∇UΤ)-23∇U(3μ1∇U+ρk);μt=Cμρk2ε;Cε1=1.44;Cε2=1.92;Cμ=0.09;σk=1.0;σε=1.3.1.4密度和流场密度自然空化状态是相变过程,为描述均质多相流中此相变过程定义传输方程:∂∂t(γaρa)+∇(γaρaUa)=Sa+Γa,(4)式中:Sa为广义源项;Γa为由于相间质量传递产生的进入a相的单位体积质量源。其中Γa=Νp∑b=1Γab(Γab:由b相到a相单位质量流率,本文中a相为水,b相为水蒸汽),Γab=˙mabAab,(5)式中:˙mab是b相到a相单位界面密度的单位质量流率;Aab为相间界面密度。根据空泡动力学知,Rayleigh-Plesset方程描述空泡与水之间的相变过程。方程如下:RBd2RBdt2+32dRBdt2+2τρLRB=pc-pρL,(6)式中:RB为气泡半径;νL为运动粘度;ρL为液体密度;τ为表面张力;p为流场参考压力;pc为水的蒸汽压。在忽略表面张力以及高次项后化简公式dRBdt=±√23pc-pρL,令蒸汽体积分数γb=VBΝB=43πRB3ΝB,其中NB、VB分别是单位体积气泡数和气泡体积。单位体积总质量传递率m˙lg=ΝBdmBdt=F3γbρgRB23|pc-p|ρLsgn(pc-p),(7)式中:F为通过实验获得的经验系数,Fvap=50,Fcond=0.01.令代γnuc(1-γb)替γb,则m˙lg=F3γnuc(1-γb)ρgRB23|pc-p|ρLsgn(pc-p),(8)式中:γnuc为成核位置的体积率。通常在数值模拟中计算空泡模型时采用人为给定单位体积气泡数NB=9×106(单位:个),而对于RB<1mm的气泡无法得到好的仿真结果。采用γnuc(1-γb)代替γb方法,可以计算小尺寸气泡半径如RB=1μm,即使气泡半径再缩小2个数量级均可以仿真。这是因为RB=1μm,其体积分数γb<10-12,根据方程(8)可知,生成的空泡由成核位置的体积率γnuc决定,这样避免由于计算时产生舍入误差影响仿真结果。1.5流域外边界的确定基于有限元的有限体积法,N-S方程中的对流项采用CFX高精度格式,它利用调整混合因子提高计算精度且平衡收敛精度。以一阶迎风为例,对流项ϕi=ϕu+ξ∇ϕΔζ˙‚(9)式中:ϕi为变量ϕ第i个节点值;ϕu为变量ϕ迎风节点值;ζ为混合因子;∇ϕ为迎风节点的节点梯度;ζ˙为从迎风节点到i节点的矢量。针对不同条件下,ζ会作出相应的调整满足计算要求获得精确的解,例如在i节点临近迎风节点的节点梯度变化剧烈,则混合因子趋近于0,反之趋近于1.边界条件为:上游边界给定来流速度和各流体体积率;下游边界给定压强;流域外边界采用速度为零物面边界,或者根据文献采用上游边界,但是根据模拟研究发现流域外边界采用上游边界时,仿真的收敛效果不如采用物面边界,即RMS残差曲线振荡较剧烈;另一方面,当v>12m/s,阻塞比(即,模型截面面积与流域截面面积之比)小于0.01时,流域外边界采用速度为零的物面边界和采用上游边界模拟结果无差别;当阻塞比大于0.01时,流域截面直径影响空泡的形态。因此,本研究工作均采用小阻塞比的物面边界作为流域的外边界;初始条件:设定初始来流速度和各流体体积分数以及流域初始压力。三维物体数值仿真具有计算量庞大,精度高等特点,这就使得进行高精度数值仿真时普通计算机无法完成任务。因此本仿真实验使用了哈尔滨工业大学能源学院的深腾1800高性能服务器。2下覆盖物的数值模拟利用商业CFD软件CFX,针对三维圆盘空化器模型进行数值模拟,所得结果将为水下超空泡水下高速航行体减阻机理研究提供重要的理论基础和依据。2.1空化器网格划分运动体几何形状如图1所示,(a)图表示无攻角圆盘空化器,圆盘空化器尺寸为R=6cm,L=1cm;(b)图表示带α攻角圆盘空化器正视图。尺寸与无攻角圆盘空化器一致。在仿真时二者表面均认为是速度为零的光滑物面。模型所在坐标系为三维直角坐标系,XY平面为正视图,XZ平面为俯视图,YZ平面为侧视图。对于本文问题,若采用非结构化网格,生成的空泡外边界会出现凸凹不平的毛刺,影响对于仿真结果的判断。若想消除毛刺现象,就要对网格进行细分,但会带来网格质量差和计算量呈指数倍增长的缺点。因此采用基于六面体结构化网格划分方法。无攻角圆盘空化器网格分布图和流场域网格分布图如图2所示。总共六面体网格数为2165124,节点数为2204421,其中圆盘空化器表面网格数为5378.流域总长为14m,空化器前面流域长度为1m,为满足阻塞比要求,流域截面积为圆形,其半径为1.5m.其中流域设定14m是因为在小空化数下空泡的长度可达十几米长。如表1,在速度为100m/s时,空泡长度达到11.04m.此外,根据软件ICEMCFD评定网格质量方法可知,结构化网格的网格质量大于等于0.2足以满足仿真要求。此空化器模型的网格质量为0.23,完全满足仿真要求。带攻角圆盘空化器所在流域的网格数量与其相差无几,只是空化器带有攻角,网格质量为0.2.因此,根据本文研究内容,采用基于六面体结构化网格划分方法,满足数值模拟的必要条件。3初始条件及研究目标的确定空化数是衡量流体空化程度的一个物理量,计算公式σ=(p∞-pc)/0.5ρv2,其中p∞,pc,ρ,v分别为某基准位置上的绝对静压力、空泡内绝对压力、液体密度、相对于物体的未扰动液体的流速。研究中设定初始条件p∞、25℃时的水的pc和ρ.雷诺数反应流体粘性的物理量,其公式Re=v∞Dn/ν,在本文研究中其变化范围为3.56×106～1.19×107.由于空化数较大情况下无法生成超空泡,反之,则由于情况比较复杂超出软件的计算能力,因此设定空化数0.019～0.23为研究目标。由于计算软件定义的上游边界条件和初始条件是以速度或者压力为参数,所以在仿真的时候需要将空化数转换成速度参数,即取速度范围在30～100m/s范围内为研究目标,并以每次仿真速度递增10m/s进行采样。3.1仿真结果与实验结果对比直径Dn=120mm无攻角圆盘空化器,速度范围为30～100m/s,仿真结果取空泡形态参数,即空泡的长度Lc和中截面最大直径Dc,见表1,空泡长度与厚度随着速度地增加而增加,即随空化数的减小而变大。表1是以初始速度为30m/s并以10m/s为步长增加至速度为100m/s的仿真结果.乌克兰专家Savchenko依据实验结果和空泡截面扩大独立性原理,利用自由流线扩展的渐近律,得到并用实验验证了:在自然空化情形下,空泡的基本几何形状是一个近似的椭球体。基本的超空泡尺寸有中截面直径Dc和长度Lc,在2%～4%的精度范围内可以由以下渐近关系式决定:DcDn=Cx0(1+σ)σ,LcDn=1σCx0(1+σ)ln1σ,(10)式中:Dc为中截面直径即超空泡最大直径;Lc为超空泡最大长度;Dn为空化器直径;σ为空化数,Cx0空化器阻力系数。由于本文采用的是圆盘空化器,从文献中可以得知Cx0=0.82.为进一步验证仿真结果的可靠性,将仿真结果与经验公式(1)进行对比,如图3、4所示,由公式所得超空泡长度和中截面最大直径变化曲线可知,其长度和中截面最大直径是随着空化数的增加而减少,仿真结果的数据点大部分落在该曲线上。即使在较大空化数下,略有不同,但是均在2%～4%精度范围内,因此二者具有很好的一致性。其中D′c=Dc/ζ和L′c/ζ分别为经过无量纲化的空泡中截面最大直径与长度,ζ为特征长度。此外,对比结果也说明了使用商业软件CFX可以进行三维自然超空化的数值模拟研究。它的计算结果可以为以后的超空泡实验研究提供参考依据。以速度为70m/s的仿真结果为例,如图5所示,空泡最大中截面直径Dc=0.55m和长度Lc=4.94m,图5为三维视图中的正视图(XY平面)。选取方法是在25℃水蒸气的等高图中选取最小含水蒸汽率为50%处定为轮廓边界,利用软件自身的测量工具测定Dc和Lc.其中设定最小含水蒸汽率为50%是以在上述边界条件下,速度范围在30～100m/s中进行仿真结果与(1)式相比时最匹配为标准而设定的。另一个方面,根据(1)式计算Lc=4.94458m和Dc=0.55286m,两者比较之后,误差在2%～4%精度范围之内,说明二者结果相吻合,其中自然超空泡所取正视图的水蒸汽体积率为50%～100%.图5中,在超空泡尾部呈现双涡管闭合方式,这与文献中关于超空泡尾部闭合方式相一致。如图6所示,是其三维视图中的三维形状的超空泡,取含水蒸气50%的等值曲面图形。速度70m/s时3-D超空泡形状是个椭球体,从其侧面图图3看到的双涡管形状是超空泡尾部形成高逆压梯度区域,迫使部分水倒流,形成回注射流。对于回注射流现象在后面会进一步证明。3.2带攻角圆空化器的模拟采用带攻角圆盘空化器可以为航行体提供必要升力,对于水下航行体航行起控制作用,因此具有很重要的研究意义。3.2.1无攻角部分空化器回注射流仿真条件带攻角圆盘空化器生成超空泡形态可由图7(a)～(b)所示,以60°攻角为例,Dn=12cm圆盘空化器,σ=0.219,航行速度v=30m/s,水蒸汽含汽率为50%～100%.从俯视图和侧视图对比可知两种角度观测到的空泡厚度是不一致的,而长度一致;空泡尾部均有凹入部分,其含汽率为0,表明此处的流体为水,因此可以证明产生了十分明显的回注射流,且形状不同。无攻角圆盘空化器生成超空泡形态可由图7(c)～(d)所示。仿真条件除无攻角外均与上面一致。与带攻角圆盘空化器生成超空泡形态比较,从侧视图和俯视图观察其空泡厚度相等,长度也一致,尾部产生的回注射流形状相同。如图8所示,带攻角圆盘空化器生成的超空泡正面投影为椭圆形,中间区域是回注射流所占位置,产生的位置不在椭圆形中心,且上下不对称,左右对称。而无攻角圆盘空化器生成的超空泡正面投影较均匀,呈正圆形,回注射流的位置在圆心处。回注射流现象的产生是根据汽水百分比含量推断出来的。为了进一步验证该区域是回注射流我们引入矢量图9.从图9中的水流速度矢量方向可以知道,水流绕过空化器到达空泡尾部低压区水流速度方向发生变化,水向相反方向流动,因此进一步证明回注射流现象的发生。如图10所示,3D带攻角圆盘空化器生成带回注射流的空泡,图11所示无攻角圆盘空化器生成带回注射流的空泡,此时二者的空化数均为0.219,所取含水蒸汽率为50%的等值面。由图7～图11可知,带攻角圆盘空化器生成的超空泡长度比无攻角圆盘空化器生成的超空泡长度短。产生回注射流现象与空化器是否带攻角无关。4数值模拟结果对比运用N-S方程及k-ε湍流模型对均质流场求解,并且利用Rayleigh-P