用Python玩转AI流体力学：从零实现PINN求解圆柱绕流问题（附完整代码）-育师

用Python玩转AI流体力学：从零实现PINN求解圆柱绕流问题（附完整代码）

流体力学，这门研究“流动”的科学，长久以来都是工程与科研领域的硬骨头。从飞机机翼的设计到心脏血液的模拟，其背后的纳维-斯托克斯方程（N-S方程）以其复杂的非线性特性，让无数工程师和研究者耗费海量计算资源。传统的计算流体力学方法，如有限体积法或有限元法，虽然成熟，但往往需要精细的网格划分、复杂的迭代求解，以及对计算硬件的高要求。有没有一种方法，能让我们用一种更“智能”、更“直接”的方式，去逼近这些复杂流动的解呢？

近年来，一个名为物理信息神经网络的概念横空出世，它巧妙地将深度学习的强大拟合能力与物理定律的约束结合起来。想象一下，你训练一个神经网络，它的损失函数不仅包含与实验数据的误差，更关键的是包含了物理控制方程本身——这意味着网络在学习数据的同时，也必须“学会”遵守物理规律。对于像我们这样有一定Python基础，喜欢动手实践的工程师或研究生来说，这无疑打开了一扇新的大门：我们不再仅仅是CFD软件的使用者，而是可以亲手构建一个“懂得物理”的AI求解器。

本文将带你从零开始，用Python和主流的深度学习框架，一步步构建一个PINN模型，去求解流体力学中那个经典的基准问题：二维圆柱绕流。我们会绕过繁复的理论推导，直接切入代码实操，在Jupyter Notebook的环境中，完成从问题定义、数据准备、模型构建、训练技巧到结果可视化的全流程。你会发现，整个过程就像在搭一个有趣的乐高模型，每一块代码都有其明确的物理意义。文末将提供完整的、可复现的代码，让你能立刻在自己的机器上运行起来，亲眼见证AI是如何“计算”出流体绕过圆柱时那优美的卡门涡街的。

1. 问题定义与环境搭建

在动手写代码之前，我们必须清晰地定义我们要解决什么问题，以及需要准备什么样的工具。圆柱绕流是一个被广泛研究的经典问题，它描述了均匀来流绕过静止圆柱体时产生的流动现象。在低雷诺数下，流动是对称的；随着雷诺数升高，流动分离，并在圆柱后方形成周期性脱落的涡旋，即卡门涡街。

我们的目标，是使用PINN来求解一个特定雷诺数下的、不可压缩流体的稳态或非稳态圆柱绕流流场。PINN的核心思想是：将流场的解（速度u, v和压力p）用一个深度神经网络来参数化，网络的输入是空间坐标(x, y)和时间t，输出就是对应的(u, v, p)。然后，我们要求这个网络满足三大类约束：

物理方程约束：网络输出的(u, v, p)必须近似满足纳维-斯托克斯方程和连续性方程。
边界条件约束：在计算域的边界上（如入口、出口、圆柱壁面、上下边界），网络输出必须满足给定的边界条件（如无滑移、压力出口等）。
初始条件约束（对于非稳态问题）：在初始时刻，网络输出必须满足给定的初始流场。

我们的损失函数就是这些约束的加权和。通过最小化这个损失函数，我们就在“训练”网络成为一个符合物理定律的流场求解器。

1.1 搭建Python环境

为了高效实现，我们需要一个配置好的Python环境。推荐使用conda创建一个独立的环境，避免包冲突。

# 创建并激活一个名为pinn_flow的新环境 conda create -n pinn_flow python=3.9 conda activate pinn_flow # 安装核心科学计算和深度学习库 pip install numpy scipy matplotlib jupyter pip install torch torchvision torchaudio --index-url https://download.pytorch.org/whl/cu118 # 根据你的CUDA版本选择 pip install tensorboard # 可选，用于训练可视化 # 安装一些好用的工具库 pip install tqdm # 进度条 pip install seaborn # 更美观的绘图

