CUMCM A题建模补漏

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也是可以的。

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

几类特殊边界条件:

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

毕奥数的数学定义式为：

Bi= \frac{h L_c}{k}

​

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

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

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

长度为 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。

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

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 三个系数。

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

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

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

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

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

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

\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，因此它是 双曲型 的。

同样，对于一维热传导方程，自变量是 t 和 x：

\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}则是速度梯度，如图所示，一个物体想要发生移动，就必须克服流体对其的水平剪切力:

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

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

一般来说\mu越大，流体越容易做规整运动（层流），\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. 无粘流动 : 忽略流体内部的粘性摩擦。粘性会消耗机械能，（导致系统内部的能量不守恒）使其转化为内能，导致伯努利方程中的常数沿流动路径减小。

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

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

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

判定准则: 若在流场某区域内处处都有\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}

伯努利方程的应用条件

总之

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

对于气体而言，当整个气体的速度较大时，在 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) 等变量，仅仅依靠上述的四个方程是无法求解得到其精准解的，因此还需要引入几个方程，理想气体状态方程和空气的萨斯兰定律。

修正之后为:

\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}

将其带入 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{纯泊肃叶流动部分}}

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