计算不可压缩流体 -- 数学基础

推荐教材： a lecture in computational fluid dynamics of incompressible flow: mathematics, algorithms and implementations.

Author: J. M. McDonough, Depart. of Mechanical Engineering and Mathematics, University of Kentucky

第一本参考，李新亮老师的CFD课程，偏压缩流，差分格式/误差分析是所有CFD课都涉及的， 间断激波处理的格式估计是特有的。课件里没有不可压的部分。所以...

第二本参考： MIT OpenCourse Numerical Fluid Mechanics， 从解线性代数方程组，差分格式，到网格，有限元CFD都有，全面不细。

第三本参考：推荐教材。差分/有限体积法CFD，入手方程性质/解性质，交错网格上的空间离散，介绍MAC， 投影法，SIMPLE等。数学比较多。看完整个CFD框架清晰很多，也有很多数学搞不定。

另外参考基于有限元的CFD，数学理论比差分更坚实, 算法精度/数值问题的解答更有据可循，单要求更多的数学基础。代表作 open source project: libmesh

--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

一 数学基础

1.1 基于Fourier展开分析的微分方程求解

斯托姆-刘维尔问题的特征方程可形成一组完备函数空间基，构成Banach 空间。其中满足内积运算（标量积）的子空间即Hilbert空间（L^2)。Hilbert函数空间内的任意函数可按Fourier展开，并具有上限（Cauchy - Schwarz 不等式），表示为函数的2-范数有界，亦即能量范数有限。

Hilbert空间的某类子集，满足其导函数仍然在Hilbert空间并满足2-范数有界，称为Sobolev空间。 可见对于需要解函数的导数仍然有界的PDE（比如，二阶以上或非线性PDE），其解必须存在于Sobolev空间。例如，NS方程的解就存在于Sobolev空间，如果存在的话。

1.2 算子空间

du(x,y,z)/dt = f(u) with u(t=0) = u_0

上式微分空间描述： C^3(x, y, z ) 空间函数映射到 C^3 * C^1(o, infinite)时空域。

考察一阶微分方程的经典解： u(x, t) = c1 * int ( u_0 * exp(c2) ) from (-infinite） to (+inifinite)

即表示：通过Fourier变换（或其他线性变换）将初值u_0空间映射到C^3 * C^1 时空。

1.3 弱解

大部分PDE不存在经典解，主要是在某些求解区域（比如边界)不协调。对于这些奇点位置采用广义函数定义解，那么结合这个奇点解就构成全域的弱解。广义函数不是真正的函数，但具有无限可导性和紧致性（奇点区域以外广义函数值为0）。 （试） 函数空间连续性定义为积分连续（弱形式），不同于函数连续要求每个点微分存在。故试探解在弱形式下具有连续性。

1.4 伴随算子

定义微分算子 L = sum ( A_k(x) d^k / dx^k )

其伴随算子 L‘ = sum (-1)^k * d^k(A_k(x)) / dx^k )

举例： L = A d^2/ dx^2 + B d/dx + C, 弱形式解表述：

< Lu, v> = int( A u'' v) + int(B u' v) + int (C uv) = int ( (Av)'' u) - int( (Bv)' u) + int ( u (Cv)) = < u, L'v>

两点注意： 1 因为试探解（v) 具有紧致性，边界自动消掉； 2 if L = L' 自伴随, 此时对应BVP矩阵具有对称性，此为大型有限元方法的基础。

定义微分算子P， 其伴随算子P‘，那么微分方程 Pu = f 的弱解 是任意满足 < u*, P'v> = < f, v> 的 u* , 其中v 来自试探解空间, f满足弱形式连续性要求， 注意这里u* 是定义在全域上的哟，即有限元方法里的Ritz法，Galerkin法将全域上定义的弱解分片到单元上满足，即更弱的形式，单元边界只满足C0连续。

1.5 Galerkin过程

1 L^2函数空间中的某类系统（NS方程），对其中各项采用Fourier展开，基函数为空间坐标的函数；

2 展开后各项与Fourier所有基函数内积，单只截取前N个基函数的乘积表达式，以此求解N个Fouerier系数。

数值方法解得系数，并带回原因变量的Fourier多项式，就是原问题的数值解。

这里采用Fouerier展开构造原方程的弱形式，此即基于Fouerier展开数值求解微分方程的基本数学原理。

1.6 Helmholtz 分解与 Leray 投影

Helmoholtz分解，有界向量场 v 可表示为 v = grad_ phi + u, 其中 phi 无旋，u 散度为0 (diver_u = 0 )

Lapalace ( phi) = grad_ v :: Possion Equation

考虑固面边界（不可穿透条件), u * n =0 --> grad_phi * n = v * n --> d phi/ dn = v * n on B.C.

上述泊松方程加边界条件即构成了描述 phi 的 Newman问题。

解存在条件要求满足， 边界积分int( v * n) = 0 <---- 对于不可穿透，无滑移边界，满足。 即对一般意义上流动模型，标量 phi总是存在的；

然后考虑 u = v - grad_ phi （*)， 即可解速度场。

(*) 式 即 Leray Projector。当向量场u散度为0，其Lp 运算等于自身。即 LP（u） = u if diver_ u = 0

1.7 NS 方程 Lp投影 / 弱解

对2D不可压缩NS方程作LP运算， u, du/dt 为其自身， 压力梯度项为0， Laplace项的LP操作得到Stokes算子 Au， 对流项的LP操作得到B(u, v)。整理即

du/dt + miu * Au + B(u) = f

令 F = f - miu * Au - B(u)， 得到 u' = F 此即NS方程的LP格式，其中压力梯度项不再了...

通过Galerkin过程给定NS弱解格式。但是这里不要求试探函数v具有C^inifite_0 无穷连续性，实际上v存在于一阶Soblev空间即可，另附加条件 v 满足 div_v = 0

d/dt < u, v> + miu * < u, v> + b(u, u, v> = <f, v>

< grad_p , v> = p * div_v + p * v *n = 0 ( div_v =0 & no flux no slip)

可见 弱解格式也不含压力梯度项。