MFEM项目中时间相关非齐次Dirichlet边界条件的波动方程实现

波动方程作为描述波动现象的经典偏微分方程，在科学计算和工程应用中具有重要地位。本文将详细介绍如何在MFEM框架下实现带有时间相关非齐次Dirichlet边界条件的波动方程求解。

问题描述

我们考虑二维区域[0,1]×[0,1]上的波动方程，边界条件设置为：

• 在x=0且1/6<y<5/6的边界上，施加随时间变化的正弦函数边界条件
• 其他边界保持零值

边界条件具体表达式为： u(x,t) = sin(4πt) ∂u/∂t(x,t) = 4πcos(4πt) ∂²u/∂t²(x,t) = -16π²sin(4πt)

数值方法实现

1. 空间离散

使用MFEM的有限元空间进行空间离散：

• 采用H1有限元空间
• 使用双线性形式构建质量矩阵M和刚度矩阵K
• 处理Dirichlet边界条件需要构建边界自由度列表

2. 时间离散

采用Newmark方法进行时间离散，波动方程离散形式为： [M + β(Δt)²K]∂²u(t+Δt)/∂t² = -K[u(t) + Δt∂u(t)/∂t + (1/2-β)(Δt)²∂²u(t)/∂t²]

其中β是Newmark方法的参数，控制数值阻尼特性。

3. 边界条件处理

关键难点在于处理时间相关的非齐次Dirichlet边界条件。实现要点包括：

1. 边界条件函数定义：创建三个函数分别表示u、∂u/∂t和∂²u/∂t²在边界上的值

2. 隐式求解处理：在ImplicitSolve方法中，需要消除边界条件对右端项的贡献，包括：

   • 当前时间步的边界值贡献
   • 时间导数项的边界值贡献
   • 二阶时间导数项的边界值贡献

3. 边界值设置：在每个时间步更新边界上的解值

实现细节

边界条件投影

使用ProjectBdrCoefficient方法将边界条件函数投影到边界自由度上：

u_bdr->SetTime(t);
u_gf.ProjectBdrCoefficient(*u_bdr, ess_bdr);

右端项修正

在隐式求解中，需要修正右端项以消除边界条件影响：

// 消除u_B的贡献
K->EliminateVDofsInRHS(ess_tdof_list, u_gf, z);

// 消除Δt*∂u/∂t_B的贡献  
K->SpMatElim().AddMult(dudt_gf, z, -current_dt);

// 消除(0.5-β)(Δt)²∂²u/∂t²_B的贡献
double fac_d2udt2 = (0.5 - beta) * (current_dt * current_dt);
K->SpMatElim().AddMult(d2udt2_gf, z, -fac_d2udt2);

时间步进控制

主循环中控制时间步进，并在每个时间步更新边界条件：

for (int ti = 1; !last_step; ti++) {
    // 更新边界条件
    u_bdr.SetTime(t);
    u_gf.ProjectBdrCoefficient(u_bdr, ess_bdr);
    
    // 时间步进
    ode_solver->Step(u, dudt, t, dt);
    
    // 后处理和可视化
    if (last_step || (ti % vis_steps) == 0) {
        // 更新可视化数据
    }
}

结果分析

通过7次均匀网格加密后，可以得到波动在区域内的传播过程。由于边界条件的正弦激励，会在x=0边界产生波动，并向区域内部传播。数值结果展示了波动在区域内的反射和干涉现象。

实现注意事项

1. 时间步长选择：需要满足CFL条件以保证数值稳定性

2. 边界条件同步：必须确保u、∂u/∂t和∂²u/∂t²的边界条件在时间上同步更新

3. 矩阵组装优化：对于固定时间步长问题，可以避免重复组装系统矩阵

4. 并行处理：当前实现主要针对串行计算，并行实现需要额外考虑边界数据的通信

本文介绍的方法可以推广到其他具有时间相关边界条件的二阶波动问题，为MFEM用户处理类似问题提供了参考实现。