Numerical Optimization of a Four-Cylinder Double-Acting Stirling Engine Based on Non-Ideal Adiabatic Thermodynamic Model and SCGM Method

: The aim of this study is to optimize a four-cylinder, double-acting α -type Stirling engine with wobble-yoke mechanism using an optimization scheme incorporated with an e ﬃ cient thermodynamic model. In this study, the non-ideal adiabatic thermodynamic model is improved by taking into account factors including pressure drops due to the sudden expansion or contraction of ﬂow cross-sectional areas in the engine, multiple nodes in the regenerator adopted to accurately capture the temperature gradient in the regenerator, and the dependence of the transport properties (thermal conductivity and dynamic viscosity) of the working ﬂuid on temperature and pressure. A parametric analysis is ﬁrstly performed to identify the designed parameters that need to be optimized. In this study, engine optimization is carried out by using the simpliﬁed conjugate-gradient method (SCGM). The e ﬀ ects of the weighting coe ﬃ cients of the objective function are studied. For a particular case considered, the optimization successfully elevates the power output from 1062.56 to 1659.72 W, and thermal e ﬃ ciency from 27.41% to 37.22%. Furthermore, the robustness of the optimization method is tested by giving di ﬀ erent sets of initial guesses. It is found that the present approach can stably lead to the same optimal design and is independent of the initial guess. [13] implemented the genetic algorithm for the optimization of a solar dish Stirling engine. This approach is a stochastic method which covers a global search of optimal geometrical parameters. However, it is less accurate and consumes more computation time. Hooshang et al. [14] describes the optimization of a gamma-type Stirling engine using neural networks as a pattern recognition tool to construct a mathematical model. This method requires a large amount of data for the neural network model to be accurate. On the other hand, Zare and Tavakolpour-Saleh [15] investigated the inﬂuence of dominant poles’ places on a free piston Stirling engine by means of the particle swarm optimization scheme.


