A Numerical Study on Blade Design and Optimization of a Helium Expander for a Hydrogen Liquefaction Plant

: A design process for cryogenic expanders that supplies 0.5 TPD of liqueﬁed hydrogen in hydrogen liquefaction plants is introduced. To improve the efﬁciency of the expander, the optimum design was conducted by adjusting two rotor shape parameters. The designed expander for hydrogen liquefaction has a target rotation speed of 75,000 rpm, and helium is applied as the working ﬂuid. Since the operating temperature of the expander is as low as 49 K, a design that reﬂects the real gas properties must be considered. For a high-efﬁciency hydrogen liquefaction plant, increasing the expander efﬁciency is one of the most critical issues. In this study, the efﬁciency of the cryogenic expander was optimized using the response surface method (RSM). The hub and shroud meridional contours and blade β angle distributions were chosen as the design parameters. As a result, through the optimized design, it was possible to improve the expander efﬁciency by up to 1.98% compared to the original expander. Flow analysis was conducted to investigate the reason for the efﬁciency improvement. Through this study, the effect of the blade meridional contour and blade β angle on the cryogenic expander efﬁciency was veriﬁed.


Introduction
A new climate regime, POST-2020, was agreed upon at the climate change conference held in Paris in December 2015. In order to implement this new climate system, all countries must strive to achieve their greenhouse gas reduction targets. To this end, many efforts are being made to develop clean energy worldwide. Representatively, there have been studies on the production and use of hydrogen energy. Recently, a full-scale study on the development of core technologies for hydrogen liquefaction plants began in South Korea.
A hydrogen liquefaction plant is a closed-loop system, and it consists of a compressor, heat exchanger, cold box, and expander. The expander is one of the core components in a cryogenic system. In the cryogenic expander, the working fluid is expanded under extremely low temperatures, e.g., below 93 K. Therefore, some of the parts, i.e., the rotor, shaft, and labyrinth seal, are operated under cryogenic conditions. In general, the cryogenic expander rotates at a high speed, from tens of thousands to hundreds of thousands of revolutions per minute, so rotordynamic stability is important.
The technologies required to develop the cryogenic expander are as follows: aerodynamic design technology with an expander loss model for cryogenic working fluids, seal and insulation technology, thermal stress analysis technology, high-speed oil-free bearing technology, clearance prediction and control technology, secondary flow stabilization technology, and expander power control technology. However, because of the information security of the industry, the introduction of these technologies in the literature or at conferences has been limited. However, it is possible to identify the latest technology trends and technical issues through the publications of a few companies that produce cryogenic expanders [1][2][3].
There are several methods to exhaust the output power generated by the expander. First, a compressor is installed on the opposite side of the expander with the same shaft. It can be applied to compress the working fluid at the cryogenic system [4]. This expandercompressor assembly is called a compander. When it is used as a compander, there is an advantage in that the high pressure required in the cryogenic system can be acquired from the compressor. However, there is a disadvantage in that the system operation is complicated.
Second, the power generated from the expander drives the brake compressor to exhaust energy [5][6][7][8]. Alternatively, there is a method of converting the power generated by the expander into heat energy using an eddy current breaker. It has an advantage in that the size of the system can be reduced and the control can be made simple. However, it dissipates the energy generated by the expander; therefore, it is not preferable for a large-capacity plant. Third, the expander rotates the generator, and electric power can be acquired from the generator. If the capacity of the expander is large enough, the application of a generator can be considered.
As described above, the method of utilizing expander energy can be determined by the characteristics of the cryogenic system. In this research, the brake compressor is connected to the expander so that the speed of the expander can be controlled by adjusting the compressor's throttle valve.
Song et al. conducted blade optimization to reduce the swirl flow and cavitation in the draft tube of a cryogenic liquid turbine [9]. Huang et al. investigated the effect of the blade loading distribution and blade lean angles on the cavitation in the draft tube of a cryogenic expander [10]. However, there has been little research on the optimization of cryogenic expander efficiency.
To develop a highly efficient cryogenic system, it is necessary to increase the efficiency of the expander. In this paper, the optimization of an expander rotor was conducted to improve expander efficiency. The expander rotor rotates at tens of thousands of revolutions per minute under cryogenic conditions. However, there is no sufficient design modeling for an expander with cryogenic fluids. Therefore, careful decisions must be made on the shape factors that determine the efficiency of the expander.
In this study, the rotor blade meridional shape and blade β angle were optimized to improve the expander efficiency. This can be used as a reference to design a high-efficiency cryogenic expander in the future.

