20190809个人笔记(川大)

官网课程视频地址

笔记主页

两相流

如果两流体密度ρ1=ρ2=1\rho_1=\rho_2=1,则能量为
E(ϕ)=λΩ12ϕ2+F(ϕ)dx E(\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(ϕ)+12u2]+Ω(Mu2+νu2)=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+1rnΔ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+1u~nΔt+unu~n+1=νΔu~n+1pn+μ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+1u~n+1Δt+(pn+1pn)un+1=0un+1nΩ=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+un2+λ(rn)2)+Δt2pn2 \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+1E~nΔt(ΩMμn+12νu~n+12dx) \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ΔtMpn+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ΔtM0λΔ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+14un+un12Δt+(2unun1)u~n+1=νΔu~n+1(2pnpn1)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+1u~n+1)2Δt+(pn+12pn+pn1)=0un+1=0un+1nΩ=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模下
不改进
uunL2(0,T,L2)+ΔtuunL2(0,T,H1)+ΔtppnL2(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)}

修正格式

{ϕn+1ϕnΔt+(un)ϕn=MΔμn+1μn+1=λ(Δϕn+1+ξn+1F(ϕn)),rn+1rn+1rnΔ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}

un=un+Δtμn+1ϕn u_\star^{n} = u^n + \Delta t\mu^{n+1} \nabla \phi^n

u~n+1u~nΔt+unu~n+1=νΔu~n+1pn+μ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+1u~n+1Δt+(pn+1pn)un+1=0un+1nΩ=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的变量

possion方程求解

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

有限差分

FFT

谱方法

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

变密度

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

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

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

动量方程
ρ0(ut+uu)=νup+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,流体11,流体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

SAV应用