工业气体是现代工业基础原料的重要组成部分，随着国内新兴产业如光纤、精细化工、电子技术等的快速发展，用气市场的分布式供应体系不断壮大，成为未来工业气体发展的趋势之一. 据QY Research^[1]针对全球低温液体市场的销售额统计结果，低温液化气体的需求是逐年增长的，预期截至2027年，年复合增长率有望达到2.67%，其中液氮作为低温液体行业最大的细分市场，占到了全球份额的49.94%. 据统计，2018—2022年，中国液氮产量复合增长率达11.51%^[2].

液氮作为一种惰性、无色无臭、无腐蚀性、不可燃且温度极低的液体，可以从空气中比较方便地提纯获取，在医疗、能源开采、科研、地层冻结加固等领域使用频繁. 当下用户使用的液氮来源主要是由气体厂家集中生产后采用低温槽车装载运输，在该过程中，由于液氮的低温属性，漏热导致的质量和能量损失是不可避免的问题. 汪青青等^[3]列举分析了120 t液氮在储运过程中一天的损失率就高达8.25%. 因此，针对科研实验室、医疗美容院、高校教学示范等具有液氮持续使用需求但用量不大的场合，长期以往，运输损耗会让液氮使用成本变得高昂；此外像一些偏远地区例如能源开采的地域，欠缺的交通条件以及使用地点的分散性，会进一步增加成本. 因此，研制新型的小型分布式液氮生产设备具有重要现实意义.

此前，已经有不少学者针对小型氮气液化系统的气体分离技术和液化技术进行了详细的研究. 在气体分离技术方面，大部分研究采用变压吸附（pressure swing adsorption, PSA）技术从空气中捕获高纯氮气，低温精馏^[4-6]则多被用于大型空分系统中；在液化技术方面，以混合工质节流循环（mixed-refrigerant J-T, MRJT）^[7-11]或者回热式制冷机^[12-19]作为冷源是目前使用较多的. Litter^[9]设计的微型氮气发生器采用PSA技术，配合MRJT，在0.8 MPa的进气压力下，可实现0.25 L/h的产量，比功耗为3.76 kW·h/L. Wang等^[10]针对不同气体分离技术对一套小型MRJT氮气液化系统进行了评估，在0.8 MPa进气工况下，采用PSA技术的系统比功耗约为2.68 kW·h/L；采用低温精馏技术的系统比功耗约为4.59 kW·h/L. 郭浩等^[11]设计了一套基于PSA技术的带预冷的MRJT制氮系统，系统整体能耗为1.77 kW·h/L. Caughley等^[14]采用制冷量为 1250 W@77 K的脉管制冷机液化由PSA技术获得的高纯氮气，在0.8 MPa进气压力工况下，产量可达12.1 L/h. 林诗燕^[15]采用PSA以及大冷量斯特林制冷机技术搭建了一套小型氮气液化系统，利用冷机未液化的氮气对气源进行预冷，最佳工况下得到7.7 L/h的液氮产量，不考虑空压机的能耗，系统比功耗约为1.46 kW·h/L.

针对采用较多的PSA技术，王颖丹^[20]对常规均压PSA制氮工艺流程进行了优化与效果测评，优化后的PSA工艺氮提取率为22%. 胡忠宽^[21]从多方面将其与低温精馏和膜分离进行对比，结果表明，面对氮气纯度需求较高的情况，PSA的氮提取率不足低温精馏的二分之一；膜分离的会更低一些并且膜的制造成本相对较高，因此在小型制氮系统中采用膜分离技术的研究不多^[7-8]. 至于液化技术，MRJT容易受制冷剂储量、配比的制约；膨胀制冷循环在小型液化系统中被采用的也较少^[22]，小产量使得系统可利用膨胀功较少，转换回收利用效益较低.

综上，为了解决分布式液氮供应体系运输过程中资源浪费、成本高昂的问题，并且改善目前小型PSA制氮系统氮提取率低、综合能耗高的状况，本研究依据不同的系统运行条件，选取不同的进气方式与压缩空气纯化方式，提出了3种采用斯特林制冷机的低温精馏制氮系统：低压进气-常压精馏系统、高压进气-常压精馏系统、高压进气-加压精馏系统. 基于Aspen HYSYS对3种流程在不同空气入口温度、压缩空气初步预冷温度、制冷机输入冷量条件下进行详细优化模拟，用于确定最佳制氮流程与最优运行工况.

