Multi-Objective Parameter Optimization of Pulse Tube Refrigerator Based on Kriging Metamodel and Non-Dominated Ranking Genetic Algorithms

Abstract

Structure parameters have an important influence on the refrigeration performance of pulse tube refrigerators. In this paper, a method combining the Kriging metamodel and Non-Dominated Sorting Genetic Algorithm II (NSGA II) is proposed to optimize the structure of regenerators and pulse tubes to obtain better cooling capacity. Firstly, the Kriging metamodel of the original pulse tube refrigerator CFD model is established to improve the iterative solution efficiency. On this basis, NSGA II was applied to the optimization iteration process to obtain the optimal and worst Pareto front solutions for cooling performance, the heat and mass transfer characteristics of which were further analyzed comparatively to reveal the influence mechanism of the structural parameters. The results show that the Kriging metamodel presents a prediction error of about 2.5%. A 31.24% drop in the minimum cooling temperature and a 31.7% increase in cooling capacity at 120 K are achieved after optimization, and the pressure drop loss at the regenerator and the vortex in the pulse tube caused by the structure parameter changes are the main factors influencing the whole cooling performance of the pulse tube refrigerators. The current study provides a scientific and efficient design method for miniature cryogenic refrigerators.

1. Introduction

The pulse tube refrigerator (PTR) proposed by Gifford and Longsworth [1] has been widely applied, especially in the field of spacecraft thermal management. Compared with other cryogenic refrigerators, PTR eliminates the moving parts at the cold end and, therefore, has the advantages of high reliability, good stability and simple construction. Afterwards, it has been continuously improved, and an orifice valve version [2], double inlet version [3], multi-bypass version [4] and inertance tube type pulse tube refrigerator [5] have appeared, and the refrigeration performance is continuously improving.

Changing the structural parameters of the PTR and finding its optimum size is an important mean to optimize the refrigeration performance of the PTR. The regenerator is one of the important components of the PTR, where the main losses of the PTR, including pressure drop losses, incomplete heat transfer losses and thermal conductivity losses are mainly concentrated. One of the main goals of optimizing the regenerator is to obtain higher thermodynamic performance by changing the thermal properties of the regenerator packing. Dingli Bao [6] mixed and filled stainless steel wire mesh materials with different mesh numbers in proportion, making the no-load refrigeration temperature of the refrigerator reach 26.7 K. Dang [7] combined experiments and simulations to discuss in detail the potential mechanisms for using mixed screens to significantly reduce regenerator losses. Wenting Wu [8] significantly improved the refrigeration performance at 20 K by coating the surface of the stainless steel wire mesh with the magnetic material Er. The second is to reasonably adjust the structural parameters of the regenerator to find the optimal size. Liu [9] simulated the oscillating flow in the regenerator by establishing a one-dimensional model; according to the numerical simulation results, the performance of regenerators with different geometric sizes was evaluated.

The pulse tube is also one of the important components of the PTR. The flow characteristics and heat transfer process inside the pulse tube are relatively complex. Losses will be caused by the gas working medium crossing the isothermal and adiabatic boundary and periodic heat transfer with the wall surface [10]. Incomplete laminarization leads to the mixing of hot and cold gases, and secondary flow losses can also occur under alternating flow [11]. By analyzing the heat and mass transfer in the pulse tube, Gu [12] found that the streaming is formed due to the generation, evolution and shedding of vortices and pressure drops, which are induced by the hydrodynamic and thermodynamic asymmetries along the refrigeration system. Zhao [13] used a combination of experiments and simulations to investigate the effect of frequency on flow and heat transfer processes in a pulse tube. Ashwin [14] used a numerical simulation method to study the influence of different length-to-diameter ratios of a pulse tube on refrigeration performance. Antao [15] used numerical simulation to study the flow field distribution when the cone angle changes in the conical pulse tube and obtained an effective method to suppress the acoustic streaming. By establishing a three-dimensional CFD model of a pulse tube, the influence of the velocity annular effect, phase angle, mass flow and other factors on pulse tube loss, such as acoustic power loss and loss caused by heat transfer between the gas and wall, were analyzed in detail by Dai [16].

The numerical simulation provides a convenient and accurate research method in comprehending the mechanism of the influence of structural parameters on the cooling effect of PTR, while it also has great limitations in the calculation efficiency, especially for the multi-parameter optimization of PTR. In order to improve the computing efficiency, researchers proposed a method combining a metamodel and optimization algorithm to design cryogenic refrigerators. By fitting the polynomial and response surface between the parameters and the target, the response surface model can better study the interaction of various parameters, predict and optimize the target, and reduce the complex and time-consuming numerical simulation process. It has been applied to the optimization research of PTR by many researchers. Rout [17] used a response surface method to optimize the size of PTR, and combined it with a genetic algorithm to obtain a series of Pareto optimal solutions. Panda [18] adopted the response surface method and particle swarm optimization method, and considered the comprehensive effects of various structure and operating parameters on the Coefficient of Performance (COP) and cooling capacity of a GM refrigerator. The optimization results are 11.39% and 6.43% higher than the cooling power and COP, respectively, of the original numerical model.

