恒容容器放气的瞬时流量的计算

有时候，你会遇到一个问题，该问题的描述如下：

你有一个已知体积的容器，设容器体积为 $V$ ，里面装有一定压力(初始压力)的气体，如空气或氢气等，设初始压力为 $P_{0}$ ，容器出口连接着一个阀门开关，开关后面接直径为 $d$ 的钢管，钢管出口为一个大气压 $P_{atm}$ 。当阀门瞬间全开时，气体出口的瞬时流量值随时间变化到底是怎么样的呢？

该问题相当于已知气体管道直径 $d$ ，即已知管道横截面积 $A=\frac{\pi d^{2}}{4}$ ，已知气体管道两端的气压差 $\Delta P=P_{0}-P_{atm}$ ，同时知道进口气体总温 $T_{0}$ 为323K，求出口瞬时质量流量 $q_{m}=\rho uA$ 或瞬时体积流量 $q_v=uA$ 随时间 $t$ 的变化关系和曲线，其中 $\rho$ 是气体密度， $u$ 为出口气体流速， $A$ 即为气体流过管道的横截面积。

1. 第一种方法：根据哈根泊谡叶方程

利用理想气体方程： $PV=\frac{m}{M}RT$ (假设放气是等温过程， $T=T_{0}$ )和哈根泊谡叶关系式： $q_v=\frac{\pi r^{4}}{8\mu L}\Delta P,\Delta P=P-P_{atm}$ ， $q_v$ 表示的是体积流量，单位为 $m^{3}/s$ ， $r$ 是管子的半径， $\mu$ 是流体的动力黏度，单位是 $kg/(m\cdot s)$ ， $L$ 是管子的长度，压强 $P$ 的单位为 $Pa$ 。两个方程联立，

$q_{m}=\frac{dm}{dt}=\frac{d(\frac{MV}{R}PT^{-1})}{dt}=\frac{MV}{RT}\frac{dP}{dt}-\frac{MVP}{RT^{2}}\frac{dT}{dt}$ ，考虑等温过程，有 $q_{m}=\frac{MV}{RT}\frac{dP}{dt}$ ，也就是说， $q_{m}$ 的变化仅与容器内压力的变化有关。进一步，根据 $q_{v}=uA=\frac{q_{m}}{\rho }=\frac{1}{\rho }\frac{dm}{dt}$ ，假设密度 $\rho$ 不变，有： $q_{v}=\frac{MV}{\rho RT}\frac{dP}{dt}=-\frac{\pi r^{4}}{8\mu L}(P-P_{atm})\Leftrightarrow \int_{P_{0}}^{P}\frac{1}{P-P_{atm}}dP=\int_{0}^{t}-\frac{\pi r^{4}}{8\mu L}\frac{\rho RT}{MV}dt$ ，积分后得到： $ln(\frac{P-P_{atm}}{P_0-P_{atm}})=-ABt,A=\frac{\pi r^{4}}{8\mu L},B=\frac{\rho RT}{MV}$ ，有 $P-P_{atm}=(P_0-P_{atm})e^{-ABt}$ ，利用该关系式，得到 $P$ 随时间 $t$ 的关系如下图所示，为一指数函数形式，而且可以通过积分，得到积分总流量为 $7.45L$ ，根据 $\Delta PV=\Delta mR_gT,\Delta m=\frac{\Delta PV}{R_gT}\Leftrightarrow \Delta q=\frac{\Delta m}{\rho }=\frac{\Delta PV}{\rho R_gT}=7.46L$ ，可见积分与差分得出的总流量非常接近。

若气体密度 $\rho$ 不是常数，则根据 $PM=\rho R T$ ，有 $q_{v}=\frac{RT}{PM}\frac{MV}{RT}\frac{dP}{dt}=\frac{V}{P}\frac{dP}{dt}$ ，进一步有：

