黏性是流体的重要属性之一,自然界中存在的流体都具有黏性。在管内流动中,如果黏性影响充满整个管道,就必须考虑整个管道的黏性流动。对于黏性影响比较小的情况,有时可以将真实流体近似地按理想流体来处理。而在某些情况下,例如求表面摩擦阻力的问题,流体的黏性起着重要的作用而不能忽略。理论和实验表明,对于气体绕物体的流动,黏性影响主要在靠近物体表面的薄层内(称为附面层)。这样,求解黏性流动的问题,既可以建立并求解黏性流动的基本万程,也可以求解附面层内的流动。因此研究附面层的目的,一方面是为了计算气流绕物体的摩擦阻力,而另一方面则是估算物体上各点的热流量,从而寻求减小摩擦阻力,减轻气动热,采取必要的设计措施。

本章首先讨论黏性流动的基本方程,由于连续方程并不涉及黏性问题,因此本章主要讨论动量方程和能量方程,然后导出湍流流动的雷诺方程,最后讨论附面层基本知识。本章内容构成了黏性流体流动的基本知识。

10.1   微分形式的动量方程 ——N一S方程

针对微元控制体(见图10.1),根据第3.5节,可以列出动量方程(式(3.45)),即

\(\sum F\text{=}\frac{\partial }{\partial t}\left( \int_{cv}{\rho Vdv} \right)+\sum {{\left( {{q}_{mi}}{{V}_{i}} \right)}_{out}}-\sum {{\left( {{q}_{mi}}{{V}_{i}} \right)}_{in}}\)                  (10.1)

同样,由于控制体为微元体,所以式(10.1)中的积分可以表示为

\(\frac{\partial }{\partial t}\left( \int_{cv}{\rho Vdv} \right)\approx \frac{\partial \left( \rho V \right)}{\partial t}dxdydz\)             (10.2)

动量通量发生在6个面上,从3个面流入,3个面流出,表10.1中给出了各个控制面上流入或流出控制体的动量。

将表中结果和式(10.2)代入式(10.1),得

\(\sum F\text{=}\left[ \frac{\partial \left( \rho V \right)}{\partial t}+\frac{\partial \left( \rho {{V}_{x}}V \right)}{\partial x}+\frac{\partial \left( \rho {{V}_{y}}V \right)}{\partial y}+\frac{\partial \left( \rho {{V}_{z}}V \right)}{\partial z} \right]dxdydz\)             (10.3)

式(10.3)为矢量方程,等号右边中括号内可以改写成

\(\frac{\partial \left( \rho V \right)}{\partial t}+\frac{\partial \left( \rho {{V}_{x}}V \right)}{\partial x}+\frac{\partial \left( \rho {{V}_{y}}V \right)}{\partial y}+\frac{\partial \left( \rho {{V}_{z}}V \right)}{\partial z}=\)V[∂ρ/∂t+(ρV)]+\(\rho \left( \frac{\partial V}{\partial t}+{{V}_{x}}\frac{\partial V}{\partial x}+{{V}_{y}}\frac{\partial V}{\partial y}+{{V}_{z}}\frac{\partial V}{\partial z} \right)\)                       (10.4)

根据连续方程,式(10.4)等号右边中括号内为零,第二大项括号内为加速度,因此方程式(10.3)可以写为

\(\sum F=\rho \frac{dV}{dt}dxdydz\)            (10.5)

式(10.5)说明,微元控制体内流体的加速度乘以控制体内流体的质量,等于控制体所受的合外力。控制体所受的外力有两大类:质量力和表面力。质量力是在某种外部场的作用下使得所有流体质量受到的力,如重力、离心力、电磁力等等。表面力是由于控制面上应力的作用而产生的力,这些应力包括压强p和流体运动而产生的黏性应力τij

\({{\sigma }_{ij}}=\left[ \begin{align}& -p+{{\tau }_{xx}},{{\tau }_{yx}},{{\tau }_{zx}} \\& {{\tau }_{xy}},-p+{{\tau }_{yy}},{{\tau }_{zy}} \\& {{\tau }_{xz}},{{\tau }_{yz}},-p+{{\tau }_{zz}} \\\end{align} \right]\)               (10.6)

σij表示在与i(x,y,z)方向垂直的面上j(x,y,z)方向的应力。

下面来分析控制体所受表面力的合力。为了简单起见,以x方向为例。图10.2给出了6个面上x方向应力作用的表面力。

将这些力进行矢量和可得出微元控制体所受表面力在x方向的分量为

\(d{{F}_{x,surf}}=\left[ \frac{\partial {{\sigma }_{xx}}}{\partial x}+\frac{\partial {{\sigma }_{yx}}}{\partial y}+\frac{\partial {{\sigma }_{zx}}}{\partial z} \right]dxdydz\)               (10.7)

将式(10.6)的第一行代入,两边同除以dv=dxdydz,得

\(\frac{d{{F}_{x}}}{dv}=-\frac{\partial p}{\partial x}+\frac{\partial {{\tau }_{xx}}}{\partial x}+\frac{\partial {{\tau }_{yx}}}{\partial y}+\frac{\partial {{\tau }_{zx}}}{\partial z}\)               (10.8a)

同理,可以得出y,z方向的合力分别为

\(\frac{d{{F}_{y}}}{dv}=-\frac{\partial p}{\partial y}+\frac{\partial {{\tau }_{xy}}}{\partial x}+\frac{\partial {{\tau }_{yy}}}{\partial y}+\frac{\partial {{\tau }_{zy}}}{\partial z}\)               (10.8b)

\(\frac{d{{F}_{z}}}{dv}=-\frac{\partial p}{\partial z}+\frac{\partial {{\tau }_{xz}}}{\partial x}+\frac{\partial {{\tau }_{yz}}}{\partial y}+\frac{\partial {{\tau }_{zz}}}{\partial z}\)               (10.8c)

将式(10.8)写成矢量形式为

\({{\left( \frac{dF}{dv} \right)}_{surf}}=-\nabla p+{{\left( \frac{dF}{dv} \right)}_{viscous}}\)               (10.9)

式(10.9)等号右边第二项为黏性力项,由9个分量组成,即

\(\begin{align}& {{\left( \frac{dF}{dv} \right)}_{viscous}}=i\left( \frac{\partial {{\tau }_{xx}}}{\partial x}+\frac{\partial {{\tau }_{yx}}}{\partial y}+\frac{\partial {{\tau }_{zx}}}{\partial z} \right)+ \\& j\left( \frac{\partial {{\tau }_{xy}}}{\partial x}+\frac{\partial {{\tau }_{yy}}}{\partial y}+\frac{\partial {{\tau }_{zy}}}{\partial z} \right)+k\left( \frac{\partial {{\tau }_{xz}}}{\partial x}+\frac{\partial {{\tau }_{yz}}}{\partial y}+\frac{\partial {{\tau }_{zz}}}{\partial z} \right) \\\end{align}\)               (10.10)

式(10.10)还可以简写成如下的散度形式:

\({{\left( \frac{dF}{dv} \right)}_{viscous}}=\)τij               (10.11)

式中

\({{\tau }_{ij}}=\left[ \begin{align}& {{\tau }_{xx}},{{\tau }_{yx}},{{\tau }_{zx}} \\& {{\tau }_{xy}},{{\tau }_{yy}},{{\tau }_{zy}} \\& {{\tau }_{xz}},{{\tau }_{yz}},{{\tau }_{zz}} \\\end{align} \right]\)               (10.12)

τij称为黏性应力张量,且为对称张量,即τij=τji,当i≠j时,因此该张量有6个独立分量。表面力的合力包含压强梯度和黏性应力散度两部分。将式(1011)、式(10.9)代入式(10.5)最后得出对于无限小微元体的微分形式动量方程为

    ρRp+τij=ρdV/dt             (10.13)

式中,ρR为单位体积所受的质量力。

用文字表示该方程的物理意义为

单位体积所受的质量力+单位体积所受的压力+

单位体积所受的黏性力=密度×加速度                      (10.14)

将方程式(10.13)写成分量形式为

\(\left. \begin{align}& \rho X-\frac{\partial p}{\partial x}+\frac{\partial {{\tau }_{xx}}}{\partial x}+\frac{\partial {{\tau }_{yx}}}{\partial y}+\frac{\partial {{\tau }_{zx}}}{\partial z}=\rho \left( \frac{\partial {{V}_{x}}}{\partial t}+{{V}_{x}}\frac{\partial {{V}_{x}}}{\partial x}+{{V}_{y}}\frac{\partial {{V}_{x}}}{\partial y}+{{V}_{z}}\frac{\partial {{V}_{x}}}{\partial z} \right) \\& \rho Y-\frac{\partial p}{\partial y}+\frac{\partial {{\tau }_{xy}}}{\partial x}+\frac{\partial {{\tau }_{yy}}}{\partial y}+\frac{\partial {{\tau }_{zy}}}{\partial z}=\rho \left( \frac{\partial {{V}_{y}}}{\partial t}+{{V}_{x}}\frac{\partial {{V}_{y}}}{\partial x}+{{V}_{y}}\frac{\partial {{V}_{y}}}{\partial y}+{{V}_{z}}\frac{\partial {{V}_{y}}}{\partial z} \right) \\& \rho Z-\frac{\partial p}{\partial z}+\frac{\partial {{\tau }_{xz}}}{\partial x}+\frac{\partial {{\tau }_{yz}}}{\partial y}+\frac{\partial {{\tau }_{zz}}}{\partial z}=\rho \left( \frac{\partial {{V}_{z}}}{\partial t}+{{V}_{x}}\frac{\partial {{V}_{z}}}{\partial x}+{{V}_{y}}\frac{\partial {{V}_{z}}}{\partial y}+{{V}_{z}}\frac{\partial {{V}_{z}}}{\partial z} \right) \\\end{align} \right\}\)               (10.15)

式(10.15)为以应力形式表示的黏性流体运动微分方程。

对于无黏性流动τij=0,因此方程式(10.13)变成

      ρRp=ρdV/dt

该式即为描述理想流动的欧拉方程(Euler ‘ s equation)。

对于牛顿流体,黏性应力与流体的变形以及黏度成正比。具体关系式为

\(\left. \begin{align}& {{\tau }_{xx}}=2\mu \frac{\partial {{V}_{x}}}{\partial x}-\frac{2}{3}\mu \left( \nabla \cdot V \right) \\& {{\tau }_{yy}}=2\mu \frac{\partial {{V}_{y}}}{\partial y}-\frac{2}{3}\mu \left( \nabla \cdot V \right) \\& {{\tau }_{zz}}=2\mu \frac{\partial {{V}_{z}}}{\partial z}-\frac{2}{3}\mu \left( \nabla \cdot V \right) \\& {{\tau }_{xy}}=\mu \left( \frac{\partial {{V}_{y}}}{\partial x}+\frac{\partial {{V}_{x}}}{\partial y} \right) \\& {{\tau }_{yz}}=\mu \left( \frac{\partial {{V}_{z}}}{\partial y}+\frac{\partial {{V}_{y}}}{\partial z} \right) \\& {{\tau }_{zx}}=\mu \left( \frac{\partial {{V}_{x}}}{\partial z}+\frac{\partial {{V}_{z}}}{\partial x} \right) \\\end{align} \right\}\)               (10.16)

式(10.16)称为广义牛顿内摩擦定律。将式(10.16)代入式(10.15),当μ=常数时,可得

\(\rho \frac{d{{V}_{x}}}{dt}=\rho X-\frac{\partial p}{\partial x}+\mu \left( \frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{z}^{2}}} \right)+\frac{\mu }{3}\frac{\partial \left( \nabla \cdot V \right)}{\partial x}\)               (10.17a)

\(\rho \frac{d{{V}_{y}}}{dt}=\rho Y-\frac{\partial p}{\partial y}+\mu \left( \frac{{{\partial }^{2}}{{V}_{y}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}{{V}_{y}}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}{{V}_{y}}}{\partial {{z}^{2}}} \right)+\frac{\mu }{3}\frac{\partial \left( \nabla \cdot V \right)}{\partial y}\)               (10.17b)

\(\rho \frac{d{{V}_{z}}}{dt}=\rho Z-\frac{\partial p}{\partial z}+\mu \left( \frac{{{\partial }^{2}}{{V}_{z}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}{{V}_{z}}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}{{V}_{z}}}{\partial {{z}^{2}}} \right)+\frac{\mu }{3}\frac{\partial \left( \nabla \cdot V \right)}{\partial z}\)               (10.17c)

式(10.17)即为描述牛顿黏性流体运动的微分方程式,又称为纳维尔一斯托克斯 (Navier-Stokes)方程,简称N一S方程。它是由C.L.M.H.Navier(1785一1836)和Sir George G Stokes(1819一1903)分别独立导出的,方程即以他们的名字联合命名。

该方程可以与成矢量形式,即

\(\rho \frac{dV}{dt}=\rho R-\nabla p+\mu {{\nabla }^{2}}V+\frac{\mu }{3}\nabla \left( \nabla \cdot V \right)\)               (10.18)

对于不可压流动,式(10.18)为

\(\frac{dV}{dt}=R-\frac{1}{\rho }\nabla p+\upsilon {{\nabla }^{2}}V\)               (10.19)

式中,v=μ/ρ称为运动黏度;

\({{\nabla }^{2}}=\frac{{{\partial }^{2}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}}{\partial {{z}^{2}}}\)

N一s方程为二阶非线性偏微分方程组。在一般情况下,从数学上精确求解此方程是不可能的。但是对于一些简单的流动,如平行平板的定常层流流动、圆管内的定常层流流动等是可以得到精确解的,而且这些精确解与实验结果符合一致。

不可压流动圆柱坐标下的N一s方程为

