Calibration method for low-temperature plate-fin heat exchanger coupled with hydrogen secondary catalytic conversion process

Abstract

A calibration method for a low-temperature plate-fin heat exchanger coupled with a hydrogen ortho- and para-catalytic conversion process comprises the following steps: first, determining the overall structure, operating conditions, and structural parameters of each fluid fin of the plate-fin heat exchanger; then, fitting the fluid physical properties, and fitting the physical properties of ortho- and para-hydrogen, respectively, for the hydrogen fluid undergoing the catalytic conversion process; then, constructing a pressure drop calculation equation group, and introducing a catalyst porous medium pressure drop source term in the fluid channel where the hydrogen ortho- and para-catalytic conversion occurs; constructing a baffle temperature node calculation equation group, and introducing a heat transfer source term caused by the expanded surface area of the catalyst if the hydrogen ortho- and para-catalytic conversion occurs in the fluid channel adjacent to the baffle; constructing a fluid temperature node calculation equation group, and introducing a hydrogen ortho- and para-conversion heat source term in the fluid channel where the hydrogen ortho- and para-catalytic conversion occurs, and introducing a component transport equation for calculating the para-hydrogen concentration; finally, iteratively solving the constructed equation group to obtain the distribution of the temperature field, pressure field, and para-hydrogen concentration field of the entire plate-fin heat exchanger. The present invention promotes the efficient use of hydrogen energy.

Description

Technical Field

[0001] The present invention relates to the technical field of low-temperature plate-fin heat exchangers coupled with hydrogen secondary catalytic conversion, and in particular to a calibration method for a low-temperature plate-fin heat exchanger coupled with a hydrogen secondary catalytic conversion process. Background Art

[0002] Hydrogen is a diatomic molecule, consisting of two hydrogen nuclei rotating about their axes. Based on the relative orientation of the two nuclear spins, hydrogen molecules can be divided into ortho-H2 and para-H2, abbreviated as o-H2 and p-H2. Ordinary hydrogen is a mixture of these two forms, and the balance between ortho-H2 and para-H2 depends solely on temperature. At temperatures above room temperature, it is generally known as normal hydrogen, consisting of 75% ortho-H2 and 25% para-H2. The parahydrogen concentration in equilibrium hydrogen decreases with decreasing temperature. When the temperature drops and hydrogen liquefies, orthohydrogen spontaneously converts to parahydrogen, releasing heat and causing significant vaporization of the stored liquid hydrogen. Even on the first day of storage, the total vaporization can reach over 20% of the total stored volume. Even if liquid hydrogen is stored in a perfectly insulated container, it will still vaporize. Within the first 24 hours, approximately 18% of the liquid hydrogen will evaporate, and after 100 hours, this loss will exceed 40%. This leads to liquid hydrogen loss and increases the pressure within the tank. Therefore, a certain margin must be left in the design of liquid hydrogen storage tanks, which increases the cost of hydrogen production, storage, and use. Under natural conditions, the rate of orthohydrogen-parahydrogen conversion is very slow. Therefore, mature hydrogen liquefaction equipment internationally generally uses multi-stage catalysis to accelerate the conversion of orthohydrogen to parahydrogen near the equilibrium concentration during the cooling process of hydrogen liquefaction, resulting in a liquid hydrogen product with a parahydrogen content of over 95%, thereby reducing the evaporation loss caused by the orthohydrogen conversion.

[0003] The plate-fin heat exchanger is an efficient, compact and lightweight heat exchange device. The fin structure is an ideal catalyst carrier and can be used in the continuous hydrogen conversion process. Continuous conversion is to remove the heat of hydrogen conversion in a timely manner through the heat exchange process of hot and cold fluids in the heat exchanger, so that the hydrogen concentration at the corresponding temperature is always maintained at a near-equilibrium hydrogen concentration during the conversion process. It is recognized as the conversion method with the lowest energy consumption. The advantages of the integrated flow heat exchange of the plate-fin heat exchanger and the hydrogen conversion are as follows: (1) the cooling heat exchange of hydrogen and the conversion reaction are carried out simultaneously, which not only accelerates the conversion rate of hydrogen, but also effectively improves the utilization efficiency of low-temperature cold energy; (2) the catalyst particles filled in the heat exchanger channel significantly enhance the heat exchange of hydrogen fluid on the high-pressure side, which can effectively reduce the heat exchange area and reduce the volume of the heat exchange equipment; (3) the heat exchanger and the conversion reactor are combined into one, eliminating the traditional conversion reactor and its supporting cooling equipment, thereby improving the compactness of the hydrogen liquefaction process system.

[0004] Plate-fin heat exchangers used in industry typically involve more than two fluid streams, potentially involving more than 20 fluids, with hundreds of fin channels. Integrating the secondary catalytic conversion of hydrogen into the heat exchanger creates a complex flow and heat transfer process, making conventional CFD numerical simulations inadequate. Mainstream heat exchanger design software such as HTRI, ASPEN MUSE, and APSEN EDR do not yet integrate the secondary catalytic conversion of hydrogen into their heat exchanger design verification models. Overall, constructing a plate-fin heat exchanger design verification model coupled with the secondary catalytic conversion of hydrogen, developing it into a systematic approach, and developing a general-purpose computational program are of strategic importance for promoting the efficient use of hydrogen energy, ensuring national energy security, and fostering new economic growth drivers. Summary of the Invention

[0005] In order to overcome the shortcomings of the above-mentioned prior art, the purpose of the present invention is to provide a calibration method for a low-temperature plate-fin heat exchanger coupled with a hydrogen paracatalytic conversion process, introduce hydrogen paracatalytic conversion into the heat exchanger calculation, and the calculation results include the temperature field, pressure field and parahydrogen concentration field distribution of the entire plate-fin heat exchanger, thereby promoting the efficient utilization of hydrogen energy.

[0006] In order to achieve the above object, the technical solution adopted by the present invention is:

[0007] A calibration method for a low-temperature plate-fin heat exchanger coupled with a hydrogen secondary catalytic conversion process comprises the following steps:

[0008] 1) Given parameters:

[0009] Select a plate-fin heat exchanger with multiple flows, and determine the overall structural parameters, operating parameters, and fin structural parameters of the plate-fin heat exchanger for each flow;

[0010] The overall structural parameters include the arrangement of the heat exchanger, the length L and width W of the heat exchanger core;

[0011] The operating parameters include the fluid of multiple streams, fluid channel layer, mass flow rate, inlet temperature, and inlet pressure;

[0012] The structural parameters of the fin are: when straight fins are selected, the parameter is fin height h i , wing span s i 、Wing thickness t i ;When selecting perforated fins, the parameter is fin height h i , wing span s i 、Wing thickness t i , porosity φ i When selecting serrated fins, the parameter is fin height h i , wing span s i 、Wing thickness t i , pitch l i ;

[0013] The parameters of the hydrogen secondary reforming catalyst particles filled in the fluid channel are: the porosity ε is 0.35, the porosity is defined as the percentage of the void volume between particles in the packing layer to the total volume, the particle diameter is 0.371 mm, and the thermal conductivity of the catalyst particles is 0.58 W·m -1 ·K -1 ;

[0014] 2) Model simplification and initialization:

[0015] The three-dimensional model was simplified to a two-dimensional plane, assuming that the fluid in the fluid channel is uniformly distributed along the z-direction. Furthermore, only steady-state conditions were considered in the calculations, and the heat exchanger was assumed to be thermally insulated from the surrounding environment. Thermal conductivity was ignored in the fluid calculations. The model construction considered the fluid domain, the solid domains of the fins and baffles, and the solid domain of the catalyst particles. In the channel where the hydrogen ortho- and para-catalytic reactions occurred, the following assumptions were made: hydrogen is a single-phase compressible gas composed of ortho- and para-hydrogen; the catalyst packing in the channel is uniform, and the particles are approximately spherical with a uniform diameter.

[0016] There are n+1 layers of partitions and n layers of fin channels. N nodes are arranged along the flow direction in the partitions, and N+1 nodes are arranged along the flow direction in the fluid channels. The temperature of each partition node T in the partition is w,i,k Define a diaphragm element, and each of the two fluid temperature nodes T in the fluid channel i,k and T i,k+1 Define a fluid unit, the fluid temperature node is also used to store pressure data P i,k and P i,k+1 If the secondary catalytic conversion of hydrogen also occurs in the fluid channel, the secondary hydrogen concentration Pa is stored at the corresponding temperature node. i,k and Pa i,k+1 , a total of N×(n+1) baffle units and n×N fluid units are set in a plate-fin heat exchanger;