Introduction
Two of the most critical issues in recent decades are the problems of global warming and energy crisis. These issues have been clearly addressed in a special report on global warming [1], which states that limiting global warming to a 1.5 • C target would require unprecedented changes in all aspects of society before unimaginable disasters happen. The report also suggested energy system transition as one significant change that needs to be made to combat climate change. In recent years, there has been a strong push to develop renewable energy options, such as solar energy and geothermal energy [2,3]. As an external combustion engine, the Stirling engine can be accommodated with various types of heat source. Egas and Clucas [4] presented a very detailed study on different engine configurations (alpha, beta and gamma) with rhombic-drive and Ross yoke mechanisms. It has been recognized that the Stirling engine is an effective power machine to harvest solar energy, geothermal energy and even waste heat in an environmentally friendly operation. For example, Sunpower Inc. developed a 9-kW free-piston Stirling engine to be used on a solar dish system, operated by Cummins Power Generation Inc., to harvest solar thermal energy [5].
In order to design a Stirling engine for delivering a high-power output, one would attempt to develop an engine with multiple cylinders. The four-cylinder double-acting α-type Stirling engine configuration shown in Figure 1 is one of the promising structures which has expansion and compression Energies 2020, 13 chambers located at different sides of each of the four pistons. The piston motion is affected by the pressures exerted on the upper and the lower faces of the piston, hence the name 'double-acting'. One working zone of the double-acting Stirling engine connects the compression chamber of one cylinder through the flow ducts to the cooler, regenerator, heater and, finally, the expansion chamber of another cylinder. If this engine configuration is compared to four conventional α-type Stirling engines, the numbers of cylinders and pistons are reduced by half. This greatly increases the power density of the engine. With four cylinders, the phase angle of oscillation between two neighboring pistons is specified as 90 • . In 1977, Clucas and Raine [6] proposed a wobble-yoke mechanism that was used to drive the moving parts of the four-cylinder double-acting α-type Stirling engine. Later, García-Canseco et al. [7] proposed a model for control of the wobble-yoke Stirling engine. Chatterton and Pennacchi [8] have conducted a study on different multi-cylinder engine configurations and their thermodynamic performances. Urieli and Berchowitz [9] presented an ideal isothermal thermodynamic model and Yu et al. [10] presented a non-ideal adiabatic thermodynamic model to predict the performance of Stirling engines. Tan and Cheng [11] performed three-dimensional computational fluid dynamic analysis of the four-cylinder, double-acting Stirling engine using the wobble-yoke mechanism. A prototype engine (DASE1) was developed by Fung [12] in Power Engine and Clean Energy Lab., National Cheng Kung University, based on the same mechanism. The photograph of DASE1 is given in Figure 2.
To the knowledge of the authors, prior to submitting the engine design for fabrication, computational optimization could help reduce the cost and time of manufacturing significantly, especially for an engine of higher power. Besides, due to the continuing improvements in computation technology, computational optimization has become an important tool to design the engines. For example, Ahmadi et al. [13] implemented the genetic algorithm for the optimization of a solar dish Stirling engine. This approach is a stochastic method which covers a global search of optimal geometrical parameters. However, it is less accurate and consumes more computation time. Hooshang et al. [14] describes the optimization of a gamma-type Stirling engine using neural networks as a pattern recognition tool to construct a mathematical model. This method requires a large amount of data for the neural network model to be accurate. On the other hand, Zare and Tavakolpour-Saleh [15] investigated the influence of dominant poles' places on a free piston Stirling engine by means of the particle swarm optimization scheme. affected by the pressures exerted on the upper and the lower faces of the piston, hence the name 'double-acting'. One working zone of the double-acting Stirling engine connects the compression chamber of one cylinder through the flow ducts to the cooler, regenerator, heater and, finally, the expansion chamber of another cylinder. If this engine configuration is compared to four conventional α-type Stirling engines, the numbers of cylinders and pistons are reduced by half. This greatly increases the power density of the engine. With four cylinders, the phase angle of oscillation between two neighboring pistons is specified as 90°.
In 1977, Clucas and Raine [6] proposed a wobble-yoke mechanism that was used to drive the moving parts of the four-cylinder double-acting α-type Stirling engine. Later, García-Canseco et al. [7] proposed a model for control of the wobble-yoke Stirling engine. Chatterton and Pennacchi [8] have conducted a study on different multi-cylinder engine configurations and their thermodynamic performances. Urieli and Berchowitz [9] presented an ideal isothermal thermodynamic model and Yu et al. [10] presented a non-ideal adiabatic thermodynamic model to predict the performance of Stirling engines. Tan and Cheng [11] performed three-dimensional computational fluid dynamic analysis of the four-cylinder, double-acting Stirling engine using the wobble-yoke mechanism. A prototype engine (DASE1) was developed by Fung [12] in Power Engine and Clean Energy Lab., National Cheng Kung University, based on the same mechanism. The photograph of DASE1 is given in Figure 2.
To the knowledge of the authors, prior to submitting the engine design for fabrication, computational optimization could help reduce the cost and time of manufacturing significantly, especially for an engine of higher power. Besides, due to the continuing improvements in computation technology, computational optimization has become an important tool to design the engines. For example, Ahmadi et al. [13] implemented the genetic algorithm for the optimization of a solar dish Stirling engine. This approach is a stochastic method which covers a global search of optimal geometrical parameters. However, it is less accurate and consumes more computation time. Hooshang et al. [14] describes the optimization of a gamma-type Stirling engine using neural networks as a pattern recognition tool to construct a mathematical model. This method requires a large amount of data for the neural network model to be accurate. On the other hand, Zare and Tavakolpour-Saleh [15] investigated the influence of dominant poles' places on a free piston Stirling engine by means of the particle swarm optimization scheme .    The SCGM method is an optimization algorithm proposed by Cheng and Chang [16], which stemmed from the traditional conjugate-gradient method (CGM). The SCGM method simplifies the CGM method by giving a fixed step size value, and meanwhile the sensitivity analysis is done by adding perturbations to the designed parameters. In this manner, the relationship between the objective function and designed parameters can be determined without overwhelming the users with mathematical derivations. This method is well known and has been widely applied in the optimization of various engineering devices, such as fuel cell [17], thermoelectric cooler [18], and micro reformers [19].
The above literature survey has revealed the necessity of optimization of the Stirling engines. However, currently there exist only a few reports which are relevant to the optimization of the Stirling engine, not to mention the four-cylinder double-acting α-type Stirling engine. Under these circumstances, it is essential to develop an efficient numerical tool for designing this type of engine. In this study, the non-ideal adiabatic thermodynamic model [10] is modified and then incorporated with the SCGM method [16] so as to build a numerical optimization tool for this purpose. It is important to note that the application of the present approach includes but is not limited to the fourcylinder, double-acting α-type Stirling engine. The optimization method developed in this study is general and could be applied for the optimal design of any Stirling engine, provided that a direct problem solver like a thermodynamic or CFD model [12,20] has been developed. Detailed information about the present thermodynamic model and the optimization method is described in the subsequent sections.

