第三章湍流模型第一节前言湍流流动模型很多,但大致可以归纳为以下三类:第一类是湍流输运系数模型,是Boussinesq于1877年针对二维流动提出的,将速度脉动的二阶关联量表示成平均速度梯度与湍流粘性系数的乘积。
根据建立模型所需要的微分方程的数目,可以分为零方程模型(代数方程模型),单方程模型和双方程模型。
(模拟大空间建筑空气流动)μt=0.03874ρvl(模拟通风空调室内的空气流动)比例系数由直接数值模拟的结果拟合而得,其中:v为当地时均速度,l为当地距壁面最近的距离。
第二类是抛弃了湍流输运系数的概念,直接建立湍流应力和其它二阶关联量的输运方程。
第三类是大涡模拟。
前两类是以湍流的统计结构为基础,对所有涡旋进行统计平均。
大涡模拟把湍流分成大尺度湍流和小尺度湍流,通过求解三维经过修正的Navier-Stokes方程,得到大涡旋的运动特性,而对小涡旋运动还采用上述的模型。
实际求解中,选用什么模型要根据具体问题的特点来决定。
参见:湍流模型的选择资料。
FLUENT提供的湍流模型包括:单方程(Spalart-Allmaras)模型、双方程模型(标准κ-ε模型、重整化群κ-ε模型、可实现(Realizable)κ-ε模型)及雷诺应力模型和大涡模拟。
如果要求解该方程,必须模拟该项以封闭方程。
如果密度是变化的流动过程如燃烧问题,我们可以用法夫雷(Favre)平均。
这样才可以求解有密度变化的流动问题。
法夫雷平均就是出了压力和密度本身以外,所有变量都用密度加权平均。
Boussinesq假设的缺点是认为湍流粘性系数tμ是各向同性标量,对一些复杂流动该条件并不是严格成立,所以具有其应用限制性。
另外的方法是求解雷诺应力各分量的输运方程。
这也需要额外再求解一个标量方程,通常是耗散率ε方程。
但是,如果湍流场各向异性很明显,如强旋流动以及应力驱动的二次流等流动中,求解雷诺应力分量输运方程无疑可以得到更好的结果。
ν~的输运方程为:ννννρννρμσνρYxCxxGDtDjbjj-+++=~~)~(1~2~3-9其中,νG是湍流粘性产生项;νY是由于壁面阻挡与粘性阻尼引起的湍流粘性的减少;νσ~和2bC是常数;ν是分子运动粘性系数。
湍流粘性系数用如下公式计算:1~ννρμft=其中,1νf是粘性阻尼函数,定义为:31331ννχχCf+=,并且ννχ~≡。
湍流粘性产生项,νG用如下公式模拟:νρν~~1SCGb=3-10其中,222~~ννfdkSS+≡,而1211ννχχff+-=。
其中,1bC和k是常数,d是计算点到壁面的距离;SijijΩΩ≡2。
ijΩ定义为:-=Ωjiijijxuxu213-11由于平均应变率对湍流产生也起到很大作用,FLUENT处理过程中,定义S为:),0min(ijijprodijSCSΩ-+Ω≡3-12其中,0.2=prodC,ijijijΩΩ≡Ω,ijijijSSS2≡,平均应变率ijS定义为:+=jiijijxuxuS213-13在涡量超过应变率的计算区域计算出来的涡旋粘性系数变小。
这适合涡流靠近涡旋中心的区域,那里只有“单纯”的旋转,湍流受到抑止。
包含应变张量的影响更能体现旋转对湍流的影响。
忽略了平均应变,估计的涡旋粘性系数产生项偏高。
湍流粘性系数减少项νY为:21~=dfCYwwνρν3-14其中,6/1636631++=wwwCgCgf3-15)(62rrCrgw-+=3-1622~~dkSrν≡3-17其中,1wC,2wC,3wC是常数,222~~ννfdkSS+≡。
在上式中,包括了平均应变率对S的影响,因而也影响用S~计算出来的r。
上面的模型常数在FLUENT中默认值为:1335.01=bC,622.02=bC,3/2~=νσ,1.71=νC,νσ~2211/)1(/bbwCkCC++=,3.02=wC,0.23=wC,41.0=k。
壁面条件在壁面,湍流运动粘性ν~设置为零。
当计算网格足够细,可以计算层流底层时,壁面切应力用层流应力-应变关系求解,即:μρττyuuu=3-18如果网格粗错不能用来求解层流底层,则假设与壁面近邻的网格质心落在边界层的对数区,则根据壁面法则:=μρττyuEkuuln13-19其中,k=0.419,E=9.793。
对流传热传质模型在FLUENT中,用雷诺相似湍流输运的概念来模拟热输运过程。
给出的能量方程为:heffijjitpiiiSuxTtckxpEuxEt+++=++)(Pr)]([)(τμρρ3-20式中,E是总能量,effij)(τ是偏应力张量,定义为:ijiieffjiijeffeffijxuxuxuδμμτ-+=32)()(3-21其中,effij)(τ表示粘性加热,耦合求解。
如果默认为分开求解,FLUENT不求解处effij)(τ。
但是可以通过变化“粘性模型”面板上的湍流普朗特数(Prt),其默认值为0.85。
Prt数:由流体物性参数组成的一个无因次数(即无量纲参数)群,表明温度边界层和流动边界层的关系,反映流体物理性质对对流传热过程的影响,它的表达式为:Pr=ν/α=cpμ/k式中,μ为动力粘度;cp为等压比热容;k为热导率;α为热扩散系数(α=λ/ρc)单位:m^2/s,v为运动粘度。
湍流质量输运与热输运类似,默认的Schmidt数是0.7,该值同样也可以在“粘性模型”面板上调节。
Schmidt数:表示动量和质量输运之间的关系:粘性系数与扩散系数的比值标量的壁面处理与动量壁面处理类似,分别选用合适的壁面法则。
综上所述,Spalart-Allmaras模型是相对简单的单方程模型,只需求解湍流粘性的输运方程,并不需要求解当地剪切层厚度的长度尺度。
该模型对于求解有壁面影响流动及有逆压力梯度的边界层问题有很好模拟效果,在透平机械湍流模拟方面也有较好结果。
Spalart-Allmaras模型的初始形式属于对低雷诺数湍流模型,这必须很好解决边界层的粘性影响区求解问题。
在FLUENT中,当网格不是很细时,采用壁面函数来解决这一问题。
当网格比较粗糙时,网格不满足精确的湍流计算要求,用壁面函数也许是最好的解决方案。
另外,该模型中的输运变量在近壁处的梯度要比ε-k中的小,这使得该模型对网格粗糙带来数值误差不太敏感。
但是,Spalart-Allmaras模型不能预测均匀各向同性湍流的耗散。
并且,单方程模型没有考虑长度尺度的变化,这对一些流动尺度变换比较大的流动问题不太适合。
比如,平板射流问题,从有壁面影响流动突然变化到自由剪切流,流场尺度变化明显。
3.3.2标准ε-k模型标准ε-k模型需要求解湍动能及其耗散率方程。
湍动能输运方程是通过精确的方程推导得到,但耗散率方程是通过物理推理,数学上模拟相似原形方程得到的。
该模型假设流动为完全湍流,分子粘性的影响可以忽略。
因此,标准ε-k模型只适合完全湍流的流动过程模拟。
标准ε-k模型的湍动能k和耗散率ε方程为如下形式:MbkiktiYGGxkxdtdk--+++=ρεσμμρ3-22kCGCGkCxxDtDbkikti2231)(ερεεσμμερεεε-+++=3-23在上述方程中,kG表示由于平均速度梯度引起的湍动能产生,bG是用于浮力影响引起的湍动能产生;MY可压速湍流脉动膨胀对总的耗散率的影响。
湍流粘性系数ερμμ2kCt=。
在FLUENT中,作为默认值常数,ε1C=1.44,ε2C=1.92,09.0=μC,湍动能k与耗散率ε的湍流普朗特数分别为kσ=1.0,εσ=1.3。
可以通过调节“粘性模型”面板来调节这些常数值。
3.3.3重整化群κ-ε模型重整化的一般思想是:减少系统的自由度,并在这个缩减的空间中,通过特定的重整化技巧,在迭代过程中保持系统的自由度数不变,并使约化系统最终收敛到真正系统的低能态中。
重整化群κ-ε模型是对瞬时的Navier-Stokes方程用重整化群的数学方法推导出来的模型。
模型中的常数与标准κ-ε模型不同,而且方程中也出现了新的函数或者项。
其湍动能与耗散率方程与标准κ-ε模型有相似的形式:()MbkieffkiYGGxkxDtDk--++=ρεμαρ3-24()RkCGCGkCxxDtDbkieffi--++=2231)(ερεεμαερεεεε3-25kG表示由于平均速度梯度引起的湍动能产生,bG是用于浮力影响引起的湍动能产生;MY可压速湍流脉动膨胀对总的耗散率的影响,这些参数与标准κ-ε模型中相同。
kα和εα分别是湍动能k和耗散率ε的有效湍流普朗特数的倒数。
湍流粘性系数计算公式为:ννννεμρ~1~~72.132dCkd--=3-26其中,μμν/~eff=,100≈νC对上面方程积分,可以精确得到有效雷诺数(涡旋尺度)对湍流输运的影响,这有助于处理低雷诺数和近壁流动问题的模拟。
对于高雷诺数,上面方程可以给出:ερμμ2kCt=,0845.0=μC。
这个结果非常有意思,和标准κ-ε模型的半经验推导给出的常数09.0=μC非常近似。
在FLUENT中,如果是默认设置,用重整化群κ-ε模型时候是针对的高雷诺数流动问题。
如果对低雷诺数问题进行数值模拟,必须进行相应的设置。
重整化群κ-ε模型有旋修正通常,平均运动有旋时候对湍流有重要影响。
FLUENT中重整化群κ-ε模型通过修正湍流粘性系数来考虑了这类影响。
湍流粘性的修正形式为:),,(0εαμμkfsttΩ=3-27其中,0tμ是不考虑有旋计算出来的湍流粘性系数;Ω是FLUENT计算出来的特征旋流数;sα是旋流常数,不同值表示有旋流动的强度不同。
流动可以是强旋或者中等旋度的。
FLUENT默认设置sα=0.05,针对中等旋度的流动问题,对于强旋流动,可以选择较大的值。
湍动能及其耗散率的有效湍流普朗特数倒数的计算公式为:effmolμμαααα=++--3679.006321.003929.23929.23929.13929.13-28式中,0α=1,在高雷诺数流动问题中,1/effmolμμ,393.1==εααk。
湍流耗散率方程右边的R为:kCR23031)/1(εβηηηρημ+-=3-29其中,εη/Sk≡,38.40=η,012.0=β。
为了更清楚体现R对耗散率的影响,我们把耗散率输运方程重写为:()kCkCGCGkCxxDtDbkieffi2*22231)(ερερεεμαερεεεεε--++=3-30则:+=εε2*2CC3031)/1(βηηηρημ+-C3-31在0ηη<的区域,R的贡献为正;*2εC大于ε2C。
以对数区为例,3≈η,0.2*2≈εC,这和标准κ-ε模型中给出的ε2C=1.92接近。
因此,对于弱旋和中等旋度的流动问题,重整化群κ-ε模型给出的结果比标准κ-ε模型的结果要大。
重整化群模型中,42.11=εC,68.12=εC。
3.3.4可实现κ-ε模型可实现κ-ε模型的湍动能及其耗散率输运方程为:MbkjktjYGGxkxDtDk--+++=ρεσμμρ3-32bjttjGCkCkCSCxxDtDεεενεερερεσμμερ31221++-++=3-33其中,+=5,43.0max1ηηC,εη/Sk=在上述方程中,kG表示由于平均速度梯度引起的湍动能产生,bG是用于浮力影响引起的湍动能产生;MY可压速湍流脉动膨胀对总的耗散率的影响。
2C和ε1C是常数;kσ,εσ分别是湍动能及其耗散率的湍流普朗特数。
在FLUENT中,作为默认值常数,ε1C=1.44,2C=1.9,kσ=1.0,εσ=1.2。
可实现κ-ε模型的湍动能的输运方程与标准κ-ε模型和重整化群κ-ε模型有相同的形式,只是模型参数不同。
但耗散率方程有较大不同。
首先耗散率产生项(方程右边第二项)不包含湍动能产生项kG,现在的形式更能体现能量在谱空间的传输。
另外的特色在于耗散率减少项中,不具有奇异性。
并不象标准κ-ε模型模型那样把K放在分母上。
该模型适合的流动类型比较广泛,包括有旋均匀剪切流,自由流(射流和混合层),腔道流动和边界层流动。
对以上流动过程模拟结果都比标准κ-ε模型的结果好,特别是可实现κ-ε模型对圆口射流和平板射流模拟中,能给出较好的射流扩张角。
湍流粘性系数公式为ερμμ2kCt=,这和标准κ-ε模型相同。
不同的是,在可实现κ-ε模型中,μC不再是个常数,而是通过如下公式计算:εμKUAACs*01+=3-34其中,ijijijijSSUΩΩ+=~~*,kijkijijωε2~-Ω=Ω,kijkijijωε-Ω=Ω,ijΩ是isthemeanrate-of–rotationtensorviewedinarotatingreferenceframewiththeangularvelocitykω。
模型常数04.40=A,φcos6=sA,而:)6arccos(31W=φ,式中W=SSSSkjjkij~,ijijSSS=~,)(21jiijijxuxuS+=我们可以发现,μC是平均应变率与旋度的函数。
在平衡边界层惯性底层,可以得到μC=0.09,与标准κ-ε模型中采用底常数一样。
双方程模型中,无论是标准κ-ε模型、重整化群κ-ε模型还是可实现κ-ε模型,三个模型有类似的形式,即都有κ和ε的输运方程,它们的区别在于:1,计算湍流粘性的方法不同;2,控制湍流扩散的湍流Prandtl数不同;3Gk关系不同。
但都包含了相同的表示由于平均速度梯度引起的湍动能产生kG,用于浮力影响引起的湍动能产生bG;可压缩湍流脉动膨胀对总的耗散率的影响MY。
对于重整化群κ-ε模型,α/1Pr=t,pCkμα/Pr/1==。
热膨胀系数pT-=ρρβ1,对于理想气体,浮力引起的湍动能产生项变为:itibxtgG-=ρρμPr3-37在FLUENT程序中,如果有重力作用,并且流场里有密度或者温度的梯度,浮力对湍动能的影响都是存在的。
浮力对耗散率的影响不是很清楚,因此,默认设置中,耗散率方程中的浮力影响不被考虑。
如果要考虑浮力对耗散率的影响,用“粘性模型”面板来控制。
浮力对耗散率影响是用ε3C来体现。
但ε3C并不是常数,而是如下的函数形式:uvCtanh3=ε3-38v是平行于重力方向的速度分量;u是垂直于重力方向的速度分量。
如果流动速度与重力方向相同的剪切流动,ε3C=1,对于流动方向与重力方向垂直的剪切流,ε3C=0。
对于高马赫数的流动问题,可压速性对湍流影响在MY中体现。
22tMMYρε=其中,tM是马赫数,定义为:2akMt=(RTaγ≡是声速)。
默认设置中,只要选择可压速理想气体,可压速效应都是考虑的。
在上述的双方程模型中,对流传热传质模型都是通过雷诺相似湍流动量输运方程得到的。
能量方程形式为:()[]heffijjieffiiiSuxTkxpEuxEt++=++)()(τρρ3-39式中,E是总的能量,effk是有效导热系数;effij)(τ是偏应力张量,定义为:ijiieffjiijeffeffijxuxuxuδμμτ-+=32)(3-40effij)(τ表示的是粘性加热,耦合求解时总是计算。
如果不是耦合求解时候,作为默认设置,并不求解该量。
如果有需要,需在“粘性模型”面板中设置。
对于重整化群κ-ε模型,有效导热系数为:effpeffckμα=3-41α用(3-28)计算,式中,pCkμα/Pr/10==。
事实上,α随着effmolμμ/的变化而变化,这是重整化群κ-ε模型的一个优点,因为实验中证明,湍流普朗特数随分子普朗特数及湍流而变化。
湍流质量输运处理过程与能量输运过程类似。
对于标准κ-ε模型和可实现的κ-ε模型,默认的Schmidt数是0.7,重整化群模型中,是通过方程3-28来计算的,其中,Sc/10=α,Sc是分子Schimidt数。
Schmidt数:表示动量和质量输运之间的关系:粘性系数与扩散系数的比值3.3.5雷诺应力模型(RSM)雷诺应力:湍流动量输送的切向应力雷诺应力模型是求解雷诺应力张量的各个分量的输运方程。
具体形式为:=+)()(jikkjiuuUxuutρρ对流项ijC[]+++-jikkjikikjkjikuuxxuupuuuxμδδρ)(湍流扩撒项TijD分子扩散LijD()θθρβρijjikikjkjkiugugxUuuxUuu+-+-应力产生项ijP浮力产生项目ijGkjkiijjixuxuxuxup-++μ2压力应变项ijΦ耗散项ijε()jkmmiikmmjkuuuuεερ+Ω-23-42系统旋转产生项ijF上面方程中,ijC,LijD,ijP,ijF不需要模拟,而TijD,ijG,ijΦ,ijε需要模拟以封闭方程。
下面简单对几个需要模拟项的模拟。
TijD可以用DelayandHarlow[L38]的梯度扩散模型来模拟,即:=ljilkksTijxuuuukxCDερ3-43但这个模型会导致数值不稳定,因此FLUENT程序中采用标量湍流扩散模型:=kjiktkTijxuuxDσμ3-44式中,湍流粘性系数用ερμμ2kCt=来计算,根据LienandLeschziner[L98],82.0=kσ,这和标准κ-ε模型中选取1.0有所不同。
根据GibsonandLaunder[L58],Fu[L55],Launder[L88,L89],压力应变项ijΦ可以分解为三项,即:wijijijijΦ+Φ+Φ=Φ2,1,3-451,ijΦ,2,ijΦ和ijwΦ分别是慢速项,快速项和壁面反射项。
--=ΦkuukCijjiijδερ3211,,常数8.11=C。
()()-+--++-=ΦCGPCGFPCijijijijijijδ3222,,60.02=C,kkPP21=,kkGG21=,kkCC21=。
壁面反射项用于重新分布近壁的雷诺正应力分布,主要是减少垂直于壁面的雷诺正应力,增加平行于壁面的雷诺正应力。
默认设置时候,FLUENT不计算ijwΦ。
如果需要计算时候,在“粘性模型”面板中设置。
这一过程只有在选择双层流模型时候,在“粘性模型”面板上调节。
ijijijijmnmnkjikijijkSbbCCbbbbCbPCCρδρερε)()31()(*332*11-+-++-=Φ)()32(54ikjkjkikijmnmnikjkjkikbbkCSbSbSbkCΩ+Ω+-++ρδρ3-48式中,ijb是雷诺应力各向异性张量,定义为:+--=kkuubijjiijρδρρ2323-49平均应变率ijS定义为:+=jiijijxuxuS21;-=Ωjiijijxuxu21;模型常数4.31=C,8.1*1=C,2.42=C,8.03=C,3.1*3=C,25.14=C,4.05=C。
二阶压力应变模型不需要考虑壁面反射影响去模拟对数区湍流边界层过程。
浮力对湍流的影响浮力引起的产生项模拟为:+=ijjitijxTgxTgtGPrμβ3-50其中,Prt是能量的湍流普朗特数,默认设置值为0.85。
对于理想气体,把热膨胀系数的定义代入上式,得:+-=ijjitijxgxgtGρρρμPr3-51耗散项ijε的模拟耗散张量ijε模拟为:)(32MijijY+=ρεδε3-52式中,22tMMYρε=,tM是马赫数;标量耗散率ε用标准k-ε模型中的采用的耗散率输运方程求解。
雷诺应力模型的边界条件在流场进口,雷诺应力模型需要各个雷诺应力分量和湍动能耗散率的值。
这些值可以直接输入,也可以湍流强度和特征长度来计算。
在壁面,雷诺应力模型通过壁面函数,给出各个雷诺应力分量和耗散率的值。
雷诺应力模型的能量与质量输运方程在雷诺应力模型中,对流传热传质模型都是通过雷诺相似湍流动量输运方程得到的。
能量方程形式为:()[]heffijjitpiiiSuxTtckxpEuxEt+++=++)()Pr()(τμρρ3-53式中,E是总的能量;effij)(τ是偏应力张量,定义为:ijiieffjiijeffeffijxuxuxuδμμτ-+=32)(3-54effij)(τ表示的是粘性加热,耦合求解时总是计算。
如果不是耦合求解时候,作为默认设置,并不求解该量,并且Prt=0.85。
3.3.6大涡模拟(LES)大涡模拟,英文简称LES(Largeeddysimulation),是近几十年才发展起来的一个流体力学中重要的数值模拟研究方法。
它区别于直接数值模拟(DNS)和雷诺平均(RANS)方法。
其基本思想是通过精确求解某个尺度以上所有湍流尺度的运动,从而能够捕捉到RANS方法所无能为力的许多非稳态,非平衡过程中出现的大尺度效应和拟序结构,同时又克服了直接数值模拟由于需要求解所有湍流尺度而带来的巨大计算开销的问题,因而被认为是最具有潜力的湍流数值模拟发展方向。
由于计算耗费依然很大,目前大涡模拟还无法在工程上广泛应用,但是大涡模拟技术对于研究许多流动机理问题提供了更为可靠的手段,可为流动控制提供理论基础,并可为工程上广泛应用的RANS方法改进提供指导。
最大长度尺度通常为平均流动的特征长度尺度。
最小尺度为Komogrov尺度。
LES的基本假设:1,动量、能量、质量及其它标量主要由大涡输运;2,流动的几何和边界条件决定了大涡的特性,而流动特性主要在大涡中体现;3,小尺度涡旋受几何和边界条件影响较小,并且各向同性;大涡模拟过程中,直接求解大涡,小尺度涡旋模拟,从而使得网格要求比DNS低。
3.3.6.1大涡模拟的控制方程LES的控制方程是对Navier-Stokes方程在波数空间或者物理空间进行过滤得到的。
过滤的过程是去掉比过滤宽度或者给定物理宽度小的涡旋,从而得到大涡旋的控制方程。
在FLUENT中,大涡模拟只能针对不可压流体(当然并非说是密度是常数)的流动。
3.3.6.2亚网格模型由于LES中亚网格应力项是未知的,并且需要模拟以封闭方程。
目前,采用比较多的亚网格模型为涡旋粘性模型,形式为:ijtijkkijSμδττ231-=-3-61式中,tμ是亚网格湍流粘性系数;ijS是求解尺度下的应变率张量,定义为:+=ijjiijxuxuS213-62求解亚网格湍流粘性系数tμ时,FLUENT提供了两种方法。
第一,Smagorinsky-Lilly模型;第二,基于重整化群的亚网格模型。
最基本的亚网格模型是Smagorinsky[L145]最早提出的,Lilly[L99]把它进行了改善,这就是今天的Smagorinsky-Lilly模型。
该模型的涡粘性计算方程为:SLst2ρμ=3-63式中,sL是亚网格的混合长度;ijijSSS2≡。
sC是Smagorinsky常数,则亚网格混合长度sL可以用下式计算。