这里我们选择PyTorch作为深度学习框架，因为它对动态计算图和自定义损失函数的支持非常灵活，非常适合PINN这种需要手动定义复杂损失项的场景。当然，使用TensorFlow的Keras API也可以实现，但PyTorch的代码通常更直观。

注意：确保你的PyTorch安装支持GPU（如果可用），这将极大加速训练过程。可以使用torch.cuda.is_available()来检查。

1.2 定义计算域与物理参数

接下来，我们在代码中定义问题的几何和物理参数。假设我们的计算域是一个矩形：x ∈ [-2, 8],y ∈ [-2, 2]，圆柱圆心位于(0, 0)，半径为0.5。来流速度U_inf = 1.0，流体密度ρ = 1.0，我们通过设定雷诺数Re = 100来动态计算动力粘度ν。

import torch import numpy as np # 物理参数 Re = 100.0 # 雷诺数 U_inf = 1.0 # 来流速度 (x方向) rho = 1.0 # 密度 nu = U_inf * 0.5 * 2 / Re # 动力粘度，特征长度取圆柱直径D=1 print(f"动力粘度 nu = {nu:.4f}") # 计算域 x_min, x_max = -2.0, 8.0 y_min, y_max = -2.0, 2.0 cylinder_center = (0.0, 0.0) cylinder_radius = 0.5 # 用于后续采样的边界和内部区域定义 # 我们将分别对内部点、边界点、圆柱壁面点、初始条件点进行采样

2. 构建物理信息神经网络模型

PINN模型的核心是一个全连接神经网络（也称为多层感知机MLP）。这个网络结构并不复杂，关键在于我们如何利用它的输出并施加物理约束。

2.1 设计网络架构

我们构建一个将坐标(x, y, t)映射到流场变量(u, v, p)的网络。网络通常有多个隐藏层，使用激活函数（如tanh或sin）引入非线性。研究表明，周期性的激活函数（如sin）有时对求解偏微分方程有奇效。

import torch.nn as nn class PINN(nn.Module): """ 物理信息神经网络模型。 输入: (x, y, t) 坐标 输出: (u, v, p) 速度分量和压力 """ def __init__(self, layers=[3, 50, 50, 50, 50, 3], activation=nn.Tanh()): super(PINN, self).__init__() self.activation = activation # 动态构建全连接层 self.linears = nn.ModuleList() for i in range(len(layers)-1): self.linears.append(nn.Linear(layers[i], layers[i+1])) # 一种常用的技巧：初始化权重，有助于训练稳定性 self._initialize_weights() def _initialize_weights(self): for m in self.linears: if isinstance(m, nn.Linear): nn.init.xavier_normal_(m.weight) nn.init.constant_(m.bias, 0.0) def forward(self, xyt): """ 前向传播。 Args: xyt: 输入张量，形状为 [batch_size, 3]，最后一维是 (x, y, t) Returns: uv: 输出张量，形状为 [batch_size, 3]，最后一维是 (u, v, p) """ a = xyt for i, linear in enumerate(self.linears[:-1]): a = self.activation(linear(a)) # 最后一层不使用激活函数，用于直接输出物理量 uv = self.linears[-1](a) return uv

这个网络结构[3, 50, 50, 50, 50, 3]意味着有4个隐藏层，每层50个神经元。你可以根据问题复杂度调整层数和神经元数量。

2.2 实现物理损失函数

这是PINN最精华的部分。我们需要计算N-S方程和连续性方程的残差。对于不可压缩流，方程如下：

连续性方程：∇·u = 0，即 ∂u/∂x + ∂v/∂y = 0
动量方程（N-S方程）： ∂u/∂t + u·∇u = - (1/ρ) ∇p + ν ∇²u （v分量类似）

我们需要利用PyTorch的自动微分功能来计算网络输出关于输入(x, y, t)的偏导数。

