主题022：RANS模型与应用

1. 引言

在计算流体力学(CFD)领域，湍流模拟一直是核心挑战之一。湍流是一种高度复杂的三维非定常流动现象，其特征是多尺度涡旋结构和能量级串。对于工程应用而言，直接数值模拟(DNS)虽然能够精确求解所有尺度的湍流运动，但其计算成本极高，难以应用于实际工程问题。大涡模拟(LES)虽然比DNS经济，但仍需要相当大的计算资源。因此，基于时间平均的雷诺平均纳维-斯托克斯(RANS)方法成为工程湍流模拟的主流方法。

RANS方法通过将流动变量分解为时均量和脉动量，对Navier-Stokes方程进行时间平均，得到描述时均流动的方程组。由于时间平均过程引入了新的未知量——雷诺应力，因此需要建立湍流模型来封闭方程组。自20世纪70年代以来，各种RANS湍流模型被提出并广泛应用于工程领域，从简单的零方程模型到复杂的多方程模型，为不同类型的流动问题提供了有效的解决方案。

在结构风振仿真中，RANS模型被广泛用于计算建筑物周围的风场分布、风压系数、以及评估风荷载特性。虽然RANS模型在预测复杂分离流动和涡脱落方面存在一定局限性，但对于定常或准定常风荷载分析、初步设计阶段的快速评估等应用场景，RANS模型仍然具有重要价值。本章将系统介绍RANS方法的基本理论、各类湍流模型的特点与应用，以及在结构风工程中的具体应用案例。

2. RANS方法的理论基础

2.1 雷诺分解

RANS方法的核心是雷诺分解，即将瞬时流动变量分解为时均量和脉动量两部分。对于任意流动变量φ，可以表示为：

$$\phi =\bar{\phi }+{\phi }^{\prime }$$

其中，$\bar{\phi }$是时均值，${\phi }^{\prime }$是脉动量。时均值的定义为：

$$\bar{\phi }=\underset{T\to \infty }{\lim }\frac{1}{T}{\int }_{t_0}^{t_0+T}\phi dt$$

在实际计算中，取足够长的时间间隔T，使得平均结果不再随T变化。对于定常流动，时间平均等价于系综平均；对于非定常流动，可以采用相平均或滤波平均。

雷诺分解具有以下重要性质：

1. 时均量的时均等于自身：
   $$\overline{\bar{\phi }}=\bar{\phi }$$

2. 脉动量的时均等于零：
   $$\overline{{\phi }^{\prime }}=0$$

3. 时均量与脉动量乘积的时均等于零：
   $$\overline{\bar{\phi }{\phi }^{\prime }}=\bar{\phi }\cdot \overline{{\phi }^{\prime }}=0$$

4. 两个脉动量乘积的时均一般不为零：
   $$\overline{{\phi }^{\prime }{\psi }^{\prime }}\ne 0$$

2.2 雷诺平均Navier-Stokes方程

对不可压缩流动的Navier-Stokes方程进行时间平均，得到RANS方程组。

连续性方程：

$$\frac{\partial {\bar{u}}_i}{\partial x_i}=0$$

动量方程：

$$\frac{\partial {\bar{u}}_i}{\partial t}+{\bar{u}}_j\frac{\partial {\bar{u}}_i}{\partial x_j}=-\frac{1}{\rho }\frac{\partial \bar{p}}{\partial x_i}+\nu \frac{{\partial }^2{\bar{u}}_i}{\partial x_j\partial x_j}-\frac{\partial }{\partial x_j}(\overline{u_i^{\prime }{u}_j^{\prime }})$$

其中，$\overline{u_i^{\prime }{u}_j^{\prime }}$是雷诺应力张量，表示湍流脉动对时均流动的影响。雷诺应力是一个对称张量，包含6个独立分量：

$${\tau }_{ij}^{(t)}=-\rho \overline{u_i^{\prime }{u}_j^{\prime }}=\left(\begin{array}{ccc}-\rho \overline{u^{\prime 2}}& -\rho \overline{u^{\prime }{v}^{\prime }}& -\rho \overline{u^{\prime }{w}^{\prime }}\\ -\rho \overline{u^{\prime }{v}^{\prime }}& -\rho \overline{v^{\prime 2}}& -\rho \overline{v^{\prime }{w}^{\prime }}\\ -\rho \overline{u^{\prime }{w}^{\prime }}& -\rho \overline{v^{\prime }{w}^{\prime }}& -\rho \overline{w^{\prime 2}}\end{array}\right)$$

2.3 封闭问题

RANS方程组中，未知量包括3个速度分量、压力和6个雷诺应力分量，共10个未知量，而方程只有4个（1个连续性方程+3个动量方程）。因此，需要建立6个额外的关系式来封闭方程组，这就是湍流模型的任务。

湍流模型的基本思想是建立雷诺应力与时均速度场之间的关系。根据建模方式的不同，湍流模型可以分为以下几类：

1. 代数模型（零方程模型）：直接建立雷诺应力与时均速度梯度的代数关系
2. 一方程模型：引入一个湍流特征量（通常是湍动能k）的输运方程
3. 两方程模型：引入两个湍流特征量的输运方程，如k-ε模型、k-ω模型
4. 雷诺应力模型（RSM）：直接求解雷诺应力的输运方程
5. 非线性涡粘模型：在涡粘假设基础上引入非线性修正

2.4 涡粘假设

大多数工程湍流模型基于Boussinesq涡粘假设，即雷诺应力与平均应变率成正比：

$$-\overline{u_i^{\prime }{u}_j^{\prime }}={\nu }_t\left(\frac{\partial {\bar{u}}_i}{\partial x_j}+\frac{\partial {\bar{u}}_j}{\partial x_i}\right)-\frac{2}{3}k{\delta }_{ij}$$

其中，${\nu }_t$是涡粘系数（湍流粘性系数），k是湍动能，${\delta }_{ij}$是Kronecker delta符号。

涡粘系数的量纲为$[L^2/T]$，可以表示为湍流特征速度$u_t$和特征长度$l_t$的乘积：

$${\nu }_t=C{u}_t{l}_t$$

其中C是常数。不同湍流模型的主要区别在于如何确定这两个特征量。

3. 零方程模型

3.1 混合长度理论

零方程模型是最简单的湍流模型，它不需要求解额外的输运方程，直接通过代数关系计算涡粘系数。Prandtl混合长度理论是最早的零方程模型。

Prandtl假设湍流脉动速度与平均速度梯度和混合长度成正比：

$$|u^{\prime }|\approx l_m\left|\frac{d\bar{u}}{dy}\right|$$

涡粘系数表示为：

$${\nu }_t=l_m^2\left|\frac{d\bar{u}}{dy}\right|$$

混合长度$l_m$需要根据流动类型确定。对于边界层流动，通常取：

$$l_m=\kappa y$$

其中，$\kappa =0.41$是von Kármán常数，y是到壁面的距离。在外层区域，混合长度趋于常数：

$$l_m=\lambda \delta$$

其中，$\lambda \approx 0.09$，$\delta$是边界层厚度。

3.2 Cebeci-Smith模型

Cebeci-Smith模型是一种改进的混合长度模型，适用于边界层流动。该模型将边界层分为内层和外层：

内层（$y\le y_c$）：

$${\nu }_t^{(i)}=l^2\left|\frac{\partial \bar{u}}{\partial y}\right|{\gamma }_{tr}$$

其中，$l=\kappa y[1-\exp (-y^{+}/A^{+})]$，$y^{+}=y{u}_{\tau }/\nu$，$A^{+}=26$，${\gamma }_{tr}$是转捩因子。

外层（$y>y_c$）：

$${\nu }_t^{(o)}=\alpha {\bar{u}}_e{\delta }^{\ast }{\gamma }_{tr}{F}_{Kleb}$$

其中，$\alpha =0.0168$，${\bar{u}}_e$是边界层外缘速度，${\delta }^{\ast }$是位移厚度，$F_{Kleb}$是Klebanoff间歇因子。

内外层在$y_c$处平滑过渡，$y_c$是使${\nu }_t^{(i)}={\nu }_t^{(o)}$的位置。

3.3 Baldwin-Lomax模型

Baldwin-Lomax模型是另一种广泛使用的零方程模型，特别适用于分离流动。该模型同样分为内外两层：

内层：

$${\nu }_t^{(i)}=l^2|\omega |$$

其中，$|\omega |$是涡量大小，$l=\kappa y[1-\exp (-y^{+}/A^{+})]$。

外层：

$${\nu }_t^{(o)}=K{C}_{cp}{F}_{wake}{F}_{Kleb}$$

其中，$F_{wake}=\min (y_{max}{F}_{max},C_{wk}{y}_{max}{U}_{dif}^2/F_{max})$，$F=y|\omega |[1-\exp (-y^{+}/A^{+})]$。

