传热学

一维稳态热传导的傅里叶定律

q=-kA\frac{dT}{dx}

其中q是热流密度,表示的是单位时间通过这个截面的热量,k是对流交换系数,表征了传热速度,一般来说对于金属而言,其传热速率较快,木头则较慢。A是截面的面积,面积越大,单位时间内传输的热量便越大。

傅里叶定律也可以表示成如下的形式:

q=-kA\nabla{T}

因为一般认为\nabla{T}等于T对x的导数。

特殊地,对于一个均匀的细杆来说,其温度是随着位移均匀变化的,于是便有:

q=-kA\frac{T_2-T_1}{L}

一维稳态热传导的傅里叶定律中,热流密度对于温度梯度之前存在一个负号,这是因为温度梯度是\frac{dT}{dx},但是由于自然界中,热量总是由高温流向低温的,如下图的说明所示,因此所有热流密度都要在温度梯度之前加上一个符号。

一维非稳态热传导方程

一维非稳态热传导方程和稳态热传导方程的区别在于温度是否是随时间均匀变化的,如果是那么就是稳态,否则都是非稳态状态。

由于在这个条件下,T不光只是x的函数,还是时间t的函数,因此不能简单使用\frac{dT}{dx},而需要从基本的能量守恒公式开始推起:

由于热能的变化率等于{\rho}c_pAdx\frac{\partial{T}}{\partial{t}}=q_x-q_{x+dx}其中,ρ是材料的密度,c_p是比热容,A是接触面积,q是热流密度,又由于q_x = -k A \left. \frac{\partial T}{\partial x} \right|_{x} q_{x+dx} = -k A \left. \frac{\partial T}{\partial x} \right|_{x+dx} ,k是上文所提到的热对流交换系数,因此上式化简为:

{\rho}c_pAdx\frac{\partial{T}}{\partial{t}}=kA(\left. \frac{\partial T}{\partial x} \right|_{x+dx}-\left. \frac{\partial T}{\partial x} \right|_{x})

又根据牛顿莱布尼兹公式有\left. \frac{\partial T}{\partial x} \right|_{x+dx}-\left. \frac{\partial T}{\partial x} \right|_{x}=\frac{\partial^2T}{\partial{x^2}},因此得到一维非稳态热传导方程:

\frac{\partial{T}}{\partial{t}}=α\frac{\partial^2{T}}{\partial{x^2}}

其中α被称为是扩散率,其等于\frac{k}{{\rho}c_p},由上式可知,对于一般的稳态条件,则有\frac{\partial{T}}{\partial{t}}=0,因此后者也是0,则有T是x的一次函数关系,T=C_1x+C_2。同样的,对于这种情况,使用一位稳态的傅里叶定律有q=-kA\frac{dT}{dx},进行化简得到T=\frac{q}{-kA}x+C.标明一维非稳态热传导方程实际上是一维稳态热传导方程的更一般形式。

一维非稳态热传导方程的求解条件

由于一维非稳态热传导方程是一个二阶空间一阶时间的PDE,不像常微分方程那样有着简单的通解公式,需要有着特殊的求解条件才能够进行求解,常见的初始条件、边界条件和特殊条件如下所述:

对于二阶空间一阶时间的PDE至少需要一个关于时间的方程和2个关于空间的方程才能求解,一般选取的初始条件T(x,0)=f(x),这个方程标明了这个系统下的初始温度关于x的分布,注意,必须是知道了t=0这个时刻各个x所对应的T的分布才能够求解,否则仍是缺条件的,t=0也并不是限定死的,对于t=t_0也是可以的。

边界条件:边界条件分为以下三种:

  1. 定温型:有T=T0,即边界的温度与外界相通,不存在热传导。

  2. 定热流:-k\frac{\partial{T}}{\partial{n}}=q_0,即外界的热流密度是恒定的,这会出现在外界的热量是恒定的时候,虽然温度是随着位移x变化的,但是热流是固定的,这一般会出现在外界的热流是已知的一个常量或者一个已知的函数的情况下。这与第三类边界条件的区别在于第三类边界条件是与外界的自然传热,是由牛顿冷却定理驱动的,这常出现在恒功率的场景下。

  3. 对流边界:-k\frac{\partial{T}}{\partial{n}}=h(T-T_{\infty}),其中h是对流换热系数,该公式的本质即为牛顿冷却定律,其中T_{\infty}指得是环境的系统温度。

几类特殊边界条件:

  1. 绝热边界条件:-k\frac{\partial{T}}{\partial{n}}=0,其中\frac{\partial{T}}{\partial{n}}指的是温度沿边界法向的导数,注意这里是对温度法向(即热流方向的偏导数,不是x方向上的)对于这种绝热边界条件,一般出现在几种情况(此边界不允许热量通过的时候):
    隔热层(无热量的传递)
    对称面:此处由于处在整个温度梯度的峰谷,因此\frac{\partial{T}}{\partial{n}}必然等于0,这种情况一般会发生在圆柱的轴线、平板的中心处等,如下图所示:

  2. 辐射边界条件:这种情况一般发生在表面温度较高,与外部空气或者真空之间的热交换非常显著的时候,一般认为\Delta{T} 越大,热辐射就会越显著,其满足斯特藩-玻尔兹曼方程,其计算公式为

    -k\frac{\partial{T}}{\partial{n}} = \varepsilon \sigma (T^4 - T_{\infty}^4)

