Analysis of Flow and Heat Transfer Characteristics and Multi-Objective Optimization for Sinusoidal PCHE

: A Printed Circuit Heat Exchanger (PCHE) is a compact heat exchanger with high temperature and pressure resistance and is considered one of the best choices for the recuperators in the Supercritical Carbon dioxide (S-CO 2 ) Brayton cycle. The ﬂow and heat transfer performance of sinusoidal channel PCHE were analyzed and a second-order regression model was established based on the response surface method to improve the performance of the continuous channel PCHE. It was found that reducing the channel diameter, increasing the channel amplitude, and reducing the channel pitch can increase the average value of the heat transfer coefﬁcient and pressure drop per unit length. Moreover, sensitivity coefﬁcient analysis was used to investigate the inﬂuence of various structural parameters on ﬂow performance, heat transfer performance, and comprehensive performance. In addition, the structure of the sinusoidal channel PCHE was optimized using a multi-objective genetic algorithm, and three sets of Pareto optimal solutions were obtained. The corresponding optimal channel diameter D , channel amplitude A , and channel pitch Lp were in the range of 1.0–1.7 mm, 2.4–3.0 mm, and 15.1–17.0 mm, respectively, which can provide theoretical basis for the design of PCHE.


Introduction
With the increasingly prominent global energy crisis, efficient energy utilization has attracted much attention.The S-CO 2 Brayton cycle has become a hotspot for research in recent years owing to the advantages of high system cycle efficiency, compact structure and wide environmental adaptability [1][2][3][4][5].The recuperators are the most significant heat transfer component in the S-CO 2 Brayton cycle, substantially affecting system cycle efficiency and compactness.It needs to operate under the conditions of high pressure and temperature [4].PCHE, with a high temperature and pressure resistance and compact structure, is considered one of the best candidates for recuperators in the S-CO 2 Brayton cycle [5].
PCHE, as a new type of plate heat exchanger, is manufactured using chemical etching and diffusion welding techniques.Due to these techniques, PCHE can withstand pressures above 90 MPa and temperatures above 900 • C, with a unit volume-to-surface area ratio of up to 2500 m 2 /m 3 and a heat transfer efficiency above 98% [6].Straight and zigzag channels are the first developed and produced structures of PCHE.Chen et al. [7] made two small-sized straight channel PCHEs and experimentally verified their feasibility in ultra-high-temperature reactors.Meshram et al. [8] compared the temperature and pressure drop changes along the high-temperature recuperator and low-temperature recuperator of zigzag channel PCHE and straight channel PCHE by numerical simulation.The results showed that the performance of heat transfer for the zigzag channel was better than that of the straight channel, but the flow performance was worse than that of the straight channel.
Energies 2023, 16, 5763 2 of 16 Finally, friction factor f of straight channel and zigzag channel and empirical formulas for Nusselt number Nu were proposed.Li et al. [9] performed a numerical simulation of zigzag channel PCHE and straight channel PCHE, indicating that the flow resistance of straight channel PCHE was much lower than that of zigzag channel PCHE.
Although zigzag channel PCHE exhibits excellent performance of heat transfer, its pressure drop is higher than other types of flow channel [10][11][12].To reduce the pressure loss, many scholars have proposed non-continuous channels such as S-shaped, airfoil-shaped, and C-shaped channels to reduce the pressure drop of PCHE [13][14][15][16][17][18][19].Although these structures can effectively reduce the pressure drop of PCHE, they are too complicated to manufacture and difficult to put into practical production [20].Some scholars have also optimized the continuous channels.Baik et al. [21] optimized the sharp corners of zigzagtype PCHE into smooth corners with a radius of 1 mm.The results showed the rounded structure can reduce the pressure drop by 40% to 65%.Abed et al. [22] researched the flow and heat transfer in a rectangular cross-section serpentine-channel PCHE by combining experimental and numerical simulation methods.It was affirmed that it enhanced the heat transfer ability and increased the pressure loss compared with the straight channel, thus obtaining improved comprehensive performance.Wen et al. [23] studied the performance of heat transfer and the characteristics of flow for sinusoidal corrugated channel PCHEs and found that increasing the amplitude and reducing the pitch can improve the heat transfer rate.The maximum rate of heat transfer calculated in the corrugated channel model was 33.74% higher than that of the straight channel.Bennett et al. [24] used a non-geometric resolution III experimental design method without replication to study the sensitivity of 20 operating conditions and geometric parameters of zigzag channel PCHE to fluid heat transfer and flow performance.It should be noted that the heat transfer and flow performance of zigzag channel PCHE were more sensitive to the changes in channel bending angle, mass flow rate, curvature radius.Jin et al. [25] selected the pitch, bending angle, and bending radius of the zigzag channel with chamfers as design variables, with heat transfer coefficient and pressure drop as optimization targets.The non-dominated sorting genetic algorithm II was used to optimize the performance of the PCHE.The results showed that both the bending angle and bending radius obviously affected the performance of PCHEs.However, the effect of pitch on the performance of the PCHE was relatively small.Du et al. [26] adopted the channel width, channel depth and aspect ratio of the rectangular section straight channel PCHE as the design variables, and the temperature drop, pressure drop and heat transfer coefficient as the objective function for multi-objective optimization by the method of the numerical simulation.It was concluded that the structure is optimal with a width of 0.29 mm and a depth of 0.39 mm.
As discussed above, many scholars have optimized their structures to improve the heat transfer and flow characteristics of continuous channel PCHEs.However, there are many performance evaluation indicators, and most studies still need to provide the optimal geometric parameters.There needs to be more research on optimizing the structure of sinusoidal PCHEs.The lack of structural parameters for sinusoidal channels makes it difficult to qualitatively and quantitatively analyze and optimize the performance of PCHEs.Therefore, this work analyzes the mechanism of the heat transfer and flow and structural optimization of the sinusoidal PCHE.

