瞬态问题中的自动时步和阶数选择

当追踪瞬态仿真日志时，您可以观察到求解器使用了不同的时步和离散化阶数。在仿真过程中，您可能会发现时步变小，求解需要很长时间才能完成。本篇博文，我们将讨论选择时步和离散化阶数的原则。在接下来的一篇博文中，我们将说明在时步较小的情况下，应该采取哪些措施来提高仿真效率。

瞬态仿真

瞬态仿真需要计算反映时间演变的离散解。从初始值开始，通过时间积分方案确定未知的瞬态自由度。如果一个时步的计算解满足给定容差下的预定义误差边界，那么该解被接受。自适应时步算法可能会使用太大的时间步长，导致已采用的步长计算结果不会通过误差测试。在这些情况下，时步会减小，并且会重复计算该时步。如果非线性求解器不能在最大迭代次数内求解代数方程，时步也会减小。这两种方案都需要额外的计算，因此会导致瞬态仿真超过最优时间。

这篇博文旨在帮助您理解求解方案背后使用的机制，并帮助解读求解器日志中提供的信息。在这篇博文的后续文章中，我们将利用对这些机制的理解来举例说明，当遇到小时步时，我们如何使时间步进更加高效和稳定。下面，我们先来看看求解器日志，看看能从中了解到什么。

瞬态求解器日志

瞬态求解器的典型求解日志如下所示:

日志显示了当前 Time 的时间积分循环迭代计数器 Step，其中 Stepsize 是当前时步的大小。接下来的三列详细说明了残差(Res)的总数，雅可比组装(Jac)的总数，以及线性代数系统解( Sol)的总数。

在这篇博文中，我们感兴趣的是观察时间步进方案的离散阶数(Order)，自适应步长选择的失败次数(Tfail)，以及代数(非线性)求解器的总失败次数(NLfail)。我们将在接下来的章节中讨论这些失败的含义。根据线性代数系统误差估计(LinErr)和线性代数系统残差( LinRes)的大小，可以看到每个时间步的最后一个线性问题的解的信息。

如果找不到求解程序日志中列出的所有时步，您可以在瞬态求解器，高级 节点的常规 部分将 日志取样( 执行时间) 设置为 0。有时通过将求解器日志 设置为详细，将单个代数求解器迭代信息添加到求解器日志输出中非常有用。仿真过程中时步的变化可以在相应的收敛图 中找到，其中输出所有时步的步长倒数。收敛图中较大的值表示较小的时步。

如果您观察到求解器的时步非常小，这可能是由几个不同的原因造成的。一个原因可能是你的模型正在接近某种奇点(物理的或非物理的)，或者解正趋于无限。另一个原因可能是模型没有给出平滑的时间变化(例如，由于网格太粗糙)，这意味着对于任何时步，误差估计都很难实现。第三个原因可能是非线性难以处理(代数求解器不收敛)。为了详细了解时步的变化意味着什么，我们必须更深入地研究时步的确定，接下来，我们来看看 COMSOL 软件使用的一些时间步进方案。

离散时间步进方案

COMSOL 软件中的瞬态求解器 提供了三种不同的时间步进方法：隐式向后差分公式（BDF)，广义 α 方法，以及显式龙格-库塔方法。

BDF 求解器是一种隐式求解器，它使用具有可变离散化阶数和自动步长选择的向后差分公式。当解的质量允许时使用高阶方案，当需要额外的鲁棒性时使用低阶方案。BDF 方法的稳定性是众所周知的，默认情况下，它们用于涉及扩散、对流和反应的问题。

广义 α 方法 是一种二阶隐式方法，在结构力学问题中很常用，但它也很适合于波传播问题。广义 α 方法 允许我们控制高频阻尼，并且通常比二阶 BDF 方案的阻尼小。在许多情况下，它比 BDF 精确，但也不太稳定。

龙格-库塔 方法是通常用于常微分方程组（ODE）的显式方法。对于使用有限元法(FEM)求解的问题，它们通常效率不高。

在这篇博文中，我们关注的是 BDF 时间步进方案。这 BDF 方案可以在瞬态求解器 中的时间步进 部分选择。根据您的模型的物理特性，它可能已经被选中。

设置时步选择

在瞬态研究 步骤，你可以通过指定一个从开始时间到结束时间和中间步骤的明确时间 列表。使用内置算子可以快速指定中间时间点列表。一个常见的误解是，这些中间时间点是求解器采用的指定时间步长，但这些时间点是输出时间，解存储在这些时间点用于后处理和评估。默认情况下，求解器采用的时间步长由算法确定，该算法试图将误差(对于每个时间步长)保持在期望的限度内。以下设置可在瞬态求解器 中时间步进 部分的求解器采用的步骤中操作实现，以作用于实际的时步选择。

