结构风振仿真-主题021-CFD数值模拟基础

kkchenkx

338人浏览 · 2026-05-09 21:17:17

kkchenkx · 2026-05-09 21:17:17 发布

第21篇：CFD数值模拟基础

摘要

计算流体力学（CFD）是现代风工程研究的核心工具，通过数值方法求解流体运动的基本方程，实现对复杂风场的精确预测。本主题系统介绍CFD的基本原理和数值方法，包括控制方程的离散化、湍流模型、边界条件处理、网格生成技术等核心内容。通过Python程序实现简单的CFD求解器，演示流场计算的基本流程，帮助读者建立CFD数值模拟的理论基础和编程能力。

关键词

计算流体力学；CFD；Navier-Stokes方程；有限体积法；湍流模型；边界条件；网格生成；数值离散；风场模拟

在这里插入图片描述

1. 引言

1.1 CFD在风工程中的重要性

计算流体力学（Computational Fluid Dynamics, CFD）是利用数值方法求解流体流动问题的学科。在结构风工程领域，CFD已经成为继理论分析和风洞试验之后的第三种重要研究手段。

CFD的优势：

成本效益：相比风洞试验，CFD可以显著降低研究成本。一个复杂的CFD模拟成本通常仅为风洞试验的10%-20%。

信息丰富：CFD可以提供流场的完整信息，包括速度、压力、湍流特性等在整个计算域内的分布，而风洞试验通常只能在有限测点获得数据。

灵活性：CFD可以轻松改变几何形状、来流条件和物理参数，便于进行参数化研究和优化设计。

极端条件：CFD可以模拟风洞试验难以实现的极端条件，如超高层建筑、复杂地形等。

CFD的局限性：

模型依赖：CFD结果的准确性高度依赖于湍流模型、壁面函数等物理模型的选择。

计算资源：高保真度的CFD模拟需要大量的计算资源，包括高性能计算机和长时间计算。

验证需求：CFD结果需要通过风洞试验或现场实测进行验证，以确保可靠性。

1.2 CFD发展历程

CFD的发展经历了几个重要阶段：

早期阶段（1960-1970年代）：

基于势流理论的简化方法
边界层方程的数值求解
计算资源极其有限

发展阶段（1980-1990年代）：

雷诺平均Navier-Stokes（RANS）方法的成熟
有限体积法成为主流
商业CFD软件的出现（如FLUENT、STAR-CD）

成熟阶段（2000-2010年代）：

大涡模拟（LES）的应用
高性能计算的普及
多物理场耦合能力增强

现代阶段（2010年代至今）：

直接数值模拟（DNS）在简单流动中的应用
机器学习与CFD的结合
云计算和GPU加速

1.3 CFD在风工程中的应用

建筑风环境：

高层建筑周围的风场分布
行人风环境评估
建筑群体干扰效应

结构风荷载：

复杂体型建筑的风压分布
大跨度屋盖的风荷载
桥梁的风致响应

风能利用：

风力机气动性能分析
风电场布局优化
尾流效应研究

环境工程：

污染物扩散
城市热岛效应
自然通风设计

2. 控制方程

2.1 连续性方程

连续性方程表达了质量守恒定律：

$∂ρ∂t+∇⋅(ρu)=0\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{u}) = 0$

对于不可压缩流动（密度 $ρ\rho$ 为常数），连续性方程简化为：

$∇⋅u=0\nabla \cdot \mathbf{u} = 0$

或展开为：

$∂u∂x+∂v∂y+∂w∂z=0\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0$

其中 $u=(u,v,w)\mathbf{u} = (u, v, w)$ 为速度矢量。

2.2 Navier-Stokes方程

Navier-Stokes方程是描述流体运动动量守恒的基本方程：

$∂(ρu)∂t+∇⋅(ρuu)=−∇p+∇⋅τ+ρg\frac{\partial (\rho \mathbf{u})}{\partial t} + \nabla \cdot (\rho \mathbf{u} \mathbf{u}) = -\nabla p + \nabla \cdot \boldsymbol{\tau} + \rho \mathbf{g}$

其中：

$p$ ：压力
$τ\boldsymbol{\tau}$ ：黏性应力张量
$g\mathbf{g}$ ：重力加速度矢量

对于牛顿流体，黏性应力张量为：

$τ=μ[∇u+(∇u)T−23(∇⋅u)I]\boldsymbol{\tau} = \mu \left[\nabla \mathbf{u} + (\nabla \mathbf{u})^T - \frac{2}{3}(\nabla \cdot \mathbf{u})\mathbf{I}\right]$

其中 $μ\mu$ 为动力黏度， $I\mathbf{I}$ 为单位张量。

对于不可压缩流动，Navier-Stokes方程简化为：

$∂u∂t+(u⋅∇)u=−1ρ∇p+ν∇2u+g\frac{\partial \mathbf{u}}{\partial t} + (\mathbf{u} \cdot \nabla)\mathbf{u} = -\frac{1}{\rho}\nabla p + \nu \nabla^2 \mathbf{u} + \mathbf{g}$

其中 $ν=μ/ρ\nu = \mu/\rho$ 为运动黏度。

2.3 能量方程

当需要考虑热效应时，需要求解能量方程：