Thermodynamic Model
In the thermodynamic model, only one working zone is considered, since these four cylinders will exhibit a periodic thermal and fluidic pattern from one another. In a typical working zone shown The SCGM method is an optimization algorithm proposed by Cheng and Chang [16], which stemmed from the traditional conjugate-gradient method (CGM). The SCGM method simplifies the CGM method by giving a fixed step size value, and meanwhile the sensitivity analysis is done by adding perturbations to the designed parameters. In this manner, the relationship between the objective function and designed parameters can be determined without overwhelming the users with mathematical derivations. This method is well known and has been widely applied in the optimization of various engineering devices, such as fuel cell [17], thermoelectric cooler [18], and micro reformers [19].
The above literature survey has revealed the necessity of optimization of the Stirling engines. However, currently there exist only a few reports which are relevant to the optimization of the Stirling engine, not to mention the four-cylinder double-acting α-type Stirling engine. Under these circumstances, it is essential to develop an efficient numerical tool for designing this type of engine. In this study, the non-ideal adiabatic thermodynamic model [10] is modified and then incorporated with the SCGM method [16] so as to build a numerical optimization tool for this purpose. It is important to note that the application of the present approach includes but is not limited to the four-cylinder, double-acting α-type Stirling engine. The optimization method developed in this study is general and could be applied for the optimal design of any Stirling engine, provided that a direct problem solver like a thermodynamic or CFD model [12,20] has been developed. Detailed information about the present thermodynamic model and the optimization method is described in the subsequent sections.