Compared with the response surface proxy model, the Kriging metamodel has the advantages of fewer samples and higher model accuracy [19]. The model was first proposed by Krige, a South African engineer, and was first applied to the prediction of gold reserves. Because of its versatility and efficiency, it was gradually widely used in various scientific fields. Zhili Sun [20] proposed an efficient method to deal with high-dimensional models, which simultaneously considered the statistical information provided by the Kriging model and the joint probability density function of input variables. Wang [21] used the two-step Design of Experiments (DOE) approach; using the Kriging model instead of finite element analysis and applying a multi-degree-of-freedom genetic algorithm to solve the optimization problem led to a significant improvement in the model correction problem he studied. Zheng [22] used the Kriging model to optimize the original Alpha Magnetic Spectrometer thermal model and showed that the accuracy of temperature prediction was significantly improved. However, the Kriging model is currently less applied in the field of PTR research, and further research is needed.

This paper proposes a novel method that combines the Kriging metamodel and NSGA-II. The aim is to obtain the optimal size parameters of the regenerator and the pulse tube corresponding to the best overall performance. The CFD model was used to simulate the sample points to provide the cooling temperature and cooling capacity for the Kriging model. NSGA-II is used to obtain the Pareto optimal parameters [23,24]. Accordingly, heat transfer and flow analyses of the optimal cases have been carried out, revealing how the structural parameters affect the performance of the regenerator and pulse tube.

2. Methodologies

2.1. Theoretical Analysis

In the PTR, the working medium gas flows alternately, with fluctuations in pressure, mass flow, temperature and volume flow all varying approximately sinusoidally, so the concept of time averaging needs to be introduced when analyzing the thermodynamic processes. Radebaugh’s enthalpy-flow theory provides an intuitive insight into the refrigeration mechanism of the pulse tube refrigerator.

The acoustic power and enthalpy flow of the working gas are calculated as follows.

$$<P\stackrel{·}{V}>=\frac{1}{\tau }{\displaystyle \underset{0}{\overset{\tau }{\int }}(p-p_0)\stackrel{·}{V}}dt=(1/2)P_d\stackrel{·}{V_d}\cos \theta$$

(1)

$$\begin{array}{l}<\stackrel{·}{H}>=\frac{1}{\tau }{\displaystyle \underset{0}{\overset{\tau }{\int }}\stackrel{·}{m}{c}_p}Tdt=(1/2)c_p\stackrel{·}{m_d}{T}_d\cos \varphi \\ \end{array}$$

(2)

where $P_d$, ${\dot{V}}_d$, ${\dot{m}}_d$, $T_d$ are the amplitudes of pressure wave, volume flow rate, mass flow rate and temperature fluctuation, respectively; $\theta$ is the phase angle between the pressure wave and the volume flow; and $\varphi$ is the phase angle between the mass flow and the temperature.

According to the first law of thermodynamics, the cooling capacity of the cold end heat exchanger is given in Equation (3) below, where $\stackrel{·}{<H_{pt-cold}>}$ is the enthalpy flow of the cold end of the pulse tube, and $\stackrel{·}{<H_{r-cold}>}$ is the enthalpy flow of the cold end of the regenerator. The ratio of the enthalpy flow and PV work at the cold end of the pulse tube is usually used to express the pulse tube expansion efficiency, ${\eta }_{pt}$, calculated as shown in Equation (4).

$$<\stackrel{·}{Q_c}>=<\stackrel{·}{H_{pt-cold}}>-<\stackrel{·}{H_{r-cold}}>$$

(3)

$${\eta }_{pt}=\frac{\stackrel{·}{<H_{pt-cold}>}}{<P{\stackrel{·}{V}}_{cold}>}$$

(4)

2.2. Mathematical Model

For the two-dimensional gas phase region, the governing equation used in the numerical simulation is:

Continuity equation:

$$\frac{\partial {\rho }_f}{\partial t}+\nabla \cdot ({\rho }_f\overrightarrow{v})=0$$

(5)

Momentum conservation equation:

$$\frac{\partial }{\partial t}({\rho }_f\overrightarrow{v})+({\rho }_f\overrightarrow{v}\overrightarrow{v})=-\nabla p+\nabla \cdot (\stackrel{=}{\tau })$$

(6)

Energy conservation equation:

$$\frac{\partial }{\partial t}({\rho }_f{E}_f)+\nabla \cdot [(\overrightarrow{v}{\rho }_f{E}_f+p)]=\nabla \cdot [k_f\nabla T+(\stackrel{=}{\tau }\cdot \overrightarrow{v})]$$