\(\left. \begin{align}& \rho \left( \frac{d{{V}_{r}}}{dt}-\frac{V_{\theta }^{2}}{r} \right)=\rho {{R}_{r}}-\frac{\partial p}{\partial r}+\mu \left( {{\nabla }^{2}}{{V}_{r}}-\frac{2}{{{r}^{2}}}\frac{\partial {{V}_{\theta }}}{\partial \theta }-\frac{{{V}_{r}}}{{{r}^{2}}} \right) \\& \rho \left( \frac{d{{V}_{\theta }}}{dt}-\frac{{{V}_{r}}{{V}_{\theta }}}{r} \right)=\rho {{R}_{\theta }}-\frac{\partial p}{r\partial \theta }+\mu \left( {{\nabla }^{2}}{{V}_{\theta }}-\frac{2}{{{r}^{2}}}\frac{\partial {{V}_{r}}}{\partial \theta }-\frac{{{V}_{\theta }}}{{{r}^{2}}} \right) \\& \rho \left( \frac{d{{V}_{z}}}{dt} \right)=\rho {{R}_{z}}-\frac{\partial p}{\partial z}+\mu \left( {{\nabla }^{2}}{{V}_{z}} \right) \\\end{align} \right\}\)               (10.20)

式中                                                       \({{\nabla }^{2}}=\frac{1}{r}\frac{\partial }{\partial r}\left( r\frac{\partial }{\partial r} \right)+\frac{1}{{{r}^{2}}}\frac{{{\partial }^{2}}}{\partial {{\theta }^{2}}}+\frac{{{\partial }^{2}}}{\partial {{z}^{2}}}\)

10.2  微分形式的能量方程

类似于第3.7节,由式( 3.71 )同样可以针对微元控制体列出能量方程

\(\overset{\cdot }{\mathop{Q}}\,-{{\overset{\cdot }{\mathop{W}}\,}_{s}}-\overset{\cdot }{\mathop{{{W}_{v}}}}\,=\frac{\partial }{\partial t}\left( \int_{cv}{e\rho dv} \right)+\int_{cs}{\left( e+\frac{p}{\rho } \right)\rho \left( V\cdot n \right)dA}\)               (10.21)

式中                                       e=u+V²/2+gz

因为在微元控制体中没有轴功,所以\({{\overset{\cdot }{\mathop{W}}\,}_{s}}=0\)。采用与导出式( 10.3 )完全相同的方法,可以得岀

\(\overset{\cdot }{\mathop{Q}}\,-{{\overset{\cdot }{\mathop{W}}\,}_{v}}=\left[ \frac{\partial \left( e\rho  \right)}{\partial t}+\frac{\partial \left( \rho {{V}_{x}}\zeta  \right)}{\partial x}+\frac{\partial \left( \rho {{V}_{y}}\zeta  \right)}{\partial y}+\frac{\partial \left( \rho {{V}_{z}}\zeta  \right)}{\partial z} \right]dxdydz\)               (10.22)

式中,ζ=e+p/ρ。类似于式(10.4),考虑到连续方程,式(10.22)成为

\(\overset{\cdot }{\mathop{Q}}\,-{{\overset{\cdot }{\mathop{W}}\,}_{v}}=\left( \rho \frac{de}{dt}+p\nabla \cdot V \right)dxdydz\)               (10.23)

传热量\(\overset{\cdot }{\mathop{Q}}\,\)可以分为两大类:一类是由于热传导对微元控制体的传热;另一类是辐射、化学反应等其他形式的热量传递。用q来表示热传导对控制体內流体在单位时间內单位面积上的传热量。下面推导由于热传导而产生的传热量。根据傅里叶热传导定律,有

q=-λ▽T

式中,λ为导热系数,与分析质量流率和动量流率相同,可以得岀6个面上由于热传导而产生的热流率,见表10.2。

将6个面上的热流量代数求和,得出

\({{\overset{\cdot }{\mathop{Q}}\,}_{k}}=-\left[ \frac{\partial {{q}_{x}}}{\partial x}+\frac{\partial {{q}_{y}}}{\partial y}+\frac{\partial {{q}_{z}}}{\partial z} \right]dxdydz=-\nabla \cdot qdxdydz\)               (10.24)

将傅里叶热传导定律代入式(10.24),得出

\({{\overset{\cdot }{\mathop{Q}}\,}_{k}}=\nabla \cdot \left( \lambda \nabla T \right)dxdydz\)               (10.25)

黏性应力所做功的功率(以下简称黏性应力的功率)等于黏性应力分量、相应的速度分量和相应的面积三项的乘积(见图10.3),与x方向垂直的左侧面上黏性应力的功率为

\({{\overset{\cdot }{\mathop{W}}\,}_{v,LF}}={{\omega }_{x}}dydz\)

式中

\({{\omega }_{x}}=-\left( {{V}_{x}}{{\tau }_{xx}}+{{V}_{y}}{{\tau }_{xy}}+{{V}_{z}}{{\tau }_{xz}} \right)\)               (10.26)

与上述分析质量流量、动量流量和热流量完全相同,可以得出在与x方向垂直的两个面上黏性应力的功率为

\[{{\overset{\cdot }{\mathop{W}}\,}_{vx}}=-\frac{\partial }{\partial x}\left( {{V}_{x}}{{\tau }_{xx}}+{{V}_{y}}{{\tau }_{xy}}+{{V}_{z}}{{\tau }_{xz}} \right)dxdydz\]

同理,可以得出另外两个方向上的功率,因此总的黏性应力的功率应为

\(\begin{align}& {{\overset{\cdot }{\mathop{W}}\,}_{v}}=-\left[ \frac{\partial }{\partial x}\left( {{V}_{x}}{{\tau }_{xx}}+{{V}_{y}}{{\tau }_{xy}}+{{V}_{z}}{{\tau }_{xz}} \right)+\frac{\partial }{\partial y}\left( {{V}_{x}}{{\tau }_{yx}}+{{V}_{y}}{{\tau }_{yy}}+{{V}_{z}}{{\tau }_{yz}} \right)+\frac{\partial }{\partial z}\left( {{V}_{x}}{{\tau }_{zx}}+{{V}_{y}}{{\tau }_{zy}}+{{V}_{z}}{{\tau }_{zz}} \right) \right]dxdydz= \\& -\nabla \cdot \left( V\cdot {{\tau }_{ij}} \right)dxdydz \\\end{align}\)               (10.27)

将式(10.27)、式(10.25)代入式(10.23),便得到微分形式的能量方程为

\(\rho \frac{de}{dt}+p\nabla \cdot V=\nabla \cdot \left( \lambda \nabla T \right)+\nabla \cdot \left( V\cdot {{\tau }_{ij}} \right)+\rho {{q}^{\prime }}\)               (10.28)

式中,q′表示由于热辐射或其他原因在单位时间内传入的热量。若忽略质量力,则

e=u+V²/2

式(10.28)中黏性力做功项还可以分解为

\(\nabla \cdot \left( V\cdot {{\tau }_{ij}} \right)=V\cdot \left( \nabla \cdot {{\tau }_{ij}} \right)+\Phi \)               (10.29)

式中,φ为黏性耗散函数。该耗散函数为

\(\begin{align}& \Phi =\mu \left[ 2{{\left( \frac{\partial {{V}_{x}}}{\partial x} \right)}^{2}}+2{{\left( \frac{\partial {{V}_{y}}}{\partial y} \right)}^{2}}+2{{\left( \frac{\partial {{V}_{z}}}{\partial z} \right)}^{2}}+{{\left( \frac{\partial {{V}_{y}}}{\partial x}+\frac{\partial {{V}_{x}}}{\partial y} \right)}^{2}}+{{\left( \frac{\partial {{V}_{z}}}{\partial y}+\frac{\partial {{V}_{y}}}{\partial z} \right)}^{2}}+{{\left( \frac{\partial {{V}_{x}}}{\partial z}+\frac{\partial {{V}_{z}}}{\partial x} \right)}^{2}} \right]- \\& \frac{2}{3}\mu {{\left( \frac{\partial {{V}_{x}}}{\partial x}+\frac{\partial {{V}_{y}}}{\partial y}+\frac{\partial {{V}_{z}}}{\partial z} \right)}^{2}} \\\end{align}\)               (10.30)

通过式(10.30)可以看出φ≥ 0,也就是说耗散项永远是正的,即黏性应力所做的功总是消耗机械能,使流体的内能增加。

将式(10.29)代入到式(10.28)中,并采用式(10.13)消去▽•τij,得到内能形式的能量方程为

\(\rho \frac{du}{dt}+p\left( \nabla \cdot V \right)=\nabla \cdot \left( \lambda \nabla T \right)+\Phi +\rho {{q}^{\prime }}\)               (10.31)

根据连续方程,有

\(p\left( \nabla \cdot V \right)=-\frac{p}{\rho }\frac{d\rho }{dt}=p\rho \frac{d}{dt}\left( \frac{1}{\rho } \right)\)               (10.32)

它表示单位时间内单位体积流体在压强p的作用下所做的膨胀(或压缩)功。

对于完全气体,由热力学公式

\(T\frac{ds}{dt}=\frac{du}{dt}+p\frac{d}{dt}\left( \frac{1}{\rho } \right)\)               (10.33)

\(\frac{dh}{dt}=\frac{du}{dt}+p\frac{d}{dt}\left( \frac{1}{\rho } \right)+\frac{1}{\rho }\frac{dp}{dt}\)               (10.34)

因此,可以将式(10.31)写成熵或焓的形式,即

\(\rho T\frac{ds}{dt}=\nabla \cdot \left( \lambda \nabla T \right)+\Phi +\rho {{q}^{\prime }}\)               (10.35)

\(\rho \frac{dh}{dt}=\frac{dp}{dt}+\nabla \cdot \left( \lambda \nabla T \right)+\Phi +\rho {{q}^{\prime }}\)               (10.36)

注意到

\(du={{c}_{v}}dT,dh={{c}_{p}}dT\)               (10.37)

设cv,cp,λ不变,则式(10.31)和式(10.36)又可以与成用温度表示的能量方程,即

\(\rho {{c}_{v}}\frac{dT}{dt}=-p\left( \nabla \cdot V \right)+\lambda {{\nabla }^{2}}T+\Phi +\rho {{q}^{\prime }}\)               (10.38)

\(\rho {{c}_{p}}\frac{dT}{dt}=\frac{dp}{dt}+\lambda {{\nabla }^{2}}T+\Phi +\rho {{q}^{\prime }}\)               (10.39)

10.3   初始条件和边界条件

通过上述的推导,得出了描述牛顿流体运动的微分方程组,共5个方程,包括连续方程(1个),动量方程(3个),能量方程(1个),而未知量有6个,即ρ,Vx,Vy,Vz,p,u(T)(以直角坐标为例,柱坐标结果一样)。因此,方程并不封闭,还需要补充一个热力学的关系式,即完全气体状态方程

p=ρRT               (10.40)

这样包括状态方程在内,基本方程组共有6个方程,构成封闭的方程组。但是要得到具体的解还要给定相应的初始和边界条件,这些条件统称为定解条件。

1.初始条件

在初始时刻,方程组的解应该等于该时刻给定的数值。在数学上可以表示为

当t=t0时,

\(\left. \begin{align}& V\left( x,y,z,{{t}_{0}} \right)={{V}_{0}}\left( x,y,z \right) \\& p\left( x,y,z,{{t}_{0}} \right)={{p}_{0}}\left( x,y,z \right) \\& \rho \left( x,y,z,{{t}_{0}} \right)={{\rho }_{0}}\left( x,y,z \right) \\& T\left( x,y,z,{{t}_{0}} \right)={{T}_{0}}\left( x,y,z \right) \\\end{align} \right\}\)               (10.41)

式中,V0(x,y,z),p0(x,y,z),ρ0(x,y,z),T0(x,y,z)均为t0时刻的己知函数。

2.边界条件

在运动流体的边界上,方程组的解所应满足的条件称为边界条件。边界条件随具体问题而定,一般来讲可能有以下几种情况:固体壁面(包括可渗透壁面)上的边界条件;不同流体分界面(包括自由液面、气液界面、液液界面)上的边界条件;无限远或管道进、出口处的边界条件等。

对于不可渗漏的固体边界速度为无滑移条件、温度为无突跃条件,即

\({{V}_{fluid}}={{V}_{wall}}{{T}_{fluid}}={{T}_{wall}}\)               (10.42)

如果固体边界为可渗漏,则边界条件要根据具体情况来确定。

对于所有的流动进、出囗截面,应给出每时刻截面上速度、压力和温度的分布。对于流体绕流物体的问题,进出囗边界变成了无穷远边界,应给出无穷远边界条件。

例10.1  用N一S方程导出圆管内的层流不可压缩流动的速度分布。

解:在第四章,我们曾经就圆管内不可压黏性层流流动,通过建立动力学基本方程,得到了该问题的解析解。这一问题还可以通过N一S方程来求解。

考虑充分发展的圆管内的不可压层流流动,圆管半径r0,直径为d,采用柱坐标。设管轴为z轴,r方向与z轴垂直,则仅有方向的分速度Vz=Vz(r)。取微元控制体,则由式(10.20),N一S方程可以写为

\(-\frac{dp}{dz}+\mu {{\nabla }^{2}}{{V}_{z}}=-\frac{dp}{dz}+\frac{\mu }{r}\frac{d}{dr}\left( r\frac{d{{V}_{z}}}{dr} \right)=0\)

对上式两边积分,考虑到ρ与r无关,则积分两次后,得

\({{V}_{z}}=\frac{1}{\mu }\frac{dp}{dz}\frac{{{r}^{2}}}{4}+{{C}_{1}}\ln r+{{C}_{2}}\)

式中的积分常数可以通过边界条件确定。当r=r0时,Vz(r0)=0;当r=0时,Vz(0)=Vmax,代入上式得

\(\begin{align}& {{C}_{1}}=0 \\& {{C}_{2}}=-\frac{dp}{dz}\frac{r_{0}^{2}}{4\mu } \\\end{align}\)

所以圆管内的速度分布为

\({{V}_{z}}=\frac{1}{4\mu }\frac{dp}{dz}\left( {{r}^{2}}-r_{0}^{2} \right)\)

例10.2   两个相距h的无限大平板,如图10.4所示,不可压缩黏性流体在其中作一维不可压定常层流流动。用N—S方程求解其速度分布规律。

解:设x方向沿流动方向,y方向垂直于平板。由于平板无限大,所以速度与坐标x,z无关,即V=Vx(y),而压强与y,z坐标无关。对该流动写出N—S方程,即由式(10.17a)得

\(-\frac{1}{\rho }\frac{dp}{dx}+\upsilon \frac{{{d}^{2}}V\left( y \right)}{d{{y}^{2}}}=0\)

由于dp/dx仅是x的函数,与y无关,因此对上式积分两次,得