其中,T_{\infty}表达的是环境温度,T是物体表面的温度,\varepsilon是材料的发射率,\sigma是斯特藩-玻尔兹曼常数,其值是5.67*10(-8)。

Biot数与集中参数方程

毕奥数的数学定义式为:

Bi= \frac{h L_c}{k}

其中h是对流换热系数,是牛顿冷却定律的参数,Lc​是特征长度,代表是一个代表物体尺寸的几何参数,通常定义为物体的体积(V)与其表面积(A)之比:L_c= \frac{V}{A}

  • 对于一个大平板,其特征长度是其厚度的一半。

  • 对于一个球体,其特征长度是其半径的三分之一。

  • 对于一个长圆柱体,其特征长度是其半径的二分之一。

k是物体的热导率。k越大,传热的速度越快,其是一维(非)稳态传热方程的重要参数。

毕奥数的一个重要用处在于其可以用来简化传热学模型:

  1. 当毕奥数小于0.1的时候说明这个系统的导热速度很快,因此可以认为温度能在很快的条件下达到一致,在这种情况下可以简化认为温度只是时间的函数,而与位移x是无关的,因此可以采用集中参数的方法简化PDE,将其转化为一个易于求解的常微分方程,且当毕奥数小于0.1的时候,这种PDE换ODE的方法可以带来极高的精度(误差一般不高于5%):
    集中参数方程:-mc_p \frac{dT}{dt} = hA(T-T_{\infty}),解这个常微分方程可以得到如下的表达式:(式中T_0表示的是0时刻的温度)

    \frac{T(t)-T_{\infty}}{T_0-T_{\infty}} = e^{-\left(\frac{hA}{\rho V c_p}\right)t}
  2. 当毕奥数大于1的时候存在极大的温度梯度,此时T与位移x之间不再是一致的,需要求解原偏微分方程。

热传导方程实例分析

长度为l的均匀杆,两端有恒定热流流入,热流的强度为q_0,写出这个热传导问题的方程的边界条件。

由于热流密度总是从高温流向低温的,其表达式如下:

q_0 = -k \frac{\partial u}{\partial n}

当x=l的时候,根据第二类边界条件,其表达式如下:

-k \left. \frac{\partial u}{\partial n} \right|_{x=l}=-k \left. \frac{\partial u}{\partial x} \right|_{x=l} = q_n = -q_0 \\

当x=0的时候,其边界条件的表达如下:

-k \left. \frac{\partial u}{\partial n} \right|_{x=0} = -k \left. \frac{\partial u}{\partial (-x)} \right|_{x=0} = q_n = -q_0

虽然最后的计算结果都是-q_0,但是其来历是完全不同的,热流密度是有高温指向低温环境的,而在这道题的情景中,棒的温度式中低于外界温度,因此两端的热流密度的方向都是由外界指向棒内,但是正是由于这个原因,左端和右端的热流密度的方向不同,导致外法线方向向量\vec{n}两者并不相同,更详细的说,是他们的方向是完全相反的。

棒的左边的\vec{n} 方向与x轴正方向相反,而右边的\vec{n} 方向则与x轴正方向是相同的,因此最后计算出来的q_0 方向也是相反的。

特别的,当棒的左侧和右侧都是绝热的时候,就讲不存在热流密度这一说,他们的热流密度都会变成0,如下式所示:

\left. \frac{\partial u}{\partial x} \right|_{x=l} = \left. \frac{\partial u}{\partial x} \right|_{x=0} = 0

热传导方程的推广

热传导方程不仅仅适用于传热学领域,其在其他物理领域,像波动学中的应用同样十分广泛,数学上一般把这些偏微分方程根据与其相似的二次代数方程进行类比,如下所述:

二维平面上的圆锥曲线的一般方程如下所示:

Ax^2 + Bxy + Cy^2 + Dx + Ey + F = 0

决定这个方程代表哪种曲线(椭圆、抛物线、双曲线)的关键在于其判别式 Δ = B² - 4AC

  1. 如果 B² - 4AC < 0,方程代表一个 椭圆 。

  2. 如果 B² - 4AC = 0,方程代表一个 抛物线。

  3. 如果 B² - 4AC > 0,方程代表一个 双曲线。

对于一个通用的二维二阶线性偏微分方程,其标准形式如下所示:

A\frac{\partial^2 u}{\partial x^2} + B\frac{\partial^2 u}{\partial x \partial y} + C\frac{\partial^2 u}{\partial y^2} + D\frac{\partial u}{\partial x} + E\frac{\partial u}{\partial y} + Fu = G

