主题006：湍流模型与应用

第一章 湍流现象与基本特征

1.1 湍流的定义与本质

湍流是流体力学中最复杂、最迷人的现象之一。从日常生活中的香烟烟雾到大气环流，从管道内的流体输送到飞机机翼周围的气流，湍流无处不在。湍流是一种高度不规则的流体运动状态，其特征是流体微团在时间和空间上呈现随机脉动。与层流的有序流动不同，湍流中充满了各种尺度的涡旋结构，能量从大尺度涡旋向小尺度涡旋传递，最终通过粘性耗散转化为热能。

湍流的本质可以从多个角度理解。从物理角度看，湍流是惯性力与粘性力相互作用的结果。当雷诺数超过临界值时，流动失稳，层流转变为湍流。从数学角度看，湍流是Navier-Stokes方程在特定条件下的解，这种解对初始条件极其敏感，表现出混沌特性。从能量角度看，湍流是一种高效的能量耗散机制，能够在较短的时空尺度内将宏观动能转化为热能。

湍流的研究具有重要的理论意义和实际价值。在理论上，湍流是经典物理学中尚未完全解决的难题之一，涉及非线性动力学、统计物理、分形几何等多个学科。在实际应用中，湍流影响着能源、环境、航空、化工等众多领域。准确预测和控制湍流对于提高能源利用效率、减少环境污染、优化工程设计具有重要意义。

1.2 湍流的基本特征

湍流具有以下基本特征：

随机性：湍流中的速度、压力等物理量随时间和空间随机变化。这种随机性不是完全无规律的，而是具有一定的统计规律。湍流脉动可以用统计方法描述，如概率密度函数、相关函数、谱分析等。

涡旋结构：湍流由各种尺度的涡旋组成。大尺度涡旋从平均流动中获取能量，通过涡旋拉伸和破碎过程将能量传递给较小尺度的涡旋，形成能量级串。在最小尺度上，粘性作用占主导，能量被耗散。这种能量传递过程是湍流的核心特征。

三维性：湍流本质上是三维的。涡旋的拉伸和扭曲需要三维速度梯度的相互作用。即使在几何上二维的流动中，湍流也会发展出三维结构。

扩散性：湍流具有强烈的混合和扩散能力。湍流扩散远强于分子扩散，这是湍流在工程应用中重要的特性。湍流传热和传质效率远高于层流。

耗散性：湍流是一种高效的能量耗散机制。能量从大尺度向小尺度传递，最终在Kolmogorov尺度上通过粘性耗散。这种能量级串过程使得湍流能够在宏观尺度上维持，同时在微观尺度上耗散能量。

间歇性：湍流并非均匀分布，而是呈现间歇性。在某些区域和时刻，湍流强度很高；在其他区域和时刻，流动可能接近层流状态。这种间歇性在湍流的尾缘和边界层外缘尤为明显。

1.3 湍流的统计描述

由于湍流的随机性，通常采用统计方法进行描述。最常用的方法是雷诺分解，将瞬时物理量分解为平均值和脉动值：

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

其中，$\bar{u}$是时间平均速度，$u^{\prime }$是速度脉动。时间平均定义为：

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

在实际应用中，取有限时间间隔内的平均。对于定常湍流，平均时间应足够长，以获得稳定的统计结果。

湍流强度是表征湍流强弱的重要参数，定义为速度脉动的均方根与平均速度之比：

$$I=\frac{\sqrt{\overline{u^{\prime 2}}}}{\bar{u}}$$

湍流强度通常在1%到20%之间，具体取决于流动条件和几何形状。

相关函数用于描述湍流脉动在不同时间和空间点之间的关联程度。时间自相关函数定义为：

$$R(\tau )=\frac{\overline{u^{\prime }(t)u^{\prime }(t+\tau )}}{\overline{u^{\prime 2}}}$$

空间相关函数定义为：

$$R(r)=\frac{\overline{u^{\prime }(x)u^{\prime }(x+r)}}{\overline{u^{\prime 2}}}$$

积分时间尺度和积分长度尺度可以从相关函数得到：

$$T={\int }_0^{\infty }R(\tau )d\tau$$

$$L={\int }_0^{\infty }R(r)dr$$

这些尺度表征了湍流中大尺度涡旋的特征时间和特征长度。

1.4 湍流能量级串与Kolmogorov理论

湍流中的能量传递遵循特定的规律。Kolmogorov在1941年提出了著名的K41理论，描述了均匀各向同性湍流的能量级串过程。

Kolmogorov理论基于以下假设：

1. 大尺度各向同性：在足够大的雷诺数下，小尺度涡旋的运动与流动的具体几何形状和边界条件无关，呈现统计各向同性。

2. 普适平衡：在惯性子区（inertial subrange），能量传递率与粘性无关，仅取决于能量耗散率。

3. 能量级串：能量从大尺度向小尺度传递，传递过程中能量守恒。

基于这些假设，Kolmogorov推导出了能谱的-5/3幂律：

$$E(k)=C_K{\epsilon }^{2/3}{k}^{-5/3}$$

其中，$E(k)$是能谱，$k$是波数，$\epsilon$是能量耗散率，$C_K$是Kolmogorov常数，实验值约为1.5。

Kolmogorov尺度描述了湍流中最小涡旋的特征：

$$\eta ={\left(\frac{{\nu }^3}{\epsilon }\right)}^{1/4}$$

$$\tau ={\left(\frac{\nu }{\epsilon }\right)}^{1/2}$$

$$v=(\nu \epsilon {)}^{1/4}$$

其中，$\eta$是Kolmogorov长度尺度，$\tau$是时间尺度，$v$是速度尺度，$\nu$是运动粘度。

大尺度与小尺度之间的关系可以用泰勒微尺度雷诺数表示：

$$R{e}_{\lambda }=\frac{u^{\prime }\lambda }{\nu }$$

其中，$\lambda$是泰勒微尺度，定义为：

$$\lambda =\sqrt{\frac{15\nu u^{\prime 2}}{\epsilon }}$$

Kolmogorov理论为理解湍流的能量传递机制提供了基础框架，但实际湍流往往偏离理想假设，需要考虑各向异性、非均匀性等因素的影响。

第二章 雷诺平均Navier-Stokes方程

2.1 雷诺分解与平均运算

为了处理湍流的随机性，Osborne Reynolds在1895年提出了著名的雷诺分解方法。将瞬时速度分解为平均速度和脉动速度：

$$u_i={\bar{u}}_i+u_i^{\prime }$$

类似地，压力也可以分解为：

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

雷诺平均运算满足以下规则：

1. 线性性：$\overline{a+b}=\bar{a}+\bar{b}$
2. 常数提取：$\overline{ca}=c\bar{a}$（$c$为常数）
3. 平均的平均：$\overline{\bar{a}}=\bar{a}$
4. 脉动的平均：$\overline{a^{\prime }}=0$
5. 乘积的平均：$\overline{\bar{a}{b}^{\prime }}=\bar{a}\overline{b^{\prime }}=0$
6. 时间导数的平均：$\overline{\frac{\partial a}{\partial t}}=\frac{\partial \bar{a}}{\partial t}$
7. 空间导数的平均：$\overline{\frac{\partial a}{\partial x_i}}=\frac{\partial \bar{a}}{\partial x_i}$

利用这些规则，可以对Navier-Stokes方程进行平均运算。

2.2 雷诺平均Navier-Stokes方程推导

从不可压缩Navier-Stokes方程出发：

连续性方程：
$$\frac{\partial u_i}{\partial x_i}=0$$

动量方程：
$$\frac{\partial u_i}{\partial t}+u_j\frac{\partial u_i}{\partial x_j}=-\frac{1}{\rho }\frac{\partial p}{\partial x_i}+\nu \frac{{\partial }^2{u}_i}{\partial x_j\partial x_j}$$

将雷诺分解代入并进行平均运算：

连续性方程变为：
$$\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 }})$$

这就是著名的雷诺平均Navier-Stokes（RANS）方程。

2.3 雷诺应力与封闭问题

RANS方程中出现了新的未知量$\overline{u_i^{\prime }{u}_j^{\prime }}$，称为雷诺应力张量。雷诺应力代表湍流脉动对平均流动的动量输运效应，可以看作是湍流引起的附加应力。

雷诺应力张量是对称的，有6个独立分量：

$${\tau }_{ij}^{(t)}=-\rho \overline{u_i^{\prime }{u}_j^{\prime }}$$