(7)

where ${\rho }_f$, $\overrightarrow{v}$, p, $E_f$, $k_f$ and T are the density, velocity, pressure, total energy, thermal conductivity and temperature of the gas, and $\stackrel{=}{\tau }$ is the stress tensor.

The above control equations are applied to the non-porous media region, including the pulse tube, inertia tube and reservoir. The other components are all porous media regions, and the continuity equations and momentum conservation equations for the porous media region are shown below, taking into account the nature of the solid matrix.

Continuity equation:

$$\frac{\partial \epsilon {\rho }_f}{\partial t}+\nabla \cdot (\epsilon {\rho }_f\stackrel{\rightharpoonup }{v})=0$$

(8)

Momentum conservation equation:

$$\frac{\partial }{\partial t}(\epsilon {\rho }_f\stackrel{\rightharpoonup }{v})+\nabla \cdot (\epsilon {\rho }_f\stackrel{\rightharpoonup }{v}\stackrel{\rightharpoonup }{v})=-\epsilon \nabla p+\nabla \cdot (\epsilon \stackrel{=}{\tau })+S_i$$

(9)

where $\epsilon$ is the porosity of the porous medium, and $S_i$ is the source term in the momentum equation, consisting of two parts, the viscous resistance term and the inertial resistance term. For isotropic porous media, the source term can be expressed as Equation (10).

$$S_i=-(\frac{\mu }{\alpha }\stackrel{\rightharpoonup }{v}+\frac{C_2}{2}{\rho }_f\overrightarrow{v}\left|\overrightarrow{v}\right|)$$

(10)

where $\mu$ is the viscosity, $\alpha$ is the permeability, and $C_2$ is the coefficient of inertia resistance.

The equation for the conservation of energy in the region of the porous medium under the heat balance model is Equation (11).

$$\frac{\partial }{\partial t}[\epsilon {\rho }_f{E}_f+(1-\epsilon ){\rho }_s{E}_s]+\nabla \cdot [\overrightarrow{v}({\rho }_f{E}_f+p)]=\nabla \cdot [k_{eff}\nabla T+(\stackrel{=}{\tau }\cdot \overrightarrow{v})]$$

(11)

where $k_{eff}$ is the thermal conductivity of the porous medium, and ${\rho }_s$ and $E_s$ are the density and total energy of the solid material, respectively.

As shown in Figure 1, a two-dimensional axisymmetric model of a pulse tube refrigerator is developed in this paper. The diameter and length of the regenerator (RD and RL) and pulse tube (PD and PL) are important parameters in this paper. The simulation is judged to have reached steady state when the average circulating temperature at the cold end heat exchanger outlet reaches steady state.

Table 1. Dimension parameters, boundary conditions and area type settings of various parts of the pulse tube refrigerator.

Component	No.	Diameter (mm)	Length (mm)	Wall Material	Boundary Condition
Compressor	A	19.08	7.5	stainless steel	Adiabatic
Transfer tube	B	3.1	101	stainless steel	Adiabatic
Aftercooler	C	8	20	Cu	293 K
Regenerator	D	/	/	stainless steel	Adiabatic
Cold heat exchanger (CHX)	E	6	5.7	Cu	Adiabatic/120 K
Pulse tube	F	/	/	Cu	Adiabatic
Warm heat exchanger (WHX)	G	8	10	Cu	293 K
Inertance tube	H	0.85	684	stainless steel	Adiabatic
Reservoir	I	26	130	stainless steel	Adiabatic

lists the geometry, material composition and boundary conditions for each component. The boundary condition of CHX is set as adiabatic when simulating the cooling temperature and 120 K when simulating the cooling capacity. The time step is set to 7.3529 × 10⁻⁴ s. The number of iterations is 40. After verification of the grid independence, the number of grids for the CFD model is determined to be about 9000, which can ensure good convergence and computational efficiency. The CFD model is solved by the commercial software ANSYS FLUENT 2020.

Using a dynamic grid to simulate the compressor operating process, the piston is modeled as a reciprocating wall having an oscillatory velocity. Using a user-defined function, the piston head motion was modeled as Equation (12).

$$A=A_0\sin (\omega t)$$

(12)

where A₀ = 4.511$\times$10⁻³ m, and $\omega$ = 2πf = 213.62 rad/s. The frequency f is 34 Hz. The charge pressure in the system is set to 3.1 MPa. In the porous media area, the three heat exchangers (C, E and G) are filled with internal wire mesh material set to copper, and the regenerator (D) is filled with internal wire mesh material of 304 stainless steel with a wire mesh number 325; both porous media areas have a porosity ε of 0.69, a permeability α of 1.06 × 10⁻¹⁰ m² and an inertia resistance coefficient $C_2$ of 76,090 m⁻¹ [25]. Running one case by a computer with a 3 GHz processor and 1 TB RAM until steady state takes a week or even more.

