有限体积方法

有限体积方法
引言
有限体积法(FVM)是在物理空间上积分形式的守恒方程进行直接离散的数值方法。

与有限差分方法相比有限体积方法更具有一般性，适用于任意形式的网格，结构网格与非结构网格均适用。

有限体积法是一种基于将CFD中最基本的量在单元内的平均值，这是与有限差分及有限元方法区别的地方，后边两种方法的数值量都取为在网格点上。

FVM方法一个重要优势是跟守恒性离散这个重要的概念联系起来，它可以自动满足具有守恒性的离散。

另一个优点就是适用于任意的网格。

5．1 守恒性离散
对于量U守恒律的一般积分形式可以由式(1.1.1)给出如下
将上式的最终表面源项合并到通量项中得到
该表达式的基本特点是存在表面积分以及在体积内U的时间变化只依赖于表面上的通量. 如图5.1.1所示可将一个体积元分解成三个亚体元,对于每个亚体元写出守恒律表达式
将这些表面积分进行加和,内部线ADB以及DE总是两次出现,但是方向相反,将三部分积分守恒律相加,这些内部的贡献量就会抵消,只剩下外边界的贡献量.例如,对于有一个通量的贡献量
而对于也有一个相似的项:
这样这两项相加就可以抵消. 故要保证格式是守恒的,通量的数值离散必须满足这样一个基本性质.下面我们以一维守恒律的情形来说明这个问题
结合图5.1.2来说明这个问题
其中f是矢量通量的x方向分量, 参考上图, 定义一个一维有限体积网格,并把中间点定义为“单元面”. 例如, 对于元(i), 单元面就是i-1/2与i+1/2的中点.
对该有限体积网格应用中心差分, 在i, i+1与i-1点处分别离散得到
将以上三个方程加和就得到了与元AB(i-3/2, i+3/2)上的守恒律相容的离散方程,即
从上式可以看出内部点的通量贡献已经抵消掉. 有时这种特性称为通量项的“telescoping property”, 对于元AB, 只考虑中间点i(不考虑i-1与i+1点),则离散形式可以直接写为
从 5.1.7的两式对比中, 我们可以看到通量部分的离散具有统一的形式, 这就是我们所要强调的守恒的特性.
如果我们要考虑方程(5.1.3)的非守恒形式, 则通量的导数就可以写为
其中, a(u)为Jacobian函数, , 故非守恒形式可以写为
利用图5.1.2所示的有限体积网格, 对非守恒形式在i点应用二阶中心差分得到
其中是的值.
同样,对于i+1点以及i-1点有,
将9式中三个离散式子进行加和得到
参考5.1.7b, 将5.1.8式直接在AB上进行离散,可得
我们发现5.1.10a右边由元AB内部点贡献的通量部分并不能互相抵消掉, 表现出源项的特点,这导致计算机程序不能将之与物理源项相区分, 故非守恒形式的离散会产生内部源.这些项被认为在网格点处展开为一项的二阶形式. 在连续流情况下可以忽略它, 但是对于计算非连续流动,比如流动中有激波产生, 就会产生巨大的误差. 数值实验显示非守恒形式比守恒形式的精度更低,尤其是在遇到梯度大的地方,由于数值源项的存在会产生更大的误差.
5.1.1 守恒的离散化的正式表示方法
对于5.1.3式,如果离散成如下的形式就可以满足守恒性要求,
为数值通量, 其为u在(2k)个邻域内点的函数.
此外, 方程5.1.11与原方程相容性要求当所有的均相等时,有
这些都可以直接推广到多维的情形, 以上条件必须分别对矢量通量的所有分量均成立.
定理: 当趋近于0时,若离散方程5.1.11的解几乎处处收敛于某个函数
值, 则是方程 5.1.3的弱解(可以存在有限个间断——Rankine-Hugoniot条件)
5.2 有限体积方法基础
有限体积方法是积分形式的守恒律方程的直接离散,这是有限体积方法与有限差分方法最大的区别,由于积分形式是守恒律的最一般的表达式,它不要求通量一定是连续的,这就是有限体积方法接近真实流动的原因.
FVM需要按以下步骤来构建:
1.划分网格,由空间离散得到有限体积,一个控制体积与每一个网格点
都相关联
2.在每一个有限体积上应用积分形式的守恒律.
有限体积选择的条件
由于具有普遍性,有限体积方法能够适用于任何类型的网格,结构与非结构.
单元居中的方法: 未知量定义在网格单元的中心,网格线定义了有限体积及表面积, 此处, 变量与单元相关,如图5.2.1a及c. 流动变量是在整个单元的平均值, 可以认为是单元内部某些有代表性点的值, 例如单元中心点.
单元顶点的方法: 未知量定义在网格角上,此处变量与网格点相关,例如单元
顶点, 如图5.2.1b,d所示
在相容的有限体积方法的体积的选择上,以下的限制条件必须得到满足:
(i)它们的总数应该覆盖整个区域
(ii)亚区域是允许重叠的,条件是表面的每一部分作为一个偶数个不同亚区域的部分而出现,这样整体的总积分守恒律就适用于任何
相邻亚区域的组合域.
(iii)通量沿单元表面必须由不依赖于当地单元的公式来计算.
(iii)确保了守恒特性的满足,因为通量的内部边界的贡献量会抵消掉(相关的有限体积相加之后)
5.2.2 有限体积离散的定义
将积分型守恒律应用到每一个控制体积, 关联到网格点J, 因此对于依附于该点或单元上未知量的离散化方程可定义为:
该方法的优点(对于无源项方程尤其有优势)是通量只在二维表面上计算,而不是三维空间中. 5.2.1可由其离散形式代替,对于
参考图5.2.1a对于单元1(i, j), 用统一表示, 是ABCD面. 通量项在4条边AB, BC, CD, DA求和.
式(5.2.2) 说明了有限体积法与有限差分及有限元法区别的一些重要特性: 1.点J的坐标是变量的准确位置,在控制体积内它将不会明确的标出.因此
在控制体内联结到一个固定点,将之看作是整个控制元上该流动变量U的
一个平均值(图5.2.1a). 5.2.2式中第一项代表在选定的控制体积上流动变量的平均值的时间变化率.
2.网格坐标只出现在确定的单元体积以及侧面上. 因此, 参照图5.2.1a, 考察点
1的控制元ABCD只有A,B,C,D的坐标将是需要用到的.
3.当不存在源项的时候,有限体积方程式表示时间间隔内U的平均值变化等
于相邻单元之间通量的交换量,对稳态流动,得到的数值解是通量进入控制体平衡的结果, 即
例子: 图5.2.1a中AB面,对于1则通量贡献量为正,而8则为负.
4. 有限体积也允许边界条件的自然引入, 例如固壁, 法向分量为0, 对连续方
程. 在固壁处. 因此对(5.2.2)及(5.2.3)的相应的贡献变为0.
5.2.3 数值格式的一般表达式
假设守恒律的积分形式(5.2.1)对于控制体积, 从到进行积分有,
引入单元平均守恒变量, 在时间的源, 单元与时间平均源, 以及
每个边上的数值通量, 分别定义如下
守恒的离散化采用如下形式:
其中与任何网格点无关, 它是整个单元上的平均. 为了在离散化的水平上实现守恒,在给定的单元面上的数值通量的估计必须独立于其所属的单元.
如果考虑空间离散完全由其数值通量来定义,时间积分项暂不处理,则以上的数值方法就会得到其一般形式. 一个一般的数值格式可以定义为对时间的常微分方程为
定义残差为整个单元上的通量平衡减去源项贡献. 5.2.6是5.2.7的时间的向前差分,也有其他的时间离散方法,例如龙格-库塔法.
守恒性条件可选择的公式
在任意数量的单元上对5.1.2进行展开, J=1-N. 对所有的单元进行加和,削去所有单元内表面的贡献项得,
定义为在整个单元的平均值,该格式的守恒性要求,在每一个时间步,如下的条件要得到满足,
边界以及源项
5.3 有限体积方法的实际应用
5.3.1 二维有限体积方法
如图5.2.1a,考察控制单元ABCD, 方程5.2.1可写为
f, g为矢量通量F的直角坐标分量,对边AB,表面矢量可定义为
对于单元,可以得到有限体积方程
ABCD展开求和包括ABCD的四条边,对于一般的四边形,面积可由对角线矢量乘积表示,如图5.3.1, 平行四边形1234的面积是ABCD的两倍,因此
, 为点A的位置矢量.
对于单元ABCD,上式右边为正.
通过单元表面通量的计算
沿侧边通量分量,如的计算
(a)对于中心离散格式以及单元中心化的有限体积方法,有以下做法:
1.通量平均
2.由于通量分量一般是U的线性函数,以下的式子与5.
3.5不等价
3.将取为通量在A及B处的平均
这里,可以在A及B处求变量的值, 例如
以及
也可以进行通量的直接平均, 例如:
可以看到, 5.3.7与5.3.10比5.3.8与5.3.11需要更少的通量计算
(b)对于中心格式及单元-顶点的有限体积方法:
5.3.7, 5.3.8是对通量的直接近似, 5.3.8是对应着对积分梯形公式的应用
通过加和在单元ABCD四个边积分的贡献量(如图5.2.1b), 可以得到
例子: 在笛卡儿网格下的中心离散格式. 在笛卡儿坐标, 均匀网格下,上述有限体积公式与有限差分的公式是一致的. 由
可以得到(此处记, 同样其他的量也采用类似的记法)
两边除以可以得到中心差分格式
将5.3.5式应用到图5.2.1a, 方程E5.3.3变为
而由5.3.8与5.3.11将推出如下的公式
(c)单元-中心化有限体积的迎风格式(利用上游点求下游)
对流通量以相关的对流速度传播方向的函数来计算,
其中
由图5.2.1a可以定义
(d)对于迎风-单元顶点的有限体积方法(图5.2.1b), 可以定义
例子:E5.3.2 “笛卡儿坐标网格中的迎风格式”
考虑二维线性对流方程的离散
如图5.2.1a所示, 在单元ABCD应用有限体积的公式:
通量定义为, 选择方程5.3.14, AB,CD为竖直边,有
对于水平边BC, DA有
故其得到的格式为一阶迎风格式的推广, 具有一阶精度
5.3.2 梯度的有限体积的计算
对于任意一个体积，由高斯定理得
此处，S是封闭的边界表面,定义平均化的梯度为
以及
对于二维控制单元,可以得到
如图5.2.1d, 在公式两边应用梯形积分公式, 得到
此处对所有的顶点求和,从1到6, 以及. 经过整理可得到
对于y同样存在这样的关系
计算单元面积: 当U=x时,方程5.3.21左侧为1. 对于任意一个单元的面积可用如下的式子进行计算,
对任意一个四边形ABCD, 如图5.3.2, 可以得到
以及
对于y方向导数有,
对于同一单元的封闭面与体有如下关系
对于二维单元, 取, 可以推出如下的公式
例子: E5.3.3 二维扩散方程
考虑二维扩散方程
对于扩散的通量分量(k为常数) 在图5.2.1a的网格上进行有限体积的离散,将整个单元ABCD的通量表示如下,
在单元的格点A,B上计算导数, 对于单元(i, j),方程5.3.3可写为
对于A点, U的导数取整个元1,6,7,8的平均值,由5.3.26得
对于B可以得到与A类似的关系
通过边AB对通量的贡献为E5.3.14与E5.3.15两式的加和, 并与相乘而得到的,同理通过BC通量的贡献为
类似的关系对于C有
最后,对于方程E5.3.13有, 可以写为
该数值离散对应的是图4.2.3中Laplace算子的离散.
更简单的办法为
这样就推出了扩散方程的标准有限差分格式(对应图4.2.2)
以上可以推广到多维的情况以及流行的结构及非结构网格上去.。