Operation Conditions
In this study, cycle analysis was conducted to develop a high-efficiency hydrogen liquefaction system. Through the cycle analysis, the operating conditions of the component were determined. The target capacity of the hydrogen liquefaction system is condensing 0.5 tons of hydrogen per day. Figure 1a shows a schematic of a 0.5-TPD class hydrogen liquefaction plant. The helium working fluid absorbs the cold heat of the vaporized nitrogen through a heat exchanger (Hx1), and by flowing through a sequence of heat exchangers and expanders, the temperature is lowered so that hydrogen can be liquefied. Hydrogen in the gaseous state is liquefied through heat exchangers Hx1−Hx4 and stored in a tank.
As shown in Figure 1a, the system consists of two expanders, one compressor, and four heat exchangers. The expanders are connected in series. In this study, design and optimization were conducted for one of the expanders. Figure 1b illustrates the temperature-entropy (T-s) diagram of the hydrogen liquefaction plant. Table 1 shows the components in each process in the T-s diagram. A heat exchanger (process 4 -5 in Figure 1) is installed between Expander 1 and Expander 2. The expander that was studied in this research was Expander 1, which relates to processes 3 -4 in the system. Figure 1b illustrates the temperature-entropy (T-s) diagram of the hydrogen liquefaction plant. Table 1 shows the components in each process in the T-s diagram. A heat exchanger (process ④-⑤ in Figure 1) is installed between Expander 1 and Expander 2. The expander that was studied in this research was Expander 1, which relates to processes ③-④ in the system.

Process @ T-s Diagram Component
Heat exchanger 3 ⑤-⑥ Expander 2   Table 2 presents the operation conditions of the expanders. The considered expansion ratios of Expander 1 and Expander 2 were 1.6 and 3.69, respectively.

Preliminary Design
A preliminary design was performed to determine the rotational speed and diameter of the expander rotor. In the preliminary design, specific velocity and specific diameter analyses were conducted. The definitions of the specific velocity and specific diameter are as follows: Specific diameter : D s = d∆h 0.25 isen Q 6 (2) Table 3 shows the preliminary design results for Expander 1. In the case of a centrifugal expander, the specific speed is generally selected between 0.5 and 0.6 in consideration of expander efficiency. However, in this study, the specific speed of the expander was selected as 0.392, considering rotordynamic performance. The design speed was 75,000 rpm. Through specific diameter analysis, it was estimated that the diameter of the rotor should be 53.77 mm.
In the case of the radial inflow expander, the relation of the specific speed N s and the specific diameter D s was N s × D s = 2 [11]. Therefore, the specific diameter D s was determined when the specific velocity N s was chosen. For this reason, U_tip in Table 3 is the same against the specific velocity N s .