[0017] The entire calculation process is iterative. Before the calculation of the first outer iteration round, an initialized temperature field, pressure field, and parahydrogen concentration field need to be given. The fluid temperature node in each layer of the fluid channel is initialized to the inlet temperature given in step 1), and the pressure value stored in the temperature node is initialized to the inlet pressure given in step 1). In the fluid channel where hydrogen positive-paracatalytic conversion occurs, the parahydrogen concentration value is also stored in the temperature node, initialized to the inlet parahydrogen concentration given in step 1).

[0018] 3) Fluid property fitting:

[0019] The required fluid properties include specific heat capacity, thermal conductivity, dynamic viscosity, density, and Prandtl number. The physical property data are extracted from the NIST database. The relationship between the fluid density, temperature, and pressure is fitted on a spline surface. The physical properties of other types of fluids are fitted with temperature on spline curves. For the thermal conductivity of solids, a thermal conductivity-temperature spline curve is fitted. Finally, all the parameters of the fitted spline curves and spline surfaces are loaded into memory.

[0020] In the fin channel where hydrogen ortho-paracatalytic conversion occurs, the physical properties of ortho-hydrogen and para-hydrogen need to be fitted separately. If the physical properties of hydrogen fluid at a certain para-hydrogen concentration need to be calculated, the physical properties of the fitted ortho-hydrogen and para-hydrogen are interpolated and calculated.

[0021] 4) Calculate the pressure field distribution in the fluid channel of the plate-fin heat exchanger:

[0022] First, the pressure drop of each fluid unit in the fluid channel is calculated. Depending on whether the secondary catalytic conversion of hydrogen occurs in the fluid channel, the pressure drop calculation method is divided into two types:

[0023]

[0024] Formula (1) is used to calculate the fluid temperature node T i,k and T i,k+1 The pressure drop in the defined fluid unit, if the fluid channel in this layer is not filled with catalyst particles, that is, no hydrogen secondary catalytic conversion occurs, then the first method is adopted to calculate the pressure drop through the friction factor f in the fluid unit. i,k To calculate, in addition G i is the mass flow rate of the fluid in each layer of fin channel, Δl x is the length of the fluid unit along the flow direction, D i The equivalent diameter of the fin structure in the fluid channel of this layer, ρ i,k+1 and ρ i,k By the fluid temperature node T i,k+1 and T i,k And the pressure value P stored in it i,k+1 and P i,k The density determined, is the density of the fluid unit, f i,k and Determined based on the qualitative temperature and qualitative pressure in the fluid unit:

[0025]

[0026] Arrange straight fins, perforated fins or serrated fins in the fluid channel. The friction factor of serrated fins is calculated using the following method:

[0027]

[0028] Where α, θ, and γ are:

[0029]

[0030] For perforated fins, the friction factor is calculated as follows:

[0031] lnf i,k =28.78906-12.31399lnRe i,k +1.565191(lnRe i,k ) 2 -0.06736098(lnRe i,k ) 3 (5)

[0032] For straight fins, the friction factor is calculated as follows:

[0033]

[0034] Re i,k is the Reynolds number:

[0035]

[0036] The fluid temperature node T i,k and T i,k+1 The dynamic viscosity of the defined fluid unit is determined according to the qualitative temperature and qualitative pressure therein;

[0037] If the fluid channel of this layer is filled with catalyst particles for hydrogen secondary catalytic conversion, the second method in formula (1) is needed to calculate the pressure drop in the fluid unit. is the superficial velocity of the fluid in the fluid unit, that is, the average velocity in the fluid unit calculated without considering the catalyst particles filled in the fluid channel, α is the permeability of the catalyst particles, and C1 is the inertial resistance coefficient of the catalyst particles. The calculation methods of the two parameters are:

[0038]

[0039]

[0040] Where d p is the diameter of the catalyst particle, ε is the porosity;

[0041] The pressure drop of each fluid unit is calculated according to formula (1). Given the inlet pressure of each layer of fluid channel, the pressure value stored in each temperature node is updated point by point along the flow direction; the updated pressure field distribution is compared with the old pressure field distribution, and the pressure field residual is calculated. If the residual does not meet the convergence condition, step 2) is repeated, that is, the pressure field is continuously updated based on the old pressure field until the pressure field residual meets the convergence condition, and step 5) is entered. The pressure field distribution updated in this step is used as the qualitative pressure of each fluid unit in the calculation of subsequent steps;

[0042] 5) Calculation of temperature field distribution of plate-fin heat exchanger baffle:

[0043] When the fluid channels adjacent to the partition are filled with catalyst particles, the heat transferred from the fluid to the partition is divided into three parts: one part is the heat transferred from the fluid through the fin surface to the fin solid domain, and then the heat is transferred from the fin solid domain to the partition solid domain; another part is the heat transferred from the fluid through the catalyst surface to the catalyst solid domain, and then through the catalyst solid domain to the partition solid domain; the third part is the heat transferred directly from the fluid through the partition surface to the partition solid domain. If the fluid channels adjacent to the partition are not filled with catalyst, the heat transferred from the fluid to the partition only includes the first and third parts of heat;

[0044] For the solid domain of fins and catalyst particles, in order to simplify the calculation, the heat conduction along the axial direction (+x or -x) is not considered, so the temperature node T i,k and T i,k+1 The one-dimensional heat conduction equation for the fin and catalyst particle solid domains in the defined fluid unit is:

[0045]

[0046] Formula (10) is used to calculate the temperature T of the fin and catalyst solid area s,i,k Temperature distribution along the y direction, λ s,i,k and λ p,i,k are the thermal conductivity of the fin and catalyst particle solid domain, is the heat transfer source term in the one-dimensional heat conduction equation, which includes the heat transfer between the fin surface and the fluid, and the heat transfer between the catalyst particle surface and the fluid:

[0047]

[0048] In formula (11), s,i,k and λ p,i,k The fluid temperature node T i,k and T i,k+1 The thermal conductivity of the solid domain of the fins and catalyst particles in the defined fluid cell, α pis the specific surface area of the catalyst particles. Assuming that the catalyst particles are uniform in size and are all perfect spheres, the calculation method is:

[0049] α p =6(1-ε) / d p (12)

[0050] h c,i,k is the convective heat transfer coefficient between the fluid and solid domain surfaces in the fluid channel. When the fluid channel is not filled with a catalyst, h c,i,k The calculation method is:

[0051]

[0052] Where Pr i,k is determined by the temperature node T i,k and T i,k+1 The Prandtl number of the defined fluid cell, λ fl,i.k is the thermal conductivity of the fluid unit, j i,k The evaluation index used in this fluid unit to evaluate the heat transfer capacity of the fin surface is the following j-factor correlation formula for serrated fins:

[0053]

[0054] For perforated fins, the correlation is:

[0055] lnj i,k =34.57583-15.92678lnRe i,k +2.137607(lnRe i,k ) 2 -0.09544151×(lnRe i,k ) 3 (15)

[0056] For straight fins, the correlation is:

[0057]

[0058] If the fluid channel is filled with catalyst particles, the convection heat transfer coefficient h on the solid domain surface in the fluid unit is c,i,k The calculation method is:

[0059]

[0060] Where Re p It is the Reynolds number calculated by taking the particle diameter as the characteristic size and the superficial velocity in the fluid unit;

[0061] According to equations (10) and (11), we can derive the temperature node T i,k and Ti,k+1 Temperature distribution of the solid domain along the y direction in the defined fluid cell:

[0062]

[0063] Where sh and ch are hyperbolic cosine and hyperbolic sine respectively, θ i,k ,θ i-i,k and θ i-(i+1),k is the excess temperature of the fin, fin top and root; m e,i,k ,θ i,k ,θ i-i,k and θ i-(i+1),j It is expressed as follows:

[0064]

[0065] When the temperature distribution in the solid domain contained in the fluid unit is obtained, the heat conduction between the solid domain contained in the fluid unit and the partition solid domain is obtained by Fourier's law of heat conduction. For a partition unit, the heat from the adjacent fluid unit also has a fourth part of heat, which is the heat conduction between the partition unit and the adjacent partition units. These four parts of heat in the partition unit satisfy the law of energy conservation. The algebraic equations for calculating the partition temperature node are derived as follows:

[0066]

[0067]

[0068] A i-1,k+1 =A i-1,k

[0069]

[0070] A i,k+1 =A i,k

[0071]

[0072]

[0073]

[0074]

[0075] A w,i,k =A i-1,k +A i-1,k+1 +A i,k +A i,k+1 +A w,i-1,k +A w,i+1,k +A w,i,k-1 +A w,i,k+1 (twenty one)

