湍流模型.pdf

24 第三章，湍流模型第三章，湍流模型 第一节， 前言 湍流流动模型很多，但大致可以归纳为以下三类 第一类是湍流输运系数模型，是 Boussinesq 于 1877 年针对二维流动提出的，将速度脉动 的二阶关联量表示成平均速度梯度与湍流粘性系数的乘积。即 2 1 21 x u uu t ∂ ∂ ′′−ρ 3－1 推广到三维问题，若用笛卡儿张量表示，即有 ij i j j i tji k x u x u uuδρρ 3 2 −         ∂ ∂ ∂ ∂ ′′− 3－2 模型的任务就是给出计算湍流粘性系数 t 的方法。根据建立模型所需要的微分方程的数 目，可以分为零方程模型（代数方程模型） ，单方程模型和双方程模型。 第二类是抛弃了湍流输运系数的概念，直接建立湍流应力和其它二阶关联量的输运方程。 第三类是大涡模拟。前两类是以湍流的统计结构为基础，对所有涡旋进行统计平均。大 涡模拟把湍流分成大尺度湍流和小尺度湍流， 通过求解三维经过修正的 Navier-Stokes 方程， 得 到大涡旋的运动特性，而对小涡旋运动还采用上述的模型。 实际求解中，选用什么模型要根据具体问题的特点来决定。选择的一般原则是精度要高， 应用简单，节省计算时间，同时也具有通用性。 FLUENT 提供的湍流模型包括单方程（Spalart-Allmaras）模型、双方程模型（标准κ- ε模型、重整化群κ- ε模型、可实现 R e a l i z a b l e κ- ε模型）及雷诺应力模型和大涡模拟。 湍流模型种类示意图 Zero-Equation Models One-Equation Models Spalart-Allmaras Two-Equation Models Standard k-ε ε RNG k-ε ε Realizable k-ε ε Reynolds-Stress Model Large-Eddy Simulation Direct Numerical Simulation 包含更多包含更多 物理机理物理机理 每次迭代每次迭代 计算量增加计算量增加 FLUENT 提 供的模型选 择 RANS-based models 25 第二节，平均量输运方程 雷诺平均就是把 Navier-Stokes 方程中的瞬时变量分解成平均量和脉动量两部分。对于速 度，有 iii uuu′ 3－3 其中， i u和 i u′ 分别是平均速度和脉动速度（i1,2,3） 类似地，对于压力等其它标量，我们也有 φφφ′ 3－4 其中，φ表示标量，如压力、能量、组分浓度等。 把上面的表达式代入瞬时的连续与动量方程，并取平均（去掉平均速度 i u上的横线） ，我 们可以把连续与动量方程写成如下的笛卡儿坐标系下的张量形式 0 ∂ ∂ ∂ ∂ i i u xt ρ ρ 3－5 ji jl l ij i j j i ji i uu xx u x u x u xx p Dt Du ′′− ∂ ∂                 ∂ ∂ − ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ −ρδρ 3 2 3－6 上面两个方程称为雷诺平均的 Navier-Stokes（RANS）方程。他们和瞬时 Navier-Stokes 方 程有相同的形式，只是速度或其它求解变量变成了时间平均量。额外多出来的项 jiu u′′−ρ是雷 诺应力，表示湍流的影响。如果要求解该方程，必须模拟该项以封闭方程。 如果密度是变化的流动过程如燃烧问题，我们可以用法夫雷（Favre）平均。这样才可以 求解有密度变化的流动问题。法夫雷平均就是出了压力和密度本身以外，所有变量都用密度加 权平均。变量的密度加权平均定义为 ρρ/ ΦΦ 3－7 符号～表示密度加权平均；对应于密度加权平均值的脉动值用 Φ′ ′ 表示，即有 Φ′ ′ ΦΦ 。很显然，这种脉动值的简单平均值不为零，但它的密度加权平均值等于零，即 0≠ Φ′ ′ ， 0 Φ′ ′ ρ Boussinesq 近似与雷诺应力输运模型近似与雷诺应力输运模型 为了封闭方程，必须对额外项雷诺应力 jiu u′′−ρ进行模拟。一个通常的方法是应用 Boussinesq 假设，认为雷诺应力与平均速度梯度成正比，即 ij i i t i j j i tji x u k x u x u uuδρρ 3 2 ∂ ∂ −         ∂ ∂ ∂ ∂ ′′− 3－8 Boussinesq 假设被用于 Spalart-Allmaras 单方程模型和ε−k双方程模型。Boussinesq 近似 的好处是与求解湍流粘性系数有关的计算时间比较少，例如在 Spalart-Allmaras 单方程模型中， 只多求解一个表示湍流粘性的输运方程；在ε−k双方程模型中，只需多求解湍动能 k 和耗散 率ε两个方程，湍流粘性系数用湍动能 k 和耗散率ε的函数。Boussinesq 假设的缺点是认为湍 26 流粘性系数 t 是各向同性标量，对一些复杂流动该条件并不是严格成立，所以具有其应用限 制性。 另外的方法是求解雷诺应力各分量的输运方程。这也需要额外再求解一个标量方程，通常 是耗散率ε方程。这就意味着对于二维湍流流动问题，需要多求解 4 个输运方程，而三维湍流 问题需要多求解 7 个方程，需要比较多的计算时间，对计算机内存也有更高要求。 在许多问题中，Boussinesq 近似方法可以得到比较好的结果，并不一定需要花费很多时间 来求解雷诺应力各分量的输运方程。但是，如果湍流场各向异性很明显，如强旋流动以及应力 驱动的二次流等流动中，求解雷诺应力分量输运方程无疑可以得到更好的结果。 第三节， 湍流模型第三节， 湍流模型 3.3.1 单方程（单方程（Spalart-Allmaras）模型）模型 Spalart-Allmaras 模型的求解变量是ν ，表征出了近壁（粘性影响）区域以外的湍流运动粘 性系数。ν 的输运方程为 ν ν ν ν ρ ν νρ σ ν ρY x C xx G Dt D j b jj −                 ∂ ∂           ∂ ∂ ∂ ∂ 1 2 3－9 其中， ν G是湍流粘性产生项； ν Y是由于壁面阻挡与粘性阻尼引起的湍流粘性的减少； ν σ 和 2b C是常数；ν是分子运动粘性系数。 湍流粘性系数用如下公式计算 1 ν νρf t 其中， 1ν f是粘性阻尼函数，定义为 3 1 3 3 1 ν ν χ χ C f ，并且 ν ν χ ≡。 湍流粘性产生项， ν G用如下公式模拟 νρ ν 1 SCG b 3－10 其中， 2 22 ν ν f dk SS≡，而 1 2 1 1 ν ν χ χ f f −。其中， 1b C和 k 是常数，d 是计算点到 壁面的距离；S ijijΩ Ω≡2。 ij Ω定义为         ∂ ∂ − ∂ ∂ Ω j i i j ij x u x u 2 1 3－11 由于平均应变率对湍流产生也起到很大作用，FLUENT 处理过程中，定义 S 为 , 0min ijijprodij SCSΩ−Ω≡ 3－12 其中，0 . 2 prod C， ijijij ΩΩ≡Ω， ijijij SSS2≡，平均应变率 ij S定义为         ∂ ∂ ∂ ∂ j i i j ij x u x u S 2 1 3－13 27 在涡量超过应变率的计算区域计算出来的涡旋粘性系数变小。 这适合涡流靠近涡旋中心的 区域，那里只有“单纯”的旋转，湍流受到抑止。包含应变张量的影响更能体现旋转对湍流的 影响。忽略了平均应变，估计的涡旋粘性系数产生项偏高。 湍流粘性系数减少项 ν Y为 2 1       d fCY ww ν ρ ν 3－14 其中， 6/1 6 3 6 6 3 1       w w w Cg C gf 3－15 6 2 rrCrg w − 3－16 22 dkS r ν ≡ 3－17 其中， 1w C， 2w C， 3w C是常数， 2 22 ν ν f dk SS≡。在上式中，包括了平均应变率对 S 的影响，因而也影响用S 计算出来的 r。 上面的模型常数在 FLUENT 中默认值为1335 . 0 1 b C，622 . 0 2 b C，3/2 ν σ， 1 . 7 1ν C， ν σ 2 2 11 / 1 / bbw CkCC，3 . 0 2 w C，0 . 2 3 w C，41. 0k。 壁面条件 在壁面，湍流运动粘性ν 设置为零。当计算网格足够细，可以计算层流底层时，壁面切应 力用层流应力－应变关系求解，即 ρ τ τ yu u u 3－18 如果网格粗错不能用来求解层流底层，则假设与壁面近邻的网格质心落在边界层的对数 区，则根据壁面法则         ρ τ τ yu E ku u ln 1 3－19 其中，k0.419，E9.793。 对流传热传质模型对流传热传质模型 在 FLUENT 中，用雷诺相似湍流输运的概念来模拟热输运过程。给出的能量方程为 heffijj i tp i i i Su x T t c k x pEu x E t         ∂ ∂         ∂ ∂ ∂ ∂ ∂ ∂ Pr ][τ ρρ 3－20 式中，E 是总能量， effij τ是偏应力张量，定义为 28 ij i i eff j i i j effeffij x u x u x u δτ ∂ ∂ − ∂ ∂ ∂ ∂ 3 2 3－21 其中， effij τ表示粘性加热， 耦合求解。 如果默认为分开求解， FLUENT 不求解处 effij τ。 但是可以通过变化“粘性模型”面板上的湍流普朗特数（Prt） ，其默认值为 0.85。 湍流质量输运与热输运类似，默认的 Schmidt 数是 0.7，该值同样也可以在“粘性模型” 面板上调节。 标量的壁面处理与动量壁面处理类似，分别选用合适的壁面法则。 综上所述，Spalart-Allmaras 模型是相对简单的单方程模型，只需求解湍流粘性的输运方 程，并不需要求解当地剪切层厚度的长度尺度。该模型对于求解有壁面影响流动及有逆压力梯 度的边界层问题有很好模拟效果，在透平机械湍流模拟方面也有较好结果。 Spalart-Allmaras 模型的初始形式属于对低雷诺数湍流模型， 这必须很好解决边界层的粘性 影响区求解问题。在 FLUENT 中，当网格不是很细时，采用壁面函数来解决这一问题。当网 格比较粗糙时，网格不满足精确的湍流计算要求，用壁面函数也许是最好的解决方案。另外， 该模型中的输运变量在近壁处的梯度要比ε−k中的小，这使得该模型对网格粗糙带来数值误 差不太敏感。 但是，Spalart-Allmaras 模型不能预测均匀各向同性湍流的耗散。并且，单方程模型没有考 虑长度尺度的变化， 这对一些流动尺度变换比较大的流动问题不太适合。 比如， 平板射流问题， 从有壁面影响流动突然变化到自由剪切流，流场尺度变化明显。 3.3.2 标准标准ε−k模型模型 标准ε−k模型需要求解湍动能及其耗散率方程。湍动能输运方程是通过精确的方程推导 得到，但耗散率方程是通过物理推理，数学上模拟相似原形方程得到的。该模型假设流动为完 全湍流，分子粘性的影响可以忽略。因此，标准ε−k模型只适合完全湍流的流动过程模拟。 标准ε−k模型的湍动能 k 和耗散率ε方程为如下形式 Mbk ik t i YGG x k xDt Dk −−       ∂ ∂         ∂ ∂ ρε σ ρ 3－22 k CGCG k C xxDt D bk ik t i 2 231 ε ρ εε σ ε ρ εεε −       ∂ ∂         ∂ ∂ 3－23 在上述方程中， k G表示由于平均速度梯度引起的湍动能产生， b G是用于浮力影响引起 的湍动能产生； M Y可压速湍流脉动膨胀对总的耗散率的影响。湍流粘性系数 ε ρ 2 k C t 。 在 FLUENT 中，作为默认值常数， ε1 C＝1.44， ε2 C1.92，09. 0 C，湍动能 k 与耗散 率ε的湍流普朗特数分别为 k σ＝1.0， ε σ＝1.3。可以通过调节“粘性模型”面板来调节这些 常数值。 29 3 . 3 . 3 重整化群κ- ε模型 3 . 3 . 3 重整化群κ- ε模型 重整化群κ- ε模型是对瞬时的 N a v i e r - S t o k e s方程用重整化群的数学方法推导出来的模 型。模型中的常数与标准κ- ε模型不同，而且方程中也出现了新的函数或者项。其湍动能与 耗散率方程与标准κ- ε模型有相似的形式 Mbk i effk i YGG x k xDt Dk −−       ∂ ∂ ∂ ∂ ρεαρ 3－24 R k CGCG k C xxDt D bk i eff i −−       ∂ ∂ ∂ ∂ 2 231 ε ρ εε α ε ρ εεεε 3－25 k G表示由于平均速度梯度引起的湍动能产生， b G是用于浮力影响引起的湍动能产生； M Y可 压速湍流脉动膨胀对总的耗散率的影响，这些参数与标准κ- ε模型中相同。 k α和 ε α分别是 湍动能 k 和耗散率ε的有效湍流普朗特数的倒数。 湍流粘性系数计算公式为 ν νν ν ε ρ 1 72 . 1 3 2 d C k d −−         3 －2 6 其中，ν/ eff ，100≈ ν C 对上面方程积分，可以精确得到有效雷诺数（涡旋尺度）对湍流输运的影响，这有助于 处理低雷诺数和近壁流动问题的模拟。 对于高雷诺数，上面方程可以给出 ε ρ 2 k C t ，0845. 0 C。这个结果非常有意 思，和标准κ- ε模型的半经验推导给出的常数09. 0 C非常近似。 在 FLUENT 中， 如果是默认设置， 用重整化群κ- ε模型时候是针对的高雷诺数流动问题。 如果对低雷诺数问题进行数值模拟，必须进行相应的设置。 重整化群κ- ε模型有旋修正 重整化群κ- ε模型有旋修正 通常，平均运动有旋时候对湍流有重要影响。F L U E N T 中重整化群κ- ε模型通过修正湍流 粘性系数来考虑了这类影响。 湍流粘性的修正形式为 ,, 0 ε α k f stt Ω 3 －2 7 其中， 0t 是不考虑有旋计算出来的湍流粘性系数； Ω是 F L U E N T 计算出来的特征旋流数； s α是 旋流常数，不同值表示有旋流动的强度不同。流动可以是强旋或者中等旋度的。F L U E N T 默认设 置 s α＝0 . 0 5 , 针对中等旋度的流动问题，对于强旋流动，可以选择较大的值。 湍动能及其耗散率的有效湍流普朗特数倒数的计算公式为 30 eff mol α α α α − − 3679. 0 0 6321. 0 0 3929. 2 3929. 2 3929. 1 3929. 1 3－28 式中， 0 α＝1 ，在高雷诺数流动问题中，1/〈〈 effmol ，393 . 1 ε ααk。 湍流耗散率方程右边的 R 为 k C R 2 3 0 3 1 /1 ε βη ηηρη − 3 －2 9 其中，εη/Sk≡，38 . 4 0 η，012. 0β。 为了更清楚体现 R 对耗散率的影响，我们把耗散率输运方程重写为 k C k CGCG k C xxDt D bk i eff i 2 * 2 2 231 ε ρ ε ρ εε α ε ρ εεεεε −−       ∂ ∂ ∂ ∂ 3－30 则 εε2 * 2 CC 3 0 3 1 /1 βη ηηρη −C 3 －3 1 在 0 ηη x x 3－66 V 是计算控制体体积；重整化群常数157. 0 rng C,而常数 C100。 对于高雷诺数流动（ t ） ， teff ≅，基于重整化群理论的亚网格模型就与 Smagorinsky-Lilly 模型相同，只是模型常数有区别。在流动场的低雷诺数区域，上面的函数就 小于零，从而只有分子粘性起作用。所以，基于重整化群理论的亚网格模型对流动转捩和近壁 流动问题有较好模拟效果。 3.3.6.3 大涡模拟的边界条件大涡模拟的边界条件 对于给定进口速度边界条件，速度等于各个方向分量与随机脉动量的和，即 uIuu ii ψ 其中，I 是脉动强度，ψ是高斯随机数，满足0ψ，1′ψ。 如果网格足够密并可以求解层流底层的流动的话， 壁面切应力采用线性应力应变关系， 即 ρ τ τ yu u u 3－67 如果网格不够细，则假定与壁面邻近网格质心落在边界层对数区内，则 ln 1 ρ τ τ yu E ku u 3－68 其中，k0.418，E9.793。 39 表 3－1，雷诺平均模型的比较 模型名字 优点 缺点 Spalart-Allmaras 计算量小，对一定复杂程度的 边界层问题有较好效果 计算结果没有被广泛测试， 缺少子模型，如考虑燃烧或 浮力问题 标准ε−k 应用多，计算量合适，有较多数据 积累和相当精度 对于流向有曲率变化， 较强压力梯度 有旋问题等复杂流动模拟效果欠缺 RNG ε−k 能模拟射流撞击，分离流，二次流， 旋流等中等复杂流动 受到涡旋粘性各向同性假设限制 Realizable ε−k 和 RNG 模型差不多，还可以模拟圆 口射流问题 受到涡旋粘性各向同性假设限制 雷诺应力模型 考虑的物理机理更仔细， 包括了湍流 各向异性影响 CPU 时间长（2～3 倍） ，动量和 湍流量高度耦合。 第四节，湍流模型算例及其设置第四节，湍流模型算例及其设置 湍流模型设置命令Define-model-viscous 40 无粘，层流和湍流        湍流模型选项          近壁处理方法选择        附加湍流选项      41 算例分析有换热的腔道流动问题算例分析有换热的腔道流动问题 步骤 1， 检查是否湍流5980Re Dh 2， 选择低雷诺数湍流模型－RNG ε−k模型；壁面处理用非平衡壁面函数，考虑压力梯 度影响 3， 网格划分 （1）四边形网格； （2）由于在靠近水平板处，垂直方向梯度较大，则近壁 网格加密，并保证第一个控制体在对数区内； （3）变化流向网格间距，用于捕捉边界 层发展； （4） ，根据计算结果，自动调节网格，用于进一步计算温度梯度。 计算结果 adiabatic wall cold air V 50 fpm T 0 F constant temperature wall T 100 F insulation 1 ft 1 ft 10 ft P Velocity contours BLs on upper lower surfaces accelerate the core flow Temperature contours Important that thermal BL was accurately resolved as well P 42 算例二，圆柱绕流算例二，圆柱绕流 步骤 1， 确定雷诺数，24600Re D 2， 钝体绕流，后面有不稳定的涡旋脱落。采用 RNG ε−k模型，壁面处理是双层区模 型； 3， 网格处理近壁网格加密，由于是双层区模型，需要网格划分到1 y wall wall 1 ft 2 ft 2 ft air V 4 fps Compute drag coefficient of the cylinder 5 ft 14.5 ft 43 计算圆柱绕流的涡旋脱落过程 Contours of effective viscosity eff t CD 0.53 Strouhal Number 0.297 U D St τ ≡ where