空间坐标下的守恒方程

质量守恒

考虑空间坐标下的控制体 \(V\), 在 \(\delta t\) 时间内, 流出控制体 \(V\) 的质量为 \[ \mathrm{d} m = \delta_t\int_{\partial V} \rho v\cdot n ~\mathrm{d}S + \mathcal{O}(\delta t^2) \] 控制体在 \(t+\delta t\) 时刻的质量为 \[ m_V(t+\delta t) = \int_{V} \rho(x,t+\delta t) ~\mathrm{d}x \] 质量守恒的表述为: 在控制体的质量改变, 只能通过质量进出控制体引起. 因此有 \[ m_{V}(t + \delta t) = m_{V}(t) - \mathrm{d} m \Rightarrow m_{V}(t + \delta t) - m_{V}(t) + \mathrm{d} m = 0 \] 两边同除 \(\delta t\), 并令 \(\delta t \rightarrow 0\), 就得到 \[ \int_{V} \frac{\partial \rho}{\partial t} ~\mathrm{d}x + \int_{\partial V} \rho v\cdot n ~\mathrm{d}S = 0 \] 再应用散度定理, 就得到 \[ \int_{V} \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho v) ~\mathrm{d}x = 0 \] 又因为 \(V\) 是任意选取的, 所以有 \[ \boxed{ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho v) = 0 } \tag{1} \] 如果材料不可压缩, 那么, 无论是物质坐标下, 还是空间坐标下的密度, 都是恒定不变的常数. 因此, 根据式 (1) 有 \[ J \equiv 1 \] 同时为了保证质量守恒, 又不能改变物质的密度, 因此必须对速度场进行约束: \[ \nabla \cdot v = 0 \]

输运定理

设想一个包裹着一团物质的曲面 \(\partial B\), 当材料变形时, 曲面也在不断地变形. 在当前构型中, 设由该曲面包裹的物质区域为 \(B\). 对该区域的物理量(的质量密度) \(\phi\) 进行积分: \[ \Phi_{B}(t) = \int_{B} \rho \phi ~\mathrm{d} x \] 与此同时, 在当前构型中考虑与 \(B\) 重合的控制体 \(V\), 对上式积分项求偏导数后积分: \[ \int_{V} \frac{\partial (\rho \phi)}{\partial t} ~\mathrm{d}x \] \(\Phi\) 在控制体 \(V\) 中的增加量, 等于物质面 \(B\) 中 \(\Phi\) 的增加量, 再加上流进 \(V\) 的 \(\Phi\), 写成公式为 \[ \int_{V} \frac{\partial (\rho \phi)}{\partial t} ~\mathrm{d}x = \frac{\mathrm{d}~}{\mathrm{d} t} \int_{B} \rho \phi ~\mathrm{d}x - \int_{\partial V} \rho \phi v \cdot n ~\mathrm{d} S \] 应用散度定理, 上式就可以写成 \[ \frac{\mathrm{d}~}{\mathrm{d} t} \int_{B} \rho \phi ~\mathrm{d}x = \int_{V} \frac{\partial (\rho \phi)}{\partial t} + \mathrm{div}(\rho \phi v) ~\mathrm{d} x \tag{2} \] 观察式 (2), 可以看到, 等式左端是物质团 \(B\) 关于时间的变化率, 等式右端是关于空间区域 \(V\) 的积分, 因此, 输运定理将物质团变化率表述为为空间坐标的形式. 需要指出的是, 输运定理并不是守恒定律, 但确实可以通过质量守恒进行化简. 式 (2) 右端项可以展开为 \[ \int_{V} \frac{\partial (\rho \phi)}{\partial t} + \mathrm{div}(\rho \phi v) ~\mathrm{d} x = \int_{V} \rho \frac{\partial \phi}{\partial t} + \rho v \cdot \nabla \phi + \phi \left( \frac{\partial \rho}{\partial t} + \mathrm{div}(\rho v) \right) ~\mathrm{d} x \\ = \int_{V} \rho \left( \frac{\partial \phi}{\partial t} + v \cdot \nabla \phi \right) ~\mathrm{d} x = \int_{V} \rho \frac{\mathrm{D} \phi}{\mathrm{D} t} ~\mathrm{d} x \] 所以, 式 (2) 可以通过质量守恒 (1) 化简为 \[ \boxed{ \frac{\mathrm{d}~}{\mathrm{d} t} \int_{B} \rho \phi ~\mathrm{d}x = \int_{V} \rho \frac{\mathrm{D} \phi}{\mathrm{D} t} ~\mathrm{d} x } \tag{3} \] 可见物质团 \(B\) 积分号外的关于时间的导数, 移动到控制体积分中对物理量 \(\phi\) 的物质导数.

