高等流体力学第2讲.doc

第二讲 流体运动微分方程
一、应力张量
作用在流体上的力可以分为两类，即质量力和表面力两大类。作用在连续介质表面上的表面力通常用作用在单位面积上的表面力——应力来表示，参见图2-1，即
（2-1）
式中　n为表面积ΔA的外法线方向；ΔP为作用在表面积ΔA上的表面力。pn除了与空间位置和时间有关外，还与作用面的取向有关。因此，有

需要特别指出，应力pn表示的是作用在以n为外法线方向的作用面上应力，其下标n并不表示应力的方向，而是受力面的外法线方向，见图2-1；一般来说，应力pn的方向并不与作用面的外法线n一致，pn除了有n方向的分量pnn外，还有τ方向的分量pnτ。只有当pnτ=0时pn才与n的方向一致；图中ΔA右侧的流体通过ΔA作用在左侧流体上的力为ΔP=pnΔA，而ΔA左侧的流体通过ΔA作用在右侧流体上的力为ΔP=p－nΔA，这两个力互为作用力和反作用力，所以有

可得
pn=-p-n （2-2）
图2-2 一点处的应力状态
z
p-y
M
x
y
B
C
A
p-z
pn
n
p-x
n
图2-1　pn与n的关系
M
n
ΔP
ΔA
－ΔP
－n
n
pn
pnn
pnτ
M
p-n

为了研究一点处微元面积上的表面力，先在流体中以M为顶点取一个微四面体，如图2-2所示。设MA=Δx，MB=Δy，MC=Δz，ΔABC的法向单位矢量为n，则

或简写为
（2-3）
设ΔABC的面积为ΔS，于是ΔMBC、ΔMCA、ΔMAB的面积可分别以ΔSx、ΔSy、ΔSz表示为
（2-4）
四面体的体积可表示为

式中h为M点到ΔABC的距离。根据达朗贝尔原理，可给出四面体受力的平衡方程为

当四面体趋近于M点时，h为一阶小量，ΔS为二阶小量，ΔV为三阶小量，略去高阶小量后可得

再考虑式（2-2）和（2-4）可得
（2-5）
上式在直角坐标系中的投影可表示为

（2-6）

上式也可以用矩阵形式表示为
（2-7）
也可以表示为

式中 P= （2-8）
称为应力张量。
这里需要着重指出的是：应力张量各分量的两个下标中，第一个下标表示的是该应力作用面的法线方向；第二个下标表示的是该应力的投影方向，例如pxy表示它是作用于外法线为x轴正向的面积元上的应力px在y轴上的投影分量。应力张量P描述的是某一点处的应力状态，过该点的任意一个曲面上的应力pn均可由式（2-7）确定。与矢量相似，张量也是客观的，正如矢量确定以后，它的大小和方向不会随着坐标系的改变而改变，所改变的只是在不同坐标系下其分量的大小。
无粘流体或静止流场中，由于不存在切向应力，即pij=0（i≠j），此时有
P===-p= -pI
式中I为单位张量，p为流体静压力。
流体力学中，常将应力张量表示为
（2-9）
式中p为静压力或平均压力，由于其作用方向与应力定义的方向相反，所以取负值；T称为偏应力张量，即
T= （2-10）
偏应力张量的分量与应力张量各分量的关系为：i＝j时，pij为法向应力，ii = pij- p；当i≠j时pij为粘性剪切应力，ij =pij。ii=0的流体称为非弹性流体或纯粘流体，ii≠0的流体称为粘弹性流体。
二、应变张量
u(M)
图2-3一点邻域的速度
M0
M
u(M0)
dr
与刚体相比，连续介质运动过程中还有可能发生变形，因此连续介质的运动比刚体的运动要复杂得多。在这里，首先回顾一下刚体运动速度分解定理。刚体的运动可以分解为随质心的平动和绕质心的转动，即