2.3. Kriging-Based Meta-Modeling Method

The Kriging model is an optimal, linear, unbiased estimation method based on variational function theory and structural analysis in a finite interval range. The global prediction of the sample space is based on a weighted analysis of the information of known points to predict the information of surrounding unknown points. Replacing the original CFD numerical model with the Kriging surrogate model can greatly improve the efficiency of solving the target response.

In the experimental design, the Kriging model was built using Latin Hypercube Sampling (LHS) with a randomly selected sample set $S={[x_1,x_2,\dots ,x_m]}^T$ in a finite interval, corresponding to an output response set $Y={[y_1,y_2,\dots ,y_m]}^T$.

The Kriging model is an unbiased estimation model that can predict the response of unknown test points by knowing the information of sample points. It can usually be expressed as Equation (13).

$$y={\displaystyle \sum _{j=1}^p{\beta }_j{f}_j}(x)+\delta (x)$$

(13)

where $\beta ={[{\beta }_1,{\beta }_2,\dots ,{\beta }_p]}^T$ is the unknown vector of the regression coefficients. The basis function $f(x)$ is a polynomial about x, providing a global approximation to the global space simulation. The $\delta (x)$ is a random error term that satisfies the following conditions.

$$E(\delta (x))=0$$

(14)

$$Var(\delta (x))={\sigma }^2$$

(15)

$$\mathrm{cov}[\delta (x_i),\delta (x_j)]={\sigma }^{2}R(\theta ,x_i,x_j)$$

(16)

$R(\theta ,x_i,x_j)$ is the spatial correlation model of $x_i,x_j$ sample points with $\theta$ as the parameter, in the form of Equation (17).

$$R(\theta ,x_i,x_j)=\exp (-{{\displaystyle \sum _{p=1}^q{\theta }_p\left|x_{ip}-x_{jp}\right|}}^2)$$

(17)

where $\theta ={[{\theta }_1,{\theta }_2,\dots ,{\theta }_q]}^T$, and the $x_{ip}$ and $x_{jp}$ are the ith and pth elements, respectively, of the sample point. The predicted response value at the unknown point x₀ can be calculated from Equation (18).

$$\widehat{y}(x_0)=f^T(x_0)\cdot \widehat{\beta }+r^T(x_0)R^{-1}(y-F\widehat{\beta })$$

(18)

where the correlation coefficient vector $r(x)$ between the sample point and unknown point x, the regression coefficient $\widehat{\beta }$, and the estimate of variance ${\widehat{\sigma }}^2$ are the following equations, respectively.

$$r(x)={[R(\theta ,x,x_1),R(\theta ,x,x_2),\dots ,R(\theta ,x,x_n)]}^T$$

(19)

$$\widehat{\beta }={(F^T{R}^{-1}F)}^{-1}{F}^T{R}^{-1}y$$

(20)

$${\widehat{\sigma }}^2=\frac{1}{n}{(y-F\widehat{\beta })}^T{R}^{-1}(y-F\widehat{\beta })$$

(21)

where $F={[f_1,f_2,\dots ,f_m]}^T$ is the estimated value of the sample set, and $R$ is the estimated value matrix of $R(\theta ,x_i,x_j)$ at the sample point, as shown in Equation (22). The ${\theta }_p$ can be obtained by the following Equation (23).

$$R=\left[\begin{array}{ccc}R(\theta ,x_1,x_2)& \cdots & R(\theta ,x_1,x_n)\\ ⋮& \ddots & ⋮\\ R(\theta ,x_n,x_1)& \cdots & R(\theta ,x_n,x_n)\end{array}\right]$$

(22)

$${\theta }_p=\underset{\theta }{\max }(-\frac{m}{2}\ln {\widehat{\sigma }}^2-\frac{1}{2}\ln \left|R\right|)$$

(23)

The accuracy of the Kriging model is evaluated by the root mean square error. The calculation formula is as follows.

$$RMSE(f,\widehat{f})=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{(f_i-{\widehat{f}}_i)}^2}}$$

(24)

2.4. Model Updating Method Based on NSGA II

In this paper, we have established the Kriging metamodel of dimensional parameters on the cooling temperature (T) and the cooling capacity (W); multiple established Kriging metamodels collectively act as the target response solver for each iteration of the genetic algorithm. In order to solve the multi-objective optimization problem, the Kriging model optimized by NSGA II based on the Pareto optimal solution is adopted. NSGA-II [24] adopts the elite strategy to expand the sample collection range. Non-dominated sorting is used to compete the parent population with the offspring, and the excellent individuals are retained to the next generation population, which improves the optimization accuracy. At the same time, the method of “crowding distance sorting” is adopted, so that the individual can fully expand to the whole Pareto calculation domain and avoid falling into the local optimal solution.

