数值传热学(课件)-1

《数值传热学(课件)-1》由会员分享，可在线阅读，更多相关《数值传热学(课件)-1（63页珍藏版）》请在金锄头文库上搜索。

1、第11章 求解二维椭圆型流动 与换热问题的通用程序介绍Presented by Professor W Q Tao2003-12-17,西安,目 录 11-1 程序所采用的数值方法及离散方程 的形式 11-2 程序的结构和各模块的功能 11-3 网格系统 11-4 程序编制的技巧 11-5 程序使用方法及主程序的说明,11-1 程序所采用的数值方法及离散方程 的形式,1111 所用算法的主要特点 原始变量法,以 作为流场求解变量;但把 作为一般 变量时也可以用于 方法；采用区域离散方法B，即先定界面，再定节点位置；采用控制容积积分法导出离散方程，离散格式是守恒的；采用交叉网格， 、 、 分别位

2、于三套网格上；, 对流扩散项采用乘方格式，但可方便地改为五种三点格式的任一种格式（中心差分，一阶迎风，混合格式），可以用延迟修正方法来引入更高阶的格式； 源项采用局部线性化方式处理，即 界面上当量扩散系数采用调和平均方式计算； 非稳态问题在时间坐标方向采用全隐格式； 边界条件均按第一类处理，第二、三类边界条件采用附加源项法处理；, 速度与压力耦合关系采用SIMPLER方式处理：在同一层次上要解两个方程：压力Poisson方程及压力修正值Poisson方程. 11.采用迭代式求解方法：用迭代法求解代数方程组(内迭代)；用亚松弛方式迭代处理非线性关系(外迭代)。,12.用 ADI 线迭代辅以块修正

3、求解代数方程.,1112 三类变量的离散方程,1. 一般 变量变量,经整理可得：,以,等代入质量守恒离散,形式，可得,除 b 外， ， 等计算式同 方程。,边界条件均采用“绝热型”与边界相应的系数为零。,质量守恒方程中无非稳态项 ，速度修正仅考虑了压力修正的影响而未顾及密度修正的影响，因而本程序仅适用于不可压流。,3 ， 变量动量方程（以 方程为例）,控制方程：,离散形式：,确定离散方程系数的方法与一般 变量相同,但注意对一般 ，,，对 则为,图111 控制容积,11-1-3 亚松弛的迭代方式为有利于非线性问题迭代的收敛,两个迭代层次之间变量的变化不宜太大,亚松弛处理可以控制这一变化速度.除了

4、 方程以外,其余、 、 及一般 变量的方程均把亚松弛处 理纳入到代数方程求解过程中,即由该代数方程求解而得的结果就是已经经过亚松弛了的结果:,故有：,最后送入求解器中求解的是,但 的亚松弛采用求解出来 新值后再进行 （为了修正速度的需要）：,112 程序的结构及各模块的功能,1121 通用程序的特点,通用程序可以分为两大类:一类是商用软件(commercial codes),一类是科研人员自己开发,自己使用的程序.对于前者由于开发人员与使用 人员不是同一批人, 因而更应考虑软件使用的方 便性.,商用软件特点简介现在全世界已有五十个左右流动与传热问题计算的商用软件,影响比较大的有：PHEONIC

5、S， FIDAP， FLUENT, CFX, STAR-CD, FLOW-3D等,除了FIDAP为有限元方法外,其余的均为有限容积法.这些商用软件的共同特点是: 有灵活的前处理与输入系统,包括输入计算条件及生成网格；, 有完善的后处理系统,使计算结果的图形显示与输出很方便； 有方便的模块接口,用户可以加入自己开发的模块； 有算例库； 有文件系统(包括在线帮助)帮助用户熟悉与操作该软件； 有完备的防错与检测系统。,科研人员自己开发自己使用的通用软件在“方便” 这一点上不及商用软件,但也要相当的适用性.为此一 般采用的方法是:,（1）采用模块式的结构来写程序，一个程序有若干个模块所组成。所谓模块（

6、module）是指一些程序语句的集合。它有输入、输出，能执行一定的功能。对于调用它的人不必知道模块具体内容是什么，只要了解它的输入输出是什么。FORTRAN语言中的子程序即为模块的一种。采用模块结构使程序可读性强，便于维护。,（2）各模块间的联系应松弛，使一个模块的变动对 其它模块的影响最小。,（3）把程序的分成不变部分及用户部分两大块，不变部分在本程序使用范围内不必变动，它只向用户部分要求一定信息，经过加工处理，把结果输出。构成了该程序的黑匣子（blackbox）。变动部分则与具体问题有关。,（4）方程的离散求解应分为两个模块（算法与格式 应分为不同模块），以适应对格式改进及代数方程 求解方

7、法的改进的研究，也适用于比较不同格式计 算结果。,（5）对常用变量应设置预置值（default values）,（6）有一定的前、后处理功能。,本课程所介绍的程序属于这一类，但前后处理模 块不完善。,1122 本程序的总体结构,一个工况的计 算只执行一次,每一个层次的迭代执行一次,主程序,用户子主程序,图112 程序结构,1123 各模块的主要功能 0MAIN （1）规定计算过程总体的流程 （2）判断是否停止执行程序 1GRID网格划分(1)选定坐标系：MODE1，2，3 代表x-y,r-x及r-(2)设置x,y方向计算区域长度：XL,YL(3)规定x,y方向节点总数L1，M1(4)设置x,y

