hanyeah

压强计算，千问给的答案。

没问题，完全按照你的思路。

我们将利用 $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. **气体部分 (**$PV=nRT \Rightarrow P = \frac{nRT}{A(H_{max}-h)}$**)**：
    - 对 $h$ 求导（体积压缩）：$\frac{\partial P_{gas}}{\partial h} = \frac{nRT}{A(H_{max}-h)^2} = \frac{P_{gas}}{H_{gas}}$
    - 对 $n$ 求导（质量增加）：$\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 + \frac{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$
*(注：这里系数调整为 *$\rho$* 以匹配量纲，核心是用 *$v_{old}$* 将方程降为一阶)*

代入压强微分公式：
$(P_j^{old} + dP_j) - (P_i^{old} + dP_i) = \rho \cdot v_{old,k} \cdot v_k$

移项，将未知量 $v$ 放在左边，已知量（$P^{old}$）放在右边：
$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 的压强反馈}}$

**求解步骤：**

1. 计算当前 $P_{gas}$ 和 $h$。
2. 计算刚度 $K = \rho g + P_{gas}/H_{gas}$。
3. 组装矩阵 $\mathbf{A}$ 和向量 $\mathbf{b}$。
4. 求解 $\mathbf{v}$。
5. 更新 $h$ 和 $n$。

新标签页打开。