Physical Model and Boundary Conditions
The flow channels of PCHE are very small, and in practical engineering, the number of flow channels in a PCHE can reach millions.Modeling and simulating the entire PCHE would require an unacceptable amount of computational resources and time, so simplification of the PCHE is necessary.Since the PCHE is made by welding a metal plate with identical flow channels, its geometric structure is periodic.Therefore, selecting one flow channel as a heat transfer unit for simulation can reflect the overall flow and heat transfer performance of the PCHE. Figure 1 shows the overall schematic diagram of the heat transfer unit of the sinusoidal channel PCHE.The inlet and outlet face of the cold and hot fluids parallel the yoz plane, and the channel extends in the x-axis direction.To ensure uniform flow velocity at the inlet and prevent backflow at the outlet, a length of 100 mm extends the inlet and outlet ends, and the total length of the channel in the x-axis direction is 440 mm.channel as a heat transfer unit for simulation can reflect the overall flow and heat transfer performance of the PCHE.
Figure 1 shows the overall schematic diagram of the heat transfer unit of the sinusoidal channel PCHE.The inlet and outlet face of the cold and hot fluids parallel the yoz plane, and the channel extends in the x-axis direction.To ensure uniform flow velocity at the inlet and prevent backflow at the outlet, a length of 100 mm extends the inlet and outlet ends, and the total length of the channel in the x-axis direction is 440 mm. Figure 2 shows the structural parameters of the channel.A and Lp represent the amplitude and pitch of the channel, respectively, where A is set to 2 mm, and Lp is set to 20 mm.The structure of the sinusoidal channel is periodic in the x-axis, with a length of one pitch for each cycle, and there are a total of 12 cycles.Figure 3 shows a cross-sectional view at location B-B, where the diameter of the cold channel D is 1.5 mm, the unit width W is 2.5 mm, and the unit height H is 3.2 mm.The structural parameters of the cold and hot channels are equal to the distance from the left and right outer walls.Figure 3 shows a cross-sectional view of section B-B.Periodic walls are set on the top and bottom outer walls and the left and right outer walls in the solid domain.The inlet condition is set as mass flow inlet, and the outlet condition is set as pressure outlet.CO2 is chosen as the operating fluid for both hot and cold streams.The working conditions selected for numerical simulation in [8] are as follows: the inlet temperature Tin, c of the cold stream is 500 K, and the outlet pressure Pout, c is 22.5 MPa; the inlet temperature Tin, h of the hot stream is 730 K, and the outlet pressure Pout, h is 9 MPa.The CO2 thermophysical Figure 2 shows the structural parameters of the channel.A and Lp represent the amplitude and pitch of the channel, respectively, where A is set to 2 mm, and Lp is set to 20 mm.The structure of the sinusoidal channel is periodic in the x-axis, with a length of one pitch for each cycle, and there are a total of 12 cycles.Figure 3 shows a cross-sectional view at location B-B, where the diameter of the cold channel D is 1.5 mm, the unit width W is 2.5 mm, and the unit height H is 3.2 mm.The structural parameters of the cold and hot channels are equal to the distance from the left and right outer walls.
channel as a heat transfer unit for simulation can reflect the overall flow and heat transfer performance of the PCHE.
Figure 1 shows the overall schematic diagram of the heat transfer unit of the sinusoidal channel PCHE.The inlet and outlet face of the cold and hot fluids parallel the yoz plane, and the channel extends in the x-axis direction.To ensure uniform flow velocity at the inlet and prevent backflow at the outlet, a length of 100 mm extends the inlet and outlet ends, and the total length of the channel in the x-axis direction is 440 mm. Figure 2 shows the structural parameters of the channel.A and Lp represent the amplitude and pitch of the channel, respectively, where A is set to 2 mm, and Lp is set to 20 mm.The structure of the sinusoidal channel is periodic in the x-axis, with a length of one pitch for each cycle, and there are a total of 12 cycles.Figure 3 shows a cross-sectional view at location B-B, where the diameter of the cold channel D is 1.5 mm, the unit width W is 2.5 mm, and the unit height H is 3.2 mm.The structural parameters of the cold and hot channels are equal to the distance from the left and right outer walls.Figure 3 shows a cross-sectional view of section B-B.Periodic walls are set on the top and bottom outer walls and the left and right outer walls in the solid domain.The inlet condition is set as mass flow inlet, and the outlet condition is set as pressure outlet.CO2 is chosen as the operating fluid for both hot and cold streams.The working conditions selected for numerical simulation in [8] are as follows: the inlet temperature Tin, c of the cold stream is 500 K, and the outlet pressure Pout, c is 22.5 MPa; the inlet temperature Tin, h of the hot stream is 730 K, and the outlet pressure Pout, h is 9 MPa.The CO2 thermophysical  Figure 2 shows the structural parameters of the channel.A an plitude and pitch of the channel, respectively, where A is set to 2 m mm.The structure of the sinusoidal channel is periodic in the x-ax pitch for each cycle, and there are a total of 12 cycles.Figure 3 show at location B-B, where the diameter of the cold channel D is 1.5 m 2.5 mm, and the unit height H is 3.2 mm.The structural paramet channels are equal to the distance from the left and right outer wa  Figure 3 shows a cross-sectional view of section B-B.Periodic and bottom outer walls and the left and right outer walls in the s condition is set as mass flow inlet, and the outlet condition is set a is chosen as the operating fluid for both hot and cold streams.T selected for numerical simulation in [8] are as follows: the inlet t Figure 3 shows a cross-sectional view of section B-B.Periodic walls are set on the top and bottom outer walls and the left and right outer walls in the solid domain.The inlet condition is set as mass flow inlet, and the outlet condition is set as pressure outlet.CO 2 is chosen as the operating fluid for both hot and cold streams.The working conditions selected for numerical simulation in [8] are as follows: the inlet temperature T in, c of the cold stream is 500 K, and the outlet pressure P out, c is 22.5 MPa; the inlet temperature T in, h of the hot stream is 730 K, and the outlet pressure P out, h is 9 MPa.The CO 2 thermophysical parameters are obtained through NIST 9.0 software and fitted using the coefficients listed in Table 1.Inconel 617 alloy is the solid material which can be used under extreme temperature and pressure conditions [8].In the subsequent simulation process, Inconel 617 alloy is treated as a constant material property, with a density of 8360 kg/m 3 , a specific heat capacity at constant pressure of 417 J/(kg•K), and a thermal conductivity of 21 W/(m•K).Based on the operating conditions of the S-CO 2 Brayton cycle high-temperature recuperator [11,12], the Reynold number at the inlet of the cold stream is selected to be 10,000, 20,000, 30,000, 40,000 and 50,000 for simulation.
In addition, the SIMPLE algorithm is used to couple and solve the velocity and pressure terms in the simulation.The PRESTO! format is used for the pressure term and the second-order upwind discretization format is used for the momentum equation, energy equation, turbulent kinetic energy equation, and specific dissipation rate equation.The simulation is considered converged as the residual of the energy equation is less than 10 −6 , the residuals of the other equations are less than 10 −3 , and the pressure drop change at the inlet and outlet of CO 2 is within 0.01%.

