《离心叶轮内流数值计算基础》根据作者多年来在叶轮机械与流体力学相关领域的积累和研究成果提炼而成。主要内容包括流体基本属性、基本方程组的推导、网格生成的代数法与微分法、网格量的计算、模型方程的分类及求解特征、差分及其稳定性分析、有限体积法的基本原理、不可压缩N-S 方程的离散计算、边界条件的实施、代数方程系统的迭代法、动-静子耦合流动模型与算法，以及并行编程基础等。《离心叶轮内流数值计算基础》注重理论体系的完整、系统和实用性，将抽象的理论与具体实例相结合、数理基础与当前热点相结合，强调研究思路与解决方法的贯通，既可作为教学用书，也可供科研参考。

更多科学出版社服务，请扫码获取。

　　第1 章流体基本概念
　　1.1 概况
　　一般而言，计算流体力学属于流体力学的一个分支领域，也习惯将其称为计算流体动力学(computational.uiddynamics)，因其所涉及的主要是流体运动问题的数值求解。计算流体力学实质上是在流体力学的基础上，通过与计算技术的融合，拓宽了我们认识新的流动现象的深度和广度。
　　从另一个视角，计算流体力学也属于偏微分方程理论和数值计算的范畴，由这些严格的数学分析和计算技术作为基础，才有可能顺利地开展流动问题的求解，在这个过程中，自然也离不开如偏微分方程的理论分析、数值算法以及程序设计等具体内容。
　　因此，首先从了解流体的物理背景开始，理解并熟悉流体运动的描述方法，熟练掌握流动控制方程的推导，由此建立正确运用流动方程并灵活进行数值计算的必需基础。
　　1.2 流体属性
　　热力学上，通常将物质分为固体、液体和气体三态，很容易通过观察由其形态判别和相互区分。在流体力学中，只有两类物质，即流体和非流体(固体)。固体能够承受外部施加的剪切力并保持静止，而流体是不能的。但这也不是完全清晰的区分，例如，在室温下的一桶沥青，看起来硬得像石头，在上面放一块砖头，也丝毫不受影响。但是，若是把砖头一直放着，则过几天沉到桶底，再要取出来就困难了。因此，沥青通常划归流体。再看金属铝，在室温下，铝是固态的，可以做成任何形状，并且只要不超过材料极限它可以承受任意的外加剪应力。然而，一旦加热以后，铝就呈现液态，在外加剪应力作用下铝液作连续的变形。但是，也不能将高温作为区分金属流体特征的标准，因为金属铅在室温下就呈现出类似的黏性蠕动流特征。同时，我们还可以看到，水银也是一种流体，而且在常见的物质里面，相对于其密度而言，其黏度(运动黏度)是最小的。
　　这里所探讨的是公认为流体的介质，例如，所有的气体都可以认为是真正的流体，还包括一些常见液体，如水、油、酒精等。
　　在真正流体的范畴内，我们来定义和解释它们的属性。大致有如下四类：
　　(1)运动属性(线速度、角速度、涡、加速度以及应变率)，严格地讲，这些是流场属性，而非流体属性；


　　(2)传输属性(黏度、导热系数、质量扩散度)；


　　(3)热力学属性(压力、密度、温度、焓、熵、比热容、普朗特数、体积模量、热膨胀系数)；