$\frac{V}{P}\frac{dP}{dt}=-A(P-P_{atm})\Leftrightarrow \frac{1}{P(P-P_{atm})}\frac{dP}{dt}=-\frac{A}{V}\Leftrightarrow (\frac{1}{P}-\frac{1}{P-P_{atm}})dP=ABdt,A=\frac{\pi r^{4}}{8\mu L},B=\frac{P_{atm}}{V}$ ，进一步积分，得到 $ln(\frac{P}{P_{atm}})-ln(\frac{P-P_{atm}}{P_{0}-P_{atm}})=ABt$ ，得到的曲线如下图所示。

通过数值计算，时间小量取 $\Delta t=10^{-6}s$ ，得到的质量流量曲线如下图，积分得到总流过的质量为 $\Delta m=0.64g$ ，与 $\Delta m=\frac{\Delta PVM}{RT}=0.669g$ 有差异，这是因为时间小量 $\Delta t$ 不够小导致的，当你取 $\Delta t=10^{-7}s$ 时，积分得到总流过的质量为 $\Delta m=0.668g$ ，此时就已经与 $0.669g$ 非常接近了。

以密度 $\rho$ 不变的解法为例，一开始的瞬时流量值非常离谱，可以去到 $4093105(L/min)$ ，根据 $\frac{dm}{dt}=\rho \cdot u\cdot A\Leftrightarrow \frac{dq}{dt}=\frac{d(m/\rho )}{dt}=u\cdot A$ ， $A=3.167\times10^{-5} m^{2}$ 可以知道出口流体平均速度 $u=2154092m/s$ ，光速是 $u_c=299792458m/s$ ，出口速度已经达到 $0.007$ 倍的光速，也超过空气声速 $\approx 314m/s$ ，妥妥是一个超音速流，而且放气过程时间非常短，不超过 $0.005s$ 。经过大量的资料查询，该结果似乎与实际测试不符。

附：关于哈根泊谡叶关系式的推导，见下图。

2. 第二种方法：根据气体动力学推算

为什么第一种方法就不符合实际呢？前人发现，收缩的出口在气流流速加速到马赫数1时，即 $Ma=1$ 时，气体流速达到上限。

假设排气过程与气体管道壁面的换热忽略不计，即壁面是绝热的，气体流体是一个准稳态问题，排气口相当于是收缩，没有扩张，根据气体动力学可知，出口气体流速只能加速到1马赫数，即 $Ma=1$ 。根据总静温关系式 $\frac{T_{0}}{T_{b}}=1+\frac{k-1}{2}Ma^{2},T_0=323K$ ，得知 $T_b\approx 269.17K$ 。再根据马赫数定义式 $Ma^{2}=\frac{v^{2}}{c^{2}}=\frac{v^{2}}{kR_gT}$ ，这里 $k$ 是气体比热容比，定义为定压比热 $C_p$ 与定容比热 $C_v$ 之比，变换后有 $v_b=Ma\cdot \sqrt{kR_gT_b}$ ， $k=1.4$ ， $T_b=269.17K$ ，比气体常数 $R_g$ 为： $R_g=\frac{R=8.314J/(mol\cdot K)}{0.002kg/mol}=4157J/(kg\cdot K)$ ，得到氢气气体流速 $v_b\approx 1295.22m/s$ 。

根据 $\frac{dm}{dt}=\rho \cdot u\cdot A$ ， $P=\rho R_gT$ , $\frac{T_{0}}{T_{b}}=1+\frac{k-1}{2}Ma^{2},T_0=323K$ ， $Ma=1$ ， $\frac{P}{P_b}=(1+\frac{k-1}{2}Ma^{2})^{\frac{k}{k-1}}$ ， $u=Ma\cdot \sqrt{kR_gT_b}$ ， $P_{b}$ 为一个大气压。在 $Ma=1$ 的壅塞流阶段，可解得 $ln\frac{P}{P_0}=\frac{\sqrt{1.2kR_gT_0}}{V}MaAt\Leftrightarrow P=f(t)=P_0e^{-47.734t}$ 。这阶段，理解为流速 $u$ 不变， $P$ 变化导致的 $\rho$ 变化，瞬时质量流量也会随之变化，但体积流量 $q_v=u\cdot A$ 不变。如下图所示，绿色曲线是瞬时流量，紫色曲线是体积流量，绿色部分面积是积分得到的总质量流量，通过积分得到壅塞流下的总质量流量为 $0.602g$ ，换算成密度为 $0.0899kg/m^{3}$ 的体积流量为 $6.69L$ 。