基于大冷量斯特林制冷机的小型低温精馏制氮系统流程原理如图1所示. 采用合适的压缩设备压缩空气；初步预冷后采用高效的空气纯化措施除去原料气中的水和CO₂；此后干燥纯净的压缩空气与精馏柱底部的富氧液态空气充分换热实现精馏系统冷量回收；预冷后的压缩空气在精馏柱内部借助高效散堆填料充分传热传质，在精馏柱顶部富集纯度大于99.5%的氮气；高纯氮气通过气泡泵结构^[23]上方夹层流道进入斯特林制冷机被进一步降温液化得到液氮；液氮返回流入气泡泵夹层结构后，一部分作为产品流出系统，其余在气泡泵中间同心管的底部换热器处受热形成气液两相流，进而借助密度差向上流动实现精馏柱回流液的泵送，由于底部换热器所需的换热量较少，该热量可从预冷后的压缩空气中获取. 本系统工艺方法已申请发明专利，申请号为2023115101875.

基于上述流程原理，依据不同的系统运行条件，具体设计了3种系统流程：

流程1：低压进气-常压精馏流程. 选用鼓风机作为压缩设备，压比仅为1.6，原料气温升较小，可采用自然冷却的初步预冷手段甚至省去该环节. 在低压工况下，难以通过吸附法去除原料气中的CO₂，因此该流程采用冷冻干燥法在预冷换热器内同步实现空气纯化，图1中压缩空气净化设备仅起干燥作用. 纯化后的干燥压缩空气在回收系统冷量实现预冷后进入精馏柱提纯，精馏柱顶部获得的高纯氮气进入制冷机实现液化.

流程2：高压进气-常压精馏流程. 选用空气压缩机作为压缩设备，压比为8.9，原料气温升明显，须借助冷水机组进行初步降温. 而后采用常温吸附法进行空气纯化. 纯化后的干燥压缩空气先在预冷换热器回收系统冷量，而后通过节流的方式进一步降压降温（如图1中 θ₄所在流程路线上），余下流程与流程1保持一致. 本流程按照实际设备结构尺寸设计有如图2所示的三维模型，撬装化后占地面积约为 5 m²，流程具备实际可行性以及高度紧凑性.

流程3：高压进气-加压精馏流程. 纯化后的干燥压缩空气在回收系统冷量后直接进入精馏柱提纯，整个系统全程处于高压操作环境. 其余过程均与流程2保持一致.

3种流程均采用Aspen HYSYS软件进行稳态模拟分析. 软件中强化的P-R方程对低温深冷体系能进行准确的相平衡计算，被广泛用于气体液化流程，很多研究者^[15,24]基于该物性方法进行液化流程模拟，理论结果与实际较吻合，故本研究选用P-R方程作为物性计算方法；将流程内气泡泵结构简化为换热器模型，其详细参数借助工程方程求解器（engineering equation solver, EES），根据质量守恒方程、动量守恒方程、漂移流模型等编程计算^[25-26].

流程内核心部件精馏塔采用RadFrac模型，以理论板为基本计算单位，并做出以下假设：每一级塔板处于相平衡状态；气、液相以理想状态混合；塔板压降为1.3 kPa；忽略雾沫和液沫夹带. 在此基础上建立平衡级MESH方程^[24].

物料平衡方程（M方程）、相平衡方程（E方程）、归一化方程（S方程）、能量衡算方程（H方程）分别如下：

(1)$ \begin{split} {M}_{ij}=&{L}_{j-1}{x}_{i,j-1}+{V}_{j+1}{y}_{i,j+1}-\left({L}_{j}+{S}_{j}\right){x}_{i,j}-\\ &\left({V}_{j}+{s}_{j}\right){y}_{i,j}+\left({F}_{j}{Z}_{ij}^{\text{L}}+{f}_{j}{Z}_{ij}^{\text{V}}\right)\text=0,\end{split} $

(2)$ {y}_{i,j}\text={K}_{i,j}{x}_{i,j}, $

(3)$ \sum\limits_i {{x_{i,j}} - 1 = 0} ,\quad\sum\limits_i {{y_{i,j}} - 1 = 0},$

(4)$ \begin{split} {M_j}{H_j}{\text{ = }}&{L_{j - 1}}{H_{j - 1}} - \left( {{L_j}+{S_j}} \right){H_j}+{F_j}H_j^{\text{F}} -{E_j}+ \\ & {V_{j+1}}{h_{j+1}} - \left( {{V_j}+{s_j}} \right){h_j}+{f_j}h_j^{\text{f}} - {e_j}{\text{ = }}0. \end{split} $

