20190808个人笔记(川大)

官网课程视频地址

笔记主页

复习回顾

{un+1u~n+1Δt+ψn+1=0un+1=0un+1nΩ=0 \begin{cases} \dfrac{u^{n+1} - \tilde u^{n+1}}{\Delta t} + \nabla \psi^{n+1} = 0 \\[1em] \nabla \cdot u ^{n+1}= 0 \\[1em] u^{n+1} \cdot n|_{\Omega} = 0 \end{cases}

pn+1=pnνu~n+1 p^{n+1} = p^n - \nu \nabla \cdot \tilde u^{n+1}

投影法BDF二阶格式

容易想到下面格式
{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}

但是上面的格式稳定性难以证明
{3u~n+14un+un12Δt+(2unun1)u~n+1=νΔu~n+1(pn+12Δt3Δ(pn+1pn))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 (p^{n+1} - \frac{2\Delta t}{3} \Delta( p^{n+1} - p^n) )\\[1em] \tilde u^{n+1} |_{\Omega} = 0 \end{cases}

{3(un+1u~n+1)2Δt+(pn+1pn)=0un+1=0un+1nΩ=0un+1τΩ=2Δt3Δ(pn+1pn)τ \begin{cases} \dfrac{3(u^{n+1} - \tilde u^{n+1})}{2\Delta t} + \nabla (p^{n+1} - p^n )= 0 \\[1em] \nabla \cdot u ^{n+1}= 0 \\[1em] u^{n+1} \cdot n|_{\Omega} = 0\\[1em] u^{n+1}\cdot \tau|_{\partial \Omega} = \frac{2\Delta t}{3} \Delta( p^{n+1} - p^n) \tau \end{cases}

一阶格式

{un+1unΔt+unun+1+12unun+1=νΔun+1pnun+1ϵΔpn+1=0,pn+1nΩ=0 \begin{cases} \dfrac{u^{n+1} - u^{n}}{\Delta t} + u^n \cdot \nabla u^{n+1} + \frac{1}{2}\nabla \cdot u^{n} u^{n+1} = \nu\Delta u^{n+1} - \nabla p^{n} \\[1em] \nabla \cdot u^{n+1} - \epsilon \Delta p^{n+1} = 0,\quad \frac{\partial p^{n+1}}{\partial \bm n}|_{\partial \Omega} = 0 \end{cases}
第一式同时与un+1u^{n+1},
注意到
(unun+1,un+1)=(unun+1,un+1)(unun+1,un+1) (u^n\cdot \nabla u^{n+1}, u^{n+1}) = -(\nabla \cdot u^{n} u^{n+1}, u^{n+1})- (u^nu^{n+1}, \nabla u^{n+1})
于是有
un+12un2=2Δtνun+122Δt(pn,un+1) ||u^{n+1}||^2 - || u^n ||^2=-2\Delta t \nu ||\nabla u^{n+1} ||^2 - 2\Delta t(\nabla p^n, u^{n+1})

2Δt(un+1,pn+1)+2Δtϵpn+12=0 2\Delta t(\nabla u^{n+1}, p^{n+1}) + 2\Delta t \epsilon || \nabla p^{n+1} ||^2 = 0
上面两式合并起来
要合理利用
(un+1un)ϵΔ(pn+1pn)=0 \nabla \cdot (u^{n+1} - u^n ) - \epsilon \Delta (p^{n+1} - p^n) = 0

ϵ(pn+1pn)2=((pn+1pn),  un+1un) \epsilon || \nabla (p^{n+1} - p^n) ||^2 = (\nabla (p^{n+1} - p^n), \; u^{n+1}-u^n)
最后可证出稳定性,要求ϵΔt\epsilon \geq \Delta t

结合(1)n+1+(2)n(1)^{n+1} + (2)^n
{u~n+1u~nΔt=νΔu~n+1pnu~n+1ΔtΔpn+1=0,pn+1nΩ=0 \begin{cases} \frac{\tilde u^{n+1} - \tilde u^n}{\Delta t} = \nu \Delta \tilde u^{n+1} - \nabla p^{n}\\[1em] \nabla \cdot \tilde u^{n+1} - \Delta t\Delta p^{n+1} = 0, \quad \frac{\partial p^{n+1}}{\partial n}|_{\partial \Omega} = 0 \end{cases}
结论
(u~n+1u(tn+1))CΔt12 || \nabla (\tilde u^{n+1} - u(t^{n+1})) || \leq C\Delta t^{\frac{1}{2}}

利用恒等式

(3a4b+c,2a)=a2+(2ab)2+(a2b+c)2(b2+(2bc)2) (3a-4b+c, 2a) = a^2 + (2a-b)^2 + (a-2b+c)^2 - (b^2 + (2b-c)^2)

要巧妙利用un+1=0\nabla \cdot u^{n+1} = 0

把SAV方法用到NS方程上

模型问题

{ut+uu=νΔupu=0,uΩ=0.\begin{cases} u_t + u \cdot \nabla u = \nu \Delta u - \nabla p\\[1em] \nabla \cdot u =0,\quad u|_{\partial \Omega} = 0. \end{cases}
并不是梯度流
SAV的实质是
μ=δEδϕ=Δϕ+ξF(ϕ) \mu = \frac{\delta E}{\delta \phi} = - \Delta \phi + \xi F'(\phi)

其中ξ=r(t)ΩFdx+C0\xi = \frac{r(t)}{\sqrt{\int_{\Omega}{ Fdx+C_0}}}
因此这里可以令
ξ=r(t)Ωu2dx+C0 \xi = \frac{r(t)}{\sqrt{\int_{\Omega}{ u^2dx+C_0}}}

方程修正

于是方程可以修改为
{ut+ξuu=νΔupξ=r(t)Ω12u2dx+C02rrt=(ut,u)=(u,ut+ξuu)u=0,uΩ=0\begin{cases} u_t + \xi u \cdot \nabla u = \nu \Delta u - \nabla p\\[1em] \xi = \frac{r(t)}{\sqrt{\int_{\Omega}{ \frac{1}{2}u^2dx+C_0}}}\\[1em] 2rr_t = (u_t, u)= (u, u_t +\xi u \cdot \nabla u)\\[1em] \nabla \cdot u =0,\quad u|_{\partial \Omega} = 0 \end{cases}

{un+1unΔt+ξn+1unun=νΔun+1pn+1ξn+1=rn+1Ω12(un+1)2dx+C02rn+1rn+1rnΔt=(un+1,un+1unΔt+ξn+1unun)un+1=0,uΩ=0\begin{cases} \frac{u^{n+1}- u^n}{\Delta t}+ \xi^{n+1} u^n \cdot \nabla u^n = \nu \Delta u^{n+1} - \nabla p^{n+1}\\[1em] \xi^{n+1} = \dfrac{r^{n+1}}{\sqrt{\int_{\Omega}{ \frac{1}{2}(u^{n+1})^2dx+C_0}}}\\[2em] 2r^{n+1}\dfrac{r^{n+1} - r^n}{\Delta t} = (u^{n+1}, \frac{u^{n+1}- u^n}{\Delta t} +\xi^{n+1} u^n \cdot \nabla u^{n})\\[2em] \nabla \cdot u^{n+1} =0,\quad u|_{\partial \Omega} = 0 \end{cases}
如何解耦
先算u~n+1\tilde u^{n+1}
{u~n+1unΔt+ξn+1unun=νΔu~n+1pnξn+1=rn+1Ω(un+1)2dx+C02rn+1rn+1rnΔt=(u~n+1,u~n+1unΔt+ξn+1unun)u~n+1Ω=0\begin{cases} \frac{\tilde u^{n+1}- u^n}{\Delta t}+ \xi^{n+1} u^n \cdot \nabla u^n = \nu \Delta\tilde u^{n+1} - \nabla p^{n}\\[1em] \xi^{n+1} = \dfrac{r^{n+1}}{\sqrt{\int_{\Omega}{ (u^{n+1})^2dx+C_0}}}\\[2em] 2r^{n+1}\dfrac{r^{n+1} - r^n}{\Delta t} = (\tilde u^{n+1}, \frac{\tilde u^{n+1}- u^n}{\Delta t} +\xi^{n+1} u^n \cdot \nabla u^{n})\\[2em] \tilde u^{n+1}|_{\partial \Omega} = 0 \end{cases}
再校正
{un+1u~n+1Δt+(pn+1pn)=0un+1=0\begin{cases} \dfrac{u^{n+1} - \tilde u^{n+1}}{\Delta t}+ \nabla(p^{n+1}- p^n)=0 \\[1em] \nabla \cdot u^{n+1} =0 \end{cases}
利用下式(不会丢精度)
un+1unΔt=un+1u~n+1Δt+u~n+1unΔt \frac{u^{n+1} - u^n}{\Delta t} = \frac{u^{n+1} - \tilde u^{n+1}}{\Delta t}+ \frac{\tilde u^{n+1} - u^n}{\Delta t}
可以证明出稳定性

两相流

两种流体混合在一起,假设这两种流体会互相排斥,例如油与水,
油滴入水中如何模拟?
引入一个指标函数
ϕ(x,t)={1,流体11,流体2 \phi(x,t) = \begin{cases} 1 ,\quad &\text{流体1}\\ -1, \quad &\text{流体2} \end{cases}
在界面处用一个光滑函数连接起来,并且界面层很薄,因此需要一个方程去描述ϕ\phi的运动,ϕ\phi的方程满足输运方程,uu是流体速度,例如
ϕt+uϕ=0 \phi_t + u\cdot \nabla \phi = 0
但是上面的格式有缺陷,容易丢失界面,导致与实验偏差很大,因此这种方法可以很好用于动画模拟。当然也可进行特殊处理,但是复杂

另外一种方法
用梯度流来表示传输性质
从一个能量出发,先引入一个能量
E(ϕ)=γ12ϕ2+1ϵ2F(ϕ)dx E(\phi) = \gamma\int \frac{1}{2}|\nabla \phi|^2 + \frac{1}{\epsilon^2}F(\phi)dx
其中F(ϕ)=14(ϕ21)2F(\phi) = \frac{1}{4}(\phi^2 - 1)^2
其平衡态的解为
Δϕ+1ϵ2F(ϕ)=0 -\Delta \phi + \frac{1}{\epsilon^2}F'(\phi) = 0
于是可以逼近
ϕt+uϕ=σδEδϕ \phi_t + u\cdot \nabla \phi = -\sigma\frac{\delta E}{\delta \phi}
边界层一直存在,相场模型可以保持界面不变,这是它的好处之一;第二个好处,如果界面有拓扑变换的话,这个方程自然而然就出来了。

两个体积保持不变,因此希望有
ddtΩϕdx=0 \frac{d}{dt}\int_{\Omega} \phi dx = 0
uu具有性质
u=0,uΩ=0 \nabla \cdot u = 0, \quad u|_{\partial \Omega} = 0
于是有
tΩϕdx=σδEδϕdx0 \partial_t \int_{\Omega} \phi dx = - \sigma \int \frac{\delta E}{\delta \phi}dx \neq 0
因此质量不守恒

改进办法为
加上变分导数边界条件
nδEδϕ=0 \frac{\partial }{\partial \bm n}\cdot \frac{\delta E}{\delta \phi}=0

于是有质量守恒

用AC方程不太好,因此用CH方程

动量方程

ρ(ut+uu)=T \rho (u_t + u\cdot \nabla u) = \nabla \cdot T
其中TT是合外力
T=μ(Du+DuT)pI+λϕϕ T = \mu (D u + Du ^T) \cdot pI + \lambda\nabla \phi \otimes\nabla \phi

这是牛顿流,上述方程右边后面一项是表面张力,它是张量

上式子同求散度\nabla \cdot可得
T=μΔupλ(ϕϕ) T = \mu \Delta u - \nabla p - \lambda \nabla \cdot (\nabla \phi \otimes\nabla \phi)

密度方程

假设 ρ=ρ0=1\rho = \rho_0 = 1

习题

习题1

SAV方法解NS方程,构造BDF-2格式,并证明稳定性

笔记主页