数学家们发现,这个方程的性质主要由其最高阶(二阶)的导数项决定,也就是 A, B, C 三个系数。

于是,他们借用了圆锥曲线的判别式,对这个偏微分方程进行分类:

  1. 椭圆型: 如果 B² - 4AC < 0

  2. 抛物型 : 如果 B² - 4AC = 0

  3. 双曲型: 如果 B² - 4AC > 0

\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = 0

这里ux 求了两次导数,所以 A=1uy 求了两次导数,所以 C=1。没有 ∂²u/∂x∂y 项,所以 B=0

判别式: B² - 4AC = 0² - 4(1)(1) = -4

结果: -4 < 0,因此它是 椭圆型 的。

波动方程的自变量是时间和空间,我们以一维波动方程为例,并将变量看作 tx

\frac{\partial^2 u}{\partial t^2} - c^2 \frac{\partial^2 u}{\partial x^2} = 0

为了对应 A, B, C,我们可以把 t 看作 yx 还是 x。方程变为 u_yy - c²u_xx = 0

所以A = -c²C = 1,没有 u_xy (即 u_xt) 项,所以 B=0

判别式: B² - 4AC = 0² - 4(-c²)(1) = 4c²

结果: 因为 c (波速) 是正实数,所以 4c² > 0,因此它是 双曲型 的。

同样,对于一维热传导方程,自变量是 tx

\frac{\partial u}{\partial t} - \alpha \frac{\partial^2 u}{\partial x^2} = 0

很显然,其是抛物型的。

对于一些常见的物理领域内的偏微分方程归纳如下,他们都有着形式雷同的特征。相同行的偏微分方程基本上是相似的。

双曲型

声波、电磁波、杆的振动

抛物型

热传导,物质扩散时的浓度变化规律,土壤力学中的渗透方程

椭圆型

稳定的浓度分布,静电场的电位,流体的势

:设弦的两端固定于x=0 和x=l,弦的初始位移如下图,初速度为零,求弦满足的定解问题。

如果设整个弦的波动的长度是u(x,t),那么便可进行如下的操作:

图片上的整个定解问题(包括波动方程、边界条件和初始条件)可以表示为如下的方程:

\frac{\partial^2 u}{\partial t^2} = a^2 \frac{\partial^2 u}{\partial x^2}, \quad 0 < x < l, \ t > 0\\ u|_{x=0} = u|_{x=l} = 0 \\ u|_{t=0} = \left\{ \begin{array}{ll} x, & 0 < x \le \frac{l}{2} \\ l-x, & \frac{l}{2} \le x < l \end{array} \right. \\ \frac{\partial u}{\partial t}\bigg|_{t=0} = 0

其中上式的第一个方程为控制方程,限制了波动方程满足弦的波动定律,第二个表达式则规定了弦的两端的初始位置都是0,第三个表达式则表达了在初筛条件下的弦形状的边界方程,最后一个表达式则是说整个弦在一开始都是固定死的(即一开始的时候整个弦是由静止开始释放的)。

高阶稳态传热方程

下面将一维传热方程向更高维进行推广,\frac{\partial^2 T}{\partial x^2} = {\alpha} \frac{\partial T}{\partial t} ,注意这个公式是在无内热源的情况下的,如果存在着内热源,上面的公式需要修正为下式,其中q指的是内部源(例如:系统内部发生化学反应、核反应等需要发生热变化的的过程),而前者均来自外界的传热与热源:

\frac{\partial^2 T}{\partial x^2} + \frac{\dot{q}}{k} = \frac{1}{\alpha} \frac{\partial T}{\partial t}

对一维传热方程进行推广,便可得到二维,三维的稳态热传导方程:

\frac{\partial^2 T}{\partial x^2} + \frac{\partial^2 T}{\partial y^2} + \frac{\partial^2 T}{\partial z^2} +\dot{q} = 0

其中当q不等于0时,方程为泊松方程,q=0时被称为拉普拉斯方程。

常见几何体的热传导方程

对于内径为r1,外径r2,高度为L的空心圆柱体来说,其满足拉普拉斯方程,

这个方程在圆柱坐标系下的通解公式如下所示:

令\begin{align*} x &= r \cos(\theta) \\ y &= r \sin(\theta) \\ z &= z \end{align*}则r = \sqrt{x^2 + y^2}, \theta = \arctan(y/x)\\

通过链式法则进行坐标变换,可以得到拉普拉斯算子在通用圆柱坐标系下的表达式。这是一个标准结果,我们直接引用:

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

其径向r上的热传导方程(不考虑其在θ方向和z方向上的热传导)可以表述为如下形式:

\frac{1}{r} \frac{d}{dr} \left( r \frac{dT}{dr} \right) = 0\\ \frac{d^2T}{dr^2} + \frac{1}{r} \frac{dT}{dr} = 0

其中其满足几个边界条件T(r1)=T1,T(r2)=T2。

这是一个典型的二阶常微分方程,求解这个方程便可得到如下通解:

T(r) = C_1 \ln(r) + C_2.

进一步展开计算,得到径向r的温度分布:

\frac{T(r) - T_1}{T_2 - T_1} = \frac{\ln(r/r_1)}{\ln(r_2/r_1)}

对于球体而言与空心圆柱体类似,有如下方程:

\frac{1}{r^2} \frac{d}{dr} \left( r^2 \frac{dT}{dr} \right) = 0\\ \frac{d^2T}{dr^2} + \frac{2}{r} \frac{dT}{dr} = 0

求解得到这个方程的通解为:

T(r) = -\frac{C_1}{r} + C_2.

进一步得到球体径向r所对应的温度分布:

\frac{T(r) - T_1}{T_2 - T_1} = \frac{\frac{1}{r_1} - \frac{1}{r}}{\frac{1}{r_1} - \frac{1}{r_2}}

由上式可知,对于球体而言,其温度分布函数与球体的半径成反比关系。

流体力学

基本概念

  • 可压缩与不可压缩的判断:一个物体能否压缩取决于其密度是否会随压强的变化而变化。

一般来说,固体、液体都是不可压缩的,对于气体需要分类讨论,对于流速较慢的气体(流速小于0.3马赫,约102m/s)是不可压缩的,但气体的速度超过0.3马赫之后则认为其是可以压缩的。

  • 粘度:所说的粘度一般指的是动力粘度,即阻止物体运动的粘度,一般来说,粘度越大,流量就会越小,其受到的阻力就会越大。一般我们讨论的都是牛顿流体,即流体的粘度是恒定的,不会随着受力的情况的改变而改变,对于这种牛顿流体,其满足牛顿粘性方程:

    \tau=\mu\frac{du}{dy}

其中\tau 是水平剪切力,\mu 是动力粘度,\frac{du}{dy}则是速度梯度,如图所示,一个物体想要发生移动,就必须克服流体对其的水平剪切力:

当将这个流体所处的环境推广到细长的管道后便可得到如下的公式:

  • 哈根-泊肃叶定律:

    Q=\frac{{\pi}R^4\Delta{P}}{8{\mu}L}

其中\Delta{P} 指的是细管前后端的压力差,Q是体积流量,但是这个公式有着较为严苛的适用前提:对于一个在又细又长的圆管中做平稳、缓慢的层流运动的流体。平稳指的是流量不会随着时间变化而发生变化,缓慢则是用来限制湍急、混乱的湍流。

一般来说\mu越大,流体越容易做规整运动(层流),\mu 越小流体越容易形成湍流(混乱运动)

  • 雷诺数:用来判断一个流体的粘度是否可以忽略不计,其计算公式为:

    Re=\frac{{\rho}vL}{\mu}

其中\rho是流体的密度,v 是流体的平流流速,L是特征长度,其取值是取决于具体的流动场景的,对于不同的场景,其特征长度也不同。具体情况需要具体分析。

一般来说当Re较小的时候,说明流体的粘性占据主导地位,物体大多采取层流的手段。
当Re值较大时,说明此时粘性是可以忽略的,流体不再受到粘度的约束,多采取湍流的手段。

流体的静力学与运动学

从如下的基本模型开始推起,如下是一个底面积为S,高为z的几何体

对这个物体进行受力分析:

取高度为dz的微元,则有\\ V=Sdz,m=S{\rho}dz\\ {\rho}Sgdz+pS=(p+dp)S\\ 化简之后即可得到流体静力学基本方程:\\ \frac{dp}{dz}={\rho}g

对于这个常微分方程,他的解是极易求解的,不难得到其解为p=p_0+{\rho}gz。这表明对于不同的流体来说,密度越大,压力变化的也会越快!

流体力学-运动学

拉格朗日描述法:观察流体运动过程中单独一个质点的运动轨迹。设t=0时刻质点所在的轨迹是(a,b,c)设在t时刻该质点的坐标表示为:

\begin{cases} x = x(a, b, c, t) \\ y = y(a, b, c, t) \\ z = z(a, b, c, t) \end{cases}

由于速度是位移的导数,那么这个质点在t时刻的速度可以表达为下式:

\begin{cases} u = \frac{\partial{x}}{\partial{t}} \\ v = \frac{\partial{y}}{\partial{t}} \\ w = \frac{\partial{z}}{\partial{t}} \end{cases}

其中对于速度而言,我们一般采用u,v,w表示质点的速度,x,y,z则用于衡量其坐标形式。拉格朗日表述法的优点在于直观清晰,尤其是对于单个质点的时候,这种方法对单个质点的分析非常有效。但是其缺点也很明显,如果想对整个流体系统都采用拉格朗日表述的话,计算量将会是非常大的。

实际上在工程建模上最常用的不是拉格朗日表述,而是欧拉表述。

欧拉法与拉格朗日表述的不同在于其不再是利用各个点处的位移,而是采用流场的思想:

v=x \vec{i} + y \vec{j} +z\vec{k}
对于x=x(t),y=y(t),z=z(t)其位移可以表达为如下形式:\\ \phi(t)=\phi({x(t),y(t),z(t)}),对此式进行求导,得到:\\ \frac{d\phi}{dt}=\frac{\partial{\phi}}{\partial{x}}\frac{dx}{dt}+\frac{\partial{\phi}}{\partial{y}}\frac{dy}{dt}+\frac{\partial{\phi}}{\partial{z}}\frac{dz}{dt}+\frac{\partial{\phi}}{\partial{t}}

由于u,v,w被表示为流体力学里的速度,因此有:

\begin{cases} \frac{dx}{dt}=u\\ \frac{dy}{dt}=v\\ \frac{dz}{dt}=w \end{cases}

代入上式可以得到如下表达式:

\frac{d\phi}{dt} = \underbrace{ \left( u\frac{\partial\phi}{\partial x} + v\frac{\partial\phi}{\partial y} + w\frac{\partial\phi}{\partial z} \right) }_{\text{迁移加速度 (Convective)}} + \underbrace{ \frac{\partial\phi}{\partial t} }_{\text{局部加速度 (Local)}}

关于上式中的\frac{\partial{\phi}}{\partial{t}},其表示的是局部加速度,即自身的加速度,关于前者迁移加速度,这是由于质点自己的运动,导致其进入了一个新的区域,在这个区域里,速度发生了突然地改变,这种变换的根源不是因为时间,而是因为空间的骤变带来了巨大影响。

对于上式中的\frac{\partial{\phi}}{\partial{t}} ,其为局部加速度,当这个局部加速度的值为0时(\frac{\partial{\phi}}{\partial{t}} =0),我们称其进入了稳态流状态,在这个情况下,速度不会随着时间的变化而变化,其只与空间存在关系,在这个情况下,流线的形状一般是固定的。

例如:在一个二维的例子中,假设速度场可以用如下的表达式表示:

假设\vec{u}(x,y,t)=(2x+t)\vec{i}-2y\vec{j}为这个系统的速度场\\ 求t=3s时,点(1,2)处的加速度大小\\ 在这个二维问题中,u(x,y,t)=2x+t,v(x,y,t)=-2y\\ a_x=\frac{\partial{u}}{\partial{t}}+u\frac{\partial{u}}{\partial{x}}+v\frac{\partial{u}}{\partial{y}}\\ a_y=\frac{\partial{v}}{\partial{t}}+u\frac{\partial{v}}{\partial{x}}+v\frac{\partial{v}}{\partial{y}}\\

分别计算u,v对x,y,t的导数,便可得到:

a_x=1+4x+2t=11\\ a_y=-4y=-8

因此最后的加速度为11\vec{i}-8\vec{j}

通过python编程求解得到最后的速度热力图如下所示,越靠近中心的速度越慢,越往边上走速度越快:

连续性方程(质量守恒定律)

从如下的物理例子中推导出连续性方程:

如图所示,存在一个顶点为(x,y,z),边长分别是\Delta{x},\Delta{y},\Delta{z},则这个长方体的体积为\Delta{V}=\Delta{x}\Delta{y}\Delta{z},当其密度为\rho 时,其总质量为\Delta{M}={\rho}\Delta{x}\Delta{y}\Delta{z},易知其质量的增加速率为\frac{\partial{M}}{\partial{t}}=\frac{\partial{\rho}}{\partial{t}}\Delta{x}\Delta{y}\Delta{z}

根据质量守恒定律,质量增加的速率必然等于净流入速率,也即流入的质量速率-流出的质量速率。

由于三个方向上是完全等价对称的结构,因此可以只分析一侧的情况来推得其他面的具体情况,操作方式如下:

对于x方向的流动对于顶点x而言\\ \dot{m_{入}}={\rho}\Delta{y}\Delta{z}u_x\\ 在x+dx处有\dot{m_{出}}={\rho}\Delta{y}\Delta{z}u_{x+dx}

由于

由于\\ \rho u_x(x) - \rho u_x(x+dx) = - \frac{\partial (\rho u_x)}{\partial x} dx\\ \underbrace{ (\rho u_x)|_x }_{\text{流入通量}} - \underbrace{ \left( (\rho u_x)|_x + \frac{\partial (\rho u_x)}{\partial x}dx \right) }_{\text{流出通量 (泰勒展开)}} \approx - \frac{\partial (\rho u_x)}{\partial x}dx\\ m_{净}=-\frac{\partial (\rho u)}{\partial x}\Delta{x}\Delta{y}\Delta{z}

同理,对于y方向,z方向的分析是等价的,因此则流入的质量速率减去流出的质量速率等于-(\frac{\partial (\rho u)}{\partial x}+\frac{\partial (\rho v)}{\partial y}+\frac{\partial (\rho w)}{\partial z})\Delta{x}\Delta{y}\Delta{z},联立之前得到的流体质量增加速率\frac{\partial{M}}{\partial{t}}=\frac{\partial{\rho}}{\partial{t}}\Delta{x}\Delta{y}\Delta{z} 并化简得到连续性方程:

\frac{\partial \rho}{\partial t} + \frac{\partial(\rho u)}{\partial x} + \frac{\partial(\rho v)}{\partial y} + \frac{\partial(\rho w)}{\partial z} = 0

对于上式,如果该流体是不可压缩流体,由于其密度不变,那么上式可以简化为:

\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = 0

由连续性方程可以得到管道情况下的一个重要推论方程(即在任何情况下,横截面积乘上对应的速度都是一个定值,这就可以解释当将一个水管的出口堵上一截后,射程会大幅变远。):

A_1v_1=A_2v_2=常数

最后给出流量的计算公式:

Q=Av\\ Q = \iint_S \mathbf{u} \cdot d\mathbf{A}

上式仅可适用于横截面积和流体的平动速度都是定值的情况下,而下式则更具有普遍性,利用曲面积分计算更广泛的解。

伯努利方程:

能量形式:\boxed{\frac{V^2}{2} + gz + \frac{p}{\rho} = \text{常数}}\\ 压力形式:\boxed{\frac{1}{2}\rho V^2 + \rho gz + p = \text{常数}}

在应用此方程前,必须首先确认流动满足以下三个适用前提:

1. 稳态流动 : 流场中任意点的流体参数(速度、压力等)(需满足\frac{\partial{v}}{\partial{t}})不随时间改变。

2. 不可压缩流动: 流体密度\rho 视为常数。

3. 无粘流动 : 忽略流体内部的粘性摩擦。粘性会消耗机械能,(导致系统内部的能量不守恒)使其转化为内能,导致伯努利方程中的常数沿流动路径减小。

为了更好阐述伯努利方程的适用前提,首先引入流线和旋度的介绍。

流线

  1. 定义: 流线是在某一瞬间,流场中一系列的假想曲线。这些曲线上任意一点的切线方向都与该点流体的速度矢量方向完全一致。
    性质1 - 无法穿越性: 流体上的所有质点只能沿着流线运动,不能横穿。流线束构成了像管道一样的“流管”。
    性质2 - 密度反映速度性: 根据质量守恒,流线密集处截面小,流速快;流线稀疏处截面大,流速慢。

  2. 流线的数学计算与求解:流线的定义“切线方向与速度矢量方向平行”直接导出了其数学方程。若速度场为\vec{V} = (u, v, w),流线上微小位移为d\vec{r} = (dx, dy, dz),则二者平行意味着各分量成正比:

\frac{dx}{u} = \frac{dy}{v} = \frac{dz}{w}

旋度

流线只描绘了运动的路径,而旋度则揭示了流体微团在前进时自身是否在转动。这是决定伯努利方程应用范围的终极判据。

  1. 直观理解:小桨轮测试
    有旋流动 : 想象中的微小桨轮在随波逐流时自身也在旋转。
    无旋流动: 桨轮只平移,不自转。(重要:流线弯曲不等于有旋!)

  2. 精确判定:旋度矢量 (\vec{\omega})
    旋度是速度场\vec{V} = (u, v, w) 的旋量,其三个分量由速度的偏导数组合而成:

    \omega_x = \frac{\partial w}{\partial y} - \frac{\partial v}{\partial z}\\ \omega_y = \frac{\partial u}{\partial z} - \frac{\partial w}{\partial x}\\ \omega_z = \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y}