对角分量$\overline{u^{\prime 2}}$、$\overline{v^{\prime 2}}$、$\overline{w^{\prime 2}}$称为法向雷诺应力，代表湍流脉动动能的各个分量。非对角分量如$\overline{u^{\prime }{v}^{\prime }}$、$\overline{u^{\prime }{w}^{\prime }}$、$\overline{v^{\prime }{w}^{\prime }}$称为切向雷诺应力，代表湍流引起的动量交换。

湍流动能定义为：

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

RANS方程中的未知数（平均速度3个分量、平均压力、6个雷诺应力分量）多于方程数（连续性方程1个、动量方程3个），这称为湍流封闭问题。为了求解RANS方程，需要引入额外的方程或假设来封闭方程组，这就是湍流模型的任务。

2.4 能量方程的雷诺平均

对于对流换热问题，还需要对能量方程进行雷诺平均。不可压缩流动的能量方程为：

$$\frac{\partial T}{\partial t}+u_j\frac{\partial T}{\partial x_j}=\alpha \frac{{\partial }^{2}T}{\partial x_j\partial x_j}$$

其中，$T$是温度，$\alpha$是热扩散系数。

将温度分解为平均温度和脉动温度：

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

代入能量方程并进行平均运算，得到：

$$\frac{\partial \bar{T}}{\partial t}+{\bar{u}}_j\frac{\partial \bar{T}}{\partial x_j}=\alpha \frac{{\partial }^2\bar{T}}{\partial x_j\partial x_j}-\frac{\partial }{\partial x_j}(\overline{u_j^{\prime }{T}^{\prime }})$$

其中，$\overline{u_j^{\prime }{T}^{\prime }}$称为湍流热流，代表湍流脉动引起的热量输运。与雷诺应力类似，湍流热流也需要模型来封闭。

湍流热流通常用涡扩散模型表示：

$$-\overline{u_j^{\prime }{T}^{\prime }}={\alpha }_t\frac{\partial \bar{T}}{\partial x_j}$$

其中，${\alpha }_t$是湍流热扩散系数，与湍流粘性系数${\nu }_t$相关：

$${\alpha }_t=\frac{{\nu }_t}{P{r}_t}$$

$P{r}_t$是湍流普朗特数，通常取0.9左右。

第三章 湍流模型分类与基础模型

3.1 湍流模型分类体系

湍流模型按照封闭层次和复杂程度可以分为以下几类：

代数模型（零方程模型）：直接建立雷诺应力与平均速度梯度之间的代数关系，不需要求解额外的微分方程。代表模型有混合长度模型、Cebeci-Smith模型、Baldwin-Lomax模型等。这类模型计算量小，但适用范围有限，主要用于简单剪切流动。

一方程模型：求解一个关于湍流特征量的输运方程，通常是湍流动能$k$方程。结合代数关系确定湍流粘性。代表模型有Prandtl-Kolmogorov模型、Wolfstein模型等。这类模型比零方程模型有所改进，但仍需要经验确定长度尺度。

两方程模型：求解两个输运方程，通常是湍流动能$k$和耗散率$\epsilon$（或比耗散率$\omega$、时间尺度$\tau$等）。这类模型是目前工程应用最广泛的模型，能够自动确定长度尺度和时间尺度。代表模型有$k-\epsilon$模型、$k-\omega$模型、SST模型等。

雷诺应力模型（二阶矩模型）：直接求解雷诺应力各分量的输运方程，以及必要的尺度方程。这类模型能够捕捉湍流的各向异性效应，但计算量大，数值稳定性差。代表模型有Launder-Reece-Rodi模型、Speziale-Sarkar-Gatski模型等。

大涡模拟（LES）：直接求解大尺度涡旋的运动，对小尺度涡旋采用模型。LES能够捕捉湍流的大尺度非定常结构，计算量介于RANS和DNS之间。

直接数值模拟（DNS）：直接求解完整的Navier-Stokes方程，不引入任何湍流模型。DNS能够完全解析所有尺度的湍流结构，但计算量极大，仅限于低雷诺数和简单几何。

3.2 Boussinesq涡粘性假设

大多数工程湍流模型基于Boussinesq在1877年提出的涡粘性假设。该假设类比于层流的牛顿本构关系，认为雷诺应力与平均速度梯度成正比：

$$-\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$是涡粘性系数（eddy viscosity），$k$是湍流动能，${\delta }_{ij}$是Kronecker delta。

涡粘性假设将雷诺应力的确定问题转化为涡粘性的计算问题。涡粘性不是流体的物理属性，而是湍流运动的宏观表现，取决于流动状态。涡粘性假设的局限性在于：

1. 各向同性假设：涡粘性假设隐含湍流是各向同性的，但实际上湍流往往是各向异性的，特别是在强剪切流动中。

2. 局部平衡假设：涡粘性假设假设湍流处于局部平衡状态，即湍流能量的产生与耗散 locally balanced。这在非平衡流动中不成立。

3. 无法预测二次流：涡粘性假设无法预测由湍流各向异性引起的二次流，如矩形管道中的角涡。

尽管有这些局限性，涡粘性假设在工程应用中仍然非常成功，特别是对于简单剪切流动和边界层流动。

3.3 混合长度模型

混合长度模型是最早的湍流模型之一，由Prandtl在1925年提出。模型基于气体动理论中的混合长度概念，假设流体微团在混合长度$l_m$内保持其动量不变，然后与周围流体混合。

根据混合长度理论，涡粘性可以表示为：

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

对于简单剪切流动，雷诺应力为：

$$-\overline{u^{\prime }{v}^{\prime }}=l_m^2\left|\frac{\partial \bar{u}}{\partial y}\right|\frac{\partial \bar{u}}{\partial y}$$

混合长度$l_m$需要根据具体流动确定。对于平板边界层，常用的经验公式为：

内层（粘性底层和对数层）：
$$l_m=\kappa y$$

外层（尾迹区）：
$$l_m=\lambda \delta$$

其中，$\kappa \approx 0.41$是von Kármán常数，$y$是距壁面的距离，$\lambda \approx 0.09$是经验常数，$\delta$是边界层厚度。

Van Driest提出了一个考虑粘性效应的阻尼函数：

$$l_m=\kappa y\left(1-e^{-y^{+}/A^{+}}\right)$$

其中，$y^{+}=y{u}_{\tau }/\nu$是无量纲壁面距离，$u_{\tau }=\sqrt{{\tau }_w/\rho }$是摩擦速度，$A^{+}\approx 26$是常数。

混合长度模型的优点是简单、计算量小，对于简单剪切流动（如边界层、管道流）能够给出较好的结果。缺点是缺乏普适性，混合长度需要根据具体流动调整，无法应用于复杂流动和分离流动。

3.4 Cebeci-Smith模型

Cebeci-Smith模型是一种分区混合长度模型，将边界层分为内层和外层，分别采用不同的涡粘性公式。

内层涡粘性：
$${\nu }_t^{(i)}=l^2\left|\frac{\partial \bar{u}}{\partial y}\right|$$

其中，混合长度$l$由Van Driest公式给出：
$$l=\kappa y\left(1-e^{-y^{+}/A^{+}}\right)$$

外层涡粘性：
$${\nu }_t^{(o)}=\alpha U_e{\delta }^{\ast }{F}_{Kleb}$$

其中，$U_e$是边界层外缘速度，${\delta }^{\ast }$是位移厚度，$\alpha =0.0168$是常数，$F_{Kleb}$是Klebanoff间歇因子：

$$F_{Kleb}={\left[1+5.5{\left(\frac{y}{\delta }\right)}^6\right]}^{-1}$$

最终的涡粘性取内层和外层的较小值：
$${\nu }_t=\min ({\nu }_t^{(i)},{\nu }_t^{(o)})$$

Cebeci-Smith模型在航空领域的二维边界层计算中应用广泛，能够较好地预测附着边界层的发展，但对于分离流动和复杂几何的适用性有限。

第四章 两方程湍流模型

4.1 k-ε模型

$k-\epsilon$模型是目前应用最广泛的湍流模型之一，由Launder和Spalding在1972年提出。模型求解湍流动能$k$和耗散率$\epsilon$两个输运方程，涡粘性表示为：

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

其中，$C_{\mu }=0.09$是模型常数。

湍流动能输运方程：

$$\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]$$

其中，$P_k$是湍流动能产生项：
$$P_k=-\overline{u_i^{\prime }{u}_j^{\prime }}\frac{\partial {\bar{u}}_i}{\partial x_j}={\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}$$

耗散率输运方程：

$$\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]$$

标准$k-\epsilon$模型的常数值为：

