Traditional cryogenic air separation units are unsuitable for distributed, small-scale liquid oxygen production. Cryocooler-based liquefaction technology offers an alternative solution, featuring a large cooling capacity, high efficiency, a compact structure, and rapid start–stop capability. In this paper, an oxygen liquefaction system based on a high-capacity Stirling cryocooler was developed. Because the heat transfer performance of cryocoolers varies significantly across different temperature ranges, heat exchanger designs must be tailored to specific operating conditions. However, research on cold-end heat exchangers for large-capacity cryocoolers used in liquefaction systems remains limited. In the liquid oxygen temperature range, factors such as liquid film formation and incomplete condensation severely affect heat transfer performance and must be considered. In this paper, numerical simulations were performed to analyze the condensation behavior of oxygen, with particular attention paid to the matching between the heat exchange structure and the cooling capacity. Subsequently, a small-scale experimental system was constructed and tested. The successful operation of the experimental system validated the feasibility of the proposed heat exchanger design. Under the conditions of 300 K and an oxygen inlet gauge pressure of 0.45 MPa, the system achieved a liquefaction capacity of 7.4 L/h, corresponding to a cooling capacity of 787 W. The specific power consumption was 0.89 kW·h/kg, with a coefficient of performance (COP) of 0.116. This performance is competitive among small-scale cryocooler-based oxygen liquefaction systems. This study provides both theoretical and experimental support for further performance optimization and engineering application of such cryocoolers in liquid oxygen production. 1. Introduction Currently, liquid oxygen is produced on a large scale primarily through cryogenic air separation [ 5]. This technology exploits differences in the boiling points of air components. Through cryogenic distillation, impurities such as nitrogen and argon are removed, yielding oxygen with a purity of at least 99.99% [ 6]. However, such systems are complex, costly, and inflexible, with slow startup profiles. Thus, they are economically viable only for large-scale, continuous industrial operations. Crucially, technical constraints currently prevent effective miniaturization of the equipment. Consequently, this method is not feasible for small-scale liquid oxygen production. To meet the demand for flexibility in small-scale applications, the hybrid approach combining non-cryogenic separation with liquefaction has become a mainstream solution. This method typically first produces oxygen-enriched gas via techniques such as pressure swing adsorption (PSA) [ 7] and membrane separation [ 8]. The gas is then condensed into liquid oxygen using a compact liquefaction system. This integrated system offers modularity and rapid start-up capability, making it suitable for small-scale scenarios. Stirling Cryogenics B.V. [ 19] has developed a medical oxygen liquefaction system based on a Stirling cryocooler, which achieves a liquefaction rate of 5–36 Nm 3/h for medical-grade oxygen (93% purity at 0.4 MPa). Essex Industries Inc. [ 20] created a portable unit for oxygen production and liquefaction with a capacity of 1 L/h. RIX Industries Inc. [ 21] designed a shipboard liquid oxygen system utilizing a pulse tube cryocooler, capable of liquefying 55 gallons per day (GPD). This system provides both liquid and gaseous oxygen, featuring high integration, autonomous operation, and low lifecycle costs. Boissin et al. [ 22] proposed a portable medical device for delivering oxygen-enriched gas, compatible with Stirling cryocoolers. Spoor et al. [ 23] deployed an oxygen liquefaction system based on a pulse tube cryocooler aboard a U.S. Navy aircraft carrier, achieving a liquefaction rate exceeding 60 GPD. To address the need for small-scale oxygen liquefaction, this study developed an oxygen liquefaction system based on a large-capacity Stirling cryocooler. In cryocooler-based liquefaction systems, the cold-end heat exchanger is a critical component responsible for transferring cooling power to an external fluid. Its performance determines whether the refrigeration capacity meets the requirements and enables efficient liquefaction. Research on large-capacity cryocoolers for liquefaction has primarily concentrated on enhancing the performance of the cryocoolers themselves. Studies focusing on the performance of the cold-end heat exchanger during the liquefaction process remain scarce. The nominal cooling capacity of a cryocooler is typically measured using the thermal opposition method [ 24]. In this method, heating resistors are placed inside the cold-end heat exchanger, and the applied heating power counterbalances the cooling effect at a steady target temperature. The cooling capacity equals the applied electrical heating power minus parasitic heat losses. In practice, optimization should be performed based on the actual application scenario. For small-scale cryocoolers, the cold-end design is less critical because their cooling capacity is typically low. However, for liquefaction or cooling processes involving kilowatt-class, high-capacity cryocoolers, the cold-end heat exchanger typically adopts structures such as fins or coils to effectively increase the heat transfer area and provide appropriate flow paths. Under actual operating conditions, the actual heat transfer performance delivered by the cold-end heat exchanger and the effective cooling capacity it provides will differ from the cryocooler’s nominal cooling capacity measured using the thermal opposition method. This study employs computational fluid dynamic simulation methods to investigate the condensation process of superheated oxygen on the outer side of the cold-end heat exchanger. Based on the simulation results, the oxygen condensation characteristics inside the heat exchanger are revealed, and its liquefaction performance is evaluated. Accordingly, the system feasibility is validated, and a gas reliquefaction scheme is proposed. Further experimental research validated the actual performance of the liquefaction system, with results demonstrating that the system exhibits excellent liquefaction efficiency. The developed oxygen liquefaction system based on the Stirling cryocooler possesses high performance and high reliability. 2.1. The Simulation Model The type and design of the cold-end heat exchanger play a critical role in the performance of a cryocooler, particularly during gas liquefaction. In the slit-type cold-end heat exchanger, heat transfer proceeds in three stages from the inside outward. First, heat is transferred between the helium and the inner slit. Second, heat conducts through the fin. Third, heat is exchanged via convection and condensation between the oxygen and the outer slit. For small-capacity cryocoolers, the cold-end heat exchanger generally requires no customized design and only needs to satisfy the heat transfer requirements of the internal helium circulation. In contrast, for large-capacity cryocoolers, the design must account for the heat transfer capability on the helium side and provide sufficient space for the external liquefaction process. illustrates a single liquefaction channel in the external region of the Stirling cryocooler’s cold-end heat exchanger. The channel features a rectangular cross-section with a width of 1.6 mm, length of 20 mm, and height of 90 mm. The cold-end heat exchanger has a total of 80 slits. A half-symmetry model is adopted to reduce computational cost. One side of the channel is defined as the cold wall, the opposite side as an adiabatic boundary. The model assumes that the inlet flow direction aligns with gravity. This assumption is widely adopted in condensation analysis for vertical narrow channels, where the main flow is dominated by gravity. In practice, however, the flow direction inevitably deviates, which hinders liquid film flow and phase change, thereby reducing the heat transfer rate. The wall where the fin is located serves as a thermally conductive wall and is made of copper. Heat conduction through its thickness is considered in the calculation, and the fin thickness is set to 1 mm. Since multiple rectangular channels are arranged in parallel, each fin is shared by the channels on both sides. Therefore, in the modeling, half of the fin thickness is taken based on symmetry. The oxygen vapor is treated as an ideal gas, allowing its temperature-dependent density variation to be modeled. Liquid oxygen is assumed to be incompressible with a constant density. The inlet was defined as a mass flow boundary with superheated oxygen at 300 K. The computational mesh is refined near the inlet, cold wall, and adiabatic boundary, with a minimum cell size of 1 µm × 1 µm × 1 µm. A non-uniform grid with a bias growth rate of 1.05 is employed to achieve a smooth transition. The mesh spacing gradually increases in the direction of gravity and with increasing distance from the wall. A preliminary grid-independence test was conducted on mesh configurations with cell counts ranging from 0.86 million to 2.44 million. Considering the trade-off between liquid film resolution and computational accuracy of heat transfer, the mesh with 1.48 million cells was selected for all subsequent calculations. With this mesh, the calculated cooling capacity differs by less than 2% from that obtained with the 2.44 million cell mesh. Oxygen cools from room temperature (300 K) to its saturation point. Over this range, its thermophysical properties vary significantly with temperature. To ensure numerical accuracy, the data from NIST are used. These data are fitted and expressed as polynomial functions of temperature to represent the thermophysical properties of gaseous oxygen, as shown in Equations (1)–(4). In addition, the thermal conductivity of copper is also temperature-dependent, and its fitted correlation is provided in Equation (5). The fitting error is indicated by R 2 > 0.99 for all correlations. At a gauge pressure of 0.45 MPa, the latent heat of vaporization at this condition is 189,443 J/kg. C p ν T = − 9.96 ୍ଠ 10 − 9 T 5 + 1.11 ୍ଠ 10 − 5 T 4 − 0.005 T 3 + T 2 − 114 T + 5926 (1) ρ ν T = − 2.68 ୍ଠ 10 − 6 T 3 + 0.002 T 2 − 0.55 T + 61 . 17 (2) k ν T = 8 . 42432 ୍ଠ 10 − 5 T + 0.0016 (3) μ ν T = 6.39349 ୍ଠ 10 − 7 T + 1.85858 ୍ଠ 10 − 5 (4) k s T = 39596.11 ୍ଠ e x p − 0.07382 T + 414.97 (5) where Cpv is the specific heat capacity at constant pressure of oxygen vapor (J·kg −1·K −1), ρv is the density of oxygen vapor (kg/m 3), kv is the thermal conductivity of oxygen (W·m −1·K −1), ks is the thermal conductivity of copper (W·m −1·K −1), and μv is the dynamic viscosity of oxygen vapor (Pa·s). Condensation occurs mainly at the gas–liquid interface, and the Volume of Fluid (VOF) method is adopted in this paper for numerical simulation. By defining the phase volume fraction within each computational cell, the VOF method enables explicit tracking of the position and morphology of the gas–liquid interface, with the sum of volume fractions of the two phases equal to unity in every control volume. The flow is laminar at a low Reynolds number. A steady-state solver is used in the numerical calculation, and the transient terms in the governing equations are eliminated. The governing equations are given as follows: The volume fraction equation is ∇ ⋅ α q u = S ρ q (6) where αq is the volume fraction of phase q (dimensionless); αvαl are the volume fractions of vapor and liquid phases; ρq is the density of phase q (kg/m 3); u is the velocity vector (m/s); and S is mass source due to phase change (kg/m 3·s). The continuity equation is ∇ ⋅ ρ u = 0 (7) The momentum equation is ∇ ⋅ ( ρ u u ) = − ∇ p + ∇ ⋅ [ μ ( ∇ u + ( ∇ u ) T ) ] + ρ g + F σ (8) where p is the pressure (Pa); μ is the molecular dynamic viscosity (Pa·s); g is the gravitational acceleration vector (m/s 2); Fσ is the surface tension force per unit volume (N/m 3). Cryogenic fluids have low surface tension, which exerts negligible effects on gravity-driven flow, so the surface tension term is omitted. The energy equation is ∇ ⋅ ρ u h = ∇ ⋅ k ∇ T + q ˙ (9) where h is the specific sensible enthalpy (J/kg), k is the effective thermal conductivity (W/(m·K)), T is the temperature (K), and q ˙ is the energy source due to phase change. To model the phase change process at the interface, this study adopts the Lee model [ 25]. The Lee model is widely used for gas–liquid phase change simulation and offers good convergence in numerical calculations. Based on the Hertz–Knudsen equation and the Clausius–Clapeyron equation [ 26], the evaporation–condensation mass flux m ˙ i on a smooth gas–liquid interface is expressed as: m ˙ i = 2 κ 2 − κ h f g 1 2 π R T s a t ρ ν ρ l ρ l − ρ ν T s a t − T T s a t (10) where m ˙ i is the net mass flux per unit area across the interface (kg/m 2·s), hfg is the latent heat (J/kg), R is the specific gas constant (J/kg·K), Tsat is the saturation temperature (K), and κ is the condensation coefficient. Its physical definition is: κ = Number of gas molecules absorbed by the liquid phase Number of gas molecules striking the liquid phase , 0 T s a t (15) The corresponding energy source term is: q ˙ = m ˙ l h f g = − m ˙ ν h f g (16) The condensation frequency coefficient r in Equations (13)–(15) is difficult to determine quantitatively due to dependencies on condensing conditions, surface geometry, and working fluid. The choice of condensation frequency involves a trade-off between physical fidelity and numerical stability. A frequency that is too low fails to maintain thermal equilibrium at the interface, causing the interface temperature to deviate unrealistically from the saturation state. Conversely, if the frequency is too high, convergence difficulties may arise, leading to solution oscillations or even divergence. Liu et al. [ 27] conducted a numerical simulation study on the filmwise condensation of water vapor between two parallel plates and determined an appropriate condensation frequency of 5000 s −1. In numerical simulations of R134a condensation in microchannels, Riva et al. [ 28] set the condensation– frequency coefficient r to 7.5 × 10 5 s −1. For condensation simulations of R410a in round and flat tubes, Zhang et al. [ 29] adopted r = ୧.୫ ୍ଠ ୧୦ 6 s −1, which limited the interfacial temperature deviation to within 0.5 K. In experimental studies on R32 condensation, Wu et al. [ 30] used r= ୩.୩୩ ୍ଠ ୧୦ 9 s −1. Shen et al. [ 31] proposed a correlation linking the condensation frequency to the interfacial temperature. To assess the applicability of the numerical model to cryogenic fluids, the condensation process of oxygen on a vertical flat plate was simulated under an absolute pressure of 0.3 MPa and a subcooling of 1 K. A series of condensation coefficients with varying magnitudes were selected for comparative analysis. Due to the limited availability of experimental data on oxygen condensation under cryogenic conditions, the comparison focuses on the agreement between the simulation results and the Nusselt theoretical solution, as shown in . To quantitatively evaluate the accuracy of the CFD model, the condensation simulation results were compared with the Nusselt theoretical solution, and a sensitivity analysis of the condensation coefficient r was conducted. When r = ୧ ୍ଠ ୧୦ 6, the simulation results are in good agreement with the theoretical solution. The CFD-predicted outlet liquid film thickness is approximately 42.28 μm, with a relative deviation of about 1.77% from the Nusselt theoretical value (43.04 μm). This deviation is within an acceptable range, indicating that the current CFD model can reasonably exhibit the theoretical trend of film condensation under simplified conditions. When r = ୧ ୍ଠ ୧୦ 5, the CFD-predicted outlet liquid film thickness is approximately 35.27 μm, with a relatively large deviation (about 18.05%); when r = ୧ ୍ଠ ୧୦ 7, the solution fails to converge. Based on the above results, r = ୧ ୍ଠ ୧୦ 6 was adopted for all subsequent simulations. It should be noted that the model calibration is based on ideal laminar condensation theory and does not account for vapor shear-induced interfacial waves or local turbulent effects. Consequently, the model may overestimate the liquid film thermal resistance and underestimate the actual heat transfer capability, which inevitably leads to certain discrepancies between the simulation results and the experimental data. In this paper, the PISO algorithm is used to handle pressure-velocity coupling. This algorithm is based on a high-order approximate relationship between pressure and velocity corrections, which effectively enhances convergence efficiency. For pressure discretization, the PRESTO scheme is selected, as it is suitable for flows involving strong swirl or significant density variations. The advection terms in the momentum and energy equations are discretized using the QUICK scheme. On structured grids, this scheme provides third-order accuracy, reducing numerical diffusion and improving solution stability. The volume fraction is discretized with the Compressive scheme to ensure accurate phase interface capture. The residual of the energy equation was set to 10 −6, while that of the momentum and volume fraction equations were set to 10 −4, and the continuity equation residual was set to 10 −3. 2.2. Simulation Results and Analysis With an inlet temperature of 300 K and an outlet saturation temperature of 110.18 K, when the cryocooler’s cooling capacity is 1000 W, the enthalpy difference is 369.941 kJ/kg, and the oxygen flow rate is approximately 2.7 × 10 −3 kg/s. The number of channels machined on the side of the cold-end heat exchanger is taken as 80. In the symmetric model, the inlet corresponds to half of the actual channel cross-sectional area, so the flow rate is approximately 1.69 × 10 −5 kg/s. Therefore, the inlet oxygen mass flow rate in the model was set to 1 × 10 −5 kg/s, 2 × 10 −5 kg/s, and 3 × 10 −5 kg/s, corresponding to conditions M1, M2, and M3, respectively. The differences between the cold wall temperature and the saturation temperature of the working fluid at a given pressure were set to 1 K, 2 K, and 3 K, and are labeled as T1, T2, and T3. illustrates the variation in the total heat transfer rate of the cold-end heat exchanger under different operating conditions. According to condensation heat transfer theory, the phase change mass flow rate is proportional to the degree of subcooling. Therefore, a larger temperature difference results in more condensate per unit time and greater latent heat release, leading to a higher total heat transfer rate. Under the same temperature difference, increasing the mass flow rate enhances the total heat transfer rate. A higher mass flow rate strengthens convective heat transfer and reduces the thermal resistance on the vapor side. It also generates greater interfacial shear stress, which thins the liquid film and lowers its conductive resistance. A lower thermal resistance thus results in a higher heat flux under the same temperature difference. The increase in total heat transfer rate due to a higher mass flow rate is consistent across all temperature differences examined, indicating that the effect of mass flow rate on reducing thermal resistance is independent of the temperature difference. Meanwhile, an increase in pressure leads to a corresponding increase in the cooling capacity. At mass flow rates of 1 × 10 −5 kg/s to 3 × 10 −5 kg/s, a temperature difference of 2 K, and a gauge pressure of 0.45 MPa, the total heat transfer rate of the 160 half-channels ranges from 1002.8 W to 1348.5 W. The cryocooler from Zhejiang Chillmaking Cryogenic Technology Co., Ltd. (Hangzhou, Zhejiang, China) used in the system has a cooling capacity of 1200 W at 110.18 K based on the thermal opposition test, with a helium operating pressure of 2.12 MPa. The total heat transfer rate range calculated from the condensation heat transfer analysis is close to the nominal cooling capacity of the cryocooler under the operating condition of 2.12 MPa, indicating that the heat exchanger is well-matched to the cryocooler. shows the average heat transfer coefficient under various working conditions. With a constant temperature difference, a higher mass flow rate produces a larger average heat transfer coefficient. At a gauge pressure of 0.45 MPa and temperature difference of 2.0 K, the values for M3, M2 and M1 are 3520.8, 2196.5 and 1918.6 W/(m 2·K) respectively. At 0.2 MPa with the same temperature difference, the corresponding values are 2116.9, 1873.2 and 1571.9 W/(m 2·K). For a given mass flow rate, the average heat transfer coefficient decreases with rising temperature difference, and this trend is more prominent at higher flow rates, indicating higher sensitivity to temperature variation. Higher flow rates enhance fluid–wall heat transfer and amplify the effect of temperature difference on boundary layer thermal resistance. In addition, the heat transfer coefficient increases with elevated pressure. shows the liquid film distribution along the gravity-driven flow in a rectangular fin channel. Near the inlet ( z = 3000 μm), the oxygen temperature is high. As cooling brings local subcooling to a threshold, a thin liquid film forms on the cold wall ( x = 20,000 μm) and fin surface ( y = 800 μm). At this stage, the channel is vapor-dominated, and heat transfer relies on single-phase vapor convection with limited efficiency. Downstream, the film gradually thickens, and at z = 80,000 μm, the fin surface is fully covered. Film formation and thickening are governed by subcooling, wall wettability, and interfacial shear stress. Initially, subcooling is high but the film is very thin, with surface tension dominating its morphology. Downstream, film thermal resistance increases, heat transfer efficiency decreases, and subcooling decays along the flow, leading to performance deterioration. Even with constant wall temperature, the thickening film shields the phase-change driving force, equivalent to an increase in interfacial thermal resistance. The condensation efficiency shifts from being limited by single-phase vapor convection to being dominated by liquid film conduction, resulting in a severalfold increase in overall thermal resistance. The cold wall temperature gradient drives condensation, fin geometry restricts lateral film spread, and gravity-driven flow generates shear stress promoting downstream film extension. The film near the cold wall rapidly thickens to 24.5 μm ( z = 80,000 μm), while the film farther from the cold wall remains thin (~9 μm). To prevent ice blockage, the slit width must exceed the maximum film thickness. The present model uses a slit width of 1.6 mm, much larger than the simulated maximum film thickness, which also facilitates uncondensed gas return to the cold end for reliquefaction, ensuring long-term stability. presents the temperature distribution at different z-sections, with the mass flow rate of 3 × 10 −5 kg/s and the temperature difference of 2 K. The isotherm of 102.02 K represents the saturation condition at 0.2 MPa and is used to identify the gas–liquid phase interface. At z = 3000 μm, a liquid film starts to form on the cold wall and the adjacent fin surfaces. In this region of the fin surface, the heat transfer is still governed by convection with the gas phase. Due to the low thermal conductivity of oxygen gas, a pronounced temperature gradient is observed near the channel inlet. Even after heat exchange along the fin, the core flow temperature remains above 200 K at z = 25,000 μm. The liquid film fully covers all fin surface at z = 80,000 μm. The latent heat is released at the film surface, the non-condensed gas transfers heat convectively to the film surface, and the heat is finally conducted to the fin surface. This heat-transfer mechanism results in a marked decrease in the local heat-transfer coefficient in regions covered by the liquid film. shows the temperature distribution of the fin surface for different flow rates. The highest temperature is located at the upper-right corner. This is the inlet region of the heat exchange channel, away from the cold wall. With the higher flow rate, the superheated vapor zone is expanded, resulting in the increased maximum fin temperature. In the presented cases, the maximum fin temperature reaches 107.02 K. This result indicates that the superheated vapor zone maintains a substantial temperature difference, from approximately 300 K down to about 102 K. Furthermore, the temperature distribution at the channel outlet cross-section confirms the presence of superheated gas. At the mass flow rate of 3 × 10 −5 kg/s and a temperature difference of 2 K, the outlet gas flow rate is 1.85 × 10 −5 kg/s, indicating that a portion of the gas remains uncondensed. This incomplete liquefaction must be considered in system design. To address this issue, the present system incorporates a gas recirculation and reliquefaction scheme to fully condense the residual gas. The simulation results provide reliable theoretical guidance for the subsequent experimental design. Based on the condensation evolution law, heat transfer performance, and flow characteristics obtained from the CFD simulations, the rationality of the structural parameters of the cold-end heat exchanger was effectively validated. Furthermore, the incomplete condensation characteristics of oxygen revealed by the simulations offer a design basis for arranging the gas backflow and reliquefaction structure in the experimental system. 3. Experiment and Analysis Based on the numerical simulation results, the experimental system was rationally designed and constructed. The structural configuration of the cold-end heat exchanger and the inlet pressure were implemented in accordance with the numerical analysis conclusions. Key indicators such as liquefaction rate, cooling capacity, and specific power consumption were obtained through experimental testing, thereby effectively verifying the overall feasibility of the system. presents the schematic diagram of an oxygen liquefaction system utilizing a high-capacity Stirling cryocooler. A Stirling cryocooler with a nominal cooling capacity of 1200 W at 110 K is used as the cold source for the liquefaction process. The cryocooler employs helium as the working fluid at an operating pressure of 2.12 MPa. Gaseous oxygen enters the liquefaction chamber of the cryocooler through a control valve, where it undergoes liquefaction. The produced liquid oxygen then flows into a storage tank via an inlet valve and is subsequently discharged through an outlet valve. The tank incorporates a self-pressurization line and associated valves, while the vacuum piping is fitted with an intake valve and a drain valve. In the system, incomplete gas liquefaction and re-evaporation of stored liquid oxygen cause gas accumulation. This raises the tank pressure, which hinders system intake and reduces production. Although periodic venting can regulate the pressure, the direct release of oxygen poses a safety hazard due to its oxidizing properties. To solve this, the system redirects vented gas from the tank back to the liquefaction chamber for reliquefaction, thereby maintaining stable pressure. This recirculation design ensures safe operation and improves overall liquefaction efficiency. Liquid oxygen is stored in a dewar vessel (model HBRF80/100H-00) from Hebei Runfeng Cryogenic Equipment Co., Ltd. (Hengshui, Hebei, China) with an effective volume of 80 L and a specified daily evaporation rate of less than 3%. The flask employs vacuum insulation and the oxygen purity is 99.99%. The liquid level is monitored in real time by a level gauge (model MD-80L/(0–0.7) m) from Suzhou Liannuo Cryogenic Equipment Technology Co., Ltd. (Kunshan, Jiangsu, China), with an expanded uncertainty of ±0.35 mm. The inlet oxygen flow rate is measured by a flow meter (model MEMS-0615A, Guangzhou Weiliang Industrial Control Technology Co., Ltd., Guangzhou, Guangdong, China) with an accuracy of ±2% of full scale. Continuous monitoring and data acquisition of temperature, pressure, and flow rate at key system locations are performed using a temperature transmitter (model WG-WZPKBS-246/TH48C, Anhui Weig Instrument Co., Ltd., Tianchang, Anhui, China) with an accuracy of ±0.5% of full scale, and pressure transmitters (model WGPT3051TG7SAM3d1, Anhui Weig Instrument Co., Ltd., Tianchang, Anhui, China) for both the storage tank pressure and the inlet pressure, also with an accuracy of ±0.5% of full scale. The entire system is vacuum-insulated and placed in an environmental chamber. Based on the full-scale accuracy of each instrument and assuming a uniform error distribution, the standard uncertainties of the measured parameters are as follows: ±3.46 SLM for flow rate, ±0.72 K for temperature, and ±0.003 MPa for pressure. Using the error propagation formula, the expanded uncertainty (with a coverage factor of k = 2) of the system’s cooling capacity is estimated to be approximately 4.24% of the measured value. The relative expanded uncertainty of the liquefaction rate is approximately 0.54%. To ensure safe operation of the liquid oxygen and high-pressure oxygen system, all oxygen-contacting components undergo a rigorous degreasing and purification process, including ultrasonic cleaning in alkaline detergent at 65–75 °C, high-purity rinsing with pH verification, and vacuum or nitrogen drying. Residual oil is verified by ultraviolet inspection and a carbon tetrachloride wipe test. The system also implements particulate control, pressure testing, oxygen concentration monitoring, and forced ventilation, with quantitative assessments of leakage, overpressure, and oxygen enrichment risks. It is equipped with an emergency relief device, a ventilation system providing at least six air changes per hour, and explosion-proof and anti-static designs. For the oxygen recirculation and reliquefaction loop, a complete safety framework is established covering oxygen-compatible materials, cleanliness, contamination control, overpressure protection, venting strategies, ignition risks, oxygen-enriched atmosphere management, and safe closed-loop operation. Oxygen is supplied from a high-pressure gas cylinder through a pressure reducing valve. During this process, the flow rate and temperature are difficult to precisely control, leading to variations in the baseline operating conditions across different test runs. The representative data from multiple experiments were selected for analysis in this paper. Based on the typical experimental data selected above, illustrates the evolution of the inlet oxygen pressure over time. Under the condition of an initial gauge pressure of 0.5 MPa at the oxygen inlet, the pressure variation exhibits two distinct operational stages. During the precooling stage, the pressure fluctuates before decreasing. This stage concludes around 90 min. The pressure fluctuations during the precooling stage are partially attributed to manual or automatic adjustments of the intake pressure during the experiment. Thereafter, the pressure stabilizes near 0.45 MPa, indicating the onset of the steady-state liquefaction phase. This reflects the thermal equilibrium at the cold head, where the condensation heat of the oxygen liquefaction process balances the cooling capacity of the Stirling cryocooler. Pressure evolution directly indicates the phase-change process during oxygen liquefaction. The initial decline results from gas cooling and partial condensation lowering the system pressure, while the subsequent steady level confirms sustained liquefaction under equilibrium conditions. shows the tank level and cryocooler power consumption over time. During the precooling stage, the cryocooler cooled from room temperature to the liquefaction temperature of oxygen within 7 min, with its power consumption increasing from 0 to 6.98 kW. Once the cold-end heat exchanger reaches the liquefaction temperature, it begins to produce a small amount of liquid oxygen. However, no liquid oxygen accumulated in the tank during this phase, since the pipelines and tank were still at an elevated temperature. Rapid accumulation of liquid oxygen begins after approximately 48 min. The steady liquefaction stage begins at 90 min with a constant liquid level rise, and the system reaches thermal equilibrium. The overall liquefaction rate during startup (including precooling) is 5.1 L/h, while the steady-state value is 7.4 L/h. The effective cooling capacity is 787 W and the specific power consumption is 0.89 kW·h/kg. At 110.18 K, the system COP is 0.116. The theoretical Carnot efficiency at this saturation temperature is 58.04%, and the measured system COP equals 11.6%, yielding a relative Carnot efficiency of about 20.0%. The actual cryocooler efficiency is higher when heat leakage and other system losses are excluded. Based on the liquid oxygen production rate of 7.4 L/h, the inlet gas mass flow rate per half-channel of the cold-end heat exchanger is estimated to be approximately 1 × 10 −5 kg/s. Simulations show that at this flow rate and a subcooling degree of 2 K, the theoretical cooling capacity is 1035.1 W, while the experimental value is 787 W, yielding a relative error of approximately 31.5%. From the perspective of model assumptions and computational methods, the deviations between simulation and experiment mainly arise from the following three aspects. First, the ideal gravity flow assumed in simulations cannot be fully realized in real single channels, leading to lower actual heat transfer rates. Second, the liquid outlet of the cryocooler cold-end heat exchanger causes uneven temperature distribution. To reduce computational cost, a single-channel model with uniform temperature approximation was adopted, which overpredicts results. Such simulation errors cannot be accurately quantified currently. Third, the simulation uses a constant cold wall temperature and ignores thermal conduction resistance between helium and the cold wall. At a heat transfer rate of 1000 W, this resistance creates a 2.16 K temperature difference. Thus, the temperature difference between helium and the saturation temperature amounts to 4.16 K. The cooling capacity decreases at a rate of approximately 13.6 W/K, resulting in a total loss of 56.6 W for a temperature drop of 4.16 K. From the perspective of system experimental testing, although thermal insulation measures have been implemented for both the cryocooler and the entire system, heat leakage remains difficult to completely eliminate. Specifically, the dewar vessel contributes approximately 6.0 W of heat leakage, the liquid outlet bellows with rubber-plastic insulation cotton contributes approximately 26.2 W, and the vacuum line contributes approximately 9.4 W. In addition, other components such as instruments and valves introduce additional heat leakage. Meanwhile, instrument measurement uncertainty also contributes to the deviation, corresponding to an absolute expanded uncertainty of approximately 33.4 W for a typical cooling capacity of 787 W. In summary, flow direction, temperature inhomogeneity, heat transfer temperature difference, heat leakage losses, and instrument measurement uncertainties collectively result in the experimentally measured cooling capacity being lower than the theoretical prediction. summarizes the deviation analysis between the predicted cooling capacity and the experimental measurement. shows that not all oxygen condenses immediately at the cold-end heat exchanger. Some uncondensed gas accumulates in the pipes and the tank, and although the system is insulated, external heat ingress still causes slight evaporation of the accumulated liquid oxygen. The resulting vapor returns to the cold end through a dedicated return line and is reliquefied. This “evaporation–return–reliquefaction” cycle plays a key role in maintaining stable pressure during continuous liquid production. compares its key performance indicators with similar small-scale oxygen liquefaction products reported in the literature. As parameters including oxygen purity and specific energy consumption are unavailable in the open literature for the majority of the systems compared, certain table entries remain incomplete. Nevertheless, a performance comparison can still be made based on available data. Compared with existing small-scale oxygen liquefaction systems, the proposed system achieves a liquid oxygen production rate of 7.4 L/h using a single Stirling cryocooler. The StirLOX series system from Stirling Cryogenics Inc. (Son, The Netherlands). Ref. [ 19] has a rated gaseous oxygen output of 5–36 Nm 3/h. If all the produced oxygen is liquefied, the corresponding liquid oxygen output is 6.2–45.1 L/h. This system produces about 6.3–12.5 L/h of liquid oxygen per cold head. The system from Essex Industries Inc. (St. Louis, MO, USA) [ 20] achieves only 1 L/h. The thermoacoustic Stirling cryocooler system from RIX Industries (Benicia, CA, USA) [ 21] achieves a production rate of 55 GPD, with a production rate of 2.9 L/h per cold head. The system reported by Spoor from Chart Industries (Troy, NY, USA) [ 23] requires three cold heads in parallel to achieve a production rate of 60 GPD, corresponding to a production rate of approximately 3.2 L/h per cold head. Based on the above comparisons, the proposed system demonstrates good liquefaction capacity. Therefore, it can be concluded that under the tested conditions, the Stirling-based liquefaction system achieves a favorable balance between efficiency and compactness, showing potential applicability for small-scale, distributed liquid oxygen production. To address the demand for small-scale, fast-startup liquid oxygen production, this study developed an oxygen liquefaction system based on a high-capacity Stirling cryocooler. Numerical simulations were performed to investigate the condensation behavior of superheated oxygen within the rectangular channels of the cold-end heat exchanger, with a focus on heat transfer rate, liquid film distribution, and backflow characteristics of uncondensed oxygen. The results show that the total heat transfer rate increases almost linearly with increasing subcooling degree and mass flow rate, while liquid film thickening degrades local heat transfer performance. These condensation characteristics and heat transfer behaviors provide a basis for the structural design of the cold-end heat exchanger and the gas return design of the system. In terms of experimental validation, under inlet oxygen conditions of 300 K and 0.45 MPa, the system achieves a liquefaction rate of 7.4 L/h, a specific power consumption of 0.89 kW·h/kg, an operating COP of 0.116, and an effective cooling capacity of 787 W. The deviation between the measured value and the numerical simulation is primarily attributable to heat leakage losses, subcooling degree deviations, and instrumental measurement uncertainties. Under the tested conditions, the proposed Stirling-based system achieves competitive energy efficiency and offers advantages in structural compactness and rapid startup response, demonstrating promising potential for small-scale liquid oxygen production applications. Figure 1. Schematic of the cold-end heat exchanger simulation model. Figure 1. Schematic of the cold-end heat exchanger simulation model. Figure 2. Nusselt film condensation of oxygen under different condensation coefficients. Figure 2. Nusselt film condensation of oxygen under different condensation coefficients. Figure 3. Total heat transfer rate under different mass flow rates and subcooling degrees. Figure 3. Total heat transfer rate under different mass flow rates and subcooling degrees. Figure 4. Average heat transfer coefficient under different mass flow rates and subcooling degrees. Figure 4. Average heat transfer coefficient under different mass flow rates and subcooling degrees. Figure 5. Liquid film thickness distribution along the flow direction in the rectangular channel. Figure 5. Liquid film thickness distribution along the flow direction in the rectangular channel. Figure 6. Temperature field distribution at different z-section positions. Figure 6. Temperature field distribution at different z-section positions. Figure 7. Fin surface temperature distribution at different oxygen mass flow rates. Figure 7. Fin surface temperature distribution at different oxygen mass flow rates. Figure 8. Schematic diagram of the small-scale oxygen liquefaction system. Figure 8. Schematic diagram of the small-scale oxygen liquefaction system. Figure 9. Variation of oxygen pressure with operating time. Figure 9. Variation of oxygen pressure with operating time. Figure 10. Liquid oxygen level and cryocooler power consumption versus time. Figure 10. Liquid oxygen level and cryocooler power consumption versus time. Energy balance deviation. Energy balance deviation. Deviation Source Valve (W) Percentage Thermal resistance 56.6 7.2% Flow direction deviation — — Non-uniform temperature distribution — — Dewar vessel (Heat leakage) 6.0 0.8% Outlet bellows (Heat leakage) 26.2 3.3% Vacuum pipe line (Heat leakage) 9.4 1.2% Other components (Heat leakage) — — Instrument measurement uncertainty 33.4 4.3% Quantified deviation 131.6 16.7% Unquantified deviation 116.5 14.8% Total 248.1 31.5% “—” indicates that the contribution of this factor cannot be reliably quantified at present. Performance of small-scale oxygen liquefaction systems. Performance of small-scale oxygen liquefaction systems. Source Cryocooler Purity Number of Cold Heads Specific Power Consumption Liquefaction Rate * Application *—indicates data that have not been publicly reported. All liquefaction rates are converted to L/h for comparison. Based on the original references, some values are estimated from the standard density of liquid oxygen at its normal boiling point.