[0076] In formula (21), t sp is the baffle thickness of the plate-fin heat exchanger, λ w is the thermal conductivity of the partition. In the formula, the subscripts (i, k) are marked outside the brackets. The subscripts of the parameters in the brackets are divided into three cases. If the parameter Δl x ,W,t sp ,ε,α p , then no subscript is required. If it is a parameter s, t, h, the subscript is i. If it is a parameter λ w ,λ s ,h c ,m e , the subscript is (i,k); if the subscript outside the brackets is (i-1,k) or (i,k-1), the parameters in the brackets are marked in the same way. In short, the parameters in the brackets are marked in the same way as those in other formulas in the text; λ w,i,k is the partition temperature node T w,i,k and T w,i,k+1 The thermal conductivity of the interface, if k = N, is directly calculated based on the partition temperature node T w,i,N To determine λ w,i,N ,λ w,i,k-1 is the partition temperature node T w,i,k-1 and T w,i,k The thermal conductivity of the interface, if k = 1, is directly calculated based on the partition temperature node T w,i,1 To determine λ w,i,1 In addition, m in formula (21) e,i,k The calculation of can be divided into two cases according to formula (19);

[0077] For each partition element defined by the partition temperature node, it is required to solve the algebraic equation shown in Equation (20). A total of N×(n+1) equations need to be solved. In order to speed up the calculation of the partition temperature nodes, the line iteration method combined with the Gauss-Seidel method is used. The partition solid domain solution is decomposed into multiple lines parallel to the y-axis. Each line contains only one layer of temperature nodes. The values of each node on the same line are obtained by directly solving the algebraic equation. The values of each node in the block formed by the partition temperature nodes on line AA' are implicitly connected, but the advancement from one line to another along the +x direction is performed iteratively:

[0078]

[0079] Where T P represents the required bulkhead temperature node, T N ,T W ,TS ,T E (north, west, south, east) are the adjacent partition temperature nodes, n0 is understood as the identifier of the outer iteration round, b represents the heat transfer source term that cannot be classified as the influence of the adjacent partition temperature node. Since the scanning direction is along the +x direction, that is, from left to right, T N ,T W ,T S Adopting the updated shelf temperature ( Figure 4 As shown), and T E The temperature value updated in the previous outer iteration round is used. For each partition temperature node on a line, a set of algebraic equations in the form of equation (22) is required to be solved. There are a total of n+1 equations on this line. The Thomas algorithm based on the Gauss elimination method is used to solve the system of equations simultaneously.

[0080] When all the partition temperature nodes are updated from left to right, go to step 6);

[0081] 6) Calculation of temperature distribution and parahydrogen concentration distribution in the fluid channel:

[0082] If the hydrogen secondary catalytic conversion process occurs in the fluid channel numbered i, the temperature node T i,k and T i,k+1 The heat of secondary catalytic conversion of hydrogen in the defined fluid unit is:

[0083]

[0084] Where S o-p,i,k The fluid temperature node T i,k and T i,k+1 The heat source term for the secondary conversion of hydrogen in the defined fluid unit, M is the molar mass of hydrogen molecules, ΔH i,k is the heat of reaction of the secondary catalytic reaction of hydrogen in the fluid unit, which is a function of temperature:

[0085]

[0086] The rate of hydrogen conversion per unit volume Calculated according to the Elovich model:

[0087]

[0088] Where, P c is the critical pressure of hydrogen, is determined by the temperature node T i,k and T i,k+1 The parahydrogen concentration in the defined fluid cell, is the reaction rate constant of the Elovich calculation model; a, b1, c, and d are the fitting coefficients of the experimental data, and their values are 1.0924, 59.7 mol·m -3 ·s -1 , -253.9 mol·m -3 ·s -1 , -11.6 mol·m -3 ·s -1 ; To balance the parahydrogen content in hydrogen, the calculation method is:

[0089]

[0090] Where T c is the critical temperature of hydrogen;

[0091] The calculation of the temperature nodes in the fluid channel is divided into two cases. If the hydrogen secondary catalytic conversion process does not occur in the fluid channel, the fluid in the fluid unit exchanges heat with the surface of the partition solid domain and the surface of the fin solid domain. If the hydrogen secondary catalytic reaction occurs in the fluid channel, the fluid in the fluid unit also exchanges heat with the surface of the catalyst solid particles, and the hydrogen secondary catalytic conversion heat also provides heat for the fluid. These heat sources in the fluid unit satisfy the energy conservation law. Combined with the previous equations (18) and (19), the equation group for calculating the fluid temperature node is derived. If the fluid flows along the +x direction, the derived equation group is in the form of:

[0092]

[0093] The influence coefficient B in the formula i,k , B w,i,k ,B w,i+1,k ,B i,k+1 The calculation method is:

[0094]

[0095]

[0096] B i,k+1 =B i,k +B w,i,k +B w,i+1,k (29)

[0097] In the formula, the subscript outside the brackets is (i,k), then the m inside the brackets is m i , represents the mass flow rate in the fluid channel numbered i, c p Writing c p,i,k , is the temperature node T i,k and T i,k+1The specific heat capacity of the fluid unit defined by the equation is: If hydrogen secondary catalytic conversion occurs in the fluid channel, then m in equation (29) e,i,k The calculation method should be used when the catalyst is filled in the fluid channel of formula (19). In this case, formula (28) is not closed, and the component transport equation of parahydrogen needs to be added:

[0098]

[0099] If the fluid flows along the +x direction and the para-catalytic conversion of hydrogen occurs in the fluid channel, the para-hydrogen concentration value stored in the fluid temperature node is determined point by point along the +x direction according to equations (28) and (30);

[0100] If the fluid channel is along the -x direction, the equations for calculating the temperature node of the fluid element are:

[0101]

[0102] The calculation method of the influence coefficient is:

[0103]

[0104]

[0105] B i,k =B i,k+1 +B w,i,k +B w,i+1,k (32)

[0106] The component transport equation of parahydrogen in the fluid unit is:

[0107]

[0108] By iterative method, equations (31) and (33) determine the temperature field and parahydrogen concentration field distribution in the fluid channel where the flow is in the -x direction;

[0109] 7) Calculation of temperature field residual and parahydrogen concentration field residual:

[0110] Compare the partition temperature field distribution calculated in step 5) with the partition temperature field distribution determined in the previous external iteration round, calculate the partition temperature field residual, compare the fluid temperature field distribution calculated in step 6) with the fluid temperature field distribution determined in the previous external iteration round, calculate the fluid temperature field residual, and take the maximum value of the partition temperature field residual and the fluid temperature field residual as the temperature field residual; compare the parahydrogen concentration field calculated in step 6) with the parahydrogen concentration field distribution determined in the previous external iteration round, and calculate the parahydrogen concentration field residual; if one of the temperature field residual and the parahydrogen concentration field residual does not meet the convergence condition, continue to repeat steps 4) to 7) as a new round of external iteration until both residuals meet the convergence condition, that is, output the calculation results, which include the overall temperature field distribution of the plate-fin heat exchanger, the pressure field distribution and the parahydrogen concentration distribution in the fluid channel; the temperature field includes the partition temperature field and the fluid temperature field.

[0111] In step 5), for the partition temperature node, successive sub-relaxations are introduced as needed:

[0112] T n0 =T n0-1 +α sur (T n0 -T n0-1 ),0<α sur ≤1 (23)

[0113] T n0 is the temperature field of the partition in the n0th outer iteration round, α sur is the relaxation factor. When this parameter is less than 1, it is a successive sub-relaxation iteration.

[0114] In step 6), the parahydrogen concentration value stored in the fluid temperature node is determined point by point along the +x direction according to equations (28) and (30). i,k and T i,k+1 The defined fluid unit is known to and T i,k , the specific steps are as follows:

[0115] a. Estimated outlet parahydrogen concentration

