• 3.5. 算子分裂法
  • 4. 压力校正法
  • 5. 习题
  • {{article.mytitle}}

    • NS方程的求解

    1. 笔记主页

    笔记主页

    小技巧：
    二阶格式需要两个初值，第一步可以用一阶格式计算，因为一阶格式截断误差是二阶所以不会影响最后的精度。

    2. N-S(Navier-Stokes)方程

    2.1. 来源

    2.1.1. 动量守恒

    ma=Fma = F

    ρdudt=∇⋅T(内力)+f(外力)\rho \frac{d \bm u}{ dt} = \nabla \cdot T \text{(内力)}+ \bm f\text{(外力)}

    2.1.2. 牛顿流体:

    T=(σij)=μ(Du+DuT)+(λ∇⋅u−p)IT = (\sigma_{ij}) = \mu(D\bm u + D \bm u^T) +(\lambda \nabla \cdot \bm u -p)I
    其中Du=∂jui.D \bm u = \partial_j u_i.
    ρ(∂u∂t+u⋅∇u)=μΔu+(μ+λ)∇∇⋅u−∇ϕ+f\rho (\frac{\partial \bm u}{\partial t} + \bm u \cdot \nabla \bm u) = \mu \Delta \bm u + (\mu + \lambda)\nabla \nabla\cdot \bm u - \nabla \phi + \bm f
    如果密度为常数ρ=ρ0=1\rho= \rho_0=1，则有不可压条件
    ∇⋅u=0\nabla \cdot \bm u =0
    于是可得方程
    {∇⋅u=0ut+u⋅∇u=γΔu−∇p+f\begin{cases} \nabla \cdot \bm u =0\\[1em] \bm u_t + \bm u \cdot \nabla \bm u = \gamma \Delta \bm u- \nabla p + \bm f \end{cases}
    其中γ=μ/ρ0\gamma = \mu /\rho_0
    边界条件
    {u∣∂ω=0T⋅n∣∂Ω=g\begin{cases} \bm u |_{\partial \omega} = 0 \\[1em] T \cdot \bm n|_{\partial \Omega} = \bm g \end{cases}
    或者周期边界条件

    2.2. 引理

    如果∇⋅u=0,  u⋅n∣∂Ω=0\nabla \cdot \bm u=0, \;\bm u \cdot \bm n |_{\partial \Omega = 0}，则有
    ∫Ωu⋅∇uvdx=0\int_{\Omega} \bm u\cdot \nabla \bm u\bm v dx = 0
    因此非线性项对能量没有贡献。

    根据
    (−∇p,u)=(p,∇⋅u)(-\nabla p, \bm u) = (p, \nabla \cdot \bm u)
    可以推出能量方程
    12ddt∣∣u∣∣2=γ∣∣∇u∣∣2+(f,u)\frac{1}{2}\frac{d}{dt}||\bm u||^2 = \gamma || \nabla \bm u ||^2 + (f, \bm u)
    可以证明出
    u∈L∞(0,T;L2)∩L2(0,T;H01)\bm u \in L^\infty(0,T;L^2) \cap L^2(0, T;H_0^1)
    强解满足

    12ddt∣∣∇u∣∣2=(u⋅∇u,Δu)≤∣∣u∣∣L∞∣∣∇u∣∣∣∣Δu∣∣\frac{1}{2}\frac{d}{dt}||\nabla \bm u||^2 = (\bm u\cdot\nabla \bm u, \Delta \bm u) \leq||\bm u||_{L^\infty} ||\nabla \bm u||||\Delta \bm u||
    非线性分析和维数有关，线性分析可以不讨论维数

    2.3. 引理

    ∣∣u∣∣L∞≤{∣∣u∣∣12∣∣Δu∣∣12,d=2{∣∣∇u∣∣12∣∣Δu∣∣12∣∣u∣∣14∣∣Δu∣∣34,d=3||\bm u||_{L^\infty} \leq \begin{cases} || \bm u ||^{\frac{1}{2}} || \Delta \bm u||^{\frac{1}{2}} , \quad d=2\\[1em] \begin{cases} || \nabla\bm u ||^{\frac{1}{2}} || \Delta \bm u||^{\frac{1}{2}} \\[1em] || \bm u ||^{\frac{1}{4}} || \Delta \bm u||^{\frac{3}{4}} \end{cases} ,\quad d = 3 \end{cases}

    2.4. 引理 Holder 不等式

    (u,v)=∣∣u∣∣p∣∣v∣∣q,1p+1q=1(u,v) = ||u||_p || v ||_q , \quad \frac{1}{p} + \frac{1}{q} =1

    2.5. 杨不等式

    ab≤ϵap+C(ϵ)bqab \leq \epsilon a^p + C(\epsilon)b^q

    对于d=2d=2情形
    ≤∣∣u∣∣12∣∣∇u∣∣∣∣Δu∣∣32≤γ2∣∣Δu∣∣2+C(γ)∣∣u∣∣∣∣∇u∣∣4\leq || \bm u||^{\frac{1}{2}} ||\nabla \bm u || || \Delta \bm u||^{\frac{3}{2}}\leq \frac{\gamma}{2}||\Delta \bm u||^2 + C(\gamma)||\bm u|| ||\nabla \bm u ||^4
    于是有
    y:=ddt∣∣∇u∣∣2≤C∣∣∇u∣∣4y:=\frac{d}{dt} || \nabla \bm u|| ^2 \leq C|| \nabla \bm u||^4
    y′≤Cy2=θ(t)yy '\leq Cy^2 = \theta(t)y
    其中
    θ(t)=C∣∣∇u∣∣2with ∫0Tθ(t)dt 有界\theta(t) = C||\nabla \bm u ||^2 \text{with } \int_{0}^T \theta(t) dt \text{ 有界}

    根据单调性
    e−∫θdty(t)≤y(0)e^{-\int \theta dt} y(t) \leq y(0)
    所以
    y(t)≤y(0)e∫θdt≤一致有界Const (t≤T)y(t) \leq y(0)e^{\int \theta dt} \leq \text{一致有界Const } (t\leq T)

    因此存在唯一强解
    u∈L∞(0,T;L2)∩L2(0,T;H01)\bm u \in L^\infty(0,T;L^2) \cap L^2(0, T;H_0^1)

    d=3d=3情形从略
    但是强解依赖于
    t≤12C0y2(0)t \leq \frac{1}{2C_0 y^2(0)}
    C0C_0为粘性系数

    3. 离散

    {un+1−unΔt+un⋅∇un=γΔun+1−∇pn+1∇⋅un+1\begin{cases} \frac{\bm u^{n+1} - \bm u^n}{\Delta t}+\bm u^n \cdot \nabla \bm u^{n} = \gamma \Delta \bm u^{n+1} - \nabla p^{n+1}\\[1em] \nabla \cdot \bm u^{n+1} \end{cases}

    3.1. 惩罚方法

    {αuϵ−Δuϵ+∇pϵ=f∇⋅uϵ+ϵpϵ=0\begin{cases} \alpha \bm u_\epsilon - \Delta\bm u_\epsilon + \nabla p_\epsilon = f\\[1em] \nabla \cdot \bm u_{\epsilon} + \epsilon p_{\epsilon} = 0 \end{cases}
    好处是算子对称正定，坏处是ϵ\epsilon要比较小，而且不好选取，并且方程耦合了

    3.1.1. 分析

    记
    eϵ=u−uϵ,qϵ=p−pϵ.e_{\epsilon} = \bm u - \bm u_{\epsilon},\quad q_{\epsilon} = p - p_{\epsilon}.
    可以得出
    α∣∣eϵ∣∣2+∣∣∇eϵ∣∣2+ϵ∣∣qϵ∣∣2=ϵ(p,qϵ)\alpha || e_{\epsilon} ||^2 + || \nabla e_{\epsilon}||^2 + \epsilon||q _{\epsilon}|| ^2 = \epsilon (p, q_{\epsilon})
    上式可以初略估计为o(ϵ12)o(\epsilon^{\frac{1}{2}})
    ≤ϵ2∣∣p∣∣2\leq \frac{\epsilon}{2}|| p ||^2
    如果加上 infsup条件
    sup⁡v∈H01(∇⋅v,q)∣∣∇v∣∣∣q∣∣≥β∣∣q∣∣>0,∀q∈L02\sup_{v\in H_0^1}\dfrac{(\nabla \cdot v, q)}{|| \nabla v|||q||} \geq \beta||q||>0,\quad \forall q \in L_0^2
    可以证得估计为o(ϵ)o(\epsilon)

    3.2. 惩罚方法(迭代)

    {αuϵn−Δuϵn+∇pϵn=f∇⋅uϵn+ϵpϵn=ϵpϵn−1\begin{cases} \alpha \bm u^n_\epsilon - \Delta\bm u^n_\epsilon + \nabla p^n_\epsilon = f\\[1em] \nabla \cdot \bm u^n_{\epsilon} + \epsilon p^n_{\epsilon} = \epsilon p^{n-1}_{\epsilon} \end{cases}

    3.2.1. 分析

    记
    eϵn=u−uϵn,qϵn=p−pϵne^n_{\epsilon} = \bm u- \bm u^n_{\epsilon},\quad q^n_{\epsilon} = p- p^n_{\epsilon}
    根据误差方程
    α∣∣eϵn∣∣2+∣∣∇eϵn∣∣2+ϵ∣∣qϵn∣∣2=ϵ(qn−1,qn)\alpha || e^n_{\epsilon} ||^2 + || \nabla e^n_{\epsilon}||^2 + \epsilon||q^n _{\epsilon}||^2 = \epsilon(q^{n-1}, q^n)
    利用infsup条件，可以证明出
    ∣∣∇en∣∣≤ϵn,∣∣qn∣∣≤ϵn−1|| \nabla e^n || \leq \epsilon^n ,\quad ||q^n|| \leq \epsilon^{n-1}

    3.3. 另一种方法1

    Au+∇p=fA \bm u + \nabla p = \bm f
    其中A=(αI−Δ)−1A = (\alpha I - \Delta )^{-1}
    上式同时取算子∇⋅A−1\nabla \cdot A^{-1}，则有

    ∇⋅A−1∇p=∇⋅A−1f\nabla \cdot A^{-1} \nabla p = \nabla \cdot A^{-1}\bm f

    当α=0\alpha =0相当于是一阶算子，条件数很好，很容易算

    但是当α≠0\alpha \neq0时比较麻烦，因为条件数大

    3.4. 另一种方法(压力稳定法)

    {αuϵ−Δuϵ+∇pϵ=f∇⋅uϵ+Δpϵ=0\begin{cases} \alpha \bm u_\epsilon - \Delta\bm u_\epsilon + \nabla p_\epsilon = f\\[1em] \nabla \cdot \bm u_{\epsilon} + \Delta p_{\epsilon} = 0 \end{cases}
    同时给pp增加边界条件
    ∂∂npϵ∣∂Ω=0\frac{\partial }{\partial \bm n}p_\epsilon |_{\partial \Omega} = 0

    这样处理之后就不是一个鞍点问题，因此可以直接求解，可以不需要infsup条件

    3.4.0.1. 理论分析

    误差方程
    α∣∣eϵn∣∣2+∣∣∇eϵn∣∣2+ϵ∣∣∇qϵn∣∣2=ϵ(∇p,∇qϵ)\alpha || e^n_{\epsilon} ||^2 + || \nabla e^n_{\epsilon}||^2 + \epsilon||\nabla q^n _{\epsilon}||^2 = \epsilon(\nabla p, \nabla q_{\epsilon})
    可以有
    ∣∣∇eϵ∣∣≤o(ϵ12)||\nabla e_\epsilon|| \leq o(\epsilon^{\frac{1}{2}})

    3.5. 算子分裂法

    ∂tu=A1u+A2u\partial_t u = A_1 u + A_2 u
    先走半步
    un+12−unΔt=A1un+1\frac{u^{n+\frac{1}{2}} - u^n}{\Delta t} = A_1 u^{n+1}
    再走半步
    un+1−un+12Δt=A2un+1\frac{u^{n+1} - u^{n+\frac{1}{2}} }{\Delta t} = A_2 u^{n+1}
    可以把这种想法应用与NS方程，这就是投影法的来源(Choin - Femam)
    对于方程
    {ut+u⋅∇u=γΔu−∇p∇⋅u=0\begin{cases} \bm u_t + \bm u \cdot \nabla \bm u = \gamma \Delta \bm u- \nabla p \\[1em] \nabla \cdot \bm u =0\\ \end{cases}
    取
    先走半步
    un+12−unΔt+un⋅∇un=γΔun+1\frac{u^{n+\frac{1}{2}} - u^n}{\Delta t} + u^n \cdot \nabla u^{n}= \gamma \Delta u^{n+1}
    边界为：un+12∣∂Ω=0u^{n+\frac{1}{2}}|_{\partial \Omega} = 0

    该方程容易求解

    再走半步
    un+1−un+12Δt+∇pn+1=0\frac{u^{n+1} - u^{n+\frac{1}{2}} }{\Delta t} +\nabla p^{n+1}= 0
    边界： un+1⋅n∣∂Ω=0u^{n+1} \cdot \bm n|_{\partial\Omega} = 0

    对方程同取散度算子，则有
    Δpn+1=1Δt∇⋅un+12\Delta p^{n+1} = \frac{1}{\Delta t} \nabla\cdot u^{n+\frac{1}{2}}
    同时边界满足
    ∂pn+1∂n∣∂Ω=0\frac{\partial p^{n+1}}{\partial \bm n}|_{\partial \Omega} = 0
    跟新un+1u^{n+1}
    un+1=un+12−Δt∇pn+1u^{n+1} = u^{n+\frac{1}{2}} - \Delta t \nabla p^{n+1}
    这是早起的投影法，该方法是收敛的，但是精度不明确

    4. 压力校正法

    第一步算一个预报
    {u~n+1−unΔt+un⋅∇un=γΔu~n+1−∇pnu~n+1∣∂Ω=0\begin{cases} \frac{\tilde u^{n+1} - u^{n}}{\Delta t} + u^n \cdot \nabla u^{n} = \gamma \Delta \tilde u^{n+1} - \nabla p^{n} \\[1em] \tilde u^{n+1} |_{\partial \Omega} = 0 \end{cases}

    第二步校正
    {un+1−u~n+1Δt+∇(pn+1−pn)=0∇⋅un+1=0un+1⋅n∣∂Ω=0\begin{cases} \frac{u^{n+1} - \tilde u^{n+1}}{\Delta t} + \nabla (p^{n+1} - p^n) = 0 \\[1em] \nabla \cdot u^{n+1} = 0 \\[1em] u^{n+1} \cdot \bm n|_{\partial \Omega} = 0 \end{cases}

    两个式子相加，可以消掉u~n+1\tilde u^{n+1}，可以看出精度

    4.1. 稳定性

    非线性项：半隐式
    {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} = \gamma \Delta \tilde u^{n+1} - \nabla p^{n} \\[1em] \tilde u^{n+1} |_{\partial \Omega} = 0 \end{cases}
    上式子
    12Δt(∣∣u~n+1∣∣2−∣∣u~n∣∣−∣∣u~n+1−u~n∣∣)=R\frac{1}{2\Delta t} \left ( ||\tilde u^{n+1}||^2 - || \tilde u^{n}|| - || \tilde u^{n+1} - \tilde u^n ||\right) = R
    其中

    R=−γ∣∣∇u~n+1∣∣2−(∇pn,u~n+1)R =- \gamma || \nabla \tilde u^{n+1}||^2 - (\nabla p^{n} , \tilde u^{n+1})
    由
    un+1+Δt∇pn+1=u~n+1+Δt∇pnu^{n+1} + \Delta t \nabla p^{n+1} = \tilde u^{n+1} + \Delta t\nabla p^n

    可以证出这种格式是无条件稳定的，
    E~n+1(u,p)−E~n(u,p)=−2Δtγ∣∣∇u~n+1∣∣2\tilde E^{n+1}(u,p) - \tilde E^n(u,p) = -2\Delta t\gamma|| \nabla \tilde u^{n+1}||^2
    其中
    E~(u,p)=12∣∣u∣∣2+Δt22∣∣∇p∣∣2\tilde E(u,p) = \frac{1}{2}||u||^{2} + \frac{\Delta t^2}{2}|| \nabla p ||^2

    这种格式具有一阶精度，但是比较难证明，pp的估计需要用到infsup条件

    这种方法的缺陷是，对压力pn+1p^{n+1}的边界条件本质上是初始边界条件
    ，因此对压力算不好

    注：公式中的γ\gamma 其实就是ν\nu

    4.2. 改进

    {u~n+1−unΔt+un⋅∇u~n+1=γΔu~n+1−∇pn+1∇⋅un+1=0un+1⋅n∣Ω=0\begin{cases} \dfrac{\tilde u^{n+1} - u^{n}}{\Delta t} + u^n \cdot \nabla \tilde u^{n+1} = \gamma \Delta \tilde u^{n+1} - \nabla p^{n+1} \\[1em] \nabla \cdot u ^{n+1}= 0 \\[1em] u^{n+1} \cdot n|_{\Omega} = 0 \end{cases}
    其中

    pn+1=pn+γΔtΔψn+1=pn+γ∇⋅u~n+1p^{n+1} = p^n + \gamma \Delta t\Delta \psi^{n+1} = p^n + \gamma \nabla \cdot \tilde u^{n+1}
    这样pn+1p^{n+1}的边界条件就会一直在变动

    5. 习题

    证明原始投影法的稳定性

    笔记主页