其中u0为刚体质心的平动速度；u为刚体内部任意一点处的运动速度；ω为刚体绕质心的旋转角速度；dr为质心至某点的微元矢量。
在t时刻的连续介质中取出包括点M0（x，y，z）的任意微元体积，同时取微元体积内的另一点M（x+dx，y+dy，z+dz），如图2-3所示。假设点M0的速度为u（x，y，z），当dr=（dx，dy，dz）为小量时，M点的速度可用M0的速度的泰勒展开式来表示，即
　　　　　　　　 （2-11）
或分量形式

显然，du或(du，dv，dw)是M点相对于M0点的相对运动速度，它可以用矩阵的形式为
　 　　　 （2-12）
上式中的方形矩阵可分解为

= R＋D （2-13）
上式中第一个矩阵R是反对称的，第二个矩阵D是对称的，这两个矩阵在流体力学中也称为二阶张量，下面就来具体分析这两个张量的物理意义。
反对称矩阵R中的九个分量中只有三个独立分量，即
，， （2-14）
这三个分量恰好就是流体微团旋转角速度矢量的三个分量，因此，将R称为旋转张量。同时ω=ω1i+ω2j+ω3k也就是速度矢量的旋度的一半，即
（2-15）
对称矩阵D中的九个分量中只有六个独立分量，
，，，
， （2-16）
Dii (i=x，y，z)恰好是流体力学中研究过的流体微团在三个坐标轴方向上的线应变速率，而Dij (i=x，y，z；j=x，y，z且i≠j)也恰好是其角变形速度。因此，流体力学中将张量D称为应变速率张量，或简称为应变张量，将R＋D称为速度梯度张量，用gradu表示。
在非牛顿流体力学中，也常用一阶Rivlin-Ericksen张量A来表述应变速率的大小，它与D的关系为
A=2D （2-17）
一阶Rivlin-Ericksen张量A的分量直角坐标系中的表达式可由式（2-16和17）得出，其在柱坐标系和球坐标系中的表达式的推导比较复杂，其结果见表2-1。
表2-1 一阶Rivlin-Ericksen张量A的分量在柱坐标系和球坐标系中的表达式
柱坐标系(r，θ，z)
柱坐标系(r，θ，φ)

由矩阵分析可知，对称张量A有三个不变量，即

（2-18）

其中最常用的是第二不变量，用来描述流场的剪切速率的大小，在简单剪切流动中常用来表示。。
例2-1试分析下板不动上板做匀速运动的两个无限大平板间的简单剪切流动
，，
式中k为常数，且k=u0/b。
b
u=ky
u0
x
y
图2-4 简单剪切流动
解：由速度分布和式(2-14、16和17)可得

再由式(2-18)可得

所以II=k=u0/b。

流动的旋转张量R的分量不全为零说明流动是有旋流动，I=trA=0表明流动为不可压缩流动，II=k表明了流场的剪切速率为常数。
三、以应力表示的运动方程
以应力表示的运动方程是将于动量定理改写为适用于控制体的形式后所得到的数学表达式。
图2-6微元体受力示意图

在连续介质中取一个如图2-6所示微元体。假设微元体中心处的应力张量为
P=
则作用在图中右、前和上三个侧面上的应力分量如图所示。作用在微元体上的所有表面力在x方向上的合力为
（在垂直于x 轴的面上的合力）
（在垂直于y 轴的面上的合力）
（在垂直于z 轴的面上的合力）
简化后可得控制体内的流体所受的表面力的合力在x方向的分量为

作用在微元体上的动量变化率在x方向上的分量为

代入动量定理后两端同除以微元体内的流体的质量ρdxdydz后，可得
（2-23a）
同理，利用y、z两个方向上的动量定理可得
（2-23b）
（2-23c）
式中的Ｘ，Ｙ和Ｚ分别是作用在单位质量流体上的质量力在x，y和z三个坐标轴方向上的分量。这便是以应力表示的粘性流体运动方程，也称为柯西应力方程，其矢量形式为
　　　　　 （2-24）
式中，为拉普拉斯算符与应力张量的内积。式(2-23)是柯西应力方程在直角坐标系中的表达式，它们在柱坐标系中的表达式见表2-2。
表2-2 柱坐标系中的连续方程和运动方程
连续方程　
动力学方程
r分量：
θ分量：
z分量：
注：表2-2的公式中的fr，fθ，和fｚ分别是作用在单位质量流体上的质量力在r，θ和z三个坐标轴方向上的分量。
柯西应力方程对任何粘性流体、任何运动状态都适用。很容易看出，柯西应力方程只包含三个方程，加上一个连续性方程也不过4个方程，而其中未知变量却有9个(6个应力分量和三个速度分量)，所以，这组方程是不封闭的。为使该方程组在理论上可解，必须进一步考虑应力和应变率之间的关系——本构方程，作为补充方程以获得封闭的方程组。

