SAV求解两流问题

1.1. 笔记主页

笔记主页

1.2. 两相流

如果两流体密度$\rho_1=\rho_2=1$，则能量为
$$E(\phi) = \lambda \int_{\Omega} \frac{1}{2}|\nabla \phi|^2 + F(\phi) dx$$

• $\phi_t +(u\cdot \nabla) \phi = M\Delta \mu, \frac{\partial u}{\partial\bm n}|_{\Omega} = 0$，$M$是mobility

• $\mu = \lambda(-\Delta \phi + F'(\phi)), \frac{\partial \phi}{\partial\bm n}|_{\Omega} = 0$

• $\nabla \cdot u = 0$

• $u_t +u\cdot\nabla u = \nu\Delta u - \nabla p - \lambda\nabla \phi \otimes\nabla\phi, u|_{\partial \Omega}=0$

利用公式
$$\nabla\cdot(\nabla \phi \otimes\nabla\phi) = \Delta \phi \nabla \phi$$
再结合原方程可以得到
$$-(\nabla\cdot(\nabla \phi \otimes\nabla\phi),\mu) = (\mu\nabla \phi, \mu )$$
最后可以得到
$$\partial_t [E(\phi)+\frac{1}{2}||u||^2] + \int_{\Omega}(M|\nabla u|^2 + \nu ||\nabla u||^2)=0$$
对于二维可以证明出任意时刻的强解的存在，对于三维情形只能证明局部存在。

相场模型的好处是从变分原理导出来的

给出一阶格式如下

$$\begin{cases} \frac{\phi^{n+1} - \phi^n}{\Delta t} +(\tilde u^{n+1}\cdot \nabla) \phi^n = M\Delta \mu^{n+1} \\[1em] \mu^{n+1} = \lambda(-\Delta \phi^{n+1} + \xi^{n+1}F'(\phi^n)), \\[1em] r^{n+1}\frac{r^{n+1} - r^n}{\Delta t} = \xi^{n+1}(F'(\phi^n), \frac{\phi^{n+1} - \phi^n}{\Delta t} ) \\[1em] \xi^{n+1} = \frac{r^{n+1}}{\sqrt{\int_{\Omega} F(\phi^n) + C_0}} \end{cases}$$

$$\frac{\tilde u^{n+1} - \tilde u^n}{\Delta t} + u^n\cdot \nabla \tilde u^{n+1} = \nu \Delta \tilde u^{n+1} - \nabla p^{n} + \mu^{n+1}\nabla \phi^n$$

$$\begin{cases} \frac{u^{n+1} - \tilde u^{n+1}}{\Delta t} +\nabla(p^{n+1} - p^n)\\[1em] \nabla \cdot u^{n+1} = 0 \\[1em] u^{n+1}\cdot \bm n|_{\partial \Omega} = 0 \end{cases}$$

可以证明出修正的能量
$$\tilde E^{n} = \frac{1}{2}\left( ||\nabla \phi||^2 + ||u^n||^2 + \lambda (r^n)^2 \right) + \Delta t^2|| \nabla p^n ||^2$$
满足

$$\tilde E^{n+1} - \tilde E^n \leq \Delta t\left( -\int_{\Omega}M|\nabla \mu^{n+1}|^2 - \nu|| \nabla \tilde u^{n+1} ||^2 dx \right)$$

可以证明出左端系数矩阵是正定的，因此该系统存在唯一。

$$\mathbb A = \left[ \begin{matrix} \frac{1}{\Delta t}& - \nabla \cdot M\nabla &\cdot \nabla p^{n+1}\\[1em] \lambda\Delta &I&0\\[1em] 0&\nabla \phi^n&\frac{1}{\Delta t}I-\nu\Delta+u^n\cdot \nabla\\[1em] \end{matrix} \right]$$

取预条件子为
$$\mathbb A = \left[ \begin{matrix} \frac{1}{\Delta t}& - \nabla \cdot M\nabla &0\\[1em] \lambda \Delta&I&0\\[1em] 0&0&\frac{1}{\Delta t}I-\nu\Delta\\[1em] \end{matrix} \right]$$

可以利用2阶NS方程的方法构造上述二阶格式
$$\begin{cases} \dfrac{3\tilde u^{n+1} -4 u^{n}+u^{n-1}}{2\Delta t} + (2u^n-u^{n-1}) \cdot \nabla \tilde u^{n+1} \\[1em]\qquad= \nu\Delta \tilde u^{n+1} - \nabla (2p^{n}-p^{n-1}) \\[1em] \tilde u^{n+1} |_{\Omega} = 0 \end{cases}$$

$$\begin{cases} \dfrac{3(u^{n+1} - \tilde u^{n+1})}{2\Delta t} + \nabla (p^{n+1} - 2p^n + p^{n-1})= 0 \\[1em] \nabla \cdot u ^{n+1}= 0 \\[1em] u^{n+1} \cdot n|_{\Omega} = 0 \end{cases}$$
上面格式对速度为2阶但是压力是1阶，如果要把压力提高半阶，
可以引入一个变量$\psi$，并且这种做法可以在$H^1$模下
不改进
$$||u-u^n||_{L^2(0,T,L^2)} +\Delta t ||u-u^n||_{L^2(0,T,H^1)} + \Delta t|| p -p^n ||_{L^2(0,T,L^2)}$$

1.3. 修正格式

$$\begin{cases} \frac{\phi^{n+1} - \phi^n}{\Delta t} +(u_{\star}^{n}\cdot \nabla) \phi^n = M\Delta \mu^{n+1} \\[1em] \mu^{n+1} = \lambda(-\Delta \phi^{n+1} + \xi^{n+1}F'(\phi^n)), \\[1em] r^{n+1}\frac{r^{n+1} - r^n}{\Delta t} = \xi^{n+1}(F'(\phi^n), \frac{\phi^{n+1} - \phi^n}{\Delta t} ) \\[1em] \xi^{n+1} = \frac{r^{n+1}}{\sqrt{\int_{\Omega} F(\phi^n) + C_0}} \end{cases}$$

$$u_\star^{n} = u^n + \Delta t\mu^{n+1} \nabla \phi^n$$

$$\frac{\tilde u^{n+1} - \tilde u^n}{\Delta t} + u^n\cdot \nabla \tilde u^{n+1} = \nu \Delta \tilde u^{n+1} - \nabla p^{n} + \mu^{n+1}\nabla \phi^n$$

$$\begin{cases} \frac{u^{n+1} - \tilde u^{n+1}}{\Delta t} +\nabla(p^{n+1} - p^n)\\[1em] \nabla \cdot u^{n+1} = 0 \\[1em] u^{n+1}\cdot \bm n|_{\partial \Omega} = 0 \end{cases}$$

注意到
$$\mu\nabla \phi + \phi \nabla \mu = \nabla(\phi\mu)$$
因此$\mu\nabla \phi$与$\phi \nabla \mu$本质是一样的，文献中可能会不一样。

但是这种格式有变系数，并且无法构造二阶格式

修正：引入两个SAV的变量

1.4. possion方程求解

在可分区域(球，圆柱等)都有快速算法

有限差分

FFT

谱方法

罗威某教授写了个python包 shenfun(高度并行)
https://github.com/spectralDNS/shenfun

1.5. 变密度

如果$\rho_1 \neq \rho_2$，但是密度差不是很大的情况

可以做一个逼近，例如Boussinessq逼近

$$\rho = \frac{\rho_1 + \rho_2}{2}$$

动量方程
$$\rho_0 (u_t + u\cdot \nabla u) = \nabla\cdot \nu \nabla u - \nabla p + f(x)$$
以及
$$f(x) = - (1+\phi)g(\rho_1 - \rho_2) - (1-\phi)g(\rho_1- \rho_2)$$

$$\phi= \begin{cases} 1,\quad \text{流体1}\\ -1,\quad \text{流体2} \end{cases}$$

但是当两种密度相差很大的情况下，例如水和气想吃大约3个数量级，不能用，所以需要另外一种格式

质量守恒
$$\rho_t + \nabla \cdot (\rho u) = 0$$

不可压
$$\nabla \cdot u = 0$$

密度
$$\rho = \frac{\phi+1}{2} \rho_1 + \frac{1-\phi}{2}\rho_1$$
与黏性系数
$$\nu = \frac{\phi+1}{2} \nu_1 + \frac{1-\phi}{2}\nu_2$$
质量守恒与不可压方程不能同时满足，因此只能取其中一种格式，每种格式的模型是不一样的。

具体参考
http://math.purdure.edu/~shen/pub/SheY15.pdf

1.6. SAV应用

• 例如图像处理，自由能
  $$E(\phi) = \int_{\Omega} |\nabla \phi| + \frac{\lambda}{2}(\phi-g)^2 dx$$
  在能量里加上高阶稳定项目，高阶项前面乘以一个很小的参数；
  因为$|\nabla \phi|$求导之后容易出现分母为$0$，因此可以做修正
  $$|\nabla \phi| = \sqrt{ |\nabla \phi| ^2+\epsilon}$$

• 如果非线性部分出现在$\nabla$前面
  $$E(\phi) = \int_{\Omega} |\gamma(\phi)\nabla \phi|^2 + F(\phi) dx$$
  则可以加上稳定项
  这种处理对强非线性很有效

• 此外还可以与其它物理方程结合起来，整个系统一定要保证能量耗时，例如增加温度方程，磁流体麦克斯韦方程。