后面非壅塞流状态下的亚声速流，原则上也是利用 $\frac{dm}{dt}=\rho \cdot u\cdot A$ ， $P=\rho R_gT_b$ ， $u=Ma\cdot \sqrt{kR_gT_b}$ ， $\frac{T_{0}}{T_{b}}=1+\frac{k-1}{2}Ma^{2},T_0=323K$ ， $\frac{P}{P_b}=(1+\frac{k-1}{2}Ma^{2})^{\frac{k}{k-1}}$ ，这5个式子得到 $\frac{dm}{dt}=f(P)$ 的关系，我用欧拉法获得解析解的近似值，得到后续的流量曲线，具体步骤是，知道压力初始条件 $P=P_{Ma=1}\approx 191801Pa$ ，初始瞬时流量为 $\frac{dm}{dt}|_{Ma=1,P=P_{Ma=1}}$ ，也就是等于壅塞流状态下最后一刻时间的流量，然后利用瞬时流量乘以时间小量，得到 $\Delta m$ ，再利用关系式 $\frac{dm}{dt}=f(P)$ ，得到 $P$ 的变化量，然后计算马赫数 $Ma$ 、速度 $u$ ，温度 $T_b$ 等参数，不断进行迭代计算，当 $P/P_{b}\approx 1$ 时结束迭代。如下图中绿色的质量流量曲线和紫色的体积流量曲线，通过积分面积算得亚声速流下总质量流量为 $0.0673g$ ，换算成密度为 $0.0899kg/m^{3}$ 的体积流量为 $0.749L$ ，因此放氢整个过程总质量流量为 $0.602g+0.067g=0.669g$ ，与 $\Delta PV=\Delta mR_gT,\Delta m=\frac{\Delta PV}{R_gT}\Leftrightarrow \Delta m=0.669g$ 算出来的基本一致。整个过程的总体积流量为 $6.69L+0.749L\approx 7.44L$ 。

3. 第三种方法：绝热小孔自由放气模型

根据伯努利方程得到气体的能量方程 $h+\frac{1}{2}u^{2}=const$ ，根据能量守恒定律，将上述方程应用于小孔，得到小孔上游和下游状态参数的关系： $h_{1}+\frac{1}{2}u_{1}^{2}=h_{2}+\frac{1}{2}u_{2}^{2}$ ，容腔内气体近似保持静止，即小孔上游速度 $u_{1}=0$ 。根据焓值的定义，上式可变为 $u_{2}=\sqrt{2C_{p}(T_{1}-T_{2})}$ ，其中 $C_{p}$ 为气体的恒压热容， $T_1$ 为上游热力学温度， $T_2$ 为下游热力学温度。

由于气流流经小孔时，与管壁接触面小，流动快，可近似认为气体流动过程是绝热过程。将绝热方程式 $PT^{\frac{k}{1-k}}=const$ ， $P=\rho R_gT,R_g=C_p-C_v$ ， $k=\frac{C_p}{C_v}$ ， $k$ 是比热容比， $C_p$ 是定压热容， $C_v$ 是定容热容，以及将 $h=C_pT$ 代入方程 $u_{2}=\sqrt{2C_{p}(T_{1}-T_{2})}$ ，可得 $u_{2}=\sqrt{2C_{p}T_{1}(1-\frac{T_{2}}{T_{1}})}$ ，根据理想气体方程式 $P=\rho R_gT$ ，可得到以下关系：

$P_1=\rho _1R_gT_1\Leftrightarrow \frac{P_1}{\rho _1}=(C_p-C_v)T_1=C_pT_1(1-\frac{C_v}{C_p})=C_pT_1(1-\frac{1}{k})\Leftrightarrow C_pT_1=\frac{k}{k-1}\cdot \frac{P_1}{\rho _1}$