def compute_physics_loss(model, points, nu, rho): """ 计算物理方程残差损失。 Args: model: PINN模型实例 points: 内部采样点，形状 [N, 3]，包含 (x, y, t) nu: 动力粘度 rho: 密度 Returns: loss_f: 动量方程残差的均方误差 loss_c: 连续性方程残差的均方误差 """ # 要求计算梯度 points.requires_grad_(True) # 网络前向传播，得到预测的u, v, p uv = model(points) u = uv[:, 0:1] # 保持维度，便于后续求导 v = uv[:, 1:2] p = uv[:, 2:3] # 计算一阶偏导：u_t, u_x, u_y, v_t, v_x, v_y, p_x, p_y u_grad = torch.autograd.grad(u.sum(), points, create_graph=True)[0] v_grad = torch.autograd.grad(v.sum(), points, create_graph=True)[0] p_grad = torch.autograd.grad(p.sum(), points, create_graph=True)[0] u_t, u_x, u_y = u_grad[:, 2:3], u_grad[:, 0:1], u_grad[:, 1:2] v_t, v_x, v_y = v_grad[:, 2:3], v_grad[:, 0:1], v_grad[:, 1:2] p_x, p_y = p_grad[:, 0:1], p_grad[:, 1:2] # 计算二阶偏导：u_xx, u_yy, v_xx, v_yy u_x_grad = torch.autograd.grad(u_x.sum(), points, create_graph=True)[0] u_xx, u_xy = u_x_grad[:, 0:1], u_x_grad[:, 1:2] u_y_grad = torch.autograd.grad(u_y.sum(), points, create_graph=True)[0] u_yy = u_y_grad[:, 1:2] v_x_grad = torch.autograd.grad(v_x.sum(), points, create_graph=True)[0] v_xx, v_xy = v_x_grad[:, 0:1], v_x_grad[:, 1:2] v_y_grad = torch.autograd.grad(v_y.sum(), points, create_graph=True)[0] v_yy = v_y_grad[:, 1:2] # 连续性方程残差：f_c = u_x + v_y f_continuity = u_x + v_y # x方向动量方程残差：f_mx = u_t + u*u_x + v*u_y + (1/rho)*p_x - nu*(u_xx + u_yy) f_momentum_x = u_t + u*u_x + v*u_y + (1/rho)*p_x - nu*(u_xx + u_yy) # y方向动量方程残差：f_my = v_t + u*v_x + v*v_y + (1/rho)*p_y - nu*(v_xx + v_yy) f_momentum_y = v_t + u*v_x + v*v_y + (1/rho)*p_y - nu*(v_xx + v_yy) # 损失为残差的均方误差 loss_f = torch.mean(f_momentum_x**2 + f_momentum_y**2) loss_c = torch.mean(f_continuity**2) return loss_f, loss_c

这段代码是PINN的“心脏”。它利用torch.autograd.grad进行高阶自动微分，精确地计算了控制方程的残差。网络训练的目标就是让这些残差尽可能接近零。

3. 数据采样与边界条件实现

PINN是“无网格”方法，但我们仍然需要在计算域内和边界上采样一系列“点”来评估损失。采样策略直接影响训练的效率和效果。

3.1 智能采样策略

我们不需要均匀的网格，而是可以更有针对性地采样。例如，在流动梯度大的区域（如圆柱壁面附近、尾流区）可以采样更密集的点。