Meanline Design
Meanline design, which is also called one-dimensional design, was conducted based on the results of the preliminary design. Meanline design is a flow path design including the blade shape and flow angle. In this study, meanline design was conducted using ConceptsNREC Inc's RITAL program. The Rodgers loss model was applied for the nozzle design, and the CETI passage loss model was applied for the rotor design. Figure 2 shows the meanline shape of the expander. The designed expander consists of a volute, nozzle, and rotor. Generally, a volute is symmetric. However, in this research, an overhang-shaped volute was selected to shorten the shaft length and increase the assembly. Cho et al. compared losses according to the shape of volute [12]. They conducted a CFD analysis of the loss in a circular volute and a rounded square volute. In this study, the volute shape was determined based on the fact that a circular volute has less loss than that of a rounded square volute. The major dimensions of the expander are shown in Table 4. The designed tip clearance between the rotor and shroud was 0.3 mm.
using ConceptsNREC Inc's RITAL program. The Rodgers loss model was applie nozzle design, and the CETI passage loss model was applied for the rotor design Figure 2 shows the meanline shape of the expander. The designed expan sists of a volute, nozzle, and rotor. Generally, a volute is symmetric. Howeve research, an overhang-shaped volute was selected to shorten the shaft length crease the assembly. Cho et al. compared losses according to the shape of vol They conducted a CFD analysis of the loss in a circular volute and a rounded volute. In this study, the volute shape was determined based on the fact that a volute has less loss than that of a rounded square volute. The major dimension expander are shown in Table 4. The designed tip clearance between the rotor and was 0.3 mm.
The velocity triangles at the inlet and outlet of the rotor are shown in Fig  this study, the incidence angle was −26.74°, and the deviation angle was 22.8°. shows the performance curves derived from the meanline design. They were p 80% (60,000 rpm), 90% (67,500 rpm), and 100% (75,000 rpm) of the design speed.    The velocity triangles at the inlet and outlet of the rotor are shown in Figure 3. In this study, the incidence angle was −26.74 • , and the deviation angle was 22.8 • . Figure 4 shows the performance curves derived from the meanline design. They were plotted at 80% (60,000 rpm), 90% (67,500 rpm), and 100% (75,000 rpm) of the design speed.    Figure 4a shows the relation between the expansion ratio and mass flow rate. Figure 4b shows the relation between the expansion ratio and efficiency. The efficiency from the meanline design was recalculated by reflecting the loss in the expander. It is defined as the ratio of the enthalpy difference and isentropic enthalpy difference between the inlet and outlet of the expander. It was verified that the target efficiency was satisfied at the design expansion ratio. Figure 4c shows the relation between the expansion ratio and the power. Like the efficiency curve, the power from the meanline design was recalculated by reflecting the loss in the expander. Finally, Figure 4d shows the relation between the expansion ratio and the expander outlet temperature. The cryogenic expander must satisfy the expander outlet temperature in the system requirements. As shown in the figure, the estimated expander outlet temperature meets the design temperature.

3D Shape Generation
Based on the results of the meanline design, a three-dimensional expander geometry was generated to conduct the numerical analysis. The three-dimensional geometry was generated using AxCent program of ConceptsNREC. Figure 5 shows the 3D shapes of the nozzle and rotor.
Considering the assembly of the expander, the length of the straight zones between the volute and nozzle inlet and between the nozzle and the rotor were determined. The trailing edge thickness was determined considering the blade manufacturer. The shape of the leading edge was rounded to prevent separation at the leading edge of the rotor.

Numerical Analysis
Numerical analysis was conducted using NUMECA's FINE/Turbo program. As a turbulence model, the Spalart-Allmaras model, a one-equation model with excellent calculation efficiency, was applied to determine the exact properties of the boundary layer in the pressure gradient.     Figure 4a shows the relation between the expansion ratio and mass flow r ure 4b shows the relation between the expansion ratio and efficiency. The e from the meanline design was recalculated by reflecting the loss in the expand defined as the ratio of the enthalpy difference and isentropic enthalpy differ tween the inlet and outlet of the expander. It was verified that the target efficie satisfied at the design expansion ratio. Figure 4c shows the relation between th sion ratio and the power. Like the efficiency curve, the power from the meanlin was recalculated by reflecting the loss in the expander. Finally, Figure 4d shows tion between the expansion ratio and the expander outlet temperature. The c expander must satisfy the expander outlet temperature in the system requirem shown in the figure, the estimated expander outlet temperature meets the des perature.

3D Shape Generation
Based on the results of the meanline design, a three-dimensional expander try was generated to conduct the numerical analysis. The three-dimensional g was generated using AxCent program of ConceptsNREC. Figure 5 shows the 3D of the nozzle and rotor.
Considering the assembly of the expander, the length of the straight zones the volute and nozzle inlet and between the nozzle and the rotor were determin trailing edge thickness was determined considering the blade manufacturer. Th of the leading edge was rounded to prevent separation at the leading edge of the