Thermodynamic Model
In the thermodynamic model, only one working zone is considered, since these four cylinders will exhibit a periodic thermal and fluidic pattern from one another. In a typical working zone shown in Figure 1, the solution domain is divided into five sub-chambers, which are the expansion chamber, heater, regenerator, cooler and compression chamber. The piston displacement can be determined from the geometrical parameters of the wobble-yoke mechanism illustrated in Figure 3. Let the angles between the bottom surface of the wobble yoke arms and the horizontal plane be represented as θ e and θ c , for the expansion and the compression chambers, respectively, and the shaft rotation angles for the expansion and compression chambers be φ e and φ c . The values of φ e and φ c can be determined in terms of θ e and θ c as where φ e = φ c + π 2 .
Energies 2020, 13, x FOR PEER REVIEW 4 of 19 for the expansion and compression chambers be e φ and c φ . The values of e φ and c φ can be determined in terms of e θ and c θ as Knowing , e c θ , one can determine the piston displacement as  Knowing θ e,c , one can determine the piston displacement as where ζ = L 1 /L 2 . Then, the lengths of expansion and compression chambers can be calculated by taking point O as the origin.
The schematic of one working zone of the engine is plotted with the geometrical parameters in Figure 4. The volume of each chamber can then be calculated based on the lengths obtained above as where V TDS and V BDS represent the dead space volumes, and ε is the regenerator's porosity which is determined from the wire mesh geometrical parameters as in which d rm is the diameter of the mesh wire and λ is the wire pitch.
Energies 2020, 13, x FOR PEER REVIEW 5 of 19 The schematic of one working zone of the engine is plotted with the geometrical parameters in Figure 4. The volume of each chamber can then be calculated based on the lengths obtained above as  Initial values for the pressure in all chambers are assumed to be equal and denoted as P 0 . The initial temperature of the expansion chamber and heater is assumed to be equal to the heater wall temperature, which is treated as the hot end temperature. The initial temperature of the compression chamber and cooler is assumed to be equal to the cooler wall temperature, which is the cold end temperature. Regarding the regenerator, multiple nodes are adopted to better predict the temperature gradient in the regenerator, as suggested in references [21,22]. The initial streamwise temperature profile in the regenerator is assumed to be a linear profile between the hot and the cold end temperatures. The use of multiple nodes can lead to a final streamwise temperature profile in the regenerator. Besides, the initial fluid temperatures in all the chambers are given as where e, h, r, k, c indicate the expansion chamber, heater, regenerator, cooler and compression chamber, respectively; n r is the number of nodes in the regenerator, and, in this study, typically n r is assigned to 40 after a node-number-independence check; and R gas is the gas constant of the working fluid, helium. The pressure variation in the working zone and the mass variations in all the chambers can be derived from the differentials of mass conservation equation, energy conservation equations and the ideal-gas equation of state. The results are given below where T he and T ck denote the temperatures at the interface at heater/expansion chamber and interface at compression chamber/cooler, respectively, which are given with upstream conditions based on the flow direction. The temperature variation in each chamber can be determined with the energy conservation equation as where δQ denotes the heat transfer between the walls and the working fluid in the control volume, dH the enthalpy change of working fluid, δW the output work by the working fluid, and dE the change in internal energy which equals mC v dT. For the regenerator, the temperature of the solid matrix is determined from the energy balance as The thermal resistances for heat transfer between the working fluid and the walls in the heater and the cooler are calculated in terms of the Nusselt numbers as where Nu h and Nu k represent the Nusselt numbers on the surfaces of the heater and the cooler, respectively, which are calculated with the correlations provided in reference [23] for the flow in the circular-tube annulus, with a constant surface temperature on one surface and the other insulated. Similarly, the thermal resistance between the working fluid and the solid matrix in the regenerator is where the Nusselt number on the heat transfer area between the working fluid and the solid matrix is calculated by the correlations suggested by Rohsenow and Cho [24]. The temperature in the engine is varied from 300 to 1000 K. Meanwhile, the pressure in the engine is remarkably changed in a cycle, particularly in a highly pressurized engine. Therefore, the dependence of the transport properties (thermal conductivity k th and dynamic viscosity µ) of the working fluid on temperature and pressure is taken into account in the present model. The transport properties are expressed as functions of temperature and pressure based on empirical equations presented by Kuehl [25], as below. The coefficients in the following equation are tabulated in Table 1. The working volume density is calculated using the ideal gas equation of the state. The constant pressure and constant volume specific heat capacity are assumed as constant values. Table 1. Coefficients for the evaluation of helium's dynamic viscosity and thermal conductivity. [25].