$P_{1}T_{1}^{\frac{k}{1-k}}=P_{2}T_{2}^{\frac{k}{1-k}}\Leftrightarrow\frac{P_2}{P_1}=(\frac{T_1}{T_2})^{\frac{k}{1-k}} \Leftrightarrow \frac{T_2}{T_1}=(\frac{P_2}{P_1})^{\frac{k-1}{k}}$ ，因此可进一步得到：

$u_{2}=\sqrt{\frac{2k}{k-1}\frac{P_1}{\rho _1}[1-(\frac{P_2}{P_1})^{\frac{k-1}{k}}]}$

根据质量流量的定义，有 $q_m=\rho uA$ ，这里由于计算的是小孔的气体流速，因此密度 $\rho =\rho _2$ ，速度即为 $u_2$ ， $A$ 为气体流过小孔管的横截面积。将绝热关系式 $\frac{P}{\rho ^{k}}=const$ ，和 $P=\rho RT$ 代入质量流量定义式，首先有 $\frac{P_1}{\rho_1^{k}}=\frac{P_2}{\rho_2^{k}}\Leftrightarrow \rho_2=\rho_1(\frac{P_2}{P_1})^{\frac{1}{k}}$ ， $P_1=\rho _1R_gT_1\Leftrightarrow \rho _1=\frac{P_1}{R_gT_1}$ ，因此有 $\rho _2=\frac{P_1}{R_gT_1}\cdot (\frac{P_2}{P_1})^{\frac{1}{k}}$ ，所以有 $q_m=\frac{P_1}{R_gT_1}\cdot (\frac{P_2}{P_1})^{\frac{1}{k}}\sqrt{\frac{2k}{k-1}\frac{P_1}{\rho _1}[1-(\frac{P_2}{P_1})^{\frac{k-1}{k}}]} \cdot A$ ，化简有：

$q_m=A\cdot P_{1}\sqrt{\frac{2k}{k-1}\frac{1}{R_gT_{1}}[(\frac{P_2}{P_1})^{\frac{2}{k}}-(\frac{P_2}{P_1})^{\frac{k+1}{k}}]}$ ，其中 $R_g=R/M$ ， $R=8.314$ 气体常数， $M$ 是气体的摩尔质量。如空气的气体常数 $R_g=286.7J/(kg\cdot K)$ 。

$q_m=A\cdot P_{1}\sqrt{\frac{2k}{k-1}\frac{1}{R_gT_{1}}[(\frac{P_2}{P_1})^{\frac{2}{k}}-(\frac{P_2}{P_1})^{\frac{k+1}{k}}]}$ ，定义流量函数 $\phi =\sqrt{\frac{2k}{R_g(k-1)}[(\frac{P_2}{P_1})^{\frac{2}{k}}-(\frac{P_2}{P_1})^{\frac{k+1}{k}}]}$ ，则有： $q_m=A\cdot P_{1}\sqrt{\frac{1}{T_1}} \cdot \phi$ ，将 $\frac{P_2}{P_1}$ 看成自变量， $\phi$ 是因变量，研究 $\phi$ 随 $\frac{P_2}{P_1}$ 的变化曲线，有以下图形关系：

通过对 $\phi$ 进行求导，令 $x=\frac{P_2}{P_1}$ ， $C=\sqrt{\frac{2k}{R_g(k-1))}}$ ， $\phi =C\sqrt{x^{\frac{2}{k}}-x^{\frac{k+1}{k}}}$ ， $\frac{d\phi }{dx}=C\frac{\frac{2}{k}x^{\frac{2}{k}-1}-\frac{k+1}{k}x^{\frac{k+1}{k}-1}}{2\sqrt{x^{\frac{2}{k}}-x^{\frac{k+1}{k}}}}=0\Leftrightarrow \frac{2}{k}x^{\frac{2-k}{k}}=\frac{k+1}{k}x^{\frac{1}{k}}\Leftrightarrow x^{\frac{1-k}{k}}=\frac{k+1}{2}\Leftrightarrow x=\frac{P_2}{P_1}=(\frac{2}{k+1})^{\frac{k}{k-1}}$