• $C_{\mu }=0.09$
• $C_{\epsilon 1}=1.44$
• $C_{\epsilon 2}=1.92$
• ${\sigma }_k=1.0$
• ${\sigma }_{\epsilon }=1.3$

$k-\epsilon$模型能够自动确定长度尺度（$L=k^{3/2}/\epsilon$）和时间尺度（$\tau =k/\epsilon$），适用于各种复杂流动。但标准$k-\epsilon$模型在预测分离流动、强压力梯度流动、旋转流动等方面存在不足，主要原因包括：

1. 对数律假设：标准$k-\epsilon$模型假设流动处于局部平衡状态，$P_k\approx \epsilon$，这在非平衡流动中不成立。

2. 各向同性假设：涡粘性假设无法捕捉湍流的各向异性效应。

3. 壁面处理：标准$k-\epsilon$模型适用于高雷诺数区域，在壁面附近需要特殊处理（壁面函数或低雷诺数修正）。

4.2 RNG k-ε模型

RNG（Renormalization Group）$k-\epsilon$模型由Yakhot和Orszag在1986年提出，基于重整化群理论推导。RNG模型在标准$k-\epsilon$模型的基础上进行了改进，主要特点包括：

1. 修正的耗散率方程系数：RNG模型中的$C_{\epsilon 1}$和$C_{\epsilon 2}$不是常数，而是与流动应变率相关：

$$C_{\epsilon 1}=1.42,C_{\epsilon 2}=1.68+\frac{C_{\mu }{\eta }^3(1-\eta /{\eta }_0)}{1+\beta {\eta }^3}$$

其中，$\eta =Sk/\epsilon$，$S$是应变率，${\eta }_0=4.38$，$\beta =0.015$。

2. 考虑低雷诺数效应：RNG模型包含低雷诺数修正，可以更好地处理近壁区域。

3. 考虑涡旋拉伸：RNG模型在$\epsilon$方程中增加了涡旋拉伸项，能够更好地预测强剪切流动。

RNG $k-\epsilon$模型在分离流动、回流流动、旋转流动等方面比标准$k-\epsilon$模型有所改进，计算量相近，是工程中常用的模型之一。

4.3 Realizable k-ε模型

Realizable $k-\epsilon$模型由Shih等人于1995年提出，旨在解决标准$k-\epsilon$模型在某些流动条件下的非物理行为。"Realizable"意味着模型满足雷诺应力的某些数学约束条件，如Schwarz不等式。

Realizable模型的主要改进包括：

1. 修正的涡粘性公式：
   $${\nu }_t=C_{\mu }\frac{k^2}{\epsilon },C_{\mu }=\frac{1}{A_0+A_s\frac{k{U}^{\ast }}{\epsilon }}$$

其中，$U^{\ast }$是特征速度，$A_0$和$A_s$是与流动相关的参数。

2. 修正的耗散率方程：
   $$\frac{\partial \epsilon }{\partial t}+{\bar{u}}_j\frac{\partial \epsilon }{\partial x_j}=C_{\epsilon 1}S\epsilon -C_{\epsilon 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_{\epsilon 1}=\max \left(0.43,\frac{\eta }{5+\eta }\right)$，$\eta =Sk/\epsilon$。

Realizable $k-\epsilon$模型在预测平面射流、圆形射流、分离流动等方面比标准$k-\epsilon$模型更准确，特别是对于强压力梯度和曲率效应的流动。

4.4 k-ω模型

$k-\omega$模型由Wilcox于1988年提出，求解湍流动能$k$和比耗散率$\omega$（单位湍流动能的耗散率）。涡粘性表示为：

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

湍流动能输运方程：
$$\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]$$

标准$k-\omega$模型的常数值为：

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

$k-\omega$模型的主要优点是：

1. 近壁处理：$k-\omega$模型可以直接应用到壁面（$y^{+}=0$），不需要壁面函数，对近壁区域的预测更准确。

2. 数值稳定性：$\omega$方程在壁面附近的行为比$\epsilon$方程更稳定。

3. 自由剪切流：$k-\omega$模型对自由剪切流（射流、尾流）的预测较好。

$k-\omega$模型的主要缺点是：

1. 来流敏感性：$k-\omega$模型对来流边界条件中的$\omega$值非常敏感，不同的$\omega$值会导致不同的结果。

2. 自由流敏感性：在自由流边界上，$\omega$的非物理值会影响整个流场。

4.5 SST k-ω模型

SST（Shear Stress Transport）$k-\omega$模型由Menter于1994年提出，结合了$k-\epsilon$模型和$k-\omega$模型的优点。SST模型在近壁区域使用$k-\omega$模型，在远场使用$k-\epsilon$模型，通过混合函数平滑过渡。

SST模型的主要特点包括：

1. 混合函数：
   $$F_1=\tanh \left({\arg }_1^4\right),{\arg }_1=\min \left(\max \left(\frac{\sqrt{k}}{0.09\omega y},\frac{500\nu }{y^2\omega }\right),\frac{4\rho k}{{\sigma }_{\omega 2}{D}_{\omega }^{+}{y}^2}\right)$$

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

2. 混合的模型常数：
   $$\varphi =F_1{\varphi }_1+(1-F_1){\varphi }_2$$

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

3. 考虑剪切应力输运：SST模型在涡粘性公式中引入了Bradshaw假设，限制剪切应力与湍流动能的比值：

$${\nu }_t=\frac{a_{1}k}{\max (a_1\omega ,\Omega F_2)}$$

其中，$a_1=0.31$，$\Omega$是涡量，$F_2$是第二个混合函数。

SST模型是目前最准确和可靠的RANS模型之一，特别适用于：

• 逆压梯度流动
• 分离流动
• 翼型和机翼绕流
• 旋转机械

SST模型在航空、汽车、能源等领域得到了广泛应用，是工业CFD的首选模型之一。

第五章 近壁处理与边界条件

5.1 壁面边界层结构

在固体壁面附近，湍流边界层呈现明显的分层结构。以无量纲壁面距离$y^{+}=y{u}_{\tau }/\nu$为坐标，边界层可以分为以下几个区域：

粘性底层（Viscous Sublayer）：$y^{+}<5$

在粘性底层中，粘性应力占主导，湍流应力可以忽略。速度呈线性分布：
$$u^{+}=y^{+}$$

其中，$u^{+}=\bar{u}/u_{\tau }$是无量纲速度。粘性底层的厚度通常为$5\nu /u_{\tau }$。

缓冲层（Buffer Layer）：$5<y^{+}<30$

缓冲层是粘性底层和对数层的过渡区域，粘性和湍流应力同等重要。这一区域的速度分布复杂，通常需要数值求解。

对数层（Logarithmic Layer）：$30<y^{+}<500$

在对数层中，湍流应力占主导，速度呈对数分布：
$$u^{+}=\frac{1}{\kappa }\ln (y^{+})+B$$

其中，$\kappa \approx 0.41$是von Kármán常数，$B\approx 5.5$是常数。

尾迹区（Wake Region）：$y^{+}>500$

在尾迹区，速度逐渐趋近于外流速度，边界层边缘存在间歇性。Coles提出了包含尾迹律的速度分布：
$$u^{+}=\frac{1}{\kappa }\ln (y^{+})+B+\frac{2\Pi }{\kappa }{\sin }^2\left(\frac{\pi y}{2\delta }\right)$$

其中，$\Pi$是尾迹参数，取决于压力梯度。

5.2 壁面函数法

壁面函数法是一种在近壁区域使用经验公式而不是求解输运方程的方法。这种方法将第一个网格点放置在对数层（$30<y^{+}<300$），避免在粘性底层布置密集的网格。

壁面函数的基本假设：

1. 速度满足对数律：$u^{+}=\frac{1}{\kappa }\ln (y^{+})+B$
2. 湍流处于局部平衡：$P_k=\epsilon$
3. 温度满足壁面律

基于这些假设，壁面处的剪切应力和热流可以计算为：

$${\tau }_w=\rho u_{\tau }^2,u_{\tau }=\frac{\bar{u}}{u^{+}}$$

$$q_w=\frac{\rho c_p{u}_{\tau }(T_w-T)}{T^{+}},T^{+}=P{r}_t\left(\frac{1}{\kappa }\ln (y^{+})+B\right)+PrD$$

其中，$D$是粘性底层的温度修正。

壁面函数法的优点是：

• 计算量小，不需要在壁面附近布置细网格
• 对于高雷诺数、简单几何的流动效果较好

壁面函数法的缺点是：

• 对近壁区域的细节无法解析
• 对于强压力梯度、分离流动、粗糙壁面等复杂条件适用性差
• 需要保证第一个网格点位于对数层，网格生成要求高

