Numerical Optimization of the β -Type Stirling Engine Performance Using the Variable-Step Simpliﬁed Conjugate Gradient Method

: This study focuses on optimizing a 100-W-class β -Type Stirling engine by combining the modiﬁed thermodynamic model and the variable-step simpliﬁed conjugate gradient (VSCGM) method. For the modiﬁed thermodynamic model, non-uniform pressure is directly introduced into the energy equation, so the indicated power and heat transfer rates can reach energy balance while the VSCGM is an updated version of the simpliﬁed conjugate gradient method (SCGM) with adaptive increments and step lengths to the optimization process; thus, it requires fewer iterations to reach the optimal solution than the SCGM. For the baseline case, the indicated power progressively raises from 88.2 to 210.2 W and the thermal efﬁciency increases from 34.8 to 46.4% before and after optimization, respectively. The study shows the VSCGM possesses robust property. All optimal results from the VSCGM are well-matched with those of the computational ﬂuid dynamics (CFD) model. Heating temperature and rotation speed have positive effects on optimal engine performance. The optimal indicated power rises linearly with the charged pressure, whereas the optimal thermal efﬁciency tends to decrease. The study also points out that results of the modiﬁed thermodynamic model with ﬁxed values of unknowns agree well with the CFD results at points far from the baseline case.


Introduction
Many configurations of Stirling engines possess dominantly longitudinal variations of temperature and velocity. This intriguing feature favors thermodynamic models, which exhibit the one-dimensional nature [1][2][3][4][5]. In addition, thermodynamic models exploit several simplified assumptions, such as uniform temperature and pressure in each chamber, ideal gas equation of state, adiabatic conditions, etc. Consequently, thermodynamic models consume much less computational time and memory storage than CFD ones. These striking merits adapt thermodynamics models to many applications, for example, investigation of thermodynamic performance of Stirling engines [6,7], study of dynamic behaviors of Stirling engines [8][9][10], and optimization of Stirling engine performance [1,11]. According to the history of theoretical model development, there exist several kinds of thermodynamic models based on the complexity level. The isothermal model uses assumptions about the constant temperature in each chamber and ignores pressure loss [12]. As a result, this kind of model generates less accurate results. The ideal adiabatic model assumes constant temperatures in heater, cooler, and regenerators and adiabatic conditions on surfaces of compression and expansion chambers [13]. The non-ideal adiabatic model introduces at least one of the following effects: convection loss, shuttle loss, gas spring hysteresis, the effectiveness of regenerator, and pressure loss [14][15][16]. The modified non-ideal adiabatic model, proposed by Yang and Cheng [7], overcomes several shortcomings of the non-ideal adiabatic thermodynamic model by adding a temporal variation of temperature and pressure loss in the heater, cooler, and regenerator and introducing mechanical loss to obtain shaft power. Recently, Cheng and Phung [17] proposed a modified thermodynamic model by completely removing the adiabatic condition on expansion and compression chambers and simultaneously introducing non-uniform pressure to the energy equation so cyclic heat transfer rates and cyclic indicated power can balance at the final cycle. This energy balance lays a firm foundation for optimization in this study. In contrast, CFD models use fewer assumptions than thermodynamic ones and can expand the transient evolution of thermal and flow fields in three-dimensional space at the expense of higher computational time and memory resource [18][19][20]. There appears no direct application of CFD models for optimizing Stirling engine performance due to these severe limitations.
The design and optimization of the Stirling engines are of strong correlation and challenging phases to approach a new prototype engine with the highest performance. Patel and Savsani [21] minimized pressure losses and maximized the power and the thermal efficiency by their proposed multi-objective tutorial training and self-learning-inspired teaching-learning-based optimization method for the Stirling engines. This method can optimize several objective functions simultaneously. Duan et al. [22] exploited the multiobjective particle swarm optimization algorithm and the Pareto optimal frontier to optimize the irreversibility, power, and thermal efficiency. They considered not only the geometry parameters, but also the temperature of working gas and the charged pressure as design variables. Ahmadi et al. [23] combined a non-dominated sorting genetic algorithm and finite speed thermodynamic analysis to gain the optimal output power and thermal efficiency and minimize the pressure losses. The results are in good agreement with the experimental data. Arora et al. [24] optimized the thermo-economic value and engine performance using NSGA-II in MATLAB Simulink for the parabolic-dish concentrator Stirling engine. Xiao et al. [25] conducted a multi-objective optimization based on pressure and volume, given from the CFD analysis.
The conjugate gradient method (CGM) is an efficient optimization algorithm based on the exploitation of first-order derivatives of the objective function and the quadraticconvergence property [26]. The SCGM, proposed by Cheng and Chang [27], is a modified version of the CGM and fixes the step length and approximates the gradient of the objective function by finite difference schemes. The VSCGM, proposed by Cheng and Lin [28], is the most updated version of the SCGM with varying design-variable increments and automatically adjustable step length, which adapts to the optimization process. Therefore, this version reduces the number of optimization iterations significantly.