def sample_points(batch_size, mode='interior', t=0.0): """ 采样不同区域的点。 Args: batch_size: 每批采样点数 mode: 采样模式，'interior', 'inlet', 'outlet', 'walls', 'cylinder', 'initial' t: 时间（对于非稳态问题） Returns: points: 采样点张量，形状 [batch_size, 3] """ points = torch.zeros(batch_size, 3) points[:, 2] = t # 设置时间 if mode == 'interior': # 内部点，排除圆柱内部区域 num_samples = 0 while num_samples < batch_size: # 在矩形域内均匀采样候选点 candidates = torch.rand(batch_size * 2, 2) candidates[:, 0] = candidates[:, 0] * (x_max - x_min) + x_min candidates[:, 1] = candidates[:, 1] * (y_max - y_min) + y_min # 判断是否在圆柱外 dist_from_center = torch.sqrt((candidates[:, 0] - cylinder_center[0])**2 + (candidates[:, 1] - cylinder_center[1])**2) mask = dist_from_center > cylinder_radius valid_candidates = candidates[mask][:batch_size - num_samples] points[num_samples:num_samples+len(valid_candidates), 0:2] = valid_candidates num_samples += len(valid_candidates) elif mode == 'inlet': # 入口边界 (x = x_min) points[:, 0] = x_min points[:, 1] = torch.rand(batch_size) * (y_max - y_min) + y_min elif mode == 'outlet': # 出口边界 (x = x_max) points[:, 0] = x_max points[:, 1] = torch.rand(batch_size) * (y_max - y_min) + y_min elif mode == 'walls': # 上下边界 (y = y_min 或 y = y_max) half = batch_size // 2 points[:half, 0] = torch.rand(half) * (x_max - x_min) + x_min points[:half, 1] = y_min points[half:, 0] = torch.rand(batch_size - half) * (x_max - x_min) + x_min points[half:, 1] = y_max elif mode == 'cylinder': # 圆柱壁面点，在圆柱表面均匀采样角度 theta = torch.rand(batch_size) * 2 * np.pi points[:, 0] = cylinder_center[0] + cylinder_radius * torch.cos(theta) points[:, 1] = cylinder_center[1] + cylinder_radius * torch.sin(theta) elif mode == 'initial': # 初始条件点（整个域，t=0） points = sample_points(batch_size, 'interior', t=0.0) return points

3.2 实现边界条件损失

边界条件通过额外的损失项来施加。例如，在入口处，我们要求速度等于来流速度；在圆柱壁面，要求无滑移条件（速度为零）；在上下边界，可以是自由滑移或对称边界。

def compute_boundary_loss(model, points, mode): """ 计算边界条件损失。 """ uv_pred = model(points) u_pred, v_pred, p_pred = uv_pred[:, 0:1], uv_pred[:, 1:2], uv_pred[:, 2:3] loss = 0.0 if mode == 'inlet': # 入口: u = U_inf, v = 0 u_target = torch.ones_like(u_pred) * U_inf v_target = torch.zeros_like(v_pred) loss = torch.mean((u_pred - u_target)**2 + (v_pred - v_target)**2) elif mode == 'outlet': # 出口: 压力p = 0 (或 Neumann 条件)，这里简单用 Dirichlet 条件 p=0 p_target = torch.zeros_like(p_pred) loss = torch.mean((p_pred - p_target)**2) elif mode == 'walls': # 上下壁面: 自由滑移条件，法向速度为零？或者无滑移？这里假设为自由滑移，v=0 # 更常见的是对称边界，但为简化，我们设 v=0 v_target = torch.zeros_like(v_pred) loss = torch.mean((v_pred - v_target)**2) elif mode == 'cylinder': # 圆柱壁面: 无滑移条件，u=0, v=0 u_target = torch.zeros_like(u_pred) v_target = torch.zeros_like(v_pred) loss = torch.mean((u_pred - u_target)**2 + (v_pred - v_target)**2) elif mode == 'initial': # 初始条件: t=0时，全流场 u=U_inf, v=0 (或更复杂的初始场) u_target = torch.ones_like(u_pred) * U_inf v_target = torch.zeros_like(v_pred) loss = torch.mean((u_pred - u_target)**2 + (v_pred - v_target)**2) return loss

通过组合物理损失和边界损失，我们得到了PINN的总损失函数。各个损失项之间的权重平衡是一个需要调整的超参数，通常需要根据具体问题反复试验。

4. 模型训练、技巧与可视化

有了模型和损失函数，我们就可以开始训练了。训练一个PINN模型比训练一个普通的分类网络更具挑战性，需要一些特别的技巧。

4.1 训练循环与自适应加权

一个基础的训练循环如下所示。我们使用Adam优化器，并可能采用学习率衰减。

