对流换热仿真-主题061_大涡模拟LES-大涡模拟LES

kkchenjj

409人浏览 · 2026-04-20 20:57:21

kkchenjj · 2026-04-20 20:57:21 发布

主题061：大涡模拟（LES）

一、引言

1.1 湍流模拟的挑战

湍流是自然界和工程中最常见的流动状态，其特征是复杂的非定常、三维、多尺度的涡旋结构。湍流包含从最大能量含能尺度到最小耗散尺度的连续谱系，这种多尺度特性使得湍流的直接数值模拟（DNS）在计算资源上极为昂贵。

对于雷诺数 $Re = 10^6$ 的湍流流动，DNS所需的网格点数约为 $Re9/4≈1013Re^{9/4} \approx 10^{13}$ ，这对于当前的计算能力而言是不可行的。因此，工程实践中需要采用湍流模型来降低计算成本。
在这里插入图片描述

1.2 雷诺平均与大涡模拟

目前主流的湍流模拟方法包括：

雷诺平均Navier-Stokes（RANS）方法：

对所有湍流尺度进行统计平均
计算成本低，适合工程应用
无法捕捉湍流的瞬态特性和大尺度结构

大涡模拟（LES）方法：

直接解析大尺度涡旋，仅对小尺度涡旋进行建模
能够捕捉湍流的瞬态特性和相干结构
计算成本介于RANS和DNS之间

直接数值模拟（DNS）：

解析所有湍流尺度
计算成本极高，仅用于基础研究

1.3 LES的优势与应用领域

大涡模拟的主要优势：

高保真度：能够捕捉湍流的瞬态特性和大尺度相干结构
物理准确性：大尺度运动包含大部分能量，对流动特性起主导作用
适用范围广：特别适用于非定常流动、分离流动和复杂几何
与实验对比：可生成与实验测量（如PIV）直接对比的瞬态数据

LES的主要应用领域：

航空航天：飞行器气动特性、发动机燃烧室
能源工程：风力机尾流、核反应堆热工水力
环境工程：大气边界层、污染物扩散
建筑工程：建筑风环境、室内通风
汽车工业：发动机缸内流动、车辆气动

二、大涡模拟的基本原理

2.1 滤波概念

大涡模拟的核心思想是通过滤波将流动变量分解为大尺度（可解）分量和小尺度（亚格子）分量：

$f(x,t)=fˉ(x,t)+f′(x,t)f(\mathbf{x},t) = \bar{f}(\mathbf{x},t) + f'(\mathbf{x},t)$

其中：

