• {{article.mytitle}}

    • SAV求解NS方程

    1.1. 笔记主页

    笔记主页

    1.2. 复习回顾

    • 最原始的投影法基于算子分裂法

    • 压力校正方法在上一步的基础上两方个程同时加上∇pn\nabla p^{n}，但是这种方法因为压力pp的边界有问题，导致压力算不好。

    • 目前投影法，在最原始的基础上假设一个新的函数∇ψn+1\nabla\psi^{n+1}，他的特点是边界上不完全为00，即法方向为00切方向不为00
      {u~n+1−unΔt+un⋅∇u~n+1=νΔu~n+1−∇pnu~n+1∣Ω=0\begin{cases} \dfrac{\tilde u^{n+1} - u^{n}}{\Delta t} + u^n \cdot \nabla \tilde u^{n+1} = \nu\Delta \tilde u^{n+1} - \nabla p^{n} \\[1em] \tilde u^{n+1} |_{\Omega} = 0 \end{cases}

    {un+1−u~n+1Δt+∇ψn+1=0∇⋅un+1=0un+1⋅n∣Ω=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+1p^{n+1} = p^n - \nu \nabla \cdot \tilde u^{n+1}

    1.3. 投影法BDF二阶格式

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

    但是上面的格式稳定性难以证明
    {3u~n+1−4un+un−12Δt+(2un−un−1)⋅∇u~n+1=νΔu~n+1−∇(pn+1−2Δt3Δ(pn+1−pn))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+1−u~n+1)2Δt+∇(pn+1−pn)=0∇⋅un+1=0un+1⋅n∣Ω=0un+1⋅τ∣∂Ω=2Δt3Δ(pn+1−pn)τ\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}

    1.4. 一阶格式

    {un+1−unΔt+un⋅∇un+1+12∇⋅unun+1=νΔun+1−∇pn∇⋅un+1−ϵΔpn+1=0,∂pn+1∂n∣∂Ω=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},
    注意到
    (un⋅∇un+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+1∣∣2−∣∣un∣∣2=−2Δtν∣∣∇un+1∣∣2−2Δ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+1∣∣2=02\Delta t(\nabla u^{n+1}, p^{n+1}) + 2\Delta t \epsilon || \nabla p^{n+1} ||^2 = 0
    上面两式合并起来
    要合理利用
    ∇⋅(un+1−un)−ϵΔ(pn+1−pn)=0\nabla \cdot (u^{n+1} - u^n ) - \epsilon \Delta (p^{n+1} - p^n) = 0
    与
    ϵ∣∣∇(pn+1−pn)∣∣2=(∇(pn+1−pn),  un+1−un)\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+1−u~nΔt=νΔu~n+1−∇pn∇⋅u~n+1−ΔtΔpn+1=0,∂pn+1∂n∣∂Ω=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+1−u(tn+1))∣∣≤CΔt12|| \nabla (\tilde u^{n+1} - u(t^{n+1})) || \leq C\Delta t^{\frac{1}{2}}

    利用恒等式

    (3a−4b+c,2a)=a2+(2a−b)2+(a−2b+c)2−(b2+(2b−c)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

    1.5. 把SAV方法用到NS方程上

    1.5.1. 模型问题

    {ut+u⋅∇u=νΔu−∇p∇⋅u=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}}}

    1.5.2. 方程修正

    于是方程可以修改为
    {ut+ξu⋅∇u=νΔu−∇pξ=r(t)∫Ω12u2dx+C02rrt=(ut,u)=(u,ut+ξu⋅∇u)∇⋅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+1−unΔt+ξn+1un⋅∇un=νΔun+1−∇pn+1ξn+1=rn+1∫Ω12(un+1)2dx+C02rn+1rn+1−rnΔt=(un+1,un+1−unΔt+ξn+1un⋅∇un)∇⋅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+1−unΔt+ξn+1un⋅∇un=νΔu~n+1−∇pnξn+1=rn+1∫Ω(un+1)2dx+C02rn+1rn+1−rnΔt=(u~n+1,u~n+1−unΔt+ξn+1un⋅∇un)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+1−u~n+1Δt+∇(pn+1−pn)=0∇⋅un+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+1−unΔt=un+1−u~n+1Δt+u~n+1−unΔ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}
    可以证明出稳定性

    1.6. 两相流

    两种流体混合在一起，假设这两种流体会互相排斥，例如油与水，
    油滴入水中如何模拟？
    引入一个指标函数
    ϕ(x,t)={1,流体1−1,流体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(ϕ)dxE(\phi) = \gamma\int \frac{1}{2}|\nabla \phi|^2 + \frac{1}{\epsilon^2}F(\phi)dx
    其中F(ϕ)=14(ϕ2−1)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δϕdx≠0\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方程

    1.6.1. 动量方程

    ρ(ut+u⋅∇u)=∇⋅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=μΔu−∇p−λ∇⋅(∇ϕ⊗∇ϕ)T = \mu \Delta u - \nabla p - \lambda \nabla \cdot (\nabla \phi \otimes\nabla \phi)

    1.6.2. 密度方程

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

    1.7. 习题

    1.7.1. 习题1

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

    笔记主页