• 自由: 基于局部误差估计选择自适应时步。基于当前误差估计与容差的比较，减少或增加时步。如果所采取的时步仍然不满足误差估计，那么该步骤用减小的时步重新计算(Tfail 记录在求解器日志中)。由于失败的步骤是无效的，因此这种步长减小的代价很高。如果非线性求解器循环没有在最大迭代次数内收敛( NLfail 记录在求解器日志中)，指定的输出时间与求解器采用的时步无关。指定输出时间的解通过求解器所采取的时步之间的插值来计算。
• 精确：与自由 设置相同，但求解程序会调整其时步，也针对指定的输出时间进行求解。
• 中级：与自由 设置相同，但求解程序会调整其时步，在由输出时间定义的每个子区间中至少求解一次。
• 手动：采用用户指定的时步，而不是自动的。时步可以是常数、全局变量表达式，或者由时间值的单调列表之间的间隔确定。因此，该时步完全由用户控制，并且该时步的误差测试被禁用。但请注意，代数求解器的误差估计仍然有效，因此容差仍然是需要考虑的重要因素。如果在最大迭代次数后没有达到代数误差估计，则时步减小。

在下面的内容中，我们将重点关注自由 时间步进。中级 和精确 时间步进可用于将自适应时步选择的优势与某些重要建模时间或建模时步的手动强制相结合。精确 时间步进还避免了用户指定时间列表的插值，这对某些应用非常重要。但请注意，输出的插值也可以在自由方法 中避免，通过选择要存储的时步:求解器采用的步长。

时间离散化

在 COMSOL Multiphysics 中，一个瞬态偏微分方程组(PDEs)通过有限元空间离散化被转化为一个隐式的常微分方程(或微分代数方程；即 DAEs)。(注意，间断伽辽金法并非这种情况。)使用时间离散化(在最简单的情况下，即向后欧拉方案)为未知自由度生成一组代数方程。在每一个时间步中，代数方程都是通过类牛顿方法来求解（全耦合 求解器或含阻尼系数牛顿法分离式 求解器)并使用自动线性化。所得到的线性方程通过直接求解器(例如，MUMPS 或 PARDISO)，或者通过迭代求解器求解。

在 COMSOL Multiphysics 中，通过有限元离散化得到的隐式 ODE 系统根据用户提供的相对和绝对容差以预期的精度要求进行求解。有两种不同类型的要求，一种针对时间步进(求解器)误差，另一种针对代数方程(求解器)误差。第一种类型适用于自由，中级 和精确 方法，但不适用于手动 方法。第二个要求都适用。

人们通常会误解，认为当使用手动方法时，不使用容差。时间步进误差要求建立在对所用方法的局部截断误差的估计。代数求解器误差使用非常相似的要求，如在求解稳态问题时，即基于解向量增量的大小。但是该误差估计的归一化适合瞬态求解器使用的归一化。这样做是为了使代数误差范数与时间步进误差范数具有可比性，从而降低代数求解器误差的标量因子(容差系数)应该与时间步进误差相比均匀地抑制代数误差。有时减少这个因子很重要，以避免用代数误差影响计算的解。事实上，如果代数误差不小于时间步进误差，它会导致其他问题，如时间步进不稳定或误差测试失败。

上述两个要求是求解 ODE 系统的核心，下面进一步描述。但是，当这个系统源于有限元离散化时，需要记住的是其他的误差来源不考虑在内。例如，有限元法的截断误差、用数值方法近似有限元积分产生的求积误差以及实际几何的多项式表示产生的几何近似误差。

局部和全局截断误差以及离散化误差

我们来仔细看看截断误差。从有限元离散化，我们得到了隐式微分方程形式

其中，U=[U_1(t), …, U_N(t)] 为未知量的矢量或瞬态自由度， L 叫做残余矢量。

这里，我们假设只包含一阶时间导数，为了简单起见，我们忽略了约束处理。具有可变步长 \tau_k 的 q 阶 BDF 方法通过下式定义

其中，符号 U^k 表示解的近似值， U(t_k) 具有步长相关系数 \alpha_{k,i} 和 \beta_k。

这个表示显示了向后微分公式的 名称。这个公式可以用来消除时间导数，得到一个代数问题

在时间 t_k 的局部截断误差被定义为 e_{k} := U^k – U(t_k;U^{k-1})，其中 U(t_k;U^{k-1}) 是确切满足 U(t_{k-1}) = U^{k-1} 的解。通过这个解的泰勒展开，我们得到