Study Aims
From the literature review, the CFD method is inappropriate for optimizing engine performance due to its limitations. In this study, the performance of the compact 100-W-class β type Stirling engine in Ref. [18] is optimized by the VSCGM, proposed by Cheng and Lin [28], where the modified thermodynamic model, proposed by Cheng and Phung [17], plays the role of a direct solver, which generates the indicated power and the thermal efficiency necessary for estimating the objective function. The optimization method exploits both the highly convergent speed of the VSCGM and the advantages of the modified thermodynamic model such as computational time, memory resource, and the cyclic balance between heat transfer and indicated power. The latter contributes the non-uniformity of pressure to the optimized engine performance and lays a firm foundation for multiobjective optimization. Results before optimization and after optimization are doubly checked with the aid of the CFD model.
However, the limitation of the optimization method is only to work with the continuous design variables and the smooth objective function. Section 3.1 describes the proposed smooth objective function, while Section 4 lists selected continuous design variables.

Numerical Methods
Section 3.1 proposes a formula of the objective function and succinctly introduces the VSCGM, which takes the role of an optimization method. Section 3.2 discusses a modified Energies 2021, 14, 7835 3 of 14 thermodynamic model which computes the engine performance to estimate values of the objective function. For the sake of simplicity of terminology, the integration of the modified thermodynamic model into the VSCGM creates an optimizer called the VSCGM optimizer. Section 3.3 describes briefly the CFD model with the role of doubly checking the results obtained from the VSCGM optimizer and determining values of unknowns in the modified thermodynamic model.

Variable-Step Simplified Conjugate Gradient Method
In this study, the multi-goal objective function depends on the engine performance as follows: where C 1 and C 2 are both non-negative and C 1 + C 2 = 1 and the respective coefficients C 1 and C 2 define how much indicated power • W and thermal efficiency η contribute to the objective function J. In this study, reference values, denoted by a subscript "ref", are identical to the initial guess so the value of the objective function is initially equal to 1. During the optimization process, the objective function is bounded between 0 and 1.
From the partial derivatives of the objective function with respect to engine performance and the triangle inequality, we can define a convergence-check function ∆J for the optimization process as follows: where δ is the convergent criterion. The convergent-check function ∆J contains rich information about the optimization process: (1) ratio of thermal efficiency and indicated power at the current optimization iteration to those of the initial guess; (2) the relative change in engine performance between two consecutive iterations; and (3) contribution coefficients of thermal efficiency and indicated power to the objective function. In this study, δ is set equal to 10 −5 for all optimization cases. If the iteration satisfied the inequality (2), the optimization process will be terminated. The VSCGM is the most updated version of the simplified conjugate gradient method (SCGM), where the design-variable increments and the step lengths in the VSCGM are automatically adjusted to favor the optimization process. Consequently, the VSCGM requires fewer optimization iterations than the SCGM [28]. Design-variable increment ∆X i depends on the initial increment in design variable ∆X (1) i , initial step length β (1) i , and previous step length β (n−1) i as follows: The gradient of the objective functions is approximately evaluated based on the central difference scheme as follows: where the vector of design variables (X 1 , X 2 , . . . , X M ) T is denoted by the symbol X and unit vector e i = (0, . . . , 1, . . . , 0) T has all components of zero except 1 at the ith component. The search direction is the linear combination of the current gradient and the previous search direction multiplied with a coefficient for modifying the current gradient direction.
where γ i is the gradient-component ratio given by: The current step length can be determined based on ratio of search directions at the current and previous step R i , the previous step length β (n−1) i , and the exponent G i as follows: β where: The design variables are then updated as follows: Table 1 lists specifications of the VSCGM used in this study, while Figure 1 shows the flowchart of the VSCGM.