式中：下标i为单一组分，j为塔内第j个非平衡级；上标字母V表示气相，L表示液相，f表示气相进料，F表示液相进料；h、H分别为气、液相焓值；V_j、L_j分别为气、液相流率；y_{i, j}、x_{i, j}分别为气、液相组分$ i $所占的摩尔分数；f_j、F_j分别为气、液相进料流率；s_j、S_j分别为气、液相侧线产率；e_j、E_j分别为气、液相热负荷；Z_ij表示进料组成中i组分的摩尔分数.

HYSYS中关于塔收敛求解策略采用Inside-out算法，通过该算法对MESH方程求解，逐板计算给定塔的操作结果. Inside-out算法表达式如下：

(5)$ \phi = {\sum\limits_k {\omega k\left( {\frac{{\widehat Gk - Gk}}{{{G^0}k}}} \right)} ^2}. $

式中：k为设计规定序号，$ \widehat G $为设计值，G为期望值，G⁰为比例，$ \omega $为权重因子.

其余压缩机、换热器一律满足质量、能量守恒方程. 在模拟计算中，各换热通道内压降取5 kPa，精馏柱内压降依据狄克松散堆填料的经验参数取9 kPa，塔顶氮气纯度要求大于等于99.5%.

通过控制变量，充分考虑空气入口温度θ₀为10~25 ℃、压缩空气初步预冷温度θ₂为5~35 ℃、制冷机输入冷量Q为0.8~1.0 kW（系统入口空气质量流量随之改变）的不同运行条件，对制氮系统的液氮产率$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $、氮提取率N以及单位产量能耗W的影响. 其中氮提取率和单位产量能耗的表达式如下：

(6)$ N = {{\left( {q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} \times 0.993\;9} \right)} \mathord{\left/ {\vphantom {{\left( {q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} \times 0.993\;9} \right)} {\left( {q_m^{{\text{Air}}} \times 0.755\;5} \right)}}} \right. } {\left( {q_m^{{\text{Air}}} \times 0.755\;5} \right)}}. $

式中：$ q_m^{{\text{Air}}} $为系统入口空气质量流量.

(7)$ W = \left( {{w_1} + {w_2} + {w_3}} \right)/q_m^{{\rm{L}}{{\rm{N}}_{\rm{2}}}}. $

式中：w₁为压缩设备功耗，w₂为初步预冷设备功耗，w₃为制冷机电功耗.

在优化计算流程中，制冷机作为系统唯一的电-冷来源，其输入冷量Q被3处流量所消耗：一是精馏柱塔顶纯化的氮气被液化所需的冷量；二是气泡泵被加热汽化部分流体重新被液化所需的冷量，该量值由EES气泡泵计算程序的相关参数与HYSYS流程中相关数据不断迭代直至两者统一而获得，具体数值为8.8~10.1 W；三是CO₂析出脱离系统所带走的冷量（仅针对流程1）. 此外，默认制冷机工作在额定状态，w₃取11.0 kW定值. 空气入口温度θ₀仅影响w₁和w₂，因此以下对变量θ₂和Q的模拟分析在θ₀=10 ℃的情况下展开.

通常对精馏柱模块进行初步模拟，以“塔板数（P）×回流比（R）”这一指标对精馏设备进行经济衡算，当该指标最小时，认为精馏设备的综合费用最低.

借助Aspen HYSYS软件，建立如图3所示的简易精馏柱模型，设置入口流体M1的质量流量，使精馏柱冷凝器负荷满足制冷机可提供制冷量1 kW，以馏出物产品M2中氮纯度达99.5%为平衡标准，对入口流体M1的温度、压力以及精馏柱塔板数进行变量分析，求解不同工况下的最佳理论塔板数.

如图4所示为常压精馏与加压精馏在不同进料工况下的经济衡算结果. 图中，x为干度，$\Delta \theta $为过热度. 可以看出，在同一进料工况下，随着理论塔板数的增加，精馏柱的经济衡算结果均呈现先减小后增大的趋势；在同一进料压力条件下，干度x或过热度$\Delta \theta $越大，经济衡算结果越大；在同一进料温度条件下，进料压力越大，最佳塔板数越大. 在常压精馏状态下，最佳理论塔板数为7；在加压精馏状态下，最佳理论塔板数为11.

