MFEM中求解非线性瞬态热传导问题的隐式方法

概述

在MFEM框架中求解非线性瞬态热传导问题时，开发者常常会遇到如何处理非线性项和时域离散的问题。本文将详细介绍如何使用MFEM的NonlinearForm类来实现完全隐式的非线性瞬态问题求解。

问题描述

考虑以下非线性热传导方程：

∂u/∂t = ∇·((κ+αu)∇u)

这是一个典型的非线性瞬态问题，其中热导率κ与温度场u相关。与线性问题不同，这里的扩散系数依赖于解本身，使得问题具有非线性特性。

传统线性化方法

在MFEM的示例ex16中，采用了显式处理非线性项的方法：

1. 在每个时间步使用上一步的温度场u计算扩散系数
2. 通过BilinearForm构建线性系统
3. 求解该线性系统

这种方法虽然实现简单，但属于半隐式方法，时间步长可能受到稳定性限制。

完全隐式方法

要实现完全隐式求解，我们需要：

1. 使用时域离散方法(如后向欧拉法)离散时间导数项
2. 将非线性扩散项在当前时间步进行隐式处理
3. 使用Newton迭代法求解非线性系统

实现关键点

1. 非线性形式定义

使用NonlinearForm类来定义非线性弱形式。对于后向欧拉离散，残差可以表示为：

F(u) = (u-uₙ)/Δt - ∇·((κ+αu)∇u)

对应的变分形式为：

<F(u),v> = ∫(u-uₙ)v/Δt dx + ∫(κ+αu)∇u·∇v dx

2. Jacobian矩阵计算

Newton方法需要计算Jacobian矩阵J=∂F/∂u。对于上述问题：

<J(u)δu,v> = ∫δuv/Δt dx + ∫[αδu∇u·∇v + (κ+αu)∇δu·∇v] dx

3. 时间步进循环

在每个时间步内：

1. 初始化Newton迭代使用上一步解
2. 构建非线性残差和Jacobian
3. 求解线性系统获得增量
4. 更新解并检查收敛性

示例参考

MFEM的示例10(p/ex10.cpp)展示了如何使用NonlinearForm求解非线性稳态问题。对于瞬态问题，可以在此基础上：

1. 添加时间导数项到残差计算中
2. 在Jacobian计算中包括1/Δt项
3. 实现时间步进循环

实现建议

1. 派生自定义NonlinearForm类，重写GetGradient和Mult方法
2. 在Mult中计算包含时间项的残差
3. 在GetGradient中计算包含时间导数和非线性项的Jacobian
4. 使用SUNDIALS包中的ARKStep或CVODE求解器(如示例16)处理时间积分

结论

通过使用MFEM的NonlinearForm类配合Newton迭代法，可以实现非线性瞬态问题的完全隐式求解。这种方法虽然计算量较大，但具有更好的数值稳定性和精度，特别适用于强非线性问题。开发者可以根据具体问题特点选择合适的线性化策略和求解方法。