Thermodynamic Model
In the modified thermodynamic model, proposed by Cheng and Phung [17], pressure losses are directly introduced into the energy equation, so heat transfer rates and indicated

Thermodynamic Model
In the modified thermodynamic model, proposed by Cheng and Phung [17], pressure losses are directly introduced into the energy equation, so heat transfer rates and indicated power meet the energy balance at the final cycle. This builds a solid foundation for multigoal optimization based on these two parameters, as mentioned in Section 3.1. The main calculations of the modified thermodynamic model can be summarized as follows: Uniform pressure in the working gas domain: Instantaneous mass in each chamber: Mass difference between the current time and the previous time: Mass flow rate across the interface of two successive sub-chambers: Non-uniform pressure in the working gas domain: where: The relative pressure of the kth sub-chamber to the central sub-chamber of the regenerator, P rel,k , is estimated as: P rel,k reg = 0 (17) The pressure drop along the kth sub-chamber is computed as: The friction factor for the compression chamber, the expansion chamber, and the heater tubes: The friction factor for the regenerator: The friction factor for the cooler fins: The instantaneous temperature can be updated by the following difference equation: where Figure 2 illustrates the way to solve the governing equations in the modified thermodynamic model. The thermal efficiency and the indicated power • from the modified thermodynamic model are employed to estimate the multi-goal objective function in the VSCGM. These two parameters are estimated as follows:

CFD Model
The CFD model consists of transport equations for mass, momentum, and energy. Additionally, equations for turbulent kinetic energy and its dissipation are included. It should be emphasized that the CFD model is used for doubly comparing optimal results from the VSCGM optimizer.
Continuity equation: Momentum equation: The thermal efficiency η and the indicated power • W from the modified thermodynamic model are employed to estimate the multi-goal objective function in the VSCGM. These two parameters are estimated as follows:

CFD Model
The CFD model consists of transport equations for mass, momentum, and energy. Additionally, equations for turbulent kinetic energy and its dissipation are included. It should be emphasized that the CFD model is used for doubly comparing optimal results from the VSCGM optimizer.
Continuity equation: Momentum equation: Energy equation: Dissipation rate of turbulent kinetic energy ε: Turbulent eddy: For the sake of brevity, readers can refer to the detailed specifications of the compact 100-W class β-type Stirling engine and validation of the CFD model in Ref. [18]. Initialboundary conditions of the CFD model are listed in Table 2. In this study, the engine speed is given, so the expressions of the piston and displacer velocities are as follows: In the CFD model, the indicated power and thermal efficiency are computed as mentioned in Equations (25) and (26). However, P e and P c are the volume-averaged pressure of the expansion chamber and the compression chamber, respectively, as mentioned in Ref. [18].
The CFD model serves two purposes: (1) completing our proposed modified thermodynamic model, as mentioned in detail in Ref. [17], and (2) verifying optimal solutions, which are obtained from the combination of VSCGM and the modified thermodynamic model. The matching process between the CFD model and the modified thermodynamic model for the baseline case generates values of unknowns listed in Table 3. These values remain unchanged for all optimization cases throughout the study. As will be seen later, results from the VSCGM optimizer with fixed values of unknowns and the CFD model are in good agreement. The schematic of the β-type Stirling engine is displayed in Figure 3. Table 3. Values of unknowns in the thermodynamic model [17].