The crowding algorithm requires that all solutions of the population be arranged in ascending order according to the value of the objective function. For each objective function, first find the boundary solution (the solution corresponding to the maximum and minimum value of the function), and set its distance parameter md value to infinity. The distance parameter values of each solution i are shown in Equation (25).

$$d_m(i)=\frac{f_m(i+1)-f_m(i-1)}{f_m^{\max }-f_m^{\min }}$$

(25)

where $f_m^{\max }$ and $f_m^{\min }$ are the maximum and minimum values of the mth objective function. The congestion value of individual i is the weighting of all the objective function distance parameters d_m. The specific process description of NSGA-II is shown in Figure 2a. The schematic diagram of the optimization iterative process is presented in Figure 2b. The set of Pareto solutions as a result of the iterations is the combined set of optimal values on the multiple Kriging metamodels.

3. Design of Experiments

3.1. Model Validation

In this section, the validity of the developed CFD model is verified by comparison with Cha’s work [25]. The dimensions of the model are 8 mm (RD), 58 mm (RL), 5 mm (PD) and 60 mm (PL) for the regenerator and pulse tube dimensions, which are consistent with the dimensions of Cha’s experimental setup. Figure 3 presents the contrasting results of the CFD simulation with Cha. The cooling trend of the CFD model is approximately the same as that of Cha’s experiment, and the average temperature of the simulated CHX wall is finally stable at 71 K, which is lower than the experimental result of 87 K. Chen [26] compares his numerical simulation results with Cha’s experimental results and explains the reasons for the difference: the model adopts the laminar flow assumption and the heat balance model, without wall thickness, and the working gas is an ideal gas. In addition, the difference between the ideal and real heat transfer in the regenerator and heat exchanger is also an important reason for the difference in results. The reliability of this paper has been verified, as it is highly consistent with Chen’s simulation.

3.2. Kriging Meta Modeling

Analyzing the effect of variations in the dimensional parameters of the regenerator and the pulse tube on the refrigeration cooling performance is the purpose of this paper; therefore, the factors are identified as the diameter and length of the regenerator (RD and RL) and the diameter and length of the pulse tube (PD and PL), and the response as the cooling temperature (T) and the cooling capacity (W) used to judge the refrigeration performance, in which the W is measured when the wall temperature of CHX is 120 K. The diameter and length range from 4 to 10 mm and 50 to 150 mm, respectively. A Design of Experiment (DOE) was carried out using Latin Hypercube Sampling (LHS) [27] with 26 sample points in the determined factor interval, among which 4 sample points were used for validation, and a CFD model was developed for each sample point, and T and W were obtained by numerical simulation. The Kriging metamodel is developed and solved based on Matlab 2020.

Table 2. LHS sample set and numerical simulation results.

No.	RD (mm)	RL (mm)	PD (mm)	PL (mm)	T (K)	W (W)
1	5.50	65.36	6.98	116.01	64.49	2.82
2	4.08	106.02	8.96	105.47	92.23	0.65
3	9.05	118.99	6.52	51.63	70.23	2.35
4	9.91	112.27	4.15	133.26	68.55	1.92
5	8.58	53.86	6.00	140.65	65.21	3.05
6	9.48	56.58	7.29	127.15	77.47	1.92
7	6.00	145.69	4.81	83.92	59.00	2.54
8	4.58	127.48	7.60	123.76	70.75	1.10
9	7.35	132.60	4.95	149.81	51.71	2.57
10	6.82	64.47	8.52	68.17	60.07	3.35
11	6.50	140.83	7.56	61.71	50.91	2.45
12	5.25	76.62	8.41	96.34	68.31	2.01
13	8.41	104.62	5.26	59.60	62.58	2.63
14	7.27	138.95	5.50	74.16	53.30	2.57
15	6.52	81.35	9.81	75.70	67.25	2.06
16	9.17	99.14	6.50	100.50	59.73	2.35
17	8.00	93.50	9.43	85.54	66.85	1.76
18	8.00	70.61	8.00	113.65	68.65	2.15
19	5.03	121.07	4.53	91.93	62.02	2.51
20	4.61	87.54	9.36	135.92	94.23	0.52
21	8.00	58	5.00	60	70.08	3.69
22	9.28	54.73	9.90	51.5	79.28	2.18
23	8.00	101.38	8.64	65.26	59.60	2.54

shows the selected sample points and numerical simulation results.

3.3. Kriging Model Validation

To verify the accuracy of the constructed Kriging model, three sample points were randomly generated as test samples within the interval, using the LHS.

Table 3. Test sample results and errors.

No.		Design Parameters				Response Values
	RD (mm)	RL (mm)	PD (mm)	PL (mm)		Kriging	CFD	Error
1	5.22	62.77	8.6	89.9	T(K)	70.48	71.00	0.74%
					W(W)	2.23	2.17	2.82%