　　(4)其他各种属性(表面张力、汽化压力、涡扩散系数、适应系数)。
　　第四项中的某些不属于真正的属性，但取决于流动条件、表面条件以及流体中的所含杂质。使用第三项中的属性也需要注意场合。严格意义上讲，处于运动中的黏性流体并不处于热平衡状态，这些属性并不适用。所幸的是，其偏离平衡态并不远，除了流体的弛豫时间很短且分子数较少，如稀薄气体超音速流。其原理可以这样理解，例如，从统计意义上，常压下的气体也是非常浓的，以1μm边长立方体内气体为例，其中大约含有108 个分子，对于这一体积的气体，当状态发生变化，甚至在发生激波变化时，它也能迅速恢复平衡态。这是因为如此多的分子在短距离内将发生大量的分子碰撞，碰撞后很快又恢复平衡，而液体更浓，所以完全可以假设它们是处于热平衡态的。
　　1.2.1 运动属性
　　对于流体，首先涉及的是流体的速度。而固体力学中，首先涉及的是质点位移，这是因为固体中，质点间以相对刚性的方式相联系。
　　由刚体动力学考虑火箭运动轨迹，只需求解出任意三个不共线的质点运动轨迹就可以了，因为任意其他质点的运动轨迹都可以由已知三点轨迹推算而得出。这种追溯单个质点轨迹的方法称为拉格朗日运动描述，在固体力学中非常有用。
　　但是，考察火箭喷管排出的流体流动时，显然不可能追溯上百万的质点轨迹。这里考察方式就很重要了，地面上的观察者看到的是复杂的非定常流，而与火箭一同飞行的观察者看到的是几乎非常规则的定常流动。因此，在流体力学中采用如下措施通常是非常有用的：①选择最便利的坐标系，使得观察时的流动是定常的；
　　② 只研究作为位置和时间函数的速度场，而不追溯任何质点轨迹。这种描述流动在每个固定点随时间变化的函数方式，称为欧拉方法。欧拉速度向量场可定义为如下直角坐标形式：
　　V(r，t)=V(x，y，z，t)=iu(x，y，z，t)+jv(x，y，z，t)+kw(x，y，z，t)(1-1)
　　求解三个随(x，y，z，t)变化的标量函数u、v、w，通常是流体力学的主要任务。与固体力学采用位移分量描述不同，流体通常采用(u，v，w)三个速度分量来表征。相对而言，位移在流体中没有太多实际用途，因此，也很少给出位移的变量描述。
　　欧拉方法，或者速度场无疑是描述流体力学问题的合适选择，但又存在一个矛盾，即三个基本物理守恒定律：质量守恒、动量守恒和能量守恒，它们都以质点为研究对象，即它们是拉格朗日属性的。这三个定律都与一个质点的某种属性随时间的变化率有关。令Q代表流体的任意属性，而dx、dy、dz、dt分别代表这四个变量的任意变化，这样，Q的总微元变化如下：
　　dQ =
　　Q
　　x
　　dx +
　　Q
　　y
　　dy +
　　Q
　　z
　　dz +
　　Q
　　t
　　dt (1-2)
　　因为我们追踪一个质点，那么空间增量如下：
　　dx=udt，dy=vdt，dz=wdt(1-3)代入式(1-2)，可以得到一个质点的Q属性的时间导数为
　　Q
　　+ u
　　Q
　　x
　　+ v
　　Q
　　y
　　+ w
　　Q
　　z
　　(1-4)dQ
　　=
　　dt
　　t
　　其中，dQ/dt的名称可以是物质导数、质点导数，名字不同，意义都是要让读者产生一种感觉，即我们是追踪一个质点的。为加强这种感觉，通常给它一个特殊符号DQ/Dt，主要是便于记忆，而无其他意义。式(1-4)右端后三项称为对流导数，因
　　Q
　　为，当速度为0，或者Q没有空间变化时，这三项就消失了。
　　项称为局部导数。
　　t
　　注意到，式(1-4)可以写为
　　Q
　　DQ +(V · .) Q (1-5)
　　=
　　Dt
　　t
　　其中，. 为梯度算子，展开如下：
　　i
　　x
　　y
　　+ k
　　z
　　(1-6)+ j
　　1.2.2 质点加速度
　　如果Q就是V本身，可以得到第一个运动属性，质点的加速度矢量为
　　V
　　DV DuDvDw+(V · .) V (1-7)= i + j + k=
　　t
　　Dt Dt Dt Dt
　　u
　　t
　　、v
　　t
　　、w
　　注意到，该加速度与u、v、w，以及12个标量导数有关，如
　　空间导数，形如
　　ui
　　xj
　　，这里，i、j代表三个坐标轴方向。
　　，还有
　　t
　　D/Dt 中的对流导数项不幸地遇到数学上的困难，因为这些是关于速度变量的非线性项。于是，存在有限对流加速度的黏性流动方程在数学上是非线性的，并且从解析的角度也令人头疼，例如，即使是稳态层流，也存在非唯一解、叠加原理不适用等问题。需要注意的是，这些非线性项是加速度，不是黏性应力。具有讽刺意味的是，黏性流动的主要障碍竟然是一个无黏项，而假设黏度为常数时，黏性应力本身反而是线性的。
　　假设流体无黏，即无摩擦，同时又无旋，这时，加速度项仍然存在，但情况会发生乐观的转变：
　　(1-8)
　　此时，
　　压力项是线性的。
　　1.2.3 其他运动属性