\(V\left( y \right)=\frac{1}{\mu }\frac{dp}{dx}\left( \frac{{{y}^{2}}}{2}+{{C}_{1}}y+{{C}_{2}} \right)\)

当y=0或y=h时,V(y)=0,代入上式得C1=-h/2,C2=0,所以库埃特黏性流动的速度分布为

\(V\left( y \right)=-\frac{1}{2\mu }\frac{dp}{dx}y\left( h-y \right)\)

如果将上式写成无量纲的形式,并引进无量纲的压力参数P,则上式写为

\(\frac{V}{{{V}_{\infty }}}=P\frac{y}{h}\left( 1-\frac{y}{h} \right)\)

式中,\(P=-\frac{{{h}^{2}}}{2\mu {{V}_{\infty }}}\frac{dp}{dx}\),V∞为特征速度。从上式可以看出,库埃特黏性流动的速度分布为抛物线规律。由于平板不动,因此在两板壁面处,速度为零;在两板中间(y/h=1/2)处,速度达到最大值。

对于可压黏性流动,由于速度、压强、密度和温度都在变化,因此求解的方程除了连续方程和动量方程外,还必须考虑能量方程和状态方程。求出速度、压强和温度分布后,可以进一步求出平板的阻力和热交换。关于可压的库埃特黏性流动,由于求解复杂,此处不再赘述。

10.4   雷诺方程和雷诺应力

在第四章中曾经对湍流流动的速度分布、流动特点和流动损失等作了简单的讨论。如果要知道湍流流场中的流动细节,即计算流场中各点的流动参数,就需要建立适合于湍流流动的基本方程。本节就是要导出湍流流动的雷诺方程。

从对湍流的研宄可知,湍流运动中任何物理量都随时间和空间不断地变化,所以要想用 N—S方程求解这种运动的瞬时速度是非常困难的。研宄表明,虽然湍流运动十分复杂,但是它仍然遵循连续介质运动的特征和一般力学规律,因此,雷诺提出用时均值概念来研究湍流运动的方法,导出了以时间平均速度场为基础的雷诺时均N—S方程。

雷诺从不可压缩流体的N—S方程导出湍流平均运动方程(后人称此为雷诺方程),并引出雷诺应力的概念。之后,人们引用时均值概念导出湍流基本方程,使湍流运动的理论分析得到了很大的发展。

10.4.1  常用的时均运算关系式

设A,B,C为湍流中物理量的瞬时值,\(\overline{A},\overline{B},\overline{C}\)为物理量的时均值,A′,B′,C′为物理量的脉动值,则具有下述时均运算规律。

(1)时均量的时均值等于原来的时均值,即

\(\overline{\overline{A}}=\overline{A}\)               (10.43)

因为在时间平均周期T内\(\overline{A}\)是个定值,所以其时均值仍为原来的值。

(2)脉动量的时均值等于零,即

\(\overline{{{A}^{\prime }}}=0\)               (10.44)

\(\overline{{{A}^{\prime }}}=\frac{1}{T}\int_{0}^{T}{{{A}^{\prime }}dt=}\frac{1}{T}\int_{0}^{T}{\left( A-\overline{A} \right)dt=}\overline{A}-\overline{A}=0\)

(3)瞬时物理量之和的时均值,等于各个物理量时均值之和,即

\(\overline{A+B}=\overline{A}+\overline{B}\)               (10.45)

\(\overline{A+B}=\frac{1}{T}\int_{0}^{T}{\left( A+B \right)dt=}\frac{1}{T}\int_{0}^{T}{Adt}+\frac{1}{T}\int_{0}^{T}{Bdt}=\overline{A}+\overline{B}\)

(4)时均物理量与脉动物理量之积的时均值等于零,即

\(\overline{\overline{A}{{B}^{\prime }}}=0\)               (10.46)

因为在平均周期内\(\overline{A}\)是个定值,所以有

\(\overline{\overline{A}{{B}^{\prime }}}=\frac{1}{T}\int_{0}^{T}{\overline{A}{{B}^{\prime }}dt=}\overline{A}\frac{1}{T}\int_{0}^{T}{{{B}^{\prime }}dt=}\overline{A}\overline{{{B}^{\prime }}}=0\)

(5)时均物理量与瞬时物理量之积的时均值等于两个时均物理量之积,即

\(\overline{\overline{A}B}=\overline{A}\overline{B}\)               (10.47)

同样,在平均周期内\(\overline{A}\)是个定值,所以

\(\overline{\overline{A}B}=\frac{1}{T}\int_{0}^{T}{\overline{A}Bdt=}\overline{A}\frac{1}{T}\int_{0}^{T}{Bdt=}\overline{A}\overline{B}\)

(6)两个瞬时物理量之积的时均值,等于两个时均物理量之积与两个脉动量之积的时均值之和,即

\(\overline{AB}=\overline{A}\overline{B}+\overline{{{A}^{\prime }}{{B}^{\prime }}}\)               (10.48)

\(\begin{align}& \overline{AB}=\frac{1}{T}\int_{0}^{T}{ABdt=}\frac{1}{T}\int_{0}^{T}{\left( \overline{A}+{{A}^{\prime }} \right)\left( \overline{B}+{{B}^{\prime }} \right)dt=}\frac{1}{T}\int_{0}^{T}{\left( \overline{A}\overline{B}+{{A}^{\prime }}\overline{B}+{{B}^{\prime }}\overline{A}+{{A}^{\prime }}{{B}^{\prime }} \right)dt=} \\& \overline{A}\overline{B}+\overline{B}\frac{1}{T}\int_{0}^{T}{{{A}^{\prime }}dt+}\overline{A}\frac{1}{T}\int_{0}^{T}{{{B}^{\prime }}dt+}\frac{1}{T}\int_{0}^{T}{{{A}^{\prime }}{{B}^{\prime }}dt=}\overline{A}\overline{B}+\overline{{{A}^{\prime }}{{B}^{\prime }}} \\\end{align}\)

推论

\(\overline{ABC}=\overline{A}\overline{B}\overline{C}+\overline{A}\overline{{{B}^{\prime }}{{C}^{\prime }}}+\overline{B}\overline{{{A}^{\prime }}{{C}^{\prime }}}+\overline{C}\overline{{{A}^{\prime }}{{B}^{\prime }}}+\overline{{{A}^{\prime }}{{B}^{\prime }}{{C}^{\prime }}}\)               (10.49)

(7)瞬时物理量对空间坐标各阶导数的时均值,等于时均物理量对同一坐标的各阶导数,即

\(\left. \begin{align}& \frac{\overline{{{\partial }^{\prime \prime }}A}}{\partial {{s}^{\prime \prime }}}=\frac{{{\partial }^{\prime \prime }}\overline{A}}{\partial {{s}^{\prime \prime }}} \\& \frac{\overline{{{\partial }^{\prime \prime }}A}}{\partial {{s}^{\prime \prime }}}=\frac{1}{T}\int_{0}^{T}{\frac{{{\partial }^{\prime \prime }}A}{\partial {{s}^{\prime \prime }}}dt=}\frac{{{\partial }^{\prime \prime }}}{\partial {{s}^{\prime \prime }}}\left( \frac{1}{T}\int_{0}^{T}{Adt} \right)=\frac{{{\partial }^{\prime \prime }}\overline{A}}{\partial {{s}^{\prime \prime }}} \\\end{align} \right\}\)               (10.50)

式中,s代表任意坐标方向,如x,y,z。

推论                脉动量对空间坐标各阶导数的时均值等于零,即

\(\frac{\overline{{{\partial }^{\prime \prime }}{{A}^{\prime }}}}{\partial {{s}^{\prime \prime }}}=0\)               (10.51)

(8)瞬时物理量对于时间导数的时均值,等于时均物理量对时间的导数,即

\(\frac{\overline{\partial A}}{\partial t}\frac{\partial \overline{A}}{\partial t}\)               (10.52a)

在准定常的条件下,

\(\frac{\partial \overline{A}}{\partial t}=0\)               (10.52b)

10.4.2  湍流运动的连续方程

由于湍流流动中各物理量都具有某种统计特征的规律,所以基本方程中任一瞬间物理量都可用平均物理量和脉动物理量之和来代替,并且可以对整个方程进行时间平均的运算。

在湍流运动中,瞬时运动的速度应满足流体运动的基本方程。其连续方程为

\(\frac{\partial \rho }{\partial t}=\frac{\partial \left( \rho {{V}_{x}} \right)}{\partial x}+\frac{\partial \left( \rho {{V}_{y}} \right)}{\partial y}+\frac{\partial \left( \rho {{V}_{z}} \right)}{\partial z}=0\)

对其进行时均运算

\(\begin{align}& \overline{\frac{\partial \rho }{\partial t}+\frac{\partial \left( \rho {{V}_{x}} \right)}{\partial x}+\frac{\partial \left( \rho {{V}_{y}} \right)}{\partial y}+\frac{\partial \left( \rho {{V}_{z}} \right)}{\partial z}}=\overline{\frac{\partial \rho }{\partial t}}+\overline{\frac{\partial \left( \rho {{V}_{x}} \right)}{\partial x}}+\overline{\frac{\partial \left( \rho {{V}_{y}} \right)}{\partial y}}+\overline{\frac{\partial \left( \rho {{V}_{z}} \right)}{\partial z}}= \\& \frac{\partial \overline{\rho }}{\partial t}+\frac{\partial \overline{\left( \rho {{V}_{x}} \right)}}{\partial x}+\frac{\partial \overline{\left( \rho {{V}_{y}} \right)}}{\partial y}+\frac{\partial \overline{\left( \rho {{V}_{z}} \right)}}{\partial z}= \\& \frac{\partial \overline{\rho }}{\partial t}+\frac{\partial \left( \overline{\rho }{{\overline{V}}_{x}}+\overline{{{\rho }^{\prime }}V_{x}^{\prime }} \right)}{\partial x}+\frac{\partial \left( \overline{\rho }\overline{{{V}_{y}}}+\overline{{{\rho }^{\prime }}V_{y}^{\prime }} \right)}{\partial y}+\frac{\partial \left( \overline{\rho }\overline{{{V}_{z}}}+\overline{{{\rho }^{\prime }}V_{z}^{\prime }} \right)}{\partial z}= \\& \frac{\partial \overline{\rho }}{\partial t}+\frac{\partial \left( \overline{\rho }\overline{{{V}_{x}}} \right)}{\partial x}+\frac{\partial \left( \overline{{{\rho }^{\prime }}V_{x}^{\prime }} \right)}{\partial x}+\frac{\partial \left( \overline{\rho }\overline{{{V}_{y}}} \right)}{\partial y}+\frac{\partial \left( \overline{{{\rho }^{\prime }}V_{y}^{\prime }} \right)}{\partial y}+\frac{\partial \left( \overline{\rho }\overline{{{V}_{z}}} \right)}{\partial z}+\frac{\partial \left( \overline{{{\rho }^{\prime }}V_{z}^{\prime }} \right)}{\partial z} \\\end{align}\)

所以可压湍流运动的连续方程为

\(\frac{\partial \overline{\rho }}{\partial t}+\frac{\partial \left( \overline{\rho }\overline{{{V}_{x}}} \right)}{\partial x}+\frac{\partial \left( \overline{{{\rho }^{\prime }}V_{x}^{\prime }} \right)}{\partial x}+\frac{\partial \left( \overline{\rho }\overline{{{V}_{y}}} \right)}{\partial y}+\frac{\partial \left( \overline{{{\rho }^{\prime }}V_{y}^{\prime }} \right)}{\partial y}+\frac{\partial \left( \overline{\rho }\overline{{{V}_{z}}} \right)}{\partial z}+\frac{\partial \left( \overline{{{\rho }^{\prime }}V_{z}^{\prime }} \right)}{\partial z}=0\)

与瞬时值的连续方程相比,多出了三项脉动量乘积的时均值的导数。

对于不可压湍流运动,ρ=C,∂ρ/∂t=0,则连续方程可化为

\(\frac{\partial \overline{{{V}_{x}}}}{\partial x}+\frac{\partial \overline{{{V}_{y}}}}{\partial y}+\frac{\partial \overline{{{V}_{z}}}}{\partial z}=0\)               (10.53a)

并可得到

\(\frac{\partial V_{x}^{\prime }}{\partial x}+\frac{\partial V_{y}^{\prime }}{\partial y}+\frac{\partial V_{z}^{\prime }}{\partial z}=0\)               (10.53b)

可见,对不可压湍流运动,时均运动和脉动运动的连续方程和瞬时运动的连续方程具有相同的形式。

10.4.3  雷诺方程

对于不可压黏性流动,在不考虑质量力的情况下,湍流运动的瞬时速度应满足不可压缩黏性流体的N一S方程,即

\(\left. \begin{align}& \frac{\partial {{V}_{x}}}{\partial t}+{{V}_{x}}\frac{\partial {{V}_{x}}}{\partial x}+{{V}_{y}}\frac{\partial {{V}_{x}}}{\partial y}+{{V}_{z}}\frac{\partial {{V}_{x}}}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\upsilon \left( \frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{z}^{2}}} \right) \\& \frac{\partial {{V}_{y}}}{\partial t}+{{V}_{x}}\frac{\partial {{V}_{y}}}{\partial x}+{{V}_{y}}\frac{\partial {{V}_{y}}}{\partial y}+{{V}_{z}}\frac{\partial {{V}_{y}}}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial y}+\upsilon \left( \frac{{{\partial }^{2}}{{V}_{y}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}{{V}_{y}}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}{{V}_{y}}}{\partial {{z}^{2}}} \right) \\& \frac{\partial {{V}_{z}}}{\partial t}+{{V}_{x}}\frac{\partial {{V}_{z}}}{\partial x}+{{V}_{y}}\frac{\partial {{V}_{z}}}{\partial y}+{{V}_{z}}\frac{\partial {{V}_{z}}}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial z}+\upsilon \left( \frac{{{\partial }^{2}}{{V}_{z}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}{{V}_{z}}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}{{V}_{z}}}{\partial {{z}^{2}}} \right) \\\end{align} \right\}\)               (10.54a)

利用不可压流动瞬时运动的连续方程

\(\frac{\partial {{V}_{x}}}{\partial x}+\frac{\partial {{V}_{y}}}{\partial y}+\frac{\partial {{V}_{z}}}{\partial z}=0\)

