下载文档原格式
附录F 二维不可压缩黏性流体方腔流动问题的有限体积算法与计算程序二维方腔流动问题是一个不可压缩黏性流动中典型流动。
虽然目前尚不能求得它的解析解，但是它常被用来作为检验各种数值算法计算精度和可靠性的算例。
文献中几乎大多数算法都对它进行过计算。
在本算例中采用有限体积算法三阶迎风型QUICK 离散格式对它进行数值求解。
同时，为了初学者入门和练习方便,这里给出了用C 语言和Fortran77语言编写的计算二维不可压缩黏性方腔流动问题计算程序，供大家学习参考。
F-1利用有限体积算法三阶迎风型QUICK 离散格式求解二维不可压缩黏性流体方腔流动问题1.二维不可压缩黏性流体方腔流动问题二维不可压缩黏性流体方腔流动（cavity flow ）：有一正方形腔室，其量纲为一的宽度为1.0，里面充满静止的不可压缩黏性流体，方腔内初始时刻压力和密度为=1.0=1.0p ρ、它周围壁面（左右壁面和底面）固定不动，上壁面以量纲为一的速度1.0沿着上壁面方向自左向右运动（图F.1）。
2. 基本方程组、初始条件和边界条件设流体是黏性流体。
二维方腔流动问题在数学上可以由二维不可压缩黏性流动N - S()()u v x y x x φρφρφμ∂∂∂∂⎛⎫+= ⎪∂∂∂∂⎝⎭⎝⎭） 其中u 为变量φ在水平x 方向的流速，v 为φ在垂直y 方向的流速，μ为黏度，S φ为源项。
源项中不仅包含压力梯度项，也包含时间导数项。
初始条件：方腔上壁面以量纲为一的速度1.0沿着上壁面方向自左向右运动。
图F.1二维不可压缩黏性方腔流问题示意图边界条件：流动速度u v 、均可采用无滑移边界条件，压力p 采用自由输出边界条件。
3.计算网格划分和控制体单元与节点定义采用交错网格，图F.2和图F.3是计算网格、控制体单元和节点示意图。
节点P 所在主控制单元如图F.2中有阴影部分所示。
在x 方向与节点P 相邻的节点为W 和E ，在y 方向与节点P 相邻的节点为S 和N ，主控制单元界面分别为w s e n 、、、。
模糊控制技术在SIMPLER算法中的应用及求解性能分析王艳宁;孙东亮;苗政;陈家庆;蔡晓君【摘要】为了提高SIMPLER算法在三维流动问题上的求解性能,引入模糊控制方法来自动调控速度亚松弛因子的大小.在数值计算过程中,将相邻两个迭代层次上的最大动量残差比值作为模糊控制输入量,速度亚松弛因子的变化量作为模糊控制输出量,基于最大动量残差的变化趋势可实现速度亚松弛因子的自动调控,从而达到加快收敛的目的.最后,通过3个经典的流动问题验证了模糊控制方法的优越性.研究表明:当初始亚松弛因子为最不利值时,模糊控制方法的收敛速度约是固定松弛因子方法的5～30倍;当初始亚松弛因子为最佳值时,模糊控制方法迭代次数与固定松弛因子方法迭代次数之比为0.7～2.0,收敛速度相差不大;采用模糊控制方法后,SIMPLER 算法在不同初始亚松弛因子下均能得到高速收敛的解,同时健壮性也显著提高.研究工作将为大幅提升SIMPLER算法在三维流动问题上的求解性能起到重要作用.【期刊名称】《西安交通大学学报》【年(卷),期】2016(050)001【总页数】7页(P78-84)【关键词】模糊控制;SIMPLER算法;三维流动问题;亚松弛因子;收敛速度;健壮性【作者】王艳宁;孙东亮;苗政;陈家庆;蔡晓君【作者单位】华北电力大学可再生能源学院,102206,北京;北京石油化工学院机械工程学院,102617,北京;华北电力大学可再生能源学院,102206,北京;北京石油化工学院机械工程学院,102617,北京;北京石油化工学院机械工程学院,102617,北京【正文语种】中文【中图分类】TK124求解流动与传热问题的压力修正算法首次由著名学者Patankar和Spalding提出,并被命名为SIMPLE算法[1]。
为了克服SIMPLE算法中初始压力和速度场单独进行设定的缺点,随后Patankar提出了改进算法SIMPLER[2]。
目前,SIMPLE系列算法已被广泛应用于求解流动与传热问题[3-7]。
数值传热学数值传热学常用的数值方法1.有限差分法历史上最早采用的数值方法，对简单几何形状中的流动与换热问题最容易实施的数值方法。
其基本点是:将求解区域中用于坐标轴平行的一系列网格的交点所组成的点的集合来代替，在每个节点上，将控制方程中每一个导数用相应的差分表达式来代替，从而在每个节点上，形成一个代数方程，每个方程中包括了本节点及其附近一些节点上的未知值，求解这些代数方程就获得了所需的数值解。
2.有限容积法将所计算的区域划分成一系列控制容积划分为一系列控制容积，每个控制容积都有一个节点做代表。
通过将守恒型的控制方程对控制容积坐积分导出离散方程。
在导出过程中，需要对界面上的被求函数本身及其一阶导数的构成做出假定，是目前流动与换热问题的数值计算中应用最广的一种方法。
3.有限元法把计算区域划分为一系列原题(在二维情况下，元体多为三角形或四边形)，由每个元体上去数个点作为节点，然后通过对控制方程做积分来获得离散方程。
有限元法最大的优点是对不规则区域的适应性较好。
但计算的工作量一般要比有限容积法大，而且在求解流动与换热问题是，对流项的离散处理方法及不可压缩流体原始变量法求解方面没有有限容积法成熟。
4.有限分析法由陈景仁教授在1981年提出。
在这种方法中，也像有限差分法那样，用一系列网格线将区域离散，所不同的是每一个节点与相邻4个网格(二维)问题组成计算单元，即一个计算单元由一个中心节点与8个l 邻点组成。
在计算单元中把控制方程中的非线性项局部线性化，并对该单元上未知函数的变化型线作出假设，把所选定型线表达式中系数和常数项用单元边界节点上位置的变量值来表示，找出其分析解。
然后利用其分析解，得到该单元中点及其边界上的位置值的代数方程，即单元中点的离散方程。
有限容积法和有限体积法有限容积法和有限体积法是计算流体力学中常用的两种数值方法，它们在流体动力学的数值计算中占有非常重要的地位。
本文将从概念、原理、特点、应用等方面，对这两种方法进行详细介绍。
一、有限容积法1.概念有限容积法（Finite Volume Method，FVM）是一种离散化的数值方法，它将连续的物理量离散化为有限个体积元，在每个体积元内计算其平均值，进而求解整个流体系统的物理量。
FVM方法的核心是质量守恒原理，即物质的进出必须平衡，这种保证了物理量在每个体积元内的守恒关系，从而保证了数值计算的准确性。
2.原理FVM方法的数值计算是基于网格的，它将流体动力学问题离散化为一个由有限体积元组成的系统，将原问题转化为流量守恒方程的求解，即$$\frac{\Delta m}{\Delta t}=\Sigma_{faces}\rho uA$$其中，$\Delta m$是在$\Delta t$时间内通过一个表面的质量变化量，$\rho$是介质的密度，$u$是速度，$A$是面积。
对于每个有限体积元，上式可以写为其中，$F_{ij}^p$和$F_{ij}^n$分别是流向有限体积元内部和外部的通量，$i,j$是有限体积元的编号。
3.特点（1）FVM方法基于质量守恒原理，具有非常强的数值稳定性和保真性；（2）FVM方法的计算结果具有局部守恒性，能够准确反映流场内部的物理现象；（3）FVM方法可以处理非结构化网格，适用范围广泛；（4）FVM方法求解的是面积分，所需的时间和空间存储相对较少。
4.应用（1）流体力学领域，如空气动力学、水力学、燃烧问题等；（2）材料科学领域，如薄膜生长、材料变形等。
有限体积法（Finite Element Method，FEM）是一种离散化的数值方法，它将求解的物理场离散化为有限个单元，然后在每个单元内进行近似计算。
相比于FVM方法，FEM方法更加精确，适用于需要高精度计算的问题。
MAC算法计算二维方腔顶盖流动李江飞;石兆东;段兴华;李岩芳;张康;逯国强;陈颖超;任亚东【摘要】二维方腔流动是不可压缩黏性的典型流动,可以用来检验各种数值算法计算精度和可靠性,目前尚不能求得它的解析解.基于Matlab编程,采用交错网格MAC算法求解二维方腔流动,计算采用控制容积积分法离散控制方程,对流项和扩散项采用中心差分格式,得到流动达到稳定状态时各物理量的分布.【期刊名称】《宜宾学院学报》【年(卷),期】2015(015)006【总页数】4页(P28-31)【关键词】数值模拟;方腔流动;控制容积积分法;MAC算法;离散【作者】李江飞;石兆东;段兴华;李岩芳;张康;逯国强;陈颖超;任亚东【作者单位】承德石油高等专科学校热能工程系,河北承德 067000;承德石油高等专科学校热能工程系,河北承德 067000;承德石油高等专科学校热能工程系,河北承德 067000;承德石油高等专科学校热能工程系,河北承德 067000;承德石油高等专科学校热能工程系,河北承德 067000;承德石油高等专科学校热能工程系,河北承德067000;承德石油高等专科学校热能工程系,河北承德 067000;承德石油高等专科学校热能工程系,河北承德 067000【正文语种】中文【中图分类】TB126Li JF,ShiZD,Duan XH,etal.Calculation of Two-dimensionalCavity Flow Based on MAC[J].Journal of Yibin Univer⁃sity,2015,15（6）：28-31.二维不可压缩黏性流体方腔流动顶盖拖动速度为utop，方腔的长度和高度均为H,流体密度为ρ、动力粘度为μ.边界条件：流动速度u、v采用无滑移边界条件，利用动量方程推导压力p的边界条件[1].流动与传热的控制方程如下：其中，p为压力，u、v分别为x、y方向速度分量.用高度H、流体密度ρ和拖动速度utop作为无量纲标尺，将控制方程无量纲化，流场初始状态为静止，Re=1000求流动达到稳定状态时，x方向中垂线（x=H/2）上的无量纲速度U，y 方向中垂线（y=H/2）上的无量纲速度V，绘制出速度分布曲线，并求出中垂线上||U、||V的平均值.1.1 涡量控制方程无量纲化以高度H、流体密度ρ和速度utop作为无量纲标尺，将控制方程无量纲化[2]：将上述无量纲量代入题中流动与传热的控制方程，得出如下的无量纲方程：1.2 边界条件边界条件为：流动速度采用无滑移边界条件，壁面处法向速度恒为0，切向速度也为零.顶盖u=1，v=0；其余u=v=0.采用均分网格，网格数80×80的交错网格来离散方程，将压力和速度放在不同位置，压力放在网格中心，以Pi,j为主节点，背离P点的U、V与P点有相同的编号，速度分量U与P在X方向位置相错半个网格，速度分量V与P在Y方向位置相错半个网格，具体如下[3]：P：X方向：0—81，左边点0，右边点81，边点与内点距离为，其余ΔX；Y方向：0—81，下边点0，上边点81，边点与内点距离，其余ΔY；U：X方向：0—80，左边点0，右边点80，相邻两点距离ΔX；Y方向：0—81，下边点0，上边点81，边点与内点距离，其余ΔY；V：X方向：0—81，左边点0，右边点81，边点与内点距离，其余ΔX；Y方向：0—80，下边点0，上边点80，相邻两点距离ΔY.对于MAC算法而言，采用交错网格，用控制容积积分法离散控制方程，对流项和扩散项采用中心差分格式.时间步长为Δτ，空间步长为ΔX、ΔY.对速度分量U进行离散[4-5]：内点处理：非稳态项：对流项：扩散项：压力项：边界点处理：对于上边界点，扩散项：对于下边界点，扩散项：内点离散后的动量方程为：对速度分量V进行离散：内点处理：非稳态项：对流项：扩散项：压力项：对于左边界点，其扩散项：对于右边界点，其扩散项：离散后的动量方程为：将，代入离散的连续性方程，如下：整理化简可得压力离散方程：其中：aP=aE+aW+aN+aS,aE=aW= aN=与上边界相邻的内点：与下边界相邻的内点：与左边界相邻的内点：与右边界相邻的内点：求解步骤如下[6-8]：①确定网格信息，如空间步长、时间步长：ΔX，ΔY，Δτ；②定义变量，给速度场和压力场赋初始值和边界值；③经过（1）、（2），可得完整的速度场离散结果，据公式求；④根据（3）求解压力泊松方程，采用Gauss-Seidel迭代求解，循环直至满足收敛条件；⑤用该时层满足收敛条件最新的压力场去更新速度场，得到下一时层的，；⑥用下一时层的，返回（4），直到稳态的解，求出速度场和压力场.程序流程如图1所示.U、V、P放置在相同的网格位置时，取中心差分，切断相邻点间关系，计算结果会出现非物理意义的振荡；通过交错网格，引入相邻两点的压差，在物理上保证了相邻点压力的相互影响，避免了振荡的压力解.稳定状态时，计算结果如图2、图3所示.【相关文献】[1] Peng Y F,Shiau Y H,Hwang R R.Transtion in a 2-D lid-driven cavity flow[J].Computer＆Fluids,2003,32（3）：337-352.[2] 陶文铨.数值传热学[M].第二版.西安：西安交通大学出版社，2010.[3] Abdallah S.Numerical solutions for the pressure poisson equation with neumann boundary conditions using a non-staggeredgrid[J].Journalofcomputationalphysics,1987,70（1）：182-192.[4] Hortmann M,PerićM,ScheuererG.Finite volumemultigrid predic⁃tion of laminar natural convection：Bench-mark solutions[J].Inter⁃national Journal for Numerical Methods in Fluids,1990,11（2）：189-207.[5] DemirdžićI,PerićM.Finite volumemethod for prediction of fluid flow in arbitrarily shaped domainswithmoving boundaries[J].Inter⁃national Journal for Numerical Methods in Fluids,1990,10（7）：771-790.[6] Brandt A.Multi-level adaptive technique（MLAT）for fast numeri⁃cal solution to boundary value problems[C].Proceedings of the Third International Conference on NumericalMethods in Fluid Me⁃chanics,Springer Berlin/Heidelberg,1973：82-89.[7] Wang J,Li JF,ChengW X,et parison of finite difference and finite volume method for numerical simulation of the incom⁃pressible viscous driven cavityflow[J].Advanced Materials Re⁃search,2013（732-733）：413-416.[8] Li JF,Long J,Yuan L A,parison of finite difference and finite volumemethod for numerical simulation of driven cavity flow based on MAC[C].Computational and Information Sciences（IC⁃CIS）,2013 Fifth International Conference on,Shiyang,2013：891-894.。
有限容积法有限容积法（Finite Volume Method）又称为控制体积法。
其基本思路是：将计算区域划分为一系列不重复的控制体积，并使每个网格点周围有一个控制体积；将待解的微分方程对每一个控制体积积分，便得出一组离散方程。
其中的未知数是网格点上的因变量的数值。
为了求出控制体积的积分，必须假定值在网格点之间的变化规律，即假设值的分段的分布的分布剖面。
从积分区域的选取方法看来，有限体积法属于加权剩余法中的子区域法；从未知解的近似方法看来，有限体积法属于采用局部近似的离散方法。
简言之，子区域法属于有限体积发的基本方法。
有限体积法的基本思路易于理解，并能得出直接的物理解释。
离散方程的物理意义，就是因变量在有限大小的控制体积中的守恒原理，如同微分方程表示因变量在无限小的控制体积中的守恒原理一样。
限体积法得出的离散方程，要求因变量的积分守恒对任意一组控制体积都得到满足，对整个计算区域，自然也得到满足。
这是有限体积法吸引人的优点。
有一些离散方法，例如有限差分法，仅当网格极其细密时，离散方程才满足积分守恒；而有限体积法即使在粗网格情况下，也显示出准确的积分守恒。
就离散方法而言，有限体积法可视作有限单元法和有限差分法的中间物。
有限单元法必须假定值在网格点之间的变化规律（既插值函数），并将其作为近似解。
有限差分法只考虑网格点上的数值而不考虑值在网格点之间如何变化。
有限体积法只寻求的结点值，这与有限差分法相类似；但有限体积法在寻求控制体积的积分时，必须假定值在网格点之间的分布，这又与有限单元法相类似。
在有限体积法中，插值函数只用于计算控制体积的积分，得出离散方程之后，便可忘掉插值函数；如果需要的话，可以对微分方程中不同的项采取不同的插值函数。
有限容积法（FVM）是计算流体力学（CFD）和计算传热学（NHT）中应用最广泛的数值离散方法。
它通常包括如下五个部分：1.? ? ? ? 网格生成2.? ? ? ? 对流项的离散化3.? ? ? ? 边界条件的离散化4.? ? ? ? 压力速度耦合5.? ? ? ? 离散方程的求解对以上五个部分的处理将直接影响到最准结果的SIMPLE算法自1972年问世以来在世界各国计算流体力学及计算传热学界得到了广泛的应用，这种算法提出不久很快就成为计算不可压流场的主要方法，随后这一算法以及其后的各种改进方案成功的推广到可压缩流场计算中，已成为一种可以计算任何流速的流动的数值方法。
利用有限体积算法三阶迎风型QUICK 离散格式求解二维不可压缩黏性流体方腔流动问题1.二维不可压缩黏性流体方腔流动问题二维不可压缩黏性流体方腔流动（cavity flow ）：有一正方形腔室，其量纲为一的宽度为1.0，里面充满静止的不可压缩黏性流体，方腔内初始时刻压力和密度为=1.0=1.0p ρ、它周围壁面（左右壁面和底面）固定不动，上壁面以量纲为一的速度1.0沿着上壁面方向自左向右运动（图F.1）。
2. 基本方程组、初始条件和边界条件设流体是黏性流体。
二维方腔流动问题在数学上可以由二维不可压缩黏性流动N - S()()u v x y x x φρφρφμ∂∂∂∂⎛⎫+= ⎪∂∂∂∂⎝⎭⎝⎭） 其中u 为变量φ在水平x 方向的流速，v 为φ在垂直y 方向的流速，μ为黏度，S φ为源项。
源项中不仅包含压力梯度项，也包含时间导数项。
初始条件：方腔上壁面以量纲为一的速度1.0沿着上壁面方向自左向右运动。
边界条件：流动速度u v 、均可采用无滑移边界条件，压力p 采用自由输出边界条件。
3.计算网格划分和控制体单元与节点定义采用交错网格，图F.2和图F.3是计算网格、控制体单元和节点示意图。
图F.1二维不可压缩黏性方腔流问题示意图节点P 所在主控制单元如图F.2中有阴影部分所示。
在x 方向与节点P 相邻的节点为W 和E ，在y 方向与节点P 相邻的节点为S 和N ，主控制单元界面分别为w s e n 、、、。
压力p 和速度u v 、分别在三套不同网格中如图F.3中有阴影部分所示。
4.有限体积算法三阶迎风型QUICK 离散格式对方程（F.1）在图F.2所示节点P 所在控制体单元内积分，有：()()V V V V V u dV v dV dV dV S dV x y x x y y φφφρφρφμμ∆∆∆∆∆⎛⎫∂∂∂∂∂∂⎛⎫+=++ ⎪ ⎪∂∂∂∂∂∂⎝⎭⎝⎭⎰⎰⎰⎰⎰（F.2） 由于二维不可压缩黏性流体方腔流动是二维问题，因此控制体单元体积V ∆仅是面积，而它的边界是长度。
设 ,w e s n A A y A A x ==∆==∆，利用Gauss 定理，可将方程（F.2）改写成如下有限体积算法离散格式：()()()()e w n s e w n s u A u A v A v A A A A A S x yx x y y φρφρφρφρφφφφφμμμμ⎡⎤⎡⎤-+-⎣⎦⎣⎦⎛⎫⎛⎫∂∂∂∂⎛⎫⎛⎫=-+-+∆∆ ⎪ ⎪ ⎪ ⎪∂∂∂∂⎝⎭⎝⎭⎝⎭⎝⎭ （F.3） 对上式中e w n sx x y y φφφφ⎛⎫⎛⎫∂∂∂∂⎛⎫⎛⎫ ⎪ ⎪ ⎪ ⎪∂∂∂∂⎝⎭⎝⎭⎝⎭⎝⎭、、、采用一阶向前差分近似，则有： 图F.2方腔流动计算网格、 控制体单元和节点示意图图F.3计算采用的交错网格示意图,,P W E P e wN P P Sn sx x x x y y y y φφφφφφφφφφφφ--∂∂⎛⎫⎛⎫== ⎪ ⎪∂∆∂∆⎝⎭⎝⎭⎛⎫⎛⎫--∂∂== ⎪ ⎪∂∆∂∆⎝⎭⎝⎭ （F.4）同时记：()()()(),,e e w we w n n s sn s F u A F u A F v A F v A ρρρρ==== （F.5）,,,e e w w e w PE PWn n s sn s PN PSA AD D x x A A D D y y μμδδμμδδ====（F.6）则可由式（F.2）写成：()()()()e e w w n n s s e E P w P W n N P s P S F F F F D D D D S x yφφφφφφφφφφφφφ-+-=---+---+∆∆ （F.7）式中P E W N S e w n s D D D D φφφφφ、、、、、、、、都是控制体单元内节点上的已知量，如果利用差分计算得到控制体单元边界上的流通量e e w w n n s s F F F F φφφφ、、、，就可以求出节点上未知量P E W N S φφφφφ、、、、。
为了便于讨论，现对一维对流扩散方程的三阶迎风型QUICK 离散格式进行分析：在三阶迎风型QUICK 离散格式中，计算主控制单元界面上流动量φ需要取主控制单元界面两侧3个节点处的流动量值进行插值计算得到，其中两个节点位于界面紧邻的两侧，第三个节点位于迎风一侧较远邻点，如图F.4所示。
图F.4三阶迎风型QUICK 离散格式示意图当0,0e w u u >>时，通过WW 、W 和P 三个节点值拟合曲线来计算主控制单元左侧界面参数w φ。
通过节点W 、P 和E 三个节点值拟合曲线来计算主控制单元右侧界面参数e φ。
当0,0e w u u <<，则分别通过节点W 、P 、E 和P 、E 、EE 三个节点值计算主控制单元左、右两侧界面参数w φ和e φ。
根据上述计算原则，可以得到界面参数w φ计算公式如下：当0w u >时，界面参数w φ计算公式为：636888ww wW wP wWW φφφφ=+- （F.8a ）当0e u >时，界面参数w φ计算公式为：636888e P E W φφφφ=+- （F.8b ）对于一维无源项一维对流扩散方程三阶迎风型QUICK 离散格式： 当0,0e w u u >>时，三阶迎风型QUICK 离散格式为：P P W W E E WW WW a a a a φφφφ=++ （F.9） 其中()31883818W w w eE e eWW wP W E WW e w a D F F a D F a F a a a a F F =++=-=-=+++- （F.9a ） 同理，若0,0e w u u <<，三阶迎风型QUICK 离散格式为：P P W W E E EE EE a a a a φφφφ=++ （F.10） 其中()3,861,8818W w w E e e w WW eP W E EE e w a D F a D F F a F a a a a F F =+=--==+++- （F.10a ） 将两种流动方向离散方程（F.9）和（F.10）合并后，可得到统一的一维对流扩散方程三阶迎风型QUICK 离散格式：P P W W E E EE EE EE EE a a a a a φφφφφ=+++（F.11） 其中()()()()()61318881836111888118W w w w e e w wWW w wE e e e e e w w EE e eP W E WW EE e w a D F F F a F a D F F F a F a a a a a F F αααααααα=+++-=-=-----=-=++++- （F.11a ）式中01000100w w w w e e e e F F F F αααα>=<=>=<=当时，；当时，；当时，；当时，。
（F.11b ）同理，可以得到带有源项的二维对流扩散方程三阶迎风型QUICK 离散格式为：P P W W E E EE EE EE EE S S N N SS SS NN NN a a a a a a a a a S Vφφφφφφφφφφ=++++++++∆ （F.12）其中S φ为有限体积算法中源项平均值。
式中各个系数为：()()()61318881836111888W w w w e e w wWW w wE e e e e e w wa D F F F a F a D F F F ααααααα=+++-=-=----- ()()()()()11861318881836111888118EE e e S s s s n n s sSS s s N n n n n n s sNN n na F a D F F F a F a D F F F a F ααααααααα=-=+++-=-=-----=-()P W E WW EE S N SS NN e w n s a a a a a a a a a F F F F =++++++++-+- （F.12a ）式中0100,,,k k k k F F k w e s nαα>=<==当时，；当时，。
（F.12b ）源项S φ为：()u pS t xφρ∂∂=--∂∂ （F.13） 若把()nu ρ表示n t 时刻动量，()1n u ρ+表示1n t +时刻动量，则可以得到源项S φ离散格式为：()()()1n nPPe w Vu u S dV x y p p y tφρρ+∆-⋅=-∆∆--∆∆⎰ （F.14）最后，得到有限体积算法二维对流扩散方程三阶迎风型QUICK 离散格式：()()()1n n n n nP P W W E E S S N N n nnn PPewa u a u a u a u a u u u x y p pytρρ+=+++--∆∆--∆∆ （F.15）式中系数k a 为一阶迎风格式中各对应系数。
5.计算结果分析利用三阶迎风型QUICK 离散格式和相应的初始条件和边界条件，求解二维不可压缩黏性流体方腔流动问题。
图F.5是不同雷诺数Re 条件下采用三阶迎风型QUICK 离散格式得到的二维不可压缩黏性流体方腔流动的计算结果。
计算结果和文献中其他高精度算法得到的计算结果进行了比较，两者计算结果十分吻合，能把方腔下壁面两个底角附近二次小涡清晰地计算出来。
这表明有限体积算法三阶迎风型QUICK 离散格式具有相当高的计算精度。
图F.5不同雷诺数Re 条件下采用三阶迎风型QUICK离散格式计算二维不可压缩黏性方腔流动的计算结果Re =100 Re =1 000Re =5 000 Re =10 000从图F.5中可以看出：二维不可压缩黏性流体方腔流动的中心大涡并不在中心位置，方腔内流动也并不对称。
这是因为，方腔上壁面以量纲为一的速度1.0沿着上壁面方向自左向右运动时，在方腔上壁面两侧的两个顶角处不再满足边界条件，这是一个带奇性的方腔流动。
计算结果表明，方腔流动和雷诺数有关，随着雷诺数Re增加，计算精度在降低。
当雷诺数Re较低时，方腔下壁面的两个底角附近的二次小涡十分清晰，随着雷诺数Re的增加，二次小涡变得越来越模糊。
由于三阶迎风型QUICK离散格式计算精度较高，因此三阶迎风型QUICK离散格式计算效果一般要比一阶迎风型离散格式相对来说好些。
using System;using System.Collections.Generic;using ponentModel;using System.Data;using System.Drawing;using System.Linq;using System.Text;using System.IO;using System.Windows.Forms;using ;using System.Collections;namespace FVM_计算方腔流场{public partial class Form1 : Form{public Form1(){InitializeComponent();}public int MX= 30;public int MY =30;public double Re =100.00;public double Ha = 10.00;public double Ri = 1.0;public double Pr = 0.733;public double dt =0.0005;public double c2 =2.25;public double[,] u, v, p, t,uo, vo, po,to, un, vn, pn,tn,streamline;//public double dx, dy;FileStream fileResult, file1,f4;StreamWriter fileResultWriter, fileW1,fw4;public double dx, dy, err, value,x,y;private void Form1_Load(object sender, EventArgs e){u = new double[MX + 1, MY + 2];v = new double[MX + 2, MY + 1];p = new double[MX + 2, MY + 2];t = new double[MX + 2, MY + 1];uo = new double[MX + 1, MY + 2];vo = new double[MX + 2, MY + 1];po = new double[MX + 2, MY + 2];to = new double[MX + 2, MY + 1];un = new double[MX + 1, MY + 2];vn = new double[MX + 2, MY + 1];pn = new double[MX + 2, MY + 2];tn = new double[MX + 2, MY + 1];streamline = new double[MX + 1, MY + 1];dx = 1.0 / MX;dy = 1.0 / MY;fileResult = new FileStream(" X_Matrix.txt", FileMode.OpenOrCreate); fileResultWriter = new StreamWriter(fileResult);file1 = new FileStream(" file1.txt", FileMode.OpenOrCreate);fileW1 = new StreamWriter(file1);f4 = new FileStream(" f4.txt", FileMode.OpenOrCreate);fw4 = new StreamWriter(f4);//double dx1, dy1,//init(u, v, p, dx, dy);//MessageBox.Show(u[5, MY + 1].ToString());}public double max(double a, double b){if (a < b)return b;elsereturn a;}public double alfa(double x){if (x >= 0)return 1.0;elsereturn 0.0;}void init(double[,] u, double[,] v, double[,] p, double[,] t, double dx, double dy) {//u = new double[MX + 1, MY + 2];//v = new double[MX + 2, MY + 1];//p = new double[MX + 2, MY + 2];int i,j;//dx=1.0/MX;//dy=1.0/MY;for(i=0;i<=MX;i++){for(j=0;j<=MY+1;j++){u[i,j]=0.0;if(j==MY+1) u[i,j]=4.0/3.0;if(j==MY) u[i,j]=2.0/3.0;}}for(i=0;i<=MX+1;i++)for(j=0;j<=MY;j++)v[i,j]=0.0;for(i=0;i<=MX+1;i++)for(j=0;j<=MY+1;j++)p[i,j]=1.0;for (i = 0; i <= MX+1; i++){for (j = 0; j <= MY; j++){t[i, j] = 0.0;if (i == MX + 1) t[i, j] = 4.0 / 3.0;if (i == MX) t[i, j] = 2.0 / 3.0;}}}//// 一阶迎风型离散格式// 二维的三阶迎风型离散格式为9点格式，因此有两层边界网格需要//处理，本程序采用一阶迎风型离散格式处理内层，用物理边界条件处理外层。
