4. 流-固耦合（FSI）简化模型

任务目标

构建一个简化流-固耦合（FSI）PINN模型：流体域为不可压缩层流，结构域为线弹性薄板。在流固界面处，通过硬约束或损失加权实现两个物理场的耦合条件——速度连续（无滑移）与力平衡（流体压力等于结构表面应力）。训练后，模型应能同时输出界面位移与流体压力，且二者在界面处满足匹配关系。

问题定义与几何

二维通道（x∈[0,1], y∈[0,0.5]），底部为刚性壁面，顶部为弹性薄板（厚度忽略，视为y=0.5处的可变形边界）。流体从左侧以抛物线速度剖面流入，右侧自由流出。弹性板在流体压力作用下发生小变形，变形后板的位置为 y = 0.5 + w(x)，其中w为垂向位移。

流体域采用欧拉描述，结构域采用拉格朗日描述。为简化，将结构位移作为界面边界条件的一部分，通过迭代或联合训练实现耦合。本案例采用联合训练方式：单一神经网络同时输出流体速度(u,v)、压力p和结构位移w，通过损失函数将界面条件耦合。

控制方程与边界条件

  • 流体（稳态不可压缩Navier-Stokes，低雷诺数简化）：
    • 动量：ν∇²u - (u·∇)u - ∇p/ρ = 0（此处取ν=0.01, ρ=1）
    • 连续性：∇·u = 0
  • 结构（Euler-Bernoulli梁，小变形）：
    • EI·w'''' = p_interface(x)（压力载荷），取EI=0.01
    • 边界：w(0)=w(1)=0, w'(0)=w'(1)=0（固支）
  • 界面条件（y=0.5处）：
    • 运动学：v(x, y=0.5+w) = 0（无穿透，线性化后v(x,0.5)=0）
    • 动力学：p(x, y=0.5) = EI·w''''(x)（压力等于结构恢复力）

PINN架构与损失函数

网络输入(x,y)，输出(u,v,p,w)。由于w仅定义在界面上，实际实现中网络输出w作为x的函数（输入x，输出w），但为统一训练，我们让网络同时接收(x,y)，当y接近0.5时w有效。更简洁的做法：使用两个子网络——流体网络输出(u,v,p)，结构网络输出w(x)。耦合通过损失项实现。

关键完整代码（PyTorch实现）

import torch
import torch.nn as nn
import numpy as np

# ---------- 网络定义 ----------
class FsiPINN(nn.Module):
    def __init__(self):
        super().__init__()
        self.f_net = nn.Sequential(
            nn.Linear(2, 50), nn.Tanh(),
            nn.Linear(50, 50), nn.Tanh(),
            nn.Linear(50, 50), nn.Tanh(),
            nn.Linear(50, 3)   # (u, v, p)
        )
        self.s_net = nn.Sequential(
            nn.Linear(1, 30), nn.Tanh(),
            nn.Linear(30, 30), nn.Tanh(),
            nn.Linear(30, 1)   # w
        )
    
    def forward(self, x, y):
        uv_p = self.f_net(torch.cat([x, y], dim=1))
        u, v, p = uv_p[:, 0:1], uv_p[:, 1:2], uv_p[:, 2:3]
        w = self.s_net(x)
        return u, v, p, w

# ---------- 自动微分工具 ----------
def grad(outputs, inputs):
    return torch.autograd.grad(outputs, inputs, grad_outputs=torch.ones_like(outputs), create_graph=True)[0]

# ---------- 训练准备 ----------
model = FsiPINN()
optimizer = torch.optim.Adam(model.parameters(), lr=1e-3)

# 采样点
N_f = 5000   # 流体内部点
N_b = 200    # 边界点
N_i = 300    # 界面点

# 流体内部点 (0<x<1, 0<y<0.5)
x_f = torch.rand(N_f, 1) * 1.0
y_f = torch.rand(N_f, 1) * 0.5

# 界面点 (y=0.5)
x_i = torch.rand(N_i, 1) * 1.0
y_i = torch.full_like(x_i, 0.5)

# 入口边界 (x=0, 0<y<0.5) 抛物线速度
y_in = torch.rand(N_b, 1) * 0.5
x_in = torch.zeros_like(y_in)
u_in = 4 * y_in * (0.5 - y_in)   # 最大速度1在y=0.25

# 壁面 (y=0)
x_w = torch.rand(N_b, 1) * 1.0
y_w = torch.zeros_like(x_w)

# 出口 (x=1) 自由流出：p=0
y_out = torch.rand(N_b, 1) * 0.5
x_out = torch.ones_like(y_out)

# 结构边界 (x=0和x=1处w=0)
x_sb = torch.tensor([[0.0], [1.0]])