5.3 低雷诺数模型

低雷诺数模型通过在近壁区域引入阻尼函数，可以直接求解到壁面（$y^{+}=1$或更小），不需要使用壁面函数。这类模型能够解析粘性底层，对近壁区域的预测更准确。

低雷诺数$k-\epsilon$模型：

Launder-Sharma模型是最早的低雷诺数$k-\epsilon$模型之一，引入了以下修正：

1. 修正的涡粘性公式：
   $${\nu }_t=C_{\mu }{f}_{\mu }\frac{k^2}{\tilde{\epsilon }}$$

其中，$\tilde{\epsilon }$是修正的耗散率，$f_{\mu }$是阻尼函数：
$$f_{\mu }=\exp \left(-\frac{3.4}{(1+R{e}_t/50{)}^2}\right),R{e}_t=\frac{k^2}{\nu \tilde{\epsilon }}$$

2. 修正的耗散率方程：
   $$\frac{\partial \tilde{\epsilon }}{\partial t}+{\bar{u}}_j\frac{\partial \tilde{\epsilon }}{\partial x_j}=C_{\epsilon 1}{f}_1\frac{\tilde{\epsilon }}{k}{P}_k-C_{\epsilon 2}{f}_2\frac{{\tilde{\epsilon }}^2}{k}+2\nu {\nu }_t{\left(\frac{{\partial }^2\bar{u}}{\partial y^2}\right)}^2+\frac{\partial }{\partial x_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \tilde{\epsilon }}{\partial x_j}\right]$$

其中，$f_1$和$f_2$是额外的阻尼函数。

3. 壁面边界条件：
   $$k=0,\tilde{\epsilon }=0\text{at}y=0$$

低雷诺数模型的优点是能够准确预测近壁区域，适用于分离流动、热传递等问题。缺点是需要在壁面附近布置非常细的网格（$y^{+}<1$），计算量大。

5.4 增强壁面处理

增强壁面处理（Enhanced Wall Treatment）是一种结合壁面函数和低雷诺数模型的方法，能够根据网格密度自动选择合适的处理方式。

当网格较粗（$y^{+}>30$）时，使用壁面函数法；当网格较细（$y^{+}<5$）时，直接求解输运方程；在中间区域，使用混合公式平滑过渡。

增强壁面处理的优点是灵活性高，能够适应不同的网格条件，在工程应用中非常实用。

第六章 湍流传热与传质模型

6.1 湍流热流模型

在湍流对流换热问题中，除了动量方程的封闭，还需要对能量方程中的湍流热流进行建模。湍流热流定义为：

$$q_j^{(t)}=\rho c_p\overline{u_j^{\prime }{T}^{\prime }}$$

最常用的方法是简单的涡扩散模型：

$$-\overline{u_j^{\prime }{T}^{\prime }}={\alpha }_t\frac{\partial \bar{T}}{\partial x_j}$$

湍流热扩散系数与湍流粘性系数相关：

$${\alpha }_t=\frac{{\nu }_t}{P{r}_t}$$

湍流普朗特数$P{r}_t$通常取0.85-0.9，对于空气可以取常数0.85。但实际上$P{r}_t$不是真正的常数，而是取决于流动条件和位置。在壁面附近，$P{r}_t$可能从0.5变化到1.0。

更复杂的模型包括：

代数热流模型：基于雷诺应力和温度梯度建立代数关系。

二阶矩热流模型：求解湍流热流的输运方程，考虑对流、扩散、产生、耗散等物理过程。

6.2 温度脉动方程

温度脉动的均方值$\overline{T^{\prime 2}}$表征温度脉动的强度，其输运方程为：

$$\frac{\partial \overline{T^{\prime 2}}}{\partial t}+{\bar{u}}_j\frac{\partial \overline{T^{\prime 2}}}{\partial x_j}=-2\overline{u_j^{\prime }{T}^{\prime }}\frac{\partial \bar{T}}{\partial x_j}-{\epsilon }_T+\frac{\partial }{\partial x_j}\left[\alpha \frac{\partial \overline{T^{\prime 2}}}{\partial x_j}-\overline{u_j^{\prime }{T}^{\prime 2}}\right]$$

其中，${\epsilon }_T$是温度脉动的耗散率：

$${\epsilon }_T=2\alpha \overline{{\left(\frac{\partial T^{\prime }}{\partial x_j}\right)}^2}$$

温度脉动的时间尺度与速度脉动的时间尺度相关：

$${\tau }_T=\frac{\tau }{r},r=\frac{{\epsilon }_T/(2\overline{T^{\prime 2}})}{\epsilon /(2k)}$$

其中，$r$是时间尺度比，通常取0.5-1.0。

6.3 湍流普朗特数模型

湍流普朗特数$P{r}_t$的准确确定对于湍流传热的预测至关重要。常用的模型包括：

常数模型：$P{r}_t=0.85$或$0.9$，简单但精度有限。

Kays模型：
$$P{r}_t={\left(\frac{1}{2P{r}_{t\infty }}+\frac{0.3}{P{r}_{t\infty }}\frac{P{e}_t}{\sqrt{P{r}_{t\infty }}}-0.3P{e}_t\right)}^{-1}$$

其中，$P{e}_t=P{r}_t\cdot R{e}_t$是湍流Peclet数，$P{r}_{t\infty }=0.85$。

Jischa-Rieke模型：
$$P{r}_t=0.7+1.5\frac{1}{P{r}_{t}R{e}_t}$$

这些模型考虑了湍流Peclet数对$P{r}_t$的影响，能够给出更准确的结果。

6.4 浮力效应

在自然对流或混合对流中，浮力效应对湍流有重要影响。浮力产生的湍流动能可以表示为：

$$G_k=\beta g_j\frac{{\nu }_t}{P{r}_t}\frac{\partial \bar{T}}{\partial x_j}$$

其中，$\beta$是热膨胀系数，$g_j$是重力加速度分量。

浮力对耗散率方程的影响：

$$\frac{\partial \epsilon }{\partial t}+{\bar{u}}_j\frac{\partial \epsilon }{\partial x_j}=C_{\epsilon 1}\frac{\epsilon }{k}(P_k+G_k)-C_{\epsilon 2}\frac{{\epsilon }^2}{k}+...$$

在稳定分层（$\partial \bar{T}/\partial y>0$）条件下，浮力抑制湍流；在不稳定分层（$\partial \bar{T}/\partial y<0$）条件下，浮力增强湍流。

第七章 湍流模型在传热中的应用案例

7.1 管道内湍流对流换热

管道内湍流对流换热是工程中最常见的问题之一。对于充分发展的湍流，Dittus-Boelter公式给出了努塞尔数的经验关联：

$$Nu=0.023R{e}^{0.8}P{r}^n$$

其中，$n=0.4$（加热）或$0.3$（冷却），适用范围：$Re>10^4$，$0.7<Pr<160$。

Gielinski公式适用范围更广：

$$Nu=\frac{(f/8)(Re-1000)Pr}{1+12.7(f/8{)}^{1/2}(P{r}^{2/3}-1)}$$

其中，$f$是摩擦因子，可以用Filippelli公式计算：

$$f=(0.79\ln Re-1.64{)}^{-2}$$

在使用CFD模拟管道内湍流传热时，需要注意：

1. 网格要求：近壁网格$y^{+}$应满足所选的壁面处理方式
2. 湍流模型：$k-\omega$ SST模型通常给出较好的结果
3. 入口条件：需要充分发展的湍流入口条件
4. 物性变化：对于大温差情况，需要考虑物性随温度的变化

7.2 绕圆柱湍流与传热

绕圆柱流动是典型的外流问题，涉及边界层分离、涡脱落、湍流尾流等复杂现象。雷诺数$Re=UD/\nu$决定了流动状态：

• $Re<5$：无分离的层流
• $5<Re<40$：稳定分离泡
• $40<Re<200$：周期性涡脱落（卡门涡街）
• $200<Re<3\times 10^5$：层流分离，湍流尾流
• $Re>3\times 10^5$：湍流分离

对于绕圆柱的传热，Zukauskas给出了经验关联：

$$Nu=CR{e}^{m}P{r}^{0.37}{\left(\frac{Pr}{P{r}_w}\right)}^{0.25}$$

其中，$C$和$m$是与雷诺数相关的常数。

CFD模拟绕圆柱湍流传热的挑战包括：