Dynamic Viscosity Thermal Conductivity
Based on the obtained volumes and fluid temperatures in all the chambers, the average pressure can be calculated by using the ideal-gas equation of state as Next, the average pressure is further corrected by considering pressure drops in each chamber. The friction coefficient in each chamber can be expressed as a function of Reynolds number. It is also noted that the influence of temperature has been taken into account in the dynamic viscosity and density terms in the calculation of Reynolds number. The empirical equations are taken from the existing report [26]. Moreover, since a real engine has connecting ducts of different sizes, additional pressure losses are expected due to a sudden expansion or contraction of the cross-sectional area of the ducts. For predictions of the pressure losses due to abrupt area change, the friction coefficients as functions of area ratios are obtained from references [27,28]. With the friction coefficients known, the pressure losses in all the chambers as well as the pressure losses due to abrupt area changes are calculated. Then, these pressure drops are introduced into the average pressure for the determination of the pressures in the expansion and the compression chambers as The computation is started from the initial conditions described in the precedent section and is advanced to a stable regime eventually. After the stable regime is reached, then the power output and the thermal efficiency of the engine are calculated.
The test case of the thermodynamic model is specified with the prototype engine DASE1, whose parameters are listed in Table 1. Note that the charged mass of the working fluid, helium, in one working zone, is 63.4 mg to yield a cycle average pressure of approximately 4 bar, which is calculated using the ideal-gas equation of state by considering the charged pressure, initial volume and initial temperature in the engine chambers.
Next, a parametric analysis is performed to select the designed parameters. The effects of each parameter are evaluated using the thermodynamic model. In this section, the size of the engine and the charged mass of the working fluid are fixed. In addition, as a parameter is investigated, all other parameters are fixed at the given values displayed in Table 2. The results of the parametric analysis are plotted in Figure 5.  Figure 5a shows the effects of the heater outer diameter on the engine indicated power output and thermal efficiency. The magnitude of thermal efficiency increases with the outer diameter and finally reaches a maximum of 33.97% at D h 2 = 0.071 m. On the other hand, the curve for the indicated power output exhibits a peak value at D h 2 = 0.0683 m. It is clear that the peak values of the two curves appear at different heater outer diameters. A peak exists with a curve because when the heater outer diameter is small, an increase in the heater outer diameter leads to an increase in the heat transfer surface area at the top of the cylinder and in the heater annular flow duct volume as well. Over the peak, the performance of the engine descends due to the excessive dead volume with larger heater outer diameter.
The effects of the regenerator outer diameter on the engine performance are depicted in Figure 5b. Again, one observes two different regenerator outer diameters: one corresponds to highest power output and another one to the highest thermal efficiency. This is because a certain amount of wire mesh is required to carry out heat transfer between the solid matrix and the working fluid. However, as D r continues to grow, the dead volume and the pressure losses due to abrupt area changes also increase.
The effects of the wire diameter and wire pitch of the wire mesh used for making the regenerator are shown in Figure 5c,d, respectively. A change in wire diameter or the wire pitch of the regenerator affects the regenerator's porosity. The porosity then affects the dead volume of the working chamber and the heat transfer between the fluid and solid matrix of the regenerator. In Figure 5c, it is seen that there is an optimal wire diameter for the maximum indicated power output, whereas the curve for thermal efficiency does not exhibit a peak value. In Figure 5d, an optimal wire pitch for maximum indicated power output is found at λ = 0.15 mm, and again the curve for thermal efficiency does not exhibit a peak value.
Finally, the effects of the cooler length on the engine performance are conveyed in Figure 5e. Since the height of the engine is fixed, a change in cooler length is accompanied by a change in the heater length. In this figure, an optimal value of the cooler length exists for the highest indicated power output, while the thermal efficiency increases with the cooler length monotonically.

Thermodynamic Model vs CFD Model
On the other hand, CFD models can lead to more detailed and complete information of the flow and thermal fields, which is not easily obtained with the thermodynamic model. However, one major concern about the CFD analysis is its time consumption in the computation. To estimate the time consumption in the CFD analysis of the present Stirling engine, a CFD computation has also been performed for the test case. Detailed information regarding the numerical schemes used was described in detail by Cheng and coauthors [12,20]. Figure 6 shows the comparison in indicated power output as a function of engine speed between the thermodynamic model and the CFD analysis. Close agreement between the two sets of data is clearly seen in this figure. Note that the readings of the power output in Figure 6 are delivered by only a working zone. With a four-cylinder double-acting engine, the total power output is four times the reading given in this figure. This implies that the engine can produce a total power output of approximately 1.2 kW at 2000 rpm.
heater length. In this figure, an optimal value of the cooler length exists for the highest indicated power output, while the thermal efficiency increases with the cooler length monotonically.