Numerical Analysis
Numerical analysis was conducted using NUMECA's FINE/Turbo progra turbulence model, the Spalart-Allmaras model, a one-equation model with e calculation efficiency, was applied to determine the exact properties of the b layer in the pressure gradient.
In order to increase the convergence of the calculation, the choice of bound dition is important. In general, pressure, temperature, and flow rate are ap boundary conditions in turbomachinery numerical analysis. Applicable press In order to increase the convergence of the calculation, the choice of boundary condition is important. In general, pressure, temperature, and flow rate are applied as boundary conditions in turbomachinery numerical analysis. Applicable pressure and temperature values are static and stagnation value. In addition, the applicable flow information is the mass flow rate and the velocity vector. In this study, numerical analysis was performed by analyzing the calculated result using the pressure condition as the main boundary condition. The applied boundary conditions were the inlet total pressure, inlet total temperature, velocity direction, turbulent viscosity, and outlet static pressure. To obtain the off-design point performance, the outlet static pressure was varied. The meanline design results were used for the flow angle information flowing into the nozzle inlet. That is, the absolute flow angle of the volute outlet was used as the nozzle inlet flow angle in the numerical analysis.
A mixing plane condition was applied between the nozzle outlet and the rotor inlet. Therefore, the average physical property value in the pitch direction was applied as the rotor inlet condition. In real operation conditions, the shape of the nozzle wake changes along with the rotor revolution, and to analyze it, unsteady calculations should be conducted. In this study, the focus was on efficiency improvement, and a steady calculation was conducted. The calculation was conducted for a single periodic passage.
To determine the appropriate number of grids, the efficiency according to the number of grids was analyzed. Efficiency was compared for the following number of grids: 199,079, 239,403, 345,415, 534,415, 1,088,147, and 1,332,851. When the number of grids was greater than 239,403, there was no significant change in efficiency. In this study, the number of grid cells in the calculation domain was 1,088,147. Figure 6 shows the relationship of efficiency according to grid number to check grid independence. A hexa mesh was applied to generate a dense grid near the blade wall, and Y+ was set to 2.
Appl. Sci. 2022, 12, x FOR PEER REVIEW sure, inlet total temperature, velocity direction, turbulent viscosity, and out pressure. To obtain the off-design point performance, the outlet static pressure ied. The meanline design results were used for the flow angle information flow the nozzle inlet. That is, the absolute flow angle of the volute outlet was used as zle inlet flow angle in the numerical analysis.
A mixing plane condition was applied between the nozzle outlet and the ro Therefore, the average physical property value in the pitch direction was appli rotor inlet condition. In real operation conditions, the shape of the nozzle wake along with the rotor revolution, and to analyze it, unsteady calculations should ducted. In this study, the focus was on efficiency improvement, and a steady ca was conducted. The calculation was conducted for a single periodic passage.
To determine the appropriate number of grids, the efficiency accordin number of grids was analyzed. Efficiency was compared for the following nu grids: 199,079, 239,403, 345,415, 534,415, 1,088,147, and 1,332,851. When the nu grids was greater than 239,403, there was no significant change in efficiency study, the number of grid cells in the calculation domain was 1,088,147. Figure the relationship of efficiency according to grid number to check grid indepen hexa mesh was applied to generate a dense grid near the blade wall, and Y+ was The cryogenic expander design was conducted using the real gas properti lium. To reduce the computation time during the numerical analysis, a property the working fluid was applied to the calculation. NUMECA's TabGen was used erate the data. TabGen creates a fluid property table in connection with REFPROP, a program that provides fluid properties. Interpolation was perform on the generated property table so that it could be used as a property value at a necessary for numerical analysis. Figure 7 shows the performance curves of the cryogenic expander derived numerical analysis. In the same way as for the meanline design results, the m rate, efficiency, power, and outlet temperature were plotted against the expans at 80%, 90%, and 100% of the design speed. The efficiency and power from the cal analysis were recalculated by taking into account the extra mechanical loss expander. The extra mechanical losses were the loss caused by the volute and caused by the windage in the rotor back disk, and in this study, the losses were to be 2% and 4%, respectively. The outlet temperature of the expander satisfies perature required by the system. The cryogenic expander design was conducted using the real gas properties of helium. To reduce the computation time during the numerical analysis, a property table of the working fluid was applied to the calculation. NUMECA's TabGen was used to generate the data. TabGen creates a fluid property table in connection with NIST's REFPROP, a program that provides fluid properties. Interpolation was performed based on the generated property table so that it could be used as a property value at a location necessary for numerical analysis. Figure 7 shows the performance curves of the cryogenic expander derived from the numerical analysis. In the same way as for the meanline design results, the mass flow rate, efficiency, power, and outlet temperature were plotted against the expansion ratio at 80%, 90%, and 100% of the design speed. The efficiency and power from the numerical analysis were recalculated by taking into account the extra mechanical losses in the expander. The extra mechanical losses were the loss caused by the volute and the loss caused by the windage in the rotor back disk, and in this study, the losses were assumed to be 2% and 4%, respectively. The outlet temperature of the expander satisfies the temperature required by the system. Appl