模型常数：$\kappa =0.4$，$A^{+}=26$，$K=0.0168$，$C_{cp}=1.6$，$C_{wk}=0.25$。

3.4 零方程模型的特点与应用

优点：

• 计算简单，不需要求解额外的输运方程
• 计算成本低，收敛性好
• 对简单剪切流动（如边界层、管道流）预测较好

缺点：

• 无法考虑湍流的历史效应和对流输运
• 对复杂流动（分离流、回流区）预测较差
• 混合长度需要针对不同流动类型调整

应用：

• 航空器外部流场初步分析
• 简单管道流动和边界层计算
• 作为更复杂模型的初始场

4. 一方程模型

4.1 湍动能方程

一方程模型通过求解湍动能k的输运方程来确定涡粘系数。湍动能定义为：

$$k=\frac{1}{2}(\overline{u^{\prime 2}}+\overline{v^{\prime 2}}+\overline{w^{\prime 2}})$$

湍动能的精确输运方程可以从Navier-Stokes方程推导得到：

$$\frac{\partial k}{\partial t}+{\bar{u}}_j\frac{\partial k}{\partial x_j}=P_k-\epsilon +\frac{\partial }{\partial x_j}\left[\nu \frac{\partial k}{\partial x_j}-\frac{1}{2}\overline{u_i^{\prime }{u}_i^{\prime }{u}_j^{\prime }}-\frac{1}{\rho }\overline{p^{\prime }{u}_j^{\prime }}\right]$$

其中，$P_k=-\overline{u_i^{\prime }{u}_j^{\prime }}\frac{\partial {\bar{u}}_i}{\partial x_j}$是湍动能生成项，$\epsilon =\nu \overline{\frac{\partial u_i^{\prime }}{\partial x_j}\frac{\partial u_i^{\prime }}{\partial x_j}}$是湍动能耗散率。

4.2 模型封闭

一方程模型假设涡粘系数与湍动能和特征长度相关：

$${\nu }_t=C_{\mu }\sqrt{k}L$$

其中，$C_{\mu }$是模型常数，L是湍流长度尺度。长度尺度通常由代数关系给出，如：

$$L=\min (\kappa y,C_L\delta )$$

湍动能耗散率通过：

$$\epsilon =C_D\frac{k^{3/2}}{L}$$

湍流扩散项采用梯度扩散假设：

$$-\frac{1}{2}\overline{u_i^{\prime }{u}_i^{\prime }{u}_j^{\prime }}-\frac{1}{\rho }\overline{p^{\prime }{u}_j^{\prime }}=\frac{{\nu }_t}{{\sigma }_k}\frac{\partial k}{\partial x_j}$$

最终得到模型化的湍动能输运方程：

$$\frac{\partial k}{\partial t}+{\bar{u}}_j\frac{\partial k}{\partial x_j}=P_k-C_D\frac{k^{3/2}}{L}+\frac{\partial }{\partial x_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]$$

4.3 一方程模型的特点

优点：

• 考虑了湍动能的对流和扩散输运
• 比零方程模型更适合非平衡流动
• 计算成本仍然较低

缺点：

• 长度尺度仍需经验确定
• 对复杂流动的适用性有限
• 不如两方程模型通用

代表模型：

• Prandtl-Kolmogorov模型
• Wolfshtein模型
• Baldwin-Barth模型

5. 标准k-ε模型

5.1 模型方程

k-ε模型是目前应用最广泛的湍流模型之一。该模型引入湍动能k和湍动能耗散率ε两个变量，通过求解各自的输运方程来确定涡粘系数。

湍动能输运方程：

$$\frac{\partial k}{\partial t}+{\bar{u}}_j\frac{\partial k}{\partial x_j}=P_k-\epsilon +\frac{\partial }{\partial x_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]$$

耗散率输运方程：

$$\frac{\partial \epsilon }{\partial t}+{\bar{u}}_j\frac{\partial \epsilon }{\partial x_j}=C_{\epsilon 1}\frac{\epsilon }{k}{P}_k-C_{\epsilon 2}\frac{{\epsilon }^2}{k}+\frac{\partial }{\partial x_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]$$

涡粘系数：

$${\nu }_t=C_{\mu }\frac{k^2}{\epsilon }$$

湍动能生成项：

$$P_k={\nu }_t\left(\frac{\partial {\bar{u}}_i}{\partial x_j}+\frac{\partial {\bar{u}}_j}{\partial x_i}\right)\frac{\partial {\bar{u}}_i}{\partial x_j}$$

5.2 模型常数

标准k-ε模型的常数通过拟合简单流动（如衰减各向同性湍流、均匀剪切流、对数律边界层）的实验数据确定：

常数	数值	物理意义
$C_{\mu }$	0.09	涡粘系数常数
$C_{\epsilon 1}$	1.44	生成项系数
$C_{\epsilon 2}$	1.92	耗散项系数
${\sigma }_k$	1.0	k的Prandtl数
${\sigma }_{\epsilon }$	1.3	ε的Prandtl数

这些常数在广泛流动中表现良好，但对于某些特定流动可能需要调整。

5.3 边界条件

入口边界：

• 给定k和ε的值，或给定湍流强度I和特征长度L：
  $$k=\frac{3}{2}(UI{)}^2,\epsilon =C_{\mu }^{3/4}\frac{k^{3/2}}{L}$$

出口边界：

• 通常采用零梯度条件：$\frac{\partial k}{\partial n}=0$，$\frac{\partial \epsilon }{\partial n}=0$

对称边界：

• 法向梯度为零，法向速度为零

固壁边界：

• 使用壁面函数或低雷诺数修正

5.4 壁面处理

壁面函数法：

在壁面附近，采用壁面函数避免求解粘性底层。对于对数律区域（$y^{+}>30$）：

$$u_t=\kappa y{u}_{\tau },u^{+}=\frac{1}{\kappa }\ln (E{y}^{+})$$

壁面剪切应力：

$${\tau }_w=\rho u_{\tau }^2=\rho C_{\mu }^{1/4}{k}^{1/2}\frac{\bar{u}}{u^{+}}$$

壁面处k和ε的边界条件：

$$\frac{\partial k}{\partial n}=0,\epsilon =\frac{C_{\mu }^{3/4}{k}^{3/2}}{\kappa y}$$

低雷诺数修正：

对于需要求解到壁面的情况（$y^{+}<1$），采用低雷诺数k-ε模型，引入衰减函数：

$${\nu }_t=C_{\mu }{f}_{\mu }\frac{k^2}{\tilde{\epsilon }},\epsilon =\tilde{\epsilon }+D$$

其中，$f_{\mu }$是粘性衰减函数，D是近壁修正项。

5.5 标准k-ε模型的特点

优点：

• 简单、稳定、计算成本低
• 对简单剪切流动和回流区预测较好
• 广泛验证，有大量应用经验

缺点：

• 对强逆压梯度流动预测不足
• 对旋转流动、曲率效应敏感
• 对圆形射流的扩展率预测过高
• 近壁区域需要特殊处理

6. RNG k-ε模型

6.1 重整化群理论

RNG（Renormalization Group）k-ε模型基于重整化群理论推导得到。RNG理论通过系统性地消除小尺度涡的影响，从大尺度方程推导湍流模型。

6.2 模型方程

RNG k-ε模型与标准k-ε模型形式相似，但系数不同：

k方程：

$$\frac{\partial k}{\partial t}+{\bar{u}}_j\frac{\partial k}{\partial x_j}=P_k-\epsilon +\frac{\partial }{\partial x_j}\left[{\alpha }_k{\nu }_{eff}\frac{\partial k}{\partial x_j}\right]$$

ε方程：

$$\frac{\partial \epsilon }{\partial t}+{\bar{u}}_j\frac{\partial \epsilon }{\partial x_j}=C_{\epsilon 1}^{\ast }\frac{\epsilon }{k}{P}_k-C_{\epsilon 2}\frac{{\epsilon }^2}{k}+\frac{\partial }{\partial x_j}\left[{\alpha }_{\epsilon }{\nu }_{eff}\frac{\partial \epsilon }{\partial x_j}\right]$$

其中，${\nu }_{eff}=\nu +{\nu }_t$，${\alpha }_k={\alpha }_{\epsilon }\approx 1.393$。

6.3 关键改进

应变率修正：

$$C_{\epsilon 1}^{\ast }=C_{\epsilon 1}-\frac{\eta (1-\eta /{\eta }_0)}{1+\beta {\eta }^3}$$

其中，$\eta =Sk/\epsilon$，$S=\sqrt{2S_{ij}{S}_{ij}}$是应变率大小，${\eta }_0=4.377$，$\beta =0.012$。

有效粘性：

RNG模型给出了低雷诺数情况下的有效粘性公式：

$${\nu }_{eff}=\nu [1+H(\frac{\tilde{\epsilon }}{{\nu }^{3/2}}-100){]}^{1/3}$$