[0116] b. Calculation (calculate Required T i,k+1 (derived from the fluid temperature field determined by the previous external iteration round), and then calculated according to equations (26) and (27):

[0117] c. Calculate according to formula (28) and (30) and T i,k+1 ;

[0118] d. Compare newly calculated and old If the residual does not meet the convergence condition, continue to perform a new round of calculation from step b until The residuals of meet the convergence conditions;

[0119] When it is determined and T i,k+1 , continue to determine the outlet of the next fluid unit and T i,k+2 , until the outlet of the fluid channel numbered i is calculated and T i,N+1 .

[0120] Compared with the prior art, the present invention has the following beneficial effects:

[0121] 1. The technology of integrating hydrogen secondary catalytic conversion with flow heat exchange in a plate-fin heat exchanger adopted in the present invention can increase heat exchange, reduce equipment volume, optimize heat exchange process, and reduce the complexity of system components.

[0122] 2. The present invention is a calibration method for plate-fin heat exchange coupled with hydrogen paracatalytic conversion. It takes into account the physical property changes of the fluid and the axial heat conduction effect of the partition. It introduces a porous medium pressure drop source term, a hydrogen paracatalytic conversion heat source term, and a parahydrogen component transport process into the channel where the catalytic conversion occurs. The calculation speed is fast and the accuracy is high.

[0123] 3. The method of the present invention can be written into a general calculation program that can perform calibration calculations for any arrangement of multi-stream countercurrent plate-fin heat exchangers coupled with hydrogen paracatalytic conversion. This overcomes the shortcomings of traditional CFD, which is incapable of numerically simulating plate-fin heat exchangers with hundreds of layers of fin channels. It also overcomes the problem that traditional plate-fin heat exchanger design software such as HTRI, ASPEN MUSE, and APSEN EDR cannot incorporate hydrogen paracatalytic conversion into heat exchanger calculations. The calculation results include the distribution of the pressure field, temperature field, and parahydrogen concentration field of the entire plate-fin heat exchanger. BRIEF DESCRIPTION OF THE DRAWINGS

[0124] Figure 1 Schematic diagram of the three-dimensional structure of a multi-stream countercurrent plate-fin heat exchanger filled with hydrogen secondary conversion catalyst (using serrated fins).

[0125] Figure 2 This is a simplified two-dimensional model of a plate-fin heat exchanger and a schematic diagram of the local temperature node layout.

[0126] Figure 3 Schematic diagram of heat conduction in the fluid channels of a plate-fin heat exchanger filled with catalyst.

[0127] Figure 4 Schematic diagram of block iterative calculation of the plate-fin heat exchanger partition temperature node.

[0128] Figure 5 It is a schematic diagram of the process of the present invention.

[0129] Figure 6 Schematic diagram of the temperature distribution of the fluid along the flow direction in the three-layer channel calculated in an embodiment of the present invention.

[0130] Figure 7 Schematic diagram of the outlet temperature and outlet parahydrogen concentration of each layer of fluid channels calculated in an embodiment of the present invention.

[0131] Figure 8 Schematic diagram of the distribution of parahydrogen concentration and equilibrium hydrogen concentration along the flow direction in a layer of fluid channels calculated according to an embodiment of the present invention. DETAILED DESCRIPTION

[0132] The present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0133] A calibration method for a low-temperature plate-fin heat exchanger coupled with a hydrogen secondary catalytic conversion process comprises the following steps:

[0134] 1) Given parameters:

[0135] Select a three-stream plate-fin heat exchanger. The heat exchanger is arranged in the form of (CBCBCAA*16)CBCBC. There are a total of 117 layers of fluid channels, such as Figure 1 As shown, the length L of the heat exchanger core is 800mm and the width W is 550mm; fluid A is hydrogen, and 32 layers of fluid channels are arranged, which are filled with catalysts for hydrogen ortho-para conversion, flowing along the +x direction, with a mass flow rate of 0.064kg / s, a para-hydrogen concentration of 45% at the inlet, a temperature of 78.2K, and a pressure of 2.4MPa. The structural parameters of the fins in fluid channel A are: perforated fins, fin height h i is 9.5mm, wing pitch s i is 1.4mm, wing thickness t i is 0.3mm, porosity φ i is 0.05; fluid B is helium, with 34 layers of fluid channels arranged, flowing along the +x direction, with a mass flow rate of 1.285 kg / s, an inlet temperature of 78.5 K, and an inlet pressure of 2.1 MPa. The structural parameters of the fins in fluid channel B are: serrated fins, fin height h i is 6.3mm, wing pitch s i is 1.2mm, wing thickness t i is 0.2mm, pitch l iThe C fluid is helium, with 51 layers of fluid channels arranged, flowing along the -x direction, with a mass flow rate of 1.285 kg / s, an inlet temperature of 67.5 K, and an inlet pressure of 0.42 MPa. The structural parameters of the fins in the C fluid channel are: serrated fins, fin height h i is 6.35mm, wing pitch s i is 1.2mm, wing thickness t i is 0.2mm, pitch l i It is 3mm;

[0136] The parameters of the hydrogen secondary reforming catalyst particles filled in fluid channel A are: porosity ε is 0.35, porosity is defined as the percentage of the void volume between particles in the packing layer to the total volume, particle diameter is 0.371 mm, and the thermal conductivity of the catalyst particles is 0.58 W·m -1 ·K -1 ;

[0137] 2) Model simplification and initialization:

[0138] In order to facilitate the construction of the model, it is assumed that the fluid in the fluid channel is uniformly distributed along the z direction, and the three-dimensional model ( Figure 1 ) is simplified to a two-dimensional plane ( Figure 2 ) within the research. In addition, only steady-state was considered in the calculations, and the heat exchanger was assumed to be thermally insulated from the surrounding environment. Heat conduction was ignored in the fluid calculations. The model construction considered the fluid domain, the solid domains of the fins and baffles, and the solid domain of the catalyst particles. In the channel where the hydrogen ortho- and para-catalytic reactions occurred, the following assumptions were made: hydrogen is a single-phase compressible gas composed of ortho- and para-hydrogen, the catalyst packing in the channel is uniform, and the particles are approximately spherical with a uniform diameter.

[0139] like Figure 2 As shown, there are n+1 layers of partitions and n layers of fin channels. Figure 2 N nodes are arranged in the x direction in the fluid channel, and N+1 nodes are arranged along the flow direction. The temperature nodes T w,i,k Define a diaphragm element, and each of the two fluid temperature nodes T in the fluid channel i,k and T i,k+1 Define a fluid unit, the fluid temperature node is also used to store pressure data P i,k and P i,k+1 If the secondary catalytic conversion of hydrogen also occurs in the fluid channel, the secondary hydrogen concentration Pa is stored at the corresponding temperature node. i,k and Pa i,k+1 In a plate-fin heat exchanger, there are a total of N×(n+1) baffle units and n×N fluid units ( Figure 2 shown);

[0140] The entire calculation process is iterative, that is, each outer iteration round calculates a new field based on the old field obtained in the previous outer iteration round. Therefore, before the calculation of the first outer iteration round, an initialized temperature field, pressure field, and parahydrogen concentration field need to be given; the fluid temperature node in each layer of the fluid channel is initialized to the inlet temperature given in step 1), and the pressure value stored in the temperature node is initialized to the inlet pressure given in step 1). In the fluid channel where hydrogen positive-paracatalytic conversion occurs, the temperature node also stores the parahydrogen concentration value, which is initialized to the inlet parahydrogen concentration given in step 1);

[0141] 3) Fluid physical property fitting:

[0142] The required fluid properties include specific heat capacity, thermal conductivity, dynamic viscosity, density, and Prandtl number. The physical property data are extracted from the NIST database. Generally, the relationship between the density of the fluid and temperature and pressure is fitted on a spline surface. The physical properties of other types of fluids are fitted on spline curves with temperature. The thermal conductivity of solids also has a functional relationship with temperature and can be fitted into a thermal conductivity-temperature spline curve. Finally, all the parameters of the fitted spline curves and spline surfaces are loaded into memory for subsequent calculations.

[0143] In the fin channel where hydrogen ortho-paracatalytic conversion occurs, the physical properties of ortho-hydrogen and para-hydrogen need to be fitted separately. In subsequent calculations, if the physical properties of hydrogen fluid at a certain para-hydrogen concentration need to be calculated, the physical properties of ortho-hydrogen and para-hydrogen that have been fitted can be interpolated and calculated.

[0144] 4) Calculate the pressure field distribution in the fluid channel of the plate-fin heat exchanger:

[0145] To obtain the pressure field distribution in the fluid channel, it is necessary to first calculate the pressure drop of each fluid unit in the fluid channel. Depending on whether the secondary catalytic conversion of hydrogen occurs in the fluid channel, there are two ways to calculate the pressure drop:

[0146]

[0147] Formula (1) is used to calculate the fluid temperature node T i,k and T i,k+1 The pressure drop in the defined fluid unit, if the fluid channel in this layer is not filled with catalyst particles, that is, no hydrogen secondary catalytic conversion occurs, then the first method is adopted to calculate the pressure drop through the friction factor f in the fluid unit. i,k To calculate, G i is the mass flow rate of the fluid in each layer of fin channel, Δl x is the length of the fluid unit along the flow direction, D i The equivalent diameter of the fin structure in the fluid channel of this layer, ρ i,k+1 and ρi,k By the fluid temperature node T i,k+1 and T i,k And the pressure value P stored in it i,k+1 and P i,k The density determined, is the density of the fluid unit, f i,k and It can be determined based on the qualitative temperature and qualitative pressure in the fluid unit:

[0148]

[0149] Straight fins, perforated fins, or serrated fins may be arranged in the fluid channel. The friction factor of serrated fins is calculated using the following method:

[0150]

[0151] Where α, θ, and γ are:

[0152]

[0153] For perforated fins, the friction factor is calculated as follows:

[0154] lnf i,k =28.78906-12.31399lnRe i,k +1.565191(lnRe i,k ) 2 -0.06736098(lnRe i,k ) 3 (5)

[0155] For straight fins, the friction factor is calculated as follows:

[0156]

[0157] Re i,k is the Reynolds number:

[0158]

[0159] The fluid temperature node T i,k and T i,k+1 The dynamic viscosity of the defined fluid unit is determined according to the qualitative temperature and qualitative pressure therein;

[0160] If the fluid channel of this layer is filled with catalyst particles for hydrogen secondary catalytic conversion, the second method in formula (1) is needed to calculate the pressure drop in the fluid unit. is the superficial velocity of the fluid in the fluid unit, that is, the average velocity in the fluid unit calculated without considering the catalyst particles filled in the fluid channel, α is the permeability of the catalyst particles, and C1 is the inertial resistance coefficient of the catalyst particles. The calculation methods of the two parameters are:

[0161]

[0162]

[0163] Where d p is the diameter of the catalyst particle, ε is the porosity;

[0164] The pressure drop of each fluid unit is calculated according to formula (1). Given the inlet pressure of each layer of fluid channel, the pressure value stored in each temperature node can be updated point by point along the flow direction; the updated pressure field distribution is compared with the old pressure field distribution, and the pressure field residual is calculated. If the residual does not meet the convergence condition, step 2) is repeated, that is, the pressure field is continuously updated based on the old pressure field until the pressure field residual meets the convergence condition, and step 5) is entered. The pressure field distribution updated in this step is used as the qualitative pressure of each fluid unit in the calculation of subsequent steps;

[0165] 5) Calculation of temperature field distribution of plate-fin heat exchanger baffle:

[0166] When catalyst particles are filled in the fluid channels adjacent to the partition, they expand the heat exchange area between the fluid and the solid surface. The heat transferred from the fluid to the partition is divided into three parts: one part is the fluid heat transferred through the fin surface to the fin solid domain, and then the heat is transferred from the fin solid domain to the partition solid domain; another part is the fluid heat transferred through the catalyst surface to the catalyst solid domain, and then through the catalyst solid domain to the partition solid domain; the third part is the fluid heat transferred directly from the partition surface to the partition solid domain. If the fluid channels adjacent to the partition are not filled with catalyst, the heat transferred from the fluid to the partition only includes the first and third parts.

[0167] For the solid domain of fins and catalyst particles, in order to simplify the calculation, the heat conduction along the axial direction (+x or -x direction) is not considered, as shown in Figure 3 As shown, then the temperature node T i,k and T i,k+1 The one-dimensional heat conduction equation for the fin and catalyst particle solid domains in the defined fluid unit is:

[0168]

[0169] Equation (10) can be used to calculate the temperature T of the fin and catalyst solid area s,i,k Temperature distribution along the y direction, λs,i,k and λ p,i,k are the thermal conductivity of the fin and catalyst particle solid domain, is the heat transfer source term in the one-dimensional heat conduction equation, which includes the heat transfer between the fin surface and the fluid, and the heat transfer between the catalyst particle surface and the fluid:

[0170]

[0171] In formula (11), s,i,k and λ p,i,k The fluid temperature node T i,k and T i,k+1 The thermal conductivity of the solid domain of the fins and catalyst particles in the defined fluid cell, α p is the specific surface area of the catalyst particles. Assuming that the catalyst particles are uniform in size and are all perfect spheres, the calculation method is:

[0172] α p =6(1-ε) / d p (12)

[0173] h c,i,k is the convective heat transfer coefficient between the fluid and solid domain surfaces in the fluid channel. When the fluid channel is not filled with a catalyst, h c,i,k The calculation method is:

[0174]

[0175] Where Pr i,k is determined by the temperature node T i,k and T i,k+1 The Prandtl number of the defined fluid cell, λ fl,i.k is the thermal conductivity of the fluid unit, j i,k The evaluation index used in this fluid unit to evaluate the heat transfer capacity of the fin surface is the following j-factor correlation formula for serrated fins:

[0176]

[0177] For perforated fins, the correlation is:

[0178] lnj i,k =34.57583-15.92678lnRe i,k +2.137607(lnRe i,k ) 2 -0.09544151×(lnRe i,k ) 3 (15)

[0179] For straight fins, the correlation is:

[0180]

[0181] If the fluid channel is filled with catalyst particles, the convection heat transfer coefficient h on the solid domain surface in the fluid unit is c,i,k The calculation method is:

[0182]

[0183] Where Re p It is the Reynolds number calculated by taking the particle diameter as the characteristic size and the superficial velocity in the fluid unit;

[0184] According to equations (10) and (11), it can be deduced that the temperature node T i,k and T i,k+1 Temperature distribution of the solid domain along the y direction in the defined fluid cell:

[0185]

[0186] Where sh and ch are hyperbolic cosine and hyperbolic sine respectively, θ i,k ,θ i-i,k and θ i-(i+1),k is the excess temperature of the fin, fin top and root; m e,i,k ,θ i,k ,θ i-i,k and θ i-(i+1),j It is expressed as follows:

[0187]

[0188] Once the temperature distribution in the solid domain contained in the fluid unit is obtained, the heat conduction between the solid domain contained in the fluid unit and the solid domain of the partition can be obtained by Fourier's law of heat conduction. This is actually the first two parts of the three heat exchange components discussed at the beginning of step 5). For a partition unit, in addition to these three parts of heat from adjacent fluid units, there is a fourth part of heat, which is the heat conduction between the partition unit and the adjacent partition units. These four parts of heat in the partition unit satisfy the law of conservation of energy. The algebraic equations for calculating the partition temperature node are derived as follows:

[0189]

[0190]

[0191] A i-1,k+1 =A i-1,k

[0192]

[0193] A i,k+1 =A i,k

[0194]

[0195]

[0196]

[0197]

[0198] A w,i,k =A i-1,k +A i-1,k+1 +A i,k +A i,k+1 +A w,i-1,k +A w,i+1,k +A w,i,k-1 +A w,i,k+1 (twenty one)

[0199] In formula (21), t sp is the baffle thickness of the plate-fin heat exchanger, λ w is the thermal conductivity of the partition. In the formula, the subscripts (i, k) are marked outside the brackets. The subscripts of the parameters in the brackets are divided into three cases. If the parameter Δl x ,W,t sp ,ε,α p If the parameters are s, t, h, the subscript is i, and if the parameters are λ, w ,λ s ,h c ,m e etc., the subscript is (i,k); if the subscript outside the brackets is (i-1,k) or (i,k-1), the parameters in the brackets are marked in the same way. In short, the marking method of the parameters in the brackets is the same as that of the parameters appearing in other formulas in the text; λ w,i,k is the partition temperature node T w,i,k and T w,i,k+1 The thermal conductivity of the interface, if k = N, is directly calculated based on the partition temperature node T w,i,N To determine λ w,i,N ,λ w,i,k-1 is the partition temperature node T w,i,k-1 and T w,i,k The thermal conductivity of the interface, if k = 1, is directly calculated based on the partition temperature node T w,i,1 To determine λ w,i,1 In addition, m in formula (21) e,i,k The calculation of can be divided into two cases according to formula (19);

[0200] For each partition element defined by the partition temperature node, it is required to solve the algebraic equation shown in Equation (20). A total of N×(n+1) equations need to be solved. In order to speed up the calculation of the partition temperature nodes, the line iteration method combined with the Gauss-Seidel method is used to decompose the partition solid domain solution into multiple lines parallel to the y-axis. Each line contains only one layer of temperature nodes. The values of each node on the same line are obtained by directly solving the algebraic equation. Figure 4 The block of diaphragm temperature nodes on line AA' is shown in Figure 1. The values of nodes within the same block are implicitly linked, but the progression from one line to another along the +x direction is performed iteratively:

[0201]