可将式( 10.54a)改写成

\(\left. \begin{align}& \frac{\partial {{V}_{x}}}{\partial t}+\frac{\partial \left( {{V}_{x}}{{V}_{x}} \right)}{\partial x}+{{V}_{y}}\frac{\partial \left( {{V}_{x}}{{V}_{y}} \right)}{\partial y}+\frac{\partial \left( {{V}_{x}}{{V}_{z}} \right)}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\upsilon \left( \frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{z}^{2}}} \right) \\& \frac{\partial {{V}_{y}}}{\partial t}+\frac{\partial \left( {{V}_{x}}{{V}_{y}} \right)}{\partial x}+{{V}_{y}}\frac{\partial \left( {{V}_{y}}{{V}_{y}} \right)}{\partial y}+\frac{\partial \left( {{V}_{z}}{{V}_{y}} \right)}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial y}+\upsilon \left( \frac{{{\partial }^{2}}{{V}_{y}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}{{V}_{y}}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}{{V}_{y}}}{\partial {{z}^{2}}} \right) \\& \frac{\partial {{V}_{z}}}{\partial t}+\frac{\partial \left( {{V}_{x}}{{V}_{z}} \right)}{\partial x}+{{V}_{y}}\frac{\partial \left( {{V}_{y}}{{V}_{z}} \right)}{\partial y}+\frac{\partial \left( {{V}_{z}}{{V}_{z}} \right)}{\partial z}=-\frac{1}{\rho }\frac{\partial p}{\partial z}+\upsilon \left( \frac{{{\partial }^{2}}{{V}_{z}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}{{V}_{z}}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}{{V}_{z}}}{\partial {{z}^{2}}} \right) \\\end{align} \right\}\)               (10.54b)

然后对式(10.54b)中的第一式进行时间平均运算,则有

\(\frac{\partial \overline{{{V}_{x}}}}{\partial t}+\frac{\partial \left( \overline{{{V}_{x}}{{V}_{x}}} \right)}{\partial x}+{{V}_{y}}\frac{\partial \left( \overline{{{V}_{x}}{{V}_{y}}} \right)}{\partial y}+\frac{\partial \left( \overline{{{V}_{x}}{{V}_{z}}} \right)}{\partial z}=-\frac{1}{\rho }\frac{\partial \overline{p}}{\partial x}+\upsilon \left( \frac{{{\partial }^{2}}\overline{{{V}_{x}}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}\overline{{{V}_{x}}}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}\overline{{{V}_{x}}}}{\partial {{z}^{2}}} \right)\)               (10.55)

由于\({{V}_{x}}=\overline{{{V}_{x}}}+V_{x}^{\prime }\),应用时均物理量与脉动物理量之积的时均值等于零的运算规则,即\(\overline{\overline{A}{{B}^{\prime }}}=0,\overline{\overline{B}{{A}^{\prime }}}=0\),可得

\(\begin{align}& \overline{{{V}_{x}}{{V}_{x}}}=\overline{{{V}_{x}}}\overline{{{V}_{x}}}+\overline{V_{x}^{\prime }V_{x}^{\prime }} \\& \overline{{{V}_{x}}{{V}_{y}}}=\overline{{{V}_{x}}}\overline{{{V}_{y}}}+\overline{V_{x}^{\prime }V_{y}^{\prime }} \\& \overline{{{V}_{x}}{{V}_{z}}}=\overline{{{V}_{x}}}\overline{{{V}_{z}}}+\overline{V_{x}^{\prime }V_{z}^{\prime }} \\\end{align}\)

这样,式(10.55)经过化简后,可表示为

\(\begin{align}& \rho \left( \frac{\partial \overline{{{V}_{x}}}}{\partial t}+\frac{\partial \overline{{{V}_{x}}}\overline{{{V}_{x}}}}{\partial x}+\frac{\partial \overline{{{V}_{x}}}\overline{{{V}_{y}}}}{\partial y}+\frac{\partial \overline{{{V}_{x}}}\overline{{{V}_{z}}}}{\partial z} \right)= \\& -\frac{\partial \overline{p}}{\partial x}+\mu \left( \frac{{{\partial }^{2}}\overline{{{V}_{x}}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}\overline{{{V}_{x}}}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}\overline{{{V}_{x}}}}{\partial {{z}^{2}}} \right)+\frac{\partial \left( -\rho \overline{V_{x}^{\prime }V_{x}^{\prime }} \right)}{\partial x}+\frac{\partial \left( -\rho \overline{V_{x}^{\prime }V_{y}^{\prime }} \right)}{\partial y}+\frac{\partial \left( -\rho \overline{V_{x}^{\prime }V_{z}^{\prime }} \right)}{\partial z} \\\end{align}\)

再应用时均运动的连续方程式(10.53a),上式可化为

\(\left. \begin{align}& \rho \left( \frac{\partial \overline{{{V}_{x}}}}{\partial t}+\overline{{{V}_{x}}}\frac{\partial \overline{{{V}_{x}}}}{\partial x}+\overline{{{V}_{y}}}\frac{\partial \overline{{{V}_{x}}}}{\partial y}+\overline{{{V}_{z}}}\frac{\partial \overline{{{V}_{x}}}}{\partial z} \right)=-\frac{\partial \overline{p}}{\partial x}+\mu \left( \frac{{{\partial }^{2}}\overline{{{V}_{x}}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}\overline{{{V}_{x}}}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}\overline{{{V}_{x}}}}{\partial {{z}^{2}}} \right)+ \\& \frac{\partial \left( -\rho \overline{V_{x}^{\prime }V_{x}^{\prime }} \right)}{\partial x}+\frac{\partial \left( -\rho \overline{V_{x}^{\prime }V_{y}^{\prime }} \right)}{\partial y}+\frac{\partial \left( -\rho \overline{V_{x}^{\prime }V_{z}^{\prime }} \right)}{\partial z} \\&  \\& \rho \left( \frac{\partial \overline{{{V}_{y}}}}{\partial t}+\overline{{{V}_{x}}}\frac{\partial \overline{{{V}_{y}}}}{\partial x}+\overline{{{V}_{y}}}\frac{\partial \overline{{{V}_{y}}}}{\partial y}+\overline{{{V}_{z}}}\frac{\partial \overline{{{V}_{y}}}}{\partial z} \right)=-\frac{\partial \overline{p}}{\partial y}+\mu \left( \frac{{{\partial }^{2}}\overline{{{V}_{y}}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}\overline{{{V}_{y}}}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}\overline{{{V}_{y}}}}{\partial {{z}^{2}}} \right)+ \\& \frac{\partial \left( -\rho \overline{V_{x}^{\prime }V_{y}^{\prime }} \right)}{\partial x}+\frac{\partial \left( -\rho \overline{V_{y}^{\prime }V_{y}^{\prime }} \right)}{\partial y}+\frac{\partial \left( -\rho \overline{V_{y}^{\prime }V_{z}^{\prime }} \right)}{\partial z} \\& \rho \left( \frac{\partial \overline{{{V}_{z}}}}{\partial t}+\overline{{{V}_{x}}}\frac{\partial \overline{{{V}_{z}}}}{\partial x}+\overline{{{V}_{y}}}\frac{\partial \overline{{{V}_{z}}}}{\partial y}+\overline{{{V}_{z}}}\frac{\partial \overline{{{V}_{z}}}}{\partial z} \right)=-\frac{\partial \overline{p}}{\partial z}+\mu \left( \frac{{{\partial }^{2}}\overline{{{V}_{z}}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}\overline{{{V}_{z}}}}{\partial {{y}^{2}}}+\frac{{{\partial }^{2}}\overline{{{V}_{z}}}}{\partial {{z}^{2}}} \right)+ \\& \frac{\partial \left( -\rho \overline{V_{x}^{\prime }V_{z}^{\prime }} \right)}{\partial x}+\frac{\partial \left( -\rho \overline{V_{y}^{\prime }V_{z}^{\prime }} \right)}{\partial y}+\frac{\partial \left( -\rho \overline{V_{z}^{\prime }V_{z}^{\prime }} \right)}{\partial z} \\\end{align} \right\}\)               (10.56)

方程组式(10.56)就是著名的不可压缩流体作湍流运动时的时均运动方程,称为雷诺方程。

将时均运动方程式(10.56)和N一S方程式(10.54a)相比可以看出,湍流中的应力,除了由于黏性所产生的应力外,还有由于湍流脉动运动所形成的附加应力,这些附加应力称为雷诺应力。雷诺方程与N一S方程在形式上是相同的,只不过在黏性应力项中多出了附加的湍流应力项。

以上导出的雷诺方程和连续方程中,除过要求解的4个变量Vx,Vy,Vz和p外,还有与脉动速度有关的如\(\overline{V_{x}^{\prime }V_{x}^{\prime }},\overline{V_{x}^{\prime }V_{y}^{\prime }}\)等6个未知数。4个方程中有10个未知数,即方程组不封闭。要使方程组封闭,必须补充其他未知量的关系式才能够进行求解。

10.4.4  雷诺应力

将雷诺方程与黏性流体应力形式的动量方程进行比较,由式( 10.56)可以看出,在湍流的时均运动方程中,除了原有的黏性应力分量外,还多出了由脉动速度乘积的时均值 \(-\rho \overline{V_{x}^{\prime }V_{x}^{\prime }},-\rho \overline{V_{x}^{\prime }V_{y}^{\prime }}\)等构成的附加项,这些附加项构成了一个对称的二阶张量G,即

\(G=\left[ \begin{align}& -\rho \overline{V_{x}^{\prime }V_{x}^{\prime }},-\rho \overline{V_{x}^{\prime }V_{y}^{\prime }},-\rho \overline{V_{x}^{\prime }V_{z}^{\prime }} \\& -\rho \overline{V_{x}^{\prime }V_{y}^{\prime }},-\rho \overline{V_{y}^{\prime }V_{y}^{\prime }},-\rho \overline{V_{y}^{\prime }V_{z}^{\prime }} \\& -\rho \overline{V_{x}^{\prime }V_{z}^{\prime }},-\rho \overline{V_{y}^{\prime }V_{z}^{\prime }},-\rho \overline{V_{z}^{\prime }V_{z}^{\prime }} \\\end{align} \right]\)               (10.57)

式(10.57)中的各项构成了所谓的雷诺应力。雷诺应力的物理意义可理解如下。

在稳定湍流中的某点M处取一微元六面体(见图10.5(a)),考察过点M取与x方向垂直的某微元面,其面积为dA1。在单位时间内通过单位面积的动量为\(\rho V_{x}^{2}\),其时均值为

\(\overline{\rho V_{x}^{2}}=\rho {{\overline{{{V}_{x}}}}^{2}}+\rho \overline{V_{x}^{\prime }V_{x}^{\prime }}\)               (10.58)

式(10.58)等号左端是单位时间内通过垂直于x方向的单位面积所传递的真实动量的平均值;等号右端第一项是同一时间内通过同一面积所传递的按时均速度计算的动量,第二项是由于x方向上速度脉动所传递的动量。根据动量定理,通过dA1面有动量传递,那么在dA1面上就有力的作用。式(10.58)中各项都具有力的量纲,从而证明了在湍流情况下,沿x方向的时均真实应力,应等于时均运动情况下x方向上的应力加上由于湍流中的x方向脉动引起的附加应力。对dA1面来说,附加应力\(\rho \overline{V_{x}^{\prime }V_{x}^{\prime }}\)与它垂直,所以是法向应力,因此称之为附加湍流正应力。

由于在点M处沿y方向上有脉动速度V’y,则在单位时间内通过微元面dA2(垂直于y方向)上的单位面积流入的质量为ρV’y,如图10.5(a)所示。这部分流体本身具有x方向的速度\({{V}_{x}}=\overline{{{V}_{x}}}+V_{x}^{\prime }\),因而随之传递的x方向上的动量为ρV’yVx,其时均值为

\(\overline{\rho V_{y}^{\prime }{{V}_{x}}}=\overline{\rho V_{y}^{\prime }\left( \overline{{{V}_{x}}}+V_{x}^{\prime } \right)}=\overline{\rho V_{y}^{\prime }\overline{{{V}_{x}}}}+\overline{\rho V_{y}^{\prime }V_{x}^{\prime }}\)

根据时均运算关系式,\(\overline{\rho V_{y}^{\prime }\overline{{{V}_{x}}}}=0\),所以

\(\overline{\rho V_{y}^{\prime }{{V}_{x}}}=\rho \overline{V_{y}^{\prime }V_{x}^{\prime }}\)               (10.59a)

图10.5(b)表示一个单位长度的流体微团因y方向的速度脉动V′y,而在单位时间内通过单位面积上增加的x方向上的动量的时均值,即

\(\overline{\rho V_{y}^{\prime }V_{x}^{\prime }}-\left( \overline{\rho V_{y}^{\prime }V_{x}^{\prime }+\frac{\partial \left( \rho V_{y}^{\prime }V_{x}^{\prime } \right)}{\partial y}\times 1} \right)=\frac{\partial \left( -\rho \overline{V_{y}^{\prime }V_{x}^{\prime }} \right)}{\partial y}\)               (10.59b)

式(10.59a)表明,在单位时间内通过垂直于y方向的dA2面的单位面积所传递出去的x方向动量为\(\rho \overline{V_{y}^{\prime }V_{x}^{\prime }}\),因而该单位面积就受到一个沿x方向的、大小为\(\rho \overline{V_{y}^{\prime }V_{x}^{\prime }}\)的作用力。式(10.59b)说明了这个力的变化量。可以理解为:当流体质点由时均速度较高的流体层向时均速度较低的流体层脉动时由于脉动引起的动量传递,使低速层被加速;反过来,如果脉动由低速层向高速层发生,高速层被减速,因此这两层流体在x方向上各受到切应力的作用。\(\rho \overline{V_{y}^{\prime }V_{x}^{\prime }}\)是湍流中流体微团的脉动造成的,称为附加湍流切应力。

湍流正应力和湍流切应力统称为雷诺应力。

10.4.5  普朗特混合长度理论

