NS方程的求解
笔记主页
小技巧:
二阶格式需要两个初值,第一步可以用一阶格式计算,因为一阶格式截断误差是二阶所以不会影响最后的精度。
ma=F
ρdtdu=∇⋅T(内力)+f(外力)
T=(σij)=μ(Du+DuT)+(λ∇⋅u−p)I
其中Du=∂jui.
ρ(∂t∂u+u⋅∇u)=μΔu+(μ+λ)∇∇⋅u−∇ϕ+f
如果密度为常数ρ=ρ0=1,则有不可压条件
∇⋅u=0
于是可得方程
⎩⎨⎧∇⋅u=0ut+u⋅∇u=γΔu−∇p+f
其中γ=μ/ρ0
边界条件
⎩⎨⎧u∣∂ω=0T⋅n∣∂Ω=g
或者周期边界条件
如果∇⋅u=0,u⋅n∣∂Ω=0,则有
∫Ωu⋅∇uvdx=0
因此非线性项对能量没有贡献。
根据
(−∇p,u)=(p,∇⋅u)
可以推出能量方程
21dtd∣∣u∣∣2=γ∣∣∇u∣∣2+(f,u)
可以证明出
u∈L∞(0,T;L2)∩L2(0,T;H01)
强解满足
21dtd∣∣∇u∣∣2=(u⋅∇u,Δu)≤∣∣u∣∣L∞∣∣∇u∣∣∣∣Δu∣∣
非线性分析和维数有关,线性分析可以不讨论维数
∣∣u∣∣L∞≤⎩⎨⎧∣∣u∣∣21∣∣Δu∣∣21,d=2⎩⎨⎧∣∣∇u∣∣21∣∣Δu∣∣21∣∣u∣∣41∣∣Δu∣∣43,d=3
(u,v)=∣∣u∣∣p∣∣v∣∣q,p1+q1=1
ab≤ϵap+C(ϵ)bq
对于d=2情形
≤∣∣u∣∣21∣∣∇u∣∣∣∣Δu∣∣23≤2γ∣∣Δu∣∣2+C(γ)∣∣u∣∣∣∣∇u∣∣4
于是有
y:=dtd∣∣∇u∣∣2≤C∣∣∇u∣∣4
y′≤Cy2=θ(t)y
其中
θ(t)=C∣∣∇u∣∣2with ∫0Tθ(t)dt 有界
根据单调性
e−∫θdty(t)≤y(0)
所以
y(t)≤y(0)e∫θdt≤一致有界Const (t≤T)
因此存在唯一强解
u∈L∞(0,T;L2)∩L2(0,T;H01)
d=3情形从略
但是强解依赖于
t≤2C0y2(0)1
C0为粘性系数
⎩⎨⎧Δtun+1−un+un⋅∇un=γΔun+1−∇pn+1∇⋅un+1
⎩⎨⎧αuϵ−Δuϵ+∇pϵ=f∇⋅uϵ+ϵpϵ=0
好处是算子对称正定,坏处是ϵ要比较小,而且不好选取,并且方程耦合了
记
eϵ=u−uϵ,qϵ=p−pϵ.
可以得出
α∣∣eϵ∣∣2+∣∣∇eϵ∣∣2+ϵ∣∣qϵ∣∣2=ϵ(p,qϵ)
上式可以初略估计为o(ϵ21)
≤2ϵ∣∣p∣∣2
如果加上 infsup条件
v∈H01sup∣∣∇v∣∣∣q∣∣(∇⋅v,q)≥β∣∣q∣∣>0,∀q∈L02
可以证得估计为o(ϵ)
⎩⎨⎧αuϵn−Δuϵn+∇pϵn=f∇⋅uϵn+ϵpϵn=ϵpϵn−1
记
eϵn=u−uϵn,qϵn=p−pϵn
根据误差方程
α∣∣eϵn∣∣2+∣∣∇eϵn∣∣2+ϵ∣∣qϵn∣∣2=ϵ(qn−1,qn)
利用infsup条件,可以证明出
∣∣∇en∣∣≤ϵn,∣∣qn∣∣≤ϵn−1
Au+∇p=f
其中A=(αI−Δ)−1
上式同时取算子∇⋅A−1,则有
∇⋅A−1∇p=∇⋅A−1f
当α=0相当于是一阶算子,条件数很好,很容易算
但是当α=0时比较麻烦,因为条件数大
⎩⎨⎧αuϵ−Δuϵ+∇pϵ=f∇⋅uϵ+Δpϵ=0
同时给p增加边界条件
∂n∂pϵ∣∂Ω=0
这样处理之后就不是一个鞍点问题,因此可以直接求解,可以不需要infsup条件
误差方程
α∣∣eϵn∣∣2+∣∣∇eϵn∣∣2+ϵ∣∣∇qϵn∣∣2=ϵ(∇p,∇qϵ)
可以有
∣∣∇eϵ∣∣≤o(ϵ21)
∂tu=A1u+A2u
先走半步
Δtun+21−un=A1un+1
再走半步
Δtun+1−un+21=A2un+1
可以把这种想法应用与NS方程,这就是投影法的来源(Choin - Femam)
对于方程
⎩⎨⎧ut+u⋅∇u=γΔu−∇p∇⋅u=0
取
先走半步
Δtun+21−un+un⋅∇un=γΔun+1
边界为:un+21∣∂Ω=0
该方程容易求解
再走半步
Δtun+1−un+21+∇pn+1=0
边界: un+1⋅n∣∂Ω=0
对方程同取散度算子,则有
Δpn+1=Δt1∇⋅un+21
同时边界满足
∂n∂pn+1∣∂Ω=0
跟新un+1
un+1=un+21−Δt∇pn+1
这是早起的投影法,该方法是收敛的,但是精度不明确
第一步算一个预报
⎩⎨⎧Δtu~n+1−un+un⋅∇un=γΔu~n+1−∇pnu~n+1∣∂Ω=0
第二步校正
⎩⎨⎧Δtun+1−u~n+1+∇(pn+1−pn)=0∇⋅un+1=0un+1⋅n∣∂Ω=0
两个式子相加,可以消掉u~n+1,可以看出精度
非线性项:半隐式
⎩⎨⎧Δtu~n+1−un+un⋅∇u~n+1=γΔu~n+1−∇pnu~n+1∣∂Ω=0
上式子
2Δt1(∣∣u~n+1∣∣2−∣∣u~n∣∣−∣∣u~n+1−u~n∣∣)=R
其中
R=−γ∣∣∇u~n+1∣∣2−(∇pn,u~n+1)
由
un+1+Δt∇pn+1=u~n+1+Δt∇pn
可以证出这种格式是无条件稳定的,
E~n+1(u,p)−E~n(u,p)=−2Δtγ∣∣∇u~n+1∣∣2
其中
E~(u,p)=21∣∣u∣∣2+2Δt2∣∣∇p∣∣2
这种格式具有一阶精度,但是比较难证明,p的估计需要用到infsup条件
这种方法的缺陷是,对压力pn+1的边界条件本质上是初始边界条件
,因此对压力算不好
注:公式中的γ 其实就是ν
⎩⎨⎧Δtu~n+1−un+un⋅∇u~n+1=γΔu~n+1−∇pn+1∇⋅un+1=0un+1⋅n∣Ω=0
其中
pn+1=pn+γΔtΔψn+1=pn+γ∇⋅u~n+1
这样pn+1的边界条件就会一直在变动
证明原始投影法的稳定性
笔记主页