Thermodynamic Model vs CFD Model
On the other hand, CFD models can lead to more detailed and complete information of the flow and thermal fields, which is not easily obtained with the thermodynamic model. However, one major concern about the CFD analysis is its time consumption in the computation. To estimate the time consumption in the CFD analysis of the present Stirling engine, a CFD computation has also been performed for the test case. Detailed information regarding the numerical schemes used was described in detail by Cheng and coauthors [12,20]. Figure 6 shows the comparison in indicated power output as a function of engine speed between the thermodynamic model and the CFD analysis. Close agreement between the two sets of data is clearly seen in this figure. Note that the readings of the power output in Figure 6 are delivered by only a working zone. With a four-cylinder double-acting engine, the total power output is four times the reading given in this figure. This implies that the engine can produce a total power output of approximately 1.2 kW at 2000 rpm.
A comparison between the temperature from the thermodynamic model and the volumeaverage temperature from the CFD analysis may be helpful. Figure 7 shows that there exists a small discrepancy between the two sets of temperature data. Nonetheless, the trend of the curves is nearly the same and the maximum difference in temperature is acceptable. A comparison between the temperature from the thermodynamic model and the volume-average temperature from the CFD analysis may be helpful. Figure 7 shows that there exists a small discrepancy between the two sets of temperature data. Nonetheless, the trend of the curves is nearly the same and the maximum difference in temperature is acceptable.
A comparison in computation time between the two models is shown in Table 3. The computations are executed on a personal computer with Intel Xeon CPU E5-2620 v3 processor. It is seen that, for a typical test case, the CFD model needs more than 52 h to finish 20 cycles in simulation, whereas the thermodynamic model requires only 1 min. In the optimization process, the number of iterations could exceed 200, and hence it is not practical to choose the CFD model as the direct problem solver. According to the comparison described above, the validity of the developed thermodynamic model is basically confirmed, and therefore the thermodynamic model is chosen as the direct problem solver for the optimization, instead of the CFD model. A comparison in computation time between the two models is shown in Table 3. The computations are executed on a personal computer with Intel Xeon CPU E5-2620 v3 processor. It is seen that, for a typical test case, the CFD model needs more than 52 h to finish 20 cycles in simulation, whereas the thermodynamic model requires only 1 min. In the optimization process, the number of iterations could exceed 200, and hence it is not practical to choose the CFD model as the direct problem solver. According to the comparison described above, the validity of the developed thermodynamic model is basically confirmed, and therefore the thermodynamic model is chosen as the direct problem solver for the optimization, instead of the CFD model.

Optimization Method
This study intended to develop an efficient approach for optimizing the engine by integrating the developed thermodynamic model and the SCGM method [16]. As the baseline design of this engine is readily available, local optimization is carried out to improve the engine performance while keeping the original engine configuration. Thus, the SCGM method is chosen for its deterministic approach and fast convergence.
In addition, the approach is modified for a multi-goal optimization, in which the objective function is defined such that the indicated power output and the thermal efficiency can both be increased at the same time. For this purpose, the objective function is defined as the reciprocal of the sum of indicated power output and thermal efficiency in the following

Optimization Method
This study intended to develop an efficient approach for optimizing the engine by integrating the developed thermodynamic model and the SCGM method [16]. As the baseline design of this engine is readily available, local optimization is carried out to improve the engine performance while keeping the original engine configuration. Thus, the SCGM method is chosen for its deterministic approach and fast convergence.
In addition, the approach is modified for a multi-goal optimization, in which the objective function is defined such that the indicated power output and the thermal efficiency can both be increased at the same time. For this purpose, the objective function is defined as the reciprocal of the sum of indicated power output and thermal efficiency in the following where in the right-hand side, the denominator of the fraction includes two terms. The magnitude of the indicated power output is represented by the first term, and the magnitude of thermal efficiency second term. The two positive weighting coefficients, c 1 and c 2 , are prescribed by the users based on the design purpose. Note that c 1 + c 2 = 1.0. The effects of different combinations of c 1 and c 2 are investigated later. In this study, a set of k designed parameters are optimized simultaneously. Let {Xj|j =1 , 2, 3, . . . , k} denote the designed parameters to optimize. The minimization of the objective function will be achieved in the iterative process. By means of the SCGM method [16], the set of the designed parameters are given initial values by the specifications of the prototype engine, and then the optimal values will be yielded in iterations. The thermodynamic model is adopted as a direct problem solver to predict the indicated power output and the thermal efficiency of the engine, corresponding to the updated designed parameters. The direct problem solutions are sent back to the SCGM code in order to calculate the objective function by Equation (33) and then to update the designed parameters. The designed parameters are updated by where the step size (β s j ), and the search direction (Π s j ) are calculated at every iteration with the following equations For the first iteration, the step size β 1 j is given a constant value and the conjugate gradient coefficient Λ 1 j is assigned to be zero. With the variable step size calculated from Equation (35), the optimization iteration can speed up as the gradient of search direction is large and slow down as the search approaches the optimal design.
The conjugate gradient coefficient Λ s j is calculated by where the gradient functions (∂J j /∂X j , j = 1, 2, 3, . . . , k) are determined by a direct numerical sensitivity analysis: First, a perturbation (∆X j ) is added to each of the designed variables, and then the change in the objective function (∆J j ) caused by ∆X j is determined. The gradient function with each designed variable is approximated by numerical differentiation: ∂J j /∂X j = ∆J j /∆X j . The optimization process is terminated as the objective function reaches its minimum value. Therefore, the process is completed as the relative error of objective functions between two consecutive iterations is decreased to be less than 10 −9 . The coupling of the optimization method with the thermodynamic model is shown in Figure 8.