# ---------- 训练循环 ----------
for epoch in range(3000):
    optimizer.zero_grad()
    
    # 流体内部残差
    x_f.requires_grad_(True)
    y_f.requires_grad_(True)
    u, v, p, _ = model(x_f, y_f)
    
    u_x = grad(u, x_f); u_y = grad(u, y_f)
    v_x = grad(v, x_f); v_y = grad(v, y_f)
    p_x = grad(p, x_f); p_y = grad(p, y_f)
    u_xx = grad(u_x, x_f); u_yy = grad(u_y, y_f)
    v_xx = grad(v_x, x_f); v_yy = grad(v_y, y_f)
    
    nu = 0.01
    # 动量方程 (忽略对流项，Stokes流)
    f_mom_x = nu * (u_xx + u_yy) - p_x
    f_mom_y = nu * (v_xx + v_yy) - p_y
    f_cont = u_x + v_y
    
    loss_fluid = torch.mean(f_mom_x**2) + torch.mean(f_mom_y**2) + torch.mean(f_cont**2)
    
    # 入口边界
    u_in_pred, v_in_pred, _, _ = model(x_in, y_in)
    loss_in = torch.mean((u_in_pred - u_in)**2) + torch.mean(v_in_pred**2)
    
    # 壁面边界
    u_w, v_w, _, _ = model(x_w, y_w)
    loss_wall = torch.mean(u_w**2) + torch.mean(v_w**2)
    
    # 出口边界
    _, _, p_out, _ = model(x_out, y_out)
    loss_out = torch.mean(p_out**2)
    
    # 界面条件
    u_i, v_i, p_i, w_i = model(x_i, y_i)
    
    # 运动学：v=0（无穿透）
    loss_kin = torch.mean(v_i**2)
    
    # 动力学：p = EI * w''''（用有限差分近似四阶导数）
    # 使用界面点x_i计算w的四阶导数
    dx = 0.01
    x_plus2 = x_i + 2*dx; x_plus1 = x_i + dx
    x_minus1 = x_i - dx; x_minus2 = x_i - 2*dx
    
    w_plus2 = model.s_net(x_plus2)
    w_plus1 = model.s_net(x_plus1)
    w_minus1 = model.s_net(x_minus1)
    w_minus2 = model.s_net(x_minus2)
    
    # 四阶中心差分
    w_4th = (w_plus2 - 4*w_plus1 + 6*w_i - 4*w_minus1 + w_minus2) / (dx**4)
    EI = 0.01
    loss_dyn = torch.mean((p_i - EI * w_4th)**2)
    
    # 结构边界条件 w(0)=w(1)=0
    w_sb = model.s_net(x_sb)
    loss_sb = torch.mean(w_sb**2)
    
    # 总损失
    loss = loss_fluid + loss_in + loss_wall + loss_out + loss_kin + loss_dyn + loss_sb
    
    loss.backward()
    optimizer.step()
    
    if epoch % 500 == 0:
        print(f'Epoch {epoch}, Loss: {loss.item():.2e}')

# ---------- 验证 ----------
# 在界面处检查压力与结构位移关系
with torch.no_grad():
    x_test = torch.linspace(0.1, 0.9, 50).view(-1, 1)
    y_test = torch.full_like(x_test, 0.5)
    _, _, p_test, w_test = model(x_test, y_test)
    
    # 计算EI*w''''
    dx = 0.01
    w_4th_test = (model.s_net(x_test+2*dx) - 4*model.s_net(x_test+dx) + 
                  6*w_test - 4*model.s_net(x_test-dx) + model.s_net(x_test-2*dx)) / (dx**4)
    
    # 打印最大残差
    residual = torch.max(torch.abs(p_test - 0.01 * w_4th_test)).item()
    print(f'界面最大压力-结构残差: {residual:.4f}')
    print('若残差<0.1，说明界面匹配良好')

运行方式

  1. 安装依赖：pip install torch numpy
  2. 直接运行上述脚本，训练约3000轮（CPU约2-3分钟）

验证方法与通过条件

  • 损失收敛：总损失降至1e-3以下，且各分量（尤其loss_dyn和loss_kin）不再显著下降
  • 界面匹配：训练结束后，程序输出界面最大压力-结构残差，若该值<0.1则通过
  • 物理合理性：入口处速度剖面近似抛物线，壁面处速度接近0，界面处v≈0

常见问题

  • 训练不收敛：尝试降低学习率至1e-4，或增加内部点数量N_f
  • 四阶导数不稳定：dx取值不宜过小（0.01较合适），若出现NaN可增大dx至0.02
  • 耦合强度不足：可对loss_dyn和loss_kin增加权重系数（如×10）以强化界面约束
入口 出口 刚性壁面 y=0 弹性板 y=0.5+w(x) p(x) w(x) 流-固耦合界面示意

该简化模型忽略了流体对流项（Stokes流）和结构剪切变形，适用于低雷诺数微流控或薄板小变形场景。实际工程中若需考虑惯性效应，可在动量方程中加入对流项(u·∇)u，代码中取消注释并增加相应自动微分即可。

← 上一章 📖 返回目录 下一章 →