当 $k=1.4$ ， $x=\frac{P_2}{P_1}\approx 0.5283$ ，进一步，选取 $R_g=286.7J/(kg\cdot K)$ ，可以计算得到此时 $\phi\approx 0.04043956 s/m$ 。将 $x=\frac{P_2}{P_1}=(\frac{2}{k+1})^{\frac{k}{k-1}}$ 代入 $\phi =\sqrt{\frac{2k}{R_g(k-1)}[(\frac{P_2}{P_1})^{\frac{2}{k}}-(\frac{P_2}{P_1})^{\frac{k+1}{k}}]}$ ，可化简为： $\phi =\sqrt{\frac{k}{R_g}\cdot \frac{2}{k-1}\cdot \frac{2}{k+1}^{\frac{k+1}{k-1}}\cdot (\frac{2}{k+1}^{\frac{1-k}{k-1}}-1)}=\sqrt{\frac{k}{R_g}\cdot \frac{2}{k-1}\cdot \frac{2}{k+1}^{\frac{k+1}{k-1}}\cdot (\frac{k-1}{2})}=\sqrt{\frac{k}{R_g}( \frac{2}{k+1})^{\frac{k+1}{k-1}}}$ ， $q_m=A\cdot P_{1}\sqrt{\frac{k}{R_gT_1}( \frac{2}{k+1})^{\frac{k+1}{k-1}}}$ 。

根据《气动系统的基础与计算特性》一书所述，雷诺首先对上述 $x=\frac{P_2}{P_1}=(\frac{2}{k+1})^{\frac{k}{k-1}}$ 时， $q_m=A\cdot P_{1}\sqrt{\frac{k}{R_gT_1}( \frac{2}{k+1})^{\frac{k+1}{k-1}}}$ 获得极值的物理现象，解释为：当马赫数等于1的状态时，气流处于声速，下游的气流信息不能向上传递，上游的气流状态不随下游压力的变化而变化，许多学者称该现象为壅塞状态。 $\frac{P_2}{P_1}$ 被称为临界压力比，小于此值时，流量达到饱和，也可称为声速流。因此有如下式子，其图形如下图红实线所示。