$∂(ρe)∂t+∇⋅(ρeu)=−∇⋅q−p(∇⋅u)+τ:∇u+Se\frac{\partial (\rho e)}{\partial t} + \nabla \cdot (\rho e \mathbf{u}) = -\nabla \cdot \mathbf{q} - p(\nabla \cdot \mathbf{u}) + \boldsymbol{\tau} : \nabla \mathbf{u} + S_e$

其中：

$e$ ：单位质量内能
$q\mathbf{q}$ ：热通量矢量
$S_e$ ：能量源项

对于理想气体：

$e = c_v T$

其中 $c_v$ 为定容比热容， $T$ 为温度。

2.4 方程的守恒形式

守恒形式的控制方程可以统一表示为：

$∂ϕ∂t+∇⋅f=ψ\frac{\partial \boldsymbol{\phi}}{\partial t} + \nabla \cdot \mathbf{f} = \boldsymbol{\psi}$

其中：

$ϕ\boldsymbol{\phi}$ ：守恒变量矢量
$f\mathbf{f}$ ：通量矢量
$ψ\boldsymbol{\psi}$ ：源项矢量

对于质量、动量和能量守恒，有：

$ϕ=[ρρuρvρwρe]\boldsymbol{\phi} = \begin{bmatrix} \rho \\ \rho u \\ \rho v \\ \rho w \\ \rho e \end{bmatrix}$

守恒形式的重要性在于：

保证物理量的严格守恒
适用于激波等间断问题
便于有限体积法离散

3. 湍流模型

3.1 湍流的基本特征

湍流是流体力学中最复杂的现象之一，具有以下特征：

随机性：湍流运动具有高度的随机性，速度和压力等物理量呈现不规则的脉动。

三维性：湍流本质上是三维的，即使在二维平均流动中也存在三维湍流结构。

涡旋结构：湍流由不同尺度的涡旋组成，大涡从平均流中获取能量，小涡通过黏性耗散能量。

能量级串：能量从大尺度涡旋向小尺度涡旋传递，形成能量级串过程。

雷诺数：湍流的产生与雷诺数密切相关，当 $Re > Re_{crit}$ 时，流动从层流转捩为湍流。

3.2 雷诺平均方法

雷诺平均Navier-Stokes（RANS）方法是目前工程应用中最广泛的湍流模拟方法。

雷诺分解：

将瞬时量分解为平均量和脉动量：

$\bar{u} + u'$

其中 $uˉ\bar{u}$ 为时均速度， $u^{'}$ 为脉动速度。

雷诺平均方程：

对Navier-Stokes方程进行雷诺平均，得到：

$∂uˉi∂t+uˉj∂uˉi∂xj=−1ρ∂pˉ∂xi+ν∂2uˉi∂xj∂xj−∂∂xj(ui′uj′‾)\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' u_j'})$

其中 $−ρui′uj′‾-\rho\overline{u_i' u_j'}$ 为雷诺应力，代表湍流对平均流的影响。

封闭问题：

雷诺应力是未知量，需要通过湍流模型来封闭方程组。

3.3 k-ε湍流模型

k-ε模型是最常用的两方程湍流模型，通过求解湍动能k和湍流耗散率ε来封闭雷诺应力。

湍动能方程：

$∂k∂t+uˉj∂k∂xj=Pk−ε+∂∂xj[(ν+νtσk)∂k∂xj]\frac{\partial k}{\partial t} + \bar{u}_j \frac{\partial k}{\partial x_j} = P_k - \varepsilon + \frac{\partial}{\partial x_j}\left[\left(\nu + \frac{\nu_t}{\sigma_k}\right)\frac{\partial k}{\partial x_j}\right]$

耗散率方程：

$∂ε∂t+uˉj∂ε∂xj=Cε1εkPk−Cε2ε2k+∂∂xj[(ν+νtσε)∂ε∂xj]\frac{\partial \varepsilon}{\partial t} + \bar{u}_j \frac{\partial \varepsilon}{\partial x_j} = C_{\varepsilon 1}\frac{\varepsilon}{k}P_k - C_{\varepsilon 2}\frac{\varepsilon^2}{k} + \frac{\partial}{\partial x_j}\left[\left(\nu + \frac{\nu_t}{\sigma_\varepsilon}\right)\frac{\partial \varepsilon}{\partial x_j}\right]$

其中：

$P_k$ ：湍动能生成项
$νt=Cμk2ε\nu_t = C_\mu \frac{k^2}{\varepsilon}$ ：湍流黏度

标准模型常数：

$Cμ=0.09,Cε1=1.44,Cε2=1.92,σk=1.0,σε=1.3C_\mu = 0.09, \quad C_{\varepsilon 1} = 1.44, \quad C_{\varepsilon 2} = 1.92, \quad \sigma_k = 1.0, \quad \sigma_\varepsilon = 1.3$

雷诺应力计算：

$−ui′uj′‾=νt(∂uˉi∂xj+∂uˉj∂xi)−23kδij-\overline{u_i' u_j'} = \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}$

3.4 k-ω SST模型

k-ω SST（Shear Stress Transport）模型结合了k-ε和k-ω模型的优点，在壁面附近使用k-ω模型，在远场使用k-ε模型。

湍动能方程：

$∂k∂t+uˉj∂k∂xj=Pk−β∗ωk+∂∂xj[(ν+σkνt)∂k∂xj]\frac{\partial k}{\partial t} + \bar{u}_j \frac{\partial k}{\partial x_j} = P_k - \beta^* \omega k + \frac{\partial}{\partial x_j}\left[\left(\nu + \sigma_k \nu_t\right)\frac{\partial k}{\partial x_j}\right]$