其中，H是Heaviside函数。

6.4 RNG k-ε模型的特点

优点：

• 对高应变率和弯曲流线流动预测更好
• 对分离流和回流区预测改进
• 对室内通风、旋转机械等应用更合适
• 包含低雷诺数效应，无需壁面函数

缺点：

• 计算成本略高于标准k-ε模型
• 对圆形射流预测仍不理想
• 某些情况下可能过度预测湍流生成

7. Realizable k-ε模型

7.1 可实现性约束

标准k-ε模型在某些流动条件下违反可实现性条件，即雷诺应力可能取非物理值。Realizable k-ε模型通过修正涡粘系数和耗散率方程来保证可实现性。

7.2 模型方程

涡粘系数：

$${\nu }_t=C_{\mu }^{\ast }\frac{k^2}{\epsilon }$$

其中，$C_{\mu }^{\ast }$是应变率的函数：

$$C_{\mu }^{\ast }=\frac{1}{A_0+A_s\frac{k{U}^{\ast }}{\epsilon }}$$

$$U^{\ast }=\sqrt{S_{ij}{S}_{ij}+{\tilde{\Omega }}_{ij}{\tilde{\Omega }}_{ij}},A_0=4.04,A_s=\sqrt{6}\cos \varphi$$

ε方程：

$$\frac{\partial \epsilon }{\partial t}+{\bar{u}}_j\frac{\partial \epsilon }{\partial x_j}=C_{1}S\epsilon -C_2\frac{{\epsilon }^2}{k+\sqrt{\nu \epsilon }}+\frac{\partial }{\partial x_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]$$

其中，$C_1=\max (0.43,\frac{\eta }{5+\eta })$，$\eta =Sk/\epsilon$。

7.3 Realizable k-ε模型的特点

优点：

• 满足可实现性条件，数值稳定性好
• 对平面射流、圆形射流、轴对称射流预测一致
• 对分离流和复杂二次流预测改进
• 对强剪切流动和旋转流动更合适

缺点：

• 计算成本略高
• 对壁面边界层需要仔细处理
• 某些情况下可能预测过度分离

8. k-ω模型

8.1 模型方程

k-ω模型由Wilcox提出，用比耗散率ω代替ε。ω的定义为：

$$\omega =\frac{\epsilon }{C_{\mu }k}=\frac{\epsilon }{0.09k}$$

k方程：

$$\frac{\partial k}{\partial t}+{\bar{u}}_j\frac{\partial k}{\partial x_j}=P_k-{\beta }^{\ast }k\omega +\frac{\partial }{\partial x_j}\left[\left(\nu +{\sigma }_k\frac{k}{\omega }\right)\frac{\partial k}{\partial x_j}\right]$$

ω方程：

$$\frac{\partial \omega }{\partial t}+{\bar{u}}_j\frac{\partial \omega }{\partial x_j}=\alpha \frac{\omega }{k}{P}_k-\beta {\omega }^2+\frac{\partial }{\partial x_j}\left[\left(\nu +{\sigma }_{\omega }\frac{k}{\omega }\right)\frac{\partial \omega }{\partial x_j}\right]$$

涡粘系数：

$${\nu }_t=\frac{k}{\omega }$$

8.2 模型常数

Wilcox k-ω模型的常数：

常数	数值
$\alpha$	5/9
$\beta$	3/40
${\beta }^{\ast }$	9/100
${\sigma }_k$	1/2
${\sigma }_{\omega }$	1/2

8.3 壁面边界条件

k-ω模型的一个优势是可以直接应用到壁面，无需壁面函数：

壁面处：

$$k=0,\omega =\frac{6\nu }{\beta y^2}$$

或采用渐进边界条件：

$$\omega =\frac{u_{\tau }^2}{\nu }\text{或}\omega =\frac{4000u_{\tau }}{y_1^2}$$

8.4 k-ω模型的特点

优点：

• 对逆压梯度流动和分离流预测更好
• 可以积分到壁面，无需壁面函数
• 对低雷诺数流动适用
• 对尾流和剪切层预测较好

缺点：

• 对自由流ω值敏感
• 在k-ω和k-ε模型转换区可能不稳定
• 对圆形射流预测仍不理想

9. SST k-ω模型

9.1 混合函数

SST（Shear Stress Transport）k-ω模型由Menter提出，结合了k-ω模型在近壁区的优势和k-ε模型在远场的优势。通过混合函数实现两种模型的平滑过渡。

9.2 模型方程

k方程：

$$\frac{\partial k}{\partial t}+{\bar{u}}_j\frac{\partial k}{\partial x_j}=P_k-{\beta }^{\ast }k\omega +\frac{\partial }{\partial x_j}\left[\left(\nu +{\sigma }_k{\nu }_t\right)\frac{\partial k}{\partial x_j}\right]$$

ω方程：

$$\frac{\partial \omega }{\partial t}+{\bar{u}}_j\frac{\partial \omega }{\partial x_j}=\frac{\alpha }{{\nu }_t}{P}_k-\beta {\omega }^2+\frac{\partial }{\partial x_j}\left[\left(\nu +{\sigma }_{\omega }{\nu }_t\right)\frac{\partial \omega }{\partial x_j}\right]+2(1-F_1)\frac{{\sigma }_{\omega 2}}{\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j}$$

涡粘系数：

$${\nu }_t=\frac{a_{1}k}{\max (a_1\omega ,S{F}_2)}$$

其中，$S=\sqrt{2S_{ij}{S}_{ij}}$是应变率大小。

9.3 混合函数

内层函数：

$$F_1=\tanh \left({\left\{\min \left[\max \left(\frac{\sqrt{k}}{{\beta }^{\ast }\omega y},\frac{500\nu }{y^2\omega }\right),\frac{4{\sigma }_{\omega 2}k}{C{D}_{k\omega }{y}^2}\right]\right\}}^4\right)$$

其中，$C{D}_{k\omega }=\max \left(2{\sigma }_{\omega 2}\frac{1}{\omega }\frac{\partial k}{\partial x_j}\frac{\partial \omega }{\partial x_j},10^{-10}\right)$。

外层函数：

$$F_2=\tanh \left({\left[\max \left(\frac{2\sqrt{k}}{{\beta }^{\ast }\omega y},\frac{500\nu }{y^2\omega }\right)\right]}^2\right)$$

9.4 模型常数

模型常数通过混合函数加权：

$$\varphi =F_1{\varphi }_1+(1-F_1){\varphi }_2$$

其中，${\varphi }_1$是k-ω模型的常数，${\varphi }_2$是k-ε模型的常数。

常数	内层（1）	外层（2）
$\alpha$	5/9	0.44
$\beta$	3/40	0.0828
${\beta }^{\ast }$	9/100	9/100
${\sigma }_k$	0.85	1.0
${\sigma }_{\omega }$	0.5	0.856

9.5 剪切应力限制器

SST模型引入了Bradshaw假设，限制剪切应力与湍动能的关系：

$${\tau }_{xy}=\rho a_{1}k,a_1=0.31$$

这通过修改涡粘系数实现：

$${\nu }_t=\frac{a_{1}k}{\max (a_1\omega ,S{F}_2)}$$

9.6 SST k-ω模型的特点

优点：

• 对逆压梯度流动和分离流预测优秀
• 可以积分到壁面
• 对翼型和航空应用特别成功
• 对自由流值不敏感
• 在广泛流动中表现一致

缺点：

• 计算成本略高于标准k-ε模型
• 需要计算到壁面的网格
• 对复杂三维流动可能需要验证

10. 雷诺应力模型（RSM）

10.1 模型思想

雷诺应力模型（Reynolds Stress Model，RSM）也称为二阶矩封闭模型，直接求解雷诺应力的输运方程，而不是基于涡粘假设。这种方法理论上更精确，能够捕捉湍流的各向异性效应。

10.2 雷诺应力输运方程

雷诺应力的精确输运方程为：

$$\frac{\partial \overline{u_i^{\prime }{u}_j^{\prime }}}{\partial t}+{\bar{u}}_k\frac{\partial \overline{u_i^{\prime }{u}_j^{\prime }}}{\partial x_k}=P_{ij}+D_{ij}-{\epsilon }_{ij}+{\Pi }_{ij}+{\Omega }_{ij}$$

各项的物理意义：

生成项：

$$P_{ij}=-\left(\overline{u_i^{\prime }{u}_k^{\prime }}\frac{\partial {\bar{u}}_j}{\partial x_k}+\overline{u_j^{\prime }{u}_k^{\prime }}\frac{\partial {\bar{u}}_i}{\partial x_k}\right)$$

扩散项：