Grid Independence Verification and Model Validation
The meshing component in the ANSYS Workbench platform is used with a multiregion approach to generate a hexahedral mesh for the geometric model.For the boundary layer mesh, the first layer thickness is set to 0.005 mm with five layers, and the boundary layer growth rate is set to 1.2.The mesh near the wall needs to satisfy y+ less than 1.The mesh details are shown in Figure 4.
Energies 2023, 16, x FOR PEER REVIEW 4 of 16 parameters are obtained through NIST 9.0 software and fitted using the coefficients listed in Table 1.Inconel 617 alloy is the solid material which can be used under extreme temperature and pressure conditions [8].In the subsequent simulation process, Inconel 617 alloy is treated as a constant material property, with a density of 8360 kg/m 3 , a specific heat capacity at constant pressure of 417 J/(kg⋅K), and a thermal conductivity of 21 W/(m⋅K).Based on the operating conditions of the S-CO2 Brayton cycle high-temperature recuperator [11,12], the Reynold number at the inlet of the cold stream is selected to be 10,000, 20,000, 30,000, 40,000 and 50,000 for simulation.

Work Condition Polynomials
500-730 K, 22.5 MPa 500-730 K, 9 MPa In addition, the SIMPLE algorithm is used to couple and solve the velocity and pressure terms in the simulation.The PRESTO! format is used for the pressure term and the second-order upwind discretization format is used for the momentum equation, energy equation, turbulent kinetic energy equation, and specific dissipation rate equation.The simulation is considered converged as the residual of the energy equation is less than 10 −6 , the residuals of the other equations are less than 10 −3 , and the pressure drop change at the inlet and outlet of CO2 is within 0.01%.