Parameter Value
Expansion chamber, R k (K/W) 0.624 Heater tubes, R k (K/W) 0.641 Cooler fins, R k (K/W) 1.087 Compression chamber, R k (K/W) 0.892 Coefficient C in Equation (22) 700 In the CFD model, the indicated power and thermal efficiency are computed as mentioned in Equations (25) and (26). However, and are the volume-averaged pressure of the expansion chamber and the compression chamber, respectively, as mentioned in Ref. [18].
The CFD model serves two purposes: (1) completing our proposed modified thermodynamic model, as mentioned in detail in Ref. [17], and (2) verifying optimal solutions, which are obtained from the combination of VSCGM and the modified thermodynamic model. The matching process between the CFD model and the modified thermodynamic model for the baseline case generates values of unknowns listed in Table 3. These values remain unchanged for all optimization cases throughout the study. As will be seen later, results from the VSCGM optimizer with fixed values of unknowns and the CFD model are in good agreement. The schematic of the β-type Stirling engine is displayed in Figure  3.

Results and Discussion
In this study, a 100-W-class β-Type Stirling engine, mentioned in Ref. [18], is optimized by the VSCGM optimizer. Readers can refer to the specifications of the engine in Table 1 of Ref. [18]. Heater length, regenerator length, cooler length, and piston diameter are herein selected as design variables, illustrated in Figure 3. Their initial values and bounds are listed in Table 4. In this study, the values of C 1 and C 2 are set equal to 0.5.  Figure 4a shows that the optimization process terminates after 17 iterations for the baseline case. The value of the objective function goes down monotonically from 1 to 0.5856. The indicated power progressively raises from 88. 2 Figure 4a shows that the optimization process terminates after 17 iterations for the baseline case. The value of the objective function goes down monotonically from 1 to 0.5856. The indicated power progressively raises from 88.2 to 210.2 W and the thermal efficiency increases from 34.8 to 46.4%. Figure 4b shows the CFD model and the VSCGM optimizer generate well-matched P-V curves before and after the optimization process. Additionally, it should be noted that there is a significant enlargement of the P-V area at the end of the optimization process in comparison to that at the beginning. This can be explained as follows. During the optimization iteration, a gradual increase in piston diameter, causing the right shift of the compression-chamber P-V area, enhances not only the swept volume but the compression ratio. It leads to an increase in pressure during the expansion process and the resulting indicated power. Moreover, it elevates the Reynolds number, so the heat transfer becomes more effective and the temperature of working gas becomes closer to wall temperatures. It also explains the fact that the thermal efficiency increases monotonically, as seen in Figure 4a.  To verify the results obtained from the VSCGM optimizer, a CFD simulation case is run with all optimal values of design variables from the VSCGM optimizer. Table 5 shows that the results of the two methods are predicted with a slight difference. Table 5. Verification of the optimal solution for the baseline case.   Figure 4b shows the CFD model and the VSCGM optimizer generate well-matched P-V curves before and after the optimization process. Additionally, it should be noted that there is a significant enlargement of the P-V area at the end of the optimization process in comparison to that at the beginning. This can be explained as follows. During the optimization iteration, a gradual increase in piston diameter, causing the right shift of the compression-chamber P-V area, enhances not only the swept volume but the compression ratio. It leads to an increase in pressure during the expansion process and the resulting indicated power. Moreover, it elevates the Reynolds number, so the heat transfer becomes more effective and the temperature of working gas becomes closer to wall temperatures. It also explains the fact that the thermal efficiency increases monotonically, as seen in Figure 4a.
To verify the results obtained from the VSCGM optimizer, a CFD simulation case is run with all optimal values of design variables from the VSCGM optimizer. Table 5 shows that the results of the two methods are predicted with a slight difference. The VSCGM is one of the local robust optimization techniques, which is inherited from the CGM and the SCGM. It means the VSCGM converges to a local minimum point with all guessed points close enough to the local optimal point. To illustrate this property of the VSCGM method, piston diameter, heater length, and regenerator length are selected. Table 6 summarizes five guessed sets of their values. Figure 5 shows that the five guessed points are spatially scattered in the design-variable space and all five trajectories from these initial points converge to a single optimal point. It proves that the VSCGM possesses a robust property.  Figure 6 indicates the effect of heating temperature on the optimal solution obtained from the VSCGM optimizer. The heating temperature varying from 673 to 1173 K is considered. Before optimization, the comparison between the modified thermodynamic model with values of unknown coefficients in Table 3 and the CFD simulation obtained from the reference [18] is performed. The engine performance run by the two models is in good agreement. From the modified thermodynamic model, the original indicated power varies monotonically from 29.1 to 116.1 W and the original thermal efficiency does from 14.6 to 41.6%. After optimization, the optimal indicated power is raised from 98.7 W to 262.0 W and the optimal thermal efficiency is from 26.4 to 53.4% over this range of heating temperature. From classical thermodynamics, a rise in heating temperature raises engine performance. It explains why the increase in heating temperature has a positive effect on optimal engine performance. To doubly check the results from the VSCGM optimizer, CFD simulation cases are performed with the optimal configuration specifications obtained from the VSCGM optimizer. The CFD results agree well with those from the VSCGM optimizer. The maximum difference in the indicated power between the two models is 5.5 W, while the maximum difference in the thermal efficiency is only 1.1%. Figure 7 shows the effect of charged pressure on the optimal engine performance over the range from 3 to 15 bar. The thermal efficiency and the indicated power before optimization are well predicted by the modified thermodynamic model with values of unknown coefficients in Table 3 and the CFD model, presented in Ref. [18]. The indicated power varies linearly from 88.2 W to 526.0 W and the thermal efficiency varies linearly from 34.8 to 37.9% with a charged pressure from 3 to 15 bar. The thermal efficiency tends to decrease when the charged pressure is over 9 bar. After optimization, the optimal indicated power linearly surges from 210.2 to 997.5 W and the optimal thermal efficiency decreases slightly from 46.4 to 43.7% as the charged pressure varies from 3 to 15 bar. The double-check from the CFD model shows that the maximum difference in the optimal thermal efficiency is 2.2% and that in the optimal indicated power is 25.2 W at a 15-bar charged pressure.   Figure 6 indicates the effect of heating temperature on the optimal solution obtained from the VSCGM optimizer. The heating temperature varying from 673 to 1173 K is considered. Before optimization, the comparison between the modified thermodynamic model with values of unknown coefficients in Table 3 and the CFD simulation obtained from the reference [18] is performed. The engine performance run by the two models is in good agreement. From the modified thermodynamic model, the original indicated power varies monotonically from 29.1 to 116.1 W and the original thermal efficiency does from 14.6 to 41.6%. After optimization, the optimal indicated power is raised from 98.7 W to 262.0 W and the optimal thermal efficiency is from 26.4 to 53.4% over this range of heating temperature. From classical thermodynamics, a rise in heating temperature raises engine performance. It explains why the increase in heating temperature has a positive effect on optimal engine performance. To doubly check the results from the VSCGM optimizer, CFD simulation cases are performed with the optimal configuration specifications obtained from the VSCGM optimizer. The CFD results agree well with those from the VSCGM optimizer. The maximum difference in the indicated power between the two  Figure 7 shows the effect of charged pressure on the optimal engine performance over the range from 3 to 15 bar. The thermal efficiency and the indicated power before optimization are well predicted by the modified thermodynamic model with values of unknown coefficients in Table 3 and the CFD model, presented in Ref. [18]. The indicated  Figure 6. Variation of the optimal engine performance with heating temperature.  Figure 8 points out the effect of the rotation speed on the optimal engine performance over the range from 900 to 2100 rpm. Before optimization, the CFD model and the modified thermodynamic model in the VSCGM optimizer show a good agreement. Their respective maximum differences in the indicated power and the thermal efficiency are 2.4 W at 1350 rpm and 1.1% at 2100 rpm. After optimization, the raise in rotation speed over this range leads to an improvement in both the optimal indicated power and the optimal thermal efficiency. The optimal thermal efficiency goes down from 51.1 to 38.7% and the optimal indicated power rises from 146.7 to 226.6 W as the rotation speed increases from 900 to 2100 rpm. The comparison between the CFD model and the VSCGM optimizer indicates that the respective maximum differences in the optimal thermal efficiency and the optimal indicated power are 1.7% and 15.9 W at 2100 rpm. Figure 6. Variation of the optimal engine performance with heating temperature. Figure 7 shows the effect of charged pressure on the optimal engine performance over the range from 3 to 15 bar. The thermal efficiency and the indicated power before optimization are well predicted by the modified thermodynamic model with values of unknown coefficients in Table 3 and the CFD model, presented in Ref. [18]. The indicated power varies linearly from 88.2 W to 526.0 W and the thermal efficiency varies linearly from 34.8 to 37.9% with a charged pressure from 3 to 15 bar. The thermal efficiency tends to decrease when the charged pressure is over 9 bar. After optimization, the optimal indicated power linearly surges from 210.2 to 997.5 W and the optimal thermal efficiency decreases slightly from 46.4 to 43.7% as the charged pressure varies from 3 to 15 bar. The double-check from the CFD model shows that the maximum difference in the optimal thermal efficiency is 2.2% and that in the optimal indicated power is 25.2 W at a 15-bar charged pressure.  Figure 7. Variation of the optimal engine performance with the charged pressure.
Energies 2021, 14, x FOR PEER REVIEW 14 of 16 Figure 8 points out the effect of the rotation speed on the optimal engine performance over the range from 900 to 2100 rpm. Before optimization, the CFD model and the modified thermodynamic model in the VSCGM optimizer show a good agreement. Their respective maximum differences in the indicated power and the thermal efficiency are 2.4 W at 1350 rpm and 1.1% at 2100 rpm. After optimization, the raise in rotation speed over this range leads to an improvement in both the optimal indicated power and the optimal thermal efficiency. The optimal thermal efficiency goes down from 51.1 to 38.7% and the optimal indicated power rises from 146.7 to 226.6 W as the rotation speed increases from 900 to 2100 rpm. The comparison between the CFD model and the VSCGM optimizer indicates that the respective maximum differences in the optimal thermal efficiency and the optimal indicated power are 1.7% and 15.9 W at 2100 rpm. This study also shows that the modified thermodynamic model in the VSCGM optimizer with fixed predicted values of unknown coefficients can generate results wellmatched with those of the CFD model at points far from the baseline case.  Figure 8. Variation of the optimal engine performance with the rotation speed. This study also shows that the modified thermodynamic model in the VSCGM optimizer with fixed predicted values of unknown coefficients can generate results wellmatched with those of the CFD model at points far from the baseline case.