Optimization Methodology
Optimization is an effective tool for selecting the best design among a variety of available options through identifying design variables and design objectives. Figure 8 shows a schematic of the design optimization process in this study. The initial design geometry was obtained via meanline design, and the performance of the design geometry was evaluated via CFD. The building of a database using design of experiments (DOE) and obtaining and optimizing a metamodel were performed using Minitab software. Metamodel optimization allows special insight into system performance and allows the system to be examined across a wide range of design conditions at a relatively modest computational cost [13].
Optimization starts with building a database using DOE. The number of regression model coefficients increases rapidly as the number of design parameters increases. Various DOEs have been proposed to obtain a highly accurate response surface function with a minimum number of experiments under given conditions. In this study, the Box-Behnken design method was adopted because it can derive quadratic regression equations, easily create orthogonal blocks, and use fewer design points than 3k-factorial design or central composite design [14]. The boundaries of the design space were normalized to the range selected for each design parameter so that the existing geometry value was set to 0 and the maximum range of change for the design variables was set to ±1. For 4 design variables, 25 test cases were designed to obtain the equation relating the design parameters and the objective function. In each test case, two design variables were fixed at 0, the existing geometry value, and the other two variables were set to ±1.

Optimization Methodology
Optimization is an effective tool for selecting the best design among a variety of available options through identifying design variables and design objectives. Figure 8 shows a schematic of the design optimization process in this study. The initial design geometry was obtained via meanline design, and the performance of the design geometry was evaluated via CFD. The building of a database using design of experiments (DOE) and obtaining and optimizing a metamodel were performed using Minitab software. Metamodel optimization allows special insight into system performance and allows the system to be examined across a wide range of design conditions at a relatively modest computational cost [13].
Optimization starts with building a database using DOE. The number of regression model coefficients increases rapidly as the number of design parameters increases. Various DOEs have been proposed to obtain a highly accurate response surface function with a minimum number of experiments under given conditions. In this study, the Box-Behnken design method was adopted because it can derive quadratic regression equations, easily create orthogonal blocks, and use fewer design points than 3k-factorial design or central composite design [14]. The boundaries of the design space were normalized to the range selected for each design parameter so that the existing geometry value was set to 0 and the maximum range of change for the design variables was set to ±1. For 4 design variables, 25 test cases were designed to obtain the equation relating the design parameters and the objective function. In each test case, two design variables were fixed at 0, the existing geometry value, and the other two variables were set to ±1. termined by the R2adj (correction coefficient of determination) of the response function. That is, if R2adj > 0.9, the response surface function is very accurate, and R2adj > 0.7, the response surface function is relatively accurate [15,16].
An optimal point of the response surface provides a new design geomet performance of the optimal design geometry is then evaluated. This optimizati cess repeats until the objective function reaches the target value.  A meta-model was built using the DOE points to construct the function F correlating the objective functions U j with the design parameters X i as follows: In the present study, the response surface method (RSM) was used to obtain an optimal design. It is widely used among approximate model-based design optimization techniques. The RSM is a technique that approximates the response surface function from the value of the objective function obtained after performing 3D CFD on a set of design variables. The response surface function U j can be expressed as a quadratic polynomial using the design parameters X i and X j as follows: where β O , β i , β ii , and β ik are the polynomial coefficients, and these undefined coefficients are determined via least square regression using the 3D CFD results. A detailed analysis of the effect of parameters and their interactions on the response was also performed through the surface plots using Minitab software. In the response surface method, the reliability of the response surface function can be checked, and this reliability can be increased by removing insignificant undefined coefficients, which is called response surface model pooling. Response surface model reduction is achieved through analysis of variance (ANOVA) and regression analysis. The analysis of variance is used to determine whether each term of the response surface function is significant, and the significance is confirmed under the condition that the significance probability p is <0.05. Insignificant terms are removed in the order of block, alternating term, and square term, and analysis of variance is performed. After the analysis of variance is finished, regression analysis is performed. Regression analysis is an analysis method that measures the fitness of how well the response surface function represents the actual value, and it is determined by the R2 adj (correction coefficient of determination) of the response surface function. That is, if R2 adj > 0.9, the response surface function is very accurate, and if 0.9 > R2 adj > 0.7, the response surface function is relatively accurate [15,16]. An optimal point of the response surface provides a new design geometry. The performance of the optimal design geometry is then evaluated. This optimization process repeats until the objective function reaches the target value.