[0202] Where T P represents the required bulkhead temperature node, T N ,T W ,T S ,T E (north, west, south, east) are the adjacent partition temperature nodes, n0 can be understood as the identifier of the outer iteration round, b represents the heat transfer source term that cannot be classified as the influence of the adjacent partition temperature node. Since the scanning direction is along the +x direction, that is, from left to right, T N ,T W ,T S The updated shelf temperature is used, and T E The temperature value updated in the previous outer iteration round is used. For each partition temperature node on a line, a system of algebraic equations in the form of equation (22) must be solved. This line contains a total of n + 1 equations. The Thomas algorithm based on the Gauss elimination method can be used to solve the system of equations simultaneously (see Chapter 4.4.1.1 of the second edition of Tao Wenquan's "Numerical Heat Transfer"). For the partition temperature nodes, it is sometimes necessary to introduce successive sub-relaxation:

[0203] T n0 =T n0-1 +α sur (T n0 -T n0-1 ),0<α sur ≤1 (23)

[0204] T n0 is the temperature field of the partition in the n0th outer iteration round, α sur is the relaxation factor. When this parameter is less than 1, it is a successive sub-relaxation iteration.

[0205] When all the partition temperature nodes are updated from left to right, you can proceed to step 6);

[0206] 6) Calculation of temperature distribution and parahydrogen concentration distribution in the fluid channel:

[0207] If the hydrogen secondary catalytic conversion process occurs in the fluid channel numbered i, the temperature node T i,k and T i,k+1 The heat of secondary catalytic conversion of hydrogen in the defined fluid unit is:

[0208]

[0209] Where S o-p,i,k The fluid temperature node T i,k and T i,k+1 The heat source term for the secondary conversion of hydrogen in the defined fluid unit, M is the molar mass of hydrogen molecules, ΔH i,k is the heat of reaction of the secondary catalytic reaction of hydrogen in the fluid unit. Generally speaking, the heat of reaction is a function of temperature:

[0210]

[0211] The rate of hydrogen conversion per unit volume Calculated according to the Elovich model:

[0212]

[0213] Where, P c is the critical pressure of hydrogen, is determined by the temperature node T i,k and T i,k+1 The parahydrogen concentration in the defined fluid cell, is the reaction rate constant of the Elovich calculation model; a, b1, c, and d are the fitting coefficients of the experimental data, and their values are 1.0924, 59.7 mol·m -3 ·s -1 , -253.9 mol·m -3 ·s -1 , -11.6 mol·m -3 ·s -1 . To balance the parahydrogen content in hydrogen, the calculation method is:

[0214]

[0215] Where T c is the critical temperature of hydrogen;

[0216] The calculation of the temperature nodes in the fluid channel is divided into two cases. If the hydrogen secondary catalytic conversion process does not occur in the fluid channel, the fluid in the fluid unit exchanges heat with the surface of the partition solid domain and the surface of the fin solid domain. If the hydrogen secondary catalytic reaction occurs in the fluid channel, in addition to the two aforementioned heat sources, the fluid in the fluid unit also exchanges heat with the surface of the catalyst solid particles, and the hydrogen secondary catalytic conversion heat also provides heat for the fluid. These heat sources in the fluid unit should satisfy the energy conservation law. Combined with the previous equations (18) and (19), the equation group for calculating the fluid temperature node can be derived. If the fluid flows along the +x direction, the equation group can be derived as follows:

[0217]

[0218] The influence coefficient B in the formula i,k , B w,i,k ,B w,i+1,k ,B i,k+1 The calculation method is:

[0219]

[0220]

[0221] B i,k+1 =B i,k +B w,i,k +B w,i+1,k (29)

[0222] In the formula, the subscript outside the brackets is (i,k), then the m inside the brackets is m i , represents the mass flow rate in the fluid channel numbered i, c p Can be written as c p,i,k , is the temperature node T i,k and T i,k+1 The specific heat capacity of the fluid unit defined by the equation is: If hydrogen secondary catalytic conversion occurs in the fluid channel, then m in equation (29) e,i,k The calculation method for the case where the fluid channel is filled with a catalyst as described in formula (19) should be used. In this case, formula (28) is not closed, and the component transport equation for parahydrogen needs to be added:

[0223]

[0224] If the fluid flows along the +x direction and the para-catalytic conversion of hydrogen occurs in the fluid channel, then according to equations (28) and (30), the para-hydrogen concentration value stored in the fluid temperature node can be determined point by point along the +x direction, so that the para-hydrogen concentration value stored in the fluid temperature node T i,k and T i,k+1 As an example of the defined fluid unit, the specific steps are as follows (it is known that and T i,k):

[0225] a. Estimated outlet parahydrogen concentration

[0226] b. Calculation (calculate Required T i,k+1 (derived from the fluid temperature field determined by the previous external iteration round), and then calculated according to equations (26) and (27):

[0227] c. Calculate according to formula (28) and (30) and T i,k+1 ;

[0228] d. Compare newly calculated and old If the residual does not meet the convergence condition, continue to perform a new round of calculation from step b until The residuals of meet the convergence conditions;

[0229] When it is determined and T i,k+1 , you can continue to determine the outlet of the next fluid unit and T i,k+2 , until the outlet of the fluid channel numbered i is calculated and T i,N+1 ;

[0230] If the fluid channel is along the -x direction, the equations for calculating the temperature node of the fluid element are:

[0231]

[0232] The calculation method of the influence coefficient is:

[0233]

[0234]

[0235] B i,k =B i,k+1 +B w,i,k +B w,i+1,k (32)

[0236] The component transport equation of parahydrogen in the fluid unit is:

[0237]

[0238] By using the general iterative method introduced above, equations (31) and (33) can determine the temperature field and parahydrogen concentration field distribution in the fluid channel where the flow is in the -x direction;

[0239] 7) Calculation of temperature field residual and parahydrogen concentration field residual:

[0240] Compare the partition temperature field distribution calculated in step 5) with the partition temperature field distribution determined in the previous external iteration round, and calculate the partition temperature field residual. Compare the fluid temperature field distribution calculated in step 6) with the fluid temperature field distribution determined in the previous external iteration round, and calculate the fluid temperature field residual. Take the maximum value of the partition temperature field residual and the fluid temperature field residual as the temperature field residual; Take the parahydrogen concentration field calculated in step 6) with the parahydrogen concentration field distribution determined in the previous external iteration round, and calculate the parahydrogen concentration field residual. If one of the temperature field residual and the parahydrogen concentration field residual does not meet the convergence condition, continue to repeat steps 4) to 7) as a new round of external iteration until both residuals meet the convergence condition, and then output the calculation results. The calculation results include the overall temperature field (partition temperature field and fluid temperature field) distribution of the plate-fin heat exchanger, the pressure field distribution in the fluid channel, and the parahydrogen concentration distribution.

[0241] The process of the present invention is shown in Figure 5 , step 4) to step 7) is a complete outer iteration round, wherein step 4) is the inner iteration calculation of the pressure field, step 5) is to use the line iteration method combined with the Gauss-Seidel method to calculate the temperature field in the partition, and step 6) is to step-by-step calculate the distribution of the fluid temperature field and the parahydrogen concentration field. In each calculation step of step 6), it is necessary to first estimate the parahydrogen concentration, and then introduce the idea of the general iterative method to calculate the outlet parahydrogen concentration and outlet temperature of each fluid unit.

[0242] Figure 6 、 Figure 7 and Figure 8 are some calculation results of this embodiment, Figure 6 The temperature distribution of the fluid in the three-layer fluid channel of the plate-fin heat exchanger is shown in Figure 1. Figure 7 is the calculated outlet temperature of each layer of fluid channel (a total of 117 layers), and the outlet parahydrogen concentration of fluid A, Figure 8 It is the distribution of parahydrogen concentration and equilibrium hydrogen concentration in the sixth layer fluid channel of the plate-fin heat exchanger along the fluid direction.

Claims

