20190807个人笔记(川大)

官网课程视频地址

笔记主页

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

N-S(Navier-Stokes)方程

来源

动量守恒

ma=F ma = F

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

牛顿流体:

T=(σij)=μ(Du+DuT)+(λup)I T = (\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.
ρ(ut+uu)=μΔ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+uu=γΔup+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ω=0TnΩ=g \begin{cases} \bm u |_{\partial \omega} = 0 \\[1em] T \cdot \bm n|_{\partial \Omega} = \bm g \end{cases}
或者周期边界条件

引理

如果u=0,  unΩ=0\nabla \cdot \bm u=0, \;\bm u \cdot \bm n |_{\partial \Omega = 0},则有
Ωuuvdx=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)
可以推出能量方程
12ddtu2=γu2+(f,u) \frac{1}{2}\frac{d}{dt}||\bm u||^2 = \gamma || \nabla \bm u ||^2 + (f, \bm u)
可以证明出
uL(0,T;L2)L2(0,T;H01) \bm u \in L^\infty(0,T;L^2) \cap L^2(0, T;H_0^1)
强解满足

12ddtu2=(uu,Δu)uLuΔ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||
非线性分析和维数有关,线性分析可以不讨论维数

引理

uL{u12Δu12,d=2{u12Δu12u14Δu34,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}

引理 Holder 不等式

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

杨不等式

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

对于d=2d=2情形
u12uΔu32γ2Δu2+C(γ)uu4 \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:=ddtu2Cu4 y:=\frac{d}{dt} || \nabla \bm u|| ^2 \leq C|| \nabla \bm u||^4
yCy2=θ(t)y y '\leq Cy^2 = \theta(t)y
其中
θ(t)=Cu2with 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 (tT) y(t) \leq y(0)e^{\int \theta dt} \leq \text{一致有界Const } (t\leq T)

因此存在唯一强解
uL(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情形从略
但是强解依赖于
t12C0y2(0) t \leq \frac{1}{2C_0 y^2(0)}
C0C_0为粘性系数

离散

{un+1unΔt+unun=γΔun+1pn+1un+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}

惩罚方法

{αuϵΔuϵ+pϵ=fuϵ+ϵ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要比较小,而且不好选取,并且方程耦合了

分析


eϵ=uuϵ,qϵ=ppϵ. 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}})
ϵ2p2 \leq \frac{\epsilon}{2}|| p ||^2
如果加上 infsup条件
supvH01(v,q)vqβq>0,qL02 \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)

惩罚方法(迭代)

{αuϵnΔuϵn+pϵn=fuϵn+ϵpϵn=ϵpϵn1 \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}

分析


eϵn=uuϵn,qϵn=ppϵn e^n_{\epsilon} = \bm u- \bm u^n_{\epsilon},\quad q^n_{\epsilon} = p- p^n_{\epsilon}
根据误差方程
αeϵn2+eϵn2+ϵqϵn2=ϵ(qn1,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ϵn1 || \nabla e^n || \leq \epsilon^n ,\quad ||q^n|| \leq \epsilon^{n-1}

另一种方法1

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

A1p=A1f \nabla \cdot A^{-1} \nabla p = \nabla \cdot A^{-1}\bm f

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

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

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

{αuϵΔuϵ+pϵ=fuϵ+Δ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条件

理论分析

误差方程
αeϵn2+eϵn2+ϵqϵn2=ϵ(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}})

算子分裂法

tu=A1u+A2u \partial_t u = A_1 u + A_2 u
先走半步
un+12unΔt=A1un+1 \frac{u^{n+\frac{1}{2}} - u^n}{\Delta t} = A_1 u^{n+1}
再走半步
un+1un+12Δt=A2un+1 \frac{u^{n+1} - u^{n+\frac{1}{2}} }{\Delta t} = A_2 u^{n+1}
可以把这种想法应用与NS方程,这就是投影法的来源(Choin - Femam)
对于方程
{ut+uu=γΔupu=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+12unΔt+unun=γΔ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+1un+12Δt+pn+1=0 \frac{u^{n+1} - u^{n+\frac{1}{2}} }{\Delta t} +\nabla p^{n+1}= 0
边界: un+1nΩ=0u^{n+1} \cdot \bm n|_{\partial\Omega} = 0

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

压力校正法

第一步算一个预报
{u~n+1unΔt+unun=γΔu~n+1pnu~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+1u~n+1Δt+(pn+1pn)=0un+1=0un+1nΩ=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},可以看出精度

稳定性

非线性项:半隐式
{u~n+1unΔt+unu~n+1=γΔu~n+1pnu~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+12u~nu~n+1u~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+12(pn,u~n+1) R =- \gamma || \nabla \tilde u^{n+1}||^2 - (\nabla p^{n} , \tilde u^{n+1})

un+1+Δtpn+1=u~n+1+Δtpn u^{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+12 \tilde E^{n+1}(u,p) - \tilde E^n(u,p) = -2\Delta t\gamma|| \nabla \tilde u^{n+1}||^2
其中
E~(u,p)=12u2+Δt22p2 \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

改进

{u~n+1unΔt+unu~n+1=γΔu~n+1pn+1un+1=0un+1nΩ=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+1 p^{n+1} = p^n + \gamma \Delta t\Delta \psi^{n+1} = p^n + \gamma \nabla \cdot \tilde u^{n+1}
这样pn+1p^{n+1}的边界条件就会一直在变动

习题

证明原始投影法的稳定性

笔记主页