比耗散率方程：

$∂ω∂t+uˉj∂ω∂xj=γνtPk−βω2+∂∂xj[(ν+σωνt)∂ω∂xj]+2(1−F1)σω2ω∂k∂xj∂ω∂xj\frac{\partial \omega}{\partial t} + \bar{u}_j \frac{\partial \omega}{\partial x_j} = \frac{\gamma}{\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}$

其中 $ω=ε/(Cμk)\omega = \varepsilon/(C_\mu k)$ 为比耗散率， $F_1$ 为混合函数。

SST模型的优势：

在逆压梯度流动中表现更好
对分离流的预测更准确
在风工程中应用广泛

3.5 大涡模拟（LES）

大涡模拟直接求解大尺度涡旋，对小尺度涡旋采用亚格子模型。

滤波操作：

$fˉ(x,t)=∫Df(x′,t)G(x−x′)dx′\bar{f}(\mathbf{x}, t) = \int_D f(\mathbf{x}', t) G(\mathbf{x} - \mathbf{x}') d\mathbf{x}'$

其中 $G$ 为滤波函数。

滤波后的Navier-Stokes方程：

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

其中 $τij=uiuj‾−uˉiuˉj\tau_{ij} = \overline{u_i u_j} - \bar{u}_i \bar{u}_j$ 为亚格子应力。

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)$ ：滤波应变率张量
$Cs≈0.1C_s \approx 0.1$ ：Smagorinsky常数
$Δ\Delta$ ：滤波宽度

LES在风工程中的应用：

建筑尾流分析
湍流风场生成
风致振动研究

4. 数值离散方法

4.1 有限差分法

有限差分法（FDM）用差商代替导数，直接在网格点上离散控制方程。

一阶导数近似：

前向差分：

$(∂u∂x)i≈ui+1−uiΔx\left(\frac{\partial u}{\partial x}\right)_i \approx \frac{u_{i+1} - u_i}{\Delta x}$

后向差分：

$(∂u∂x)i≈ui−ui−1Δx\left(\frac{\partial u}{\partial x}\right)_i \approx \frac{u_i - u_{i-1}}{\Delta x}$

中心差分：

$(∂u∂x)i≈ui+1−ui−12Δx\left(\frac{\partial u}{\partial x}\right)_i \approx \frac{u_{i+1} - u_{i-1}}{2\Delta x}$

二阶导数近似：

$(∂2u∂x2)i≈ui+1−2ui+ui−1Δx2\left(\frac{\partial^2 u}{\partial x^2}\right)_i \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{\Delta x^2}$

截断误差：

一阶差分： $O(Δx)O(\Delta x)$
中心差分： $O(Δx2)O(\Delta x^2)$

4.2 有限体积法

有限体积法（FVM）是目前CFD中最常用的离散方法，基于控制体的积分形式守恒方程。

基本思想：

将计算域划分为若干控制体，在每个控制体上积分控制方程：

$∫V∂ϕ∂tdV+∫Sf⋅dS=∫VψdV\int_V \frac{\partial \phi}{\partial t} dV + \int_S \mathbf{f} \cdot d\mathbf{S} = \int_V \psi dV$

离散形式：

$ϕPn+1−ϕPnΔtVP+∑facesff⋅Sf=ψPVP\frac{\phi_P^{n+1} - \phi_P^n}{\Delta t} V_P + \sum_{faces} \mathbf{f}_f \cdot \mathbf{S}_f = \psi_P V_P$

其中下标 $P$ 表示控制体中心， $f$ 表示控制面。

FVM的优势：

严格保证守恒性
适用于任意网格
物理意义清晰

4.3 时间离散方法

显式格式：

欧拉显式：

$ϕn+1=ϕn+Δt⋅f(ϕn)\phi^{n+1} = \phi^n + \Delta t \cdot f(\phi^n)$

优点：计算简单，易于并行
缺点：时间步长受稳定性限制

隐式格式：

欧拉隐式：

$ϕn+1=ϕn+Δt⋅f(ϕn+1)\phi^{n+1} = \phi^n + \Delta t \cdot f(\phi^{n+1})$

优点：无条件稳定
缺点：需要求解方程组

Crank-Nicolson格式：

$ϕn+1=ϕn+Δt2[f(ϕn)+f(ϕn+1)]\phi^{n+1} = \phi^n + \frac{\Delta t}{2} \left[f(\phi^n) + f(\phi^{n+1})\right]$

二阶精度，无条件稳定。

4.4 对流项离散

对流项的离散对计算稳定性和精度有重要影响。

中心差分格式：

$\frac{\partial \phi}{\partial x} \approx u_e \frac{\phi_E - \phi_P}{\Delta x}$

二阶精度，但在高雷诺数下可能产生振荡。

迎风格式：

根据流速方向选择上游值：

$\frac{\partial \phi}{\partial x} \approx \begin{cases} u_e \frac{\phi_P - \phi_W}{\Delta x} & u_e > 0 \\ u_e \frac{\phi_E - \phi_P}{\Delta x} & u_e < 0 \end{cases}$

一阶精度，稳定性好但数值耗散大。

QUICK格式：

