计算流体力学在目前工程计算领域应用越来越广泛,其中有限元法是计 算流体力学的主要计算方法之一。国内大量相关领域的科研工作都是依托国 外公司大型商业软件进行计算的,虽然能够得到各种各样问题的计算结果, 但是科研人员无法从核心算法角度来考察计算结果的正确性,也无法根据自 己的需要对这些软件进行二次开发,这限制了国内相关领域基础理论研究的 发展。
作者毕超,北京化工大学一位年轻有为的副教授,中国塑料加工工业协 会专家,曾荣获2013年北京市青年英才荣誉称号;主要从事聚合物混炼成 型加工装备的开发研制,利用计算流体力学有限元方法分析聚合物在混炼、 成型装备内的流动和混合机理。他将多年科研和教学工作内容进行了总结, 完成了《计算流体力学有限元方法及其编程详解》 一书,从基本理论和编 程方法两个方面对有限元方法在计算流体力学和传热学方面的应用进行了详 细的介绍。
该书采用新颖的理论讲解和实例编程相结合的撰写模式,讲述了计算流体 力学有限元方法的基本理论。书中内容难度由浅入深,将计算实例、理论推 导、编程逻辑、程序编写及结果分析有机结合,归纳简化了使用有限元方法求 解计算流体力学和传热学问题的复杂烦琐过程,讲述了理想流体、牛顿流体、 非牛顿流体流动问题的有限元求解方法,特别是还包括了考虑惯性项影响、非 定常流动以及流热耦合等多种复杂非线性问题的求解方法。该书内容丰富、理 论深入、逻辑清晰,有利于读者更加清晰地了解计算流体力学有限元方法的基 本理论,不仅可以作为计算流体力学领域本科生和研究生的课程教材,而且书 中内容还可以为开发具有自主知识产权的大型工程计算软件提供理论基础。该 书是近些年来我国计算流体力学领域为数不多的理论和实践并重的专业性著作。
毕超博士的这本书,在继承与创新、理论与实践、主要问题与科学体系等 诸多方面都取得了较为理想的成果。因此,我很高兴为本书作序,并谨向广大 科技工作者和工程技术人员推荐此书。
有限元方法是计算流体力学的一个重要分支,在工程计算领域中的应用越来越广泛。笔 者从2005年博士在读阶段就开始从事这一领域的研究,目前已经有8年了,且留校任教后 开设的研究生课程也是围绕这个内容。笔者与很多同学一样,在刚开始接触计算流体力学 时,发现虽然计算流体力学的教材很多,但是大多数教材中的理论都很高深,需要费很大工 夫才能理解其中的意思,编程就更是消耗精力了。于是笔者就想写一本能够让读者易懂的有 关计算流体力学及其相关编程的书。经过8年的学习和教学积累,笔者在2008年开始着手 撰写本书,目前已经有将近5年的时间了,终于在2013年年初完成了本书的书稿。
本书内容从最基本的有限元基础知识讲起,难度逐渐加深,每一章都是针对一个计算实 例进行理论讲解和公式推导的,在此基础上,每个实例都配置有十分清晰的程序代码。本书 共分8章,第1章以有限元方法求解常微分方程为例,讲解有限元方法求解微分方程的基础 知识;第2章以理想流体流动为例,介绍有限元方法求解Laplace 方程的方法;第3章讲述 速度-压力有限元法和罚函数有限元法求解牛顿流体Navier-Stocks ( 简记为N-S) 方程组的方 法,为后续章节奠定基础;第4章讲述非牛顿流体问题的求解方法;第5章讲解考虑惯性项 时N-S方程组的求解方法;第6章讲述与时间有关的流体流动问题的求解方法;第7章讲述 与时间有关的热传导问题的求解方法;第8章讲述速度与温度耦合问题的有限元求解方法。
本书采用MATLAB语言编写计算程序,以便于读者阅读。本书可作为本科生或研究生 计算流体力学课程教材,也可作为相关课程的辅导教材。
本书的原始资料来自于笔者参加北京化工大学江波教授开设的“聚合物加工数值计算” 课程时所记的随堂笔记,在此特别感谢江波教授当年上课时对笔者的指导。
中国矿业大学(北京)朱红青教授担任本书主审,他提出的宝贵意见对本书的撰写具 有很大的指导和帮助,在此表示感谢。
本书献给我的爱妻赵昕楠和刚刚出生的女儿安亿。
由于编者水平有限,书中难免存在不足之处,欢迎广大读者对书中内容多提宝贵意见, 可发邮件至笔者信箱bichao@mail.buct.edu.cn 或bichao812@sohu.com。
本书所有章节的配套程序可在诚睿研究室新浪博客( http://blog.sina.com.cn lul 3296676545)中下载,也欢迎广大读者通过诚睿研究室新浪微博( weibo.comlu13296676545) 与笔者沟通交流。
毕 超
2013 年于北京
本章将以简单的常微分方程为例,介绍有限元方法数值求解过程的基本思路,从而让读 者对有限元方法的基本概念和求解过程具有一定了解。本章分为两节,分别为一维一次常微 分方程的有限元数值求解和一维二次常微分方程的有限元数值求解。1.1节中的实例计算相 对简单,希望读者能亲手逐步推导,体会计算过程。1.2节的内容较1.1节复杂,讲述以编 程思路为主,希望读者能把程序逐行输入MATLAB, 体会编程思路。
1.1 一维一次常微分方程的有限元数值求解
1.1.1 方程及精确解
求解如下微分方程:

边界条件为
x=0 时 ,u=0 (1.2)
值得注意的是,微分方程的边界条件主要分为两类:第一类边界条件,又称为本质边界条 件,即已知边界处场量的数值,如式(1.2)所示;第二类边界条件,又称自然边界条件,
即已知边界处场量导数的数值,如 x=0 时, 带有第二类边界条件的微分方程的求 解,我们将在1.2节中介绍。
式(1.1)和式(1.2)构成的微分方程的精确解为
u =x (1.3)
稍后我们会对比有限元方法计算得到的数值解与该精确解的差别。
1.1.2 有限元方法求解
1. 计算区域的离散
(1)单元类型的选择
在进行有限元计算之前,需要将计算区域进行离散化,也就是通常所说的网格划分。本 节研究的是一维问题,使用一维线性单元对计算区域进行等间距离散,离散结果如图1-1所 示。图中n 为结点序号,e为单元序号,共有N 个结点、E 个单元,且N=E+1。

图1-1 一维线性单元离散结果
(2)结点序号
结点序号共分为两类,即单元内部结点序号和总体结点序号。单元内部结点序号就是各个结点在单元内部的序号,与离散时所选单元类型及插值函数构造有关;总体结点序号是结 点在整个离散区域内统一排列的结点序号。同一个结点可能属于多个单元,该结点在不同的 单元内总体结点序号是唯一的,而内部结点序号可能是不同的。
(3)离散数据的存储
选择单元类型后,便可以建立离散数据,供有限元计算使用。离散结果通常包括以下数 据:
单元信息数据JM: 存储单元所包含的结点序号,其行数等于总单元数,列数与单元内 结点个数一致,本节实例取2,第i 行存储内容为第i 个单元所包含的所有结点序号。
结点坐标数据JX: 存储结点的坐标,其行数等于总结点数,列数与所研究问题的维数 一致,本节实例取1,其第i 行存储内容为第i 个结点的坐标数据。
第一类边界条件数据JB1: 存储第一类边界条件的相关数据。其行数等于处于第一类边 界条件的总结点数,列数与边界条件的分类及所研究问题的维数有关,本节实例取2,第一 列为处于第一类边界上的结点序号,第二列为处于第一类边界上的结点对应的边界值。
第二类边界条件数据JB2: 存储第二类边界条件的相关数据。其行数等于处于第二类边 界条件的总单元边数,列数与所研究问题的维数及离散单元类型密切相关,应包括边界单元 序号、单元内边界边序号、边界边的外法线方向余弦(二维、三维问题)及边界上各个结 点的边界值。本节不涉及该类边界条件,相关内容见1.2节。
以图1-1所示一维线性单元离散结果为例取单元数E=5、 见表1-1~表1-3。

2. 插值函数和权函数
网格离散完成后,单元的插值函数也就随之确定。对于一维线性单元,插值函数包括 Φ1和Φ₂,具体表达式为

式中,ξ为局部坐标,取值范围为ξ∈[0,11。以第二个单元为例(图1-2).分析局部坐 标ξ与笛卡儿坐标x 的对应关系:

在已知单元结点处函数值u 的前提下,可通过与插值函数Φ的结合,近似地描述单元 内部各结点处的函数值u°, 即

用向量形式表示为

式 中 ,I 为单元内结点总数。
插值函数具有两个特性,简述如下:
1)插值函数的个数与单元结点的个数一致;
2)单元内,插值函数只在其对应结点局部坐标处取值为1,其余结点处取值为0。 如果需要详细了解上述内容请参阅参考文献[1]。
3. 加权余量方程的建立
设û为方程(1.1)的一个近似解,即有