从雷诺方程可以看出,由于湍流运动采用了时均方法,在运动方程中出现了雷诺应力,从而增加了方程中的未知量,因此需要补充新的关系式才能求解。如果补充的关系式是一个代数方程,而不需要补充任何附加的微分方程来求解时均流场,则称这种模型为零方程模型;若补充的关系式是一个微分方程(如湍流脉动动能方程),则称为一方程模型;若是两个微分方程,则称为双方程模型;等等。本节所讨论的普朗特混合长度理论即是所谓的代数模型(零方程模型)。

混合长度理论是基于经验性的一个经过实验验证的理论模型。在许多问题中得到了较好的应用。其基本思想是如果能够找出湍流应力与其他流场参数之间的关系,即找到了这些物理量的补充关系式,就可以使方程组封闭。为此,普朗特把湍流脉动与气体分子运动相比拟,认为雷诺应力是由流体微团的脉动引起的,它和分子运动引起黏性应力的情况十分相似。在定常层流直线运动中,由分子动量输运而引起的黏性切应力\({{\tau }_{t}}=\mu d{{V}_{x}}/dy\),与此相对应,当湍流时均流动的流线为直线时,认为脉动引起的雷诺切应力(湍流应力)也可以表示成上述类似的形式,即

\({{\tau }_{t}}={{\mu }_{t}}\frac{d\overline{{{V}_{x}}}}{dy}\)               (10.60)

式中,μt称为湍流黏度。这就是混合长度理论的基本思想。

另外,湍流应力与脉动速度有关。为了确定这种关系,普朗特做出了第一个假设,即流体微团x方向脉动速度V′x近似等于两层流体的时均速度之差,即

\(V_{x}^{\prime }\approx \left( \overline{{{V}_{x}}}+l\frac{d\overline{{{V}_{x}}}}{dy} \right)-\overline{{{V}_{x}}}=l\frac{d\overline{{{V}_{x}}}}{dy}\)

这一假设的基础是认为流体微团在y方向脉动,从这一层跳入另一层时,要经过一段与其他流体微团不相碰撞的距离I(见图10.6),在这段距离上速度保持不变。这个距离l称为混合长度,它是流体微团在湍流运动中的自由行程的平均值。经过l距离后,流体微团以自己原来的动量进入另一层和周围流体相掺混。

从图10.6上可以看出,(y+l)层上的流体质点脉动到y层时,其速度比y层上的流体时均速度大\(l\frac{d\overline{{{V}_{x}}}}{dy}\)。它引起y层上流体速度有一个正的脉动,其值为\(V_{x}^{\prime }=l\frac{d\overline{{{V}_{x}}}}{dy}\) 。同理,当流体微团从y层脉动到(y+l)层时,使(y+l)层的流体有一个负的脉动速度,其大小也是\(l\frac{d\overline{{{V}_{x}}}}{dy}\)。

普朗特又做出第二个假设。他认为y方向的脉动速度V′y和V′x成正比。其根据可用图10.6说明。两层流体混合时,由于上、下两层流体的速度差为\(l\frac{d\overline{{{V}_{x}}}}{dy}\),且两层流体的质点间相互作用,从而引起横向脉动,其速度为V′y,因此该速度必然与速度差\(l\frac{d\overline{{{V}_{x}}}}{dy}\)有关,且两者具有相同的数量级。考虑到V′y与V′x的符号相反,因此有

\(-V_{y}^{\prime }=CV_{x}^{\prime }\approx Cl\frac{d\overline{{{V}_{x}}}}{dy}\)

普朗特引入了混合长度的概念,确定了脉动速度V′x,V′y的大小与时均速度梯度之间的关系,从而确定湍流切应力的大小,即

\({{\tau }_{\text{t}}}=-\rho \overline{V_{x}^{\prime }V_{y}^{\prime }}=-\rho \frac{1}{T}\int_{0}^{T}{V_{x}^{\prime }V_{y}^{\prime }dt=}\rho \frac{1}{T}\int_{0}^{T}{C{{\left( l\frac{d\overline{{{V}_{x}}}}{dy} \right)}^{2}}dt=\rho C{{l}^{2}}}{{\left( \frac{d\overline{{{V}_{x}}}}{dy} \right)}^{2}}\)

式中,混合长度l尚未确定,因此可取C=1。这样湍流切应力就可以写为

\({{\tau }_{\text{t}}}=-\rho \overline{V_{x}^{\prime }V_{y}^{\prime }}=\rho {{l}^{2}}{{\left( \frac{d\overline{{{V}_{x}}}}{dy} \right)}^{2}}\)

考虑到湍流切应力τt的符号应与黏性切应力τt的符号相同,为标出符号,上式可写成

\({{\tau }_{\text{t}}}=\rho {{l}^{2}}\left| \frac{d\overline{{{V}_{x}}}}{dy} \right|\frac{d\overline{{{V}_{x}}}}{dy}={{\mu }_{t}}\frac{d\overline{{{V}_{x}}}}{dy}\)               (10.61)

式中,\({{\mu }_{t}}=\rho {{l}^{2}}\left| \frac{d\overline{{{V}_{x}}}}{dy} \right|\),混合长度l一般需要实验确定。

10.5  附面层基本知识

10.5.1  附面层的概念

1.附面层厚度及流动阻力

黏性是流体的重要属性。根据流体黏性的特点,在靠近物体表面处,流体将黏附在物面上而流速为零,即满足无滑移条件。而沿物而的法线方向上,流速逐渐增加,到某一距离处,流速与外边界速度近似相等。通常定义靠近物体表面,存在较大速度梯度的簿层称为附面层或边界层。当V=0.99V0(V0为附面层外边界的速度)时的垂直物面的法向距离为附面层厚度,用δ表示。在航空上,有实际意义的问题大多属于大雷诺数下的流动问题,且靠近物体表面速度梯度很大的这一层都是很薄的,因此附面层厚度δ是个小量。气流流过物体表面的距离越长,附面层厚度也越大,即附面层厚度随气流流过物体的距离而增加。黏性影响较大的另一种情况是流体在物体后面形成所谓的尾迹区。由于附面层和尾迹区内黏性的作用较强,黏性切应力作用较大,因而形成流动阻力。显然,该阻力产生的根源是流体与物体表面之间的摩擦以及附面层分离引起的。之外,在由于附面层脱离后形成的尾迹区中,还会导致物体表面上产生流动方向的压力差,因而形成所谓的压差阻力。

在附而层外边界,流速接近于外边界速度,因此附面层外边界的速度梯度很小。而空气的黏度系数也很小,所以在附面层之外,可以忽略黏性的影响,而作为理想流动来处理。总之,在靠近物体表面的附面层内以及在物体之后的尾迹区内,黏性都有显著的影响。

2.附面层中沿物面的法向压强保持近似不变

在附面层内,除了速度梯度∂V/∂y很大外,还有另外一个重要的特点,对于物面曲率半径比较大,即物而不太弯曲的情况,沿着其物面的法线方向流体压强保持近似不变。如果测量流体流过平板的附面层内沿y方向的压强梯度,则可以得到在附面层内压强p沿y方向不变,即∂p/∂y=0。该结论非常重要,它可以使附面层运动方程大大简化;同时,它还使得理想流体的结论具有实际意义。按理想流体理论计算附面层外边界的压强分布后,即可得到物面上对应点的压强。

3.位移厚度δ*(流量损失厚度或排移厚度)

所谓的位移厚度δ*就是由于附面层内速度降低(流量有损失)而要求流道加宽的厚度。

设物体上某点处的附面层厚度为δ,如图10.7所示,垂直纸面方向为单位宽度,则由于附面层内的流动速度减小,使得实际流过附面层内的质量流量比没有附面层(理想流体时)减少了,所减少的质量流量为

\(\int_{0}^{\delta }{\left( {{\rho }_{0}}{{V}_{0}}-\rho {{V}_{x}} \right)}dy\)

式中,ρ0,V0分别是附面层外边界的理想流体的密度和速度;ρ,Vx分别是附面层内流体的密度和速度。这些减少的质量流量要在主流中挤出δ*的距离才能流过去。因此,它应等于以理想流体(ρ0,V0)流过δ*距离上的质量流量,即

\[{{\rho }_{0}}{{V}_{0}}{{\delta }^{*}}=\int_{0}^{\delta }{\left( {{\rho }_{0}}{{V}_{0}}-\rho {{V}_{x}} \right)}dy\]

所以,得

\({{\delta }^{*}}=\int_{0}^{\delta }{\left( 1-\frac{\rho {{V}_{x}}}{{{\rho }_{0}}{{V}_{0}}} \right)}dy\)               (10.62)

由此可见,在质量流量相等的条件下,犹如将理想流体的流动区域自物面向外移动了一个δ*的距离。它表示了由于黏性的作用,附面层内流体质量流量相对理想流体减小的程度。

对于不可压缩流体,式(10.62)可改写为

\({{\delta }^{*}}=\int_{0}^{\delta }{\left( 1-\frac{{{V}_{x}}}{{{V}_{0}}} \right)}dy\)               (10.63)

根据以上的分析,如果按理想流体设计的型面,为了使相同质量流量的黏性流体能够通过,则物面应向外移动一个δ*的距离。

位移厚度的概念,对于流动方向要求严格的流道设计具有重要的意义。特别是对于管道内出现声速截面时,附面层会使实际流过通道的流量减少。这时,应对按理想流体设计出来的管道壁面进行修止。

例如,设计喷管时,通常先把喷管中的流动看成是无黏性的,求喷管的理想型线,然后考虑黏性的影响,把理想型线上的各点都增加当地位移厚度δ*。

4.动量损失厚度δ**

由于附面层内的流速小于理想流体的流速,因此附面层内流体的动量也会减小。通过附面层厚度δ的流体实际具有的动量为\(\int_{0}^{\delta }{\rho V_{x}^{2}}dy\),此部分流体若以附面层外边界上理想流体速度V0运动时,所具有的动量为\(\int_{0}^{\delta }{\rho {{V}_{x}}{{V}_{0}}}dy\)。因此其动量损失应等于单位时间内以速度V0、密度ρ0流过一面积为δ**(厚度)× 1(宽度)的流体所具有的动量,即

\({{\rho }_{0}}V_{0}^{2}{{\delta }^{**}}=\int_{0}^{\delta }{\rho {{V}_{x}}{{V}_{0}}}dy-\int_{0}^{\delta }{\rho V_{x}^{2}}dy\)

δ**称为动量损失厚度,即

\({{\delta }^{**}}=\int_{0}^{\delta }{\frac{\rho {{V}_{x}}}{{{\rho }_{0}}{{V}_{0}}}\left( 1-\frac{{{V}_{x}}}{{{V}_{0}}} \right)}dy\)               (10.64)

对不可压缩流体,ρ0=ρ,则

\({{\delta }^{**}}=\int_{0}^{\delta }{\frac{{{V}_{x}}}{{{V}_{0}}}\left( 1-\frac{{{V}_{x}}}{{{V}_{0}}} \right)}dy\)               (10.65)

由式(10.63)和式(10.65)可以看出,位移厚度和动量损失厚度与附面层内的速度分布和附面层厚度有关,通常用比值H=δ*/δ**表示附面层内速度分布形状的参数,H称为形状因子。H越大,说明附面层内的速度分布越呈现凹形状,如图10.8(a)所示;H越小,速度分布越饱满,如图10.8(b)所示。因此可见,湍流附面层的形状因子比层流的小,这是由于湍流流动中流体质点横向脉动的动量变换,使附面层内的速度分布更加均匀。工程中常用形状因子来判定附面层是否出现分离。一般认为,当H ≥3.5时,层流附面层会发生分离。而当H ≥ 1.4~1.75时,湍流附面层会发生分离。

10.5.2  附面层的转捩

根据雷诺实验,黏性流体存在着两种流态,即层流和湍流。附面层流动和管流一样有层流附面层和湍流附面层之分。实验观察表明,流体从物体前缘开始,先形成层流附面层。层流附面层的存在有一个极限情况,超过此极限时,层流处于不稳定状态,并逐渐过渡为湍流附面层。图10.9所示是均匀来流流过平板时的流动图形。图中,O—A称为层流附面层段,A—B称为转捩段,转捩起点A距平板前缘的距离用XT表示,对应于转捩点A的雷诺数称为临界雷诺数,即\({{\operatorname{Re}}_{cr}}=\frac{{{V}_{\infty }}{{V}_{T}}}{\upsilon }\),通常转捩雷诺数的大小要由实验确定。一般地,对于绕平板的流动,Recr= 5 ×105 ~3 × 106。经过转捩段A—B 后,即Re> Recr,附面层转变为湍流。由Recr可以得到转捩点的位置,即

\({{X}_{T}}=\frac{\mu {{\operatorname{Re}}_{cr}}}{\rho {{V}_{\infty }}}\)               (10.66)

由式(10.66)可见,转捩点的位置与流体的黏度、密度、来流速度和临界雷诺数有关。

参考文献[ 5 ]引用了米歇尔(Michel)基于实验提出的转捩点位置XT处的雷诺数和相应的动量损失厚度为参考长度的雷诺数之间的关系为

\({{\operatorname{Re}}_{{{\delta }^{**}}}}\approx 2.9{{\operatorname{Re}}_{{{X}_{T}}}}\)               (10.67a)

参考文献[ 6 ]给出了经过改进的半经验公式为

\(\begin{align}& {{\operatorname{Re}}_{{{\delta }^{**}}}}\approx 1.718\operatorname{Re}_{{{X}_{T}}}^{0.435} \\& 0.3R{{e}_{{{X}_{T}}}}\times {{10}^{-6}}20 \\\end{align}\)               (10.67b)

只要速度分布光滑和表面光滑,式(10.67)即是确定转捩点位置较好的方法。

10.6  附面层微分方程

附面层概念的提出,可以将黏性流动的求解简化为求解附面层内的流动和附面层外的理想流动;可以根据附面层的特点对N—S方程进行简化,得到求解附面层的微分方程。

10.6.1  层流附面层微分方程

为了简化推导,考虑平壁面二维不可压层流流动,取平行于物面的方向为x方向,垂直于物面的为y方向。如果忽略壁面曲率和质量力的影响,则连续方程和N—S方程可表示为

