第五章--有限体积法

第五章有限体积法前面一章我们推求出了在旋转坐标系下的RANS方程组，即**223210221(,)321iijijieijkjkjijjietktikitDkkkkkikjjuxuuuupuxxxxxppkrkkuuuukkctxxxxxxluCtx　　　　　　2221ikittkikjjjiiktjijtuuuCkxxxxxkuuuGxxxkc120.09,1.44,1.92,1.0,1.3kCCC显然，这一组方程在数学上过于复杂，求解其解析解还十分困难。因此，当前对流体流动问题的研究除了采用实验测量、实验模拟观测察的方法外，在计算上主要采取两种措施：一，根据具体问题对方程进行简化，如无粘流、稳态流、不可压缩流等；二，采用数值计算的方法求解方程。然而，即使经过简化，相当多的流动方程依然无法用解析的方法得到理论解。如无粘、无旋、不可压缩流体所满足的拉普拉斯方程2222220xyz是相当简单的偏微分方程，但当求解区域比较复杂时仍无法得到解析解。因此，数值计算方法成了求解流体流动问题的主要方法。5.1求解流体流动问题的常用数值计算方法数值计算是将描述物理现象的偏微分方程在一定的网格系统中离散，用网格节点处的场变量值近似描述微分方程方程中各项所表示的数学关系，按照一定的物理定律或数学原理构造与微分方程相关的离散代数方程组。引入边界条件后求解离散的代数方程组，得到各网格节点处的场变量分布，用这一离散的场变量分布近似代替原微分方程的解析解。目前，求解流体流动的数值计算方法比较多，有有限差分法、有限元法、有限体积法、边界元法、特征线法、谱方法、有限分析法和格子类方法等。每一种数值计算方法都有各自的特点和适用范围，其中通用性比较好、应用比较广泛的是前4种。1、有限差分法有限差分法是一种比较古老的算法，有很长的应用历史和众多的研究人员，曾经是求解各种复杂偏微分方程的最主要的数值计算方法。有限差分法用差商代替微商，用计算区域网格节点值构成差商，近似表示微分方程中的各阶导数。例如：11,,,nnnniiiiininuuuuuuttxx（5.1.1a）1121122,2nnnniiiinnniiiinuuuuuuuuxxxxxxxx（5.1.1b）式(5.1-1)是用一阶向前差分所表示的一阶导数或二阶导数，类似的可以有一阶向后差分，或中央差分等，当然也可以用二阶差分（三点差分）来表示差商。将表示场变量一阶导数和二阶导数的差商近似式代入微分方程，就可以得到关于各网格点处的差分方程。求解这一组代数方程，可得各节点处的场变量数值解。有限差分法形式简单，对任意复杂的偏微分方程都可以写出其对应的差分方程。但有限差分方程的获得只是用差商代替微分方程中的微商（导数），而微分方程中各项的物理意义和微分方程所反映的物理定律（守恒定律）在微分方程中并无体现。因此差分方程只能认为是对微分方程的数学近似，基本没有反映其物理特征。差分方程的结果有可能表现出某些不合理现象。2、有限元法有限元法是20世纪60年代出现的一种数值计算方法。它最初用于固体力学问题的数值计算，如杆系结构，梁系结构，板、壳、体结构的受力和变形问题。20世纪70年代，在英国科学家ZienkiewiczO.C.等人的努力下，将它推广到各类场问题的数值求解，如温度场、电磁场，也包括流场。有限元法离散方程的获得方法主要有直接刚度法、虚功原理推导、泛函变分原理推导或加权余量法推导。直接刚度法是直接从问题的物理定律、物理公司中得到有限元离散方程。它只适用于比较简单的问题，如梁单元受力变形的有限元离散方程。虚功原理一般只适用于推导弹性力学中物体受力和变形问题的计算过程。变分原理是将微分方程求解问题转换为某泛函求极值问题，再对泛函的表达式进行一定的运算得到有限元离散方程。它可以被用于各类场问题有限元离散方程的推导，但是首先要找到与所求解问题的微分方程对应的泛函。这不是一件容易的事情。在许多情况下所要求解的微分方程没有对应的泛函，例如前述流体流动和传热的控制方程组就没有对应的泛函，因此变分原理推导法不能应用。这时一般采用加权余量法推导。加权余量法的思想很简单，设某些物理问题的控制微分方程及其边界条件分别为()0(f在域内）（5.1.2a）()0(g在域边界上）（5.1.2b）为待求函数。首先选定一个试探函数1=niiic（5.1.3）其中，ic为待定常数，i为试探函数项。将试探函数代入式（5.1.2a）和（5.1.2b），一般来讲不可能正好满足，在域Ω内和边界S上会产生误差，即()fR(5.1.4a)()bgR(5.1.4b)式中，R和Rb称为余量。加权余量法的基本思想是在域Ω和边界S上寻找n个线性无关的函数iW(i=1,2,…,n)，使余量R和Rb在加权平均的意义上等于零，即0iRWd（5.1.5a）0biRWd（5.1.5b）这里称为权函数。式（5.1.5a）和（5.1.5b）表明了这样一个思想：尽管本身不能满足微分方程式（5.1.2a）和式（5.1.2b），但是当其余量与许多线性无关的权函数相乘并积分时，这个余量在总体上接近于零，也就是说在积分的意义上满足微分方程式（5.1.2a）和式（5.1.2b）。当n足够大时，就趋近于真解。例如前述一维稳态对流扩散方程式，按照加权余量法的思想，其解函数在有限元网格系统中的每个单元内应满足()()eddduWdxdxdxdx=0(5.1.6)也就是说解函数在积分意义上满足原微分方程。在选定权函数iW后（有限元法常取单元形状函数为权函数），将对流扩散方程经过一定的运算可得有限元离散方程。有限元法的优点是解题能力强，可以比较精确的模拟各种复杂曲线或曲面边界，网格的划分比较随意，可以统一处理多种边界条件，离散方程的形式规范，便于编制通用的程序。因此，它在固体力学方程的数值计算方面取得了巨大的成功。但在应用于流体流动的数值计算方面却遇到了一些困难。原因是按加权余量法推导出的有限元离散方程也只是对原微分方程的数学近似，有限元离散方程中各项还无法给出合理的物理解释，对计算中出现的一些误差还难以改进。因此，有限元法在流体力学中的应用还存在一些问题。5.2有限体积法的基本思想和特点5.2.1通用变量方程在前面建立的连续方程、运动方程、能量方程、湍动能输运方程、湍动耗散率方程中，尽管变量不同，但它们都有相似的形式，若我们用一个通用的变量来表示，以上方程都可以写成统一的形式，即()()()divdivgradStu（5.2.1）以上方程称为通用变量方程或通用输运方程，统一表示各变量在流体输运过程中的守恒关系。5.2.2有限体积法的基本思想有限体积法与有限元法和有限差分法一样，也对求解区域进行离散，将其分割成有限大小的离散网格。每一网格节点按一定的方式形成一个包围该节点的控制容积V，有限体积法的关键步骤就是将控制微分方程在控制容积内进行积分，即式（4.2-2）。()()()VVVVdVdivdVdivgraddVSdVtu（5.2.2）利用奥高公式，将式等号左端第二项（对流项）和等号右端第一项（扩散项）的体积分转换为关于控制容积V表面A上的面积分。()VAdivdVdAana（5.2.3）利用式（5.2-3）可将式（5.2-2）改写成()()()VAAVdVndAngraddASdVtu（5.2.4）这里我们将等号左端第一项中积分和微分的顺序变换了一下，以方便说明其物理意义。这一项表明特征变量的总量在控制容积V内随时间的变化量。而在左端第二项中左端第二项中()nu意为特征变量由于对流流动沿控制容积表面外法线方向n的流动率（流出）。因此方程左端第二项表示在控制容积中由于边界对流引起的的净减少量。等式右端第一项是扩散项的积分。扩散流的正方向应为的负梯度方向。例如热量是沿着负的温度梯度方向传导的。而n为控制容积表面外法线方向，因此()ngrad是向控制容积外的扩散率。所以()ngrad就表示向控制容积内的扩散率。从而等式右端第一项的物理意义为控制容积内特征变量由于边界扩散流动引起的净增加量。积分后控制变量在控制容积内的守恒关系可表述为：+=+由于边界对流引起的流出量由于内源引起的增加量随时间变化量由于边界扩散引起的增加量5.2.3有限体积法的特点（1）有限体积法的出发点是积分形式的控制方程，只一点不同于有限差分法；同时积分方程又表现出了特征变量在控制容积内的守恒特性，这又有别于有限元法。（2）积分方程中每一项都有明确的物理意义，从而使得方程离散时，对各离散项可以给出一定的物理解释。这一点对于流动和传热问题的其他数值计算方法还不能做到。（3）区域离散的节点网格与进行积分的控制容积分立。如图5.2.1所示二维问题离散系统，实心圆点表示节点，实线表示由节点构成的网格，图中阴影面积表示节点P的控制容积。一般来讲各节点有互不重叠的控制容积。从而整个求解域中场变量的守恒可以由各个控制容积中特征变量的守恒来保证。正是由于有限体积法的这些特点，使其成为当前求解流动和传热问题的数值计算中图5.2.1P节点周围控制容积最成功的方法，己经被绝大多数工程流体和传热计算软件采用。5.3稳态对流扩散问题的有限体积法5.3.1一维稳态对流扩散问题的有限体积法计算格式无源项一维稳态对流扩散问题对一般场变量应满足以下控制方程：()()ddddxdxdxu（5.3.1）式中u为在x方向的流动速度。这样流动必须是连续的，即满足连续方程()0ddxu（5.3.2）式中，为流体密度。此时可认为u为已知值。控制容积如图（5.2.2），对控制容积积分得()()VVddddxdxdxdxdxu（5.3.3）()0Vddxdxu（5.3.4）由奥高公式，控制容积内积分的对流扩散方程式可写成()()AAddAdAdxnunA为控制容积边界面积即()()()()ewewddAAAAdxdxuu（5.3.5）同理，连续方程可写为()0AdAnu即()()0ewAAuu（5.3.6）图5.3.1P节点周围控制容积其中控制容积边界处的值e、w、e、w、eddx、wddx采用相邻节点处值的线性插值运算求出：2PEe，2WPw，e2PE，w2WPEPeePEddxxx，P为书写方便，定义两个参数：Fu——通过单位控制容积界面的对流流量；/Dx——单位界面上扩散阻力的倒数（扩导）。从而wwFu，eeFu（5.3.7a），eePEDx（5.3.7b）设weAAA，将式（5.3.7）表示的参数代人式（5.3.5）和式（5.3.6），有()()()()22ewPEWPeEPwPWFFDD（5.3.8a）0ewFF（5.3.8b）式（5.3.8a）可写成按节点场变量排列的形式，即22