$$D_{ij}=-\frac{\partial }{\partial x_k}\left(\overline{u_i^{\prime }{u}_j^{\prime }{u}_k^{\prime }}+\frac{1}{\rho }\overline{p^{\prime }(u_i^{\prime }{\delta }_{jk}+u_j^{\prime }{\delta }_{ik})}-\nu \frac{\partial \overline{u_i^{\prime }{u}_j^{\prime }}}{\partial x_k}\right)$$

耗散项：

$${\epsilon }_{ij}=2\nu \overline{\frac{\partial u_i^{\prime }}{\partial x_k}\frac{\partial u_j^{\prime }}{\partial x_k}}$$

压力-应变项：

$${\Pi }_{ij}=\frac{1}{\rho }\overline{p^{\prime }\left(\frac{\partial u_i^{\prime }}{\partial x_j}+\frac{\partial u_j^{\prime }}{\partial x_i}\right)}$$

旋转项：

$${\Omega }_{ij}=-2{\Omega }_k\left(\overline{u_j^{\prime }{u}_m^{\prime }}{e}_{ikm}+\overline{u_i^{\prime }{u}_m^{\prime }}{e}_{jkm}\right)$$

10.3 封闭模型

RSM需要模型化扩散项、耗散项和压力-应变项。

扩散项模型化：

$$D_{ij}=\frac{\partial }{\partial x_k}\left(C_s\frac{k}{\epsilon }\overline{u_k^{\prime }{u}_l^{\prime }}\frac{\partial \overline{u_i^{\prime }{u}_j^{\prime }}}{\partial x_l}+\nu \frac{\partial \overline{u_i^{\prime }{u}_j^{\prime }}}{\partial x_k}\right)$$

耗散项模型化：

对于高雷诺数流动：

$${\epsilon }_{ij}=\frac{2}{3}\epsilon {\delta }_{ij}$$

压力-应变项模型化（LRR模型）：

$${\Pi }_{ij}=-C_1\frac{\epsilon }{k}\left(\overline{u_i^{\prime }{u}_j^{\prime }}-\frac{2}{3}k{\delta }_{ij}\right)-C_2(P_{ij}-\frac{2}{3}{P}_k{\delta }_{ij})$$

其中，$C_1=1.8$，$C_2=0.6$。

10.4 尺度方程

RSM还需要求解湍动能和耗散率的输运方程来确定湍流时间尺度：

$$\frac{\partial k}{\partial t}+{\bar{u}}_j\frac{\partial k}{\partial x_j}=\frac{1}{2}{P}_{ii}-\epsilon +\frac{\partial }{\partial x_j}\left(C_k\frac{k}{\epsilon }\overline{u_j^{\prime }{u}_k^{\prime }}\frac{\partial k}{\partial x_k}\right)$$

$$\frac{\partial \epsilon }{\partial t}+{\bar{u}}_j\frac{\partial \epsilon }{\partial x_j}=C_{\epsilon 1}\frac{P_k\epsilon }{k}-C_{\epsilon 2}\frac{{\epsilon }^2}{k}+\frac{\partial }{\partial x_j}\left(C_{\epsilon }\frac{k}{\epsilon }\overline{u_j^{\prime }{u}_k^{\prime }}\frac{\partial \epsilon }{\partial x_k}\right)$$

10.5 雷诺应力模型的特点

优点：

• 自动考虑湍流各向异性
• 对强旋流、曲率效应、浮力效应预测更好
• 理论上更严格，适用范围更广
• 对二次流和复杂三维流动预测改进

缺点：

• 计算成本高（7个额外的输运方程）
• 数值稳定性较差，收敛困难
• 模型常数多，标定复杂
• 对简单流动优势不明显

11. 非线性涡粘模型

11.1 模型思想

非线性涡粘模型（Nonlinear Eddy Viscosity Models，NLEVM）在Boussinesq线性涡粘假设基础上，引入应变率和旋转张量的高阶项，以捕捉湍流的各向异性效应。

11.2 应力-应变关系

非线性涡粘模型的一般形式为：

$$\overline{u_i^{\prime }{u}_j^{\prime }}=\frac{2}{3}k{\delta }_{ij}-2{\nu }_t{S}_{ij}+\sum _{n=1}^N{f}_n(S_{ij},{\Omega }_{ij})$$

其中，$S_{ij}=\frac{1}{2}(\frac{\partial {\bar{u}}_i}{\partial x_j}+\frac{\partial {\bar{u}}_j}{\partial x_i})$是应变率张量，${\Omega }_{ij}=\frac{1}{2}(\frac{\partial {\bar{u}}_i}{\partial x_j}-\frac{\partial {\bar{u}}_j}{\partial x_i})$是旋转张量。

11.3 二次模型

Speziale提出的二次非线性模型：

$$\overline{u_i^{\prime }{u}_j^{\prime }}=\frac{2}{3}k{\delta }_{ij}-2{\nu }_t{S}_{ij}+C_1\frac{k^3}{{\epsilon }^2}(S_{ik}{S}_{kj}-\frac{1}{3}{S}_{kl}{S}_{kl}{\delta }_{ij})+C_2\frac{k^3}{{\epsilon }^2}(S_{ik}{\Omega }_{kj}+S_{jk}{\Omega }_{ki})$$

11.4 三次模型

Craft等人提出的三次模型包含更多非线性项：

$$\overline{u_i^{\prime }{u}_j^{\prime }}=\frac{2}{3}k{\delta }_{ij}-2{\nu }_t{S}_{ij}+\sum _{n=1}^5{C}_n{T}_{ij}^{(n)}$$

其中，$T_{ij}^{(n)}$是应变率和旋转张量的各种组合。

11.5 显式代数应力模型（EASM）

EASM将雷诺应力模型中的对流和扩散项模型化，得到显式的应力-应变关系：

$$\overline{u_i^{\prime }{u}_j^{\prime }}=k\left(\frac{2}{3}{\delta }_{ij}+\sum _{\lambda }\frac{{\alpha }_{\lambda }}{{\beta }_{\lambda }}{T}_{ij}^{(\lambda )}\right)$$

其中，系数${\alpha }_{\lambda }$和${\beta }_{\lambda }$是应变率和旋转张量不变量的函数。

11.6 非线性涡粘模型的特点

优点：

• 计算成本低于RSM
• 能捕捉湍流各向异性的主要效应
• 对二次流和分离流预测改进
• 数值稳定性比RSM好

缺点：

• 模型复杂，常数标定困难
• 对某些流动的预测仍不如RSM
• 理论严格性不如RSM

12. RANS模型在结构风工程中的应用

12.1 建筑物风荷载分析

RANS模型广泛用于计算建筑物表面的风压分布。对于规则形状的建筑物，稳态RANS计算可以提供平均风压系数，用于结构设计。

应用流程：

1. 建立几何模型和计算域
2. 生成合适的计算网格（重点关注壁面和分离区）
3. 选择湍流模型（通常SST k-ω或Realizable k-ε）
4. 设置边界条件（入口速度剖面、出口压力、壁面无滑移）
5. 进行稳态计算直至收敛
6. 提取表面风压系数分布
7. 计算等效静力风荷载

注意事项：

• RANS模型主要提供平均风压，脉动风压需要结合湍流强度估计
• 对强分离流动（如钝体绕流），RANS可能低估脉动分量
• 对于重要工程，建议结合风洞试验或LES验证

12.2 风环境评估

在城市规划和建筑设计中，需要评估行人高度的风环境舒适性。RANS模型可以计算建筑群周围的平均风速分布。

评估指标：

• 行人高度（1.5-2.0m）的风速比
• 风环境舒适度等级
• 风加速区域识别

湍流模型选择：

• 标准k-ε模型：计算快，适合初步评估
• RNG k-ε模型：对复杂建筑群预测更好
• SST k-ω模型：对分离流和尾流区预测更准确

12.3 桥梁风荷载

大跨度桥梁的风荷载分析是RANS模型的重要应用。主要关注：

• 主梁断面表面风压分布
• 静力三分力系数（阻力、升力、力矩）
• 桥面风环境

计算策略：

• 采用二维断面计算获取静力系数
• 结合抖振理论计算动力响应
• 对于颤振分析，需要结合CFD与结构动力学

12.4 冷却塔和烟囱

大型冷却塔和烟囱等高耸结构的风荷载分析需要考虑：

• 结构表面的风压分布
• 尾流区的涡脱落效应
• 群体效应（多塔干扰）

RANS模型的局限性：

• 稳态RANS无法捕捉涡脱落的非定常特性
• 对脉动风压的预测不足
• 建议采用非定常RANS（URANS）或DES/LES方法

13. RANS模型的选择指南

13.1 模型选择原则

选择合适的RANS模型需要考虑以下因素：

1. 流动特性：

   • 是否存在分离流
   • 逆压梯度的大小
   • 旋转或曲率效应
   • 自由剪切层还是壁面约束流

2. 计算资源：

   • 网格规模
   • 计算时间限制
   • 精度要求