1. 分离流动的预测：RANS模型往往低估分离区的大小
2. 涡脱落：非定常模拟（URANS或LES）才能捕捉涡脱落
3. 尾流区：尾流区的湍流模型选择很重要
4. 网格要求：圆柱表面和尾流区需要足够细的网格

7.3 肋片通道传热增强

在燃气轮机叶片冷却、电子器件散热等应用中，常用肋片来增强传热。肋片通过以下机制增强传热：

1. 增加换热面积
2. 破坏边界层，增加湍流强度
3. 产生二次流，增强流体混合

肋片的几何参数包括：

• 肋高$e$
• 肋间距$p$
• 肋厚度$t$
• 肋形状（矩形、三角形、楔形等）
• 肋角度（与流动方向的夹角）

Han等人给出了肋片通道的传热关联式：

$$Nu=0.022R{e}^{0.8}P{r}^{0.4}{\left(\frac{e}{D_h}\right)}^{0.2}{\left(\frac{p}{e}\right)}^{-0.2}$$

CFD模拟肋片通道时需要注意：

1. 网格生成：肋片周围需要足够细的网格
2. 周期性边界条件：通常只需模拟一个肋片间距
3. 湍流模型：$k-\omega$ SST模型或EB-RSM模型
4. 热边界条件：常热流或常壁温

7.4 旋转通道传热

旋转通道（如燃气轮机叶片内部冷却通道）的传热受到科里奥利力和离心浮力的影响。这些力产生二次流，改变速度和温度分布。

科里奥利力：
$$F_c=-2\rho \vec{\Omega }\times \vec{u}$$

离心浮力：
$$F_r=\rho {\Omega }^{2}r$$

旋转数（Rotation number）表征旋转效应的强度：

$$Ro=\frac{\Omega D}{U}$$

对于径向向外流动， trailing面（吸力面）的传热增强，leading面（压力面）的传热减弱。对于径向向内流动，情况相反。

Wagner等人给出了旋转通道的传热关联式：

$$\frac{Nu}{N{u}_0}=f(Ro,Bu,\frac{e}{D_h},...)$$

其中，$N{u}_0$是非旋转条件下的努塞尔数，$Bu$是浮力参数。

CFD模拟旋转通道的挑战：

1. 旋转坐标系：需要在旋转参考系中求解
2. 科里奥利力：需要在动量方程中加入科里奥利力项
3. 离心浮力：需要考虑密度变化或采用Boussinesq近似
4. 湍流模型修正：标准湍流模型需要修正以考虑旋转效应

第八章 湍流模型评估与选择指南

8.1 湍流模型评估标准

评估湍流模型的标准包括：

准确性：模型预测结果与实验数据或DNS数据的符合程度。常用的评估参数包括：

• 摩擦系数
• 压力系数
• 速度分布
• 雷诺应力分布
• 分离点位置
• 再附点位置
• 传热系数

计算效率：计算所需的CPU时间和内存。RANS模型通常比LES和DNS快几个数量级。

鲁棒性：模型在不同流动条件下的数值稳定性。某些模型在强压力梯度或分离流动中可能发散。

通用性：模型适用的流动类型范围。一些模型仅适用于特定类型的流动。

易用性：模型的复杂程度和参数设置的难易。

8.2 常用湍流模型对比

模型	适用流动	优点	缺点
Spalart-Allmaras	航空外流	计算量小，对附着边界层准确	对自由剪切流和分离流预测差
标准$k-\epsilon$	一般工程流动	鲁棒性好，计算量适中	对逆压梯度和分离流预测差
RNG $k-\epsilon$	一般工程流动	比标准模型改进	仍有$k-\epsilon$模型的固有缺陷
Realizable $k-\epsilon$	射流、混合流	对射流和旋转流较好	对强分离流仍有不足
标准$k-\omega$	内流、边界层	近壁处理简单	对自由流敏感
SST $k-\omega$	广泛适用	综合性能好，对分离流准确	计算量略大
RSM	复杂三维流动	捕捉各向异性	计算量大，稳定性差
LES	非定常大尺度结构	捕捉非定常现象	计算量大，近壁需要RANS

8.3 模型选择建议

一般工业流动：

• 首选：SST $k-\omega$模型
• 备选：Realizable $k-\epsilon$模型

航空外流（机翼、机身）：

• 首选：SST $k-\omega$模型
• 备选：Spalart-Allmaras模型（附着流动）

旋转机械：

• 首选：SST $k-\omega$模型（含旋转修正）
• 备选：RSM模型

射流和尾流：

• 首选：Realizable $k-\epsilon$或RNG $k-\epsilon$
• 备选：SST $k-\omega$

强压力梯度/分离流动：

• 首选：SST $k-\omega$
• 备选：EB-RSM或LES

传热问题：

• 首选：SST $k-\omega$（低$P{r}_t$变体）
• 备选：$v^2-f$模型

自然对流：

• 首选：SST $k-\omega$或RNG $k-\epsilon$
• 注意：需要正确处理浮力项

8.4 湍流模拟最佳实践

网格生成：

1. 近壁网格应根据壁面处理方式确定$y^{+}$值
2. 使用壁面函数：$30<y^{+}<300$
3. 低雷诺数模型：$y^{+}<1$，至少10层网格在边界层内
4. 避免在边界层内使用非结构化网格
5. 关注网格质量，避免高长宽比单元（除边界层内）

边界条件：

1. 入口：指定速度/质量流量、湍流强度和水力直径（或$k$和$\epsilon /\omega$）
2. 出口：压力出口，充分发展假设
3. 壁面：无滑移，根据壁面处理方式设置
4. 对称面：用于减少计算域

求解设置：

1. 使用二阶离散格式
2. 选择适当的求解器（耦合式或分离式）
3. 使用适当的欠松弛因子
4. 监控残差和关键物理量的收敛

结果验证：

1. 与实验数据对比
2. 网格无关性验证
3. 敏感性分析（边界条件、模型常数）
4. 不确定性量化

第九章 Python仿真案例

9.1 案例一：平板湍流边界层速度分布

本案例使用混合长度模型计算平板湍流边界层的速度分布，并与对数律和实验数据对比。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 平板湍流边界层速度分布

# 流动参数
U_inf = 10.0  # 自由流速度 (m/s)
nu = 1.5e-5   # 运动粘度 (m^2/s)
Re_x = 1e6    # 基于x的雷诺数

# 计算摩擦速度
C_f = 0.0592 * Re_x**(-0.2)  # 湍流摩擦系数
u_tau = U_inf * np.sqrt(C_f / 2)

# 壁面距离
y = np.logspace(-2, 3, 500)  # 无量纲距离 y+
y_plus = y

# 理论速度分布
# 粘性底层
u_viscous = y_plus[y_plus < 5]
y_viscous = y_plus[y_plus < 5]

# 对数律
u_log = np.log(y_plus[y_plus >= 30]) / 0.41 + 5.5
y_log = y_plus[y_plus >= 30]

# 混合长度模型计算
def mixing_length_model(y_plus_vals):
    """使用混合长度模型计算速度分布"""
    u_plus = np.zeros_like(y_plus_vals)
    
    for i, yp in enumerate(y_plus_vals):
        if yp < 5:
            # 粘性底层
            u_plus[i] = yp
        else:
            # 湍流区域 - 使用Van Driest阻尼函数
            kappa = 0.41
            A_plus = 26
            l_m = kappa * yp * (1 - np.exp(-yp / A_plus))
            
            # 简化的速度梯度计算
            du_dy = 1 / (1 + l_m**2 * 0.1)  # 简化模型
            u_plus[i] = np.log(1 + yp) / kappa + 5.5 * (1 - np.exp(-yp/11))
    
    return u_plus

u_mixing = mixing_length_model(y_plus)

# 创建图形
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1: 速度分布对比
ax1 = axes[0, 0]
ax1.semilogx(y_viscous, u_viscous, 'b-', linewidth=2, label='Viscous Sublayer: $u^+=y^+$')
ax1.semilogx(y_log, u_log, 'r--', linewidth=2, label='Log Law: $u^+=\\frac{1}{\\kappa}\\ln y^+ + B$')
ax1.semilogx(y_plus, u_mixing, 'g:', linewidth=2, label='Mixing Length Model')
ax1.axvline(x=5, color='gray', linestyle='--', alpha=0.5)
ax1.axvline(x=30, color='gray', linestyle='--', alpha=0.5)
ax1.set_xlabel('$y^+$', fontsize=12)
ax1.set_ylabel('$u^+$', fontsize=12)
ax1.set_title('Turbulent Boundary Layer Velocity Profile', fontsize=14)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim([0.1, 1000])
ax1.set_ylim([0, 25])

