SAV求解两流问题

Author : zbzhen,        Modified : Tue Feb 21 15:55:26 2023

1.1. 笔记主页

笔记主页

1.2. 两相流

如果两流体密度ρ1=ρ2=1\rho_1=\rho_2=1，则能量为
E(ϕ)=λ∫Ω12∣∇ϕ∣2+F(ϕ)dxE(\phi) = \lambda \int_{\Omega} \frac{1}{2}|\nabla \phi|^2 + F(\phi) dx

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

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

给出一阶格式如下

{ϕn+1−ϕnΔt+(u~n+1⋅∇)ϕn=MΔμn+1μn+1=λ(−Δϕn+1+ξn+1F′(ϕn)),rn+1rn+1−rnΔt=ξn+1(F′(ϕn),ϕn+1−ϕnΔt)ξn+1=rn+1∫ΩF(ϕn)+C0\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}

u~n+1−u~nΔt+un⋅∇u~n+1=νΔu~n+1−∇pn+μn+1∇ϕ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

{un+1−u~n+1Δt+∇(pn+1−pn)∇⋅un+1=0un+1⋅n∣∂Ω=0\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}

可以证明出修正的能量
E~n=12(∣∣∇ϕ∣∣2+∣∣un∣∣2+λ(rn)2)+Δt2∣∣∇pn∣∣2\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
满足

E~n+1−E~n≤Δt(−∫ΩM∣∇μn+1∣2−ν∣∣∇u~n+1∣∣2dx)\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)

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

A=[1Δt−∇⋅M∇⋅∇pn+1λΔI00∇ϕn1ΔtI−νΔ+un⋅∇]\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]

取预条件子为
A=[1Δt−∇⋅M∇0λΔI0001ΔtI−νΔ]\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方程的方法构造上述二阶格式
{3u~n+1−4un+un−12Δt+(2un−un−1)⋅∇u~n+1=νΔu~n+1−∇(2pn−pn−1)u~n+1∣Ω=0\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}

{3(un+1−u~n+1)2Δt+∇(pn+1−2pn+pn−1)=0∇⋅un+1=0un+1⋅n∣Ω=0\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，并且这种做法可以在H1H^1模下
不改进
∣∣u−un∣∣L2(0,T,L2)+Δt∣∣u−un∣∣L2(0,T,H1)+Δt∣∣p−pn∣∣L2(0,T,L2)||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. 修正格式

{ϕn+1−ϕnΔt+(u⋆n⋅∇)ϕn=MΔμn+1μn+1=λ(−Δϕn+1+ξn+1F′(ϕn)),rn+1rn+1−rnΔt=ξn+1(F′(ϕn),ϕn+1−ϕnΔt)ξn+1=rn+1∫ΩF(ϕn)+C0\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⋆n=un+Δtμn+1∇ϕnu_\star^{n} = u^n + \Delta t\mu^{n+1} \nabla \phi^n

u~n+1−u~nΔt+un⋅∇u~n+1=νΔu~n+1−∇pn+μn+1∇ϕ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

{un+1−u~n+1Δt+∇(pn+1−pn)∇⋅un+1=0un+1⋅n∣∂Ω=0\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. 变密度

如果ρ1≠ρ2\rho_1 \neq \rho_2，但是密度差不是很大的情况

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

ρ=ρ1+ρ22\rho = \frac{\rho_1 + \rho_2}{2}

动量方程
ρ0(ut+u⋅∇u)=∇⋅ν∇u−∇p+f(x)\rho_0 (u_t + u\cdot \nabla u) = \nabla\cdot \nu \nabla u - \nabla p + f(x)
以及
f(x)=−(1+ϕ)g(ρ1−ρ2)−(1−ϕ)g(ρ1−ρ2)f(x) = - (1+\phi)g(\rho_1 - \rho_2) - (1-\phi)g(\rho_1- \rho_2)

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

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

质量守恒
ρt+∇⋅(ρu)=0\rho_t + \nabla \cdot (\rho u) = 0

不可压
∇⋅u=0\nabla \cdot u = 0

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

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

1.6. SAV应用