我们将利用 \(P=\rho gh\) 和 \(PV=nRT\) 计算压强微分 \(dP\)，并将伯努利方程中的 \(v^2\) 线性化为 \(v_{old} \cdot v\)，从而构建一个关于流速 \(\mathbf{v}\) 的一阶线性矩阵方程。
1. 变量定义
- 未知量：\(m\) 根管子在下一时刻的流速向量 \(\mathbf{v} = [v_1, v_2, \dots, v_m]^T\)。
- 状态量：
- \(h_i\)：容器 \(i\) 的液位。
- \(n_i\)：容器 \(i\) 的气体摩尔数。
- 几何量：
- \(A_i\)：容器 \(i\) 的截面积。
- \(S_k\)：管子 \(k\) 的截面积。
2. 压强微分推导 \(dP\)
对于任意一根管子 \(k\)，连接容器 \(i\)（尾）和容器 \(j\)（头）。
流动方程基于压强差：\(\Delta P_k = P_j - P_i\)。
我们需要计算由于液体和气体转移引起的压强差变化 \(d(\Delta P_k)\)。
容器 \(i\) 的总压强变化 \(dP_i\)
总压强 \(P_i = P_{gas,i} + \rho g h_i\)。
根据全微分公式：
\[
dP_i = \frac{\partial P_i}{\partial h_i} dh_i + \frac{\partial P_i}{\partial n_i} dn_i
\]
1. 液体部分 \((\rho g h)\)：
\[
\frac{\partial P_i}{\partial h_i} = \rho g
\]
2. 气体部分 \(\big(PV=nRT \Rightarrow P = \dfrac{nRT}{A(H_{max}-h)}\big)\)：
\[
\frac{\partial P_{gas}}{\partial h} = \frac{nRT}{A(H_{max}-h)^2} = \frac{P_{gas}}{H_{gas}}
\]
\[
\frac{\partial P_{gas}}{\partial n} = \frac{RT}{V_{gas}} = \frac{P_{gas}}{n}
\]
综上，容器 \(i\) 的压强微分方程为：
\[
dP_i = \left( \rho g + \frac{P_{gas,i}}{H_{gas,i}} \right) dh_i + \frac{P_{gas,i}}{n_i} dn_i
\]
令 \(K_i = \rho g + \dfrac{P_{gas,i}}{H_{gas,i}}\)（总刚度），则：
\[
dP_i = K_i \cdot dh_i + \frac{P_{gas,i}}{n_i} dn_i
\]
3. 关联状态变化与流速 \(dh, dn \rightarrow v\)
在时间步长 \(\Delta t\) 内，状态的变化由流过管子的流体决定。
假设流速 \(v_k\) 定义为从 \(i\) 流向 \(j\)。
1. 液体转移 \((dh)\)：
容器 \(i\) 失去液体，容器 \(j\) 获得液体。
\[
dh_i = - \frac{\Delta t}{A_i} \sum_{k \in \text{out}(i)} S_k v_k + \frac{\Delta t}{A_i} \sum_{k \in \text{in}(i)} S_k v_k
\]
简写为：
\[
dh_i = \sum_k C_{ik} v_k
\]
（\(C\) 是几何系数矩阵）
2. 气体转移 \((dn)\)：
假设气体随液体流动或通过专门管道流动（此处假设气体摩尔数变化与流速 \(v\) 有关，例如气液两相流或纯气管）：
\[
dn_i = \sum_k D_{ik} v_k
\]
(如果是纯液体管网，\(dn=0\)，此项忽略；如果是气液耦合，此项由气体流速决定)
4. 伯努利方程线性化 (二阶变一阶)
伯努利方程（忽略阻力，考虑动能项）：
\[
P_i + \frac{1}{2} \rho v_i^2 = P_j + \frac{1}{2} \rho v_j^2
\]
在管子 \(k\) 中，我们关注压差驱动流速。简化模型通常写作：
\[
\Delta P_k = \frac{1}{2} \rho v_k^2
\]
线性化 \(v^2 \approx v_{old} \cdot v\)
\[
P_j - P_i = \rho \cdot v_{old,k} \cdot v_k
\]
代入压强微分公式：
\[
(P_j^{old} + dP_j) - (P_i^{old} + dP_i) = \rho \cdot v_{old,k} \cdot v_k
\]
移项，将未知量 \(v\) 放在左边，已知量放在右边：
\[
dP_j - dP_i - \rho \cdot v_{old,k} \cdot v_k = P_i^{old} - P_j^{old}
\]
5. 最终矩阵方程
将 \(dP\) 展开为 \(v\) 的函数：
\[
dP_j - dP_i = \sum_{m} \left( K_j C_{jm} + \frac{P_{gas,j}}{n_j} D_{jm} - K_i C_{im} - \frac{P_{gas,i}}{n_i} D_{im} \right) v_m
\]
令系数矩阵元素为 \(M_{km}\)，则对于每一根管子 \(k\)：
\[
\sum_{m} M_{km} v_m - \rho \cdot v_{old,k} \cdot v_k = \Delta P_{static}
\]
线性方程组 \(\mathbf{A}\mathbf{v} = \mathbf{b}\)：
\[
\begin{bmatrix}
M_{11} - \rho v_{old,1} & M_{12} & \cdots & M_{1m} \\
M_{21} & M_{22} - \rho v_{old,2} & \cdots & M_{2m} \\
\vdots & \vdots & \ddots & \vdots \\
M_{m1} & M_{m2} & \cdots & M_{mm} - \rho v_{old,m}
\end{bmatrix}
\begin{bmatrix}
v_1 \\
v_2 \\
\vdots \\
v_m
\end{bmatrix}
=
\begin{bmatrix}
P_{i,1}^{old} - P_{j,1}^{old} \\
P_{i,2}^{old} - P_{j,2}^{old} \\
\vdots \\
P_{i,m}^{old} - P_{j,m}^{old}
\end{bmatrix}
\]
矩阵元素 \(M_{km}\) 的物理意义：
\[
M_{km}
=
\underbrace{\left( K_j \frac{\partial h_j}{\partial v_m} + \frac{P_{gas,j}}{n_j} \frac{\partial n_j}{\partial v_m} \right)}_{\text{管头 j 的压强反馈}}
-
\underbrace{\left( K_i \frac{\partial h_i}{\partial v_m} + \frac{P_{gas,i}}{n_i} \frac{\partial n_i}{\partial v_m} \right)}_{\text{管尾 i 的压强反馈}}
\]
求解步骤：
- 计算当前 \(P_{gas}\) 和 \(h\)。
- 计算刚度 \(K = \rho g + P_{gas}/H_{gas}\)。
- 组装矩阵 \(\mathbf{A}\) 和向量 \(\mathbf{b}\)。
- 求解 \(\mathbf{v}\)。
- 更新 \(h\) 和 \(n\)。