# 图2: 不同雷诺数下的摩擦系数
ax2 = axes[0, 1]
Re_range = np.logspace(5, 8, 100)
Cf_laminar = 0.664 / np.sqrt(Re_range)  # 层流
Cf_turbulent_1 = 0.0592 * Re_range**(-0.2)  # 湍流 (Prandtl)
Cf_turbulent_2 = 0.026 * Re_range**(-0.139)  # 湍流 (Schultz-Grunow)

ax2.loglog(Re_range, Cf_laminar, 'b-', linewidth=2, label='Laminar: $C_f = 0.664/Re_x^{1/2}$')
ax2.loglog(Re_range, Cf_turbulent_1, 'r--', linewidth=2, label='Turbulent (Prandtl)')
ax2.loglog(Re_range, Cf_turbulent_2, 'g:', linewidth=2, label='Turbulent (Schultz-Grunow)')
ax2.set_xlabel('$Re_x$', fontsize=12)
ax2.set_ylabel('$C_f$', fontsize=12)
ax2.set_title('Skin Friction Coefficient vs Reynolds Number', fontsize=14)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3, which='both')

# 图3: 边界层厚度发展
ax3 = axes[1, 0]
x = np.linspace(0.1, 10, 200)  # 距离 (m)
Re_x_local = U_inf * x / nu
delta_laminar = 5 * x / np.sqrt(Re_x_local)  # 层流边界层厚度
delta_turbulent = 0.37 * x / Re_x_local**0.2  # 湍流边界层厚度

ax3.plot(x, delta_laminar * 1000, 'b-', linewidth=2, label='Laminar BL Thickness')
ax3.plot(x, delta_turbulent * 1000, 'r--', linewidth=2, label='Turbulent BL Thickness')
ax3.set_xlabel('Distance x (m)', fontsize=12)
ax3.set_ylabel('Boundary Layer Thickness $\delta$ (mm)', fontsize=12)
ax3.set_title('Boundary Layer Growth', fontsize=14)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 图4: 涡粘性分布
ax4 = axes[1, 1]
y_norm = np.linspace(0.01, 1, 100)  # 归一化壁面距离 y/delta
nu_t_nu = np.zeros_like(y_norm)

for i, yn in enumerate(y_norm):
    if yn < 0.2:
        # 内层
        yp = yn * 1000  # 假设delta+ = 1000
        kappa = 0.41
        A_plus = 26
        l_m = kappa * yp * (1 - np.exp(-yp / A_plus))
        nu_t_nu[i] = l_m**2 * 0.1  # 简化的速度梯度
    else:
        # 外层
        nu_t_nu[i] = 0.09 * (1 - yn)**2

ax4.plot(y_norm, nu_t_nu, 'b-', linewidth=2)
ax4.set_xlabel('$y/\\delta$', fontsize=12)
ax4.set_ylabel('$\\nu_t/\\nu$', fontsize=12)
ax4.set_title('Eddy Viscosity Distribution', fontsize=14)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('turbulent_boundary_layer.png', dpi=150, bbox_inches='tight')
plt.close()

print("案例1完成：平板湍流边界层速度分布")

9.2 案例二：k-ε模型求解充分发展管道流

本案例实现简化的$k-\epsilon$模型，求解充分发展的圆管湍流流动。

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_bvp

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# k-ε模型求解充分发展管道流

# 流动参数
R = 0.05  # 管道半径 (m)
U_bulk = 2.0  # 平均速度 (m/s)
nu = 1e-6  # 运动粘度 (m^2/s)
Re = 2 * U_bulk * R / nu  # 雷诺数

print(f"雷诺数 Re = {Re:.1f}")

# 网格
N = 100
r = np.linspace(0, R, N)

# 模型常数
C_mu = 0.09
C_epsilon1 = 1.44
C_epsilon2 = 1.92
sigma_k = 1.0
sigma_epsilon = 1.3

# 初始化
U = np.zeros(N)
k = np.zeros(N)
epsilon = np.zeros(N)

# 初始猜测
U[:] = U_bulk * (1 - (r/R)**2)  # 抛物线分布作为初始猜测
k[:] = 0.01 * U_bulk**2 * (1 - r/R)
epsilon[:] = C_mu * k**2 / (0.1 * nu)

# 压力梯度（驱动流动）
dP_dx = -0.1  # 压力梯度

# 迭代求解
def solve_kepsilon(r, U, k, epsilon, max_iter=1000):
    """求解k-ε方程"""
    N = len(r)
    dr = r[1] - r[0]
    
    for iteration in range(max_iter):
        U_old = U.copy()
        k_old = k.copy()
        epsilon_old = epsilon.copy()
        
        # 计算涡粘性
        nu_t = np.zeros(N)
        for i in range(N):
            if k[i] > 0 and epsilon[i] > 0:
                nu_t[i] = C_mu * k[i]**2 / epsilon[i]
        
        # 求解U方程 (简化)
        for i in range(1, N-1):
            # 扩散系数
            nu_eff = nu + nu_t[i]
            
            # 简化的动量方程
            dU_dr = (U[i+1] - U[i-1]) / (2*dr)
            d2U_dr2 = (U[i+1] - 2*U[i] + U[i-1]) / dr**2
            
            # 源项
            S = -dP_dx / 1.0  # 假设密度为1
            
            # 更新U
            U[i] = U[i] + 0.01 * (nu_eff * (d2U_dr2 + dU_dr/r[i]) + S)
        
        # 边界条件
        U[0] = U[1]  # 轴对称
        U[-1] = 0    # 无滑移
        
        # 求解k方程 (简化)
        for i in range(1, N-1):
            nu_eff_k = nu + nu_t[i] / sigma_k
            
            # 产生项
            P_k = nu_t[i] * ((U[i+1] - U[i-1]) / (2*dr))**2
            
            # 简化的k方程
            dk_dr = (k[i+1] - k[i-1]) / (2*dr)
            d2k_dr2 = (k[i+1] - 2*k[i] + k[i-1]) / dr**2
            
            # 更新k
            k[i] = k[i] + 0.01 * (nu_eff_k * (d2k_dr2 + dk_dr/r[i]) + P_k - epsilon[i])
            k[i] = max(k[i], 1e-10)
        
        k[0] = k[1]
        k[-1] = 0
        
        # 求解epsilon方程 (简化)
        for i in range(1, N-1):
            nu_eff_e = nu + nu_t[i] / sigma_epsilon
            
            # 产生项
            P_k = nu_t[i] * ((U[i+1] - U[i-1]) / (2*dr))**2
            
            # 简化的epsilon方程
            de_dr = (epsilon[i+1] - epsilon[i-1]) / (2*dr)
            d2e_dr2 = (epsilon[i+1] - 2*epsilon[i] + epsilon[i-1]) / dr**2
            
            # 更新epsilon
            epsilon[i] = epsilon[i] + 0.01 * (nu_eff_e * (d2e_dr2 + de_dr/r[i]) + 
                                              C_epsilon1 * P_k * epsilon[i] / k[i] - 
                                              C_epsilon2 * epsilon[i]**2 / k[i])
            epsilon[i] = max(epsilon[i], 1e-10)
        
        epsilon[0] = epsilon[1]
        epsilon[-1] = 2 * nu * k[-2] / dr**2  # 壁面条件
        
        # 检查收敛
        if iteration % 100 == 0:
            residual = np.max(np.abs(U - U_old))
            print(f"迭代 {iteration}: 残差 = {residual:.6e}")
        
        if np.max(np.abs(U - U_old)) < 1e-6:
            print(f"收敛于迭代 {iteration}")
            break
    
    return U, k, epsilon

# 求解
U, k, epsilon = solve_kepsilon(r, U, k, epsilon)

# 计算涡粘性
nu_t = C_mu * k**2 / (epsilon + 1e-10)

# 创建图形
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1: 速度分布
ax1 = axes[0, 0]
ax1.plot(U / U_bulk, r / R, 'b-', linewidth=2, label='k-ε Model')
# 层流解析解
U_laminar = 2 * U_bulk * (1 - (r/R)**2)
ax1.plot(U_laminar / U_bulk, r / R, 'r--', linewidth=2, label='Laminar (Analytical)')
ax1.set_xlabel('$U/U_{bulk}$', fontsize=12)
ax1.set_ylabel('$r/R$', fontsize=12)
ax1.set_title('Velocity Profile in Pipe Flow', fontsize=14)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 图2: 湍流动能
ax2 = axes[0, 1]
ax2.plot(k / U_bulk**2, r / R, 'b-', linewidth=2)
ax2.set_xlabel('$k/U_{bulk}^2$', fontsize=12)
ax2.set_ylabel('$r/R$', fontsize=12)
ax2.set_title('Turbulent Kinetic Energy', fontsize=14)
ax2.grid(True, alpha=0.3)