Grid Independence Verification and Model Validation
The meshing component in the ANSYS Workbench platform is used with a multiregion approach to generate a hexahedral mesh for the geometric model.For the boundary layer mesh, the first layer thickness is set to 0.005 mm with five layers, and the boundary layer growth rate is set to 1.2.The mesh near the wall needs to satisfy y+ less than 1.The mesh details are shown in Figure 4.The mesh independence is verified by selecting the results of cold stream heat transfer coefficient h and unit pressure drop Δp as the number of mesh nodes changes.Figure 5 shows the variation in cold stream heat transfer coefficient h and unit pressure drop Δp with increased mesh nodes.As the number of mesh nodes increases from 1,154,175 to 1,520,620, from 1,520,620 to 1,910,230, from 1,910,230 to 2,374,084, and from 2,374,084 to The mesh independence is verified by selecting the results of cold stream heat transfer coefficient h and unit pressure drop ∆p as the number of mesh nodes changes.Figure 5 shows the variation in cold stream heat transfer coefficient h and unit pressure drop ∆p with increased mesh nodes.As the number of mesh nodes increases from 1,154,175 to 1,520,620, from 1,520,620 to 1,910,230, from 1,910,230 to 2,374,084, and from 2,374,084 to 3,129,565, the Energies 2023, 16, 5763 5 of 16 variation rates of heat transfer coefficient are 2.14%, 1.24%, 0.52%, and 0.15%, respectively.The variation rates of unit pressure drop are 1.97%, 1.06%, 0.41%, and 0.18%, respectively.It can be seen that when the number of mesh nodes reaches 2,374,084, the changes in cold stream heat transfer coefficient and unit pressure drop are minimal.Increasing the number of mesh nodes beyond this point will not affect the simulation results.Therefore, the number of mesh nodes is controlled at around 2,374,084.
Energies 2023, 16, x FOR PEER REVIEW 5 of 16 3,129,565, the variation rates of heat transfer coefficient are 2.14%, 1.24%, 0.52%, and 0.15%, respectively.The variation rates of unit pressure drop are 1.97%, 1.06%, 0.41%, and 0.18%, respectively.It can be seen that when the number of mesh nodes reaches 2,374,084, the changes in cold stream heat transfer coefficient and unit pressure drop are minimal.
Increasing the number of mesh nodes beyond this point will not affect the simulation results.Therefore, the number of mesh nodes is controlled at around 2,374,084.To validate the reliability of the numerical simulation, this section conducts a simulation verification of the experiment conducted by Ishizuka et al. [27].The experiment tested the flow and heat transfer of zigzag-type PCHE using S-CO2 as the working fluid.The periodic boundary conditions are imposed on the top, bottom, left and right walls.The experimental results and the numerical simulation results are shown in Table 2.The maximum error between the numerical simulation results and the experimental results for pressure drop is 9.8%, and for temperature difference is 3.4%.The error between the numerical simulation results and the experimental results is within 10%, proving the method's reliability.Table 3 shows the validation of the simulation conducted by Meshran et al. [8].The maximum error between the numerical simulation results and the literature for pressure drop is 1.9%, and for temperature difference is 1.3%.The reliability of the method is further proved.To validate the reliability of the numerical simulation, this section conducts a simulation verification of the experiment conducted by Ishizuka et al. [27].The experiment tested the flow and heat transfer of zigzag-type PCHE using S-CO 2 as the working fluid.The periodic boundary conditions are imposed on the top, bottom, left and right walls.The experimental results and the numerical simulation results are shown in Table 2.The maximum error between the numerical simulation results and the experimental results for pressure drop is 9.8%, and for temperature difference is 3.4%.The error between the numerical simulation results and the experimental results is within 10%, proving the method's reliability.Table 3 shows the validation of the simulation conducted by Meshran et al. [8].The maximum error between the numerical simulation results and the literature for pressure drop is 1.9%, and for temperature difference is 1.3%.The reliability of the method is further proved.