Optimization of the Expander
Since the aim of the optimization in the present work is to maximize the total-to-total efficiency for the given design specification, the isentropic total-to-total efficiency of the expander at the design point was considered the objective function, as Equation (5) shows: where h T0 is the expander inlet total enthalpy, h T6 is the expander outlet total enthalpy, and h isen T6 is the expander outlet total enthalpy in isentropic conditions of the expander inlet. The hub and shroud meridional contours in the z-r plane and the blade β angle distributions at the hub and shroud were chosen as the design parameters of the 3D radial expander rotor. The hub and shroud meridional contours and the control points of each Bezier curve are shown in Figure 9. The locations of the inlet and outlet points were fixed. To maximize the shape profiles with the minimum number of test cases, and Z H2 , R H3 , Z S2 , and R S3 were selected as the design parameters for the optimization of the meridional contours. Figure 10 shows the blade β angle distributions at the hub and shroud and the control points of each Bezier curve. β H3 , β H4 , β S3 , and β S4 were selected as the design parameters for the optimization of the blade angle distributions.

Optimization of the Expander
Since the aim of the optimization in the present work is to maximize the total-to-total efficiency for the given design specification, the isentropic total-to-total efficiency of the expander at the design point was considered the objective function, as Equation (5) shows: is the expander inlet total enthalpy, is the expander outlet total enthalpy, and is the expander outlet total enthalpy in isentropic conditions of the expander inlet.
The hub and shroud meridional contours in the z-r plane and the blade angle distributions at the hub and shroud were chosen as the design parameters of the 3D radial expander rotor. The hub and shroud meridional contours and the control points of each Bezier curve are shown in Figure 9. The locations of the inlet and outlet points were fixed. To maximize the shape profiles with the minimum number of test cases, and 2 , 3 , 2 , and 3 were selected as the design parameters for the optimization of the meridional contours. Figure 10 shows the blade angle distributions at the hub and shroud and the control points of each Bezier curve.
3 , 4 , 3 , and 4 were selected as the design parameters for the optimization of the blade angle distributions.  For the 4 design variables, 25 test cases were designed using the Box-Behnken design method using Minitab software. The design of 25 test cases alternated between the shape of the impeller meridional contours and the blade β angle. In this study, the optimization process for the blade β angle distributions was first conducted with the initially given geometry for the hub and shroud meridional contours. An optimal design for the blade β angle distributions was achieved using the results of the 25 test cases. The next optimization process was performed for the hub and shroud meridional contours with the obtained geometry of the blade β angle distributions from the first optimization process. The 25 CFD test cases presented an optimized geometry for the impeller meridional contours. As such, the optimization process alternated between optimization of the blade β angle distributions and that of the impeller hub and shroud meridional contours. In this study, eight optimization cycles were executed. The metal models of optimizations are summarized in Table 5, and the enhancement of the isentropic total-to-total efficiency of the expander is presented in Figure 11. Figures 12 and 13 show comparisons of the blade β angle distributions and meridional contours obtained by optimization. For the 4 design variables, 25 test cases were designed using the B sign method using Minitab software. The design of 25 test cases alterna shape of the impeller meridional contours and the blade angle. In th timization process for the blade angle distributions was first conduc tially given geometry for the hub and shroud meridional contours. An for the blade angle distributions was achieved using the results of t The next optimization process was performed for the hub and shroud tours with the obtained geometry of the blade angle distributions fro mization process. The 25 CFD test cases presented an optimized geom peller meridional contours. As such, the optimization process alternate mization of the blade angle distributions and that of the impeller h meridional contours. In this study, eight optimization cycles were exec models of optimizations are summarized in Table 5, and the enhancem tropic total-to-total efficiency of the expander is presented in Figure 11 13 show comparisons of the blade angle distributions and meridion tained by optimization.   Table 5. Design parameters and metamodel in each optimization cycle.