$fˉ\bar{f}$ 是滤波后的变量（大尺度分量）
$f^{'}$ 是亚格子尺度变量（小尺度分量）

2.1.1 盒式滤波器

盒式滤波器（Box/Top-hat filter）是最常用的滤波器之一：

$G(x−x′)={1Δ3∣x−x′∣≤Δ20其他G(\mathbf{x}-\mathbf{x}') = \begin{cases} \frac{1}{\Delta^3} & |\mathbf{x}-\mathbf{x}'| \leq \frac{\Delta}{2} \\ 0 & \text{其他} \end{cases}$

滤波运算定义为：

$fˉ(x,t)=∫−∞∞G(x−x′)f(x′,t)dx′\bar{f}(\mathbf{x},t) = \int_{-\infty}^{\infty} G(\mathbf{x}-\mathbf{x}') f(\mathbf{x}',t) d\mathbf{x}'$

2.1.2 高斯滤波器

高斯滤波器具有更好的数学性质：

$G(x−x′)=(6πΔ2)3/2exp⁡(−6∣x−x′∣2Δ2)G(\mathbf{x}-\mathbf{x}') = \left(\frac{6}{\pi \Delta^2}\right)^{3/2} \exp\left(-\frac{6|\mathbf{x}-\mathbf{x}'|^2}{\Delta^2}\right)$

高斯滤波器的优点是光滑性好，不会产生高频振荡。

2.1.3 谱截断滤波器

在谱方法中常用的滤波器：

$fˉ^(k)={f^(k)∣k∣≤kc0∣k∣>kc\hat{\bar{f}}(k) = \begin{cases} \hat{f}(k) & |k| \leq k_c \\ 0 & |k| > k_c \end{cases}$

其中 $kc=π/Δk_c = \pi/\Delta$ 是截断波数。

2.2 滤波后的控制方程

对不可压缩Navier-Stokes方程进行滤波，得到LES控制方程：

连续性方程：

$∂uˉi∂xi=0\frac{\partial \bar{u}_i}{\partial x_i} = 0$

动量方程：

$∂uˉi∂t+∂(uˉiuˉj)∂xj=−1ρ∂pˉ∂xi+ν∂2uˉi∂xj∂xj−∂τij∂xj\frac{\partial \bar{u}_i}{\partial t} + \frac{\partial (\bar{u}_i \bar{u}_j)}{\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 \tau_{ij}}{\partial x_j}$

能量方程：

$∂Tˉ∂t+∂(uˉjTˉ)∂xj=α∂2Tˉ∂xj∂xj−∂qj∂xj\frac{\partial \bar{T}}{\partial t} + \frac{\partial (\bar{u}_j \bar{T})}{\partial x_j} = \alpha \frac{\partial^2 \bar{T}}{\partial x_j \partial x_j} - \frac{\partial q_j}{\partial x_j}$

2.3 亚格子应力与热流

滤波过程引入了新的未知项：

亚格子应力张量：

$τij=uiuj‾−uˉiuˉj\tau_{ij} = \overline{u_i u_j} - \bar{u}_i \bar{u}_j$

亚格子热流：

$qj=ujT‾−uˉjTˉq_j = \overline{u_j T} - \bar{u}_j \bar{T}$

这些项代表小尺度涡旋对大尺度运动的影响，需要通过亚格子模型来封闭。

2.4 亚格子应力的分解

亚格子应力可以进一步分解为：

$τij=Lij+Cij+Rij\tau_{ij} = L_{ij} + C_{ij} + R_{ij}$

其中：

Leonard应力： $Lij=uˉiuˉj‾−uˉiuˉjL_{ij} = \overline{\bar{u}_i \bar{u}_j} - \bar{u}_i \bar{u}_j$ （可计算部分）
交叉应力： $Cij=uˉiuj′‾+ui′uˉj‾C_{ij} = \overline{\bar{u}_i u'_j} + \overline{u'_i \bar{u}_j}$
Reynolds应力： $Rij=ui′uj′‾R_{ij} = \overline{u'_i u'_j}$

在大多数LES实现中，这三部分不单独处理，而是统一建模。

三、亚格子尺度模型

3.1 涡粘模型

3.1.1 Smagorinsky模型

Smagorinsky模型是最早的亚格子模型：

$τij−13τkkδij=−2νsgsSˉij\tau_{ij} - \frac{1}{3}\tau_{kk}\delta_{ij} = -2\nu_{sgs}\bar{S}_{ij}$

其中亚格子涡粘系数：

$νsgs=(CsΔ)2∣Sˉ∣\nu_{sgs} = (C_s \Delta)^2 |\bar{S}|$

应变率张量：

$Sˉij=12(∂uˉi∂xj+∂uˉj∂xi)\bar{S}_{ij} = \frac{1}{2}\left(\frac{\partial \bar{u}_i}{\partial x_j} + \frac{\partial \bar{u}_j}{\partial x_i}\right)$

应变率张量的模：

$∣Sˉ∣=2SˉijSˉij|\bar{S}| = \sqrt{2\bar{S}_{ij}\bar{S}_{ij}}$

Smagorinsky常数：

对于各向同性湍流： $Cs≈0.1−0.2C_s \approx 0.1-0.2$
对于壁面流动： $Cs≈0.065C_s \approx 0.065$

Smagorinsky模型的局限性：

需要经验确定 $C_s$
在层流区域产生过多的耗散
无法自动适应流动状态的变化

3.1.2 动态Smagorinsky模型

Germano等人提出的动态模型通过双滤波技术动态计算 $C_s$ ：

$Cs2=⟨LijMij⟩⟨MijMij⟩C_s^2 = \frac{\langle L_{ij}M_{ij} \rangle}{\langle M_{ij}M_{ij} \rangle}$

其中：

测试滤波与网格滤波之间的应力： $Lij=uˉiuˉj^−uˉ^iuˉ^jL_{ij} = \widehat{\bar{u}_i \bar{u}_j} - \hat{\bar{u}}_i \hat{\bar{u}}_j$
$Mij=2Δ2∣Sˉ∣Sˉij^−2Δ^2∣Sˉ^∣Sˉ^ijM_{ij} = 2\Delta^2 \widehat{|\bar{S}|\bar{S}_{ij}} - 2\hat{\Delta}^2 |\hat{\bar{S}}|\hat{\bar{S}}_{ij}$

动态模型的优点：

自动适应局部流动状态
在层流区域 $Cs→0C_s \rightarrow 0$
能够模拟能量反向传递（backscatter）

动态模型的缺点：

计算成本增加
$C_s^2$ 可能出现负值，导致数值不稳定
需要平均或限制处理

3.1.3 WALE模型

WALE（Wall-Adapting Local Eddy-viscosity）模型由Nicoud和Ducros提出：

$νsgs=(CwΔ)2(SijdSijd)3/2(SˉijSˉij)5/2+(SijdSijd)5/4\nu_{sgs} = (C_w \Delta)^2 \frac{(S_{ij}^d S_{ij}^d)^{3/2}}{(\bar{S}_{ij}\bar{S}_{ij})^{5/2} + (S_{ij}^d S_{ij}^d)^{5/4}}$

其中 $S_{ij}^d$ 是应变率张量的无迹部分：

$Sijd=12(gˉij2+gˉji2)−13δijgˉkk2S_{ij}^d = \frac{1}{2}\left(\bar{g}_{ij}^2 + \bar{g}_{ji}^2\right) - \frac{1}{3}\delta_{ij}\bar{g}_{kk}^2$

这里 $gˉij=∂uˉi/∂xj\bar{g}_{ij} = \partial \bar{u}_i / \partial x_j$ 。

WALE模型的优点：

在壁面附近自动趋于零
不需要动态计算
对近壁流动建模较好

3.2 结构函数模型

3.2.1 尺度相似模型

Bardina等人提出的尺度相似模型：

$τij=Css(uˉiuˉj^−uˉ^iuˉ^j)\tau_{ij} = C_{ss}(\widehat{\bar{u}_i \bar{u}_j} - \hat{\bar{u}}_i \hat{\bar{u}}_j)$

该模型基于假设：小尺度运动的统计特性与大尺度运动相似。

3.2.2 混合模型

结合涡粘模型和尺度相似模型：

$τij=uˉiuˉj^−uˉ^iuˉ^j⏟尺度相似−2CmixΔ2∣Sˉ^∣Sˉ^ij⏟涡粘修正\tau_{ij} = \underbrace{\widehat{\bar{u}_i \bar{u}_j} - \hat{\bar{u}}_i \hat{\bar{u}}_j}_{\text{尺度相似}} - \underbrace{2C_{mix}\Delta^2|\hat{\bar{S}}|\hat{\bar{S}}_{ij}}_{\text{涡粘修正}}$

3.3 一方程模型

3.3.1 亚格子动能模型

求解亚格子动能输运方程：

$∂ksgs∂t+uˉj∂ksgs∂xj=−τijSˉij−Cϵksgs3/2Δ+∂∂xj(νsgsσk∂ksgs∂xj)\frac{\partial k_{sgs}}{\partial t} + \bar{u}_j\frac{\partial k_{sgs}}{\partial x_j} = -\tau_{ij}\bar{S}_{ij} - C_\epsilon \frac{k_{sgs}^{3/2}}{\Delta} + \frac{\partial}{\partial x_j}\left(\frac{\nu_{sgs}}{\sigma_k}\frac{\partial k_{sgs}}{\partial x_j}\right)$

涡粘系数：

$νsgs=CkΔksgs\nu_{sgs} = C_k \Delta \sqrt{k_{sgs}}$

3.4 亚格子热流模型

3.4.1 梯度扩散模型

最简单的亚格子热流模型：

$qj=−νsgsPrsgs∂Tˉ∂xjq_j = -\frac{\nu_{sgs}}{Pr_{sgs}}\frac{\partial \bar{T}}{\partial x_j}$

其中 $Pr_{sgs}$ 是亚格子Prandtl数，通常取 $0.4 - 0.6$ 。

3.4.2 动态热流模型

类似于动态Smagorinsky模型，动态计算 $Pr_{sgs}$ ：

$1Prsgs=⟨FjHj⟩⟨HjHj⟩\frac{1}{Pr_{sgs}} = \frac{\langle F_j H_j \rangle}{\langle H_j H_j \rangle}$

其中 $Fj=uˉjTˉ^−uˉ^jTˉ^F_j = \widehat{\bar{u}_j \bar{T}} - \hat{\bar{u}}_j \hat{\bar{T}}$ 。

四、近壁面处理

4.1 壁面附近湍流特性

壁面附近存在粘性底层，其特征：

无量纲距离 $y+=yuτ/ν<5y^+ = y u_\tau / \nu < 5$ 为粘性底层
$5 < y^+ < 30$ 为缓冲层
$y^+ > 30$ 为对数律层

在粘性底层内，速度呈线性分布： $u^+ = y^+$ 。

4.2 壁面解析LES

要求：

网格在壁面法向足够细， $Δy+<1\Delta y^+ < 1$
壁面切向分辨率： $Δx+<50\Delta x^+ < 50$ ， $Δz+<20\Delta z^+ < 20$

优点：直接解析近壁湍流结构
缺点：计算成本极高

4.3 壁面模型LES（WMLES）

4.3.1 对数律壁面模型

假设壁面剪切应力：

$τw=ρuτ2=ρ(κuˉln⁡(y/y0))2\tau_w = \rho u_\tau^2 = \rho \left(\frac{\kappa \bar{u}}{\ln(y/y_0)}\right)^2$

其中 $uτu_\tau$ 是摩擦速度， $κ=0.41\kappa = 0.41$ 是Karman常数。

4.3.2 代数壁面模型

通过壁面函数计算壁面剪切应力：

$u+=1κln⁡(Ey+)u^+ = \frac{1}{\kappa}\ln(E y^+)$

其中 $E = 9.8$ 是粗糙度参数。

4.3.3 双层模型

结合RANS和LES：

近壁区域使用RANS模型
远离壁面使用LES

五、LES的数值方法

5.1 空间离散格式

5.1.1 高阶有限差分

LES需要高阶空间离散来准确解析湍流结构：

4阶中心差分： $O(Δx4)O(\Delta x^4)$
6阶紧致格式： $O(Δx6)O(\Delta x^6)$
谱方法：指数收敛

5.1.2 谱方法

对于周期性边界条件，谱方法具有最高的精度：

$∂u∂x=∑kiku^(k)eikx\frac{\partial u}{\partial x} = \sum_{k} ik \hat{u}(k) e^{ikx}$

5.2 时间推进方法

5.2.1 显式方法

Runge-Kutta方法：

3阶Runge-Kutta：稳定性好，适合LES
4阶Runge-Kutta：精度高，但计算量大

Adams-Bashforth方法：

$un+1=un+Δt(32fn−12fn−1)u^{n+1} = u^n + \Delta t \left(\frac{3}{2}f^n - \frac{1}{2}f^{n-1}\right)$

5.2.2 隐式方法

对于壁面解析LES，近壁区域需要较小的时间步长，隐式方法可以提高稳定性。

5.3 压力-速度耦合

5.3.1 投影法

预测步（忽略压力梯度）：

$u∗−unΔt=−N(un)+L(un)\frac{u^* - u^n}{\Delta t} = -N(u^n) + L(u^n)$

压力修正：

$∇2pn+1=ρΔt∇⋅u∗\nabla^2 p^{n+1} = \frac{\rho}{\Delta t}\nabla \cdot u^*$

修正速度：

$un+1=u∗−Δtρ∇pn+1u^{n+1} = u^* - \frac{\Delta t}{\rho}\nabla p^{n+1}$

5.3.2 SIMPLE算法

对于低马赫数流动，可以使用SIMPLE算法处理压力-速度耦合。

5.4 边界条件

5.4.1 入口边界

合成湍流入口：

生成符合目标湍流统计特性的入口速度场
方法：傅里叶合成、数字滤波、回收/映射方法

回收/映射方法：

从下游某位置回收湍流数据
调整平均速度和湍流强度

5.4.2 出口边界

对流边界条件：

$∂u∂t+Uc∂u∂x=0\frac{\partial u}{\partial t} + U_c \frac{\partial u}{\partial x} = 0$

其中 $U_c$ 是对流速度。

特征边界条件：

基于特征分析，允许信息传出

5.4.3 周期性边界

对于充分发展的湍流，可以使用周期性边界条件：

$u (x + L) = u (x)$

5.5 计算网格要求

5.5.1 网格分辨率准则

壁面解析LES：

$Δx+<50\Delta x^+ < 50$ （流向）
$Δy+<1\Delta y^+ < 1$ （法向，第一层网格）
$Δz+<20\Delta z^+ < 20$ （展向）

壁面模型LES：

$Δx+≈100−200\Delta x^+ \approx 100-200$
$Δy+≈30−100\Delta y^+ \approx 30-100$
$Δz+≈50−100\Delta z^+ \approx 50-100$

5.5.2 网格质量

避免网格扭曲和过度拉伸
保持网格正交性
在关键区域（如剪切层）加密网格

六、LES的验证与确认

6.1 统计收敛性

湍流是随机过程，需要足够长的时间平均：

$⟨u⟩=1T∫0Tu(t)dt\langle u \rangle = \frac{1}{T}\int_0^T u(t) dt$

收敛准则：

一阶统计量（平均速度）收敛较快
二阶统计量（雷诺应力）需要更长时间
通常需要 $10 - 100$ 个流动通过时间

6.2 网格无关性验证

通过网格细化验证LES结果：

粗网格： $Δ=2Δfine\Delta = 2\Delta_{fine}$
中等网格： $Δ=Δfine\Delta = \Delta_{fine}$
细网格： $Δ=0.5Δfine\Delta = 0.5\Delta_{fine}$

6.3 与实验对比

6.3.1 平均量对比

平均速度分布
雷诺应力分量
湍流强度

6.3.2 瞬态特性对比

功率谱密度
两点相关函数
条件平均

6.4 与DNS/RANS对比

与DNS对比：验证亚格子模型的准确性
与RANS对比：展示LES的优势

# -*- coding: utf-8 -*-
"""
主题061：大涡模拟（LES）
Large Eddy Simulation (LES) Demonstration
"""

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.gridspec import GridSpec
from matplotlib.animation import FuncAnimation, PillowWriter
from matplotlib.patches import Circle, Rectangle
import matplotlib
matplotlib.use('Agg')
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False

import os
output_dir = '仿真结果'
os.makedirs(output_dir, exist_ok=True)

print("="*70)
print("主题061：大涡模拟（LES）")
print("="*70)

# ==================== 案例1：滤波函数与湍流能谱 ====================
print("\n【案例1：滤波函数与湍流能谱】")

def turbulence_spectrum_analysis():
    """分析湍流能谱和滤波效果"""
    
    # 波数范围
    k = np.linspace(0.1, 100, 1000)
    
    # Kolmogorov能谱 E(k) ~ k^(-5/3)
    E_k = k**(-5/3)
    E_k[0] = 0  # 避免无穷大
    
    # 截断波数
    k_c = 20
    
    # 不同滤波器的传递函数
    # 谱截断滤波器
    G_spectral = np.where(k <= k_c, 1.0, 0.0)
    
    # 高斯滤波器
    G_gaussian = np.exp(-(k/k_c)**2 / 2)
    
    # Box滤波器（近似）
    G_box = np.sinc(k/k_c)
    G_box = np.abs(G_box)
    
    # 滤波后的能谱
    E_filtered_spectral = E_k * G_spectral
    E_filtered_gaussian = E_k * G_gaussian**2
    E_filtered_box = E_k * G_box**2
    
    # 创建图形
    fig = plt.figure(figsize=(16, 10))
    gs = GridSpec(2, 3, figure=fig, hspace=0.3, wspace=0.3)
    
    # 子图1：原始湍流能谱
    ax1 = fig.add_subplot(gs[0, 0])
    ax1.loglog(k, E_k, 'b-', linewidth=2.5, label='Energy Spectrum E(k)')
    ax1.axvline(x=k_c, color='red', linestyle='--', linewidth=2, alpha=0.7, label=f'Cutoff k_c={k_c}')
    ax1.fill_between(k, E_k, alpha=0.3, color='blue')
    
    # 标注不同区域
    ax1.axvspan(0.1, k_c, alpha=0.1, color='green', label='Resolved (LES)')
    ax1.axvspan(k_c, 100, alpha=0.1, color='orange', label='Subgrid (SGS)')
    
    ax1.set_xlabel('Wavenumber k', fontsize=11)
    ax1.set_ylabel('E(k)', fontsize=11)
    ax1.set_title('Turbulence Energy Spectrum', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=9)
    ax1.grid(True, alpha=0.3, which='both')
    ax1.set_ylim([1e-4, 10])
    
    # 子图2：滤波器传递函数
    ax2 = fig.add_subplot(gs[0, 1])
    ax2.semilogx(k, G_spectral, 'b-', linewidth=2.5, label='Spectral Cutoff')
    ax2.semilogx(k, G_gaussian, 'r--', linewidth=2.5, label='Gaussian')
    ax2.semilogx(k, G_box, 'g:', linewidth=2.5, label='Box (Top-hat)')
    ax2.axvline(x=k_c, color='black', linestyle='-.', linewidth=1.5, alpha=0.5)
    ax2.set_xlabel('Wavenumber k', fontsize=11)
    ax2.set_ylabel('Filter Transfer Function G(k)', fontsize=11)
    ax2.set_title('Filter Transfer Functions', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=9)
    ax2.grid(True, alpha=0.3)
    ax2.set_ylim([0, 1.1])
    
    # 子图3：滤波后的能谱对比
    ax3 = fig.add_subplot(gs[0, 2])
    ax3.loglog(k, E_k, 'k-', linewidth=2, alpha=0.5, label='Original')
    ax3.loglog(k, E_filtered_spectral, 'b-', linewidth=2.5, label='Spectral Filtered')
    ax3.loglog(k, E_filtered_gaussian, 'r--', linewidth=2.5, label='Gaussian Filtered')
    ax3.loglog(k, E_filtered_box, 'g:', linewidth=2.5, label='Box Filtered')
    ax3.axvline(x=k_c, color='red', linestyle='--', linewidth=2, alpha=0.7)
    ax3.set_xlabel('Wavenumber k', fontsize=11)
    ax3.set_ylabel('Filtered E(k)', fontsize=11)
    ax3.set_title('Filtered Energy Spectra', fontsize=12, fontweight='bold')
    ax3.legend(fontsize=9)
    ax3.grid(True, alpha=0.3, which='both')
    ax3.set_ylim([1e-5, 10])
    
    # 子图4：能量累积分布
    ax4 = fig.add_subplot(gs[1, 0])
    E_cumulative = np.cumsum(E_k) / np.sum(E_k)
    E_cum_filtered = np.cumsum(E_filtered_gaussian) / np.sum(E_filtered_gaussian)
    
    ax4.plot(k, E_cumulative, 'b-', linewidth=2.5, label='Original')
    ax4.plot(k, E_cum_filtered, 'r--', linewidth=2.5, label='Filtered (Gaussian)')
    ax4.axvline(x=k_c, color='green', linestyle='--', linewidth=2, alpha=0.7)
    ax4.axhline(y=0.8, color='orange', linestyle=':', linewidth=2, alpha=0.7, label='80% Energy')
    ax4.set_xlabel('Wavenumber k', fontsize=11)
    ax4.set_ylabel('Cumulative Energy', fontsize=11)
    ax4.set_title('Cumulative Energy Distribution', fontsize=12, fontweight='bold')
    ax4.legend(fontsize=9)
    ax4.grid(True, alpha=0.3)
    ax4.set_xlim([0, 100])
    ax4.set_ylim([0, 1.1])
    
    # 子图5：湍流尺度示意图
    ax5 = fig.add_subplot(gs[1, 1])
    ax5.set_xlim(0, 10)
    ax5.set_ylim(0, 10)
    ax5.axis('off')
    ax5.set_title('Turbulent Scales Hierarchy', fontsize=12, fontweight='bold')
    
    # 绘制不同尺度的涡旋
    circle_large = Circle((2, 7), 1.5, fill=True, facecolor='red', alpha=0.3, edgecolor='red', linewidth=2)
    ax5.add_patch(circle_large)
    ax5.text(2, 7, 'L', ha='center', va='center', fontsize=20, fontweight='bold', color='red')
    ax5.text(2, 5, 'Energy\nContaining', ha='center', fontsize=9, color='red')
    
    circle_medium = Circle((5, 7), 1.0, fill=True, facecolor='blue', alpha=0.3, edgecolor='blue', linewidth=2)
    ax5.add_patch(circle_medium)
    ax5.text(5, 7, 'M', ha='center', va='center', fontsize=16, fontweight='bold', color='blue')
    ax5.text(5, 5, 'Inertial\nSubrange', ha='center', fontsize=9, color='blue')
    
    circle_small = Circle((8, 7), 0.6, fill=True, facecolor='green', alpha=0.3, edgecolor='green', linewidth=2)
    ax5.add_patch(circle_small)
    ax5.text(8, 7, 'S', ha='center', va='center', fontsize=12, fontweight='bold', color='green')
    ax5.text(8, 5, 'Dissipation\nRange', ha='center', fontsize=9, color='green')
    
    ax5.annotate('', xy=(3.5, 7), xytext=(2.5, 7), arrowprops=dict(arrowstyle='->', lw=2, color='purple'))
    ax5.annotate('', xy=(6.5, 7), xytext=(5.5, 7), arrowprops=dict(arrowstyle='->', lw=2, color='purple'))
    ax5.text(5, 2.5, 'Energy Cascade →', ha='center', fontsize=11, fontweight='bold', color='purple')
    
    # 子图6：LES概念说明
    ax6 = fig.add_subplot(gs[1, 2])
    ax6.axis('off')
    concept_text = """
    LES Concept:
    
    1. FILTERING:
       Separate large scales (resolved)
       from small scales (modeled)
    
    2. RESOLVED SCALES:
       - Carry most energy
       - Flow-dependent
       - Solved directly on grid
    
    3. SUBGRID SCALES:
       - Universal behavior
       - Isotropic (usually)
       - Modeled via SGS model
    
    4. CUTOFF:
       Grid size Delta determines
       the separation
       
    k_c = pi/Delta
    """
    ax6.text(0.5, 0.5, concept_text, transform=ax6.transAxes, fontsize=10,
            verticalalignment='center', horizontalalignment='center',
            bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8))
    
    plt.savefig(f'{output_dir}/case1_spectrum_filtering.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"  截断波数 k_c = {k_c}")
    print(f"  解析能量比例: {np.sum(E_filtered_gaussian)/np.sum(E_k)*100:.1f}%")
    print(f"  亚格子能量比例: {(1-np.sum(E_filtered_gaussian)/np.sum(E_k))*100:.1f}%")
    
    return {'k_c': k_c, 'energy_ratio': np.sum(E_filtered_gaussian)/np.sum(E_k)}

result_case1 = turbulence_spectrum_analysis()
print("✓ 案例1完成")

# ==================== 案例2：亚格子模型对比 ====================
print("\n【案例2：亚格子模型对比】")

def subgrid_models_comparison():
    """对比不同亚格子模型的性能"""
    
    # 创建简单的速度场（模拟剪切层）
    nx, ny = 64, 64
    x = np.linspace(0, 2*np.pi, nx)
    y = np.linspace(0, 2*np.pi, ny)
    X, Y = np.meshgrid(x, y)
    
    # 基础速度场（平均剪切）
    U_mean = 1.0
    u_base = U_mean * np.tanh(2*(Y - np.pi))
    v_base = np.zeros_like(X)
    
    # 添加湍流脉动
    np.random.seed(42)
    n_modes = 10
    u_fluct = np.zeros_like(X)
    v_fluct = np.zeros_like(X)
    
    for i in range(n_modes):
        kx = np.random.uniform(1, 5)
        ky = np.random.uniform(1, 5)
        amp = 1.0 / (kx**2 + ky**2)**(5/6)
        phase = np.random.uniform(0, 2*np.pi)
        u_fluct += amp * np.sin(kx*X + ky*Y + phase)
        v_fluct += amp * np.cos(kx*X + ky*Y + phase)
    
    # 归一化脉动
    u_fluct = 0.3 * u_fluct / (np.std(u_fluct) + 1e-10)
    v_fluct = 0.3 * v_fluct / (np.std(v_fluct) + 1e-10)
    
    u = u_base + u_fluct
    v = v_base + v_fluct
    
    # 计算应变率张量
    dx = x[1] - x[0]
    dy = y[1] - y[0]
    
    dudx = np.gradient(u, dx, axis=1)
    dudy = np.gradient(u, dy, axis=0)
    dvdx = np.gradient(v, dx, axis=1)
    dvdy = np.gradient(v, dy, axis=0)
    
    S11 = dudx
    S12 = 0.5 * (dudy + dvdx)
    S22 = dvdy
    
    S_mag = np.sqrt(2*(S11**2 + 2*S12**2 + S22**2))
    
    # 滤波宽度（网格尺寸）
    Delta = np.sqrt(dx**2 + dy**2)
    
    # 1. Smagorinsky模型
    Cs = 0.1
    nu_sgs_smag = (Cs * Delta)**2 * S_mag
    
    # 2. WALE模型
    Cw = 0.5
    g11 = dudx
    g12 = dudy
    g21 = dvdx
    g22 = dvdy
    
    Sd11 = 0.5*(g11**2 + g12*g21) - (1/3)*(g11**2 + g12*g21 + g21*g12 + g22**2)
    Sd12 = 0.5*(g11*g12 + g12*g22)
    Sd21 = 0.5*(g21*g11 + g22*g21)
    Sd22 = 0.5*(g21*g12 + g22**2) - (1/3)*(g11**2 + g12*g21 + g21*g12 + g22**2)
    
    Sd_mag = np.sqrt(Sd11**2 + Sd12**2 + Sd21**2 + Sd22**2)
    
    nu_sgs_wale = (Cw * Delta)**2 * Sd_mag**2.5 / (S_mag**2.5 + Sd_mag**2.5 + 1e-10)
    
    # 3. 动态模型（简化实现）
    from scipy.ndimage import uniform_filter
    u_test = uniform_filter(u, size=3)
    v_test = uniform_filter(v, size=3)
    
    dudx_test = np.gradient(u_test, dx, axis=1)
    dudy_test = np.gradient(u_test, dy, axis=0)
    dvdx_test = np.gradient(v_test, dx, axis=1)
    dvdy_test = np.gradient(v_test, dy, axis=0)
    
    S11_test = dudx_test
    S12_test = 0.5 * (dudy_test + dvdx_test)
    S22_test = dvdy_test
    S_mag_test = np.sqrt(2*(S11_test**2 + 2*S12_test**2 + S22_test**2))
    
    Cs_dyn = 0.1 * (S_mag_test / (S_mag + 1e-10))
    Cs_dyn = np.clip(Cs_dyn, 0, 0.2)
    nu_sgs_dyn = (Cs_dyn * Delta)**2 * S_mag
    
    # 创建图形
    fig = plt.figure(figsize=(16, 10))
    gs = GridSpec(2, 3, figure=fig, hspace=0.3, wspace=0.3)
    
    # 子图1：速度场
    ax1 = fig.add_subplot(gs[0, 0])
    speed = np.sqrt(u**2 + v**2)
    strm = ax1.streamplot(X, Y, u, v, color=speed, cmap='viridis', density=2, linewidth=1.5)
    plt.colorbar(strm.lines, ax=ax1, label='Velocity Magnitude')
    ax1.set_xlabel('x', fontsize=11)
    ax1.set_ylabel('y', fontsize=11)
    ax1.set_title('Velocity Field', fontsize=12, fontweight='bold')
    ax1.set_aspect('equal')
    
    # 子图2：Smagorinsky涡粘系数
    ax2 = fig.add_subplot(gs[0, 1])
    im2 = ax2.contourf(X, Y, nu_sgs_smag, levels=20, cmap='hot')
    plt.colorbar(im2, ax=ax2, label='nu_sgs')
    ax2.set_xlabel('x', fontsize=11)
    ax2.set_ylabel('y', fontsize=11)
    ax2.set_title(f'Smagorinsky Model (Cs={Cs})', fontsize=12, fontweight='bold')
    ax2.set_aspect('equal')
    
    # 子图3：WALE涡粘系数
    ax3 = fig.add_subplot(gs[0, 2])
    im3 = ax3.contourf(X, Y, nu_sgs_wale, levels=20, cmap='hot')
    plt.colorbar(im3, ax=ax3, label='nu_sgs')
    ax3.set_xlabel('x', fontsize=11)
    ax3.set_ylabel('y', fontsize=11)
    ax3.set_title(f'WALE Model (Cw={Cw})', fontsize=12, fontweight='bold')
    ax3.set_aspect('equal')
    
    # 子图4：动态Smagorinsky
    ax4 = fig.add_subplot(gs[1, 0])
    im4 = ax4.contourf(X, Y, nu_sgs_dyn, levels=20, cmap='hot')
    plt.colorbar(im4, ax=ax4, label='nu_sgs')
    ax4.set_xlabel('x', fontsize=11)
    ax4.set_ylabel('y', fontsize=11)
    ax4.set_title('Dynamic Smagorinsky Model', fontsize=12, fontweight='bold')
    ax4.set_aspect('equal')
    
    # 子图5：涡粘系数对比（沿中心线）
    ax5 = fig.add_subplot(gs[1, 1])
    mid = ny // 2
    ax5.plot(x, nu_sgs_smag[mid, :], 'b-', linewidth=2.5, label='Smagorinsky')
    ax5.plot(x, nu_sgs_wale[mid, :], 'r--', linewidth=2.5, label='WALE')
    ax5.plot(x, nu_sgs_dyn[mid, :], 'g:', linewidth=2.5, label='Dynamic')
    ax5.set_xlabel('x', fontsize=11)
    ax5.set_ylabel('nu_sgs', fontsize=11)
    ax5.set_title('Eddy Viscosity Comparison (Centerline)', fontsize=12, fontweight='bold')
    ax5.legend(fontsize=9)
    ax5.grid(True, alpha=0.3)
    
    # 子图6：模型参数表格
    ax6 = fig.add_subplot(gs[1, 2])
    ax6.axis('off')
    
    stats_data = [
        ['Model', 'Mean nu_sgs', 'Max nu_sgs', 'Model Constant'],
        ['Smagorinsky', f'{np.mean(nu_sgs_smag):.4f}', f'{np.max(nu_sgs_smag):.4f}', f'Cs={Cs}'],
        ['WALE', f'{np.mean(nu_sgs_wale):.4f}', f'{np.max(nu_sgs_wale):.4f}', f'Cw={Cw}'],
        ['Dynamic', f'{np.mean(nu_sgs_dyn):.4f}', f'{np.max(nu_sgs_dyn):.4f}', 'Variable Cs'],
    ]
    
    table = ax6.table(cellText=stats_data[1:], colLabels=stats_data[0], 
                     cellLoc='center', loc='center', bbox=[0, 0.3, 1, 0.4])
    table.auto_set_font_size(False)
    table.set_fontsize(10)
    table.scale(1, 2)
    for i in range(len(stats_data[0])):
        table[(0, i)].set_facecolor('#2196F3')
        table[(0, i)].set_text_props(weight='bold', color='white')
    ax6.set_title('SGS Model Comparison', fontsize=12, fontweight='bold', y=0.85)
    
    model_desc = """
    Model Features:
    
    Smagorinsky:
    - Simple and robust
    - Constant Cs
    - Too dissipative in laminar regions
    
    WALE:
    - Automatically vanishes at walls
    - No dynamic computation needed
    - Better near-wall behavior
    
    Dynamic:
    - Adapts to local flow
    - Can model backscatter
    - More expensive
    """
    ax6.text(0.5, 0.15, model_desc, transform=ax6.transAxes, fontsize=9,
            verticalalignment='top', horizontalalignment='center',
            bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.5))
    
    plt.savefig(f'{output_dir}/case2_sgs_models.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"  Smagorinsky平均涡粘: {np.mean(nu_sgs_smag):.4f}")
    print(f"  WALE平均涡粘: {np.mean(nu_sgs_wale):.4f}")
    print(f"  Dynamic平均涡粘: {np.mean(nu_sgs_dyn):.4f}")
    
    return {'nu_smag': nu_sgs_smag, 'nu_wale': nu_sgs_wale, 'nu_dyn': nu_sgs_dyn}

result_case2 = subgrid_models_comparison()
print("✓ 案例2完成")

# ==================== 案例3：2D衰减湍流模拟 ====================
print("\n【案例3：2D衰减湍流模拟】")

def decaying_turbulence_2d():
    """模拟2D衰减湍流 - 使用简化的谱方法"""
    
    # 网格参数
    nx, ny = 64, 64
    Lx, Ly = 2*np.pi, 2*np.pi
    dx = Lx / nx
    dy = Ly / ny
    x = np.linspace(0, Lx, nx, endpoint=False)
    y = np.linspace(0, Ly, ny, endpoint=False)
    X, Y = np.meshgrid(x, y)
    
    # 时间参数
    dt = 0.01
    nt = 300
    t = np.arange(nt) * dt
    
    # 波数
    kx = 2*np.pi * np.fft.fftfreq(nx, dx)
    ky = 2*np.pi * np.fft.fftfreq(ny, dy)
    KX, KY = np.meshgrid(kx, ky)
    k2 = KX**2 + KY**2
    k2[0, 0] = 1.0  # 避免除零
    k_mag = np.sqrt(k2)
    
    # 初始条件 - 使用衰减湍流能谱
    np.random.seed(42)
    
    # 生成初始速度场（物理空间）
    u = np.zeros((ny, nx))
    v = np.zeros((ny, nx))
    
    # 使用多个涡旋叠加
    n_vortices = 15
    for _ in range(n_vortices):
        x0 = np.random.uniform(0, Lx)
        y0 = np.random.uniform(0, Ly)
        strength = np.random.uniform(-0.3, 0.3)
        sigma = np.random.uniform(0.4, 0.8)
        
        r2 = (X - x0)**2 + (Y - y0)**2
        u += -strength * (Y - y0) / (sigma**2) * np.exp(-r2 / (2*sigma**2))
        v += strength * (X - x0) / (sigma**2) * np.exp(-r2 / (2*sigma**2))
    
    # 归一化
    u_max = np.max(np.sqrt(u**2 + v**2))
    if u_max > 0:
        u = u / u_max * 0.3  # 限制最大速度
        v = v / u_max * 0.3
    
    # 转换到谱空间
    u_hat = np.fft.fft2(u)
    v_hat = np.fft.fft2(v)
    
    # 物理参数
    nu = 0.03  # 粘性系数
    
    # 存储历史
    kinetic_energy = np.zeros(nt)
    enstrophy = np.zeros(nt)
    
    # 简化的衰减模拟（纯粘性衰减）
    for n in range(nt):
        # 物理空间速度（先计算当前状态）
        u = np.real(np.fft.ifft2(u_hat))
        v = np.real(np.fft.ifft2(v_hat))
        
        # 计算涡量
        omega_hat = 1j * KX * v_hat - 1j * KY * u_hat
        omega = np.real(np.fft.ifft2(omega_hat))
        
        # 计算动能和涡量
        kinetic_energy[n] = 0.5 * np.mean(u**2 + v**2)
        enstrophy[n] = 0.5 * np.mean(omega**2)
        
        # 纯粘性衰减
        damping = np.exp(-nu * k2 * dt)
        u_hat = u_hat * damping
        v_hat = v_hat * damping
    
    # 最终场
    u_final = np.real(np.fft.ifft2(u_hat))
    v_final = np.real(np.fft.ifft2(v_hat))
    omega_final = np.real(np.fft.ifft2(1j * KX * v_hat - 1j * KY * u_hat))
    
    # 创建图形
    fig = plt.figure(figsize=(16, 10))
    gs = GridSpec(2, 3, figure=fig, hspace=0.3, wspace=0.3)
    
    # 子图1：最终速度场
    ax1 = fig.add_subplot(gs[0, 0])
    speed = np.sqrt(u_final**2 + v_final**2)
    strm = ax1.streamplot(X, Y, u_final, v_final, color=speed, cmap='viridis', density=2, linewidth=1.5)
    plt.colorbar(strm.lines, ax=ax1, label='|u|')
    ax1.set_xlabel('x', fontsize=11)
    ax1.set_ylabel('y', fontsize=11)
    ax1.set_title(f'Final Velocity Field (t={t[-1]:.2f})', fontsize=12, fontweight='bold')
    ax1.set_aspect('equal')
    
    # 子图2：涡量场
    ax2 = fig.add_subplot(gs[0, 1])
    im2 = ax2.contourf(X, Y, omega_final, levels=20, cmap='RdBu_r')
    plt.colorbar(im2, ax=ax2, label='Vorticity omega')
    ax2.set_xlabel('x', fontsize=11)
    ax2.set_ylabel('y', fontsize=11)
    ax2.set_title('Final Vorticity Field', fontsize=12, fontweight='bold')
    ax2.set_aspect('equal')
    
    # 子图3：动能衰减
    ax3 = fig.add_subplot(gs[0, 2])
    ax3.semilogy(t, kinetic_energy, 'b-', linewidth=2.5, label='Kinetic Energy')
    # 理论衰减率
    t_theory = t[t > 0.1]
    k_theory = kinetic_energy[10] * np.exp(-2*nu*4*(t_theory - t[10]))
    ax3.semilogy(t_theory, k_theory, 'r--', linewidth=2, label='Theory: exp(-2nu*k^2*t)')
    ax3.set_xlabel('Time', fontsize=11)
    ax3.set_ylabel('Kinetic Energy', fontsize=11)
    ax3.set_title('Kinetic Energy Decay', fontsize=12, fontweight='bold')
    ax3.legend(fontsize=9)
    ax3.grid(True, alpha=0.3)
    
    # 子图4：涡量衰减
    ax4 = fig.add_subplot(gs[1, 0])
    ax4.semilogy(t, enstrophy, 'g-', linewidth=2.5, label='Enstrophy')
    ax4.set_xlabel('Time', fontsize=11)
    ax4.set_ylabel('Enstrophy', fontsize=11)
    ax4.set_title('Enstrophy Decay', fontsize=12, fontweight='bold')
    ax4.legend(fontsize=9)
    ax4.grid(True, alpha=0.3)
    
    # 子图5：能谱
    ax5 = fig.add_subplot(gs[1, 1])
    E_k = 0.5 * (np.abs(u_hat)**2 + np.abs(v_hat)**2)
    
    # 径向平均
    k_bins = np.linspace(0, np.max(k_mag), 30)
    E_k_avg = np.zeros(len(k_bins)-1)
    for i in range(len(k_bins)-1):
        mask = (k_mag >= k_bins[i]) & (k_mag < k_bins[i+1])
        if np.sum(mask) > 0:
            E_k_avg[i] = np.mean(E_k[mask])
    
    k_centers = 0.5 * (k_bins[:-1] + k_bins[1:])
    ax5.loglog(k_centers[1:], E_k_avg[1:], 'b-o', linewidth=2, markersize=4)
    # 添加-5/3斜率参考线
    k_ref = k_centers[5:15]
    if len(k_ref) > 0 and E_k_avg[5] > 0:
        E_ref = E_k_avg[5] * (k_ref / k_centers[5])**(-5/3)
        ax5.loglog(k_ref, E_ref, 'r--', linewidth=2, label='k^(-5/3)')
    ax5.set_xlabel('Wavenumber k', fontsize=11)
    ax5.set_ylabel('E(k)', fontsize=11)
    ax5.set_title('Energy Spectrum', fontsize=12, fontweight='bold')
    ax5.legend(fontsize=9)
    ax5.grid(True, alpha=0.3, which='both')
    
    # 子图6：模拟参数
    ax6 = fig.add_subplot(gs[1, 2])
    ax6.axis('off')
    param_text = f"""
    Simulation Parameters:
    
    Grid: {nx} x {ny}
    Domain: {Lx:.2f} x {Ly:.2f}
    Time step: {dt}
    Final time: {t[-1]:.3f}
    
    Physical Parameters:
    Molecular viscosity nu: {nu}
    
    Results:
    Initial KE: {kinetic_energy[0]:.4f}
    Final KE: {kinetic_energy[-1]:.4f}
    Decay ratio: {kinetic_energy[-1]/kinetic_energy[0]:.4f}
    
    Model: 2D Viscous Decay
    Method: Spectral + Projection
    """
    ax6.text(0.5, 0.5, param_text, transform=ax6.transAxes, fontsize=10,
            verticalalignment='center', horizontalalignment='center',
            bbox=dict(boxstyle='round', facecolor='lightblue', alpha=0.8))
    
    plt.savefig(f'{output_dir}/case3_decaying_turbulence.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"  网格: {nx} x {ny}")
    print(f"  初始动能: {kinetic_energy[0]:.4f}")
    print(f"  最终动能: {kinetic_energy[-1]:.4f}")
    print(f"  衰减比例: {kinetic_energy[-1]/kinetic_energy[0]:.4f}")
    
    return {'kinetic_energy': kinetic_energy, 'enstrophy': enstrophy, 't': t, 
            'u_hat': u_hat, 'v_hat': v_hat, 'kx': kx, 'ky': ky}

result_case3 = decaying_turbulence_2d()
print("✓ 案例3完成")

# ==================== 案例4：混合层不稳定性 ====================
print("\n【案例4：混合层不稳定性】")

def mixing_layer_instability():
    """模拟混合层Kelvin-Helmholtz不稳定性 - 使用伪谱方法"""
    
    # 网格参数
    nx, ny = 256, 128
    Lx, Ly = 4*np.pi, 2*np.pi
    dx = Lx / nx
    dy = Ly / ny
    x = np.linspace(0, Lx, nx, endpoint=False)
    y = np.linspace(-Ly/2, Ly/2, ny, endpoint=False)
    X, Y = np.meshgrid(x, y)
    
    # 时间参数
    dt = 0.01
    nt = 100
    t = np.arange(nt) * dt
    
    # 波数
    kx = 2*np.pi * np.fft.fftfreq(nx, dx)
    ky = 2*np.pi * np.fft.fftfreq(ny, dy)
    KX, KY = np.meshgrid(kx, ky)
    k2 = KX**2 + KY**2
    k2[0, 0] = 1e-10
    
    # 初始条件 - 混合层
    U1, U2 = 1.0, -1.0
    delta = 0.5
    
    # 平均速度（双曲正切剖面）
    U_mean = 0.5*(U1 + U2) + 0.5*(U1 - U2) * np.tanh(Y / delta)
    V_mean = np.zeros_like(X)
    
    # 添加扰动
    np.random.seed(456)
    k_unstable = 1.0 / delta
    u_pert = 0.05 * np.sin(k_unstable * X) * np.exp(-Y**2 / delta**2)
    v_pert = 0.05 * np.cos(k_unstable * X) * np.exp(-Y**2 / delta**2)
    
    u = U_mean + u_pert
    v = V_mean + v_pert
    
    # 转换到谱空间
    u_hat = np.fft.fft2(u)
    v_hat = np.fft.fft2(v)
    
    # 投影到无散空间
    def project_divergence_free(u_hat, v_hat):
        div_hat = 1j * KX * u_hat + 1j * KY * v_hat
        p_hat = div_hat / (1j * k2)
        u_hat_new = u_hat - 1j * KX * p_hat
        v_hat_new = v_hat - 1j * KY * p_hat
        return u_hat_new, v_hat_new
    
    u_hat, v_hat = project_divergence_free(u_hat, v_hat)
    
    # 物理参数
    nu = 0.005
    
    # 存储快照
    snapshots = [0, nt//4, nt//2, 3*nt//4, nt-1]
    u_snapshots = []
    v_snapshots = []
    
    # 时间推进
    for n in range(nt):
        # 存储快照
        if n in snapshots:
            u_snapshots.append(np.real(np.fft.ifft2(u_hat)))
            v_snapshots.append(np.real(np.fft.ifft2(v_hat)))
        
        # 物理空间速度
        u = np.real(np.fft.ifft2(u_hat))
        v = np.real(np.fft.ifft2(v_hat))
        
        # 对流项
        dudx = np.real(np.fft.ifft2(1j * KX * u_hat))
        dudy = np.real(np.fft.ifft2(1j * KY * u_hat))
        dvdx = np.real(np.fft.ifft2(1j * KX * v_hat))
        dvdy = np.real(np.fft.ifft2(1j * KY * v_hat))
        
        conv_u = u * dudx + v * dudy
        conv_v = u * dvdx + v * dvdy
        
        conv_u_hat = np.fft.fft2(conv_u)
        conv_v_hat = np.fft.fft2(conv_v)
        
        # 时间推进
        u_hat = u_hat - dt * conv_u_hat - dt * nu * k2 * u_hat
        v_hat = v_hat - dt * conv_v_hat - dt * nu * k2 * v_hat
        
        # 投影
        u_hat, v_hat = project_divergence_free(u_hat, v_hat)
    
    # 创建图形
    fig = plt.figure(figsize=(16, 10))
    gs = GridSpec(2, 3, figure=fig, hspace=0.3, wspace=0.3)
    
    # 绘制不同时刻的涡量
    for idx, (i, u_snap, v_snap) in enumerate(zip(snapshots, u_snapshots, v_snapshots)):
        if idx >= 5:
            break
        ax = fig.add_subplot(gs[idx//3, idx%3])
        
        omega = np.gradient(v_snap, dx, axis=1) - np.gradient(u_snap, dy, axis=0)
        im = ax.contourf(X, Y, omega, levels=20, cmap='RdBu_r', vmin=-3, vmax=3)
        plt.colorbar(im, ax=ax, label='omega')
        ax.set_xlabel('x', fontsize=10)
        ax.set_ylabel('y', fontsize=10)
        ax.set_title(f't = {t[i]:.2f}', fontsize=11, fontweight='bold')
        ax.set_aspect('equal')
    
    plt.savefig(f'{output_dir}/case4_mixing_layer.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"  混合层厚度: {delta}")
    print(f"  速度差: {U1 - U2}")
    print(f"  模拟时间: {t[-1]:.2f}")
    
    return {'snapshots': snapshots, 'u_snap': u_snapshots, 'v_snap': v_snapshots}

result_case4 = mixing_layer_instability()
print("✓ 案例4完成")

# ==================== 案例5：LES性能分析 ====================
print("\n【案例5：LES性能与分辨率分析】")

def les_resolution_analysis():
    """分析LES的分辨率要求和计算成本"""
    
    # 雷诺数范围
    Re_range = np.logspace(3, 6, 50)
    
    # DNS网格要求: N ~ Re^(9/4)
    N_dns = Re_range**(9/4) / 1e6
    
    # LES网格要求: N ~ Re^(3/2) (壁面模型LES)
    N_les_wm = Re_range**(3/2) / 1e6
    
    # LES网格要求: N ~ Re^(2) (壁面解析LES)
    N_les_wr = Re_range**(2) / 1e6
    
    # RANS网格: 与Re无关，取决于几何
    N_rans = np.ones_like(Re_range) * 0.1
    
    # 计算时间估计（相对单位）
    T_dns = N_dns * Re_range**(3/4) / 1e3
    T_les_wm = N_les_wm * Re_range**(1/2) / 1e3
    T_les_wr = N_les_wr * Re_range**(1/2) / 1e3
    T_rans = N_rans * 100 / 1e3
    
    # 创建图形
    fig = plt.figure(figsize=(16, 10))
    gs = GridSpec(2, 3, figure=fig, hspace=0.3, wspace=0.3)
    
    # 子图1：网格数对比
    ax1 = fig.add_subplot(gs[0, 0])
    ax1.loglog(Re_range, N_dns, 'b-', linewidth=2.5, label='DNS')
    ax1.loglog(Re_range, N_les_wr, 'r--', linewidth=2.5, label='LES (Wall-Resolved)')
    ax1.loglog(Re_range, N_les_wm, 'g:', linewidth=2.5, label='LES (Wall-Modeled)')
    ax1.loglog(Re_range, N_rans, 'k-.', linewidth=2.5, label='RANS')
    ax1.set_xlabel('Reynolds Number', fontsize=11)
    ax1.set_ylabel('Grid Points (Millions)', fontsize=11)
    ax1.set_title('Grid Requirements vs Reynolds Number', fontsize=12, fontweight='bold')
    ax1.legend(fontsize=9)
    ax1.grid(True, alpha=0.3, which='both')
    
    # 子图2：计算时间对比
    ax2 = fig.add_subplot(gs[0, 1])
    ax2.loglog(Re_range, T_dns, 'b-', linewidth=2.5, label='DNS')
    ax2.loglog(Re_range, T_les_wr, 'r--', linewidth=2.5, label='LES (Wall-Resolved)')
    ax2.loglog(Re_range, T_les_wm, 'g:', linewidth=2.5, label='LES (Wall-Modeled)')
    ax2.loglog(Re_range, T_rans, 'k-.', linewidth=2.5, label='RANS')
    ax2.set_xlabel('Reynolds Number', fontsize=11)
    ax2.set_ylabel('Relative Compute Time', fontsize=11)
    ax2.set_title('Computational Cost vs Reynolds Number', fontsize=12, fontweight='bold')
    ax2.legend(fontsize=9)
    ax2.grid(True, alpha=0.3, which='both')
    
    # 子图3：特定Re下的对比
    ax3 = fig.add_subplot(gs[0, 2])
    Re_target = 1e5
    methods = ['RANS', 'LES-WM', 'LES-WR', 'DNS']
    grids = [0.1, Re_target**(3/2)/1e6, Re_target**(2)/1e6, Re_target**(9/4)/1e6]
    colors = ['green', 'orange', 'red', 'purple']
    
    x_pos = np.arange(len(methods))
    bars = ax3.bar(x_pos, grids, color=colors, alpha=0.7, edgecolor='black', linewidth=1.5)
    ax3.set_ylabel('Grid Points (Millions)', fontsize=11)
    ax3.set_title(f'Grid Comparison @ Re={Re_target:.0e}', fontsize=12, fontweight='bold')
    ax3.set_xticks(x_pos)
    ax3.set_xticklabels(methods)
    ax3.set_yscale('log')
    ax3.grid(True, alpha=0.3, axis='y')
    for bar, val in zip(bars, grids):
        ax3.text(bar.get_x() + bar.get_width()/2, bar.get_height()*1.2, 
                f'{val:.1e}', ha='center', fontsize=9, fontweight='bold')
    
    # 子图4：分辨率要求（y+）
    ax4 = fig.add_subplot(gs[1, 0])
    y_plus_dns = np.ones_like(Re_range) * 0.5
    y_plus_les_wr = np.ones_like(Re_range) * 1.0
    y_plus_les_wm = np.ones_like(Re_range) * 30
    y_plus_rans = np.ones_like(Re_range) * 30
    
    ax4.semilogx(Re_range, y_plus_dns, 'b-', linewidth=2.5, label='DNS')
    ax4.semilogx(Re_range, y_plus_les_wr, 'r--', linewidth=2.5, label='LES-WR')
    ax4.semilogx(Re_range, y_plus_les_wm, 'g:', linewidth=2.5, label='LES-WM')
    ax4.semilogx(Re_range, y_plus_rans, 'k-.', linewidth=2.5, label='RANS')
    ax4.fill_between(Re_range, 0, 5, alpha=0.2, color='blue', label='Viscous Sublayer')
    ax4.fill_between(Re_range, 5, 30, alpha=0.2, color='orange', label='Buffer Layer')
    ax4.fill_between(Re_range, 30, 500, alpha=0.2, color='green', label='Log Layer')
    ax4.set_xlabel('Reynolds Number', fontsize=11)
    ax4.set_ylabel('First Grid Point y+', fontsize=11)
    ax4.set_title('Wall Resolution Requirements', fontsize=12, fontweight='bold')
    ax4.legend(fontsize=8, loc='upper right')
    ax4.grid(True, alpha=0.3)
    ax4.set_ylim([0, 100])
    
    # 子图5：适用性图表
    ax5 = fig.add_subplot(gs[1, 1])
    ax5.set_xlim(0, 10)
    ax5.set_ylim(0, 10)
    ax5.axis('off')
    ax5.set_title('Method Selection Guide', fontsize=12, fontweight='bold')
    
    rect_rans = Rectangle((0.5, 7), 4, 2.5, fill=True, facecolor='lightgreen', 
                          edgecolor='green', linewidth=2, alpha=0.7)
    ax5.add_patch(rect_rans)
    ax5.text(2.5, 8.25, 'RANS', ha='center', fontsize=14, fontweight='bold', color='green')
    ax5.text(2.5, 7.5, 'Steady flows\nSimple geometry\nQuick results', ha='center', fontsize=9)
    
    rect_les = Rectangle((5.5, 7), 4, 2.5, fill=True, facecolor='lightyellow', 
                         edgecolor='orange', linewidth=2, alpha=0.7)
    ax5.add_patch(rect_les)
    ax5.text(7.5, 8.25, 'LES', ha='center', fontsize=14, fontweight='bold', color='orange')
    ax5.text(7.5, 7.5, 'Unsteady flows\nComplex physics\nHigh fidelity', ha='center', fontsize=9)
    
    rect_dns = Rectangle((0.5, 3.5), 4, 2.5, fill=True, facecolor='lightcoral', 
                         edgecolor='red', linewidth=2, alpha=0.7)
    ax5.add_patch(rect_dns)
    ax5.text(2.5, 4.75, 'DNS', ha='center', fontsize=14, fontweight='bold', color='red')
    ax5.text(2.5, 4.0, 'Low Re only\nResearch purpose\nNo modeling', ha='center', fontsize=9)
    
    rect_hybrid = Rectangle((5.5, 3.5), 4, 2.5, fill=True, facecolor='lightblue', 
                            edgecolor='blue', linewidth=2, alpha=0.7)
    ax5.add_patch(rect_hybrid)
    ax5.text(7.5, 4.75, 'Hybrid RANS-LES', ha='center', fontsize=12, fontweight='bold', color='blue')
    ax5.text(7.5, 4.0, 'Best of both\nDES, SAS, PANS\nIndustrial use', ha='center', fontsize=9)
    
    ax5.text(5, 1.5, 'Decision Flow:', ha='center', fontsize=11, fontweight='bold')
    ax5.text(5, 0.8, 'Need unsteady? -> LES\nLow Re & simple? -> DNS\nQuick results? -> RANS\nComplex + affordable? -> Hybrid', 
            ha='center', fontsize=9, family='monospace')
    
    # 子图6：LES最佳实践
    ax6 = fig.add_subplot(gs[1, 2])
    ax6.axis('off')
    best_practice = """
    LES Best Practices:
    
    1. GRID RESOLUTION:
       - Wall-resolved: Dx+ < 50, Dy+ < 1
       - Wall-modeled: Dx+ < 200
       - Isotropic cells in bulk flow
    
    2. TIME STEP:
       - CFL < 1 (preferably < 0.5)
       - Resolve turbulence time scales
    
    3. SGS MODEL:
       - Wall-resolved: Dynamic
       - Wall-modeled: WALE or WMLES
       - Avoid Smagorinsky near walls
    
    4. STATISTICS:
       - Run 10-20 flow-through times
       - Check convergence of 2nd moments
    
    5. VALIDATION:
       - Compare with DNS when possible
       - Validate mean flow and fluctuations
       - Check energy spectrum
    """
    ax6.text(0.5, 0.5, best_practice, transform=ax6.transAxes, fontsize=9,
            verticalalignment='center', horizontalalignment='center',
            bbox=dict(boxstyle='round', facecolor='lightyellow', alpha=0.8))
    
    plt.savefig(f'{output_dir}/case5_resolution_analysis.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print(f"  DNS网格数 @ Re=1e5: {1e5**(9/4)/1e6:.1f}M")
    print(f"  LES-WR网格数 @ Re=1e5: {1e5**(2)/1e6:.1f}M")
    print(f"  LES-WM网格数 @ Re=1e5: {1e5**(3/2)/1e6:.1f}M")
    print(f"  RANS网格数 @ Re=1e5: 0.1M")
    
    return {'Re': Re_range, 'N_dns': N_dns, 'N_les_wr': N_les_wr, 'N_les_wm': N_les_wm}

result_case5 = les_resolution_analysis()
print("✓ 案例5完成")

# ==================== GIF动画创建 ====================
print("\n【创建GIF动画】")

def create_animation():
    try:
        fig, axes = plt.subplots(2, 2, figsize=(14, 10))
        fig.suptitle('LES Simulation: Decaying 2D Turbulence', fontsize=14, fontweight='bold')
        
        # 使用案例3的数据
        t_data = result_case3['t']
        ke_data = result_case3['kinetic_energy']
        
        # 重新运行简化的模拟来获取动画帧
        nx, ny = 64, 64
        Lx, Ly = 2*np.pi, 2*np.pi
        dx = Lx / nx
        dy = Ly / ny
        x = np.linspace(0, Lx, nx, endpoint=False)
        y = np.linspace(0, Ly, ny, endpoint=False)
        X, Y = np.meshgrid(x, y)
        
        kx = 2*np.pi * np.fft.fftfreq(nx, dx)
        ky = 2*np.pi * np.fft.fftfreq(ny, dy)
        KX, KY = np.meshgrid(kx, ky)
        k2 = KX**2 + KY**2
        k2[0, 0] = 1e-10
        k_mag = np.sqrt(k2)
        
        # 初始条件
        np.random.seed(123)
        u_hat = np.zeros((ny, nx), dtype=complex)
        v_hat = np.zeros((ny, nx), dtype=complex)
        
        for i in range(1, 6):
            for j in range(1, 6):
                k_mag_val = np.sqrt(i**2 + j**2)
                amp = 1.0 / k_mag_val**2
                phase_u = np.random.uniform(0, 2*np.pi)
                phase_v = np.random.uniform(0, 2*np.pi)
                u_hat[j, i] += amp * np.exp(1j * phase_u)
                v_hat[-j, -i] += amp * np.exp(-1j * phase_u)
                v_hat[j, i] += amp * np.exp(1j * phase_v)
                v_hat[-j, -i] += amp * np.exp(-1j * phase_v)
        
        # 投影到无散空间
        def project_divergence_free(u_hat, v_hat):
            div_hat = 1j * KX * u_hat + 1j * KY * v_hat
            p_hat = div_hat / (1j * k2)
            u_hat_new = u_hat - 1j * KX * p_hat
            v_hat_new = v_hat - 1j * KY * p_hat
            return u_hat_new, v_hat_new
        
        u_hat, v_hat = project_divergence_free(u_hat, v_hat)
        
        # 归一化
        u = np.real(np.fft.ifft2(u_hat))
        scale = 1.0 / (np.max(np.abs(u)) + 1e-10)
        u_hat *= scale
        v_hat *= scale
        
        nu = 0.02
        dt = 0.02
        
        n_frames = 30
        
        def init():
            for ax in axes.flat:
                ax.clear()
            return []
        
        def update(frame):
            nonlocal u_hat, v_hat
            
            # 物理空间速度
            u = np.real(np.fft.ifft2(u_hat))
            v = np.real(np.fft.ifft2(v_hat))
            
            # 对流项
            dudx = np.real(np.fft.ifft2(1j * KX * u_hat))
            dudy = np.real(np.fft.ifft2(1j * KY * u_hat))
            dvdx = np.real(np.fft.ifft2(1j * KX * v_hat))
            dvdy = np.real(np.fft.ifft2(1j * KY * v_hat))
            
            conv_u = u * dudx + v * dudy
            conv_v = u * dvdx + v * dvdy
            
            conv_u_hat = np.fft.fft2(conv_u)
            conv_v_hat = np.fft.fft2(conv_v)
            
            # 时间推进
            u_hat = u_hat - dt * conv_u_hat - dt * nu * k2 * u_hat
            v_hat = v_hat - dt * conv_v_hat - dt * nu * k2 * v_hat
            
            # 投影
            u_hat, v_hat = project_divergence_free(u_hat, v_hat)
            
            # 更新显示
            current_time = frame * dt
            
            # 子图1：速度场
            ax1 = axes[0, 0]
            ax1.clear()
            speed = np.sqrt(u**2 + v**2)
            ax1.contourf(X, Y, speed, levels=15, cmap='viridis', vmin=0, vmax=1)
            ax1.set_title(f'Velocity Field @ t={current_time:.2f}')
            ax1.set_xlabel('x')
            ax1.set_ylabel('y')
            ax1.set_aspect('equal')
            
            # 子图2：涡量场
            ax2 = axes[0, 1]
            ax2.clear()
            omega = np.gradient(v, dx, axis=1) - np.gradient(u, dy, axis=0)
            ax2.contourf(X, Y, omega, levels=15, cmap='RdBu_r', vmin=-2, vmax=2)
            ax2.set_title('Vorticity Field')
            ax2.set_xlabel('x')
            ax2.set_ylabel('y')
            ax2.set_aspect('equal')
            
            # 子图3：动能衰减
            ax3 = axes[1, 0]
            ax3.clear()
            ke = 0.5 * np.mean(u**2 + v**2)
            ax3.text(0.5, 0.5, f'Kinetic Energy\n{ke:.4f}', transform=ax3.transAxes,
                    ha='center', va='center', fontsize=16)
            ax3.set_xlim([0, 1])
            ax3.set_ylim([0, 1])
            ax3.axis('off')
            
            # 子图4：能谱
            ax4 = axes[1, 1]
            ax4.clear()
            E_k = 0.5 * (np.abs(u_hat)**2 + np.abs(v_hat)**2)
            
            k_bins = np.linspace(0, np.max(k_mag), 15)
            E_k_avg = np.zeros(len(k_bins)-1)
            for i in range(len(k_bins)-1):
                mask = (k_mag >= k_bins[i]) & (k_mag < k_bins[i+1])
                if np.sum(mask) > 0:
                    E_k_avg[i] = np.mean(E_k[mask])
            
            k_centers = 0.5 * (k_bins[:-1] + k_bins[1:])
            ax4.loglog(k_centers[1:], E_k_avg[1:], 'b-o', markersize=4)
            ax4.set_xlabel('Wavenumber k')
            ax4.set_ylabel('E(k)')
            ax4.set_title('Energy Spectrum')
            ax4.grid(True, alpha=0.3, which='both')
            
            return []
        
        anim = FuncAnimation(fig, update, frames=n_frames, init_func=init, blit=False)
        writer = PillowWriter(fps=5)
        anim.save(f'{output_dir}/case_les_animation.gif', writer=writer)
        plt.close()
        
        print("✓ GIF动画创建成功")
        return True
    except Exception as e:
        print(f"  GIF动画创建失败: {e}")
        import traceback
        traceback.print_exc()
        return False

animation_success = create_animation()

# ==================== 仿真总结 ====================
print("\n" + "="*70)
print("仿真完成总结")
print("="*70)
print("\n【仿真结果文件】")
print(f"  1. case1_spectrum_filtering.png - 湍流能谱与滤波函数")
print(f"  2. case2_sgs_models.png - 亚格子模型对比")
print(f"  3. case3_decaying_turbulence.png - 2D衰减湍流模拟")
print(f"  4. case4_mixing_layer.png - 混合层不稳定性")
print(f"  5. case5_resolution_analysis.png - LES分辨率分析")
if animation_success:
    print(f"  6. case_les_animation.gif - LES湍流动画")

print("\n【关键结果】")
print(f"  截断波数 k_c: {result_case1['k_c']}")
print(f"  解析能量比例: {result_case1['energy_ratio']*100:.1f}%")
print(f"  Smagorinsky平均涡粘: {np.mean(result_case2['nu_smag']):.4f}")
print(f"  WALE平均涡粘: {np.mean(result_case2['nu_wale']):.4f}")
print(f"  衰减湍流初始KE: {result_case3['kinetic_energy'][0]:.4f}")
print(f"  衰减湍流最终KE: {result_case3['kinetic_energy'][-1]:.4f}")
print(f"  衰减比例: {result_case3['kinetic_energy'][-1]/result_case3['kinetic_energy'][0]:.4f}")
print(f"  DNS网格数 @ Re=1e6: {1e6**(9/4)/1e6:.0f}M")
print(f"  LES-WR网格数 @ Re=1e6: {1e6**(2)/1e6:.0f}M")
print(f"  LES-WM网格数 @ Re=1e6: {1e6**(3/2)/1e6:.0f}M")

print("\n【LES要点】")
print("  1. 滤波分离大尺度（解析）和小尺度（建模）")
print("  2. Smagorinsky模型简单但过于耗散")
print("  3. 动态模型自适应但计算成本高")
print("  4. WALE模型近壁行为更好")
print("  5. 壁面解析LES需要 y+ < 1")
print("  6. 壁面模型LES可大幅降低计算成本")

print("\n" + "="*70)
print("所有仿真案例已完成！")
print("="*70)

AtomGit开源社区

AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。

更多推荐

为搜索引擎 Agent 设计 Harness 查询改写与扩展

本文将带你深入探讨如何为搜索引擎 Agent 设计一套强大的Harness（框架），专门用于查询改写（Query Rewriting）与查询扩展（Query Expansion）。我们将从理论基础出发，一步步构建一个实用的系统，包括同义词扩展、意图识别、错误纠正、上下文理解等核心功能。架构设计：我们设计了一个从预处理到候选生成再到排序的完整流程查询预处理：实现了文本清洗、分词、去停用词和词形还原查

AtomGit开源社区

2026年川菜文化推广大使，你猜是谁？

AtomGit开源社区

基于SpringBoot+Vue的在线家具商城设计与实现管理系统设计与实现【Java+MySQL+MyBatis完整源码】

AtomGit开源社区

所有评论(0)

查看更多评论

kkchenjj

@weixin_42749425

已为社区贡献264条内容

对流换热仿真-主题061_大涡模拟LES-大涡模拟LES

kkchenjj

主题061：大涡模拟（LES）

一、引言

1.1 湍流模拟的挑战

1.2 雷诺平均与大涡模拟

1.3 LES的优势与应用领域

二、大涡模拟的基本原理

2.1 滤波概念

2.1.1 盒式滤波器

2.1.2 高斯滤波器

2.1.3 谱截断滤波器

2.2 滤波后的控制方程

2.3 亚格子应力与热流

2.4 亚格子应力的分解

三、亚格子尺度模型

3.1 涡粘模型

3.1.1 Smagorinsky模型

3.1.2 动态Smagorinsky模型

3.1.3 WALE模型

3.2 结构函数模型

3.2.1 尺度相似模型

3.2.2 混合模型

3.3 一方程模型

3.3.1 亚格子动能模型

3.4 亚格子热流模型

3.4.1 梯度扩散模型

3.4.2 动态热流模型

四、近壁面处理

4.1 壁面附近湍流特性

4.2 壁面解析LES

4.3 壁面模型LES（WMLES）

4.3.1 对数律壁面模型

4.3.2 代数壁面模型

4.3.3 双层模型

五、LES的数值方法

5.1 空间离散格式

5.1.1 高阶有限差分

5.1.2 谱方法

5.2 时间推进方法

5.2.1 显式方法

5.2.2 隐式方法

5.3 压力-速度耦合

5.3.1 投影法

5.3.2 SIMPLE算法

5.4 边界条件

5.4.1 入口边界

5.4.2 出口边界

5.4.3 周期性边界

5.5 计算网格要求

5.5.1 网格分辨率准则

5.5.2 网格质量

六、LES的验证与确认

6.1 统计收敛性

6.2 网格无关性验证

6.3 与实验对比

6.3.1 平均量对比

6.3.2 瞬态特性对比

6.4 与DNS/RANS对比

所有评论(0)

温馨提示：您尚未绑定手机号

kkchenjj