判定准则: 若在流场某区域内处处都有\vec{\omega} = \vec{0},则为无旋流场。否则为有旋流场。

3.旋度矢量的深刻物理意义:

方向 - 旋转轴: 旋度矢量\vec{\omega} 的方向指出了流体微团旋转轴的方向(遵循右手定则:大拇指指向\vec{\omega},四指弯曲方向为旋转方向)。例如,\omega_x \neq 0意味着流体微团正在绕着一个平行于x轴的轴旋转。

大小 - 角速度的两倍: 旋度矢量的模长|\vec{\omega}| 是流体微团真实角速度|\vec{\Omega}| 的两倍。即:|\vec{\omega}| = 2|\vec{\Omega}|

证明: 考虑一个以角速度\vec{\Omega} 做刚性旋转的流场,其速度场为\vec{V} = \vec{\Omega} \times \vec{r}。将该速度场的各分量u = \Omega_y z - \Omega_z y,v = \Omega_z x - \Omega_x z,w = \Omega_x y - \Omega_y x 代入旋度计算公式,将得到\omega_x = 2\Omega_x,\omega_y = 2\Omega_y,\omega_z = 2\Omega_z。因此,\vec{\omega} = 2\vec{\Omega}

伯努利方程的应用条件

  1. 第一步:检查基本前提。
    确认流动是稳态、不可压、无粘的。

  2. 第二步:判断流动性质,确定应用范围。

  3. 规则 1:一般情况 (有旋流动 |\vec{\omega} \neq \vec{0})
    规则: 必须严格沿同一条流线 选取两个点应用伯努利方程。
    原因: 每一条流线有其独立的伯努利常数,不同流线上的常数通常不相等。
    p_1 + \frac{1}{2}\rho V_1^2 + \rho gz_1 = p_2 + \frac{1}{2}\rho V_2^2 + \rho gz_2 = C_{\text{沿某特定流线}}

  4. 规则 2:特殊情况 (无旋流动 |\vec{\omega} = \vec{0})
    规则: 可以在流场中任意选取两个点 应用伯努利方程,无需考虑流线。
    原因: 在无旋流场中,伯努利常数在整个流场中是全局唯一的。
    p_A + \frac{1}{2}\rho V_A^2 + \rho gz_A = p_B + \frac{1}{2}\rho V_B^2 + \rho gz_B = C_{\text{全场唯一}}