四、本构方程
流变学中对本构方程的定义是：在某些假定条件下，材料力学行为的数学描述，著名的牛顿内摩擦定律、虎克定律便是最原始形式的流体及弹性体的本构方程。在非牛顿流体力学中，本构方程特指应力张量与应变张量之间的关系方程，非牛顿流体不同于牛顿流体，它没有一个单一的本构方程，因为非牛顿流体类型繁多，特性各异，不能指望用一个本构方程来描述，下面就简单介绍几种常用的本构方程。
（一）纯粘流体的本构方程
在流体力学中的本构方程，是在某些假设下，连续介质的力学行为的数学描述，如弹性力学中的虎克定律和流体力学中的牛顿内摩擦定律等。在非牛顿流体力学中，本构方程给出的是偏应力张量和应变张量之间的函数关系，即间接地给出了偏应力张量与速度之间的关系，本构方程、连续性方程和以应力表示的运动方程便可构成封闭的方程组。本构方程作为流变学的重要研究内容，所涉及的内容浩繁复杂。在这一小节中主要着重介绍石油工程中常用的几种本构方程。
1.平均压力的引入
将流体静力学中的静压力引入到应力张量后，得到了以式(2－10)表示的应力张量，其中的偏应力张量可表示为
P= （2-25）
严格意义上讲，只有在静止流体和理想流体中才存在静压力p，在实际流体的流场中不存在静压力，这是因为在实际流体的流场中pxx、pyy和pzz会随着坐标系方向的不同而变化。幸运的是，对于各项同性的流体，应力张量的第一不变量为常数，讲这个不变量定义为平均压力，即

因此，可以把平均压力看作是静压力。引入平均压力后，式（2-21）可写做
　　　　　 （2-26）
至此，连续性方程（2-22）和运动方程（2-26）中，一共有u、p、T三个未知数，如果再加上本构方程T=f（u）后便可构成一个封闭的方程组。
2.纯粘流体的本构方程
假设或直接说希望存在一种流体，它满足
　　　 　 （2-27）
或以分量的形式表示为
　　　 　 （2-28）
图2－7库特流
x
o
y
u(t,y)

b
u0
式中μ是一标量函数。上式所定义的一类流体确定存在，比如牛顿流体。
下面仍以图2－7所示的稳定库特流为例来研究对象来研究粘性流体的本构方程，流场的边界条件为