3. 验证数据：

   • 是否有实验数据验证
   • 模型在类似流动中的应用经验

13.2 不同流动的模型推荐

流动类型	推荐模型	备选模型
简单边界层	标准k-ε	SST k-ω
逆压梯度分离流	SST k-ω	Realizable k-ε
旋转机械	SST k-ω	RNG k-ε
强旋流/曲率流	RSM	SST k-ω
射流流动	Realizable k-ε	RNG k-ε
复杂三维流动	SST k-ω	RSM
航空翼型	SST k-ω	k-ω

13.3 工程实践建议

1. 从简单模型开始：

   • 先用标准k-ε模型进行初步计算
   • 评估流动特性，识别关键区域

2. 网格敏感性分析：

   • 检查近壁面y+值
   • 确保分离区和尾流区网格足够密

3. 模型验证：

   • 与实验数据对比
   • 评估不同模型的预测差异

4. 结合多种方法：

   • RANS用于初步设计和参数研究
   • LES/DES用于详细分析和验证
   • 风洞试验用于最终确认

14. RANS模型的局限性与发展趋势

14.1 主要局限性

1. 无法捕捉大尺度非定常结构：

   • 涡脱落的周期性
   • 大尺度相干结构
   • 湍流的瞬态特性

2. 对分离流的预测不足：

   • 分离点位置偏差
   • 再附长度预测不准
   • 脉动分量低估

3. 模型常数的普适性问题：

   • 不同流动需要调整常数
   • 缺乏严格的物理基础

4. 各向同性假设的局限：

   • 涡粘假设忽略了湍流各向异性
   • 对强旋流和曲率效应预测不足

14.2 改进方向

1. 混合RANS-LES方法：

   • DES（Detached Eddy Simulation）
   • DDES（Delayed DES）
   • IDDES（Improved DDES）
   • 在壁面区用RANS，在分离区用LES

2. 非定常RANS（URANS）：

   • 捕捉周期性流动特征
   • 计算成本低于LES
   • 对涡脱落流动有一定改善

3. 机器学习辅助模型：

   • 数据驱动的湍流模型
   • 基于DNS/LES数据的模型训练
   • 提高模型的适应性和精度

4. 壁面模型LES：

   • 减少近壁面网格需求
   • 结合RANS壁面模型与LES外区

15. Python仿真案例

15.1 案例一：平板边界层流动模拟

本案例使用简化的RANS方法模拟平板边界层发展，比较不同湍流模型的效果。