2	8.64	141.56	5.2	122.59	T(K)	54.61	53.33	2.40%
					W(W)	2.08	2.02	2.79%
3	6.8	121.49	4.74	132.96	T(K)	51.61	50.49	2.21%
					W(W)	2.85	3.00	4.85%

shows the results of the CFD simulations of the test samples compared to their predicted values in the Kriging model. The maximum error between the CFD results and the predictions for cooling temperature and cooling capacity is 4.85%, with the remaining error values being within 2.5%. T and W can be accurately predicted by the Kriging model within the selected parameter interval.

4. Results and Discussion

4.1. Analysis of Influencing Factors of T and W

Figure 4a shows the response surface of the regenerator sizes RD and RL with respect to the cooling temperature, where the pulse tube size PD = 6.4 mm and PL = 80 mm. With the increase of RD, the cooling temperature first decreases and then increases. The lowest cooling temperature appears between 6–8 mm in diameter, which changes with the change of RL. With the increase of RL, the cooling temperature also shows a trend of decreasing first and then increasing, and the lowest cooling temperature is reached at RL = 120 mm.

Figure 4b shows the response surface of pulse tube sizes PD and PL with respect to the cooling temperature, where the regenerator size RD = 6.9 mm and RL= 80 mm. With the increase of PD, the cooling temperature first decreases and then increases. There is an optimal PD value to minimize the cooling temperature, and the optimal PD value continues to decrease with the increase of PL. The influence of the PL change on the cooling temperature is greatly affected by PD, showing a completely opposite trend under different PD values. When PD is small, the increase of PL will cause the cooling temperature to drop, while when PD is large, the increase of PL will cause the cooling temperature to rise.

Figure 4c shows the response surface of the regenerator sizes RD and RL with respect to the cooling capacity, where the pulse tube size PD = 6.4 mm and PL = 80 mm. With the increase of RL, the cooling capacity decreases continuously. With the change of RD, the cooling capacity first increases and then decreases. According to the difference of RL, the optimal cooling capacity is between 6–7 mm in diameter.

Figure 4d shows the response surface of the pulse tube sizes PD and PL with respect to the cooling capacity, where the regenerator size RD = 6.9 mm and RL = 80 mm. With the increase of PD, the cooling capacity first increases and then decreases. When PD is different, the influence of PL on the cooling capacity shows different trends. When PD is small, the influence of PL on the cooling capacity is small, and when PD is large, the increase of PL will significantly reduce the cooling capacity.

In Figure 4, PD and PL interact with each other when the dimensional parameters of the pulse tube affect the cooling performance. A larger PD causes a larger dead space at the ends of the pulse tube, deteriorating the cooling performance, which sometimes results in the exact opposite trend of the influence of PL. The working gas in the regenerator has to undergo periodic heat exchange with the packing while transferring as much acoustic power as possible from the hot end to the cold end. A smaller RL or RD is, therefore, not conducive to heat exchange in the regenerator, while a larger RL and RD increase the losses in the regenerator and lead to a deterioration in cooling performance. The four dimensional parameters, therefore, need to be taken into account in order to optimize the optimum size.

4.2. Analysis of NSGA II Optimization Results

NSGA II is used to optimize the two Kriging models constructed for the cooling temperature and the cooling capacity with multiple objectives. If single-objective optimization is used, the two Kriging models need to be optimized separately, and it is not possible to obtain an optimal solution that satisfies both objectives. Equation (26) is the multiple objectives function and constraints for NSGA II, where $T(RD,RL,PD,PL)$ and $W(RD,RL,PD,PL)$ are constructed Kriging models of T and W, respectively.