将以式（2-28）表示的应力分量Tij代入柯西应力方程（2-23），并利用连续方程可得到
（2-29）
（2-30）
式（2-30）表明p不是y的函数。再假设体系在x方向压力梯度为0，则可得到
（2-31）
解之可得
（2-32）
由式（2-35、15和17）可知A的非零分量只有：
（2-33）
因为应力ij正比于Aij对应的分量，所以由式（2-28）可以得到，应力张量T的非零分量为
（2-34）
这就是牛顿内摩擦定律在初等流体力学中一般采用的形式。比较严格的定义是通过方程（2-27），其中使μ等于常数。
现在再回到方程式（2-27），μ是一个标量系数的情形。特别是这意味着η可能是某个函数，而不像牛顿流体那样μ为常数。在这种情况下，我们把μ称为非牛顿视粘度或粘度系数。
对于其它流体而言，其视粘度μ究竟是哪些变量的函数呢？这里不准备去验证某些想法，而是假设μ只依赖于由应变张量A的分量定义的流变学量。由于μ是一个标量，而A的却是一个张量，所以标量μ不能直接表示成张量A的函数，若想在两者之间建立起函数关系，必然要先引入一种变换或者运算将张量A的转换成标量。
首先讨论一个比较熟悉的问题，质点的动能依赖于它的速度矢量的大小。但动能是一个标量，而速度却是矢量。这个问题的解决方法在于，存在着一个矢量的标量函数，即
（2-35）
速度的平方u2是一个标量，它与坐标系的选择无关，这类与坐标系无关的矢量函数称为不变量。任何矢量V都可以组成唯一的标量函数V·V。由此可以推广到，任何张量V也都可以组成唯一的标量函数V·V。由此不禁使我们想到了应变张量的三个不变量，再回想前面已经分析过的平板间的简单剪切流动，其三个不变量中只有第二不变量II不为零，由此不难想到，对于不满足式（2-27）的非牛顿流体，其视粘度函数必然是应变张量A的第二不变量II的函数，即
（2-36）
因此，对于满足式（2-27）的非牛顿流体简单剪切流动而言，偏应力张量中不为0的两个分量为：
（2-37）
由于满足式（2-7）或（2-28）的流体的法向应力ii =0，称之为纯粘流体或广义牛顿流体。
（二）几种常用的本构模型
前面讨论的纯粘流体或广义牛顿流体的视粘度为常数，其本构关系比较简单，即粘性切应力与应变张量成正比例关系。但是，石油工程中常见的非牛顿流体的本构方程比纯粘流体复杂，大致可划分为塑性型、假塑型、屈服假塑型和膨胀型等不同的类型，下面就来探讨一下这几种类型流体的本构方程。
1.宾汉模型
尤金·库克·宾汉，美国化学家。曾任拉法耶特学院化学系教授与系主任。宾汉和马尔克斯·雷纳一起开创了流变学这一学科，并于1929年创造了“rheology”（流变学）这一词汇。宾汉在流变学的理论研究和实际应用方面均有较多贡献，最著名的是提出宾汉模型，在该模型中，当应力小于某一临界值时，流体呈现弹性体特征，应变速率为零。当应力大于临界值时，流体呈现为牛顿流体。符合这一描述的流体称为宾汉流体。1948年宾汉逝世后，流变学协会每年颁发宾汉奖章，为流变学领域的最高奖赏之一。
由图2-8所示不同类型非牛顿流体的流变曲线中，单独将塑性流体流变曲线拿剪切速率
图2-8塑性流体的流变特性曲线
粘性切应力τ
出来分析，由图2-8所示的流变曲线可以看出，当所流体受到的切应力低于某一数值时不会流动，而当所受到的切应力超过这一数值时才开始流动，这一数值称为塑性流体的极限静切应力；当应变速率较小时，粘性切应力与剪切速率成曲线关系；当应变速率较大时，粘性切应力与剪切速率成线性关系。流变曲线的直线段的斜率，称为塑性粘度或结构粘度，该直线段的延长线与切应力轴的交点，称为屈服应力或极限动切应力。因此，塑性流体符合宾汉模型，其本构方程可以用宾汉公式表示为
（2-38）
或
（2-39）
式中 τ0为屈服应力或极限动切应力，Pa；ηp为塑性粘度或结构粘度，Pa·s。
宾汉流体的视粘度为
（2-40）
上式表明，视粘度η是剪切速率的单调减函数，且随着剪切速率的增大逐渐趋于塑性粘度ηp。
2.幂律模型
n=1
n>1
n<1
剪切速率
图2-9 幂律流体的流变特性曲线
粘性切应力τ
符合幂律模型的流体简称幂律流体，其流变曲线如图2-9所示，幂律模型常被用来描述拟塑性流体和膨胀性流体的流变性，其本构方程为
（2-41）
式中 k为稠度系数，Pa·sn；n为流性指数，无量纲指数。对于拟塑性型（剪切稀化）流体来说，n＜1；对于膨胀型（剪切稠化）流体来说，n＞1；当n=1时则为牛顿流体；n与1的差值越大则流体的非牛顿特性越强。
幂律流体的视粘度为
（2-42）
由此可见，幂律流体的视粘度与稠度系数k成正比，当n＞1时为剪切速率的单增函数，当n＜1时为剪切速率的单减函数。
3.卡森模型
图2-9 幂律流体的流变曲线
石油工程中的另一个常用的带屈服值的流变模型为卡森模型，是卡森于1959年提出的，其本构方程为
（1-23）
式中 ηc为卡森粘度，Pa·s；τc为卡森屈服值，Pa。满足卡森模型的流体称为卡森流体，卡森流体的流变曲线如图2-10所示。
卡森流体的视粘度为
（2-43）
图2-9 Herschel—Buckley流体
的流变特线
粘性切应力τ
剪切速率
4.带屈服应力的幂律模型
上述三种模式均是采用双参数来描述非牛顿流体的流变特性的，另外还有一些三参数的流变模型。这里首先介绍图2-10所示的带屈服应力的幂律模型，又称为Herschel—Buckley模型，满足这一模型的流体称为Herschel—Buckley流体，其本构方程为
（2-44）
该模型相当于宾汉模型和幂律模型的综合，因此本构方程中各量的意义等同于宾汉模型和幂律模型中各参数的意义。
带屈服应力的幂律流休的视粘度为
（2-45）
5.罗伯逊一斯蒂夫模型
罗伯逊-斯蒂夫模型的形式为
（2-46）
式中 A为稠度系数，Pa·sB；B为流性指数，无因次；C为变形速度校正值，s-1。
罗伯逊-斯蒂夫模型也是一个三参数描述带屈服应力的非牛顿流体流变性的模型。
罗伯逊-斯蒂夫流体的视粘度为
（2-47）
五 初始条件和边界条件
柯西应力方程是一个与时间有关的非线性偏微分方程组，理论上只要方程组封闭便可得到方程组的通解，且在给定的初边值条件下还可以得到方程组的特解。工程实际中的很多的场都是有限的，初始条件和边界条件也都是给定的。因此，理论上可以得到任何流动的速度场、应变场和应力场。
解析方法是流体力学各种研究方法中最为准确的和最为理想的方法。解析方法首先要详细分析问题的物理学本质，通过适当的简化建立物理模型，之后运用物理定律建立数学模型，通常是建立起微分方程或微分方程组，确定流动方程边界条件和初始条件，最后运用数学方法求解出流动方程的解析解并与实验方法所得的结果进行比较，以检验物理模型和数学模型的合理性。
前面已经介绍了流动方程，后两种方法中，都离不开流场的初始条件和边界条件，下面就来简单介绍一下非牛顿流体力学中所涉及的初始条件和边界条件。
（一）初始条件
初始条件是指流动在t=t0的初始时刻，流体运动应满足的初始状态，或是给出流场中各物理量及其对时间的导数值。如速度、压强、密度和温度的初始值为
（2-48）
其中u0(x, y, z)、p0(x, y, z)和ρ0(x, y, z)为已知函数。
（二）边界条件
边界条件是指流体运动的边界上方程组的解应满足的条件。边界条件有很多种形式，如固壁边界条件、无穷远条件及对称边界条件等。
1. 固壁边界
固壁条件是最常用的边界条件，存在于固体和流体之间。考虑到流体不能流入固体也不能沿固体壁面滑移，因此固壁条件可表述为
（2-49）
式中uf和us分别表示流体和固体在边界处的速度。这一条件也称为粘附条件或无滑移条件。固壁静止时可表示为
（2-50）
2. 无穷远条件
无穷远条件指在无限大的流场中（如飞机所处的流场），无穷远处的边界条件为无穷远处的相应的流动参数值。即当r→∞时，
（2-51）
其中u∞(x, y, z, t)、p∞(x, y, z, t)和ρ∞(x, y, z, t)为已知函数。
3. 对称边界条件
在求解流体力学问题时经常会遇到流场对称的情况，为了减小工作量常常仅仅求解其中的半个流场，而没有必要对整个流场进行求解。在运用数值方法求解时这种简化会带来更高的计算效率。假设研究对象是平面上关于y轴对称的流动，这样我们可以仅仅研究左半个或者右半个流场，得到其中任何半个流场则另一半的流动情况也就不用在进行研究了。

习 题
1.试由纯粘流体的本构方程和柯西方程推导纳维尔-斯托克斯方程（N-S方程）。