def train_pinn(model, epochs, batch_sizes, lr=1e-3, device='cuda'): """ 训练PINN模型。 Args: model: PINN模型 epochs: 训练轮数 batch_sizes: 字典，包含各区域每批采样点数，如 {'interior': 1000, 'inlet': 200, ...} lr: 初始学习率 device: 训练设备 """ model.to(device) optimizer = torch.optim.Adam(model.parameters(), lr=lr) scheduler = torch.optim.lr_scheduler.StepLR(optimizer, step_size=5000, gamma=0.9) # 损失权重（超参数，需要调优） lambda_phys = 1.0 lambda_bc = 1.0 lambda_ic = 1.0 history = {'total_loss': [], 'phys_loss': [], 'bc_loss': []} for epoch in range(epochs): model.train() optimizer.zero_grad() total_loss = 0.0 # 1. 物理方程损失 interior_points = sample_points(batch_sizes['interior'], mode='interior', t=0.0).to(device) loss_f, loss_c = compute_physics_loss(model, interior_points, nu, rho) phys_loss = lambda_phys * (loss_f + loss_c) total_loss += phys_loss # 2. 边界条件损失 bc_loss_sum = 0.0 for bc_type in ['inlet', 'outlet', 'walls', 'cylinder']: points = sample_points(batch_sizes.get(bc_type, 100), mode=bc_type, t=0.0).to(device) loss_bc = compute_boundary_loss(model, points, bc_type) bc_loss_sum += loss_bc bc_loss = lambda_bc * bc_loss_sum total_loss += bc_loss # 3. 初始条件损失（对于非稳态问题，需要采样t=0的点） # initial_points = sample_points(batch_sizes['initial'], mode='initial', t=0.0).to(device) # loss_ic = lambda_ic * compute_boundary_loss(model, initial_points, 'initial') # total_loss += loss_ic # 反向传播与优化 total_loss.backward() optimizer.step() scheduler.step() # 记录损失 history['total_loss'].append(total_loss.item()) history['phys_loss'].append(phys_loss.item()) history['bc_loss'].append(bc_loss.item()) if epoch % 1000 == 0: print(f'Epoch {epoch:05d} | Total Loss: {total_loss.item():.4e} | Phys Loss: {phys_loss.item():.4e} | BC Loss: {bc_loss.item():.4e}') return history

训练PINN的一个常见难题是损失项尺度不平衡。物理残差、边界条件损失的数值量级可能相差很大，导致优化过程被大数值的损失项主导。解决方案之一是采用自适应加权，即在训练过程中动态调整各损失项的权重。

# 一种简单的自适应加权策略示例 class AdaptiveWeightedLoss: def __init__(self, num_terms, initial_weights=None, alpha=0.9): self.num_terms = num_terms self.weights = torch.ones(num_terms) if initial_weights is None else torch.tensor(initial_weights) self.loss_history = [] self.alpha = alpha # 平滑因子 def update_weights(self, current_losses): # current_losses 是一个列表，包含当前各个损失项的值 if not self.loss_history: self.loss_history.append(current_losses) else: # 指数移动平均 smoothed = [] for i in range(self.num_terms): new_val = self.alpha * self.loss_history[-1][i] + (1-self.alpha) * current_losses[i] smoothed.append(new_val) self.loss_history.append(smoothed) # 根据历史损失调整权重（例如，让损失大的项权重变小） # 这里是一个简化策略：权重与平均损失的倒数成正比 avg_losses = np.mean(self.loss_history[-10:], axis=0) if len(self.loss_history) >= 10 else self.loss_history[-1] self.weights = torch.tensor([1.0 / (l + 1e-8) for l in avg_losses]) # 归一化权重，使其和为固定值（如损失项数量） self.weights = self.weights / torch.sum(self.weights) * self.num_terms

4.2 结果可视化与验证

训练完成后，我们需要评估模型的效果。最直观的方式就是将网络预测的流场可视化出来，并与经典CFD结果或实验数据进行对比。