$$\{\begin{array}{l}Find RD,4mm\le RD\le 10mm\\ \hspace{1em}\hspace{1em}RL,50mm\le RL\le 150mm\\ \hspace{1em}\hspace{1em}PD,4mm\le PD\le 10mm\\ \hspace{1em}\hspace{1em}PL,50mm\le PL\le 150mm\\ Min[T(RD,RL,PD,PL)]\\ Max[W(RD,RL,PD,PL)]\end{array}$$

(26)

NSGA II was run using MATLAB 2020 software, the initial population size was 40, the number of iterations was 200, the probability of variation was set to 0.05, and the crossover probability was 0.8.

Figure 5a,b are the Pareto optimal solution and the worst solution sets of the constructed Kriging model, respectively. The optimal solution is found to optimize the constructed Kriging metamodel in order to determine the optimal cooling performance and the corresponding dimensional parameters of the pulse tube refrigerator. The worst solution is used to better compare the analysis with the optimal solution and to investigate the mechanism of the effect of the excellent dimensions on the physical quantities of the pulse tube refrigerator. The black square points in Figure 5 are two representative solutions selected from each of the two sets of Pareto fronts solved—solution A1 with the lowest cooling temperature, solution A2 with the highest cooling capacity, solution A3 with the lowest cooling capacity, and solution A4 with the highest cooling temperature—and their dimensional parameters and error comparisons are shown in

Table 4. Results and errors of the selected Pareto solutions.

No.		Design Parameters				Response Values
	RD (mm)	RL (mm)	PD (mm)	PL (mm)		Kriging	CFD	Error
A1	7.02	122	6.34	107.92	T(K)	43.73	48.19	4.46
					W(W)	2.65	2.64	0.01
A2	6.8	50	6.6	50	T(K)	65.2	62.77	2.43
					W(W)	4.65	4.86	0.21
A3	4	150	10	129.19	T(K)	106.36	107.23.	0.87
					W(W)	0.06	0.21	0.15
A4	4	114.94	10	150	T(K)	108.2	112.79	4.59
					W(W)	0.184	0.05	0.134

. This table shows a maximum difference of 4.59 K in cooling temperature and 0.21 W in cooling capacity between the predicted and simulated values, indicating that the Pareto solution set is highly accurate. The CFD simulation was used to solve for the initial size of the pulsed tube refrigerator (RD = 8 mm, RL = 58 mm, PD = 5 mm, PL = 60 mm), with T being 70.08 K and W being 3.69 W in the results. After optimization, A1 reduces the minimum cooling temperature by 31.24%, and A2 increases the cooling capacity at 120 K by 31.7%. In Section 5, the four solutions are compared and analyzed in detail from the loss of the regenerator and the pulse tube.

5. Analysis of CFD Simulation Results

From the previous optimization results, it is clear that different combinations of regenerator and pulse tube sizes can lead to significant differences in the refrigeration performance of the PTR. Separate analyses for the regenerator and pulse tube are carried out in the subsequent subsections to evaluate the advantages and disadvantages of the refrigeration performance of each Pareto solution from a loss perspective.

5.1. Analysis of Regenerator Losses

Figure 6a shows the variation of acoustic power in different regenerators along an equidistant cross-section in the axial direction. In the following, the regenerator and pulse tube are divided into 6 radial sections along the axial symmetry axis and numbered 1–6, so as to intuitively observe the change rule of each physical quantity along the axial direction. The change trend of acoustic power in A1 and A2 in the figure is relatively similar. Compared with A3 and A4, the acoustic power in the regenerator inlet section is smaller; the acoustic power in the outlet section is larger; and the loss along the way of acoustic power is 37.49 W and 35.13 W, respectively, while the losses along the way of A3 and A4 are 50.51 W and 56.58 W, respectively. The loss of acoustic power in the regenerator is an important reason for the difference in their refrigeration performance.

According to Equation (1), pressure amplitude is an important factor influencing acoustic power. The change of pressure amplitude can reflect the change trend of pressure drop loss. It can be seen from Figure 6b that the pressure amplitude changes along the axial equidistant section in the regenerator, and the trend is similar to Figure 6a. The working medium gas in the regenerator will exchange heat with the packing under the viscous effect, resulting in resistance loss, causing the pressure drop along the axial direction of the regenerator in the figure. The pressure amplitude of A1 and A2 is significantly smaller than that of A3 and A4, making their acoustic power loss significantly smaller than that of A3 and A4. At the same time, the pressure amplitude drops in A2 and A4 are always smaller than in A1 and A3, respectively, so that the increase in pressure drop loss in the regenerator due to an increase in RL can be clearly observed with little difference in other dimensional parameters.

In order to obtain the desired regenerator cold end pressure ratio and pulse tube enthalpy flow, A3 and A4 with a large pressure drop loss of the regenerator need to provide a larger hot end pressure ratio. Figure 6c shows the distribution of the cross-sectional pressure ratios at the cold and hot ends of the four Case regenerators mentioned above. The presence of pressure drop forces a larger pressure ratio at the hot end to obtain a higher cold end pressure ratio and pulse tube enthalpy flow to improve the cooling performance, which increases the demand on the compressor. Because the same UDF is used to simulate the reciprocating motion of the compressor in this CFD model, there is an upper limit to its performance. Therefore, even if the pressure ratio at the hot end of A3 and A4 rises to around 1.358, the pressure ratio at the cold end still only increases to around 1.025, which is a large gap from A1 and A2, because the ideal enthalpy flow of the pulse tube cannot be obtained, which is an important reason for the poor cooling performance of A3 and A4.

5.2. Analysis of Pulse Tube Losses

Different from the laminar flow in the regenerator, the velocity of the working medium gas in the pulse tube is fast, which is easy to develop into vortex. At the same time, the radial difference of the axial velocity caused by the existence of the boundary layer will form a secondary flow in the pulse tube. The complex flow of these working medium gases will restrict the performance of the refrigerator.

The expansion efficiency is usually used to evaluate the loss in the pulse tube. Figure 7 calculates the expansion efficiency of the pulse tube of 4 Cases, in which the expansion efficiency of A1 and A2 are 0.75 and 0.69, respectively, while the expansion efficiency of A3 and A4 are significantly lower than the former, which are 0.53 and 0.54, respectively, indicating that the loss of working medium gas in the pulse tube is also an important factor causing the difference in refrigeration performance.

Figure 8 shows the distribution of the mass flow rate amplitude in the pulse tube along the axial equidistant section. A1–A4 all showed a trend in which the amplitude of the mass flow rate gradually decreased along the axial direction. At the same time, the mass flow rate amplitude of A3 and A4 is significantly smaller than that of A1 and A2, and the difference is more obvious at the cold end of the pulse tube. The amplitude of A1 is 227.53% and 215.99% higher than that of A3 and A4, respectively, and the amplitude of A2 is 298.3% and 283.17% higher than that of A3 and A4, respectively, which will lead to a significant reduction in the fluid oscillation intensity in the pulse tube, especially at the cold end of the pulse tube, which is not conducive to refrigeration.

The distribution of the temperature field in the vasculature and the change of the gas piston can be observed to further verify the previous inference. Figure 9 shows the temperature contours of the four Cases after reaching steady state. Figure 9a shows the temperature distribution when the compression process ends and the expansion state is about to begin in a cycle, and Figure 9b shows the temperature distribution when the expansion process ends and the compression state is about to begin in a cycle. These two moments have been chosen to give a good reflection of the movement of the gas piston. A1 and A2 can observe the obvious movement process of the gas piston, and the reciprocating motion of the gas piston of A2 is more intense, while the reciprocating motion of the gas piston of A3 and A4 is very weak, which cannot effectively develop its refrigeration performance. This is the same as the variation trend of the amplitude of the mass flow rate in the previous article, which can better confirm the previous conclusion. There is a radial temperature gradient in A1 and A2 near the wall of the pulse tube, which is greater in A2 than in A1, as can be seen from the temperature contours. A3 and A4 have a relatively straight isotherm inside, with no significant radial temperature gradient. The radial temperature gradient in the pulse tube is mainly caused by the viscous dissipation between the boundary layer and the wall of the pulse tube, and because of the presence of the boundary layer, A1 and A2 have a significant radial velocity gradient at the wall. This can be further analyzed by observing the flow field and vortex volumes.

Figure 10 shows the distribution of the flow field and vorticity at four moments in the same cycle of the pulse tube when the four Cases reach steady state. The distribution of the vortex fields in A1 and A2 is similar, with the vortices mainly located close to the wall of the pulse tube and developing axially with the compression and expansion process, mainly due to the radial temperature gradient caused by the viscous effect of the boundary layer in A1 and A2, which leads to the formation of the Rayleigh flow. As the radial temperature gradient is greater in A2, the vortex of A2 occupies more volume in the pulse tube, and its corresponding vortex value is significantly higher than that of A1, thus further confirming that the pulse tube loss is greater in A2 than in A1, resulting in a lower expansion efficiency in A2 than in A1. It can be seen that the vortex fields of A3 and A4 are approximately the same, as the pulse tube diameters of A3 and A4 are both 10 mm, while the diameters of the CHX and WHX connected to them are 6 mm and 8 mm, respectively. This radial difference in size can lead to a larger dead volume at the inlet and outlet ends of the pulse tube, and fully developed vortices are found in both Figure 10c,d in the above-mentioned areas. In addition to the significant gradient in mass flow rate between the cases, the vortex generation conditions are also an important reason for the difference. As can be seen from Figure 7, the dead volume of the pulse tube has a more significant effect on the efficiency than the dissipation at the boundary layer.

6. Conclusions

This paper proposed an approach of combining the Kriging metamodel and genetic algorithms for the optimal design of a pulse tube refrigerator. Conclusions can be drawn as follows:

(I)

The two Kriging models established with four dimensional parameters have been verified to be highly reliable. The results show that the established Kriging model presents a prediction error of about 2.5%. Therefore, replacing the CFD model with it and using it for model updates can significantly improve the iterative efficiency.

(II)

Combining the Kriging method modeling and NSGA II can solve multi-objective optimization problems. In this paper, the Pareto front solution with the optimal and worst cooling temperature and cooling capacity is solved, and a 31.24% drop in the minimum cooling temperature and a 31.7% increase in the cooling capacity at 120 K are achieved after optimization.

(III)

A1–A4 are analyzed from the perspective of the regenerator and the pulse tube, respectively. The pressure drop loss at the regenerator is an important reason for the deterioration of the refrigeration performance, the unreasonable pulse tube structure will cause a large dead space, and the existence of eddy viscosity will seriously restrict the expansion efficiency of the pulse tube.

(IV)

After verification, the Kriging model combined with NSGA-II optimization is highly feasible in the field of pulse tube research, and the current study provides a scientific and efficient design method for miniature cryogenic refrigerators.