任意拉格朗日法是一种用于处理动边界问题的方法，尤其在流体问题中得到应用，以这个方法的学习作为这一系列笔记的结束。

参考文献如下

  • Donea, Jean, Antonio Huerta, and A Rodrıguez-Ferran. “Chapter 14 Arbitrary Lagrangian–Eulerian Methods,” n.d.
  • Fehn, Niklas, Johannes Heinz, Wolfgang A. Wall, and Martin Kronbichler. “High-Order Arbitrary Lagrangian-Eulerian Discontinuous Galerkin Methods for the Incompressible Navier-Stokes Equations.” Journal of Computational Physics 430 (April 2021): 110040. https://doi.org/10.1016/j.jcp.2020.110040.

1. 不可压 NS 方程

在区域 \(\Omega\) 上的不可压 NS 方程如下（欧拉描述）

\[ \begin{aligned} \frac{\partial u}{\partial t} + u \cdot \nabla u - \nu \nabla^2 u + \nabla p &= f \\ \nabla \cdot u &= 0 \end{aligned} \]

其中 \(u=u(x,t)\)，时间导数是固定空间坐标 \(x\) 的偏导， \(|_x\) 默认省略。

\[ \frac{\partial u}{\partial t} = \frac{\partial u}{\partial t}|_x \]

或者利用随体导数 \(|_X\)，物质质点 \(X\) 和空间坐标 \(x\) 之间存在对应关系 \(x= x(X,t)\)

\[ \begin{aligned} \frac{\partial u}{\partial t}|_X &= \frac{\partial u}{\partial t} + u \cdot \nabla u \end{aligned} \]

得到拉格朗日描述

\[ \begin{aligned} \frac{\partial u}{\partial t}|_X - \nu \nabla^2 u + \nabla p &= f \\ \nabla \cdot u &= 0 \end{aligned} \]

拉格朗日描述面临的问题是，流体的剧烈形变会破坏网格，计算难以继续。

2. ALE 描述的 NS 方程

假设底层网格也存在一个独立运动，引入网格坐标 \(\chi\)，现在有三个坐标体系：

  • 物质坐标\(X\)
  • 空间坐标 \(x\)
  • 网格坐标 \(\chi\)

存在三个坐标之间的对应转换，例如可以根据物质坐标和时间，算出当前物质坐标的空间位置，流体速度就是对应的偏导数

\[ \begin{aligned} x &= x(X,t)\\ u &= \frac{\partial \,x(X,t)}{\partial t}|_X \end{aligned} \]

还可以根据当前网格坐标和时间，算出当前网格坐标的空间位置，网格速度就是对应的偏导数

\[ \begin{aligned} x &= x(\chi,t)\\ w &= \frac{\partial \,x(\chi,t)}{\partial t}|_{\chi} \end{aligned} \]

注意到这里的流体速度 \(u\) 和网格速度 \(w\) 都与时间和空间坐标有关。

\[ \begin{aligned} u &= u(X,t) = u(X(x,t),t)\\ w &= w(\chi,t) = w(\chi(x,t),t)\\ \end{aligned} \]

对于关于空间坐标和时间的物理量 \(\phi=\phi(x,t)\)，可以视作关于物质坐标和时间的物理量

\[ \begin{aligned} \phi & = \phi(x(X,t),t)\\ \end{aligned} \]

或者视作关于网格坐标和时间的物理量

\[ \begin{aligned} \phi & = \phi(x(\chi,t),t)\\ \end{aligned} \]

固定不同坐标可以得到不同的时间偏导数，存在如下转换关系

\[ \begin{aligned} \frac{\partial \phi}{\partial t}|_X &= \frac{\partial \phi(x(X,t),t)}{\partial t}|_X = \frac{\partial \phi}{\partial t}|_x + (u \cdot \nabla) \phi \\ \frac{\partial \phi}{\partial t}|_\chi &= \frac{\partial \phi(x(\chi,t),t)}{\partial t}|_\chi = \frac{\partial \phi}{\partial t}|_x + (w \cdot \nabla) \phi \end{aligned} \]

因此有

\[ \frac{\partial\phi}{\partial t} = \frac{\partial \phi}{\partial t}|_x = \frac{\partial \phi}{\partial t}|_\chi - (w\cdot \nabla)\phi \]

取 \(\phi=u\) 代入 NS 方程即得到

\[ \begin{aligned} \frac{\partial u}{\partial t}|_\chi + (u-w) \cdot \nabla u - \nu \nabla^2 u + \nabla p &= f \\ \nabla \cdot u &= 0 \end{aligned} \]

这里的空间导数全都是关于空间坐标 \(x\) 的偏导，时间导数是 \(|_{\chi}\) 即固定网格坐标时的偏导。 ALE 引入的网格在空间中的运动是独立的，并且上式可以兼容欧拉描述和拉格朗日描述。

ALE 介于欧拉描述和拉格朗日描述之间：网格的运动既不是固定的，也不是完全跟随流体的。网格会跟随流体的运动进行一个平缓的形变，在动边界处尽可能与流体位移匹配，在固定边界处保持不动，在内部尽可能保证网格的质量。

3. ALE 推出欧拉描述

从下式出发

\[ \begin{aligned} \frac{\partial u}{\partial t}|_\chi + (u-w) \cdot \nabla u - \nu \nabla^2 u + \nabla p &= f \\ \nabla \cdot u &= 0 \end{aligned} \]

如果我们将网格固定，那么网格坐标 \(\chi\) 和空间坐标 \(x\) 的一一对应就与时间无关：\(x = x(\chi)\) ，从而网格速度恒为零