式中,r 为残差。用任意权函数u* 乘以残差r, 并在方程定义域[0,1]内进行积分:

数值方法求解方程(1.1)的根本出发点就是要使得残差r→0。式(1 . 9)中u* 又称为 权函数。有限元方法中,权函数的选择具有任意性。本书介绍内容为Galerkin 有限元方法, 该方法规定权函数等于插值函数,即

为了简化表达式,将式(1.9)中用u 代替得到

即得到式(1.1)的加权余量方程。
4. 单元有限元方程的建立
将式(1.11)的积分区域转换到单元内,有

写成局部坐标形式:

由式(1.5)求导,有

则

将式(1.7)、式(1. 14)和式(1. 15)代入式(1. 13),有

5. 总体方程的组合
建立单元有限元方程后,进行总体方程的组合。总体方程形式为

式中,K 为 N×N 矩阵;u 为N×1 向量;b 为 N×1 向量;N 为总体结点数。
(1) K 矩阵的组合
就本节研究实例来讲,第i-1,i 和i+1 个单元的系数矩阵K 可写成:

三个系数矩阵中第二项表达式内K 的下标为单元内部结点序号,第三项表达式内K°的下标 为总体结点序号。将每个单元的单元方程系数矩阵K, 按照总体结点序号对应的位置放置 到总体方程的系数矩阵K 中。如图1-1所示,第i 个结点共用于第i-1 和 i 个单元,所以 Ke(i-1)和 K(i) 中均出现下标为 ( ii) 的项,在K 的 ( ii) 位置要进行K(i-1) 和 K“(i)的求和, 即
K:=K(i-1)+K(i) (1.29)
对于编写程序完成上述组合时,需要利用JM数据进行对位求和。对于第i 个单元来说, 进行双层循环求和,第一层循环指标m 等于1到JM的列数(本例为2),第二层循环指标n 也是由1到JM的列数(本例为2),对于一组m 和 n, 查找JM(i,m) 和JM(i,n) 对应 数值,并完成如下累加计算:
K(JM(i,m),JM(i,n))=K(JM(i,m),JM(i,n))+K(m,n) (1.30)
值得注意的是,m 和 n 就是第i 个单元系数矩阵各个元素对应的单元内部结点序号,JM(i, m) 和 JM(i,n) 对应数值即为第i 个单元系数矩阵各个元素对应的总体结点序号。组合 后,总体方程系数矩阵K 可写成:

(2) b 向量的组合
对于第i-1,i 和i+1 个单元,单元右边向量可写成b°:


三个b° 向量的第二项表达式内b° 的下标为单元内部结点序号,第三项表达式内b° 的下标为 总体结点序号。同样,对于第i 个结点共用于i-1 和i 个单元来说,b 向量相应位置也要体 现两个结点的共同作用:
b₁=b:(i-¹)+b:(3) (1.35)
对于编写程序完成上述组合时,需要利用JM 数据进行对位求和。对于第i 个单元来说, 进行循环求和时,循环指标m 等于1到JM 的列数(本例为2),查找 JM(i,m) 对应数 值,并完成如下累加计算:
b(JM(i,m),1)=b(JM(i,m),1)+b(,1) (1.36)
值得注意的是,m 为第i 个单元右边向量各个元素对应的单元内部结点序号,JM(i,m) 对应数值即为第i 个单元右边向量各个元素对应的总体结点序号。组合后得到

原资料见附件
免责声明:
本页面/内容部分素材来源于互联网公 开 信 息,旨在传递更多信息,不代表本平台立场。
版权归原作者或机构所有,如涉及侵权,请通过平台联系我们,我们将在核实后第一时间处理。
本平台对转载内容的真实性、准确性不作任何保证,用户需自行判断并承担使用风险。