三阶迎风格式：

$ϕe=3ϕE+6ϕP−ϕW8\phi_e = \frac{3\phi_E + 6\phi_P - \phi_W}{8}$

精度高但可能产生过冲。

TVD格式：

总变差减小格式，结合高阶精度和稳定性。

5. 边界条件

5.1 入口边界条件

速度入口：

指定入口速度分布：

$u=uin(y,z)\mathbf{u} = \mathbf{u}_{in}(y, z)$

对于大气边界层流动，通常采用对数律或指数律剖面。

压力入口：

指定入口总压：

$p0=p+12ρu2p_0 = p + \frac{1}{2}\rho u^2$

湍流参数：

需要指定入口湍流强度 $I$ 和湍流尺度 $L$ ：

$kin=32(UrefI)2k_{in} = \frac{3}{2}(U_{ref} I)^2$

$εin=Cμ3/4k3/2L\varepsilon_{in} = C_\mu^{3/4} \frac{k^{3/2}}{L}$

5.2 出口边界条件

充分发展出口：

$∂ϕ∂x=0\frac{\partial \phi}{\partial x} = 0$

适用于下游无回流的情况。

压力出口：

指定出口静压：

$p = p_{out}$

对流出口：

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

其中 $U_c$ 为对流速度。

5.3 壁面边界条件

无滑移条件：

$u=0\mathbf{u} = 0$

壁面函数：

对于高雷诺数流动，采用壁面函数避免解析黏性底层。

对数律区：

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

其中：

$u+=u/uτu^+ = u/u_\tau$ ：无量纲速度
$y+=yuτ/νy^+ = y u_\tau / \nu$ ：无量纲距离
$κ=0.41\kappa = 0.41$ ：卡门常数
$C = 5.5$ ：常数

壁面剪切应力：

$τw=ρuτ2=ρuPu+uτ\tau_w = \rho u_\tau^2 = \rho \frac{u_P}{u^+} u_\tau$

5.4 对称边界条件

对称面上：

$∂ϕ∂n=0\frac{\partial \phi}{\partial n} = 0$

法向速度为零，切向速度梯度为零。

5.5 周期性边界条件

对于具有周期性的流动：

$ϕ(x)=ϕ(x+L)\phi(x) = \phi(x + L)$

常用于模拟无限长圆柱绕流等问题。

6. 网格生成技术

6.1 网格类型

结构化网格：

网格点排列规则，每个内部点有相同数量的相邻点。

优点：计算效率高，精度好
缺点：对复杂几何适应性差

非结构化网格：

网格点排列不规则，通常采用三角形或四面体单元。

优点：对复杂几何适应性强
缺点：计算开销大，需要更多内存

混合网格：

结合结构化和非结构化网格的优点，在关键区域使用结构化网格，在复杂区域使用非结构化网格。

6.2 网格质量

网格质量对CFD计算结果有重要影响。

关键指标：

长宽比：单元最长边与最短边之比，应小于100。

偏斜度：单元形状与理想形状的偏离程度，应小于0.85。

正交性：网格线与控制体面的正交程度，应大于0.1。

分辨率：

壁面附近： $y^+ < 1$ （解析黏性底层）或 $30 < y^+ < 300$ （使用壁面函数）
尾流区：需要足够分辨率捕捉涡旋结构
剪切层：需要加密网格

6.3 自适应网格

自适应网格根据流场特征动态调整网格分布。

自适应准则：

速度梯度
涡量大小
压力梯度
误差估计

自适应方法：

网格细化（h-refinement）：将单元细分为更小的单元。

网格粗化（h-coarsening）：合并小单元为大单元。

p-refinement：提高单元的多项式阶数。

7. 求解算法

7.1 SIMPLE算法

SIMPLE（Semi-Implicit Method for Pressure-Linked Equations）算法是求解不可压缩流动的经典算法。

基本步骤：

步骤1：速度预测

求解动量方程，获得中间速度 $u^*$ ：

$aeue∗=∑anbunb∗+b+(pP−pE)Aea_e u_e^* = \sum a_{nb} u_{nb}^* + b + (p_P - p_E)A_e$

步骤2：压力修正

求解压力修正方程：