Results and Discussion
In the parametric analysis shown in Figure 5, only one parameter is changed, while the others are fixed in each plot. In fact, multi-parameter optimization is more practical from a designer's viewpoint. The five parameters considered in Figure 5 exhibit peaks on the performance curves of the engine, and hence, as a test of the present approach, they are all selected to be the designed parameters to optimize simultaneously. Their original values are listed in Table 4.

Results and Discussion
In the parametric analysis shown in Figure 5, only one parameter is changed, while the others are fixed in each plot. In fact, multi-parameter optimization is more practical from a designer's viewpoint. The five parameters considered in Figure 5 exhibit peaks on the performance curves of the engine, and hence, as a test of the present approach, they are all selected to be the designed parameters to optimize simultaneously. Their original values are listed in Table 4. The optimization is carried out under some geometrical constraints. For example, the outer diameters of heater and cooler cannot exceed a critical value, 0.08 m. The sum of the lengths of the heater and cooler is fixed at 0.144 m. The wire diameter of the wire mesh used in the regenerator should always be less than the wire pitch. The regenerator's porosity, which can be expressed in terms of wire diameter and wire pitch, should always be between 0 and 1. These constraints are tabulated in Table 5. Table 5. Geometrical constraints.

Parameter Constraint
Heater outer diameter D h 2 < 0.08 m Regenerator outer diameter D r < 0.08 m Heater and cooler lengths L h 2 + L k 1 = 0.144 m Wire diameter and wire pitch d rm < λ Regenerator's porosity It is important to mention that in the present approach, all the parameters of the engine are categorized into designed, fixed and dependent parameters. The relationship between the dependent and the designed parameters must be clearly defined and expressed as the constraint equations shown in Table 5. If, for example, parameter a is actually another form of parameter b, one simply needs to figure out the constraint (relationship) between a and b. Therefore, other than the highly independent parameters, this model works for inter-dependent parameters by incorporating constraints among the parameters.
Firstly, the magnitudes of the weighting coefficients are set to be c 1 = 0.3 and c 2 = 0.7. The optimal search for the combination of the design parameters is in the direction of decreasing objective function. Figure 9 shows the variation in the objective function in the optimization process. It is seen that the value of the objective function indeed is continuously decreased. For this case, the iteration number to reach the optimal designs is less than 250, and the lowest value of the objective function is approximately 0.0067.
The comparison of P-V diagrams between the original and optimal designs is illustrated in Figure 10. It is seen that the areas of the cycles on P-V diagrams in the both the expansion and the compression chambers are indeed enlarged remarkably. The expansion in the areas of the cycles reveals the improvement in the power output of the engine by optimization. The results of the optimal set of the designed parameters and the engine performance before and after the optimization are tabulated in Table 6. Using the present optimization approach, the indicated power output can be increased from 1062.56 to 1659.72 W, and the thermal efficiency from 27.41% to 37.22%. function. Figure 9 shows the variation in the objective function in the optimization process. It is seen that the value of the objective function indeed is continuously decreased. For this case, the iteration number to reach the optimal designs is less than 250, and the lowest value of the objective function is approximately 0.0067. The comparison of P-V diagrams between the original and optimal designs is illustrated in Figure 10. It is seen that the areas of the cycles on P-V diagrams in the both the expansion and the compression chambers are indeed enlarged remarkably. The expansion in the areas of the cycles reveals the improvement in the power output of the engine by optimization. The results of the optimal set of the designed parameters and the engine performance before and after the optimization are tabulated in Table 6. Using the present optimization approach, the indicated power output can be increased from 1062.56 to 1659.72 W, and the thermal efficiency from 27.41% to 37.22%.  It is seen in Table 6 that the heater outer diameter tends to increase for a larger heat transfer surface area, but it is limited by the increase in dead volume. The regenerator outer diameter approaches the outer diameter of both heater and cooler so that the pressure drop due to the sudden expansion and contraction of the flow is minimized. This can be seen also in the literature [29], where  It is seen in Table 6 that the heater outer diameter tends to increase for a larger heat transfer surface area, but it is limited by the increase in dead volume. The regenerator outer diameter approaches the outer diameter of both heater and cooler so that the pressure drop due to the sudden expansion and contraction of the flow is minimized. This can be seen also in the literature [29], where the pressure drop due to the sudden expansion and contraction of flow is minimized when the heater, regenerator and cooler diameters are almost the same. The wire mesh parameters are optimized so that this leads to the optimum porosity. The cooler length increases, while the heater length reduces, because the thermal resistance at the cold side is higher than that of the hot side.
Influence of changing the weighting coefficients c 1 and c 2 on the optimal design's performance is evaluated and the results are displayed in Table 7. The main purpose of the optimization program is to maximize the indicated power output without reducing the thermal efficiency. Three different combinations of c 1 and c 2 are considered here. As stated earlier, the coefficients may be selected depending on the design purposes. By elevating the value of c 1 , the emphasis on the power output improvement is stronger, so that the power output is expected to become higher. As shown in this table, as the value of c 1 is increased from 0.15 to 0.5, the indicated power output is elevated from 1615.56 to 1663.12 W, while the thermal efficiency is slightly reduced from 37.54% to 36.37%. The results provided in Table 7 agree with this expectation. With this approach, the initial guess can be varied arbitrarily to a certain extent. It is important to know whether the optimization method leads to a unique optimal point as the initial guess is changed. Thus, three additional initial guess sets are adopted in addition to the previous case (initial guess 1). Figure 11 shows the trajectories of the optimization process from the four sets of initial guesses in the space of (D h 2 , D r , L k 1 ). The values of the parameters associated with the four initial guesses and the optimal point are tabulated in Table 8. It is clearly seen that, even though the initial guesses are separated in the space, the approach always finds the same optimal point. It means that the approach is robust and independent of the initial guess. Note that the approach is not just limited to the present group of designed parameters. When necessary, more designed parameters may be added in the group. Using more designed parameters actually means that the optimization can be done in a bigger search space so that a better solution can be expected. If the group of designed parameters are changed, the optimization requires a separate convergence study. The approach presented in this report is of great potential for the practical design of the Stirling engine.
It is important to note that the present optimization method is capable of modifying an existing original design so as to reach a higher performance. However, the side effect of manufacturing, mechanical properties may be possible. One possible manufacturing problem would be that the optimized parameters of the parts may not be commercially available. If that is the case, to reduce the cost, one can pick commercialized parts close enough to the optimized geometrical parameters. It is important to note that the present optimization method is capable of modifying an existing original design so as to reach a higher performance. However, the side effect of manufacturing, mechanical properties may be possible. One possible manufacturing problem would be that the optimized parameters of the parts may not be commercially available. If that is the case, to reduce the cost, one can pick commercialized parts close enough to the optimized geometrical parameters.

Conclusions
In this study, an optimization approach based on thermodynamic modelling incorporated with the SCGM method has been successfully developed to optimize the geometrical parameters of a four-cylinder double-acting α-type Stirling engine using the wobble-yoke mechanism.
The thermodynamic model is built by modifying an existing non-ideal adiabatic thermodynamic model. The validity of the developed thermodynamic model has been confirmed with a three-dimensional computational fluid dynamics analysis. The approach can continuously adjust the combination of the designed parameters in order to minimize the objective function.
After optimization, the indicated power output can be increased from 1062.56 to 1659.72 W, and the thermal efficiency from 27.41% to 37.22%.
Finally, the robustness of the optimization method is tested using different initial guesses. The results show that the approach is robust and independent of the initial guess for a typical case.
In conclusion, the approach presented in this report is of great potential for the practical design of the Stirling engine.