如图5所示为流程1、2在空气入口温度θ₀为10 ℃时，在不同初步预冷温度θ₂以及制冷机输入冷量Q的工况下，系统液氮产率$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $的变化情况. 以θ₂为10 ℃为例（Q-$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $平面投影曲线），Q越大，$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $相应越大. 以Q为1 kW为例（θ₂-$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $平面投影曲线），在模拟范围内，流程1的$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $随θ₂的减小，先增加后保持不变；流程2的$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $随θ₂的减小持续增加. 此外，如图5(a)所示，在θ₂＞25 ℃时，预冷换热器内冷流体冷量已被充分回收，但精馏柱进料温度仍大于饱和温度，属于过热气进料状态. 根据能量守恒定律，塔釜进料温度越高，等冷负荷下被液化的量就越少. 当θ₂＜25 ℃时，预冷换热器内冷流体的冷量被充分回收后仍有富余，此时精馏柱进料状态恒定，系统伴随有能量浪费，在该范围内$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $达到最大值，当Q=1 kW时，$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $可达8.349 kg/h. 如图5(b)所示，在流程2的结果中，由于模拟范围内预冷换热器处冷量都被回收完全，在θ₂=10 ℃，Q=1 kW时，$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $最大可达9.008 kg/h；在θ₂=25 ℃，Q=1 kW时，$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $=8.548 kg/h.

如图6所示为流程3液氮产率$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $和氮提取率N的变化曲线. 由于该流程在θ₂模拟范围内，预冷换热器处冷量均被充分回收且有冷量富余，这使得精馏柱进料状态在模拟范围内一直保持不变，此时$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $与θ₂失去关联，$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $随Q的增大而增大，在Q=1 kW时$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $=6.734 kg/h.

综上，从液氮产率的角度分析，3种流程在θ₂和Q相同的工况下，流程2由于等焓节流从压力势能转换得到额外冷量，液氮产率相对最大.

如图7所示为流程1、2在空气入口温度θ₀为10 ℃时，在不同初步预冷温度θ₂以及制冷机输入冷量Q的工况下，系统氮提取率N的变化情况.

如图7(a)所示为流程1的结果，以θ₂=35 ℃为例（Q-N平面投影曲线），Q越大，而系统入口流量$q_m^{{\mathrm{Air}}} $也增大，但与液氮产率$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $的增幅不一致，系统精馏柱进料温度存在小幅度上升，精馏回流比增大，N存在小幅度提升的现象；以Q=1 kW为例（θ₂-N平面投影曲线），在模拟范围内，N随θ₂的减小，先降低后保持不变，其转折点位于25 ℃处，该处系统氮提取率N=62.19%. 当θ₂＞25 ℃时，系统处于过热气进料状态，在精馏柱内回流的液氮被大量汽化返回顶部，塔釜废液减少，进而θ₂越大，N越高；当θ₂＜25 ℃时，进料状态恒定，回流比基本无波动，因此N也基本不变.

如图7(b)所示为流程2的结果，以θ₂=25 ℃为例（Q-N平面投影曲线），Q的改变对N基本不产生影响，这是因为节流措施的存在，原先应有的温度波动变得更小，精馏柱进料温度基本恒定；以Q=1 kW为例（θ₂-N平面投影曲线），在模拟范围内，预冷换热器中冷量回收充分，θ₂越低，精馏柱进料温度越低，所需回流比越小，导致N越小. 在θ₂=25 ℃、Q=1 kW的工况下，N=58.97%.

流程3的氮提取率结果见图6，在模拟范围内，N仅随Q的增加小幅度上升，因为系统操作压力大，模拟范围内精馏柱均以气液两相状态进料，无再沸器的精馏结构下存在液相浪费，N相较其余两流程会低很多. 在θ₂=25 ℃、Q=1 kW的工况下，N=38.54%.

综上，从氮提取率的角度分析，对比文献调研，3种流程的氮提取率均远大于常温提纯的PSA手段，且在同工况下，流程1的氮提取率最高.

如图8所示为3种流程在θ₀=10 ℃时，在不同θ₂以及Q的工况下，系统单位产量能耗W的变化情况. 3种流程在θ₂一致的情况下，系统单位产量能耗W均随Q的增加而减小；针对变量θ₂，在模拟范围内可以看到，在流程1、2中，存在W随θ₂的减小呈现先减后增的趋势；在流程3中，在模拟范围内W随θ₂的减小不断增加，这是由于流程3的θ₂的改变只影响了初步预冷设备功耗w₂. 由图可得，在Q=1 kW的工况下，当θ₂=25 ℃时，3种流程均处于最节能状态，单位产量能耗W分别为 4946、6067、8076 kJ/kg（比功耗分别为1.105、1.356、1.534 kW·h/L）.