\[ \begin{aligned} w &= \frac{\partial x(\chi)}{\partial t} = 0 \\ \frac{\partial u}{\partial t}|_\chi &= \frac{\partial u}{\partial t} + (w\cdot \nabla)u = \frac{\partial u}{\partial t} \end{aligned} \]

因此得到

\[ \begin{aligned} \frac{\partial u}{\partial t} + u \cdot \nabla u - \nu \nabla^2 u + \nabla p &= f \\ \nabla \cdot u &= 0 \end{aligned} \]

即欧拉描述下的 NS 方程。

4. ALE 推出拉格朗日描述

从下式出发

\[ \begin{aligned} \frac{\partial u}{\partial t}|_\chi + (u-w) \cdot \nabla u - \nu \nabla^2 u + \nabla p &= f \\ \nabla \cdot u &= 0 \end{aligned} \]

如果我们令网格速度总是等于流体速度，也就是网格的形变时刻对应着流体的形变，网格坐标\(\chi\) 和物质坐标 \(X\) 的一一对应就与时间无关：\(X=X(\chi)\)，由于 \(u=w\)，我们可知

\[ \frac{\partial u}{\partial t}|_X = \frac{\partial u}{\partial t}|_x + (u \cdot \nabla) u = \frac{\partial u}{\partial t}|_x + (w \cdot \nabla) u = \frac{\partial u}{\partial t}|_\chi \]

因此得到

\[ \begin{aligned} \frac{\partial u}{\partial t}|_X - \nu \nabla^2 u + \nabla p &= f \\ \nabla \cdot u &= 0 \end{aligned} \]

即拉格朗日描述下的 NS 方程。

5 . NS 方程全离散

在 \(t^n \to t^{n+1}\) 的时间步，已有 \(t^n\) 时刻的三角网格 \(\Omega^n\)，其上的某个三角单元 \(K^n\)，以及定义在 \(\Omega^n\) 上的速度 \(u^n\)。

基于 \(u^n\) 我们可以构造相应的网格速度 \(w^n\) （在每一个顶点处需要一个移动速度），要求：

  • 在动边界的节点处，\(w^n\) 根据动边界的位移确定
  • 在固定边界的节点处，\(w^n=0\)
  • 在内部节点处，\(w^n\) 平缓过渡，保证三角网格的质量

然后，将 \(\Omega^n\) 的每一个顶点基于网格速度 \(w^n\) 进行移动，得到新的三角网格 \(\Omega^{n+1}\)，三角单元 \(K^n\) 变形为 \(K^{n+1}\)。

将定义在 \(\Omega^n\) 上的 \(u^n\) 转换为定义在 \(\Omega^{n+1}\) 上的 \(\tilde{u}^n\)：记参考单元 \(K\)，存在两个映射

\[ \begin{aligned} &\phi_{K^n}: K \to K^n\\ &\phi_{K^{n+1}}: K \to K^{n+1}\\ \end{aligned} \]

取

\[ \begin{aligned} \tilde{u}^{n}(\phi_{K^{n+1}}(x)) = u^n(\phi_{K^n}(x)) \\ \tilde{u}^n = u^n \circ \phi_{K^n} \circ \phi_{K^{n+1}}^{-1} \end{aligned} \]

这样可以体现出时间导的 \(|_\chi\) ，即固定网格坐标。以上的转换并不需要额外的计算。

从下式出发，向后欧拉离散。

\[ \begin{aligned} \frac{\partial u}{\partial t}|_\chi + (u-w) \cdot \nabla u - \nu \nabla^2 u + \nabla p &= f \\ \nabla \cdot u &= 0 \end{aligned} \]

任取定义在 \(K^{n+1}\) 上的测试函数 \(v,q\)，得到

\[ \begin{aligned} \frac{1}{\Delta t}((u^{n+1}-\tilde{u}^n),v)_{K^{n+1}} + ((u^{n+1}-w^n) \cdot \nabla u^{n+1},v)_{K^{n+1}} \\- \nu (\nabla^2 u^{n+1},v)_{K^{n+1}} + (\nabla p^{n+1},v)_{K^{n+1}} &= (f^{n+1},v)_{K^{n+1}} \\ (\nabla \cdot u^{n+1},q)_{K^{n+1}} &= 0 \end{aligned} \]

可以求解得到定义在三角网格 \(\Omega^{n+1}\) 的 \(u^{n+1},p^{n+1}\)。

Remark:

在欧拉描述下的 NS 方程中，时间导是 \(|_x\)，即固定空间坐标，计算网格基本保持不动。如果 \(t^n\) 时刻需要改变计算网格，我们就需要在不同三角形单元之间进行插值构造：根据 \(\Omega^n\) 上的 \(u^n\) 来得到 \(\Omega^{n+1}\) 的 \(\tilde{u}^n\) ，满足

\[ \tilde{u}^n(x) = u^n(x) \]

ALE 和欧拉描述下，定义在 \(K^{n+1}\) 的 \(\tilde{u}^n\) 是不一样的，差异体现在了网格移动多出来的对流项 \(w\cdot\nabla u\) 。

6 . 小结

从上述推导可以看出，对于一般的方程（欧拉描述，\(x\) 为空间坐标）

\[ \frac{\partial u}{\partial t}|_x = F(u) \]

转换为拉格朗日描述（\(X\) 为物质坐标）

\[ \frac{\partial u}{\partial t}|_X - \frac{\partial x}{\partial t}|_X\cdot \nabla u = F(u) \]

转换为 ALE 描述（\(\chi\) 为网格坐标）

\[ \frac{\partial u}{\partial t}|_\chi - \frac{\partial x}{\partial t}|_\chi\cdot \nabla u = F(u) \]