总之

先保证流动“稳、不压、无粘”;再判断是否“无旋”:无旋是特权,可跨线通用;有旋是常规,须在线内专用。

对于气体而言,当整个气体的速度较大时,在0.3-0.7马赫的速度区间内,可以采用Prandtl-Glauert法则(可压缩流体力学方程)。

此公式是流体力学中最为重要,也是最为出名的方程,其是牛顿第二定律在整个流体力学领域内的一个迁移。

将整个流体的加速度进行分解到x,y,z轴上,于是得到:

a_x = \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + w\frac{\partial u}{\partial z}\\ a_y = \frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} + w\frac{\partial v}{\partial z}\\ a_z = \frac{\partial w}{\partial t} + u\frac{\partial w}{\partial x} + v\frac{\partial w}{\partial y} + w\frac{\partial w}{\partial z}

在这个公式中,(u,v,w)分别表示了在这个点的速度情况,接着对这个质元受到的力进行修正,其受到来自流体压力的梯度力,自身的重力以及由于前后速度差而带来的粘性力。

粘性力的说明:粘性力来自于速度差产生的剪切力({\mu}\frac{\partial{{\tau}_{u}}}{{\partial}y}),由于微元两面的剪切力大小不同,存在差值,计算其差值并采用牛顿莱布尼兹公式进行化简,得到:其对应的粘性合力大小为{\mu}\frac{\partial^2u}{\partial{y^2}},对于x,z方向上的粘性合力的计算方法雷同,于是得到x,y,z方向上的N-S方程:

\rho \left( \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + w\frac{\partial u}{\partial z} \right) = -\frac{\partial p}{\partial x} + \mu \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} \right) + \rho g_x\\ \rho \left( \frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} + w\frac{\partial v}{\partial z} \right) = -\frac{\partial p}{\partial y} + \mu \left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} + \frac{\partial^2 v}{\partial z^2} \right) + \rho g_y\\ \rho \left( \frac{\partial w}{\partial t} + u\frac{\partial w}{\partial x} + v\frac{\partial w}{\partial y} + w\frac{\partial w}{\partial z} \right) = -\frac{\partial p}{\partial z} + \mu \left( \frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2} + \frac{\partial^2 w}{\partial z^2} \right) + \rho g_z

由于上述的方程里含有了4个变量(u,v,w,p),因此求解起来必须有4个方程,考虑到对于不可压缩流体,还有连续性方程,因此得到了针对不可压缩流体的四个控制方程:

\left\{ \begin{aligned} % x-momentum \rho \left( \frac{\partial u}{\partial t} + u\frac{\partial u}{\partial x} + v\frac{\partial u}{\partial y} + w\frac{\partial u}{\partial z} \right) &= -\frac{\partial p}{\partial x} + \mu \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} + \frac{\partial^2 u}{\partial z^2} \right) + \rho g_x \\ \\ % y-momentum \rho \left( \frac{\partial v}{\partial t} + u\frac{\partial v}{\partial x} + v\frac{\partial v}{\partial y} + w\frac{\partial v}{\partial z} \right) &= -\frac{\partial p}{\partial y} + \mu \left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2} + \frac{\partial^2 v}{\partial z^2} \right) + \rho g_y \\ \\ % z-momentum \rho \left( \frac{\partial w}{\partial t} + u\frac{\partial w}{\partial x} + v\frac{\partial w}{\partial y} + w\frac{\partial w}{\partial z} \right) &= -\frac{\partial p}{\partial z} + \mu \left( \frac{\partial^2 w}{\partial x^2} + \frac{\partial^2 w}{\partial y^2} + \frac{\partial^2 w}{\partial z^2} \right) + \rho g_z \\ \\ % Continuity \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} &= 0 \end{aligned} \right.

上述的所有N-S方程都是仅针对于不可压缩流体而言的,对于可压缩流体,其计算过程将会更加复杂,因为涉及到动力粘度\mu的修正,以及多个变量的引入意味着唯一求解需要更多的方程,对于可压缩流体,其NS方程如下所示:

由于上述方程组中引入了(u,v,w,p,ρ,\mu)等变量,仅仅依靠上述的四个方程是无法求解得到其精准解的,因此还需要引入几个方程,理想气体状态方程和空气的萨斯兰定律。