Governing Equations and Parameter Definition
When using ANSYS Fluent 2021 software for solving, three fundamental equations must be satisfied: Continuity equation (mass conservation equation): Conservation of momentum equation: Conservation of energy equation: where k f is the thermal conductivity of fluid; kt is the turbulent thermal conductivity; δ ij is the Kronecker delta; ν is the kinematic viscosity; µ is the dynamic viscosity.Palko [28] has shown that the SST k-ω model can accurately compute flow field information of supercritical fluids when a sufficient number of mesh nodes are provided, and this turbulence model has been widely used in various simulation calculations.
The equations for the turbulent kinetic energy and specific dissipation rate in the SST k-ω model are given by: where G k and G w are the generation terms for turbulent kinetic energy and specific dissipation rate, respectively; Y k and Y w are the dissipation terms for turbulent kinetic energy and specific dissipation rate, respectively; D w is the cross-diffusion term; S k and S w are the source terms defined according to the actual problem; Γ k and Γ w are the effective diffusion coefficients.The effective diffusion coefficients Γ k and Γ w are expressed by the following formula: where µ t is the turbulent viscosity; σ k and σ ω are constants in the turbulence model.The formula for hydraulic diameter D h calculation is: where S is the cross-sectional area of the channel, in m 2 ; L h is the perimeter of the crosssectional area of the channel, in m.The formula for calculating Reynolds number Re of the fluid is: where u is the inlet velocity of the fluid, in m/s; µ is the dynamic viscosity of the fluid, in kg/(m•s).
The formula for calculating the logarithmic mean temperature difference ∆T m is: where T w is the average temperature of the wall, in K; T in is the inlet temperature of the fluid, in K; T out is the outlet temperature of the fluid, in K.
The formula for calculating the heat transfer coefficient h is: where q is the heat flux density, in W/(m 2 •K).
The formula for calculating the Nusselt number Nu is: The formula for calculating the friction factor f is: where ∆p is the pressure drop between the inlet and outlet of the fluid, in Pa.
The formula for calculating the overall performance parameter PEC [29] is: (14)

