有限体积法基本原理 课件 chap06-2D输运方程及有限体积法

大连理工大学工业装备结构分析优化与CAE软件全国重点实验室有限体积法基本原理夏广庆,张军军,陈冲,

韩亚杰,鹿畅6.12D输运方程通用温度输运的对流扩散方程：忽略时间导数项和对流项，可以得到：在直角坐标系下将梯度算子与点积算子展开，可以得到：6.2控制方程离散控制方程中各项在有限体积内积分：

6.2控制方程离散根据散度定理将控制方程中体积积分项进行变换：单元具有有限数量的面（𝑁），对于四边形单元，𝑁=4，而对于三角形单元，𝑁=3。因此，对整个表面的曲面积分可以用对单元的𝑁个面积分来代替：在二阶有限体积法中，流动量在整个面上都是线性变化的，因此，整个面的积分可以简化为面心的值（常数）：6.2控制方程离散——内部单元离散格式为进一步简化方程，参考图6–1中的2D四边形单元，该单元有左面（𝑙）、右面（𝑟）、顶面（𝑡）和底面（𝑏）。单位法向量总是指向单元外，𝑥方向从左到右为正，𝑦方向从下到上为正。这种离散方法适用于网格中所有单元形状和方向与图6–1中四边形单元（正规网格）相同，为进一步简化，接下来将对内部单元和边界单元的离散格式分别考虑。6.2控制方程离散——内部单元离散格式所有扩散项都将使用中心差法，使用图6–2中的符号表示单元中心点之间的距离𝑑，上式可以表示为：

按照单元中心温度（𝑇𝑃、𝑇𝑅、𝑇𝐿、𝑇𝑇和𝑇𝐵）对上式进行排序，可以得到：6.2控制方程离散——内部单元离散格式所以对于内部单元，2D扩散方程的有限体积离散格式及系数如下：使用求和符号对离散格式进行表示：6.2控制方程离散——边界单元（Dirichlet）右、上和下表面与内部单元相连，因此，通过这些面的热扩散可采用与内部单元相同的处理方法来进行计算：本章将在左面应用定值（Dirichlet）边界条件。左边界面上的温度为𝑇𝑤𝑎𝑙𝑙，到壁面的距离为到下一个单元中心点距离的一半。因此，通过左边界面的热扩散为：将上式代入边界单元对应离散格式中：6.2控制方程离散——边界单元（Dirichlet）为与内部单元离散格式保持一致，将上式变换为：需要注意的是，使用求和符号时对于边界面𝑎𝑁

=0，并且𝑎𝑃的求和不包括边界面。6.2控制方程离散——边界单元（Neumann）边界施加固定通量（Neumann）时的边界单元如图6-4所示，该单元右、上和下表面与内部单元相连，因此，通过这些面的热扩散可采用与内部单元相同的处理方法来进行计算：

6.2控制方程离散——边界单元（Neumann）为与内部单元离散格式保持一致，将上式变换为：需要注意的是，使用求和符号时对于边界面𝑎𝑁

=0，并且𝑎𝑃的求和不包括边界面。6.3控制方程各单元系数与线性系统根据前述推导，可以得到控制方程对应单元各系数(Dirichlet边界与内部单元)：Neumann边界与内部单元系数表：单元类型Interior0Left0Right0Bottom0Top0单元类型Interior0Left00Right00Bottom00Top006.3控制方程各单元系数与线性系统将前述单元系数用求和符号表示，如下表：单元面类型Interior0Dirichlet0Neumann00

6.2控制方程离散——多边界界单元对于2D情况，一些边界单元可能具有不止一个边界面，且不同边界面的边界条件也可能不同，如图6-5所示，计算域左上角具有2个边界面，该单元的右面与下面与内部单元相连，因此，通过这些面的热扩散通量可采用与内部单元相同的处理方法来进行计算：：如果边界面的左侧为固定温度𝑇𝑤𝑎𝑙𝑙,𝑙，顶部为固定温度𝑇𝑤𝑎𝑙𝑙,t，那么也可以在边界面上使用中心差分法。6.2控制方程离散——多边界界单元为与内部单元离散格式保持一致，将上式变换为：使用求和符号表示：6.3Example1——带热源的热传导根据问题模型，可以知道杆件中温度满足如下控制方程：6.3Example1——带热源的热传导Step-1：网格划分Step-2：添加材料属性Step-3：计算方程中各系数单元1010100-4071006021010100-2051005031010100-205100504100100-409100605

0101010-202100506

101010100100407

101010100100408

1001010-204100506.3Example1——带热源的热传导Step-3：计算方程中各系数单元90101010-2021005010

101010100100401110101010010040121001010-2041005013

010010-40510060141010010-2031005015