如果也对残差进行泰勒展开，我们得到

其中每个定义的第一项都是相同的零。

在质量矩阵 D = -\frac{\partial{L}}{\partial \dot{U}} 可逆的情况下，对于 阶方案，有

对于不可逆的质量矩阵，我们有一个 DAE，误差分析更详细，但自适应时间步进的原理是相同的。这里，我们还假设不同的时间步长相互之间有些关联，因此时间步长不会随步长的不同而任意改变。全局截断误差 E^k := U^k-U(t_k) 才是最重要的。这个误差从所有局部误差中获得贡献，从概念上来说，就是将所有的局部误差加起来。对于 q 阶瞬态方法，并且具有恒定的时步 \tau_k = \tau，我们发现(通过适当的规范化) 在适当的假设下，

这里，常数 C 与问题密切相关。一个实际的常微分方程问题的局部误差既可以被方程的性质放大，也可以被方程的性质抑制，因此对局部误差简单求和的概念应该持保留态度。例如，对于误差被放大的情况，常数可能会非常大。此外，在仿真过程中，以恒定时间步长获得的局部误差可能有很大差异。这些观察结果表明，选择时间步长的一个好策略是在时间步进期间将局部误差保持在同一水平。还应该可以调整局部误差的水平，以便可以控制全局误差，并且在必要时可以减小全局误差。

对于给定的离散格式，误差估计的显式表达式可以按照下式给出

其中， C_k 是关于 q 和过去的步长 \tau_j,\, j\le k 的可计算函数。

例如，BDF 的显示模拟可以用来计算预测值

渐进分析给出

其中，常数 \bar{C}_k 是可计算的。

因此，我们得到了一个可计算的表达式，作为局部截断误差的上界

这个可计算的表达式在 COMSOL Multiphysics 中用作局部截断误差的近似值，该近似值又可用于控制或调整时间步长。

BDF 时间步进方案的时间步进控制

对于绝对容差 A，和相对容差 R，如果 |\bar{e}_{k}| \le A + R |U^{k}|，那么 BDF 时间步进方案中的时间步长可以被接受。在自由度为 N_j 的 M 阶情况下，如果误差估计的加权均方根范数

小于 1，则时间步长可以被接受。这里， \bar{e}_{k,i} 是在时间 t_k 时刻、缩放或未缩放的绝对容差 A_i 下，缩放或未缩放的场分量 U_i^k 的局部截断误差的估计值。你可以参阅 COMSOL Multiphysics Reference Manual 中的“隐式瞬态求解器算法”章节查看有关缩放和未缩放版本的详细信息。

我们发现每个时步的误差取决于两个因素:

1. 时步的大小，\tau_k
2. 离散化的阶数 q

由于高阶离散需要进行更多的工作，因此高阶格式以更高的计算成本给出更高的求解精度。时步越短，精度越高，成本也越高，因为总的来说需要计算更多的时步。自适应的 BDF 时间步进方法在保证所需求解精度的前提下，同时调整阶数和时步，尽量最小化成本。它对时步和阶数进行了复杂的预测，但仍然有可能当前时步的结果不符合误差限制。在这种情况下，时步将以较小的时步重复 \tau_k’，并且步进故障 Tfail 被记录下来。

控制代数误差的影响

请注意，对于非线性问题，需要在每个时间步长中求解一组非线性方程。这通常是通过类似牛顿迭代的方法来完成的，这种方法会产生额外的代数误差。此外，牛顿法中的线性化方程可以用线性迭代求解器来求解，而这些求解器只能在一定程度上近似线性化系统的精确解。

对于给定预定义的迭代次数，非线性求解器可能无法收敛。在这种情况下，如果需要的话，可以更新雅可比矩阵并且减小时间步长。在日志中，你会发现计数器 NLfail 增加了。通过检查时步 k 的第 m 次迭代的增量 \delta^{k,m} := U^{k,m}-U^{k,m-1} 来估计代数求解器误差，其中，U^{k,m} 是 U^k 的 m 阶牛顿迭代的近似值，(例如，参见日晷:非线性和微分/代数方程解算器套件)。

代数误差 |U^k-U^{k,m}| 可能会影响截断误差和稳定性。你需要确保代数误差足够小。容差因子 (在全耦合 和分离 求解器中的终止技术 部分设置)是一个安全系数，与截断误差相比，它可用于抑制代数误差。用于截断误差的加权范数也可同样用于代数误差减少容差因子 通常会导致非线性求解器中迭代次数增加(或相同)。