"""
平板边界层流动的RANS模拟
使用简化的涡粘模型求解边界层发展
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import matplotlib
matplotlib.use('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

class BoundaryLayerRANS:
    """
    平板边界层RANS求解器
    使用积分方法求解边界层发展
    """
    
    def __init__(self, U_inf, nu, L, Nx=200):
        """
        初始化求解器
        
        参数:
            U_inf: 自由流速度 [m/s]
            nu: 运动粘度 [m^2/s]
            L: 平板长度 [m]
            Nx: 流向网格数
        """
        self.U_inf = U_inf
        self.nu = nu
        self.L = L
        self.Nx = Nx
        
        # 流向网格
        self.x = np.linspace(0.001, L, Nx)  # 从0.001开始避免奇点
        
        # 雷诺数
        self.Re_L = U_inf * L / nu
        
    def blasius_solution(self, eta_max=10, N=100):
        """
        Blasius方程数值解（层流边界层）
        
        参数:
            eta_max: 相似变量最大值
            N: 网格数
        
        返回:
            eta, f, f', f''
        """
        eta = np.linspace(0, eta_max, N)
        
        # Blasius方程: f''' + 0.5*f*f'' = 0
        # 转化为方程组:
        # y0' = y1  (f' = u)
        # y1' = y2  (f'' = du/deta)
        # y2' = -0.5*y0*y2
        
        def blasius_ode(y, eta):
            return [y[1], y[2], -0.5*y[0]*y[2]]
        
        # 边界条件: f(0)=0, f'(0)=0, f'(inf)=1
        # 使用打靶法求解 f''(0)
        
        def shoot(s):
            y0 = [0, 0, s]
            sol = odeint(blasius_ode, y0, eta)
            return sol[-1, 1] - 1  # f'(inf) - 1 = 0
        
        # 二分法求解
        s_low, s_high = 0.1, 1.0
        while abs(s_high - s_low) > 1e-6:
            s_mid = (s_low + s_high) / 2
            if shoot(s_low) * shoot(s_mid) < 0:
                s_high = s_mid
            else:
                s_low = s_mid
        
        s_optimal = (s_low + s_high) / 2
        y0 = [0, 0, s_optimal]
        sol = odeint(blasius_ode, y0, eta)
        
        return eta, sol[:, 0], sol[:, 1], sol[:, 2]
    
    def turbulent_velocity_profile(self, y_delta, Re_theta, model='mixing_length'):
        """
        湍流边界层速度剖面
        
        参数:
            y_delta: 无量纲高度 y/δ
            Re_theta: 动量厚度雷诺数
            model: 湍流模型 ('mixing_length', 'k_epsilon')
        
        返回:
            u_Uinf: 无量纲速度 u/U_inf
        """
        if model == 'mixing_length':
            # 混合长度模型 - 使用壁面律
            # 内层: u+ = y+ (粘性底层)
            # 外层: u+ = (1/κ)ln(y+) + B (对数律)
            
            kappa = 0.41
            B = 5.0
            
            # 假设壁面摩擦速度 u_tau
            # 对于零压力梯度边界层: Cf ≈ 0.026/Re_theta^0.25
            Cf = 0.026 * Re_theta**(-0.25)
            u_tau = np.sqrt(Cf / 2) * self.U_inf
            
            # 边界层厚度估计
            delta = self.nu * Re_theta / u_tau
            
            y = y_delta * delta
            y_plus = y * u_tau / self.nu
            
            u_plus = np.zeros_like(y_plus)
            
            # 粘性底层
            mask_viscous = y_plus < 11.6
            u_plus[mask_viscous] = y_plus[mask_viscous]
            
            # 对数律层
            mask_log = y_plus >= 11.6
            u_plus[mask_log] = (1/kappa) * np.log(y_plus[mask_log]) + B
            
            # 尾流律修正
            # Coles尾流函数
            Pi = 0.55  # 尾流强度参数
            wake = (2/kappa) * np.sin(np.pi/2 * y_delta)**2
            u_plus = u_plus + Pi * wake
            
            # 限制最大速度
            u_plus = np.minimum(u_plus, self.U_inf / u_tau)
            
            return u_plus * u_tau / self.U_inf
        
        elif model == 'power_law':
            # 幂律速度剖面
            n = 7  # 幂律指数 (随Re变化)
            if Re_theta > 10000:
                n = 8
            elif Re_theta > 100000:
                n = 9
            
            return y_delta**(1/n)
        
        else:
            raise ValueError(f"未知模型: {model}")
    
    def calculate_boundary_layer_params(self, model='mixing_length'):
        """
        计算边界层参数沿平板的发展
        
        参数:
            model: 湍流模型
        
        返回:
            delta, delta_star, theta, Cf: 边界层厚度、位移厚度、动量厚度、摩擦系数
        """
        delta = np.zeros(self.Nx)
        delta_star = np.zeros(self.Nx)
        theta = np.zeros(self.Nx)
        Cf = np.zeros(self.Nx)
        
        for i, x in enumerate(self.x):
            Re_x = self.U_inf * x / self.nu
            
            if Re_x < 5e5:  # 层流
                # Blasius解
                delta[i] = 5.0 * x / np.sqrt(Re_x)
                delta_star[i] = 1.72 * x / np.sqrt(Re_x)
                theta[i] = 0.664 * x / np.sqrt(Re_x)
                Cf[i] = 0.664 / np.sqrt(Re_x)
            else:  # 湍流
                # 经验公式 (基于1/7幂律)
                delta[i] = 0.37 * x / (Re_x**0.2)
                delta_star[i] = 0.046 * x / (Re_x**0.2)
                theta[i] = 0.036 * x / (Re_x**0.2)
                Cf[i] = 0.026 / (Re_x**0.2)
        
        return delta, delta_star, theta, Cf
    
    def compare_velocity_profiles(self):
        """
        比较不同湍流模型的速度剖面
        """
        fig, axes = plt.subplots(2, 2, figsize=(14, 12))
        
        # 1. 层流Blasius解
        ax1 = axes[0, 0]
        eta, f, fp, fpp = self.blasius_solution()
        ax1.plot(fp, eta, 'b-', linewidth=2, label='Blasius解')
        ax1.set_xlabel('u/U_∞', fontsize=11)
        ax1.set_ylabel('η = y√(U_∞/νx)', fontsize=11)
        ax1.set_title('层流边界层速度剖面 (Blasius解)', fontsize=12)
        ax1.grid(True, alpha=0.3)
        ax1.legend()
        
        # 2. 湍流速度剖面比较
        ax2 = axes[0, 1]
        y_delta = np.linspace(0.001, 1, 100)
        
        Re_theta_values = [1000, 5000, 10000]
        colors = ['blue', 'green', 'red']
        
        for Re_theta, color in zip(Re_theta_values, colors):
            u = self.turbulent_velocity_profile(y_delta, Re_theta, 'mixing_length')
            ax2.plot(u, y_delta, color=color, linewidth=2, 
                    label=f'Re_θ = {Re_theta}')
        
        # 幂律剖面
        u_power = self.turbulent_velocity_profile(y_delta, 5000, 'power_law')
        ax2.plot(u_power, y_delta, 'k--', linewidth=2, label='1/7幂律')
        
        ax2.set_xlabel('u/U_∞', fontsize=11)
        ax2.set_ylabel('y/δ', fontsize=11)
        ax2.set_title('湍流边界层速度剖面比较', fontsize=12)
        ax2.grid(True, alpha=0.3)
        ax2.legend()
        
        # 3. 边界层厚度发展
        ax3 = axes[1, 0]
        delta, delta_star, theta, Cf = self.calculate_boundary_layer_params()
        
        ax3.plot(self.x, delta, 'b-', linewidth=2, label='δ (边界层厚度)')
        ax3.plot(self.x, delta_star, 'r--', linewidth=2, label='δ* (位移厚度)')
        ax3.plot(self.x, theta, 'g:', linewidth=2, label='θ (动量厚度)')
        
        # 标记转捩点
        Re_trans = 5e5
        x_trans = Re_trans * self.nu / self.U_inf
        ax3.axvline(x=x_trans, color='k', linestyle='-.', alpha=0.5, label='转捩点')
        
        ax3.set_xlabel('x [m]', fontsize=11)
        ax3.set_ylabel('厚度 [m]', fontsize=11)
        ax3.set_title('边界层厚度沿平板发展', fontsize=12)
        ax3.grid(True, alpha=0.3)
        ax3.legend()
        
        # 4. 壁面摩擦系数
        ax4 = axes[1, 1]
        
        # 层流理论
        Re_x_lam = np.linspace(1e4, 5e5, 100)
        Cf_lam = 0.664 / np.sqrt(Re_x_lam)
        
        # 湍流经验公式
        Re_x_turb = np.linspace(5e5, self.Re_L, 100)
        Cf_turb = 0.026 / (Re_x_turb**0.2)
        
        ax4.semilogy(Re_x_lam, Cf_lam, 'b-', linewidth=2, label='层流 (Blasius)')
        ax4.semilogy(Re_x_turb, Cf_turb, 'r-', linewidth=2, label='湍流 (经验公式)')
        
        # 计算结果
        Re_x_calc = self.U_inf * self.x / self.nu
        ax4.semilogy(Re_x_calc, Cf, 'g--', linewidth=2, alpha=0.7, label='混合计算')
        
        ax4.axvline(x=5e5, color='k', linestyle='-.', alpha=0.5)
        ax4.set_xlabel('Re_x', fontsize=11)
        ax4.set_ylabel('C_f', fontsize=11)
        ax4.set_title('壁面摩擦系数', fontsize=12)
        ax4.grid(True, alpha=0.3)
        ax4.legend()
        
        plt.tight_layout()
        plt.savefig('boundary_layer_rans.png', dpi=150, bbox_inches='tight')
        plt.close()
        
        print("边界层RANS分析图已保存为 'boundary_layer_rans.png'")
    
    def analyze_turbulence_models(self):
        """
        分析不同湍流模型的边界层特性
        """
        fig, axes = plt.subplots(2, 2, figsize=(14, 12))
        
        # 计算位置
        x_eval = self.L / 2
        Re_x = self.U_inf * x_eval / self.nu
        
        # 1. 不同模型的速度剖面
        ax1 = axes[0, 0]
        y_delta = np.linspace(0.001, 1, 100)
        
        # 混合长度模型
        Re_theta = 5000
        u_ml = self.turbulent_velocity_profile(y_delta, Re_theta, 'mixing_length')
        ax1.plot(u_ml, y_delta, 'b-', linewidth=2, label='混合长度模型')
        
        # 幂律模型
        u_pl = self.turbulent_velocity_profile(y_delta, Re_theta, 'power_law')
        ax1.plot(u_pl, y_delta, 'r--', linewidth=2, label='1/7幂律')
        
        ax1.set_xlabel('u/U_∞', fontsize=11)
        ax1.set_ylabel('y/δ', fontsize=11)
        ax1.set_title(f'湍流模型速度剖面比较 (Re_x={Re_x:.2e})', fontsize=12)
        ax1.grid(True, alpha=0.3)
        ax1.legend()
        
        # 2. 形状因子比较
        ax2 = axes[0, 1]
        
        Re_range = np.logspace(4, 7, 100)
        
        # 层流形状因子
        H_lam = np.ones_like(Re_range) * 2.59  # Blasius解
        
        # 湍流形状因子 (经验公式)
        H_turb = 1.3 + 0.5 * (Re_range/1e6)**(-0.1)
        H_turb = np.clip(H_turb, 1.28, 1.6)
        
        ax2.semilogx(Re_range, H_lam, 'b-', linewidth=2, label='层流 (H=2.59)')
        ax2.semilogx(Re_range, H_turb, 'r-', linewidth=2, label='湍流 (经验)')
        ax2.axvline(x=5e5, color='k', linestyle='-.', alpha=0.5, label='转捩')
        
        ax2.set_xlabel('Re_x', fontsize=11)
        ax2.set_ylabel('形状因子 H', fontsize=11)
        ax2.set_title('边界层形状因子', fontsize=12)
        ax2.grid(True, alpha=0.3)
        ax2.legend()
        ax2.set_ylim([1.2, 3.0])
        
        # 3. 动量厚度雷诺数发展
        ax3 = axes[1, 0]
        
        _, _, theta, _ = self.calculate_boundary_layer_params()
        Re_theta_calc = self.U_inf * theta / self.nu
        
        ax3.semilogy(self.x, Re_theta_calc, 'b-', linewidth=2)
        ax3.axvline(x=5e5*self.nu/self.U_inf, color='k', linestyle='-.', alpha=0.5, label='转捩点')
        
        ax3.set_xlabel('x [m]', fontsize=11)
        ax3.set_ylabel('Re_θ', fontsize=11)
        ax3.set_title('动量厚度雷诺数发展', fontsize=12)
        ax3.grid(True, alpha=0.3)
        ax3.legend()
        
        # 4. 湍流模型特性对比
        ax4 = axes[1, 1]
        
        models = ['零方程\\n(混合长度)', '标准k-ε', 'RNG k-ε', 'Realizable\\nk-ε', 'SST k-ω', 'RSM']
        accuracy = [3, 6, 7, 7.5, 8.5, 9]  # 分离流预测精度
        cost = [2, 4, 4.5, 5, 5.5, 9]  # 计算成本
        
        x_pos = np.arange(len(models))
        width = 0.35
        
        bars1 = ax4.bar(x_pos - width/2, accuracy, width, label='分离流预测精度', color='steelblue')
        bars2 = ax4.bar(x_pos + width/2, cost, width, label='计算成本', color='coral')
        
        ax4.set_ylabel('评分 (1-10)', fontsize=11)
        ax4.set_title('湍流模型特性对比', fontsize=12)
        ax4.set_xticks(x_pos)
        ax4.set_xticklabels(models, fontsize=9)
        ax4.legend()
        ax4.grid(True, alpha=0.3, axis='y')
        
        # 添加数值标签
        for bar in bars1:
            height = bar.get_height()
            ax4.annotate(f'{height:.1f}',
                        xy=(bar.get_x() + bar.get_width() / 2, height),
                        xytext=(0, 3),
                        textcoords="offset points",
                        ha='center', va='bottom', fontsize=8)
        
        for bar in bars2:
            height = bar.get_height()
            ax4.annotate(f'{height:.1f}',
                        xy=(bar.get_x() + bar.get_width() / 2, height),
                        xytext=(0, 3),
                        textcoords="offset points",
                        ha='center', va='bottom', fontsize=8)
        
        plt.tight_layout()
        plt.savefig('turbulence_model_analysis.png', dpi=150, bbox_inches='tight')
        plt.close()
        
        print("湍流模型分析图已保存为 'turbulence_model_analysis.png'")


# 运行分析
if __name__ == "__main__":
    # 创建求解器实例
    solver = BoundaryLayerRANS(U_inf=10.0, nu=1.5e-5, L=5.0, Nx=200)
    
    print("=" * 60)
    print("平板边界层RANS模拟分析")
    print("=" * 60)
    print(f"自由流速度: {solver.U_inf} m/s")
    print(f"运动粘度: {solver.nu} m²/s")
    print(f"平板长度: {solver.L} m")
    print(f"雷诺数: Re_L = {solver.Re_L:.2e}")
    print("=" * 60)
    
    # 边界层分析
    solver.compare_velocity_profiles()
    
    # 湍流模型分析
    solver.analyze_turbulence_models()
    
    print("\n分析完成!")

15.2 案例二：方柱绕流RANS模拟

本案例使用简化的RANS方法模拟方柱绕流，比较不同湍流模型对分离流和阻力的预测。

"""
方柱绕流的简化RANS模拟
使用涡粘模型模拟湍流效应
"""

import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve
import matplotlib
matplotlib.use('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

class SquareCylinderRANS:
    """
    方柱绕流RANS求解器（简化模型）
    基于势流理论修正+涡粘模型
    """
    
    def __init__(self, D=1.0, U_inf=1.0, Re=1e4, turbulence_model='k_epsilon'):
        """
        初始化求解器
        
        参数:
            D: 方柱边长
            U_inf: 来流速度
            Re: 雷诺数
            turbulence_model: 湍流模型类型
        """
        self.D = D
        self.U_inf = U_inf
        self.Re = Re
        self.turbulence_model = turbulence_model
        
        # 计算域
        self.Lx = 20 * D
        self.Ly = 10 * D
        
        # 网格参数
        self.Nx = 200
        self.Ny = 100
        
        self.dx = self.Lx / (self.Nx - 1)
        self.dy = self.Ly / (self.Ny - 1)
        
        # 方柱位置
        self.cylinder_x = 5 * D
        self.cylinder_y = self.Ly / 2
        
        # 初始化场
        self.u = np.ones((self.Nx, self.Ny)) * U_inf
        self.v = np.zeros((self.Nx, self.Ny))
        self.p = np.zeros((self.Nx, self.Ny))
        self.k = np.ones((self.Nx, self.Ny)) * 0.001 * U_inf**2
        self.epsilon = np.ones((self.Nx, self.Ny)) * 0.001 * U_inf**3 / D
        
        # 标记方柱区域
        self.create_cylinder_mask()
        
    def create_cylinder_mask(self):
        """
        创建方柱掩码
        """
        self.cylinder_mask = np.zeros((self.Nx, self.Ny), dtype=bool)
        
        for i in range(self.Nx):
            for j in range(self.Ny):
                x = i * self.dx
                y = j * self.dy
                
                if (abs(x - self.cylinder_x) <= self.D/2 and 
                    abs(y - self.cylinder_y) <= self.D/2):
                    self.cylinder_mask[i, j] = True
        
        # 在方柱内部设置速度为零
        self.u[self.cylinder_mask] = 0
        self.v[self.cylinder_mask] = 0
    
    def calculate_nu_t(self):
        """
        计算湍流涡粘系数
        """
        if self.turbulence_model == 'k_epsilon':
            # 标准k-ε模型
            C_mu = 0.09
            nu_t = C_mu * self.k**2 / (self.epsilon + 1e-10)
        elif self.turbulence_model == 'mixing_length':
            # 混合长度模型
            # 计算涡量
            dvdx = np.gradient(self.v, self.dx, axis=0)
            dudy = np.gradient(self.u, self.dy, axis=1)
            omega = np.abs(dudy - dvdx)
            
            # 混合长度
            lm = 0.09 * self.D  # 简化处理
            nu_t = lm**2 * omega
        else:
            # 常数涡粘模型
            nu_t = np.ones_like(self.u) * 0.01 * self.U_inf * self.D
        
        return np.clip(nu_t, 0, 10 * self.U_inf * self.D)
    
    def solve_potential_flow(self, iterations=1000):
        """
        求解势流场（简化模型）
        使用迭代方法求解速度场
        """
        nu_t = self.calculate_nu_t()
        nu_eff = 1.5e-5 / (self.U_inf * self.D) + nu_t  # 无量纲有效粘度
        
        for _ in range(iterations):
            u_old = self.u.copy()
            v_old = self.v.copy()
            
            # 内部点更新（不包括边界和方柱）
            for i in range(1, self.Nx-1):
                for j in range(1, self.Ny-1):
                    if self.cylinder_mask[i, j]:
                        continue
                    
                    # 简化的速度更新（考虑对流和扩散）
                    # 使用迎风格式处理对流项
                    if u_old[i,j] > 0:
                        dudx = (u_old[i,j] - u_old[i-1,j]) / self.dx
                        dvdx = (v_old[i,j] - v_old[i-1,j]) / self.dx
                    else:
                        dudx = (u_old[i+1,j] - u_old[i,j]) / self.dx
                        dvdx = (v_old[i+1,j] - v_old[i,j]) / self.dx
                    
                    if v_old[i,j] > 0:
                        dudy = (u_old[i,j] - u_old[i,j-1]) / self.dy
                        dvdy = (v_old[i,j] - v_old[i,j-1]) / self.dy
                    else:
                        dudy = (u_old[i,j+1] - u_old[i,j]) / self.dy
                        dvdy = (v_old[i,j+1] - v_old[i,j]) / self.dy
                    
                    # 扩散项（Laplacian）
                    d2udx2 = (u_old[i+1,j] - 2*u_old[i,j] + u_old[i-1,j]) / self.dx**2
                    d2udy2 = (u_old[i,j+1] - 2*u_old[i,j] + u_old[i,j-1]) / self.dy**2
                    d2vdx2 = (v_old[i+1,j] - 2*v_old[i,j] + v_old[i-1,j]) / self.dx**2
                    d2vdy2 = (v_old[i,j+1] - 2*v_old[i,j] + v_old[i,j-1]) / self.dy**2
                    
                    # 更新速度（松弛迭代）
                    omega_relax = 0.1
                    self.u[i,j] = (1-omega_relax) * u_old[i,j] + omega_relax * (
                        u_old[i,j] - 0.1 * (u_old[i,j]*dudx + v_old[i,j]*dudy - 
                        nu_eff[i,j] * (d2udx2 + d2udy2))
                    )
                    self.v[i,j] = (1-omega_relax) * v_old[i,j] + omega_relax * (
                        v_old[i,j] - 0.1 * (u_old[i,j]*dvdx + v_old[i,j]*dvdy - 
                        nu_eff[i,j] * (d2vdx2 + d2vdy2))
                    )
            
            # 应用边界条件
            self.apply_boundary_conditions()
    
    def apply_boundary_conditions(self):
        """
        应用边界条件
        """
        # 入口：均匀来流
        self.u[0, :] = self.U_inf
        self.v[0, :] = 0
        
        # 出口：零梯度
        self.u[-1, :] = self.u[-2, :]
        self.v[-1, :] = self.v[-2, :]
        
        # 上下边界：对称或自由滑移
        self.u[:, 0] = self.u[:, 1]
        self.v[:, 0] = 0
        self.u[:, -1] = self.u[:, -2]
        self.v[:, -1] = 0
        
        # 方柱表面：无滑移
        self.u[self.cylinder_mask] = 0
        self.v[self.cylinder_mask] = 0
    
    def calculate_pressure_coefficient(self):
        """
        计算压力系数分布
        """
        # 基于速度计算压力（伯努利方程简化）
        velocity_magnitude = np.sqrt(self.u**2 + self.v**2)
        Cp = 1 - (velocity_magnitude / self.U_inf)**2
        return Cp
    
    def calculate_drag_lift_coefficients(self):
        """
        计算阻力系数和升力系数（简化估算）
        """
        Cp = self.calculate_pressure_coefficient()
        
        # 提取方柱表面的压力系数
        # 前表面
        front_mask = np.zeros_like(self.cylinder_mask)
        i_front = int(self.cylinder_x / self.dx - self.D/2/self.dx)
        j_min = int((self.cylinder_y - self.D/2) / self.dy)
        j_max = int((self.cylinder_y + self.D/2) / self.dy)
        front_mask[i_front, j_min:j_max] = True
        
        # 后表面
        back_mask = np.zeros_like(self.cylinder_mask)
        i_back = int(self.cylinder_x / self.dx + self.D/2/self.dx)
        back_mask[i_back, j_min:j_max] = True
        
        # 计算平均压力系数
        Cp_front = np.mean(Cp[front_mask]) if np.any(front_mask) else 1.0
        Cp_back = np.mean(Cp[back_mask]) if np.any(back_mask) else -0.5
        
        # 阻力系数估算
        Cd = (Cp_front - Cp_back)  # 简化计算
        
        return Cd, 0.0  # 升力系数在此简化模型中为0
    
    def visualize_results(self):
        """
        可视化计算结果
        """
        fig, axes = plt.subplots(2, 3, figsize=(16, 10))
        
        # 创建网格
        X = np.linspace(0, self.Lx, self.Nx)
        Y = np.linspace(0, self.Ly, self.Ny)
        X, Y = np.meshgrid(X, Y)
        
        # 1. 速度矢量图
        ax1 = axes[0, 0]
        skip = 5
        ax1.quiver(X[::skip, ::skip], Y[::skip, ::skip], 
                  self.u[::skip, ::skip].T, self.v[::skip, ::skip].T,
                  scale=20, alpha=0.7)
        
        # 绘制方柱
        rect = plt.Rectangle((self.cylinder_x - self.D/2, self.cylinder_y - self.D/2),
                            self.D, self.D, fill=True, color='gray', alpha=0.8)
        ax1.add_patch(rect)
        
        ax1.set_xlabel('x [m]', fontsize=10)
        ax1.set_ylabel('y [m]', fontsize=10)
        ax1.set_title('速度矢量场', fontsize=11)
        ax1.set_aspect('equal')
        ax1.set_xlim([0, self.Lx])
        ax1.set_ylim([0, self.Ly])
        
        # 2. 流线图
        ax2 = axes[0, 1]
        speed = np.sqrt(self.u.T**2 + self.v.T**2)
        strm = ax2.streamplot(X, Y, self.u.T, self.v.T, 
                             color=speed, cmap='viridis', density=2, linewidth=1)
        plt.colorbar(strm.lines, ax=ax2, label='速度 [m/s]')
        
        rect2 = plt.Rectangle((self.cylinder_x - self.D/2, self.cylinder_y - self.D/2),
                             self.D, self.D, fill=True, color='gray', alpha=0.8)
        ax2.add_patch(rect2)
        
        ax2.set_xlabel('x [m]', fontsize=10)
        ax2.set_ylabel('y [m]', fontsize=10)
        ax2.set_title('流线图', fontsize=11)
        ax2.set_aspect('equal')
        
        # 3. 压力系数分布
        ax3 = axes[0, 2]
        Cp = self.calculate_pressure_coefficient()
        contour = ax3.contourf(X, Y, Cp.T, levels=20, cmap='RdBu_r', vmin=-1.5, vmax=1.5)
        plt.colorbar(contour, ax=ax3, label='Cp')
        
        rect3 = plt.Rectangle((self.cylinder_x - self.D/2, self.cylinder_y - self.D/2),
                             self.D, self.D, fill=True, color='gray', alpha=0.8)
        ax3.add_patch(rect3)
        
        ax3.set_xlabel('x [m]', fontsize=10)
        ax3.set_ylabel('y [m]', fontsize=10)
        ax3.set_title('压力系数分布', fontsize=11)
        ax3.set_aspect('equal')
        
        # 4. 湍动能分布
        ax4 = axes[1, 0]
        k_plot = np.log10(self.k.T + 1e-10)
        contour_k = ax4.contourf(X, Y, k_plot, levels=20, cmap='hot')
        plt.colorbar(contour_k, ax=ax4, label='log₁₀(k)')
        
        rect4 = plt.Rectangle((self.cylinder_x - self.D/2, self.cylinder_y - self.D/2),
                             self.D, self.D, fill=True, color='gray', alpha=0.8)
        ax4.add_patch(rect4)
        
        ax4.set_xlabel('x [m]', fontsize=10)
        ax4.set_ylabel('y [m]', fontsize=10)
        ax4.set_title('湍动能分布 (log尺度)', fontsize=11)
        ax4.set_aspect('equal')
        
        # 5. 涡粘系数分布
        ax5 = axes[1, 1]
        nu_t = self.calculate_nu_t()
        contour_nut = ax5.contourf(X, Y, nu_t.T, levels=20, cmap='plasma')
        plt.colorbar(contour_nut, ax=ax5, label='νt [m²/s]')
        
        rect5 = plt.Rectangle((self.cylinder_x - self.D/2, self.cylinder_y - self.D/2),
                             self.D, self.D, fill=True, color='gray', alpha=0.8)
        ax5.add_patch(rect5)
        
        ax5.set_xlabel('x [m]', fontsize=10)
        ax5.set_ylabel('y [m]', fontsize=10)
        ax5.set_title('涡粘系数分布', fontsize=11)
        ax5.set_aspect('equal')
        
        # 6. 中心线速度分布
        ax6 = axes[1, 2]
        j_center = int(self.cylinder_y / self.dy)
        ax6.plot(X[:, j_center], self.u[:, j_center], 'b-', linewidth=2, label='u/U∞')
        ax6.axvline(x=self.cylinder_x - self.D/2, color='r', linestyle='--', alpha=0.5, label='方柱前缘')
        ax6.axvline(x=self.cylinder_x + self.D/2, color='r', linestyle='--', alpha=0.5, label='方柱后缘')
        ax6.set_xlabel('x [m]', fontsize=10)
        ax6.set_ylabel('u/U∞', fontsize=10)
        ax6.set_title('中心线速度分布', fontsize=11)
        ax6.grid(True, alpha=0.3)
        ax6.legend()
        
        plt.tight_layout()
        plt.savefig('square_cylinder_rans.png', dpi=150, bbox_inches='tight')
        plt.close()
        
        print("方柱绕流RANS分析图已保存为 'square_cylinder_rans.png'")
    
    def compare_turbulence_models(self):
        """
        比较不同湍流模型的结果
        """
        models = ['constant', 'mixing_length', 'k_epsilon']
        model_names = ['常数涡粘', '混合长度', 'k-ε模型']
        
        fig, axes = plt.subplots(2, 3, figsize=(16, 10))
        
        for idx, (model, name) in enumerate(zip(models, model_names)):
            self.turbulence_model = model
            self.u = np.ones((self.Nx, self.Ny)) * self.U_inf
            self.v = np.zeros((self.Nx, self.Ny))
            
            # 求解
            self.solve_potential_flow(iterations=500)
            
            X = np.linspace(0, self.Lx, self.Nx)
            Y = np.linspace(0, self.Ly, self.Ny)
            X, Y = np.meshgrid(X, Y)
            
            # 速度云图
            ax = axes[0, idx]
            speed = np.sqrt(self.u.T**2 + self.v.T**2)
            contour = ax.contourf(X, Y, speed, levels=20, cmap='viridis', vmin=0, vmax=self.U_inf*1.2)
            plt.colorbar(contour, ax=ax, label='速度 [m/s]')
            
            rect = plt.Rectangle((self.cylinder_x - self.D/2, self.cylinder_y - self.D/2),
                                self.D, self.D, fill=True, color='gray', alpha=0.8)
            ax.add_patch(rect)
            
            ax.set_title(f'{name} - 速度场', fontsize=11)
            ax.set_xlabel('x [m]', fontsize=10)
            ax.set_ylabel('y [m]', fontsize=10)
            ax.set_aspect('equal')
            
            # 压力系数
            ax2 = axes[1, idx]
            Cp = self.calculate_pressure_coefficient()
            contour2 = ax2.contourf(X, Y, Cp.T, levels=20, cmap='RdBu_r', vmin=-1.5, vmax=1.5)
            plt.colorbar(contour2, ax=ax2, label='Cp')
            
            rect2 = plt.Rectangle((self.cylinder_x - self.D/2, self.cylinder_y - self.D/2),
                                 self.D, self.D, fill=True, color='gray', alpha=0.8)
            ax2.add_patch(rect2)
            
            ax2.set_title(f'{name} - 压力系数', fontsize=11)
            ax2.set_xlabel('x [m]', fontsize=10)
            ax2.set_ylabel('y [m]', fontsize=10)
            ax2.set_aspect('equal')
        
        plt.tight_layout()
        plt.savefig('turbulence_model_comparison.png', dpi=150, bbox_inches='tight')
        plt.close()
        
        print("湍流模型比较图已保存为 'turbulence_model_comparison.png'")


# 运行分析
if __name__ == "__main__":
    print("=" * 60)
    print("方柱绕流RANS模拟")
    print("=" * 60)
    
    # 创建求解器实例
    solver = SquareCylinderRANS(D=1.0, U_inf=10.0, Re=1e4, turbulence_model='k_epsilon')
    
    print(f"方柱边长: {solver.D} m")
    print(f"来流速度: {solver.U_inf} m/s")
    print(f"雷诺数: Re = {solver.Re:.2e}")
    print("=" * 60)
    
    # 求解流场
    print("求解流场...")
    solver.solve_potential_flow(iterations=1000)
    
    # 计算气动力系数
    Cd, Cl = solver.calculate_drag_lift_coefficients()
    print(f"\n估算阻力系数: Cd ≈ {Cd:.3f}")
    print(f"估算升力系数: Cl ≈ {Cl:.3f}")
    print("(注: 简化模型结果，仅供参考)")
    
    # 可视化结果
    print("\n生成可视化结果...")
    solver.visualize_results()
    
    # 比较不同湍流模型
    print("\n比较不同湍流模型...")
    solver.compare_turbulence_models()
    
    print("\n分析完成!")