最后，对不同空气入口温度θ₀进行模拟分析，鉴于所有流程涉及到的精馏环节直接相关的温度都是初步预冷后温度θ₂，θ₀只影响了压缩设备功耗w₁和初步预冷设备功耗w₂. 在θ₂=25 ℃时，流程1、2的w₁和w₂随θ₀和Q的变化曲线如图9所示，可见w₁、w₂与θ₀、Q均呈单调关系，θ₀、Q越小，两设备能耗越小. 并且对比图9(a)、(b)，在相同工况下，流程2设备能耗约为流程1的6~7倍.

综上，从单位产量能耗的角度分析，流程1的系统能耗大幅度低于其余2种流程的，并且结合液氮产率，其单位产量能耗在θ₀=10 ℃，θ₂=25 ℃，Q=1 kW的工况下存在最小值，此时系统运行最为高效节能.

如图10所示给出了系统最佳运行温度θ₀=10 ℃、θ₂=25 ℃下的流程曲线，3种流程的液氮产率$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $、单位产量能耗W均与制冷机输入冷量Q呈单调相关且趋势一致，Q越大，$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $越大、W越小. 综合对比发现，流程1仅在液氮产率上略逊一筹，在能耗层面远低于其余2种. 并且对于流程3而言，加压精馏使得精馏柱塔顶饱和氮气的温度远大于77 K的冷头温度，制冷机冷头换热不可逆损失大的问题没有得到改善.

综上，在3种低温精馏制氮系统中，选用压比为1.614的低压风机进行空气压缩，采用冷冻干燥法纯化空气，最后进行常压精馏的方案是综合最优的. 该系统方案在θ₀=10 ℃、θ₂=25 ℃、Q=1 kW的运行条件下，可以实现$ q_m^{{\text{L}}{{\text{N}}_{\text{2}}}} $=8.349 kg/h，N=62.19%，W=4946 kJ/kg（比功耗为1.105 kW·h/L）的最佳性能.

此外，将设计流程与现有小型液氮生产系统性能相比较，如表1所示归纳总结了性能对比结果. 表中，$q_V^{{\mathrm{LN}}_2} $为LN₂产率 (体积流量)，SEC为比功耗，P为纯度. 对比发现，本研究所涉及的液氮生产系统在液氮产率以及比功耗方面都占有明显优势. 表中，带“*”的数据通常只考虑了液化过程所需设备，忽略了空压机、冷水机组这些设备的能耗；带“**”的数据忽略了空压机、冷水机组的电转换效率以及吸附纯化设备的能耗. 为了方便对比，表中对LN₂产率数值均进行了体积单位换算.

本研究提出3种采用大冷量斯特林制冷机的低温精馏制氮系统. 由流程模拟结果可知，采用鼓风机压缩空气，在预冷换热器内以冷冻干燥的方法，在回收系统塔釜冷流体冷量的同时纯化空气，而后进行常压精馏，是最优的系统方案.

模型综合考虑了空气入口温度、压缩空气初步预冷温度、制冷机输入冷量3个运行条件对各流程的性能影响，其中压缩空气初步预冷温度对流程的性能影响最大，存在最优解使得系统单位产量能耗达到最小值. 结果表明，“低压进气-常压精馏”是最优系统方案，在空气入口温度为10 ℃，初步预冷温度为25 ℃，制冷机输入冷量为1 kW（该状态下系统入口空气质量流量为17.66 kg/h）的运行条件下，系统液氮产率（体积流量）为 10.38 L/h，氮提取率高达62.19%，系统比功耗为 1.105 kW·h/L. 与现有小型液氮生产系统相比，在产品纯度需求较高的情况下，系统氮提取率是常温提纯PSA技术的3倍多，在液氮产率、比功耗方面具有明显优势；并且，由于精馏技术在塔顶获得的是低温饱和氮气，在制冷机冷头处换热温差小，换热不可逆损失大幅度降低. 除此之外，该流程在实际运用上，随着仪器仪表的减少，可以降低进气压比，改用翅片管自然对流散热的形式进行初步预冷至25 ℃，进一步实现系统节能.

当前研究仅进行了理论验证，为实验提供了前期指导，后续须进一步进行实验佐证，完善具体流程设计. 其中，对于进行空气纯化的板式换热器，如何将CO₂在换热通道内高效析出并去除是后续设计难点之一；低压大流量进气风机、系统低温绝热技术、系统启动初期空气纯化用冷量来源等均须以高效节能为目的慎重选取设计.