对于单个时步，用较小的容差因子求解方程的计算成本更高。但是，请注意，使用过大的容差因子 可能导致代数误差，干扰时间步进方法并导致不正确的时步选择、不准确的方法或不稳定的方法。最后两种影响也可能是手动的 时间步进法导致的。

时间步长选择和阶数选择

接下来，我们将描述 BDF 方案中时间步长的选择方法。对于瞬态求解器，绝对容差和相对容差控制每个积分步骤中的误差。相对公差 R 通常由物理场设置决定，但我们可以在瞬态研究 步骤中转换为手动控制。默认值 R=0.01。绝对容差 A 可以在瞬态求解器的绝对容差 部分设置。对于公差方法 中的因子选项，绝对容差是相对容差乘以指定系数。手动的 选项允许我们在默认值为 A=0.001 的情况下指定一个详细的值。通过使用变量列表，可以将绝对容差单独用于每个变量。

具体来说，BDF 计算 k^th 时步的局部截断误差的估计值，并且要求误差估计来满足不等式 \|\bar{e}_{k}\|_{\rm WRMS} < 1。在时间步进方案中，有一个初始阶段是经过特殊处理的。对于开始的几步，在每一步中步长 \tau_k 被加倍，阶数被提高，q=1,2,3, … 直到局部误差测试 \|\bar{e}_{k}\|_{\rm WRMS} < 1 失败；通过检查解的平滑性的机制来降低阶数 q；或者阶数 q 达到最大数量级，q=5。对于 q 附近的 k，BDF 方案中阶数的选择是基于比例导数范数 |\tau^k U^{(k)}| 随 k 单调递减的要求。这些范数再次使用以下事实进行评估

步长和阶数的选择从 T(q) 的单调性测试开始。 T(q) 递减，表示解是光滑的，阶数可以增加。当 T(q) 随着 q 增加是，阶数需要减小。否则，保持阶数 q。接下来，执行局部误差测试，如果失败，则按阶数 q’=q 和新的步长 \tau’ 重新执行该步骤。后者基于 \bar{e}_k 的 O(\tau^{q+1}) 的渐近行为。对于 \bar{e}_{k,\tau} > A + R |U|，我们要求对于校正的时步 \tau’，至少 \bar{e}_{k,\tau’} = A+R|U| 是成立的，这就产生了新的时步 \tau’。通过

我们衍生出

其中，在实施中使用了额外的安全系数。

如果误差小于允许值，则存在一个所谓的“死区”，在该区域中，时步不增加，直到误差小于允许值的 16 倍。然后，时间步长增加两倍。请注意，时步的确定主要由限制器调节。调节器仅在线性斜率范围内基于给定的公式调节。当估计的局部误差大于它们应有的值，并且需要减少时步时，该机制就起主导作用。在下图中(类似于参考文献1中的相应图)，可以确定斜率 -1/(q+1) 在线性斜率范围内，其中 \log(\tau_{k+1}/\tau_{k}) 被印在 \log(|\bar{e}_k|/(A+R|U^k|)) 。

如果通过局部误差测试，则调整下一步的步长和阶数。如果在上一步中进行了更改，则阶数不会改变。在阶数 q<5 和恒定步长时，如果最后的 q+1 步骤是按固定阶数进行的，BDF 考虑提高或降低阶数。如果 T(q) 随着 q 而减少，阶数 q 增加。如果
T(q) 随着 q 增加，阶数 q 减少。更多细节和启发可以在下列文档中找到：IDA solver in the Suite of Nonlinear and Differential/Algebraic Equation Solvers (SUNDIALS)。

关于自动时间步长和阶数选择的总结性思考

在这篇博文中，我们讨论了 BDF 时间步进方法中自动时间步长和阶数选择的原则。这些原则根据规定的容差，以较低的计算成本和最大的鲁棒性保证求解精度。在实践中，我们发现这种原则可以导致非常小的时间步骤和很长的模拟时间。时步过小以及依赖时间和非线性求解器步骤的重复失败表明可能需要进一步调整模型配置和求解器设置。我们还应该始终确保底层模型有一个合理的解。

在接下来的博客文章中，我们将介绍一些示例，并讨论在时步过小的情况下可以调整哪些求解器设置。敬请关注！

延伸阅读

通过浏览 COMSOL 知识库获得更多关于瞬态模型的详细信息:

• BDF 方法、α 法和龙格-库塔方法
• 求解初始值不一致的瞬态模型
• 控制瞬态求解器时步

参考文献

1. G. Söderlind, L. Wang, “Adaptive time-stepping and computational stability“, J. Comp. and Appl. Math., vol. 185, pp. 225–243, 2006.