\(\left. \begin{align}& \frac{\partial {{V}_{x}}}{\partial x}+\frac{\partial {{V}_{y}}}{\partial y}=0 \\& \frac{\partial {{V}_{x}}}{\partial t}+{{V}_{x}}\frac{\partial {{V}_{x}}}{\partial x}+{{V}_{y}}\frac{\partial {{V}_{x}}}{\partial y}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\upsilon \left( \frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{y}^{2}}} \right) \\& \frac{\partial {{V}_{y}}}{\partial t}+{{V}_{x}}\frac{\partial {{V}_{y}}}{\partial x}+{{V}_{y}}\frac{\partial {{V}_{y}}}{\partial y}=-\frac{1}{\rho }\frac{\partial p}{\partial y}+\upsilon \left( \frac{{{\partial }^{2}}{{V}_{y}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}{{V}_{y}}}{\partial {{y}^{2}}} \right) \\\end{align} \right\}\)               (10.68)

为了简化式(10.68),对它进行无量纲化。根据附面层流动的特点,选取附面层外边界速度V0、物体的特征长度L、附面层厚度δ及密度ρ为特征量,对式(10.68)进行无量纲化,即令

\(\left. \begin{align}& \overline{x}=\frac{x}{L},\overline{y}=\frac{y}{L}=\frac{y}{L/\sqrt{\operatorname{Re}}},\overline{t}=\frac{t}{L/{{V}_{0}}} \\& \overline{{{V}_{x}}}=\frac{{{V}_{x}}}{{{V}_{0}}},\overline{p}=\frac{p}{\rho V_{0}^{2}},\overline{{{V}_{y}}}=\frac{{{V}_{y}}}{{{V}_{0}}/\sqrt{\operatorname{Re}}} \\\end{align} \right\}\)               (10.69)

式中,\(\operatorname{Re}=\frac{\rho {{V}_{0}}L}{\mu }\)。将式(10.69)代入基本方程式(10.68),可得

\(\left. \begin{align}& \frac{\partial \overline{{{V}_{x}}}}{\partial \overline{x}}+\frac{\partial \overline{{{V}_{y}}}}{\partial \overline{y}}=0 \\& \frac{\partial \overline{{{V}_{x}}}}{\partial \overline{t}}+\overline{{{V}_{x}}}\frac{\partial \overline{{{V}_{x}}}}{\partial \overline{x}}+\overline{{{V}_{y}}}\frac{\partial \overline{{{V}_{x}}}}{\partial \overline{y}}=-\frac{\partial \overline{p}}{\partial \overline{x}}+\frac{1}{\operatorname{Re}}\frac{{{\partial }^{2}}\overline{{{V}_{x}}}}{\partial {{\overline{x}}^{2}}}+\frac{{{\partial }^{2}}\overline{{{V}_{x}}}}{\partial {{\overline{y}}^{2}}} \\& \frac{1}{\operatorname{Re}}\left( \frac{\partial \overline{{{V}_{y}}}}{\partial \overline{t}}+\overline{{{V}_{x}}}\frac{\partial \overline{{{V}_{y}}}}{\partial \overline{x}}+\overline{{{V}_{y}}}\frac{\partial \overline{{{V}_{y}}}}{\partial \overline{y}} \right)=-\frac{\partial \overline{p}}{\partial \overline{y}}+\frac{1}{{{\operatorname{Re}}^{2}}}\frac{{{\partial }^{2}}\overline{{{V}_{y}}}}{\partial {{\overline{x}}^{2}}}+\frac{1}{\operatorname{Re}}\frac{{{\partial }^{2}}\overline{{{V}_{y}}}}{\partial {{\overline{y}}^{2}}} \\\end{align} \right\}\)               (10.70)

式中,符号上带“—’的物理量的数量级均为1,因此各项的量级取决于相应的系数的量级。由于在附面层中Re≥1,所以方程中带有1/Re²,1/Re系数的项可以忽略。方程变为

\(\left. \begin{align}& \frac{\partial \overline{{{V}_{x}}}}{\partial \overline{x}}+\frac{\partial \overline{{{V}_{y}}}}{\partial \overline{y}}=0 \\& \frac{\partial \overline{{{V}_{x}}}}{\partial \overline{t}}+\overline{{{V}_{x}}}\frac{\partial \overline{{{V}_{x}}}}{\partial \overline{x}}+\overline{{{V}_{y}}}\frac{\partial \overline{{{V}_{x}}}}{\partial \overline{y}}=-\frac{\partial \overline{p}}{\partial \overline{x}}+\frac{{{\partial }^{2}}\overline{{{V}_{x}}}}{\partial {{\overline{y}}^{2}}} \\& \frac{\partial \overline{p}}{\partial \overline{y}}=0 \\\end{align} \right\}\)               (10.71)

利用式(10.69),可将式(10.71)还原为有量纲形式的方程,即

\(\left. \begin{align}& \frac{\partial {{V}_{x}}}{\partial x}+\frac{\partial {{V}_{y}}}{\partial y}=0 \\& \frac{\partial {{V}_{x}}}{\partial t}+{{V}_{x}}\frac{\partial {{V}_{x}}}{\partial x}+{{V}_{y}}\frac{\partial {{V}_{x}}}{\partial y}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\upsilon \frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{y}^{2}}} \\& \frac{\partial p}{\partial y}=0 \\\end{align} \right\}\)               (10.72)

式(10.72)即为平面壁的二维不可压层流附面层方程。由式(10.72)的最后一个方程可以看出,对于直壁,沿垂直于壁面方向,压强近似保持不变,即附面层内横向截面上的压强近似等于附面层外边界处的主流压强。因此,在求解绕平面物体(或物面曲率半径比较大的物体)的流动时,第三个方程可以去掉,而压强可以用附面层外边界的压强p0代替。因此,平面壁的二维不可压附面层方程为

\(\frac{\partial {{V}_{x}}}{\partial x}+\frac{\partial {{V}_{y}}}{\partial y}=0\)               (10.73a)

\(\frac{\partial {{V}_{x}}}{\partial t}+{{V}_{x}}\frac{\partial {{V}_{x}}}{\partial x}+{{V}_{y}}\frac{\partial {{V}_{x}}}{\partial y}=-\frac{1}{\rho }\frac{\partial {{p}_{0}}}{\partial x}+\upsilon \frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{y}^{2}}}\)               (10.73b)

对于曲面物体,采用沿曲面壁方向作为x方向,y方向与x方向垂直并从壁面算起,采用正交曲线坐标系,并采用与上述同样的分析方法。考虑到物面的曲率半径为r,经数量级分析后,得到曲线坐标系中的附面层方程为

\(\left. \begin{align}& \frac{\partial {{V}_{x}}}{\partial x}+\frac{\partial {{V}_{y}}}{\partial y}=0 \\& \frac{\partial {{V}_{x}}}{\partial t}+{{V}_{x}}\frac{\partial {{V}_{x}}}{\partial x}+{{V}_{y}}\frac{\partial {{V}_{x}}}{\partial y}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\upsilon \frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{y}^{2}}} \\& \frac{1}{\rho }\frac{\partial p}{\partial y}=\frac{V_{x}^{2}}{r} \\\end{align} \right\}\)               (10.74)

由式(10.74)可以看出,对于曲壁的情况,由于壁面弯曲产生的离心力,使得横向的压强梯度不为零,显然这是由于壁面弯曲造成的。

求解附面层方程式(10.73)或式(10.74),必须根据具体问题提出相应的边界条件和初始条件。下面给出初始条件和附面层内、外边界上的边界条件。

初始条件:

当t=t0时,

Vx=Vx(x,y,t0),Vy=Vy(x,y,t0)

边界条件:

(1)在物面上,满足无滑移条件,即当y=0时,Vx=0 ,Vy=0;

(2)在附面层外边界,满足外边界条件,即当y=δ时,Vx=V0(x),V0(x)是附面层外边界上的理想流体的速度,可以通过附面层外的无黏性流动求出。

10.6.2  湍流附面层微分方程

对于二维不可压湍流附面层,N—S方程式(10.68)中的动量方程中存在有湍流切应力的附加应力项,对二维流动应用式(10.53a)和式(10.56),省略各项时均化参数的记号,则有

\(\left. \begin{align}& \frac{\partial {{V}_{x}}}{\partial x}+\frac{\partial {{V}_{y}}}{\partial y}=0 \\& \frac{\partial {{V}_{x}}}{\partial t}+{{V}_{x}}\frac{\partial {{V}_{x}}}{\partial x}+{{V}_{y}}\frac{\partial {{V}_{x}}}{\partial y}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\upsilon \left( \frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{y}^{2}}} \right)- \\& \left( \frac{\partial V_{x}^{\prime }V_{x}^{\prime }}{\partial x}+\frac{\partial V_{x}^{\prime }V_{y}^{\prime }}{\partial y} \right) \\& \frac{\partial {{V}_{y}}}{\partial t}+{{V}_{x}}\frac{\partial {{V}_{y}}}{\partial x}+{{V}_{y}}\frac{\partial {{V}_{y}}}{\partial y}=-\frac{1}{\rho }\frac{\partial p}{\partial y}+\upsilon \left( \frac{{{\partial }^{2}}{{V}_{y}}}{\partial {{x}^{2}}}+\frac{{{\partial }^{2}}{{V}_{y}}}{\partial {{y}^{2}}} \right)- \\& \left( \frac{\partial V_{y}^{\prime }V_{x}^{\prime }}{\partial x}+\frac{\partial V_{y}^{\prime }V_{y}^{\prime }}{\partial y} \right) \\\end{align} \right\}\)               (10.75)

经过数量级的分析,湍流附面层方程可以与成如下形式:

\(\left. \begin{align}& \frac{\partial {{V}_{x}}}{\partial x}+\frac{\partial {{V}_{y}}}{\partial y}=0 \\& \frac{\partial {{V}_{x}}}{\partial t}+{{V}_{x}}\frac{\partial {{V}_{x}}}{\partial x}+{{V}_{y}}\frac{\partial {{V}_{x}}}{\partial y}=-\frac{1}{\rho }\frac{\partial p}{\partial x}+\upsilon \frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{y}^{2}}}-\left( \frac{\partial V_{x}^{\prime }V_{x}^{\prime }}{\partial x}+\frac{\partial V_{x}^{\prime }V_{y}^{\prime }}{\partial y} \right) \\& \frac{1}{\rho }\frac{\partial p}{\partial y}=\frac{V_{x}^{2}}{r} \\\end{align} \right\}\)               (10.76)

式中,各物理量指的是时均值,如Vx指\(\overline{{{V}_{x}}}\),V′xV′y指\(\overline{V_{x}^{\prime }V_{y}^{\prime }}\)。

求解方程组式(10.76)的边界条件与式(10.74)相同。

10.7  附面层积分方程

虽然附面层微分方程比较N—S方程有了很大的简化,但是要求解这一组偏微分方程,其计算工作量仍然很大,需要借助于计算机进行数值求解。求解附面层问题的另一种方法是附面层积分法。这种方法的基本思想是使流动参数在总体上满足附面层基本方程。在求解时,近似地给定一个只依赖于x坐标的单参数速度分布来代替附面层内真实的速度分布。解法的精确度取决于所选定的速度分布。

10.7.1  附面层的动量积分方程

附面层积分方程可以由两种方法导出。一种是将附面层微分方程在整个附面层厚度δ的区间上积分,另一种是在附面层内取一微元段,运用基本方程。前者主要是从数学上推导,而后者的物理概念比较清楚。下面采用后一种推导方法来得出附面层动量积分方程。

图10.10表示了附面层内流体沿某一壁面的流动。设流动为定常的二维平面不可压流动。在附面层中取一微元控制体ABDCA,其中AB和CD是垂直于壁面的两个控制面,相距为dx,BD是壁面,AC是附面层外边界。垂直于纸面控制体的宽度取单位宽度。对控制体运用动量定理。

表10.3给出了通过各控制面上的质量流量和相应的动量。

需要强调一点,由于dx是无限小量,所以将AC边界上的流体速度都看做是V0,实际上,V0是x的函数,V0(x)由壁面形状决定。

在单位时间内,通过控制面流出与流入控制体的动量的差值为

\(\frac{\partial }{\partial x}\left( \int_{0}^{\delta }{\rho V_{x}^{2}dy} \right)dx-{{V}_{0}}\frac{\partial }{\partial x}\left( \int_{0}^{\delta }{\rho {{V}_{x}}dy} \right)dx\)

进一步分析作用在控制体上的力。因为对平板在附面层内∂p/∂y=0,所以在AB,CD面上的压强沿y方向没有变化。于是沿x方向作用在控制体上的力,见表10.4。

在表10.4中,AC面上的压强取点A和点C压强的平均值。AC面积在x方向的投影面积大小为dδ x 1=dδdx/dx。-τw表示壁面上的摩擦应力。在CD和BD面上的作用力方向与x方向相反,所以都带有负号。

作用在控制体上沿x方向上的合力,经过化简整理后,得

\(p\delta -\left[ p\delta +\frac{\partial \left( p\delta  \right)}{\partial x}dx \right]+\left( p+\frac{\partial p}{\partial x}\frac{dx}{2} \right)\frac{d\delta }{dx}dx-{{\tau }_{w}}dx\approx -\left( {{\tau }_{w}}+\frac{\partial p}{\partial x}\delta  \right)dx\)

根据动量定理,作用在控制体上所有作用力的合力等于单位时间流出和流入控制体动量之差,即

\(-\left( {{\tau }_{w}}+\frac{\partial p}{\partial x}\delta  \right)dx=\frac{\partial }{\partial x}\left( \int_{0}^{\delta }{\rho V_{x}^{2}dy} \right)dx-{{V}_{0}}\frac{\partial }{\partial x}\left( \int_{0}^{\delta }{\rho {{V}_{x}}dy} \right)dx\)

\(-{{\tau }_{w}}-\frac{\partial p}{\partial x}\delta =\frac{\partial }{\partial x}\left( \int_{0}^{\delta }{\rho V_{x}^{2}dy} \right)-{{V}_{0}}\frac{\partial }{\partial x}\left( \int_{0}^{\delta }{\rho {{V}_{x}}dy} \right)\)               (10.77)

式(10.77)称为附面层积分方程。该方程对于层流附面层和湍流附面层都适用。对于后一种情况,可直接将附面层连续方程和动量方程相加后沿附面层积分得到,积分时注意到在壁面上及附面层外边界处湍流应力等于零。

对不可压缩流体,式(10.77)化为