import matplotlib.pyplot as plt def visualize_results(model, device='cuda'): """ 可视化预测的流场。 """ model.eval() # 创建用于评估的网格点 nx, ny = 200, 100 x = np.linspace(x_min, x_max, nx) y = np.linspace(y_min, y_max, ny) X, Y = np.meshgrid(x, y) # 将网格点转换为模型输入格式 points = np.hstack([X.reshape(-1,1), Y.reshape(-1,1), np.zeros((nx*ny, 1))]) # t=0 points_tensor = torch.tensor(points, dtype=torch.float32).to(device) with torch.no_grad(): uv_pred = model(points_tensor).cpu().numpy() U_pred = uv_pred[:, 0].reshape(ny, nx) V_pred = uv_pred[:, 1].reshape(ny, nx) P_pred = uv_pred[:, 2].reshape(ny, nx) # 创建一个遮挡数组，将圆柱内部的点设为NaN，以便在图中隐藏 dist_from_center = np.sqrt((X - cylinder_center[0])**2 + (Y - cylinder_center[1])**2) mask = dist_from_center < cylinder_radius U_pred[mask] = np.nan V_pred[mask] = np.nan P_pred[mask] = np.nan # 绘制速度大小云图 speed = np.sqrt(U_pred**2 + V_pred**2) fig, axes = plt.subplots(1, 3, figsize=(18, 5)) im0 = axes[0].contourf(X, Y, speed, levels=50, cmap='jet') axes[0].set_title('Predicted Velocity Magnitude') axes[0].set_xlabel('x') axes[0].set_ylabel('y') axes[0].add_patch(plt.Circle(cylinder_center, cylinder_radius, color='k', fill=True)) plt.colorbar(im0, ax=axes[0]) # 绘制压力云图 im1 = axes[1].contourf(X, Y, P_pred, levels=50, cmap='RdBu_r') axes[1].set_title('Predicted Pressure Field') axes[1].set_xlabel('x') axes[1].set_ylabel('y') axes[1].add_patch(plt.Circle(cylinder_center, cylinder_radius, color='k', fill=True)) plt.colorbar(im1, ax=axes[1]) # 绘制流线图 axes[2].streamplot(X, Y, U_pred, V_pred, color=speed, linewidth=1, cmap='jet', density=2) axes[2].set_title('Streamlines') axes[2].set_xlabel('x') axes[2].set_ylabel('y') axes[2].add_patch(plt.Circle(cylinder_center, cylinder_radius, color='k', fill=True)) axes[2].set_xlim([x_min, x_max]) axes[2].set_ylim([y_min, y_max]) plt.tight_layout() plt.show() # 计算并打印一些定量指标，例如圆柱表面的阻力系数（需要积分） # 这里仅作示意 print("可视化完成。可以观察尾流区是否对称，以及圆柱前后的压力分布是否合理。")

为了验证结果的正确性，我们可以将PINN预测的圆柱表面压力系数分布与经典文献（如Taira和Colonius的“浸入边界法”结果）或高精度CFD仿真结果进行对比。压力系数Cp的定义为：Cp = (p - p_inf) / (0.5 * ρ * U_inf^2)，其中p_inf是参考压力（通常取出口压力或无穷远处压力）。

def plot_pressure_coefficient(model, device='cuda'): """ 绘制圆柱表面的压力系数分布。 """ model.eval() # 在圆柱表面均匀采样角度 theta = np.linspace(0, 2*np.pi, 200) x_surface = cylinder_center[0] + cylinder_radius * np.cos(theta) y_surface = cylinder_center[1] + cylinder_radius * np.sin(theta) t_surface = np.zeros_like(theta) points_surface = np.vstack([x_surface, y_surface, t_surface]).T points_tensor = torch.tensor(points_surface, dtype=torch.float32).to(device) with torch.no_grad(): uv_pred = model(points_tensor).cpu().numpy() p_surface = uv_pred[:, 2] # 假设出口压力为参考压力 p_inf = 0 Cp = (p_surface - 0) / (0.5 * rho * U_inf**2) plt.figure(figsize=(8, 5)) plt.plot(theta * 180 / np.pi, Cp, 'b-', linewidth=2, label='PINN Prediction') # 这里可以添加经典CFD或实验数据作为对比，例如： # theta_ref, Cp_ref = load_reference_data() # plt.plot(theta_ref, Cp_ref, 'ro', markersize=4, label='Reference (CFD)') plt.xlabel('Angle (degree)') plt.ylabel('Pressure Coefficient Cp') plt.title('Pressure Coefficient on Cylinder Surface') plt.grid(True, linestyle='--', alpha=0.7) plt.legend() plt.tight_layout() plt.show()