[!NOTE]

应用输运定理推导质量守恒的空间表达形式

围绕在物质曲面 \(\partial B\) 的物理量 \(\phi\) 如果等于恒定的常数 \(1\), 那么 \(\phi\) 在空间或者物质坐标下的表达形式是相同的. 定义关于时间的函数 \(m(t)\) 为 \[ m(t) \triangleq \int_{B} \rho(x,t) ~\mathrm{d}x \] 质量守恒定律声明, 物质团的质量, 无论物质如何变形, 总是恒定不变, 因此 \[ m(t) \equiv m(0) = \int_{B_0} \rho(X,0) ~\mathrm{d}X \] 因此 \[ \frac{\mathrm{d}~}{\mathrm{d} t} \int_{B} \rho ~\mathrm{d}x \equiv 0 \] 代入式 (2) 中得到 \[ \int_{V} \frac{\partial \rho}{\partial t} + \mathrm{div}(\rho v) ~\mathrm{d} x = 0 \]

动量守恒定律

粒子所满足的动量守恒定律为: 粒子动量的变化率, 等于粒子所受的外力: \[ \frac{\mathrm{d}~}{\mathrm{d} t} (mv) = f \] 将粒子换成物质团 \(B\), 就得到 \[ \frac{\mathrm{d}~}{\mathrm{d} t} \int_{B} \rho v ~\mathrm{d}x = f \] 对上式左端应用质量守恒简化后的输运定理 (3) : \[ \frac{\mathrm{d}~}{\mathrm{d} t} \int_{B} \rho v ~\mathrm{d}x = \int_{V} \rho a ~\mathrm{d}x \] 在当前构型下, 选取与 \(B\) 重合的控制体 \(V\), 外力可以表示为两部分: \(V\) 内的体力和 \(\partial V\) 上的面力: \[ f = \int_{V} \rho b ~\mathrm{d}x + \int_{\partial V} \sigma \cdot n ~\mathrm{d}S \] 再应用散度定理, 将面积分转换为体积分, 就得到 \[ \int_{V} \sigma \cdot \nabla + \rho b - \rho a ~\mathrm{d}V = 0 \] 又因为 \(V\) 的任意性, 所以得到空间坐标下的动量守恒定律的微分形式 \[ \boxed{ \sigma \cdot \nabla + \rho b - \rho a = 0 } \tag{4} \]

如果系统满足平衡方程, 那么加速度 \(a=0\), 因此式 (4) 为 \[ \boxed{ \sigma \cdot \nabla + \rho b= 0 } \tag{5} \]

角动量守恒定律

粒子所满足的角动量守恒定律为: 粒子关于坐标系原点 \(O\) 角动量的变化率, 等于粒子所受的外力矩: \[ \frac{\mathrm{d}~}{\mathrm{d} t} (x \times mv) = x \times f \] 我们应用分量形式进行推导, 上述方程的左端为 \[ \frac{\mathrm{d}~}{\mathrm{d} t} \int_{B} \rho x \times v ~\mathrm{d}x = \int_{V} \rho e_{ijk} \frac{\mathrm{D} (x_j v_k)}{\mathrm{D} t} ~\mathrm{d}x = \int_{V} \underbrace{\rho e_{ijk} v_j v_k}_{=0} + \rho e_{ijk} x_j a_k ~\mathrm{d}x = \int_{V} \rho e_{ijk} x_j a_k ~\mathrm{d}x \] 方程的右端为 \[ x \times f = \int_{V} \rho e_{ijk} x_j b_k ~\mathrm{d}x + \int_{\partial V} e_{ijk} x_j \sigma_{kl} n_{l} ~\mathrm{d}S \] 并应用散度定理, 得到 \[ \int_{\partial V} e_{ijk} x_j \sigma_{kl} n_{l} ~\mathrm{d}S = \int_{V} e_{ijk} (x_j \sigma_{kl})_{,l} ~\mathrm{d}x \] 式中的积分项为 \[ e_{ijk} (x_j \sigma_{kl})_{,l} = e_{ijk} (\sigma_{kj} + x_{j} \sigma_{ki,i}) \] 代入到角动量守恒方程中, 整理得到 \[ \int_{V} e_{ijk} x_j (\rho a_k - \rho b_k - \sigma_{ki,i}) ~\mathrm{d}x = \int_{V} e_{ijk} \sigma_{kj} ~\mathrm{d}x \] 根据平衡方程 (4), 上式左端恒等于 0, 所以 \[ \boxed{ e_{ijk} \sigma_{kj} = 0 \Rightarrow \sigma_{jk} = \sigma_{kj} } \tag{6} \]

能量守恒定律