1010010-2031005016100010-40710060Step-3：计算方程中各系数(采用求和符号表示各系数)单元面类型Interior10100Dirichlet(L)2002000Dirichlet(B)2003000Dirichlet(R)2004000Dirichlet(T)20050006.4Example1——带热源的热传导OpenFOAM求解器构造——求解器结构及各文件功能简介(复习)files文件主要完成两个工作：1）指定需要编译的文件，这里的文件是代码文件，即.C文件，而不包括.H文件；2）指定编译类型及对应类型名称，如编译成库类型，其关键字为LIB=···；编译成可执行程序类型，其关键字为EXE=···。options文件也主要完成两个工作：1）指定编译的头文件所在目录，如果所用头文件在当前文件夹或为标准C++头文件，无需指定，EXE_INC和关键字-I用于指定头文件所在目录；2）指定编译当前程序用到的库，EXE_LIBS和关键字-L/-l用于指定用到的库。6.4Example1——带热源的热传导createFields.H——初始化基本物理场待求温度场速度场温度扩散系数热源OpenFOAM求解器构造OpenFOAM求解器构造6.4Example1——带热源的热传导convectionAndDiffusion.C——求解器源文件网格读取线性方程构造方程求解及输出算例文件路径正确性检测时间类初始化网格类初始化基本物理场初始化构造线性方程：各项与算例问题中对应线性方程求解结果文件输出6.4Example1——带热源的热传导OpenFOAM求解器构造编辑完成OpenFOAM各求解器后，在求解器下执行wmake对求解器进行编译，得到可执行程序。求解器源文件编译得到求解器名称g++之后的指令参数依次指定了采用的C++标准版本、编译器参数、机器位数、整数位数等信息，之后以“-I”开头的表示使用到的头文件所在的目录，以“-L”开头的表示指定的库，以“-1”开头的表示链接到指定的库，“-oapp”表示生成的程序名为app。6.4Example1——带热源的热传导OpenFOAM

求解器对应算例构造——带热源热传导算例构造算例文件较为方便方法修改已有的、与问题模型相似的算例，这里以修改cavity为例进行介绍，对应算例文件路径：openfoam10/tutorials/incompressible/icoFoam/cavity6.4Example1——带热源的热传导OpenFOAM求解器对应算例构造——带热源热传导算例6.4Example1——带热源的热传导OpenFOAM求解器对应算例构造——带热源热传导算例6.4Example1——带热源的热传导OpenFOAM求解器对应算例构造算例文件编辑完成后以此执行：blockMesh、../app、paraFoam分别进行网格划分、算例求解及后处理，对应各步骤过程分别如下：6.4Example1——带热源的热传导OpenFOAM求解器对应算例构造算例文件编辑完成后以此执行：blockMesh、../diffusionEquationSolver、paraFoam

分别进行网格划分、算例求解及后处理，对应各步骤过程分别如下：6.4Example1——带热源的热传导自编程序简介——2D有热源热传导通过自编程序方式进行问题求解，程序主要包含模型参数设置、网格划分、矩阵系数计算、线性系统构造、求解及后处理。程序各部分定义参考如下：6.4Example1——带热源的热传导自编程序简介——2D有热源热传导通过自编程序方式进行问题求解，程序主要包含模型参数设置、网格划分、矩阵系数计算、线性系统构造、求解及后处理。程序各部分定义参考如下：6.4Example1——带热源的热传导算例结果在进行后处理时较为方便的绘图方法是将结果保存在网格节点上而非网格单元中心，如图6–10所示，因为边界上的温度值是分布在网格节点上（图中红色部分），如果采用单元中心（图中蓝色部分）进行后处理则无法得到边界上的温度分布。常见的CFD程序采用单元中心进行计算，因此，内部节点值（图6–10中以白色显示）需要通过单元中心点值（蓝色显示）之间的插值计算得到。6.5热平衡对于2D热平衡检验时，内部单元一般有限体积离散格式如下：根据傅里叶定律通过每个单元面的热通量为：将上式代入对应离散格式中可以得到：负号是为了确保热量流动方向与温度梯度方向相反对于2D四边形单元而言，热平衡误差表示为：6.5热平衡计算前述算例问题中各单元的热平衡误差，结果如下表所示：6.5热平衡单元热通量平衡表还可用于计算平板每个面（跨越所有边界）流出的总热通量：将每个边界面（上表中用红色标出的部分）流出的热通量相加，结果如表6–8所示：正热通量表示热量穿过平板从平板边界流出，负热通量则表示热量通过边界进入平板内。左边界作为温度最低的面（𝑇𝑙

=100℃），该平面的热通量最大。相反，上表面作为最热面（𝑇𝑟

=250℃），热量通过该平面进入平板内，从而产生