4.3 完整代码整合与运行

将以上所有模块整合到一个脚本或Jupyter Notebook中，你就可以运行完整的PINN求解流程了。下面是一个简化的主程序结构：

def main(): # 0. 设置设备 device = torch.device('cuda' if torch.cuda.is_available() else 'cpu') print(f"Using device: {device}") # 1. 初始化模型 model = PINN(layers=[3, 50, 50, 50, 50, 3], activation=nn.Tanh()).to(device) # 2. 定义批次大小 batch_sizes = { 'interior': 2000, 'inlet': 200, 'outlet': 200, 'walls': 200, 'cylinder': 300, 'initial': 1000 } # 3. 训练模型 print("开始训练PINN模型...") history = train_pinn(model, epochs=20000, batch_sizes=batch_sizes, lr=1e-3, device=device) # 4. 可视化训练损失曲线 plt.figure(figsize=(10, 4)) plt.plot(history['total_loss'], label='Total Loss') plt.plot(history['phys_loss'], label='Physics Loss') plt.plot(history['bc_loss'], label='BC Loss') plt.yscale('log') plt.xlabel('Epoch') plt.ylabel('Loss (log scale)') plt.title('Training Loss History') plt.legend() plt.grid(True) plt.show() # 5. 可视化预测结果 visualize_results(model, device) # 6. 绘制压力系数 plot_pressure_coefficient(model, device) print("PINN圆柱绕流求解完成！") if __name__ == '__main__': main()

运行这段代码，你将看到训练损失逐渐下降，最终得到预测的流场。对于雷诺数100的稳态流动，你应该能看到一个对称的、稳定的尾流区。如果想模拟非稳态的卡门涡街，你需要引入时间t作为输入，并在损失函数中加入初始条件项，同时在时间维度上进行采样。这会使问题变得更复杂，但基本框架是一致的。

5. 进阶技巧与挑战

当你成功运行了基础版本的PINN后，可能会发现结果并不完美：收敛慢、精度不够、或者对于高雷诺数问题直接失败。这引出了PINN在实际应用中的几个核心挑战和进阶技巧。

5.1 处理高雷诺数与复杂流动

高雷诺数流动具有更强的非线性和更薄的边界层，这对PINN的逼近能力提出了更高要求。以下是一些改进策略：

网络架构优化：尝试更深的网络、更多的神经元，或者使用残差连接（ResNet块）来缓解梯度消失问题。也可以尝试傅里叶特征编码，将输入坐标(x, y)映射到高频空间，有助于网络学习高频信号。

class FourierFeatureEncoding(nn.Module): """傅里叶特征编码，将低维输入映射到高维空间。""" def __init__(self, input_dim, mapping_size=256, scale=10.0): super().__init__() self.B = torch.randn((input_dim, mapping_size)) * scale def forward(self, x): x_proj = 2 * torch.pi * x @ self.B return torch.cat([torch.sin(x_proj), torch.cos(x_proj)], dim=-1)

将编码后的特征再输入到MLP中。

损失函数改进：除了MSE损失，可以尝试加权残差，在物理梯度大的区域（如边界层、剪切层）赋予更高的权重。这需要先进行一轮粗略训练来估计残差分布，或者采用自适应采样策略，在训练过程中动态地在残差大的区域增加采样点。
多尺度训练与课程学习：不要一开始就用高雷诺数训练。可以从低雷诺数（如Re=10）开始训练一个模型，然后将这个模型的权重作为高雷诺数问题的预训练权重，再进行微调。这类似于“课程学习”，让网络先学会简单的流动，再逐步挑战复杂的。

5.2 加速训练与提升精度

PINN训练慢是出了名的。除了使用GPU，还有以下方法可以加速：