物质的能量由动能 \(K\) 和内能 \(E\) 构成, 定义为 \[ K \triangleq \frac{1}{2} \int_{B} \rho |v|^2 ~\mathrm{d}x, \quad E \triangleq \int_{B} \rho e ~\mathrm{d}x \] 式中, \(|v|^2 = v \cdot v\), \(e\) 为材料内能的质量密度. 物质团的能量来源有两种: 外力做功与其它形式热量的流入流出. 因此能量守恒定律就是在陈述这些能量的平衡关系.

[!NOTE]

如果我们不能够确定内能 \(e\) 或 \(E\) 的形式, 那么能量守恒定律只能是一个文字的陈述, 或者作为内能的定义. 能量守恒只有在确定内能形式 (也就是本构) 的时候才有用.

如果考虑外部能量的来源只有热流 \(h\), 那么能量守恒定律可以表示为 \[ \frac{\mathrm{d}~}{\mathrm{d} t} \int_{B} \rho \left( \frac{1}{2}|v^2| + e \right) ~\mathrm{d}x = \int_{V} \rho b \cdot v ~\mathrm{d}x + \int_{\partial V} (v \cdot \sigma - h) \cdot n ~\mathrm{d}S \] 对上式左端应用输运定理得到 \[ \frac{\mathrm{d}~}{\mathrm{d} t} \int_{B} \rho \left( \frac{1}{2}|v^2| + e \right) ~\mathrm{d}x = \int_{V} \rho a \cdot v + \rho \dot{e} ~\mathrm{d}x \] 再由 \(V\) 的任意性, 得到微分形式的能量守恒定律 \[ \rho a \cdot v + \rho \dot{e} = \rho b \cdot v + (v \cdot \sigma) \cdot \nabla - h \cdot \nabla \] 式中, \[ (v \cdot \sigma) \cdot \nabla = (v_{i} \sigma_{ij})_{,j} = (\sigma \cdot \nabla) \cdot v + \sigma : v \nabla \] 代入得到 \[ [\rho a - \rho b - (\sigma \cdot \nabla)] \cdot v + \rho \dot{e} - \sigma : v \nabla + h \cdot \nabla =0 \] 也即 \[ \boxed{ \rho \dot{e} - \sigma : v \nabla+ h \cdot \nabla = 0 } \tag{7} \]

[!NOTE]

不可压缩牛顿流体的能量守恒方程

流体的本构方程为 \[ \sigma = -p I + 2\mu D \] 热流满足 Fourier 传热定律: \[ h = - \kappa \nabla T \] 流体的内能密度只与温度相关 \[ e = C_m T \] 代入到方程 (7) 中得到 \[ \rho C_m \dot{T} + p \nabla \cdot v= 2\mu D : D + \nabla^2 T \] 又因为不可压缩性, 所以 \(\nabla \cdot v \equiv 0\), 因此 \[ \rho C_m \dot{T} = 2\mu D : D + \nabla^2 T \]

虚功原理

考虑满足平衡方程 \[ \sigma \cdot \nabla + \rho b = 0 \] 的应力场 \(\sigma\), 以及独立的速度场 \(v\), 和相应的变形率张量 \[ D = \frac{1}{2} ( \nabla v + v \nabla ) \] 需要强调的是, 应力场 \(\sigma\) 和速度场 \(v\) 是相互独立的, 没有任何的关联. 在 \(V\) 上对应力和变形率张量的双点积积分得到 \[ \int_{V} \sigma : D ~\mathrm{d} x = \int_{V} \sigma : v\nabla ~\mathrm{d}x = \int_{\partial V} v \cdot \sigma \cdot n ~\mathrm{d}S + \int_{V} \rho b \cdot v ~\mathrm{d}x \] 注意到 \(\sigma \cdot n = t^{(n)}\) 为作用在边界 \(\partial V\) 的面力, 因此 \[ \boxed{ \int_{V} \sigma : D ~\mathrm{d} x = \int_{\partial V} t^{(n)} \cdot v ~\mathrm{d}S + \int_{V} \rho b \cdot v ~\mathrm{d}x } \tag{8} \]

[!NOTE]

小变形下的虚功原理

在小变形假设下, \[ \delta u = v \delta t, \quad \delta \varepsilon \triangleq \frac{1}{2}(\nabla \delta u + \delta u \nabla) = D \delta t \] 因此, 对式 (8) 两边同乘 \(\delta t\), 得到 \[ \int_{V} \sigma : \delta \varepsilon ~\mathrm{d} x = \int_{\partial V} t^{(n)} \cdot \delta u ~\mathrm{d}S + \int_{V} \rho b \cdot \delta u ~\mathrm{d}x \]