1. A calibration method for a low-temperature plate-fin heat exchanger coupled with a hydrogen secondary catalytic conversion process, characterized in that: The following steps are involved: 1) Given parameters: Select a plate-fin heat exchanger with multiple flows, and determine the overall structural parameters, operating parameters, and fin structural parameters of the plate-fin heat exchanger for each flow; The overall structural parameters include the arrangement of the heat exchanger, the length L and width W of the heat exchanger core; The operating parameters include the fluid of the multi-stream, the fluid channel layer, the mass flow rate, the inlet temperature, the inlet pressure, and the inlet parahydrogen concentration; The structural parameters of the fin are: when straight fins are selected, the parameter is fin height h i , wing span s i 、Wing thickness t i ;When selecting perforated fins, the parameter is fin height h i , wing span s i 、Wing thickness t i , porosity φ i When selecting serrated fins, the parameter is fin height h i , wing span s i 、Wing thickness t i , pitch l i ; The parameters of the hydrogen secondary reforming catalyst particles filled in the fluid channel are: the porosity ε is 0.35, the porosity is defined as the percentage of the void volume between particles in the packing layer to the total volume, the particle diameter is 0.371 mm, and the thermal conductivity of the catalyst particles is 0.58 W·m -1 ·K -1 ; 2) Model simplification and initialization: The three-dimensional model was simplified to a two-dimensional plane, assuming that the fluid in the fluid channel is uniformly distributed along the z-direction. Furthermore, only steady-state conditions were considered in the calculations, and the heat exchanger was assumed to be thermally insulated from the surrounding environment. Thermal conductivity was ignored in the fluid calculations. The model construction considered the fluid domain, the solid domains of the fins and baffles, and the solid domain of the catalyst particles. In the channel where the hydrogen ortho- and para-catalytic reactions occurred, the following assumptions were made: hydrogen is a single-phase compressible gas composed of ortho- and para-hydrogen; the catalyst packing in the channel is uniform, and the particles are approximately spherical with a uniform diameter. There are n+1 layers of partitions and n layers of fin channels. N nodes are arranged along the flow direction in the partitions, and N+1 nodes are arranged along the flow direction in the fluid channels. The temperature of each partition node T in the partition is w,i,k Define a diaphragm element, and each of the two fluid temperature nodes T in the fluid channel i,k and T i,k+1 Define a fluid unit, the fluid temperature node is also used to store pressure data P i,k and P i,k+1 If the secondary catalytic conversion of hydrogen also occurs in the fluid channel, the secondary hydrogen concentration is stored at the corresponding temperature node. and In a plate-fin heat exchanger, a total of N×(n+1) baffle units and n×N fluid units are set; The entire calculation process is iterative. Before the calculation of the first outer iteration round, an initialized temperature field, pressure field, and parahydrogen concentration field need to be given. The fluid temperature node in each layer of the fluid channel is initialized to the inlet temperature given in step 1), and the pressure value stored in the temperature node is initialized to the inlet pressure given in step 1). In the fluid channel where hydrogen positive-paracatalytic conversion occurs, the parahydrogen concentration value is also stored in the temperature node, initialized to the inlet parahydrogen concentration given in step 1). 3) Fluid physical property fitting: The required fluid properties include specific heat capacity, thermal conductivity, dynamic viscosity, density, and Prandtl number. The physical property data are extracted from the NIST database. The relationship between the fluid density, temperature, and pressure is fitted on a spline surface. The physical properties of other types of fluids are fitted with temperature on spline curves. For the thermal conductivity of solids, a thermal conductivity-temperature spline curve is fitted. Finally, all the parameters of the fitted spline curves and spline surfaces are loaded into memory. In the fin channel where hydrogen ortho-paracatalytic conversion occurs, the physical properties of ortho-hydrogen and para-hydrogen need to be fitted separately. If the physical properties of hydrogen fluid at a certain para-hydrogen concentration need to be calculated, the physical properties of the fitted ortho-hydrogen and para-hydrogen are interpolated and calculated. 4) Calculate the pressure field distribution in the fluid channel of the plate-fin heat exchanger: First, the pressure drop of each fluid unit in the fluid channel is calculated. Depending on whether the secondary catalytic conversion of hydrogen occurs in the fluid channel, the pressure drop calculation method is divided into two types: Formula (1) is used to calculate the fluid temperature node T i,k and T i,k+1 The pressure drop in the defined fluid unit, if the fluid channel in this layer is not filled with catalyst particles, that is, no hydrogen secondary catalytic conversion occurs, then the first method is adopted to calculate the pressure drop through the friction factor f in the fluid unit. i,k To calculate, G i is the mass flow rate of the fluid in each layer of fin channel, Δl x is the length of the fluid unit along the flow direction, D i The equivalent diameter of the fin structure in the fluid channel of this layer, ρ i,k+1 and ρ i,k By the fluid temperature node T i,k+1 and T i,k And the pressure value P stored in it i,k+1 and P i,k The density determined, is the density of the fluid unit, f i,k and Determined based on the qualitative temperature and qualitative pressure in the fluid unit: Arrange straight fins, perforated fins or serrated fins in the fluid channel. The friction factor of serrated fins is calculated using the following method: Where α, δ, and γ are: For perforated fins, the friction factor is calculated as follows: lnf i,k =28.78906-12.31399lnRe i,k +1.565191(lnRe i,k ) 2 -0.06736098(lnRe i,k ) 3 (5) For straight fins, the friction factor is calculated as follows: Re i,k is the Reynolds number: The fluid temperature node T i,k and T i,k+1 The dynamic viscosity of the defined fluid unit is determined according to the qualitative temperature and qualitative pressure therein; If the fluid channel of this layer is filled with catalyst particles for hydrogen secondary catalytic conversion, the second method in formula (1) is needed to calculate the pressure drop in the fluid unit. is the superficial velocity of the fluid in the fluid unit, that is, the average velocity in the fluid unit calculated without considering the catalyst particles filled in the fluid channel, α c is the permeability of the catalyst particles, C1 is the inertial resistance coefficient of the catalyst particles, and the calculation methods of the two parameters are: Where d p is the diameter of the catalyst particle, ε is the porosity; The pressure drop of each fluid unit is calculated according to formula (1). Given the inlet pressure of each layer of fluid channel, the pressure value stored in each temperature node is updated point by point along the flow direction; the updated pressure field distribution is compared with the old pressure field distribution, and the pressure field residual is calculated. If the residual does not meet the convergence condition, step 2) is repeated, that is, the pressure field is continuously updated based on the old pressure field until the pressure field residual meets the convergence condition, and step 5) is entered. The pressure field distribution updated in this step is used as the qualitative pressure of each fluid unit in the calculation of subsequent steps; 5) Calculation of temperature field distribution of plate-fin heat exchanger baffle: When the fluid channels adjacent to the partition are filled with catalyst particles, the heat transferred from the fluid to the partition is divided into three parts: one part is the heat transferred from the fluid to the fin solid domain through the fin surface, and then the heat is transferred to the partition solid domain through the fin solid domain; the other part is the heat transferred from the fluid to the catalyst solid domain through the catalyst surface, and then transferred to the partition solid domain through the catalyst solid domain; The third part is the heat transferred from the fluid directly to the solid domain of the separator through the surface of the separator. If the fluid channels adjacent to the separator are not filled with catalyst particles, the heat transferred from the fluid to the separator only includes the first and third parts. For the solid domain of fins and catalyst particles, in order to simplify the calculation, the heat conduction along the axial direction (+x or -x) is not considered, so the temperature node T i,k and T i,k+1 The one-dimensional heat conduction equation for the fin and catalyst particle solid domains in the defined fluid unit is: Formula (10) is used to calculate the temperature T of the fin and catalyst solid area s,i,k Temperature distribution along the y direction, λ s,i,k and λ p,i,k are the thermal conductivity of the fin and catalyst particle solid domain, is the heat transfer source term in the one-dimensional heat conduction equation, which includes the heat transfer between the fin surface and the fluid, and the heat transfer between the catalyst particle surface and the fluid: In formula (11), s,i,k and λ p,i,k The fluid temperature node T i,k and T i,k+1 The thermal conductivity of the solid domain of the fins and catalyst particles in the defined fluid cell, α p is the specific surface area of the catalyst particles. Assuming that the catalyst particles are uniform in size and are all perfect spheres, the calculation method is: a p =6(1-e) / d p (12) h c,i,k is the convective heat transfer coefficient between the fluid and solid domain surfaces in the fluid channel. When the fluid channel is not filled with a catalyst, h c,i,k The calculation method is: Where Pr i,k is determined by the temperature node T i,k and T i,k+1 The Prandtl number of the defined fluid cell, λ fl,i.k is the thermal conductivity of the fluid unit, j i,k The evaluation index used in this fluid unit to evaluate the heat transfer capacity of the fin surface is the following j-factor correlation formula for serrated fins: For perforated fins, the correlation is: lnj i,k =34.57583-15.92678lnRe i,k +2.137607(lnRe i,k ) 2 -0.09544151×(lnRe i,k ) 3 (15) For straight fins, the correlation is: If the fluid channel is filled with catalyst particles, the convection heat transfer coefficient h on the solid domain surface in the fluid unit is c,i,k The calculation method is: Where Re p It is the Reynolds number calculated by taking the particle diameter as the characteristic size and the superficial velocity in the fluid unit; According to equations (10) and (11), we can derive the temperature node T i,k and T i,k+1 Temperature distribution of the solid domain along the y direction in the defined fluid cell: Where sh and ch are hyperbolic cosine and hyperbolic sine respectively, θ i,k ,θ i-i,k and θ i-(i+1),k is the excess temperature of the fin, fin top and root; m e,i,k ,θ i,k ,θ i-i,k and θ i-(i+1),k It is expressed as follows: When the temperature distribution in the solid domain contained in the fluid unit is obtained, the heat conduction between the solid domain contained in the fluid unit and the partition solid domain is obtained by Fourier's law of heat conduction. For a partition unit, the heat from the adjacent fluid unit also has a fourth part of heat, which is the heat conduction between the partition unit and the adjacent partition units. These four parts of heat in the partition unit satisfy the law of energy conservation. The algebraic equations for calculating the partition temperature node are derived as follows: A i-1,k+1 =A i-1,k A i,k+1 =A i,k A w,i,k =A i-1,k +A i-1,k+1 +A i,k +A i,k+1 +A w,i-1,k +A w,i+1,k +A w,i,k-1 +A w,i,k+1 (21) In formula (21), t sp is the baffle thickness of the plate-fin heat exchanger, λ w is the thermal conductivity of the partition. In the formula, the subscripts (i, k) are marked outside the brackets. The subscripts of the parameters in the brackets are divided into three cases. If the parameter Δl x ,W,t sp ,ε,α p , then no subscript is required. If it is a parameter s, t, h, the subscript is i. If it is a parameter λ w ,λ s ,h c ,m e , the subscript is (i,k); if the subscript outside the brackets is (i-1,k) or (i,k-1), the parameters in the brackets are marked in the same way. In short, the parameters in the brackets are marked in the same way as those in other formulas in the text; λ w,i,k is the partition temperature node T w,i,k and T w,i,k+1 The thermal conductivity of the interface, if k = N, is directly calculated based on the partition temperature node T w,i,N To determine λ w,i,N ,λ w,i,k-1 is the partition temperature node T w,i,k-1 and T w,i,k The thermal conductivity of the interface, if k = 1, is directly calculated based on the partition temperature node T w,i,1 To determine λ w,i,1 In addition, m in formula (21) e,i,k The calculation of can be divided into two cases according to formula (19); For each partition element defined by the partition temperature node, it is required to solve the algebraic equation shown in Equation (20). A total of N×(n+1) equations need to be solved. In order to speed up the calculation of the partition temperature nodes, the line iteration method combined with the Gauss-Seidel method is used. The partition solid domain solution is decomposed into multiple lines parallel to the y-axis. Each line contains only one layer of temperature nodes. The values of each node on the same line are obtained by directly solving the algebraic equation. The values of each node in the block formed by the partition temperature nodes on line AA' are implicitly connected, but the advancement from one line to another along the +x direction is performed iteratively: Where T P represents the required bulkhead temperature node, T N ,T W ,T S ,T E is the adjacent partition temperature node, n0 is understood as the identifier of the outer iteration round, b represents the heat transfer source term that cannot be classified as the influence of the adjacent partition temperature node. Since the scanning direction is along the +x direction, that is, from left to right, T N ,T W ,T S The updated shelf temperature is used, and T E The temperature value updated in the previous outer iteration round is used. For each partition temperature node on a line, a set of algebraic equations in the form of equation (22) is required to be solved. There are a total of n+1 equations on this line. The Thomas algorithm based on the Gauss elimination method is used to solve the system of equations simultaneously. When all the partition temperature nodes are updated from left to right, go to step 6); 6) Calculation of temperature distribution and parahydrogen concentration distribution in the fluid channel: If the hydrogen secondary catalytic conversion process occurs in the fluid channel numbered i, the temperature node T i,k and T i,k+1 The heat of secondary catalytic conversion of hydrogen in the defined fluid unit is: Where S o-p,i,k The fluid temperature node T i,k and T i,k+1 The heat source term for the secondary conversion of hydrogen in the defined fluid unit, M is the molar mass of hydrogen molecules, ΔH i,k is the heat of reaction of the secondary catalytic reaction of hydrogen in the fluid unit, which is a function of temperature: The rate of hydrogen conversion per unit volume Calculated according to the Elovich model: Where, P c is the critical pressure of hydrogen, is determined by the temperature node T i,k and T i,k+1 The parahydrogen concentration in the defined fluid cell, is the reaction rate constant of the Elovich calculation model; a, b1, c, and d are the fitting coefficients of the experimental data, and their values are 1.0924, 59.7 mol·m -3 ·s -1 , -253.9 mol·m -3 ·s -1 , -11.6 mol·m -3 ·s -1 ; To balance the parahydrogen content in hydrogen, the calculation method is: Where T c is the critical temperature of hydrogen; The calculation of the temperature nodes in the fluid channel is divided into two cases. If the hydrogen secondary catalytic conversion process does not occur in the fluid channel, the fluid in the fluid unit exchanges heat with the surface of the partition solid domain and the surface of the fin solid domain. If the hydrogen secondary catalytic reaction occurs in the fluid channel, the fluid in the fluid unit also exchanges heat with the surface of the catalyst solid particles, and the hydrogen secondary catalytic conversion heat also provides heat for the fluid. These heat sources in the fluid unit satisfy the energy conservation law. Combined with the previous equations (18) and (19), the equation group for calculating the fluid temperature node is derived. If the fluid flows along the +x direction, the derived equation group is in the form of: The influence coefficient B in the formula i,k , B w,i,k ,B w,i+1,k ,B i,k+1 The calculation method is: B i,k+1 =B i,k +B w,i,k +B w,i+1,k (29) In the formula, the subscript outside the brackets is (i, k), and the subscript inside the brackets is m i , represents the mass flow rate in the fluid channel numbered i, c p Writing c p,i,k , is the temperature node T i,k and T i,k+1 The specific heat capacity of the defined fluid unit; If the secondary catalytic conversion of hydrogen occurs in the fluid channel, then m in Eq. (29) e,i,k The calculation method should be used when the catalyst is filled in the fluid channel of formula (19). In this case, formula (28) is not closed, and the component transport equation of parahydrogen needs to be added: If the fluid flows along the +x direction and the para-catalytic conversion of hydrogen occurs in the fluid channel, the para-hydrogen concentration value stored in the fluid temperature node is determined point by point along the +x direction according to equations (28) and (30); If the fluid channel is along the -x direction, the equations for calculating the temperature node of the fluid element are: The calculation method of the influence coefficient is: B i,k =B i,k+1 +B w,i,k +B w,i+1,k (32) The component transport equation of parahydrogen in the fluid unit is: By iterative method, equations (31) and (33) determine the temperature field and parahydrogen concentration field distribution in the fluid channel where the flow is in the -x direction; 7) Calculation of temperature field residual and parahydrogen concentration field residual: Compare the partition temperature field distribution calculated in step 5) with the partition temperature field distribution determined in the previous external iteration round, calculate the partition temperature field residual, compare the fluid temperature field distribution calculated in step 6) with the fluid temperature field distribution determined in the previous external iteration round, calculate the fluid temperature field residual, and take the maximum value of the partition temperature field residual and the fluid temperature field residual as the temperature field residual; compare the parahydrogen concentration field calculated in step 6) with the parahydrogen concentration field distribution determined in the previous external iteration round, and calculate the parahydrogen concentration field residual; if one of the temperature field residual and the parahydrogen concentration field residual does not meet the convergence condition, continue to repeat steps 4) to 7) as a new round of external iteration until both residuals meet the convergence condition, that is, output the calculation results, which include the overall temperature field distribution of the plate-fin heat exchanger, the pressure field distribution and the parahydrogen concentration distribution in the fluid channel; the temperature field includes the partition temperature field and the fluid temperature field.

2. The method according to claim 1, wherein: In step 5), for the partition temperature node, successive sub-relaxations are introduced as needed: T n0 =T n0-1 +α sur (T n0 -T n0-1 ),0<α sur ≤1 (23) T n0 is the temperature field of the partition in the n0th outer iteration round, α sur is the relaxation factor. When this parameter is less than 1, it is a successive sub-relaxation iteration.

3. The method according to claim 1, wherein: In step 6), the parahydrogen concentration value stored in the fluid temperature node is determined point by point along the +x direction according to equations (28) and (30). i,k and T i,k+1 The defined fluid unit is known to and T i,k , the specific steps are as follows: a. Estimated outlet parahydrogen concentration b. Calculation calculate Required T i,k+1 It comes from the fluid temperature field determined by the previous external iteration round, and is calculated according to equations (26) and (27): c. Calculate according to formula (28) and (30) and T i,k+1 ; d. Compare newly calculated and old If the residual does not meet the convergence condition, continue to perform a new round of calculation from step b until The residuals of meet the convergence conditions; When it is determined and T i,k+1 , continue to determine the outlet of the next fluid unit and T i,k+2 , until the outlet of the fluid channel numbered i is calculated and T i,N+1 .

Images