混合方法：纯粹的PINN（无数据）有时很难收敛。可以引入少量高精度CFD数据或实验数据作为监督信号。这形成了“数据驱动”与“物理驱动”的混合损失，能极大地引导训练方向，加快收敛。这就是所谓的数据增强的PINN。
领域分解：对于大尺度或复杂几何问题，可以将计算域划分为多个子域，为每个子域训练一个PINN，并在子域交界处施加连续性条件。这能降低单个网络的建模难度，也便于并行训练。

使用更高效的优化器：Adam是默认选择，但对于一些病态问题，L-BFGS二阶优化器可能表现出更好的收敛性，尤其是在训练后期。可以尝试先用Adam训练几千轮，再用L-BFGS进行微调。

# 在训练后期切换优化器示例 if epoch > 10000: optimizer = torch.optim.LBFGS(model.parameters(), lr=1, max_iter=20, history_size=50) def closure(): optimizer.zero_grad() # ... 重新计算损失 ... loss.backward() return loss optimizer.step(closure)

5.3 从稳态到非稳态：捕捉卡门涡街

要模拟非稳态的圆柱绕流（卡门涡街），关键变化在于：

输入维度：网络输入变为(x, y, t)。
采样策略：需要在时空域(x, y, t)内采样。时间t的范围从0到某个总时长T。
损失函数：必须包含初始条件损失（t=0时刻整个流场的条件），并且在计算物理损失时，需要包含时间导数项u_t,v_t。
输出分析：训练完成后，固定空间坐标(x, y)，让时间t变化，网络就能输出该点的速度随时间变化的曲线，从而得到涡脱落的频率（斯特劳哈尔数）。

这无疑增加了问题的维度和难度，但框架是通用的。一个实用的技巧是，可以先训练一个稳态求解器（忽略时间项），将其结果作为非稳态问题的初始猜测，这能显著改善非稳态训练的收敛性。

6. 工程实践中的考量与扩展

将PINN从玩具问题推向实际工程应用，还需要考虑更多现实因素。

几何复杂性：本文的圆柱是简单的隐式几何（通过距离函数判断）。对于任意复杂几何，需要一种通用的方式将几何信息融入采样和损失计算。一种方法是使用符号距离函数或水平集函数来描述几何，并在SDF=0的等值面上施加边界条件。另一种更工程化的方法是与CAD软件或网格生成工具结合，直接读取表面网格点坐标作为边界采样点。

多物理场耦合：真实的流体问题往往耦合了热传导、物质输运、化学反应等多物理场。PINN的优势在于可以很自然地将多个控制方程同时纳入损失函数。例如，在能量方程和N-S方程之间，只需要在网络的输出中增加温度场T，并在损失函数中加入能量方程的残差即可。这种“端到端”的耦合求解避免了传统方法中复杂的迭代耦合流程。

逆问题与参数辨识：PINN特别擅长解决逆问题。假设我们有一些流场的部分观测数据（如PIV实验测得的某些点的速度），但不知道流体的粘度或边界条件。我们可以将未知参数（如粘度ν）也作为可训练参数，与网络权重一起优化。损失函数由“物理方程残差”和“观测数据误差”共同构成。通过训练，网络不仅能重构出完整流场，还能反推出未知的物理参数。这在实验流体力学中具有巨大潜力。

高性能计算：虽然PINN声称“无网格”，但其训练过程本身计算量巨大，尤其是需要大量采样点和自动微分。利用多GPU数据并行来分发不同的采样批次，或者探索模型并行来分解大型网络，是走向大规模应用的必经之路。此外，针对物理方程的特定结构，可以设计更高效的自动微分策略，避免计算海森矩阵等不必要的开销。

最后，我想分享一点个人在调试PINN时的体会：耐心和细致的调参是关键。损失权重的平衡、学习率策略、网络深度和宽度、激活函数的选择、采样点的分布，每一个环节都可能决定训练的成败。它不像训练一个图像分类网络那样有标准配方，更像是在进行一场科学实验，需要你根据物理直觉和训练反馈不断调整。当看到自己编写的网络最终成功复现出优美的物理流场时，那种成就感是无可比拟的。希望本文提供的代码和思路，能成为你探索AI+流体力学这个迷人交叉领域的一块坚实跳板。