Effects of Structural Parameters on the Performance of PCHE
In this section, when evaluating the influence of channel diameter on the performance of the PCHE, the channel amplitude is kept at 2 mm, and the channel pitch is 20 mm.When evaluating the effect of channel amplitude on the performance of the PCHE, the channel diameter is kept at 1.5 mm, and the channel pitch is 20 mm.When evaluating the impact of channel pitch on the performance of the PCHE, the channel diameter is kept at 1.5 mm, and the channel amplitude is 2 mm.
As shown in Figures 6 and 7, with the increase in the Reynolds number, the influence of structural parameters on heat transfer coefficient and pressure drop gradually increases.Reducing the channel diameter, increasing the channel amplitude and decreasing the channel pitch of sinusoidal PCHE can increase the average value of heat transfer coefficient by up to 55.04%, 28.25%, and 16.56%, respectively, and increase the average value of unit pressure drop by up to 87.38%, 30.36%, and 35.08%, respectively.That indicates that the fluid in the channel with a smaller diameter has a higher velocity, higher turbulence kinetic energy, and more substantial heat transport capacity at the same inlet Reynolds number, which enhances the heat transfer performance but weakens the flow performance.The position sketch of the monitoring surface is shown in Figure 8. Plane 1, Plane 2 and Plane 3 are located at the fifth pitch, Plane 4 and Plane 5 are located at the sixth pitch.Axial velocity distribution under different channel amplitude in Figure 9.The larger amplitude of the channel is, the greater offset of the fluid velocity at the corner and the higher the degree of fluid disturbance is.Meanwhile, the flow separation phenomenon become more obvious, and the boundary layer is more significantly destroyed, which improves the heat transfer performance but weakens the flow performance.Axial velocity distribution under different channel pitch in Figure 10.For the channel with a smaller pitch, the frequency of fluid velocity offset is higher, the vortexes in the flow are more frequently destroyed and recombined, and the heat transfer rate is faster.However, the energy loss is more severe, which improves the heat transfer performance but weakens the flow performance.distribution under different channel pitch in Figure 10.For the channel with a smaller pitch, the frequency of fluid velocity offset is higher, the vortexes in the flow are more frequently destroyed and recombined, and the heat transfer rate is faster.However, the energy loss is more severe, which improves the heat transfer performance but weakens the flow performance.
(a) (b) (c)   distribution under different channel pitch in Figure 10.For the channel with a smaller pitch, the frequency of fluid velocity offset is higher, the vortexes in the flow are more frequently destroyed and recombined, and the heat transfer rate is faster.However, the energy loss is more severe, which improves the heat transfer performance but weakens the flow performance.
(a) (b) (c)     As analyzed above, channels of different diameters have different velocities, which cause changes in flow and heat transfer performance.In order to understand the changes in flow and heat transfer performance of fluids under channels of different wave amplitudes and pitches more intuitively, the field synergy principle is used to explain the heat transfer improvement mechanism.According to the field synergy principle, the smaller the synergy angle, the better the synergy between velocity and temperature [30].Fortynine equidistant cross-sections are created along the direction of the cold fluid flow, with the first and last cross-sections corresponding to the inlet and outlet cross-sections of the cold fluid, respectively.The area-weighted average velocity and temperature synergy angle are used for calculation, with its formula given in Equation ( 15).(15) where  is the velocity vector, in m/s;  is the temperature gradient vector, in K.
As shown in Figure 11, the fluid velocity varies periodically along the flow channel, resulting in periodic changes in the synergy angle, which varies in a "W" shape.The smaller synergy angle is located at the peak or trough of the channel.Compared with As analyzed above, channels of different diameters have different velocities, which cause changes in flow and heat transfer performance.In order to understand the changes in flow and heat transfer performance of fluids under channels of different wave amplitudes and pitches more intuitively, the field synergy principle is used to explain the heat transfer improvement mechanism.According to the field synergy principle, the smaller the synergy angle, the better the synergy between velocity and temperature [30].Fortynine equidistant cross-sections are created along the direction of the cold fluid flow, with the first and last cross-sections corresponding to the inlet and outlet cross-sections of the cold fluid, respectively.The area-weighted average velocity and temperature synergy angle are used for calculation, with its formula given in Equation ( 15).(15) where  is the velocity vector, in m/s;  is the temperature gradient vector, in K.
As shown in Figure 11, the fluid velocity varies periodically along the flow channel, resulting in periodic changes in the synergy angle, which varies in a "W" shape.The smaller synergy angle is located at the peak or trough of the channel.Compared with As analyzed above, channels of different diameters have different velocities, which cause changes in flow and heat transfer performance.In order to understand the changes in flow and heat transfer performance of fluids under channels of different wave amplitudes and pitches more intuitively, the field synergy principle is used to explain the heat transfer improvement mechanism.According to the field synergy principle, the smaller the synergy angle, the better the synergy between velocity and temperature [30].Forty-nine equidistant cross-sections are created along the direction of the cold fluid flow, with the first and last cross-sections corresponding to the inlet and outlet cross-sections of the cold fluid, respectively.The area-weighted average velocity and temperature synergy angle are used for calculation, with its formula given in Equation (15).
where U is the velocity vector, in m/s; ∇T is the temperature gradient vector, in K.As shown in Figure 11, the fluid velocity varies periodically along the flow channel, resulting in periodic changes in the synergy angle, which varies in a "W" shape.The smaller synergy angle is located at the peak or trough of the channel.Compared with channels with wave amplitudes of 1 mm and 2 mm, the velocity and temperature synergy angle in the channel with a wave amplitude of 3 mm decreased by 8.08% and 5.02%, respectively.Compared with channels with a pitch of 20 mm and 25 mm, the velocity and temperature synergy angle in the channel with a pitch of 15 mm decreased by 3.51% and 1.95%, respectively.
Energies 2023, 16, x FOR PEER REVIEW 10 of 16 channels with wave amplitudes of 1 mm and 2 mm, the velocity and temperature synergy angle in the channel with a wave amplitude of 3 mm decreased by 8.08% and 5.02%, respectively.Compared with channels with a pitch of 20 mm and 25 mm, the velocity and temperature synergy angle in the channel with a pitch of 15 mm decreased by 3.51% and 1.95%, respectively. (

Sensitivity Analysis
The central composite face-centered design method in response surface analysis is used to design experiments for a three-factor, three-level scheme.The design results in a total of 15 experimental points.Each factor has three levels, with coded levels of −1, 0, and 1, respectively.The coding scheme and level values are shown in Table 4.Under the operating condition of the inlet Reynolds number of 30,000 for the hot and cold fluids, the ANSYS Fluent 2021 was used to simulate the 15 design points, and the corresponding objective function values were obtained.The specific experimental design scheme and results are shown in Table 5.

Sensitivity Analysis
The central composite face-centered design method in response surface analysis is used to design experiments for a three-factor, three-level scheme.The design results in a total of 15 experimental points.Each factor has three levels, with coded levels of −1, 0, and 1, respectively.The coding scheme and level values are shown in Table 4.Under the operating condition of the inlet Reynolds number of 30,000 for the hot and cold fluids, the ANSYS Fluent 2021 was used to simulate the 15 design points, and the corresponding objective function values were obtained.The specific experimental design scheme and results are shown in Table 5.
Based on Table 4, the quadratic regression models for h, ∆p, and PEC are expressed in Equations ( 16)-( 18  In order to study the sensitivity of the three response factors (D, A, Lp) to the objective function, the value levels of the three response factors are uniformly ratioed, as shown in Equation (19).The ratioed response variable values can correspond to their original numerical values, and the specific correspondence can be referred to in Table 6.

Ratio y =
x − x 0 x 0 (19) where x represents the numerical value of the response factor, and x 0 represents the value of the response factor at the 0 level.To study the sensitivity of each of the three objective functions to the response factors, the ratioed D, A, and Lp are used as the horizontal axis, and the three objective functions (h, ∆p, and PEC) are used as the vertical axis, as shown in Figure 12.When studying the sensitivity of the objective function to a certain response factor, the other response factors should be kept at the 0 level, that is, D = 1.5 mm, A = 2 mm, Lp = 20 mm.The sensitivity of the response factor to the objective function depends on the slope of the line.If the slope is positive, the response factor is positively correlated with the objective function.If the slope is negative, the response factor is negatively correlated with the objective function.The larger the absolute value of the slope of the line, the higher the sensitivity of the response factor to the objective function.To quantify the sensitivity of each response factor, the sensitivity coefficient is used as a measure of its size [31].SCy is used to represent the sensitivity coefficient of the objective function y (h, Δp, PEC) to the response factor x (D, A, Lp), and the calculation formula is shown in Equation (20).(20) where y0 represents the objective function value when encoded as 0, and x0 represents the design variable value when encoded as 0.
The results are shown in Table 7.The positive or negative sensitivity coefficient represents the correlation between the objective function and the response factor, positive for a positive correlation and unfavorable for a negative correlation, and the absolute value indicates the sensitivity level.For h, the sensitivity to D is the highest, 3.1 times and 2.8 times that of A and Lp, respectively.For Δp, the sensitivity to D is also the highest, 7.3 times and 3.1 times that of A and Lp, respectively.For PEC, the sensitivity to A is the highest, 7.5 times and 1.2 times that of D and Lp, respectively.From Figure 12, it can be seen that h has the highest sensitivity to D and the lowest sensitivity to Lp.In addition, D and Lp are negatively correlated with the heat transfer coefficient, while A and h are positively correlated.The relationships between h and A, Lp are linear, and the absolute value of the slope of the curve between h and D decreases with the increase of D, indicating that the sensitivity of h to D decreases with the increase of D. The sensitivity of ∆p to D is significantly higher than that to A and Lp.Moreover, D and Lp are negatively correlated with ∆p, while A is positively correlated.The sensitivity of ∆p to D decreases with the increase of D. PEC has the highest sensitivity to A and is positively correlated, while it is negatively correlated with D and Lp.The slope between PEC and D is close to 0, indicating that the influence of D on PEC is small.The absolute value of the slope of the curve between PEC and Lp increases with the increase of Lp, indicating that the sensitivity of PEC to Lp increases with the increase of Lp.
To quantify the sensitivity of each response factor, the sensitivity coefficient is used as a measure of its size [31].SC y is used to represent the sensitivity coefficient of the objective function y (h, ∆p, PEC) to the response factor x (D, A, Lp), and the calculation formula is shown in Equation (20).

SC y =
y − y 0 /y 0 x − x 0 /x 0 (20) where y 0 represents the objective function value when encoded as 0, and x 0 represents the design variable value when encoded as 0.
The results are shown in Table 7.The positive or negative sensitivity coefficient represents the correlation between the objective function and the response factor, positive for a positive correlation and unfavorable for a negative correlation, and the absolute value indicates the sensitivity level.For h, the sensitivity to D is the highest, 3.1 times and 2.8 times that of A and Lp, respectively.For ∆p, the sensitivity to D is also the highest, 7.3 times and 3.1 times that of A and Lp, respectively.For PEC, the sensitivity to A is the highest, 7.5 times and 1.2 times that of D and Lp, respectively.

Multi-Objective Optimization
Based on the MATLAB 2018 platform, the non-dominated sorting genetic algorithm-II with elite strategy is used to optimize and solve the three objective functions (h, ∆p, and PEC).The algorithm parameters are set as follows: the population size is 150, the maximum evolution generation is 300, the crossover probability is 80%, and the mutation probability is 5%.
Figure 13 shows the three-dimensional Pareto optimal solution diagram of the heat transfer performance, flow performance, and comprehensive performance of the sinusoidal PCHE.All points on the Pareto optimal boundary are the optimal solutions of multiobjective functions.This means that there is no second solution that can improve one function's expected value without reducing the expected values of other functions.In other words, no solution can replace the solution on the Pareto optimal boundary.

Multi-Objective Optimization
Based on the MATLAB 2018 platform, the non-dominated sorting genetic a II with elite strategy is used to optimize and solve the three objective functions ( PEC).The algorithm parameters are set as follows: the population size is 150, mum evolution generation is 300, the crossover probability is 80%, and the muta ability is 5%.
Figure 13 shows the three-dimensional Pareto optimal solution diagram o transfer performance, flow performance, and comprehensive performance of th dal PCHE.All points on the Pareto optimal boundary are the optimal solutions objective functions.This means that there is no second solution that can improve tion's expected value without reducing the expected values of other functions words, no solution can replace the solution on the Pareto optimal boundary.The black dots in the Figure 13 represent the projection of the Pareto optimal boundary on each coordinate plane, as shown in Figure 14. Figure 14a represents the Pareto optimal boundary of h and ∆p.It can be seen that point C is the Pareto optimal solution that balances the two objective functions of h and ∆p.

Conclusions
The sinusoidal channel PCHEs have the potential as a high-temperature recuperator of the S-CO 2 recompression Brayton cycle.In this paper, a three-dimensional numerical simulation was conducted to investigate the thermo-hydraulic performance of the sinusoidal channel PCHEs.Response surface analysis and multi-objective optimization were conducted on the structural parameters.The main conclusions are as follows: (1) Reducing the channel diameter, increasing the channel amplitude, and reducing the channel pitch can increase the average value of the heat transfer coefficient by a maximum of 55.04%, 28.25%, and 16.56%, respectively.It can increase the average value of the pressure drop per unit length by a maximum of 87.38%, 30.36%, and 35.08%, respectively.(2) The sensitivity of h to D is the highest, with a sensitivity coefficient equivalent to 3.1 times and 2.8 times that of A and Lp, respectively.The sensitivity of ∆p to D is the highest, with a sensitivity coefficient equivalent to 7.3 times and 3.1 times that of A and Lp, respectively.The sensitivity of PEC to A is the highest, with a sensitivity coefficient equivalent to 7.5 times and 1.2 times that of D and Lp, respectively.(3) Based on the multi-objective genetic algorithm, a three-dimensional Pareto optimal solution set was obtained with h, ∆p, and PEC as the optimization objectives for the sinusoidal PCHE.From the perspective of balancing each objective function, three sets of Pareto optimal solutions were obtained, corresponding to the optimal h, ∆p, and PEC intervals of 5574.

Figure 10 .
Figure 10.Axial velocity distribution under different channel pitch.

Figure 11 .
Figure 11.Synergy angle for the PCHE with different wave amplitudes and pitches.(a) Different amplitude; (b) different pitch.

Figure 11 .
Figure 11.Synergy angle for the PCHE with different wave amplitudes and pitches.(a) Different amplitude; (b) different pitch.

Figure 12 .
Figure 12.Sensitivity analysis of h, Δp, and PEC to each response factor in sinusoidal PCHE.From Figure12, it can be seen that h has the highest sensitivity to D and the lowest sensitivity to Lp.In addition, D and Lp are negatively correlated with the heat transfer coefficient, while A and h are positively correlated.The relationships between h and A, Lp are linear, and the absolute value of the slope of the curve between h and D decreases with the increase of D, indicating that the sensitivity of h to D decreases with the increase of D. The sensitivity of Δp to D is significantly higher than that to A and Lp.Moreover, D and Lp are negatively correlated with Δp, while A is positively correlated.The sensitivity of Δp to D decreases with the increase of D. PEC has the highest sensitivity to A and is positively correlated, while it is negatively correlated with D and Lp.The slope between PEC and D is close to 0, indicating that the influence of D on PEC is small.The absolute value of the slope of the curve between PEC and Lp increases with the increase of Lp, indicating that the sensitivity of PEC to Lp increases with the increase of Lp.To quantify the sensitivity of each response factor, the sensitivity coefficient is used as a measure of its size[31].SCy is used to represent the sensitivity coefficient of the objective function y (h, Δp, PEC) to the response factor x (D, A, Lp), and the calculation formula is shown in Equation(20).

Figure 12 .
Figure 12.Sensitivity analysis of h, ∆p, and PEC to each response factor in sinusoidal PCHE.

Figure 13 .
Figure 13.Pareto-optimal solution set of three objectives in sinusoidal PCHE displayed dimensional plot.
Figure 14b represents t optimal boundary of Δp and PEC.It can be seen that point B is the Pareto optima that balances the two objectives of Δp and PEC.
Figure 14c represents the Paret boundary of h and PEC.It can be seen that point D is the ideal point for both h where both expected objectives are achieved optimally.Point A is the non-idea both h and PEC, where both expected objectives are the worst.The objective fun ues and design variable values for points A, B, C, and D are shown in Table8.

Figure 13 .
Figure 13.Pareto-optimal solution set of three objectives in sinusoidal PCHE displayed in a threedimensional plot.
Figure 14b represents the Pareto optimal boundary of ∆p and PEC.It can be seen that point B is the Pareto optimal solution that balances the two objectives of ∆p and PEC.
Figure 14c represents the Pareto optimal boundary of h and PEC.It can be seen that point D is the ideal point for both h and PEC, where both expected objectives are achieved optimally.Point A is the non-ideal point for both h and PEC, where both expected objectives are the worst.The objective function values and design variable values for points A, B, C, and D are shown in Table8.

Figure 14 .
Figure 14.Pareto-optimal solution set displayed in a two-dimensional plot.(a) The Pareto-optimal frontier of h and ∆p; (b) the Pareto-optimal frontier of ∆p and PEC; (c) the Pareto-optimal frontier of h and PEC.

Table 1 .
CO 2 physical properties polynomial fitting results.

Table 2 .
Comparison between numerical calculation and experimental data.

Table 3 .
Validation results of the simulation by Meshran.

Table 2 .
Comparison between numerical calculation and experimental data.

Table 3 .
Validation results of the simulation by Meshran.

Table 4 .
Coding scheme and value levels table.

Table 4 .
Coding scheme and value levels table.

Table 6 .
Response factors and ratioed response variable values.

Table 7 .
Sensitivity coefficients of h, ∆p, and PEC in sinusoidal PCHE.