弹性和流体的一点基础知识

本文$\mathbf{v}, \boldsymbol{v}$混用。

有关弹性力学

1.应变

设想一个弹性体发生了形变。那么原本坐标是$\mathbf{x}$的点现在是$\mathbf{x+u}$。

实际上描述具体的形变方式，比$\mathbf{u}$更方便的是采用应变张量${u}_{ij}$。为了得到它，考虑微元$dl$在形变后的变化：

$dl'^2=dl^2+2\frac{\partial u_i}{\partial x_k}dx_idx_k+\frac{\partial u_i\partial u_i}{\partial x_k\partial x_k}dx_kdx_j$

写成

$dl'^2-dl^2=2u_{ik}dx_idx_j$

丢掉高阶项，强加一个对称性，得到应变张量的常用德行

$u_{ik}=1/2(\partial_ku_i+\partial_iu_k)$

从另一种观点来看，这个量代表了形变前后坐标系的度规的差。

2.应力

现在再考虑弹性体里面的点因形变而受到的力。一片小区域里，这个力的合力必然通过包围着它的面积与外界平衡。这样，$\int F_idV=\int\nabla\cdot(?)dS$，在问号处定义应力张量$\sigma_{ik}$:

$F_i=\partial_k\sigma_{ik}$

在这个地方，著名的知乎用户赵永峰指出可以这么做的原因还包括小形变下近似认为对应形变前后坐标系的$\partial_k$和$\partial_{k'}$可以近似相等。

应力张量的对称性没有那么自然，但是注意到，以上的式子必须对力矩也同样成立，计算

$M_{ik}=\int(x_k\partial_l\sigma_{il}-x_i\partial_l\sigma_{kl})dV=\int\partial_l(\sigma_{il}x_k-\sigma_{kl}x_i)dV-\int(\sigma_{ik}-\sigma_{ki})dV$

为了保证力矩只由边界上的量决定，应力张量必然是对称的了。

一般外力$\mathbf{P}$会引起满足$\sigma_{ik}n_k=P_i$的应力。

3.胡克定律

假设应力做功$\delta R$。

$\int dV\delta R=\int dV\partial_k\sigma_{ik}\delta u_i=-\int dV\sigma_{ik}\partial_k\delta u_i=-\int dV\sigma_{ik}\delta u_{ik}$

第二个等号的地方是因为分部积分时可以把化成的面积分的面取在没有形变的无穷远处（弹性体本身的尺寸即可认为有那么大），从而让那一项为0。

能量可以写成$dE=TdS-dR=TdS+\sigma_{ik}du_{ik}$，记 Helmholtz 自由能为$F$，热力势为$\Phi$，那么

$\sigma_{ik}=\frac{\partial F}{\partial u_{ik} };u_{ik}=-\frac{\partial {\Phi} }{\partial \sigma_{ik} }$

对于自由能可以展开成

$F=F_0+\frac{1}{2}\lambda u_{ii}^2+\mu u_{ik}^2$

没有线性的项，因为自由能最小的时候，它对应变的导数（即应力）应该是0。$\lambda,\mu$是“$\lambda\alpha\mu\epsilon$”常量。只展到两阶是合理的，因为当更高阶的项开始起作用的时候，一般来讲弹性形变已经转为了塑性形变。为了物理上的原因，把应变拆成 bulk 和 shear 两部分：

$u_{ik}=(u_{ik}-\frac{1}{3}\delta_{ik}u_{ll})+\frac{1}{3}\delta_{ik}u_{ll}$

用这种写法表示出自由能($K=\lambda+2\mu/3$)：

$F=F_0+\mu(u_{ik}-\frac{1}{3}\delta_{ik}u_{ll})^2+\frac{1}{2}Ku_{ll}^2$

对应变求导，得到应力的表达式

$\sigma_{ik}=(Ku_{ll}\delta_{ik}+2\mu(u_{ik}-\frac{1}{3}\delta_{ik}u_{ll}))$

这个就是胡克定律。

4.使用杨氏模量和泊松比

一般不用拉梅常量来表示物体的弹性性质，而使用更加易于测量的杨氏模量和泊松比。用一个拉伸杆的例子来定义它们：假设一根圆柱形杆两头受拉力$p=\sigma_{zz}$，则可以知道

$u_{xx}=u_{yy}=(1/9K-1/6\mu)p,u_{zz}=(1/9K+1/3\mu)p$

定义杨氏模量$E=p/u_{zz}$，泊松比$\sigma=-u_{zz}/u_{xx}$。

用它们表示的自由能可以写为

$F=\frac{E}{2+2\sigma}(u_{ik}^2+\frac{\sigma}{1-2\sigma}u_{ll}^2)$

用这两个量把胡克定律再写一遍以便直接用：

$\begin{aligned} \sigma_{i k} &=\frac{E}{1+\sigma}\left(u_{i k}+\frac{\sigma}{1-2 \sigma} u_{l l} \delta_{i k}\right) \\ u_{i k} &=\left[(1+\sigma) \sigma_{i k}-\sigma \sigma_{l l} \delta_{i k}\right] / E \end{aligned}$

$\begin{aligned} \sigma_{x x} &=\frac{E}{(1+\sigma)(1-2 \sigma)}\left[(1-\sigma) u_{x x}+\sigma\left(u_{y y}+u_{z z}\right)\right] \\ \sigma_{y y} &=\frac{E}{(1+\sigma)(1-2 \sigma)}\left[(1-\sigma) u_{y y}+\sigma\left(u_{x x}+u_{z z}\right)\right] \\ \sigma_{z z} &=\frac{E}{(1+\sigma)(1-2 \sigma)}\left[(1-\sigma) u_{z z}+\sigma\left(u_{x x}+u_{y y}\right)\right] \\ \sigma_{x y} &=\frac{E}{1+\sigma} u_{x y}, \sigma_{x z}=\frac{E}{1+\sigma} u_{x z}, \sigma_{y z}=\frac{E}{1+\sigma} u_{y z} \end{aligned}$

$\begin{aligned} u_{x x} &=\frac{1}{E}\left[\sigma_{x x}-\sigma\left(\sigma_{y y}+\sigma_{z z}\right)\right] \\ u_{y y} &=\frac{1}{E}\left[\sigma_{y y}-\sigma\left(\sigma_{x x}+\sigma_{z z}\right)\right] \\ u_{z z} &=\frac{1}{E}\left[\sigma_{z z}-\sigma\left(\sigma_{x x}+\sigma_{y y}\right)\right] \\ u_{x y}=\frac{1+\sigma}{E} \sigma_{x y}, & u_{x z}=\frac{1+\sigma}{E} \sigma_{x z}, u_{y z}=\frac{1+\sigma}{E} \sigma_{y z} \end{aligned}$

5.例子1

当我们求解一个均匀、各向同性的物体的弹性时，写出受力

$\partial_k\sigma_{ik}+\rho g_i=0$

用胡克定律写开了，是

$\frac{E}{2(1+\sigma)} \frac{\partial^{2} u_{i} }{\partial x_{k}^{2} }+\frac{E}{2(1+\sigma)(1-2 \sigma)} \frac{\partial^{2} u_{l} }{\partial x_{i} \partial x_{l} }+\rho g_{i}=0$

或者

$\nabla^2 \mathbf{u}+\frac{1}{1-2 \sigma} \nabla(\nabla\cdot \mathbf{u})=-\rho \mathbf{g}\frac{2(1+\sigma)}{E}$

当造成形变的主要因素不是场，而是只施加于表面的力的时候，式子的右侧是0，而外力只体现在边界条件上。这个时候再取一次$\nabla^2$，得到形变满足 biharmonic equation

$\nabla^2\nabla^2\mathbf{u}=0$

对于平面形变的问题，即一个方向（不妨说是z方向）没形变，这个时候达到平衡，则可以写下

$\frac{\partial \sigma_{x x} }{\partial x}+\frac{\partial \sigma_{x y} }{\partial y}=0, \quad \frac{\partial \sigma_{y x} }{\partial x}+\frac{\partial \sigma_{y y} }{\partial y}=0$

一个常规套路是取一个 stress function $\chi$，规定

$\sigma_{x x}=\partial^{2} \chi / \partial y^{2}, \quad \sigma_{x y}=-\partial^{2} \chi / \partial x \partial y, \quad \sigma_{y y}=\partial^{2} \chi / \partial x^{2}$

因为$u_{zz}$是0，这个时候用胡克定律可以写出来

$\nabla^2\chi=\sigma_{x x}+\sigma_{y y}=E\left(u_{x x}+u_{y y}\right) /(1+\sigma)(1-2 \sigma)\propto\nabla\cdot\mathbf{u}$

这个说明$\chi$也是biharmonic的。把问题转成解一个biharmonic的标量场。

6.例子2

还有一个问题是值得关注的，这就是无限大的平面介质受到压力而平衡时的形变。Boussinesq 在19世纪做过一个技巧性很高的解法，也许与朗道口味相投，于是我看到的也是这一版解法。但是，我并没有信心把这个解法的来源讲清楚。

从这里开始：

$\nabla^2 \mathbf{u}+\frac{1}{1-2 \sigma} \nabla(\nabla\cdot \mathbf{u})=0$

取一个解$\mathbf{u}=\mathbf{f}+\nabla\phi$，其中 f 是 harmonic 的，类似于通解的角色。认为平面是xy平面，那么外力P从z方向来。把$\mathbf{f}$表示成$f_{x}=\partial g_{x} / \partial z, \quad f_{y}=\partial g_{y} / \partial z$，其中g的分量可以知道也是 harmonic 的。现在原方程就成了

$2(1-\sigma) \Delta \phi=-\frac{\partial}{\partial z}\left(\frac{\partial g_{x} }{\partial x}+\frac{\partial g_{y} }{\partial y}+f_{z}\right)$

作者说“容易”看出来解是

$\phi=-\frac{z}{4(1-\sigma)}\left(f_{z}+\frac{\partial g_{x} }{\partial x}+\frac{\partial g_{y} }{\partial y}\right)+\psi$

其中$\psi$是 harmonic 的。由于$\sigma_{iz}=P_i$，并用胡克定律而将$\sigma$表示出来，就是

$\begin{array}{c}{\left[\frac{\partial^{2} g_{x} }{\partial z^{2} }\right]_{z=0}+\left[\frac{\partial}{\partial x}\left\{\frac{1-2 \sigma}{2(1-\sigma)} f_{z}-\frac{1}{2(1-\sigma)}\left(\frac{\partial g_{x} }{\partial x}+\frac{\partial g_{y} }{\partial y}\right)+2 \frac{\partial \psi}{\partial z}\right\}\right]_{z=0} }=-2(1+\sigma)P_x/E \\ {\left[\frac{\partial^{2} g_{y} }{\partial z^{2} }\right]_{z=0}+\left[\frac{\partial}{\partial y}\left\{\frac{1-2 \sigma}{2(1-\sigma)} f_{z}-\frac{1}{2(1-\sigma)}\left(\frac{\partial g_{x} }{\partial x}+\frac{\partial g_{y} }{\partial y}\right)+2 \frac{\partial \psi}{\partial z}\right\}\right]_{z=0} }=-2(1+\sigma)P_y/E \\ {\left[\frac{\partial}{\partial z}\left\{f_{z}-\left(\frac{\partial g_{x} }{\partial x}+\frac{\partial g_{y} }{\partial y}\right)+2 \frac{\partial \psi}{\partial z}\right\}\right]_{z=0}=-2(1+\sigma) P_{z} / E}\end{array}$

出于某种机械降神的考虑，直接令前两式最长的两段为0：

$(1-2 \sigma) f_{z}-\left(\frac{\partial g_{x} }{\partial x}+\frac{\partial g_{y} }{\partial y}\right)+4(1-\sigma) \frac{\partial \psi}{\partial z}=0$

这样以上四式可解四个 harmonic function：$g_{ \{x,y\} }, f_z, \psi$。

考虑单点的压力下造成的形变满足

$u_{i}=G_{i k}(x, y, z) F_{k}$

然后对一片压力积分。对于在无穷远处为0的 harmonic function f，如果约束它在每一点的$\partial_zf$，就可以把它的形式写成

$f(x, y, z)=-\frac{1}{2 \pi} \iint\left[\frac{\partial f\left(x^{\prime}, y^{\prime}, z\right)}{\partial z}\right]_{z=0} \frac{\mathrm{d} x^{\prime} \mathrm{d} y^{\prime} }{r}$

其中

$r=\sqrt{\{ \left(x-x^{\prime}\right)^{2}+\left(y-y^{\prime}\right)^{2}+z^{2} \} }$

就得到

$\begin{aligned}f_{z}-\left(\frac{\partial g_{x} }{\partial x}+\frac{\partial g_{y} }{\partial y}\right)+2 \frac{\partial \psi}{\partial z}=\frac{1+\sigma}{\pi E} \cdot \frac{F_{z} }{r}\\ \frac{\partial g_{x} }{\partial z}=\frac{1+\sigma}{\pi E} \cdot \frac{F_{x} }{r}\\ \frac{\partial g_{y} }{\partial z}=\frac{1+\sigma}{\pi E} \cdot \frac{F_{y} }{r}\end{aligned}$

从这里就可以较容易地解出形变了。

所以用应变的定义和胡克定律代入到受力平衡中去，定好边界条件，原则上就可以解出整个问题。比较让人头疼的地方有两个：对于稍微复杂一点的情形，一定要谨慎地寻找和确定边界条件；对于几乎所有的情形，都要解很复杂的方程。不过，好在这些方程都是线性的。

有关流体力学

1.连续性方程

流体的密度和速度写成$\rho(x, y, z, t), \mathbf{v}(x, y, z, t)$。这里，速度指的是出现在时空的固定点上的流体微元的速度，而不是某个特定的流体微元的速度。别的量也是一样的。

认为空间中某个区域中流体质量的变化等同于边界上流入和流出的和，得到连续性方程

$\partial_t\rho+\nabla\cdot(\rho\mathbf{v})=0$

2.Euler 方程

由牛顿第二定律，

$-\int\nabla pdV=\int\rho\operatorname{d}_t\mathbf{v}dV$

但是，此处的这个$\operatorname{d}_t\mathbf{v}$实际上是特定微元的速度，所以为了与语境相适应，考虑

$d\mathbf{v}=\partial_t\mathbf{v}dt+\partial_i\mathbf{v}dx_i$

即得到

$\frac{D}{Dt}\mathbf{v}=\partial_t\mathbf{v}+(\mathbf{v}\cdot\nabla)\mathbf{v}=-\nabla p/\rho$

$D/Dt$是协变导数，也叫物质导数。

对于单位质量的焓$w$，有$dw=Tds+(1/\rho )dp$，那么对于等熵的流，$\nabla p/\rho$就写成$\nabla w$。

3.定常流与 Bernoulli 方程

定常流的意思是$\partial_t\mathbf{v}=0$。对 Euler 方程，利用结论 $\nabla v^2=2(\mathbf{v}\times(\nabla\times\mathbf{v})+(\mathbf{v}\cdot\nabla)\mathbf{v})$，写成

$\partial_t\mathbf{v}+0.5\nabla v^2-\mathbf{v}\times(\nabla\times\mathbf{v})=-\nabla w-\mathbf{g}$

考虑到左边第一项是0，而第三项始终与速度方向垂直，故将此方程在流线（方向$\mathbf{l}$）上投影，得到 Bernoulli 方程（重力在z方向）

$\partial_\mathbf{l}(v^2/2+w+gz)=0$

4.能流和动量流

流体的体积能是$\rho \frac{v^{2} }{2}+\rho \varepsilon$。通过利用以上的方程做一些操作，最后可以得到能流密度

$\rho v\left(\frac{v^{2} }{2}+w\right)$

这个地方是质量焓而不是质量内能，原因是还需要考虑压力做的功。

流体的体积动量是$\rho\mathbf{v}$。同样，可以知道动量流密度是

$\Pi_{i k}=p \delta_{i k}+\rho v_{i} v_{k}$

5.势流与速度势

在等熵流中，沿着一条封闭的物质线（注意，它的空间位形并不是固定着的）取速度环量$\Gamma=\oint \boldsymbol{v} \cdot \mathrm{d} \boldsymbol{l}$，其对时间的全微分可以知道是0。因此，等熵流中速度环量守恒。这时，如果来流在无穷远处均匀，那么在全空间都有$\nabla\times\mathbf{v}=0$。这样的流动叫做势流。此时，可以定义速度势$\phi$，令$v=\nabla\phi$。

对于绕给定物体流动的问题，由于在物体表面难以取速度环量，实际上在物体的表面$\nabla\times\mathbf{v}=0$是不一定正确的，在流线上，则体现为在一些地方流线会中断在物体的表面而如同进入了物体的内部一样。像这样的“间断解”破坏了问题的解的唯一性，好在现实中由于粘度的存在而使得人们（或大自然）能够通过考察物体表面上薄薄的一层流体的情况来确定问题的唯一解。

6.不可压缩流

流体的不可压缩意味着其密度是一个常量。对于采用这种近似的判据，在非定常流中需要两个条件：

$v<<c, \tau>>l/c$

$\tau,l$是速度发生显著变化所需要的时间和距离的量级。c是声速。在定常流中，只需要第一个。这两个条件的意义可以这样解释：流体内部的压强（和速度、密度相关）达不到显著调整密度的水平；流体内部的相互作用传播相对流体的流动几乎是瞬时的。

对于不可压缩流，各方程可以简化为

$\begin{aligned}\nabla\cdot\mathbf{v}=0\\ \partial_t(\nabla\times v)=\nabla\times(\mathbf{v}\times(\nabla\times\mathbf{v}))\\\frac{v^{2} }{2}+\frac{p}{\rho}+g z=\mathrm{const}\end{aligned}$

容易知道速度势满足 laplace 方程。

8.黏性流体的N-S方程

流体的黏性的体现，在于其动量输运中不可逆的部分。连续性方程仍然成立，而 Euler 方程需要改一改。用动量流的形式将 Euler 方程写成

$\frac{\partial}{\partial t}\left(\rho v_{i}\right)=-\frac{\partial \Pi_{i k} }{\partial x_{k} }$

在动量流$\Pi$中加入不可逆输运的部分，

$\Pi_{i k}=p \delta_{i k}+\rho v_{i} v_{k}-\sigma_{i k}^{\prime}=-\sigma_{i k}+\rho v_{i} v_{k}$

$\sigma'_{ik},\sigma_{ik}$分别叫做黏性应力张量和应力张量。这里的应力的说法和弹性体当中的应力是相似的。没有黏度的流体不可能有剪切应力，因为两层流体之间可以无障碍地滑动。另一方面，除了剪切应力之外，$\sigma_{ik}'$还包含正比于$\delta_{ik}$的一项，对应了体积黏度。为了确定黏性应力张量的形式，首先应该知道它只和速度对空间的导数有关，因为只要速度均匀则没有这种应力；丢掉高阶的项，则只和一阶导线性相关。出于力矩平衡的考虑，它也是对称的。如此，写出它的一般形式

$\sigma_{i k}^{\prime}=\eta\left(\frac{\partial v_{i} }{\partial x_{k} }+\frac{\partial v_{k} }{\partial x_{i} }-\frac{2}{3} \delta_{i k} \frac{\partial v_{l} }{\partial x_{l} }\right)+\zeta \delta_{i k} \frac{\partial v_{l} }{\partial x_{l} }$

参数是剪切和体积黏度；

Euler 方程就写成

$\rho\left(\frac{\partial v_{i} }{\partial t}+v_{k} \frac{\partial v_{i} }{\partial x_{k} }\right)=-\frac{\partial p}{\partial x_{i} }+\frac{\partial}{\partial x_{k} }\left[\eta\left(\frac{\partial v_{i} }{\partial x_{k} }+\frac{\partial v_{k} }{\partial x_{i} }-\frac{2}{3} \delta_{i k} \frac{\partial v_{l} }{\partial x_{l} }\right)\right]+\frac{\partial}{\partial x_{i} }\left(\zeta \frac{\partial v_{l} }{\partial x_{l} }\right)$

尽管黏度一般是热力学量的函数，但在多数情况下，它们随空间的变化很缓慢，因此将它们提出来。

$\rho\left[\frac{\partial \boldsymbol{v} }{\partial t}+(\boldsymbol{v} \cdot \nabla) \boldsymbol{v}\right]=-\nabla p+\eta \nabla^2 \boldsymbol{v}+\left(\zeta+\frac{\eta}{3}\right) \nabla(\nabla\cdot \boldsymbol{v})$

这就是黏性流体的 Navier-Stokes 方程。朗道的书中列出了球坐标和柱坐标系下NS方程的形式以方便计算。

对于不可压缩的流体，情况会很大的简化掉，具体来说，N-S 方程简化为

$\frac{\partial \boldsymbol{v} }{\partial t}+(\boldsymbol{v} \cdot \nabla) \boldsymbol{v}=-\frac{1}{\rho} \nabla p+\frac{\eta}{\rho} \nabla^2 \boldsymbol{v}$

而应力张量也成了

$\sigma_{i k}=-p \delta_{i k}+\eta\left(\frac{\partial v_{i} }{\partial x_{k} }+\frac{\partial v_{k} }{\partial x_{i} }\right)$

体积黏度实际上是流体突然压缩或膨胀时在热力学上响应的结果。不可压缩流体自然没有体积黏度。

这个时候可以消去压强，方法是对NS方程取旋度。这时候方程就只和速度有关系。然后可以反解出压强。这个时候的做法一般是对NS方程取散度，这样会得到关于压强的泊松方程。

考虑边界条件。首先，在固体的表面，存在无滑移条件$\Delta\mathbf{v}=0$。然后对于受力，有$P_{i}=-\sigma_{i k} n_{k}=p n_{i}-\sigma_{i k}^{\prime} n_{k}$，这n对于流体是朝外的。

9.一个例子

考虑任意静止柱形管道中黏性流体的平行流，即 Poiseuille 流。不妨令管子的方向是x方向。平行流意味着压强的梯度只在x方向；而速度可表示为$(v(y,z),0,0)$。连续性方程自动满足，而NS方程给出$\partial_yp=\partial_zp=0$，$(\partial_{yy}+\partial_{zz})v=\partial_xp/\eta$。第三个方程左边只和yz有关，右边只和x有关，所以它们都只能是常数了。所以压强梯度是线性的，而速度则满足一个二维的泊松方程。当管道截面是圆形时，此泊松方程的解给出 Hagen-Poiseuille 定律。

10.雷诺数

做量纲分析就可以知道，用 Reynolds 数

$R e=\frac{\rho u l}{\eta}=\frac{u l}{\nu}$

可以衡量惯性和黏性的主导地位。其中，$\nu$是运动黏度$\eta/\rho$。

选取这种数的标准是用决定流动的最少数目的特征量凑成无量纲的量。比如在非定常流，还要用到表征非定常性的时间尺度的量$\tau$，这下$u,l,\tau,\nu$能凑出两个无量纲量，故这时候描述流动除了雷诺数还要用到 Strouhal 数

$S r=\frac{u \tau}{l}$

在低雷诺数的定常流下，NS方程简化为

$\eta \nabla^2 v-\nabla p=0$