$∇⋅(1ae∇p′)=∇⋅u∗\nabla \cdot \left(\frac{1}{a_e}\nabla p'\right) = \nabla \cdot \mathbf{u}^*$

步骤3：速度修正

$ue=ue∗+Aeae(pP′−pE′)u_e = u_e^* + \frac{A_e}{a_e}(p_P' - p_E')$

步骤4：压力更新

$p^* + \alpha_p p'$

其中 $αp\alpha_p$ 为压力松弛因子。

步骤5：迭代

重复步骤1-4直至收敛。

7.2 PISO算法

PISO（Pressure Implicit with Splitting of Operators）算法适用于瞬态流动计算。

特点：

不需要迭代
时间精度高
计算效率高

步骤：

预测步：求解动量方程
第一修正步：求解压力修正方程，修正速度和压力
第二修正步：再次求解压力修正方程

7.3 耦合求解

对于可压缩流动或强耦合问题，采用耦合求解方法。

全隐式耦合：

同时求解所有变量：

$[AuuAupApuApp][up]=[bubp]\begin{bmatrix} A_{uu} & A_{up} \\ A_{pu} & A_{pp} \end{bmatrix} \begin{bmatrix} \mathbf{u} \\ p \end{bmatrix} = \begin{bmatrix} \mathbf{b}_u \\ b_p \end{bmatrix}$

优点：稳定性好，收敛快
缺点：内存需求大

8. 数值仿真案例

8.1 二维方柱绕流

以下Python程序实现二维方柱绕流的简化CFD模拟：

"""
CFD数值模拟基础 - 二维方柱绕流仿真
使用涡量-流函数方法求解
"""

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from matplotlib.patches import Rectangle
import warnings
warnings.filterwarnings('ignore')

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

class CylinderFlowSolver:
    """二维方柱绕流求解器（涡量-流函数方法）"""
    
    def __init__(self, nx=200, ny=100, Re=100):
        """
        初始化求解器
        
        参数:
            nx, ny: 网格数
            Re: 雷诺数
        """
        self.nx = nx
        self.ny = ny
        self.Re = Re
        
        # 计算域尺寸
        self.Lx = 20.0  # 长度
        self.Ly = 10.0  # 高度
        
        # 网格间距
        self.dx = self.Lx / (nx - 1)
        self.dy = self.Ly / (ny - 1)
        
        # 方柱位置和尺寸
        self.cylinder_x = 5.0
        self.cylinder_y = self.Ly / 2
        self.cylinder_size = 1.0
        
        # 初始化场变量
        self.omega = np.zeros((nx, ny))  # 涡量
        self.psi = np.zeros((nx, ny))    # 流函数
        self.u = np.zeros((nx, ny))      # x方向速度
        self.v = np.zeros((nx, ny))      # y方向速度
        
        # 时间步长
        self.dt = min(self.dx, self.dy) / 4
        
    def is_inside_cylinder(self, i, j):
        """判断点是否在方柱内部"""
        x = i * self.dx
        y = j * self.dy
        return (abs(x - self.cylinder_x) <= self.cylinder_size/2 and 
                abs(y - self.cylinder_y) <= self.cylinder_size/2)
    
    def apply_boundary_conditions(self):
        """应用边界条件"""
        # 入口边界（左边界）
        self.psi[0, :] = self.psi[1, :]
        self.omega[0, :] = 0
        self.u[0, :] = 1.0  # 均匀来流
        self.v[0, :] = 0
        
        # 出口边界（右边界）
        self.psi[-1, :] = self.psi[-2, :]
        self.omega[-1, :] = self.omega[-2, :]
        self.u[-1, :] = self.u[-2, :]
        self.v[-1, :] = self.v[-2, :]
        
        # 上下边界（无滑移）
        self.psi[:, 0] = 0
        self.psi[:, -1] = self.Ly
        self.omega[:, 0] = -2 * self.psi[:, 1] / self.dy**2
        self.omega[:, -1] = -2 * (self.psi[:, -2] - self.Ly) / self.dy**2
        self.u[:, 0] = 0
        self.u[:, -1] = 0
        self.v[:, 0] = 0
        self.v[:, -1] = 0
        
        # 方柱边界（无滑移）
        for i in range(self.nx):
            for j in range(self.ny):
                if self.is_inside_cylinder(i, j):
                    self.psi[i, j] = 0
                    self.omega[i, j] = 0
                    self.u[i, j] = 0
                    self.v[i, j] = 0
    
    def solve_poisson_psi(self, max_iter=100, tol=1e-6):
        """求解流函数的泊松方程"""
        for _ in range(max_iter):
            psi_old = self.psi.copy()
            
            # 内部点
            for i in range(1, self.nx-1):
                for j in range(1, self.ny-1):
                    if not self.is_inside_cylinder(i, j):
                        self.psi[i, j] = (
                            (self.psi[i+1, j] + self.psi[i-1, j]) / self.dx**2 +
                            (self.psi[i, j+1] + self.psi[i, j-1]) / self.dy**2 +
                            self.omega[i, j]
                        ) / (2/self.dx**2 + 2/self.dy**2)
            
            # 检查收敛
            error = np.max(np.abs(self.psi - psi_old))
            if error < tol:
                break
    
    def update_velocity(self):
        """从流函数计算速度"""
        # 内部点
        for i in range(1, self.nx-1):
            for j in range(1, self.ny-1):
                if not self.is_inside_cylinder(i, j):
                    self.u[i, j] = (self.psi[i, j+1] - self.psi[i, j-1]) / (2*self.dy)
                    self.v[i, j] = -(self.psi[i+1, j] - self.psi[i-1, j]) / (2*self.dx)
    
    def update_vorticity(self):
        """更新涡量（对流-扩散方程）"""
        omega_new = self.omega.copy()
        
        for i in range(1, self.nx-1):
            for j in range(1, self.ny-1):
                if not self.is_inside_cylinder(i, j):
                    # 对流项（迎风格式）
                    if self.u[i, j] > 0:
                        domega_dx = (self.omega[i, j] - self.omega[i-1, j]) / self.dx
                    else:
                        domega_dx = (self.omega[i+1, j] - self.omega[i, j]) / self.dx
                    
                    if self.v[i, j] > 0:
                        domega_dy = (self.omega[i, j] - self.omega[i, j-1]) / self.dy
                    else:
                        domega_dy = (self.omega[i, j+1] - self.omega[i, j]) / self.dy
                    
                    conv = self.u[i, j] * domega_dx + self.v[i, j] * domega_dy
                    
                    # 扩散项（中心差分）
                    diff = (
                        (self.omega[i+1, j] - 2*self.omega[i, j] + self.omega[i-1, j]) / self.dx**2 +
                        (self.omega[i, j+1] - 2*self.omega[i, j] + self.omega[i, j-1]) / self.dy**2
                    ) / self.Re
                    
                    # 时间推进（欧拉显式）
                    omega_new[i, j] = self.omega[i, j] + self.dt * (-conv + diff)
        
        self.omega = omega_new
    
    def step(self):
        """单步推进"""
        # 更新涡量
        self.update_vorticity()
        
        # 求解流函数
        self.solve_poisson_psi()
        
        # 更新速度
        self.update_velocity()
        
        # 应用边界条件
        self.apply_boundary_conditions()
    
    def solve(self, n_steps=5000, plot_interval=1000):
        """
        求解流动
        
        参数:
            n_steps: 总步数
            plot_interval: 绘图间隔
        """
        print(f"开始求解: Re = {self.Re}, 网格 = {self.nx}×{self.ny}")
        
        # 初始条件
        self.apply_boundary_conditions()
        
        for step in range(n_steps):
            self.step()
            
            if (step + 1) % plot_interval == 0:
                print(f"  完成 {step+1}/{n_steps} 步")
        
        print("求解完成")
    
    def plot_results(self, filename='flow_results.png'):
        """绘制结果"""
        x = np.linspace(0, self.Lx, self.nx)
        y = np.linspace(0, self.Ly, self.ny)
        X, Y = np.meshgrid(x, y)
        
        fig, axes = plt.subplots(2, 2, figsize=(14, 10))
        
        # 图1：流线图
        ax1 = axes[0, 0]
        levels = np.linspace(-2, 2, 50)
        cs = ax1.contour(X, Y, self.psi.T, levels=levels, colors='black', linewidths=0.5)
        ax1.clabel(cs, inline=True, fontsize=8)
        ax1.set_xlabel('x', fontsize=11)
        ax1.set_ylabel('y', fontsize=11)
        ax1.set_title('Streamlines', fontsize=12, fontweight='bold')
        ax1.set_aspect('equal')
        
        # 绘制方柱
        rect = Rectangle(
            (self.cylinder_x - self.cylinder_size/2, self.cylinder_y - self.cylinder_size/2),
            self.cylinder_size, self.cylinder_size,
            facecolor='gray', edgecolor='black'
        )
        ax1.add_patch(rect)
        
        # 图2：涡量云图
        ax2 = axes[0, 1]
        omega_plot = np.clip(self.omega.T, -5, 5)
        im = ax2.contourf(X, Y, omega_plot, levels=50, cmap='RdBu_r')
        plt.colorbar(im, ax=ax2, label='Vorticity')
        ax2.set_xlabel('x', fontsize=11)
        ax2.set_ylabel('y', fontsize=11)
        ax2.set_title('Vorticity Contours', fontsize=12, fontweight='bold')
        ax2.set_aspect('equal')
        
        # 绘制方柱
        rect = Rectangle(
            (self.cylinder_x - self.cylinder_size/2, self.cylinder_y - self.cylinder_size/2),
            self.cylinder_size, self.cylinder_size,
            facecolor='gray', edgecolor='black'
        )
        ax2.add_patch(rect)
        
        # 图3：速度矢量图
        ax3 = axes[1, 0]
        skip = 5
        ax3.quiver(X[::skip, ::skip], Y[::skip, ::skip], 
                   self.u.T[::skip, ::skip], self.v.T[::skip, ::skip],
                   scale=20, width=0.003)
        ax3.set_xlabel('x', fontsize=11)
        ax3.set_ylabel('y', fontsize=11)
        ax3.set_title('Velocity Vectors', fontsize=12, fontweight='bold')
        ax3.set_aspect('equal')
        ax3.set_xlim(0, self.Lx)
        ax3.set_ylim(0, self.Ly)
        
        # 绘制方柱
        rect = Rectangle(
            (self.cylinder_x - self.cylinder_size/2, self.cylinder_y - self.cylinder_size/2),
            self.cylinder_size, self.cylinder_size,
            facecolor='gray', edgecolor='black'
        )
        ax3.add_patch(rect)
        
        # 图4：速度大小云图
        ax4 = axes[1, 1]
        velocity_magnitude = np.sqrt(self.u**2 + self.v**2)
        im = ax4.contourf(X, Y, velocity_magnitude.T, levels=50, cmap='viridis')
        plt.colorbar(im, ax=ax4, label='Velocity Magnitude')
        ax4.set_xlabel('x', fontsize=11)
        ax4.set_ylabel('y', fontsize=11)
        ax4.set_title('Velocity Magnitude', fontsize=12, fontweight='bold')
        ax4.set_aspect('equal')
        
        # 绘制方柱
        rect = Rectangle(
            (self.cylinder_x - self.cylinder_size/2, self.cylinder_y - self.cylinder_size/2),
            self.cylinder_size, self.cylinder_size,
            facecolor='gray', edgecolor='black'
        )
        ax4.add_patch(rect)
        
        plt.tight_layout()
        plt.savefig(filename, dpi=150, bbox_inches='tight')
        plt.close()
        print(f"✓ 结果图已保存: {filename}")


# ==================== 主程序 ====================

if __name__ == "__main__":
    print("="*70)
    print("CFD数值模拟基础 - 二维方柱绕流")
    print("="*70)
    
    # 创建求解器
    solver = CylinderFlowSolver(nx=150, ny=75, Re=100)
    
    # 求解
    solver.solve(n_steps=3000, plot_interval=500)
    
    # 绘制结果
    solver.plot_results('cylinder_flow_re100.png')
    
    print("\n" + "="*70)
    print("仿真完成！")
    print("="*70)

8.2 边界层流动

以下程序模拟平板边界层的发展：

"""
边界层流动模拟
使用Blasius方程求解层流边界层
"""

import numpy as np
import matplotlib
matplotlib.use('Agg')
import matplotlib.pyplot as plt
from scipy.integrate import odeint
import warnings
warnings.filterwarnings('ignore')

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


def blasius_equation(f, eta):
    """
    Blasius方程
    f''' + 0.5*f*f'' = 0
    
    转化为方程组:
    f0' = f1
    f1' = f2
    f2' = -0.5*f0*f2
    """
    f0, f1, f2 = f
    return [f1, f2, -0.5*f0*f2]


def solve_blasius():
    """求解Blasius方程"""
    # 边界条件
    # f(0) = 0, f'(0) = 0, f'(inf) = 1
    
    # 使用打靶法求解
    eta_max = 10
    eta = np.linspace(0, eta_max, 1000)
    
    # 初始猜测 f''(0)
    f2_0 = 0.332  # 已知精确值
    
    # 求解
    f0 = [0, 0, f2_0]
    sol = odeint(blasius_equation, f0, eta)
    
    return eta, sol


def plot_blasius_solution():
    """绘制Blasius解"""
    eta, sol = solve_blasius()
    
    f = sol[:, 0]  # 流函数
    df = sol[:, 1]  # 速度 u/U_inf
    d2f = sol[:, 2]  # 剪切应力
    
    fig, axes = plt.subplots(1, 3, figsize=(15, 4))
    
    # 图1：速度剖面
    axes[0].plot(df, eta, 'b-', linewidth=2)
    axes[0].set_xlabel('u/U∞', fontsize=11)
    axes[0].set_ylabel('η = y√(U∞/νx)', fontsize=11)
    axes[0].set_title('Velocity Profile', fontsize=12, fontweight='bold')
    axes[0].grid(True, alpha=0.3)
    axes[0].set_xlim(0, 1.2)
    axes[0].set_ylim(0, 10)
    
    # 图2：剪切应力分布
    axes[1].plot(d2f, eta, 'r-', linewidth=2)
    axes[1].set_xlabel("f''(η)", fontsize=11)
    axes[1].set_ylabel('η', fontsize=11)
    axes[1].set_title('Shear Stress Distribution', fontsize=12, fontweight='bold')
    axes[1].grid(True, alpha=0.3)
    axes[1].set_ylim(0, 10)
    
    # 图3：边界层厚度定义
    # 位移厚度
    delta_star = np.trapz(1 - df, eta)
    # 动量厚度
    theta = np.trapz(df * (1 - df), eta)
    
    axes[2].plot(eta, 1 - df, 'g-', linewidth=2, label='1 - u/U∞')
    axes[2].fill_between(eta, 0, 1 - df, alpha=0.3, color='green')
    axes[2].axvline(x=delta_star, color='r', linestyle='--', label=f'δ* = {delta_star:.3f}')
    axes[2].set_xlabel('η', fontsize=11)
    axes[2].set_ylabel('1 - u/U∞', fontsize=11)
    axes[2].set_title('Displacement Thickness', fontsize=12, fontweight='bold')
    axes[2].legend(loc='upper right')
    axes[2].grid(True, alpha=0.3)
    axes[2].set_xlim(0, 10)
    
    plt.tight_layout()
    plt.savefig('blasius_boundary_layer.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("Blasius解分析:")
    print(f"  壁面剪切应力 f''(0) = {d2f[0]:.6f}")
    print(f"  位移厚度 δ* = {delta_star:.6f}")
    print(f"  动量厚度 θ = {theta:.6f}")
    print(f"  形状因子 H = {delta_star/theta:.4f}")
    
    print("✓ Blasius边界层图已保存")


def simulate_boundary_layer_development():
    """模拟边界层沿平板的发展"""
    # 参数
    U_inf = 1.0  # 来流速度
    nu = 1e-5    # 运动黏度
    L = 1.0      # 平板长度
    
    # x位置
    x = np.linspace(0.01, L, 100)
    
    # 边界层厚度 (δ ≈ 5√(νx/U))
    delta = 5 * np.sqrt(nu * x / U_inf)
    
    # 位移厚度
    delta_star = 1.72 * np.sqrt(nu * x / U_inf)
    
    # 动量厚度
    theta = 0.664 * np.sqrt(nu * x / U_inf)
    
    # 壁面剪切应力
    tau_w = 0.332 * rho * U_inf**2 / np.sqrt(U_inf * x / nu)
    
    # 局部摩擦系数
    cf = 0.664 / np.sqrt(U_inf * x / nu)
    
    fig, axes = plt.subplots(2, 2, figsize=(12, 10))
    
    # 图1：边界层厚度发展
    axes[0, 0].plot(x, delta, 'b-', linewidth=2, label='δ (99% thickness)')
    axes[0, 0].plot(x, delta_star, 'r--', linewidth=2, label='δ* (displacement)')
    axes[0, 0].plot(x, theta, 'g-.', linewidth=2, label='θ (momentum)')
    axes[0, 0].set_xlabel('x (m)', fontsize=11)
    axes[0, 0].set_ylabel('Thickness (m)', fontsize=11)
    axes[0, 0].set_title('Boundary Layer Growth', fontsize=12, fontweight='bold')
    axes[0, 0].legend(loc='upper left')
    axes[0, 0].grid(True, alpha=0.3)
    
    # 图2：雷诺数
    Re_x = U_inf * x / nu
    axes[0, 1].semilogy(x, Re_x, 'b-', linewidth=2)
    axes[0, 1].set_xlabel('x (m)', fontsize=11)
    axes[0, 1].set_ylabel('Re_x', fontsize=11)
    axes[0, 1].set_title('Reynolds Number Distribution', fontsize=12, fontweight='bold')
    axes[0, 1].grid(True, alpha=0.3, which='both')
    
    # 图3：壁面剪切应力
    rho = 1.225  # 空气密度
    tau_w = 0.332 * rho * U_inf**2 / np.sqrt(Re_x)
    axes[1, 0].plot(x, tau_w, 'r-', linewidth=2)
    axes[1, 0].set_xlabel('x (m)', fontsize=11)
    axes[1, 0].set_ylabel('τ_w (Pa)', fontsize=11)
    axes[1, 0].set_title('Wall Shear Stress', fontsize=12, fontweight='bold')
    axes[1, 0].grid(True, alpha=0.3)
    
    # 图4：摩擦系数
    cf = 0.664 / np.sqrt(Re_x)
    axes[1, 1].semilogy(Re_x, cf, 'b-', linewidth=2)
    axes[1, 1].set_xlabel('Re_x', fontsize=11)
    axes[1, 1].set_ylabel('c_f', fontsize=11)
    axes[1, 1].set_title('Local Friction Coefficient', fontsize=12, fontweight='bold')
    axes[1, 1].grid(True, alpha=0.3, which='both')
    
    plt.tight_layout()
    plt.savefig('boundary_layer_development.png', dpi=150, bbox_inches='tight')
    plt.close()
    
    print("✓ 边界层发展图已保存")


# ==================== 主程序 ====================

if __name__ == "__main__":
    print("="*70)
    print("边界层流动模拟")
    print("="*70)
    
    print("\n[1/2] 求解Blasius方程...")
    plot_blasius_solution()
    
    print("\n[2/2] 模拟边界层发展...")
    simulate_boundary_layer_development()
    
    print("\n" + "="*70)
    print("仿真完成！")
    print("="*70)

8.3 结果分析

运行上述程序将生成以下结果：

方柱绕流（Re=100）：

流线图显示方柱后方的回流区
涡量云图显示剪切层和尾涡结构
速度矢量图显示流动分离和再附着
速度大小云图显示加速区和尾流区

边界层流动：

Blasius速度剖面符合理论解
边界层厚度沿平板逐渐增大
壁面剪切应力随距离减小
摩擦系数随雷诺数减小

AtomGit开源社区

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

更多推荐

中国开源软件的崛起与困境：贡献者生态的建立之难

AtomGit开源社区

技术创业的冷思考：从风口追赶到价值创造的转变

AtomGit开源社区

张量分解核心模型：从Tucker分解、CPD到耦合CPD与块项分解

本文系统介绍了四种张量分解核心模型：Tucker分解、CP分解、耦合CP分解和块项分解(BTD)。Tucker分解通过核张量与因子矩阵乘积表示高维数据；CP分解作为其特例，采用秩-1分量之和的紧凑形式；耦合CP分解通过共享因子矩阵实现多源数据融合；BTD则通过低多线性秩项组合处理复杂结构。文章详细阐述了各模型的数学表达、参数特性和适用场景，并对比了它们的优缺点：Tucker灵活性高但参数多，CP具

AtomGit开源社区

所有评论(0)

查看更多评论

kkchenkx

@2401_87715305

已为社区贡献156条内容

结构风振仿真-主题021-CFD数值模拟基础

kkchenkx

第21篇：CFD数值模拟基础

摘要

关键词

1. 引言

1.1 CFD在风工程中的重要性

1.2 CFD发展历程

1.3 CFD在风工程中的应用

2. 控制方程

2.1 连续性方程

2.2 Navier-Stokes方程

2.3 能量方程

2.4 方程的守恒形式

3. 湍流模型

3.1 湍流的基本特征

3.2 雷诺平均方法

3.3 k-ε湍流模型

3.4 k-ω SST模型

3.5 大涡模拟（LES）

4. 数值离散方法

4.1 有限差分法

4.2 有限体积法

4.3 时间离散方法

4.4 对流项离散

5. 边界条件

5.1 入口边界条件

5.2 出口边界条件

5.3 壁面边界条件

5.4 对称边界条件

5.5 周期性边界条件

6. 网格生成技术

6.1 网格类型

6.2 网格质量

6.3 自适应网格

7. 求解算法

7.1 SIMPLE算法

7.2 PISO算法

7.3 耦合求解

8. 数值仿真案例

8.1 二维方柱绕流

8.2 边界层流动

8.3 结果分析

所有评论(0)

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

kkchenkx