SAV求解两流问题
笔记主页
如果两流体密度ρ1=ρ2=1,则能量为
E(ϕ)=λ∫Ω21∣∇ϕ∣2+F(ϕ)dx
-
ϕt+(u⋅∇)ϕ=MΔμ,∂n∂u∣Ω=0,M是mobility
-
μ=λ(−Δϕ+F′(ϕ)),∂n∂ϕ∣Ω=0
-
∇⋅u=0
-
ut+u⋅∇u=νΔu−∇p−λ∇ϕ⊗∇ϕ,u∣∂Ω=0
利用公式
∇⋅(∇ϕ⊗∇ϕ)=Δϕ∇ϕ
再结合原方程可以得到
−(∇⋅(∇ϕ⊗∇ϕ),μ)=(μ∇ϕ,μ)
最后可以得到
∂t[E(ϕ)+21∣∣u∣∣2]+∫Ω(M∣∇u∣2+ν∣∣∇u∣∣2)=0
对于二维可以证明出任意时刻的强解的存在,对于三维情形只能证明局部存在。
相场模型的好处是从变分原理导出来的
给出一阶格式如下
⎩⎨⎧Δtϕn+1−ϕn+(u~n+1⋅∇)ϕn=MΔμn+1μn+1=λ(−Δϕn+1+ξn+1F′(ϕn)),rn+1Δtrn+1−rn=ξn+1(F′(ϕn),Δtϕn+1−ϕn)ξn+1=∫ΩF(ϕn)+C0rn+1
Δtu~n+1−u~n+un⋅∇u~n+1=νΔu~n+1−∇pn+μn+1∇ϕn
⎩⎨⎧Δtun+1−u~n+1+∇(pn+1−pn)∇⋅un+1=0un+1⋅n∣∂Ω=0
可以证明出修正的能量
E~n=21(∣∣∇ϕ∣∣2+∣∣un∣∣2+λ(rn)2)+Δt2∣∣∇pn∣∣2
满足
E~n+1−E~n≤Δt(−∫ΩM∣∇μn+1∣2−ν∣∣∇u~n+1∣∣2dx)
可以证明出左端系数矩阵是正定的,因此该系统存在唯一。
A=Δt1λΔ0−∇⋅M∇I∇ϕn⋅∇pn+10Δt1I−νΔ+un⋅∇
取预条件子为
A=Δt1λΔ0−∇⋅M∇I000Δt1I−νΔ
可以利用2阶NS方程的方法构造上述二阶格式
⎩⎨⎧2Δt3u~n+1−4un+un−1+(2un−un−1)⋅∇u~n+1=νΔu~n+1−∇(2pn−pn−1)u~n+1∣Ω=0
⎩⎨⎧2Δt3(un+1−u~n+1)+∇(pn+1−2pn+pn−1)=0∇⋅un+1=0un+1⋅n∣Ω=0
上面格式对速度为2阶但是压力是1阶,如果要把压力提高半阶,
可以引入一个变量ψ,并且这种做法可以在H1模下
不改进
∣∣u−un∣∣L2(0,T,L2)+Δt∣∣u−un∣∣L2(0,T,H1)+Δt∣∣p−pn∣∣L2(0,T,L2)
⎩⎨⎧Δtϕn+1−ϕn+(u⋆n⋅∇)ϕn=MΔμn+1μn+1=λ(−Δϕn+1+ξn+1F′(ϕn)),rn+1Δtrn+1−rn=ξn+1(F′(ϕn),Δtϕn+1−ϕn)ξn+1=∫ΩF(ϕn)+C0rn+1
u⋆n=un+Δtμn+1∇ϕn
Δtu~n+1−u~n+un⋅∇u~n+1=νΔu~n+1−∇pn+μn+1∇ϕn
⎩⎨⎧Δtun+1−u~n+1+∇(pn+1−pn)∇⋅un+1=0un+1⋅n∣∂Ω=0
注意到
μ∇ϕ+ϕ∇μ=∇(ϕμ)
因此μ∇ϕ与ϕ∇μ本质是一样的,文献中可能会不一样。
但是这种格式有变系数,并且无法构造二阶格式
修正:引入两个SAV的变量
在可分区域(球,圆柱等)都有快速算法
有限差分
FFT
谱方法
罗威某教授写了个python包 shenfun(高度并行)
https://github.com/spectralDNS/shenfun
如果ρ1=ρ2,但是密度差不是很大的情况
可以做一个逼近,例如Boussinessq逼近
ρ=2ρ1+ρ2
动量方程
ρ0(ut+u⋅∇u)=∇⋅ν∇u−∇p+f(x)
以及
f(x)=−(1+ϕ)g(ρ1−ρ2)−(1−ϕ)g(ρ1−ρ2)
ϕ={1,流体1−1,流体2
但是当两种密度相差很大的情况下,例如水和气想吃大约3个数量级,不能用,所以需要另外一种格式
质量守恒
ρt+∇⋅(ρu)=0
不可压
∇⋅u=0
密度
ρ=2ϕ+1ρ1+21−ϕρ1
与黏性系数
ν=2ϕ+1ν1+21−ϕν2
质量守恒与不可压方程不能同时满足,因此只能取其中一种格式,每种格式的模型是不一样的。
具体参考
http://math.purdure.edu/~shen/pub/SheY15.pdf
-
例如图像处理,自由能
E(ϕ)=∫Ω∣∇ϕ∣+2λ(ϕ−g)2dx
在能量里加上高阶稳定项目,高阶项前面乘以一个很小的参数;
因为∣∇ϕ∣求导之后容易出现分母为0,因此可以做修正
∣∇ϕ∣=∣∇ϕ∣2+ϵ
-
如果非线性部分出现在∇前面
E(ϕ)=∫Ω∣γ(ϕ)∇ϕ∣2+F(ϕ)dx
则可以加上稳定项
这种处理对强非线性很有效
-
此外还可以与其它物理方程结合起来,整个系统一定要保证能量耗时,例如增加温度方程,磁流体麦克斯韦方程。