8、方向主控制体界面位置，XU(i),YV(j), i=2, L1, j=2,M1(方法B)。,(5)当MODE1，规定计算区域径向起始半径R(1) (6)网格均匀时，调用UGRID生成界面位置,图113 极坐标系统,2UGRID 根据XL，YL及L1，M1生成均分网格的界面位置。 3SETUP 1（1）设置在计算过程中保持不变的一维几何参数；（2）设置节点坐标，X(i), Y(i), i=1L1, j=1M1（3）生成全控制体宽度：XCV(i), YCV(j),i=2L2, j=1M2（4）确定相邻两节点间的距离：XDIF(i), YDIF(j),i=2L1, j=2M1，XDIF(i)=X(i

9、)X(i-1),YDIF(j)=Y(j)Y(i-1),（4）生成U，V各自控制容积宽度：XCVS(i), i=3L2， YCVS(j), j=3M2 （5）设置Y方向半径R(j), X方向 scaling factor SX(j) （6）计算与X方向相垂直的面积： ARX(j) （7）设置 XCV(i)， XCVIP(i) （8）建立插值函数FX(i), FXM(i)等。,（9） 设置 ， ， ， ， ， 的初值。除 外其余均为0值。,4START（1）对非稳态问题规定初始条件；（2）对稳态问题规定迭代的初场；固定不变的边界条件也可在此引入。以上四个模块在一个工况计算中知执行一次。 5DENS

10、E规定流体的密度场；对常物性问题可不写任何语句，但应保留空块。,6BOUND 设置各变量的边界条件。,7OUTPUT（1）每做一个层次的迭代（代数方程系数变换一 次，即ITER之值加1）输出几个代表性的变量以观察收敛程度；（2）计算一些特殊的二次量，如： ， ， ， 等；（3）调用PRINT，输出二维的场。,8PRINT输出二维的场。,9SETUP2 （主程序之核心模块）（1）调用GAMSOR，以确定 ， ， 等；（2）调用DIFLOW，以确定所用格式的 ；（3）建立离散方程系数，设为 ， ， ， ， 及 ；,（4）调用SOLVE求解代数方程；,（5）迭代次数增值ITERITER+1。,10G

11、AMSOR（1）确定各变量的 值： 、 为 ， 为 ；,（2）确定各变量 ， ，存入 ， ；（3）确定第二，第三类边界条件的边界控制容积的， ， 并累加于 ， 中。,11DIFLOW据 及 之值计算相应的 。 12SOLVE用交替方向线迭代+块修正求解方程，采用迭代 方式，迭代次数由 来控制，一般在1-10 之间。,1124 本程序功能及限制 功能：（1）可以计算三种二维正交坐标的不可压缩流动及传热问题；（2）可同时求解10个变量，打印13个变量；（3）控制方程可为有量纲方程亦可为无量纲方程。 限制:(1)不宜计算非线性强烈的非稳态问题；(2)不宜计算几何形状复杂的问题；(3)不能计算可压缩流

12、。,11-3 网格系统 11-3-1 三种坐标系中的有关规定 1. 直角坐标系,XL,YL,（1）MODE1； （2）Z 方向为单位厚度； （3）坐标原点位于计算区域的左下方。,2. 圆柱轴对称坐标系,Y,X,X,XL,YL,R,R(1),（1）MODE2；（2）计算对 1弧度进行；（3）R(J) 从对称周起算；（4）R（1）应给定。,3. 圆柱轴对称坐标系,（1）MODE3；（2）Z 方向取单位厚度 （3）R(J) 从对圆心起 （4）R（1）应给定；（5） 角应小于2 。,11-3-2 主要几何参数及变量编号方法,1.控制容积界面：XU(i)，i=2,L1, YV(j)，j=2,M12.主节

13、点： 最后三个节点X方向为L1,L2,L3,Y方向为M1,M2,M3,3.主控制体宽度： XCV(i)，i2,L2; YCV(j), j=2,M2,Y,4.节点间距:,5.速度控制体宽度 : XCVS(i), i=3,L2; YCVS(j), j=3,M2,6.速度编号： 箭头所指节点号为速度的编号。,7.代数方程求解区域起始节点编号:,由于代数节点的求解限在内点进行，而离散方程求解区域边界的第一内节点编号对u,v,p是不一样的，因而程序中设置了两个FORTRAN简单变量IST，JST,表示X,Y方向迭代开始的内点编号，其值如下表所示：,变 量,IST,JST,11-3-3 坐标系统综合图示,

14、Y方向的几何图示与此类似 。,1134 关于压力场计算的几点说明：1边界上的压力值是在获得收敛的解后采用外推法而获得的；2与边界相邻的速度控制体中压差与作用面积的 计算：与边界相邻的速度控制体比内部控制体要长，但在该速度的离散方解中，如对u(3,j)，压力差均取为p(2.j)-p(3,j)。因为边界压力未知，采用增加压力作用面积的方法来弥补缩小了的差压，即在u(3,j)中压力差项为：,在内部XCVS(i)=XDIF(i)，因此VOL(i,j)/XDIF(i)，就是压力作用面积。而现在把作用面积扩大了，相当于把P(2,j)-P(3,j)线性外推来获得P(1,j)-P(3,j)的结果。,3计算区域四个角顶在计算中并不用到，为输出时的需要，采用以下插值方法计算之：,不可压缩流体压力仅是相对值，输出时以 (IPREF,JPREF)点的压力为基准，即：,程序中的IPREF，JPREF的预置值均取为1。,114 本程序编制的技巧,