$\phi =\left\{\begin{matrix} \sqrt{\frac{k}{R_g}( \frac{2}{k+1})^{\frac{k+1}{k-1}}} \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \; \frac{P_2}{P_1}\leq b=0.5283 \\ \sqrt{\frac{2k}{R_g(k-1)}[(\frac{P_2}{P_1})^{\frac{2}{k}}-(\frac{P_2}{P_1})^{\frac{k+1}{k}}]} \, \, \, \; \frac{P_2}{P_1}>b=0.5283 \end{matrix}\right.$

$q_m=\left\{\begin{matrix} A\cdot P_{1}\sqrt{\frac{k}{R_gT_1}( \frac{2}{k+1})^{\frac{k+1}{k-1}}}\, \, \, \, \,\, \, \, \, \, \, \, \, \, \! \! \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \,\frac{P_2}{P_1}\leq b=0.5283 \\ A\cdot P_{1}\sqrt{\frac{2k}{k-1}\frac{1}{R_gT_{1}}[(\frac{P_2}{P_1})^{\frac{2}{k}}-(\frac{P_2}{P_1})^{\frac{k+1}{k}}]}\, \,\, \, \, \, \, \, \, \, \frac{P_2}{P_1}>b=0.5283\end{matrix}\right.$

工程上，经常用 $\phi_{subsonic}=K_{G}\times 2\sqrt{\frac{P_2}{P_1}(1-\frac{P_2}{P_1}) }\, \, \frac{P_2}{P_1}>b=0.5, K_G=0.04043$ 来代替 $\phi =\sqrt{\frac{2k}{R_g(k-1)}[(\frac{P_2}{P_1})^{\frac{2}{k}}-(\frac{P_2}{P_1})^{\frac{k+1}{k}}]} \, \, \, \; \frac{P_2}{P_1}>b=0.5283$ ，曲线如下图的蓝色虚线所示，从图形上看，两者曲线基本吻合，且当 $\frac{P_2}{P_1}=0.5$ 时， $\phi_{subsonic} \approx 0.04043$ ，因此 $b=0.5$ 。两者曲线如下图所示，两者最大误差为3%。

这样得到小孔的一维等熵流动的质量流量的近似公式：

$q_{m}=\left\{\begin{matrix} AP_{1}\frac{K_{G}}{\sqrt{T_{1}}}\, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \frac{P_2}{P_1}\leqslant b=0.5\\ AP_{1}\frac{K_{G}}{\sqrt{T_{1}}}\times 2\sqrt{\frac{P_2}{P_1}(1-\frac{P_2}{P_1}) }\, \, \, \, \, \, \, \frac{P_2}{P_1}> b=0.5\end{matrix}\right.$ ， $K_G=\sqrt{\frac{k}{R_g}( \frac{2}{k+1})^{\frac{k+1}{k-1}}}$ 。

由于 $q_m=\frac{dm}{dt}, PV=mR_gT$ ，可得 $q_m=\frac{dm}{dt}=\frac{d(\frac{PV}{R_gT)})}{dt}$ ，当 $V$ 和 $T_1$ 恒定时，有：

$q_m=\frac{dm}{dt}=\frac{d(\frac{PV}{R_gT_1)})}{dt}=\frac{V}{R_gT_1}\frac{dP}{dt}$

$q_m=\frac{V}{R_gT_1}\frac{dP}{dt}=\left\{\begin{matrix} A\cdot P_{1}\sqrt{\frac{1}{T_1}}\sqrt{\frac{k}{R_g}( \frac{2}{k+1})^{\frac{k+1}{k-1}}} \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \; \frac{P_2}{P_1}\leq b=0.5283 \\ A\cdot P_{1}\sqrt{\frac{1}{T_1}}\sqrt{\frac{2k}{R_g(k-1)}[(\frac{P_2}{P_1})^{\frac{2}{k}}-(\frac{P_2}{P_1})^{\frac{k+1}{k}}]} \, \, \, \; \frac{P_2}{P_1}>b=0.5283 \end{matrix}\right.$ ，有

$\frac{dP}{dt}=\left\{\begin{matrix} \frac{A\cdot P_{1}\sqrt{T_1R_g}}{V}\sqrt{k( \frac{2}{k+1})^{\frac{k+1}{k-1}}} \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \; \frac{P_2}{P_1}\leq b=0.5283 \\ \frac{A\cdot P_{1}\sqrt{T_1R_g}}{V}\sqrt{\frac{2k}{(k-1)}[(\frac{P_2}{P_1})^{\frac{2}{k}}-(\frac{P_2}{P_1})^{\frac{k+1}{k}}]} \, \, \, \; \frac{P_2}{P_1}>b=0.5283 \end{matrix}\right.$ ， $\frac{P_2}{P_1}\leq b=0.5283$ 时，移项得 $p_1dP=(\frac{A\sqrt{T_1R_g}}{V}\sqrt{k( \frac{2}{k+1})^{\frac{k+1}{k-1}}})dt$ ，积分得 $\int_{P_0}^{P}P_1dP=\int_{0}^{t}(\frac{A\sqrt{T_1R_g}}{V}\sqrt{k( \frac{2}{k+1})^{\frac{k+1}{k-1}}})dt\Leftrightarrow P=P_0e^{(\frac{A\sqrt{T_1R_g}}{V}\sqrt{k( \frac{2}{k+1})^{\frac{k+1}{k-1}}})\cdot t}$ ，

也就是说，壅塞状态下(声速流)质量流量曲线是指数形式下降的。亚声速流下的质量流量计算则比较复杂，根据状态方程的微分形式，并按绝热过程处理，有以下公式：

$\left\{\begin{matrix} \frac{dP}{dt}=k\frac{q_mR_gT}{V}\, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \, \\ \frac{dT}{dt}=(k-1)\frac{q_mR_gT^{2}}{PV} \end{matrix}\right.$ ，其中将质量流量公式 $q_m=A\cdot P_{1}\sqrt{\frac{2k}{k-1}\frac{1}{R_gT_{1}}[(\frac{P_2}{P_1})^{\frac{2}{k}}-(\frac{P_2}{P_1})^{\frac{k+1}{k}}]}$ 代入，使用差分形式，取时间微小量为 $\Delta t=10^{-5}s$ ，取 $R_g=4157J/(kg\cdot K)$ ， $A=\frac{\pi d^{2}}{4}$ ，取 $d=6.35mm$ ， $V=1L$ ， $T_0=323K$ ，并进行数值计算，得到的质量流量曲线初看真的很像一条直线，但细看还是有点往上凸的。数值计算结果表明，亚声速流时间约为0.0287s，气体从初始323K降低至最后的269.17K。

若按简化后公式 $\phi_{subsonic}=K_G\times 2\sqrt{\frac{P_2}{P_1}(1-\frac{P_2}{P_1}) }\, \, \frac{P_2}{P_1}>b=0.5283$ ，这里 $R_g=4157J/(kg\cdot K)$ ，所以有 $K_G=\sqrt{\frac{k}{R_g}( \frac{2}{k+1})^{\frac{k+1}{k-1}}}=0.01062s/m$ ，根据此式计算质量流量，同样进行数值计算，亚声速流时间约为0.0361s，气体从初始323K降低至最后的257.32K。两者曲线偏差还不小。

下面是压力、温度、压力比、质量流量随时间的变化曲线，该图计算过程如下，先知道初始各状态，如初始压力 $1MPa=1000000Pa$ ，初始温度 $T_0=323K$ ，初始压力比 $P_2/P_1=101325/1000000=0.101325$ 。用流量公式 $q_m=A\cdot P_{1}\sqrt{\frac{k}{R_gT_1}( \frac{2}{k+1})^{\frac{k+1}{k-1}}}$ ，计算当前时刻的质量流量 $q_m$ ，将 $q_m$ 乘以时间小量 $\Delta t=10^{-5}s$ ，得到 $\Delta m$ ，再用理想气体方程 $\Delta PV=\Delta mR_gT$ 计算出 $\Delta P$ ， $T$ 取当前时刻温度。由于这里当成是绝热放气过程，因此容器内气体温度也会变化，按公式 $\frac{dT}{dt}=(k-1)\frac{q_mR_gT^{2}}{PV}$ 的差分形式，将上述计算好的 $q_m$ ，当前温度 $T$ ，当前压力 $P$ ，代入后获得温度变化量 $\Delta T$ ，放气是对外界做功，因此气体温度下降。用当前压力减去 $\Delta P$ ，获得下一时刻压力 $P-\Delta P=P_{next}$ ；同理，获得下一刻温度 $T-\Delta T=T_{next}$ 。这样，壅塞(声速)状态下，各物理量的曲线就出来了。当压力比 $P_2/P_1>0.5283$ 时，变为亚声速流，此式流量公式变为 $q_m=A\cdot P_{1}\sqrt{\frac{2k}{k-1}\frac{1}{R_gT_{1}}[(\frac{P_2}{P_1})^{\frac{2}{k}}-(\frac{P_2}{P_1})^{\frac{k+1}{k}}]}$ ，按同样的数值计算方法，获得亚声速流下的各物理量随时间的变化曲线。

从曲线可以看出，压力曲线是类似指数型下降的，温度指数型下降，压力比有点类似S型曲线，质量流量在声速阶段是一下凹的曲线，亚声速是上凸的曲线，图中中间连接点不太连续的情况是因为数值计算中，受时间小量的限制，计算到临界压力比时，并不会恰好等于0.5283，压力变化会有个台阶，显得不是那么“连续”。

另外，整个放气过程时间为 $0.08747s$ ，对质量流量曲线进行积分，可以知道总流过的质量为 $6.69\times 10^{-4}kg$ ，即 $0.669g$ ，与 $\Delta PV=\Delta mR_gT,\Delta m=\frac{\Delta PV}{R_gT}\Leftrightarrow \Delta m=0.669g$ 算出来的还是一样。