Cycle Design Parameters Metamodel
meridional contours U = 0.875270 − 0.003063 Z H2 − 0.001480 R H3 + 0.003210 Z S2 + 0.001197 R S3 − 0.000112 Z H2 Z H2 − 0.000010 R H3 R H3 − 0.000145 Z S2 Z S2 − 0.000180 R S3 R S3 + 0.000207 Z H2 R H3 + 0.000153 Z H2 Z S2 − 0.000100 Z H2 R S3 − 0.000222 R H3 Z S2 − 0.000020 R H3 R S3 + 0.000170 Z S2 R S3 3 blade β angle meridional contours meridional contours     Figure 14 shows the locations of Planes 1 and 2 in the calculation domain. Plane 1 is the mid-span of the rotor and Plane 2 is the passage area between the blade pressure side and blade suction side. Figures 15 and 16 show the entropy distribution of the rotor,    Figure 14 shows the locations of Planes 1 and 2 in the calculation domain. Plane 1 is the mid-span of the rotor and Plane 2 is the passage area between the blade pressure side and blade suction side. Figures 15 and 16 show the entropy distribution of the rotor, and Figures 17 and 18 show the relative velocity distributions at Planes 1 and 2. As shown in the figures, the entropy of the original model at the suction side of the rotor was larger than that of the optimized model.   Through the optimization, the blade height was increased and the blade β angle was distributed smoothly from the blade leading edge to the trailing edge. The increased rotor blade height induced an increase in the passage area. Therefore, the flow velocity of the optimized model was slightly lower than that of the original model.

Optimization Results at Design Point
As shown in Figure 17, in the original model, the flow area which is not the main passage flow component near the suction blade surface was formed wider than the optimized model. That is, through optimization, the flow of the main passage flow could be smoother than the original model.  As shown in Figure 17, in the original model, the flow area which is not the main passage flow component near the suction blade surface was formed wider than the optimized model. That is, through optimization, the flow of the main passage flow could be smoother than the original model.
As shown in Figure 18, in the case of the original model, a strong tip leakage flow occurs from the pressure blade side to the suction blade side. In the original model, a As shown in Figure 17, in the original model, the flow area which is not the main passage flow component near the suction blade surface was formed wider than the optimized model. That is, through optimization, the flow of the main passage flow could be smoother than the original model.
As shown in Figure 18, in the case of the original model, a strong tip leakage flow occurs from the pressure blade side to the suction blade side. In the original model, a stronger tip leakage flow occurred in the relatively narrow passage area than in the op-  Figure 18, in the case of the original model, a strong tip leakage flow occurs from the pressure blade side to the suction blade side. In the original model, a stronger tip leakage flow occurred in the relatively narrow passage area than in the optimized model, resulting in increased losses and consequently affecting efficiency. Figure 19 and Table 6 show a comparison of the performance of the original and optimized models. In Figure 19a, at the design point (DP), it can be seen that the efficiency of the optimized model improved 1.98% compared to that of the original model. Under the condition of an expansion ratio of 1.34 or higher, it was confirmed that the efficiency of the optimization model was improved compared to that of the original model ( Table 6). The extra mechanical losses which occur at the volute and the rotor back disk were considered. Figure 19b shows the temperature trends at the expander outlet relative to the expansion ratio. The outlet temperature of the optimized model was 0.32% lower than that of the original model at the design point. Figure 19c shows the relationship between the expansion ratio and the expander output power. At the design point, the power of the optimized model was approximately 2% higher than that of original model.

Conclusions
In this study, we introduced a design method for expanders for hydrogen liquefaction plants. Blade optimization was conducted to increase the expander efficiency. The capacity of the developed expander for a hydrogen liquefaction system is 0.5 tons of hydrogen liquefied per day. Based on the cycle analysis results of the hydrogen liquefaction system, the meridional and three-dimensional design of the expander were conducted. To improve the efficiency of the expander, the rotor blade meridional shape and blade β angle were optimized. Through the optimization design, the blade height was increased between the blade leading edge and the trailing edge, and the distribution of the blade β angle was rearranged smoothly. Finally, the optimized expander's efficiency increased by 1.98% compared to that of the original model. To investigate the reason for this efficiency improvement, a flow analysis was conducted. The entropy distribution on the suction blade surface was smaller in the optimized model than that of the original model. Through the relative velocity distribution, it was found that the optimized model had a smoother main passage flow and smaller tip leakage flow than the original model. Due to the passage area of the optimized model, which is larger than the original model, the effect of the tip leakage flow on the main passage flow could be reduced. It was confirmed that the optimization of the blade β angle, rotor hub, and shroud meridional contour had an effect on improving the efficiency of the cryogenic expander. The efficiency optimization method in this paper can be used as a reference in designing high-efficiency cryogenic expanders in the future.

Conflicts of Interest:
The authors declare no conflict of interest.