N-S方程的几种简化形式

  • 哈根-泊肃叶简化(圆筒内的层流运动)

  • 建立圆柱形极坐标系,在这个坐标系下原拉普拉斯方程现变化为方程:

    \nabla^2 T = \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial T}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 T}{\partial \theta^2} + \frac{\partial^2 T}{\partial z^2}
  • 基本假设:
    1.该流体是稳定流,即\frac{\partial{u}}{\partial{t}}=0,
    2.轴对称性:即在圆筒内这个流体只能朝着一个方向进行运动,向别的方向不能移动,即有u=u(r),且满足v=w=0
    3.充分发展性:即这个流体在运动的过程中,在运动方向的截面上的形状会近似保持不变的状态。即其u对x的偏导数\frac{\partial{u}}{\partial{x}} 恒为0,
    4.水平管道内部流体的重力可以忽略不计。(因为重力g是沿着竖直方向上的,而流体的运动方向u恰好与力的方向正交,因此重力忽略不计。)

  • 因此原不可压缩流体的N-S方程修正前后如下所示:

    \underbrace{\frac{\partial u}{\partial t}}_{=0 \text{ (稳定流)}} + u_r \frac{\partial u}{\partial r} + \frac{u_\theta}{r} \frac{\partial u}{\partial \theta} + u \underbrace{\frac{\partial u}{\partial x}}_{=0 \text{ (充分发展)}} = -\frac{1}{\rho}\frac{\partial p}{\partial x} + \nu \left[ \frac{1}{r}\frac{\partial}{\partial r}\left(r\frac{\partial u}{\partial r}\right) + \frac{1}{r^2}\frac{\partial^2 u}{\partial \theta^2} + \underbrace{\frac{\partial^2 u}{\partial x^2}}_{=0 \text{ (充分发展)}} \right]

修正之后为:

\frac{dp}{dx} = \mu \frac{1}{r}\frac{d}{dr}\left(r\frac{du}{dr}\right)

由此原来一个极其复杂的PDE就被转化成了一个极易求解的常微分方程,求解这个微分方程得到其解析解(速度分布函数)为:

u(r) = \frac{\Delta p}{4\mu L}(R^2 - r^2)

下面计算其流量:

dQ = u(r) \cdot 2\pi r \, dr\\ Q = \int_{0}^{R} u(r) \cdot 2\pi r \, dr

最终得到哈根-泊肃叶方程:

Q = \frac{\pi \Delta p R^4}{8\mu L}
  • 库埃特流动
    适用于平板之间的剪切运动,假设:两平行板件(此两平行板都是可以无限延伸的)是稳定的层流,且上板是匀速直线运动状态,而下板是静止的一种层流状态。
    库埃特流动相较于泊肃叶流动,其特点在于驱动力不同。泊肃叶运动的驱动力在于两端的压力差,而彼此之间是不存在剪切力的,而库埃特流动的驱动力则是类似多米勒骨牌,通过上一层的移动拖动下一层的运动,因此库埃特流动的驱动力源自上板的剪切运动。

  • 根据库埃特运动的以上特征,我们可以归纳出这种流动具有如下性质:

  • 稳定流 : \frac{\partial u}{\partial t} = 0

  • 一维流动: 流体只在x方向运动,没有y方向的速度v=0

  • 充分发展性: 平板是可以无限延伸的。流动在x方向上是充分发展的,意味着速度分布在任何x位置都相同\frac{\partial u}{\partial x} = 0

  • 无压力梯度: 流动仅由板的运动引起,没有外加的压力驱动\frac{\partial p}{\partial x} = 0

将其带入x方向上的N-S方程\rho \left( \frac{\partial u}{\partial t} + u \frac{\partial u}{\partial x} + v \frac{\partial u}{\partial y} \right) = -\frac{\partial p}{\partial x} + \mu \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} \right) 得到库埃特流动的常微分方程:

\frac{d^2 u}{d y^2} = 0

结合初值条件,求解得到这个方程的通解是u(y) = U \frac{y}{h} ,结合牛顿粘性定律得到各个位置所对应的牛顿剪切力等于:

\tau = \mu \frac{du}{dy} = \mu \frac{U}{h}

由于库埃特流动和泊肃叶流动的驱动力不同,二者完全是毫不相干的,因此他们在流体力学的各个属性上也存在着天壤之别,经过计算前者的速度与位移成正比,而后者与位移则是抛物线关系,如下图所示:

由于库埃特流动的驱动力与泊肃叶完全不同,当库埃特流动的场景下除了自身存在运动的平板以外,流体两端也存在着压力差的时候,这就变成了库埃特-泊肃叶耦合流动模型(广义泊肃叶流动),针对以上所述的公式进行整合和求解,不难得到在综合考虑多个情况下的速度关系式为:

\frac{d^2 u}{d y^2} = \frac{1}{\mu} \frac{dp}{dx}
u(y) = \underbrace{U \frac{y}{h}}_{\text{纯库埃特流动部分}} + \underbrace{\frac{1}{2\mu} \left(\frac{dp}{dx}\right) (y^2 - hy)}_{\text{纯泊肃叶流动部分}}

上式的解较为复杂,需要借助数值分析。