\(\frac{d}{dx}\left( \int_{0}^{\delta }{V_{x}^{2}dy} \right)-{{V}_{0}}\frac{d}{dx}\left( \int_{0}^{\delta }{{{V}_{x}}dy} \right)=-\frac{\delta }{\rho }\frac{dp}{dx}-\frac{{{\tau }_{w}}}{\rho }\)               (10.78)

式(10.78)右端的压强梯度可以根据附面层外边界的理想流动得出。根据伯努利方程

\({{p}_{0}}+\frac{1}{2}\rho V_{0}^{2}=\)常数

对x求导后得                                                   \(\frac{d{{p}_{0}}}{dx}=-\rho {{V}_{0}}\frac{d{{V}_{0}}}{dx}\)

注意到\(\delta =\int_{0}^{\delta }{dy}\),则式(10.78)等号右端第一项写为

\(-\frac{\delta }{\rho }\frac{dp}{dx}={{V}_{0}}\frac{d{{V}_{0}}}{dx}\delta ={{V}_{0}}\frac{d{{V}_{0}}}{dx}\int_{0}^{\delta }{dy}=\frac{d{{V}_{0}}}{dx}\int_{0}^{\delta }{{{V}_{0}}dy}\)            (a)

式(10.78)等号左端第二项,按两数乘积的求导法则,有

\(\begin{align}& {{V}_{0}}\frac{d}{dx}\left( \int_{0}^{\delta }{{{V}_{x}}dy} \right)=\frac{d}{dx}\left( {{V}_{0}}\int_{0}^{\delta }{{{V}_{x}}dy} \right)-\frac{d{{V}_{0}}}{dx}\left( \int_{0}^{\delta }{{{V}_{x}}dy} \right)= \\& \frac{d}{dx}\left( \int_{0}^{\delta }{{{V}_{0}}{{V}_{x}}dy} \right)-\frac{d{{V}_{0}}}{dx}\left( \int_{0}^{\delta }{{{V}_{x}}dy} \right) \\\end{align}\)            (b)

将式(a)、式(b)带入式(10.78),可得

\(\frac{d{{V}_{0}}}{dx}\int_{0}^{\delta }{\left( {{V}_{0}}-{{V}_{x}} \right)dy}+\frac{d}{dx}\int_{0}^{\delta }{{{V}_{x}}\left( {{V}_{0}}-{{V}_{x}} \right)dy}=\frac{{{\tau }_{w}}}{\rho }\)

根据δ*和δ**的定义式,上式可进一步化成

\(\frac{d{{V}_{0}}}{dx}{{V}_{0}}{{\delta }^{*}}+\frac{d}{dx}\left( V_{0}^{2}{{\delta }^{**}} \right)=\frac{{{\tau }_{w}}}{\rho }\)

展开合并同类项,最后得到

\(\frac{d{{\delta }^{**}}}{dx}+\frac{1}{{{V}_{0}}}\frac{d{{V}_{0}}}{dx}\left( 2{{\delta }^{**}}+{{\delta }^{*}} \right)=\frac{{{\tau }_{w}}}{\rho V_{0}^{2}}\)               (10.79)

式(10.79)即为附面层动量积分方程。

式(10.79)中,一共有4个未知数,即V0,τw,δ*和δ**,其中,V0是由理想流体计算获得的,而δ**和δ*则由Vx和δ决定,因此方程尚有三个未知量τw,δ和Vx。在求解式(10.79)时,通常补充附面层内速度分布Vx=f(x)和壁面摩擦切应力τw的表达式。

10.7.2  速度分布应满足的边界条件

用积分法求解附面层时,需要补充附面层内的速度分布。虽然所选定的速度分布不能精确地表示附面层内的流动,但是可以精确地满足边界上的条件。

在附面层外边界上,黏性流体可以近似地看做理想流体。因此在外边界上,它们的速度等于外边界速度,速度的各阶导数都等于零,即

当y=δ时,

\({{V}_{x}}={{V}_{0}},\frac{{{\partial }^{n}}{{V}_{x}}}{\partial {{y}^{n}}}=0\)        (n=1,2,3…)             (10.80)

当y=0时,应满足无滑移条件

Vx=0   ,   Vy=0               (10.81)

如果将此条件用于附面层动量微分方程,即

\({{V}_{x}}\frac{\partial {{V}_{x}}}{\partial x}+{{V}_{y}}\frac{\partial {{V}_{x}}}{\partial y}=-\frac{1}{\rho }\frac{dp}{dx}+\upsilon \frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{y}^{2}}}\)

则可得到另一个边界条件,即

当y=0时,

\(\frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{y}^{2}}}=\frac{1}{\mu }\frac{dp}{dx}=-\frac{{{V}_{0}}}{\upsilon }\frac{d{{V}_{0}}}{dx}\)               (10.82)

再把动量方程对y求导,有

\(\frac{\partial {{V}_{x}}}{\partial y}\left( \frac{\partial {{V}_{x}}}{\partial x}+\frac{\partial {{V}_{y}}}{\partial y} \right)+{{V}_{x}}\frac{{{\partial }^{2}}{{V}_{x}}}{\partial y\partial x}+{{V}_{y}}\frac{{{\partial }^{2}}{{V}_{x}}}{\partial {{y}^{2}}}=\upsilon \frac{{{\partial }^{3}}{{V}_{x}}}{\partial {{y}^{3}}}\)

根据连续方程和无滑移条件,又可得到一个边界条件,即

当y=0时,

\(\frac{{{\partial }^{3}}{{V}_{x}}}{\partial {{y}^{3}}}=0\)               (10.83)

只要选定的速度分别满足边界条件,则表明它在物体表面和边界层外部附近都和真实速度分布接近。在附面层中间部分虽然可能有一定的误差,但是在应用积分法时,由于总体上满足动量积分方程,因此可以得到满足工程需要的结果。

在上述边界条件中,无滑移条件式(10.81)和压强梯度条件式(10.82)反映了物面及物面形状对速度分布的影响。因此在附面层计算中,为了保证一定计算精度,应满足这些条件。

10.7.3  平板不可压层流附面层计算

有一直匀流,速度为V∞,密度为ρ,流过如图10.9所示的平板。假设整个平板上全部为层流流动且平板的厚度无限薄,平板长度为l,宽度为b,下面用10.7.2小节介绍的附面层积分法对其进行求解。求解的内容有:速度近似分布;附面层厚度;切应力;摩擦阻力因数等。

根据假设,可以认为平板不影响附面层外的流动,仍然可以将附面层以外的流动看成是与平板平行的理想流动。于是,附面层外的流速V0=V∞,且沿平板V0=常数。将其代入动量积分关系式(10.79),则方程简化为

\(\frac{d{{\delta }^{**}}}{dx}=\frac{{{\tau }_{w}}}{\rho V_{\infty }^{2}}\)               (10.84)

为了求解式(10.84),需要补充两个关系式,即附面层内的速度分布和壁面上的摩擦应力关系式。

求速度分布的步骤:

首先假设速度分布为y的幂函数,即

\({{V}_{x}}={{a}_{0}}+{{a}_{1}}y+{{a}_{2}}{{y}^{2}}+{{a}_{3}}{{y}^{3}}+…+{{a}_{n}}{{y}^{n}}\)

式中,待定系数a0,a1,a3,…是未知的,它们必须由速度分布应遵循的边界条件确定;幂指数n可根据具体要求选取。实验证明,取n=2,即可与实验得到的速度分布曲线吻合很好,即

\({{V}_{x}}={{a}_{0}}+{{a}_{1}}y+{{a}_{2}}{{y}^{2}}\)

然后,确定式中的三个系数,它们必须由三个边界条件确定。这些边界条件如下:

(1)在物面上y=0,Vx=0代入上式,得a0=0。

(2)在附面层外边界上y=δ,Vx=V∞,可得V∞=a1δ+a2δ²。

(3)在附面层外边界上y=δ,∂Vx/∂y= 0可得a1+2a2δ=0。

由以上各式,可以确定

\({{a}_{0}}=0,{{a}_{1}}=\frac{2{{V}_{\infty }}}{\delta },{{a}_{2}}=-\frac{{{V}_{\infty }}}{{{\delta }^{2}}}\)

于是,速度分布为

\({{V}_{x}}=\frac{2{{V}_{\infty }}}{\delta }y-\frac{{{V}_{\infty }}}{{{\delta }^{2}}}{{y}^{2}}\)

或                                          \(\frac{{{V}_{x}}}{{{V}_{\infty }}}=\frac{2y}{\delta }-\frac{{{y}^{2}}}{{{\delta }^{2}}}\)               (10.85)

需要补充的第二个关系式是牛顿内摩擦定律,它提供了τw的关系式,即

\({{\tau }_{w}}=\mu {{\left( \frac{\partial {{V}_{x}}}{\partial y} \right)}_{y=0}}=2\mu \frac{{{V}_{\infty }}}{\delta }\)               (10.86)

利用补充方程式(10.85)、式(10.86)和动量积分方程式(10.84),联立求解即可得到附面层内所需要的有关结果。

由速度分布可求得动量损失厚度

\({{\delta }^{**}}=\int_{0}^{\delta }{\frac{{{V}_{x}}}{{{V}_{\infty }}}\left( 1-\frac{{{V}_{x}}}{{{V}_{\infty }}} \right)}dy=\int_{0}^{\delta }{\left( \frac{2y}{\delta }-\frac{{{y}^{2}}}{{{\delta }^{2}}} \right)\left( 1-\frac{2y}{\delta }+\frac{{{y}^{2}}}{{{\delta }^{2}}} \right)}dy=\frac{2}{15}\delta \)

于是,有                                \(\frac{d{{\delta }^{**}}}{dx}=\frac{2}{15}\frac{d\delta }{dx}\)               (10.87)

将式(10.86)、式(10.87)代入式(10.84),得

\(\frac{2}{15}\frac{d\delta }{dx}=\frac{2\mu {{V}_{\infty }}}{\delta \rho V_{\infty }^{2}}=\frac{2\mu }{\delta \rho {{V}_{\infty }}}\)

整理上式后,得

\(\int_{0}^{\delta }{\delta }d\delta =\int_{0}^{x}{\frac{15\mu }{\rho {{V}_{\infty }}}}dx\)

积分结果为                                       \(\frac{{{\delta }^{2}}}{2}=\frac{15\mu }{\rho {{V}_{\infty }}}x\)

故得附面层厚度随x的变化关系为

\(\delta =5.477\sqrt{\frac{\upsilon x}{{{V}_{\infty }}}}\)               (10.88a)

或                                          \(\frac{\delta }{x}=\frac{5.477}{\sqrt{{{\operatorname{Re}}_{x}}}}\)               (10.88b)

式中,\({{\operatorname{Re}}_{x}}={{V}_{\infty }}x/\upsilon \),是距平板前缘为x处的当地雷诺数。由式(10.88a)和式(10.88b)可见,层流附面层厚度与x的平方根成正比,与当地雷诺数的平方根成反比。

将δ(x)代回式(10.86),经化简后可得平板表面上的切应力分布为

\({{\tau }_{w}}=\mu {{\left( \frac{\partial {{V}_{x}}}{\partial y} \right)}_{y=0}}=2\mu \frac{{{V}_{\infty }}}{\delta }=0.365\sqrt{\frac{\rho \mu V_{\infty }^{3}}{x}}\)               (10.89)

当地摩擦阻力因数Cf定义为\({{C}_{f}}=\frac{{{\tau }_{w}}}{\frac{1}{2}\rho V_{\infty }^{2}}\),将式(10.89)代入,即可得到层流附面层的Cf表达式,即

\({{C}_{f}}=\frac{0.73}{\sqrt{{{\operatorname{Re}}_{x}}}}\)               (10.90)

作用在宽度为b的平板的上表面摩擦附力为

\({{X}_{f}}=\int_{0}^{t}{{{\tau }_{w}}bdx=0.73}V_{\infty }^{\frac{3}{2}}b\sqrt{\rho l\mu }\)               (10.91)

整个平板的上表面摩擦阻力因数定义为

\({{C}_{D}}=\frac{{{X}_{f}}}{\frac{1}{2}\rho V_{\infty }^{2}bl}=\frac{1.46}{\sqrt{{{\operatorname{Re}}_{t}}}}\)               (10.92)

式中                                        \({{\operatorname{Re}}_{t}}=V_{\infty }^{{}}l/\upsilon \)

10.7.4  光滑平板不可压湍流附面层计算

一般情况,如果绕物体的附面层不发生严重的脱体现象,曲壁附面层的摩擦阻力与平板情形相差不大,因此可以简化计算。

1.光滑平板湍流附面层

当流动雷诺数足够大时,在靠近平板前缘一段是层流附面层,而靠近平板后一段是湍流附面层。下面讨论假设平板从前缘开始就是湍流附面层的情况。

为了求解湍流附面层,根据普朗特的假设:沿平板的附面层流动与管流的情况没有显著的差别。因此对于充分发展的平板湍流附面层,可以把它看做管流。其中,附面层厚度已达到管道半径,管中心的最大速度Vmax相当于附面层外边界的速度V∞。实验证明,当Re=V∞x/v < 106时,平板湍流附面层的速度分布与管流的速度分布一致。切应力的关系也可采用圆管的结果。

湍流流动的速度分布可以根据半经验的对数分布规律,也可以根据经验的幂次方的分布规律。现采用后者作为第一个补充方程,即

\(\frac{{{V}_{x}}}{{{V}_{\infty }}}={{\left( \frac{y}{\delta } \right)}^{\frac{1}{7}}}\)               (10.93)

代入附面层动量损失厚度的表达式,可得

\({{\delta }^{**}}=\int_{0}^{\delta }{\frac{{{V}_{x}}}{{{V}_{\infty }}}\left( 1-\frac{{{V}_{x}}}{{{V}_{\infty }}} \right)}dy=\int_{0}^{\delta }{{{\left( \frac{y}{\delta } \right)}^{\frac{1}{7}}}\left[ 1-{{\left( \frac{y}{\delta } \right)}^{\frac{1}{7}}} \right]}dy=\frac{7}{72}\delta \)

故                                          \(\frac{d{{\delta }^{**}}}{dx}=\frac{7}{72}\frac{d\delta }{dx}\)               (10.94)