　　写
　　。我们来
　　y = vy
　　图1-1流体微团的变形
　　分析这里所发生的变化。首先，由于平移运动，四边形的角点由B点移到B. 点。其次，由于旋转，对角线从BD变到B.D.。再次，由于膨胀，微元体变得稍微大一点。最后，由于剪应变，正方形变成菱形。
　　下面从定量角度进一步讨论。这里为帮助理解，我们给出一个简单的函数变化率，具体如下：
　　ui(x+Δx). ui(x)=
　　ui
　　x
　　Δx =
　　ui
　　x
　　dx (1-9)
　　其中，ui表示u的任一坐标分量。
　　下面进入正题，平移由B点的位移量udt和vdt来表示，平移速率分别为u和v。三维情形，分别为u、v和w。
　　对角线BD的旋转率为dΩz=φ+dα. 45.。注意到，2φ+dα+dβ=90.，代入前式消去φ，可得
　　1
　　dΩz=(dα. dβ)(1-10)
　　2
　　其中，下标z表示旋转轴平行于z轴。还可以推断旋转是逆时针的，因为从图中观察dα略大于dβ。下面给出dα和dβ的具体计算：
　　. .
　　. . .
　　.. .
　　. . .
　　.
　　v
　　x
　　dxdt
　　dx +
　　u
　　x
　　dxdt
　　v
　　x
　　dt
　　v
　　v
　　dα = arctan dt dt= arctan ≈ arctan=
　　1+
　　u
　　x
　　x
　　x
　　dt (1-11)
　　. .
　　. . .
　　. . .
　　. . .
　　.
　　u
　　y
　　dydt
　　dy +
　　v
　　y
　　dydt
　　u
　　y
　　dt
　　u
　　u
　　dβ = arctan dt dt= arctan ≈ arctan=
　　1+
　　v
　　y
　　y
　　y
　　dt
　　(1-12)
　　将式(1-11)和式(1-12)代入式(1-10)，可得
　　v
　　u
　　1dΩz/dt=(1-13)
　　x .
　　y
　　2
　　类似地，沿x和y轴的旋转率为
　　1
　　v
　　1
　　u
　　w
　　w
　　dΩx/dt=dΩy/dt=(1-14)
　　2 y . zz . x
　　这刚好是角速度向量dΩ/dt的三个分量。因为系数12 容易令人费解，习惯上常用一个2倍于它的量ω来表示，如下：dΩ
　　ω = 2 (1-15)
　　dt
　　，
　　2
　　这个新的量ω，在流体力学中是有重要物理意义的，称为流体的涡量。将(1-13)和式(1-14)代入式(1-15)，这样，就将速度向量与涡量构建了关联：ω=rotV=. × V(1-16)因此，涡量的散度为0，即divω=. · ω=div(rotV)=0(1-17)因此，纯数学角度的涡向量代表无源场。如果ω=0，那么流动就是无旋的。下面讨论剪应变，通常将它定义为初始两垂线之间夹角的平均增量，如图1-1所示，AB线和BC线为初始两垂线，那么平均角度增量，即剪应变率为
　　1.dαdβ1
　　v
　　u
　　(1-18)+ +
　　εxy=
　　=
　　x
　　y
　　2dt dt 2
　　类似地，另两个剪应变率分量分别为
　　εzx
　　=
　　w
　　+
　　v
　　z
　　u
　　w
　　1 1 (1-19)+
　　εyz
　　=
　　y
　　z
　　x
　　2 2
　　同样，类似固体力学原理，剪应变率是对称的，即εij=εji。第四个也是最后一个微团运动为膨胀或者拉伸应变。同样参考图1-1，沿x方向的拉伸定义为微团水平方向边长的增长率：
　　u
　　x
　　dxdtdx +
　　. dx
　　u
　　εxx
　　dt = dt (1-20)=
　　x
　　dx
　　同样可得另两个分量，这样，所有三个分量为
　　εxx
　　=
　　u
　　x
　　，
　　εyy
　　=
　　v
　　w
　　z
　　(1-21)
　　εzz
　　=
　　，
　　y
　　将应变率(包括拉伸和剪切)集中到一个对称二阶张量，可写为
　　εij
　　=
　　. .
　　.
　　εxxεxyεxz
　　εyxεyyεyz
　　εzxεzyεzz
　　. .
　　.
　　(1-22)
　　尽管各分量大小随所取坐标轴x、y、z而不同，应变率张量与弹性体应力张量和应变率张量类似，遵循对称张量的坐标变换原理。尤其是无论坐标轴怎么取，或
　　……