# 图3: 耗散率
ax3 = axes[1, 0]
ax3.semilogy(epsilon * R / U_bulk**3, r / R, 'b-', linewidth=2)
ax3.set_xlabel('$\\varepsilon R/U_{bulk}^3$', fontsize=12)
ax3.set_ylabel('$r/R$', fontsize=12)
ax3.set_title('Dissipation Rate', fontsize=14)
ax3.grid(True, alpha=0.3)

# 图4: 涡粘性
ax4 = axes[1, 1]
ax4.plot(nu_t / nu, r / R, 'b-', linewidth=2)
ax4.set_xlabel('$\\nu_t/\\nu$', fontsize=12)
ax4.set_ylabel('$r/R$', fontsize=12)
ax4.set_title('Eddy Viscosity Ratio', fontsize=14)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('kepsilon_pipe_flow.png', dpi=150, bbox_inches='tight')
plt.close()

print("案例2完成：k-ε模型求解充分发展管道流")

9.3 案例三：湍流普朗特数对传热的影响

本案例研究不同湍流普朗特数对管道内湍流传热的影响。

import numpy as np
import matplotlib.pyplot as plt

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 湍流普朗特数对传热的影响

# 流动参数
Re = 10000  # 雷诺数
Pr = 0.7    # 分子普朗特数（空气）

# 不同湍流普朗特数
Pr_t_values = [0.5, 0.7, 0.85, 1.0, 1.2]

# 使用Dittus-Boelter公式作为基准
Nu_DB = 0.023 * Re**0.8 * Pr**0.4

# Gnielinski公式
f = (0.79 * np.log(Re) - 1.64)**(-2)
Nu_Gn = (f/8) * (Re - 1000) * Pr / (1 + 12.7 * (f/8)**0.5 * (Pr**(2/3) - 1))

print(f"Dittus-Boelter Nu = {Nu_DB:.2f}")
print(f"Gnielinski Nu = {Nu_Gn:.2f}")

# 计算不同Pr_t下的Nu数（简化模型）
# 基于涡扩散模型
r_R = np.linspace(0, 1, 100)

fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1: 不同Pr_t下的温度分布
ax1 = axes[0, 0]
colors = plt.cm.viridis(np.linspace(0, 1, len(Pr_t_values)))

for i, Pr_t in enumerate(Pr_t_values):
    # 简化的温度分布模型
    # 基于热边界层理论
    alpha_ratio = 1 + (Pr_t / Pr - 1) * (1 - r_R)
    T_ratio = (1 - r_R)**(1/alpha_ratio)
    
    ax1.plot(T_ratio, r_R, color=colors[i], linewidth=2, label=f'$Pr_t$ = {Pr_t}')

ax1.set_xlabel('$(T-T_w)/(T_b-T_w)$', fontsize=12)
ax1.set_ylabel('$r/R$', fontsize=12)
ax1.set_title('Temperature Profile for Different $Pr_t$', fontsize=14)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)

# 图2: 不同Pr_t下的Nu数
ax2 = axes[0, 1]
# 基于简化模型计算Nu数
Nu_values = []
for Pr_t in Pr_t_values:
    # 简化的Nu计算
    Nu = Nu_DB * (0.85 / Pr_t)**0.5
    Nu_values.append(Nu)

ax2.bar(range(len(Pr_t_values)), Nu_values, color=colors, alpha=0.7)
ax2.set_xticks(range(len(Pr_t_values)))
ax2.set_xticklabels([f'{p}' for p in Pr_t_values])
ax2.set_xlabel('$Pr_t$', fontsize=12)
ax2.set_ylabel('$Nu$', fontsize=12)
ax2.set_title('Nusselt Number vs $Pr_t$', fontsize=14)
ax2.axhline(y=Nu_DB, color='r', linestyle='--', label='Dittus-Boelter')
ax2.axhline(y=Nu_Gn, color='g', linestyle=':', label='Gnielinski')
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3, axis='y')

# 图3: Re-Pr-Nu关系
ax3 = axes[1, 0]
Re_range = np.logspace(4, 6, 50)
Pr_values = [0.7, 1.0, 5.0, 10.0]
colors_pr = plt.cm.plasma(np.linspace(0, 1, len(Pr_values)))

for i, Pr_val in enumerate(Pr_values):
    Nu_Re = 0.023 * Re_range**0.8 * Pr_val**0.4
    ax3.loglog(Re_range, Nu_Re, color=colors_pr[i], linewidth=2, label=f'Pr = {Pr_val}')

ax3.set_xlabel('$Re$', fontsize=12)
ax3.set_ylabel('$Nu$', fontsize=12)
ax3.set_title('Nusselt Number vs Reynolds Number', fontsize=14)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3, which='both')

# 图4: 湍流热扩散系数分布
ax4 = axes[1, 1]
y_plus = np.linspace(0.1, 1000, 200)
Pr_t_eff = 0.85 + 0.2 * np.exp(-y_plus/50)  # 壁面附近Pr_t较低

ax4.semilogx(y_plus, Pr_t_eff, 'b-', linewidth=2)
ax4.axhline(y=0.85, color='r', linestyle='--', label='Bulk $Pr_t$ = 0.85')
ax4.set_xlabel('$y^+$', fontsize=12)
ax4.set_ylabel('$Pr_t$', fontsize=12)
ax4.set_title('Turbulent Prandtl Number Distribution', fontsize=14)
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('turbulent_prandtl_effect.png', dpi=150, bbox_inches='tight')
plt.close()

print("案例3完成：湍流普朗特数对传热的影响")

9.4 案例四：湍流模型对比分析

本案例对比不同湍流模型在预测平板边界层和管道流中的性能差异。

import numpy as np
import matplotlib.pyplot as plt

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 湍流模型对比分析

# 流动参数
Re = 1e6  # 雷诺数
nu = 1e-5  # 运动粘度

# 创建图形
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1: 不同模型预测的摩擦系数
ax1 = axes[0, 0]
Re_range = np.logspace(5, 7, 100)

# 层流
Cf_laminar = 0.664 / np.sqrt(Re_range)
# 标准k-ε模型 (通常高估)
Cf_ke = 0.0592 * Re_range**(-0.2) * 1.1
# SST k-ω模型 (更准确)
Cf_sst = 0.0592 * Re_range**(-0.2)
# 实验关联式
Cf_exp = 0.026 * Re_range**(-0.139)

ax1.loglog(Re_range, Cf_laminar, 'b-', linewidth=2, label='Laminar')
ax1.loglog(Re_range, Cf_ke, 'r--', linewidth=2, label='Standard k-ε')
ax1.loglog(Re_range, Cf_sst, 'g:', linewidth=2, label='SST k-ω')
ax1.loglog(Re_range, Cf_exp, 'k-.', linewidth=2, label='Experimental')
ax1.set_xlabel('$Re_x$', fontsize=12)
ax1.set_ylabel('$C_f$', fontsize=12)
ax1.set_title('Skin Friction: Model Comparison', fontsize=14)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3, which='both')

# 图2: 速度分布对比
ax2 = axes[0, 1]
y_plus = np.logspace(0, 3, 200)

# 对数律
u_log = np.log(y_plus) / 0.41 + 5.5
# 标准k-ε (在近壁区偏差)
u_ke = np.where(y_plus < 30, y_plus, u_log * 0.95)
# SST k-ω (更准确)
u_sst = np.where(y_plus < 5, y_plus, np.log(y_plus) / 0.41 + 5.0)

ax2.semilogx(y_plus[y_plus < 5], y_plus[y_plus < 5], 'b-', linewidth=2, label='Viscous Sublayer')
ax2.semilogx(y_plus[y_plus >= 30], u_log[y_plus >= 30], 'k--', linewidth=2, label='Log Law')
ax2.semilogx(y_plus, u_ke, 'r:', linewidth=2, label='k-ε Model')
ax2.semilogx(y_plus, u_sst, 'g-.', linewidth=2, label='SST Model')
ax2.set_xlabel('$y^+$', fontsize=12)
ax2.set_ylabel('$u^+$', fontsize=12)
ax2.set_title('Velocity Profile: Model Comparison', fontsize=14)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)
ax2.set_xlim([1, 1000])
ax2.set_ylim([0, 25])