Conclusions
In this study, the compact 100-W-class β-Type Stirling engine was optimized by the combination of the modified thermodynamic model, proposed by Cheng and Phung [17], and the VSCGM, proposed by Cheng and Lin [28]. All optimal results from the VSCGM optimizer are compared with those of the CFD model. The following remarking conclusions are drawn: 1.
The objective function reduces gradually from 1 to 0.5856 after 17 iterations. The indicated power raises from 88.2 W to 210.2 W and the thermal efficiency increases from 34.8 to 46.4%. These results from the VSCGM optimizer are doubly checked by the CFD model with the same optimal configuration. Their P-V diagrams before and after optimization are well matched. Furthermore, the difference in the optimal thermal efficiency is 0.3% and the difference in the optimal indicated power is 3.2 W.

2.
The VSCGM optimizer gives the same optimal solution for scattered initial guesses. It proves that the VSCGM possesses robustness.

3.
The heating temperature has a positive effect on optimal engine performance. The optimal indicated power raises from 98.7 to 262.0 W and the optimal thermal efficiency raises from 26.4 to 53.4% as the temperature varies from 673 to 1173 K. The maximum difference in the optimal indicated power between the CFD model and the VSCGM optimizer is 5.5 W, while that in the optimal thermal efficiency is only 1.1%.

4.
The optimal thermal efficiency decreases slightly from 46.4 to 43.7%, while the optimal indicated power raises from 210.2 to 997.5 W as the charged pressure varies from 3 to 15 bar. The double-check from the CFD model shows that the case of 15-bar charged pressure gives the maximum difference between the CFD model and the VSCGM optimizer. These maximum differences are 2.2% for the optimal thermal efficiency and 25.2 W for the optimal indicated power.

5.
The optimal thermal efficiency goes down from 51.1 to 38.7% and the optimal indicated power surges from 146.7 to 226.6 W as the rotation speed increases from 900 to 2100 rpm. The comparison shows that the maximum difference in the optimal indicated power is 15.9 W, while that in the optimal thermal efficiency is 1.7% at 2100 rpm between the CFD model and the VSCGM optimizer. Funding: This research received no external funding.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declared no conflict of interest to this paper.