第二个补充方程为τw的关系式。对于光滑圆管中的湍流流动,当Re≤105时,沿程损失因数

\(\lambda =\frac{0.3164}{{{\operatorname{Re}}^{\frac{1}{4}}}}\)               (10.95)

式中,\(\operatorname{Re}=\rho {{V}_{av}}d/\mu ,{{V}_{av}}\)是平均速度,当用1/7次方速度分布时,它与圆管轴线上的速度Vmax的关系是Vav=0.817Vmax。当将圆管中的结果用于附面层计算时,要用附面层厚度去代替管径,即δ=d/2:用附面层外边界上的速度V∞去代替Vmax。这样,应用壁面切应力τw与△p的关系,并应用式(10.95)就可得到τw的表达式,即

\(\begin{align}& {{\tau }_{w}}=\frac{\Delta p{{r}_{0}}}{2l}=\frac{\lambda }{8}\rho V_{av}^{2}=\frac{1}{8}\times 0.3164{{\operatorname{Re}}^{-0.25}}\rho {{\left( 0.817{{V}_{\infty }} \right)}^{2}}= \\& \frac{1}{8}\times 0.3164{{\left( \frac{\rho \times 0.817{{V}_{\infty }}\times 2\delta }{\mu } \right)}^{-0.25}}\rho {{\left( 0.817{{V}_{\infty }} \right)}^{2}}= \\& 0.0233\rho V_{\infty }^{2}{{\left( \frac{\mu }{\rho {{V}_{\infty }}\delta } \right)}^{0.25}} \\\end{align}\)               (10.96)

将式(10.94)、式(10.96)代入附面层积分关系式(10.84),得

\(0.0233\rho V_{\infty }^{2}{{\left( \frac{\mu }{\rho {{V}_{\infty }}\delta } \right)}^{0.25}}=\frac{7}{72}\rho V_{\infty }^{2}\frac{d\delta }{dx}\)

简化后得到

\(0.24{{\left( \frac{\mu }{\rho {{V}_{\infty }}} \right)}^{0.25}}\int_{0}^{x}{dx}=\int_{0}^{\delta }{{{\delta }^{0.25}}d\delta }\)

积分后得到附面层厚度随x的变化关系为

\(\delta =0.382{{\left( \frac{\mu }{\rho {{V}_{\infty }}x} \right)}^{\frac{1}{5}}}x\)               (10.97)

或                                                       \(\frac{\delta }{x}=\frac{0.382}{\operatorname{Re}_{x}^{\frac{1}{5}}}\)

应用式(10.96)、式(10.97),可以得到平板湍流附面层当地摩擦因数为

\({{C}_{f}}=\frac{{{\tau }_{w}}}{\frac{1}{2}\rho V_{\infty }^{2}}=\frac{0.0592}{\operatorname{Re}_{x}^{\frac{1}{5}}}\)               (10.98)

平板上部的摩擦阻力因数及摩擦阻力分别为

\({{C}_{D}}=\frac{{{X}_{f}}}{\frac{1}{2}\rho V_{\infty }^{2}bl}=\frac{0.074}{\operatorname{Re}_{t}^{\frac{1}{5}}}\)               (10.99)

\({{X}_{f}}=\int_{0}^{t}{{{\tau }_{w}}bdx=}\frac{0.074}{{{\left( \frac{{{V}_{\infty }}l}{\upsilon } \right)}^{1/5}}}\frac{1}{2}\rho V_{\infty }^{2}bl\)               (10.100)

上面的公式是应用1/7次方速度分布得出的结果,一般认为在5 × 105 < Ret < 107范围内较合适,随着Ret的增加,偏差也增大。

2.湍流附面层与层流附面层的比较

湍流附面层与层流附面层在基本特性上有较大差别。

(1)湍流附面层的速度分布曲线比层流速度分布曲线要饱满得多,湍流附面层内流体平均动量比层流的大,因此湍流附面层不易分离。

(2)湍流附面层的厚度比层流附面层的厚度增长得快,因为湍流附面层的厚度δ与\({{x}^{0.8}}\)(见式(10.97))成正比,而层流附面层的δ与\({{x}^{0.5}}\)成正比,可见湍流附面层比层流附面层要厚得多。

(3)对于湍流附面层来说,作用在平板上的摩擦附力\({{X}_{f}}\)与\(V_{\infty }^{1.8}\)及\({{l}^{0.8}}\)成正比,对于层流附面层来说,作用在平板上的摩擦附力\({{X}_{f}}\)与\(V_{\infty }^{1.5}\)及\({{l}^{0.5}}\)成正比。因此,从减小摩擦阻力来看,层流附面层将优于湍流附面层。

10.7.5  光滑平板混合附面层计算

在高雷诺数的情况下,绕物体流动的附面层往往是混合附面层,即从平板前缘开始先是一段层流附面层,经过过渡段再变为湍流附面层,如图10.11所示。在计算中忽略过渡段,即认为从转捩点开始,都是湍流附面层。混合附面层的摩擦阻力计算方法如下:

令    l——平板总长度;

       \({{x}_{T}}\)——平板上层流附面层长度;

       \({{C}_{D}}\)——从前缘开始平板上全为湍流附面层时的摩擦阻力因数;

       \({{C}_{D\text{t}}}\)——\({{x}_{T}}\)段上湍流附面层的摩擦阻力因数;

       \({{C}_{Dt}}\)——\({{x}_{T}}\)段上层流附面层的摩擦阻力因数;

       \({{C}_{D\text{m}}}\)——混合附面层的摩擦阻力因数。

根据摩擦阻力因数的定义,可得混合附面层的摩擦阻力为

\(\frac{1}{2}\rho V_{\infty }^{2}l{{C}_{D\text{m}}}=\frac{1}{2}\rho V_{\infty }^{2}\left( l{{C}_{D}}-{{x}_{T}}{{C}_{D\text{t}}}+{{x}_{T}}{{C}_{Dt}} \right)\)

故                                                     \({{C}_{D\text{m}}}={{C}_{D}}-\frac{{{x}_{T}}}{l}\left( {{C}_{D\text{t}}}-{{C}_{Dt}} \right)\)

又                                                      \({{\operatorname{Re}}_{t}}=\frac{{{V}_{\infty }}l}{\upsilon },{{\operatorname{Re}}_{cr}}=\frac{{{V}_{\infty }}{{x}_{T}}}{\upsilon }\)

因此                                      \({{C}_{D\text{m}}}={{C}_{D}}-\frac{{{\operatorname{Re}}_{cr}}}{{{\operatorname{Re}}_{t}}}\left( {{C}_{D\text{t}}}-{{C}_{Dt}} \right)={{C}_{D}}-\frac{C}{{{\operatorname{Re}}_{t}}}\)               (10.101)

式中                                                \(C=0.074\operatorname{Re}_{cr}^{4/5}-1.328\operatorname{Re}_{cr}^{1/2}\)

当Ret很大时,                                      \({{C}_{D\text{m}}}\approx {{C}_{D}}\)

以上几节讨论了压强梯度为零(平板)的情况,对于有压强梯度时的曲壁附面层计算可参看黏性流体力学的教材及有关文献。

10.8  附面层分离与控制

10.8.1  附面层分离

在流体流动中,一方面,对于平板附面层,在平板上的压强为常数,即∂p/∂x=0;而对于曲壁附面层,沿物体表面可能存在压强梯度∂p/∂x > 0,即流体是在逆压梯度下流动的,因而速度迅速衰减。另一方面,流体沿壁面流动时,附面层厚度逐渐增加,由于黏性摩擦影响,靠近壁面处动能有很大损失。在这双重的作用下,使得靠近壁面某点S处的流体停止流动,结果使S点之后的附面层内部的流体出现倒流的现象。S点称为分离点。由图10. 12可见,在M点之前,∂p/∂x<0,∂Vx/∂y>0,在M点之后,∂p/∂x>0,∂Vx/∂y<0,在分离点S处,\({{\left[ \frac{\partial {{V}_{x}}}{\partial y} \right]}_{y=0}}=0\)。由以上分析可知,只有在逆压梯度下,即∂p/∂x>0,才有可能出现分离。因此, ∂p/∂x>0是分离的必要条件。在逆压梯度下,满足\({{\left[ \frac{\partial {{V}_{x}}}{\partial y} \right]}_{y=0}}=0\)的条件,是判别分离的准则。总之,附面层分离只可能在逆压梯度的条件下发生。

附面层分离点S的准确计算非常困难,因为点S是在附面层厚度很小之处,并按外部势流场的压力分布求出的。在分离点之前的流动阻力可按前述的附面层方法计算得到。而一旦附面层出现分离,附面层厚度增加很快,即δ<<L(物体特征尺寸)的条件不再满足,此时附面层理论失效,因此要用完整的黏性流动方程来求解。

附面层分离后,流动中出现了旋涡,结果是流动附力急剧增加。这是由于物体表面上产生流动方向的压力差所致,即所谓的压差阻力。流动阻力包括了压差阻力和摩擦阻力。摩擦阻力主要取决于附面层的流动状态(层流或湍流),压差阻力则主要与附面层的分离有关。

10.8.2  附面层控制

附面层分离会使流体的一部分机械能损失,流体绕物体的阻力急剧增加,发动机各部件效率降低;有时甚至产生不稳定流动,以致造成发动机的损坏。因此在设计时,应尽量避免大范围内的附面层分离。预防和推迟附面层分离是工程设计中应关注的问题。

附面层分离是流体质点在运动中,由于黏性摩擦和逆向压差的共同作用所造成的。有许多方法可以控制和防止附面层分离,这里介绍几种有效的方法。

1.高速气流喷入附面层

附面层分离的原因之一是黏性使得附面层内的流体速度降低,动能减小。因此防止附面层分离的方法之一是向附面层内注入高速气流,使得附面层内的流体质点重新获得能量。方法是在物体内部设置气源,将高速射流从附面层将要分离之处喷入,使其避免分离。

2.附面层的吸入

在附面层容易分离的物面上设置狭缝,通过吸气装置把靠近物体表面的低能气流吸入物体内,使附面层厚度变薄,靠近物体表面处的气流具有较大的流速,可以有效地消除附面层分离。

3.安装涡流发生器

湍流附面层比层流附面层的速度分布较饱满,在物体壁面附近流体质点的动能较大,能够承受较大的逆压梯度,因此湍流附面层比层流附面层不易分离。由此可见,如果在有逆压梯度的通道里或物面上安装一些涡流发生器(能产生诱导涡的小翼形叶片),使附面层提前变成湍流,可有效防止或推迟分离。

4.设计措施

在设计亚声速通道时,为了减小过大的逆压梯度,应避免扩张通道的扩张角过大。扩张角过大,附面层容易分离,扩张角过小,达到同样面积比的管道长度大,因而摩擦阻力大。一般亚声速的扩压器的扩张角在6°~12°的范围内。等压强梯度的扩压器比锥形扩压器分离较迟。对于较短的管道,采用等压强梯度的扩压器很有利。

小结

本章讨论黏性流体流动基础及有关基本概念。主要讨论了N一S方程,湍流流动的雷诺方程和附面层的基本方程。

(1)雷诺方程与雷诺应力。

(2)附面层微分方程、积分方程及有关附面层的特点、分离和控制。

(3)速度附面层是靠近物体表面速度梯度很大的簿层。在附面层内,若物面曲率半径较大,则压强沿附面层法线方向近似不变。

思考与练习题

10.1  湍流流动用什么方法研究?如何描述湍流脉动现象?

10.2  附面层分离的物理原因是什么?如何判别流动分离?

10.3  附面层内的流动有何特点?

10.4  某一黏性流体绕物体作定常的三维流动,已知速度分布为Vx=3y+ 2z,Vy=x+2z,Vz=3x+2y,单位为m/ s。设流体的黏度系数μ=0.01Pa • s,求流场中的切应力τxy,τyz和τzx。

10.5   己知不可压流动速度场为Vx=2x ,Vy=-2y。求流场中的正应力τxx,τyy和τzz。

10.6   定常二维不可压黏性流动,黏度为μ。若己知流函数ψ=-xy,忽略质量力,试通过积分N一S方程,找出压强与x,y之间的关系。

10.7   水在两块无限大平行平板之间流动,下板静止,上板以V∞=0.5 m/s的恒定速度沿x方向运动,设两板间距h=2.5 mm,Vy=Vz=0,且黏度为μ=0.01Pa • s。试求:

(1)通过某截面体积流量为零时的压强梯度;

(2)设沿x方向无压强变化,即p=const,试求两板之间的速度分布。

10.8   在10.7题中,假设上板以V∞=0.3 m/s,压强梯度dp/dx=-200Pa/m,试求水流的最大速度及其所在的位置。

10.9   两块无限大平行平板之间,充满不可压缩绝热黏性流体,己知两板间距h=025m,压强p=const,黏度为μ=1 × 106Pa • s,上板以速度V∞=0.15 m/ s运动,下板静止不动。求流体单位体积的内能增加率。

10.10   已知平板层流附面层的速度分布分别为

\[\begin{align}& \left( 1 \right)\frac{{{V}_{x}}}{{{V}_{\infty }}}=\frac{y}{\delta } \\& \left( 2 \right)\frac{{{V}_{x}}}{{{V}_{\infty }}}=\frac{2y}{\delta }-\frac{{{y}^{2}}}{{{\delta }^{2}}} \\& \left( 3 \right)\frac{{{V}_{x}}}{{{V}_{\infty }}}=\sin \left[ \frac{\pi y}{2\delta } \right] \\\end{align}\]

试求附面层特性δ/δ*和δ**。

10.11   空气以速度V∞=75m/ s,压强p=1 × 105Pa,温度t=20 ℃流过水平放置的光滑平板,平板长2 m,宽1 m,空气的运动黏度为υ=1.5 × 10-5m2/s。

(1)若平板分别为层流和湍流附面层,求平板末端的附面层厚度和总阻力,并比较之。

(2)设临界雷诺数Recr=3 × 105,求平板总阻力。

11.12   在海平面高度,一列长120 m,宽和高均为3 m的火车,以V=100 m/s的速度行驶,顶面和两边侧面可看做光滑平板。求这三个面上所受的总摩擦阻力。

发表回复

您的电子邮箱地址不会被公开。 必填项已用 * 标注

You cannot copy content of this page