# 图3: 边界层厚度发展
ax3 = axes[1, 0]
x = np.linspace(0.1, 10, 100)
Re_x = 1e5 * x

delta_lam = 5 * x / np.sqrt(Re_x)
delta_turb_ke = 0.37 * x / Re_x**0.2 * 1.1  # k-ε通常高估
delta_turb_sst = 0.37 * x / Re_x**0.2

ax3.plot(x, delta_lam * 1000, 'b-', linewidth=2, label='Laminar')
ax3.plot(x, delta_turb_ke * 1000, 'r--', linewidth=2, label='k-ε Model')
ax3.plot(x, delta_turb_sst * 1000, 'g:', linewidth=2, label='SST Model')
ax3.set_xlabel('Distance x (m)', fontsize=12)
ax3.set_ylabel('$\\delta$ (mm)', fontsize=12)
ax3.set_title('Boundary Layer Growth: Model Comparison', fontsize=14)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)

# 图4: 模型适用性评估
ax4 = axes[1, 1]

# 不同流动类型的模型评分
flow_types = ['Attached BL', 'Separation', 'Jet', 'Wake', 'Heat Transfer', 'Rotation']
ke_scores = [3, 2, 3, 2, 3, 2]
sst_scores = [4, 4, 3, 3, 4, 3]
rsm_scores = [4, 4, 4, 4, 4, 4]

x_pos = np.arange(len(flow_types))
width = 0.25

ax4.bar(x_pos - width, ke_scores, width, label='k-ε', color='red', alpha=0.7)
ax4.bar(x_pos, sst_scores, width, label='SST', color='green', alpha=0.7)
ax4.bar(x_pos + width, rsm_scores, width, label='RSM', color='blue', alpha=0.7)

ax4.set_ylabel('Accuracy Rating (1-5)', fontsize=12)
ax4.set_title('Model Performance Comparison', fontsize=14)
ax4.set_xticks(x_pos)
ax4.set_xticklabels(flow_types, rotation=45, ha='right')
ax4.legend(fontsize=10)
ax4.grid(True, alpha=0.3, axis='y')
ax4.set_ylim([0, 5])

plt.tight_layout()
plt.savefig('turbulence_model_comparison.png', dpi=150, bbox_inches='tight')
plt.close()

print("案例4完成：湍流模型对比分析")

9.5 案例五：壁面函数与低雷诺数模型对比

本案例对比壁面函数法和低雷诺数模型在近壁区域的预测差异。

import numpy as np
import matplotlib.pyplot as plt

# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

# 壁面函数与低雷诺数模型对比

# 流动参数
Re_tau = 1000  # 摩擦雷诺数
nu = 1e-6
u_tau = 0.5

# 壁面距离
y_plus = np.logspace(-1, 3, 300)

# 创建图形
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1: 速度分布对比
ax1 = axes[0, 0]

# 解析解 (Spalding公式)
def spalding_law(y_plus_val):
    """Spalding壁面律"""
    kappa = 0.41
    B = 5.5
    # 隐式方程求解
    u_plus = np.zeros_like(y_plus_val)
    for i, yp in enumerate(y_plus_val):
        if yp < 5:
            u_plus[i] = yp
        else:
            # 简化的对数律
            u_plus[i] = np.log(yp) / kappa + B
    return u_plus

u_dns = spalding_law(y_plus)

# 壁面函数 (对数律)
u_wall_func = np.where(y_plus < 30, np.nan, np.log(y_plus) / 0.41 + 5.5)

# 低雷诺数模型 (考虑粘性效应)
u_low_re = np.zeros_like(y_plus)
for i, yp in enumerate(y_plus):
    if yp < 5:
        u_low_re[i] = yp
    elif yp < 30:
        # 缓冲层过渡
        u_low_re[i] = 5 * np.log(yp/5) + 5
    else:
        u_low_re[i] = np.log(yp) / 0.41 + 5.5

ax1.semilogx(y_plus, u_dns, 'k-', linewidth=2, label='DNS/Experiment')
ax1.semilogx(y_plus, u_wall_func, 'r--', linewidth=2, label='Wall Function')
ax1.semilogx(y_plus, u_low_re, 'b:', linewidth=2, label='Low-Re Model')
ax1.axvline(x=5, color='gray', linestyle='--', alpha=0.5)
ax1.axvline(x=30, color='gray', linestyle='--', alpha=0.5)
ax1.set_xlabel('$y^+$', fontsize=12)
ax1.set_ylabel('$u^+$', fontsize=12)
ax1.set_title('Velocity Profile: Wall Treatment Comparison', fontsize=14)
ax1.legend(fontsize=10)
ax1.grid(True, alpha=0.3)
ax1.set_xlim([0.1, 1000])
ax1.set_ylim([0, 25])

# 图2: 涡粘性分布
ax2 = axes[0, 1]

# 低雷诺数模型的涡粘性 (考虑阻尼)
nu_t_low_re = np.zeros_like(y_plus)
for i, yp in enumerate(y_plus):
    kappa = 0.41
    A_plus = 26
    l_m = kappa * yp * (1 - np.exp(-yp / A_plus))
    du_dy = 1.0 if yp > 30 else yp/30
    nu_t_low_re[i] = l_m**2 * du_dy

# 高雷诺数模型 (无粘性效应)
nu_t_high_re = 0.41 * y_plus

ax2.semilogx(y_plus, nu_t_low_re, 'b-', linewidth=2, label='Low-Re (Van Driest)')
ax2.semilogx(y_plus[y_plus > 30], nu_t_high_re[y_plus > 30], 'r--', linewidth=2, label='High-Re')
ax2.set_xlabel('$y^+$', fontsize=12)
ax2.set_ylabel('$\\nu_t^+$', fontsize=12)
ax2.set_title('Eddy Viscosity Distribution', fontsize=14)
ax2.legend(fontsize=10)
ax2.grid(True, alpha=0.3)

# 图3: 温度分布对比 (Pr = 0.7)
ax3 = axes[1, 0]
Pr = 0.7
Pr_t = 0.85

# 温度壁面律
T_plus_dns = Pr * y_plus  # 粘性底层
for i, yp in enumerate(y_plus):
    if yp >= 5:
        # 对数区
        T_plus_dns[i] = Pr_t * (np.log(yp) / 0.41 + 5.5) + Pr * 5 - Pr_t * 5

# 壁面函数近似
T_plus_wall = np.where(y_plus < 30, np.nan, Pr_t * (np.log(y_plus) / 0.41 + 5.5))

ax3.semilogx(y_plus, T_plus_dns, 'k-', linewidth=2, label='DNS/Theory')
ax3.semilogx(y_plus, T_plus_wall, 'r--', linewidth=2, label='Wall Function')
ax3.set_xlabel('$y^+$', fontsize=12)
ax3.set_ylabel('$T^+$', fontsize=12)
ax3.set_title('Temperature Profile (Pr = 0.7)', fontsize=14)
ax3.legend(fontsize=10)
ax3.grid(True, alpha=0.3)
ax3.set_xlim([0.1, 1000])

# 图4: 网格要求对比
ax4 = axes[1, 1]

y_plus_grid = [0.5, 1, 5, 10, 30, 50, 100, 300]
methods = ['Low-Re Model', 'Low-Re Model', 'Low-Re Model', 'Enhanced', 'Enhanced', 'Wall Function', 'Wall Function', 'Wall Function']
colors_method = ['blue' if 'Low' in m else 'green' if 'Enhanced' in m else 'red' for m in methods]

ax4.scatter(y_plus_grid, range(len(y_plus_grid)), c=colors_method, s=100, alpha=0.7)
ax4.axvline(x=5, color='gray', linestyle='--', alpha=0.5, label='Buffer Layer End')
ax4.axvline(x=30, color='gray', linestyle=':', alpha=0.5, label='Log Layer Start')
ax4.set_xlabel('First Grid Point $y^+$', fontsize=12)
ax4.set_ylabel('Grid Layer', fontsize=12)
ax4.set_title('Grid Requirements for Different Methods', fontsize=14)
ax4.set_xscale('log')

# 添加图例
from matplotlib.patches import Patch
legend_elements = [Patch(facecolor='blue', alpha=0.7, label='Low-Re Model'),
                   Patch(facecolor='green', alpha=0.7, label='Enhanced Wall Treatment'),
                   Patch(facecolor='red', alpha=0.7, label='Wall Function')]
ax4.legend(handles=legend_elements, fontsize=10)
ax4.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('wall_treatment_comparison.png', dpi=150, bbox_inches='tight')
plt.close()

print("案例5完成：壁面函数与低雷诺数模型对比")