Exchanging data between computational fluid dynamic and thermodynamic models for improving numerical analysis of Stirling engines

In a numerical analysis of a Stirling engine, the thermal and flow fields in the cylinder are three‐dimensional and periodically varying. Therefore, a computational fluid dynamic (CFD) analysis may provide detailed solution but it takes a long computation time to simulate enough number of cycles toward the final stable oscillatory regime. On the contrary, a thermodynamic model is relatively time‐effective; however, it is not sufficiently accurate because it is built under some ideal assumptions. In this study, a method of exchanging data between the two models is proposed to facilitate the simulation of Stirling engines. In this method, the CFD model provides some coefficients required by the thermodynamic model. The simulation is then advanced to the stable oscillatory regime by the thermodynamic model. The information of the stable oscillatory regime predicted by the thermodynamic model is introduced back to the CFD model to renew its initial guesses. Then, the CFD model corrects the required coefficients for the thermodynamic model. The iteration of prediction‐and‐correction continues until the exchanged information converges. Since the predicted initial guesses of the CFD analysis by the thermodynamic model are actually very close to the solutions in the stable oscillatory regime, the CFD model can achieve detailed solutions for the stable oscillatory regime in just a finite number of cycles. In this manner, the overall time consumption can be dramatically reduced and accuracy of the CFD analysis can be improved as well.

biomass energy, and geothermal energy is environmentally friendly and may offer solutions to the aforementioned energy crisis. Using more efficient devices may significantly reduce waste energy and save available energy. Stirling engines have received great attention of the researchers since they could be used as alternative power devices that solve the energy crisis. In particular, Stirling engines are efficient, environmentally friendly, and compatible with a variety of external heat sources, including biomass energy, solar energy, and waste heat. 2,3 It is widely recognized that theoretical models for simulation of the Stirling engine can be categorized into three levels based on the complexity or the number of physical assumptions with these models. 4,5 Schmidt 6 first proposed a primitive first-level theoretical model for predicting the Stirling engine performance in 1871. This first-level model was characterized by five ideal assumptions, including idealgas working fluid, uniform pressure throughout the engine, isothermal hot and cold chambers, no leakage of working gas, and sinusoidal motion of the displacer and the piston. These ideal assumptions make Schmidt's model overestimate engine performance by about 30%. 7 Second-level models are the upgraded versions of the first-level model by either modifying some ideal assumptions to build up more realistic models or correcting models through pressure drops, hysteresis losses, regenerator losses, radiation, convection losses, and so on. Urieli and Berchowitz 8 introduced an adiabatic model in which adiabatic boundary conditions were prescribed on surfaces of the expansion and the compression chambers, while temperatures of the working fluid in other chambers are assumed to be uniform and fixed. Parlak et al 9 built a quasisteady model stemming from the model of 8 by considering dissipation and convection losses. Cheng and coworkers 10,11 modified the adiabatic condition by distinguishing temperatures of thermal sources and working gas and reckoning pressure loss through the regenerator. Udeh et al 12 investigated a β-type Stirling engine by using a non-ideal second-level model dealing with the shuttle heat losses and mass leakage between the expansion chamber and the compression chamber and between the workspace and the crankcase.
The CFD models are referred to as the third-level models in which major assumptions with the second-level models are removed in the expense of huge computational resources caused by highly geometric complexity, nonlinearly mathematical models, and coupling between the thermal and flow fields. Nonetheless, these models improve the accuracy of predictions and give more detailed information of the thermal and flow fields. 4 For example, Cheng and Chen 13 performed a three-dimensional CFD simulation of a β-type Stirling engine. Predictions of the thermal efficiency and indicated power, as functions of the rotation speed, had been verified with experiments. El-Ghafour et al 14 carried out the CFD simulation of the Stirling engine model GPU3. Average differences between the numerical and experimental results were found to be within 4%. Authors found that the temperature profile in the regenerator is a logarithmic profile in streamwise direction and observed significant oscillation in the temperature of the regenerator.
The thermal and flow fields in the cylinders of the Stirling engines are periodically varying. As the rotation speed of the engine is 1000 rpm, the engine will go through 1000 cycles in 60 s to reach stable operation from the cold start. As a matter of fact, the temperature of the solid matrix and the working fluid in the regenerator requires considerable time to develop as compared to other chambers. Therefore, the CFD analysis may take a huge computation time to cover enough number of cycles to reach the final stable oscillatory regime, particularly in the regenerator. According to Cheng and Tan, 15 for a particular Stirling engine analysis, the CFD model required more than 52 hours to finish only 20 cycles in simulation, whereas the thermodynamic model needed no longer than 1 minute. The computation time with CFD analysis can be extremely long to cover more than 1000 cycles. On the other hand, a lower-level thermodynamic model is relatively time-effective. However, the solutions of the thermodynamic model are not sufficiently detailed because it is built under some ideal assumptions. Besides, if some required coefficients, such as thermal resistances and friction factors in the working chambers, are lacking, the model cannot yield accurate predictions of engine performance.
In this study, a method of exchanging data between the 3rd-level CFD and the 2nd-level thermodynamic models is proposed to facilitate the simulation of Stirling engines. In this method, the CFD model, with initial guesses of the thermal and fluid conditions in the chambers, provides the coefficients required by the thermodynamic model. The simulation is then rapidly advanced to the stable oscillatory regime by using the thermodynamic model. The information of the stable oscillatory regime predicted by the thermodynamic model is introduced back to the CFD model to renew its initial guesses. In general, the computation can lead to final stable solutions in just a finite number of iterations of data exchanging. It will be shown later that the overall time consumption can be dramatically reduced and accuracy of the CFD analysis can be improved as well.

| 2179
CHENG aNd PHUNG with existing experimental data for validation of the numerical approach. The physical model of the Stirling engine is plotted in Figure 1.

| CFD model
In this section, mass, momentum, and energy transport equations are formulated with the CFD model. The realizable k -turbulence model including equations of turbulent kinetic energy and dissipative rate of turbulent kinetic energy is adapted to determine the properties of turbulent flow. The governing equations are briefly listed as follows: Continuity equation: Momentum equation: where Reynolds stresses R can be expressed in terms of mean velocity components u i , eddy viscosity t , and turbulent kinetic energy k as and the viscous stresses can be expressed as where is the laminar viscosity. Energy transport equation: The transport properties, such as laminar viscosity and thermal conductivity of the working gas, are dependent on both temperature and pressure, as described by Kuehl. 17 Ideal-gas equation of state: Turbulent kinetic energy k: Dissipation rate of turbulent kinetic energy : Eddy viscosity t : The solid matrix in the regenerator causes momentum losses. The momentum equation (2) is modified by adding a Darcy-Forchheimer sink term for the porous-medium flow in the regenerator Meanwhile, local thermal-equilibrium assumption is introduced for the porous-medium flow in the regenerator. Thus, the energy equation is written as Velocities of the piston and the displacer in the Stirling engine are prescribed as After solving the above set of governing equations for the solutions of three-dimensional, transient velocity, pressure, and temperature distributions in the cylinder, one can determine the indicated power and the thermal efficiency by.
where P c and P e are volume-averaged pressures associated with volumes of compression and expansion chambers, respectively, and Q • in is the input heat transfer rate on the surfaces of heater tubes and expansion chamber.
The initial guesses and boundary conditions for the test cases are summarized in Table 1.

| Thermodynamic model
With the thermodynamic model, there are several assumptions made: 1. Temperature and pressure of the working fluid are transiently varying but spatially uniform in each subchamber. 2. The changes in kinetic and the potential energies are negligible. 3. Specific heats of the working gas and the solid material in the regenerator are assumed constant. 4. The thermal field in the regenerator is assumed to be in local thermal equilibrium. 5. The working gas obeys the ideal-gas equation of state.
It is important to mention that the thermodynamic model needs to be modified so that its accuracy can match the numerical predictions from the CFD model. These modifications includes the following: 1. Global and local mass conservation equations are applied in determination of the mass flow rates so that some fundamental assumptions used in most of the existing models, for example, constant temperatures in heat exchangers and adiabatic conditions on surfaces of the expansion and the compression chambers, are eliminated. 2. The effects of pressure losses are introduced into the energy equations for all the chambers so that the indicated power and heat transfer rates are able to meet the requirement for energy balance. 3. To capture spatial variations in temperature and pressure, the engine is divided into five main chambers, namely, compression chamber, cooler fins, regenerator, heater tubes, and expansion chamber. Among the five main chambers, the cooler fins, regenerator, and heater tubes are further divided into a number of subchambers.
The main chambers and subchambers are shown in Figure 2. Note that in total, there are n subchambers. The 1st one is the expansion chamber, and the nth one is the compression chamber as illustrated in Figure 2. The positive direction of the mass flow rates is chosen to be from the expansion chamber to the compression chamber. Note that volumes of subchambers in each main chamber are equally divided.
The local and global mass equations, the energy equation, and the ideal-gas equation of state are expressed for the kth subchamber as Local mass conservation equations: Global mass conservation equation:

Energy conservation equation:
where i w is used to specify the thermal boundary conditions. If the wall is isothermal, i w = 1, and if the wall is adiabatic, i w = 0. Once the pressure losses are evaluated, they are introduced into the pressure term on the right-hand side of the above equation. Ideal-gas equation of state: where R is the gas constant. Note that quantities with superscript "*" in Equations (15)(16)(17)(18) are exact solutions to the ordinary differential equations. Quantities without superscript "*" are approximate solutions to the discretized equations in the subsequent sections. Equations (15) and (17) for the kth subchamber can be discretized using the implicit Euler method In the existing papers, 7,19-24 uniform pressure throughout the workspace is determined from a combination of the global mass conservation equation and the ideal-gas equations of state. Then, the pressure difference between two consecutive time levels is evaluated based on two assumptions: (1) constant temperatures in heat exchangers and (2) adiabatic conditions on surfaces of the expansion and the compression chambers. Next, the mass difference between two consecutive time levels is calculated based on the obtained pressure difference. 7,22,25,26 However, since the two assumptions are impractical, a new method for determining the mass difference between two consecutive time levels is attempted. In this study, the mass difference between two consecutive time levels is directly determined from the known values of masses at the current and latest time levels. Thus, the two assumptions are no longer needed.
Accordingly, Equations (16) and (18) are used for determination of the uniform pressure and gas mass in each subchamber as Then, mass difference between two consecutive time levels is.
Mass flow rates are determined by the local mass conservation, Equation (19), as mentioned in. 7,25,26 where m • 1 = 0 and m • n + 1 = 0. The gas flow direction, or in other words, the mass flow rate, is directly evaluated from the mass difference between two consecutive time levels as mentioned in the Equations (23) and (24). In addition, pressure losses in the workspace are caused by several factors, for example, the friction between the flow and the stationary walls, the sudden change in the cross-sectional area, and the flow through the regenerator, etc. 11,26 Pressure losses cause nonuniform distribution of absolute pressure and degrade engine performance. In this study, the pressure losses are expressed by where f k is the friction factor; L k is the longitudinal length of the kth subchamber, and D h,k is the hydraulic diameter. Besides, averaged densities are calculated with k = m k ∕V k and velocities are calculated The friction factors are determined as follows: where the Reynolds number is Re k = k v k D h,k ∕ k and viscosity k = T k , P k is calculated based on the transport properties equation presented by Kuehl. 17 Relative pressures in the central and two neighboring subchambers inside the regenerator are To avoid the instability in the solution procedure, the nonuniform pressure should satisfy the global mass conservation by adjusting the pressure of the central subchamber of the regenerator P k reg .
Since (19) Re k for compression, expansion chamber, and heater tubes CHENG aNd PHUNG one has The nonuniform absolute pressure in all subchambers can be computed by The enthalpy flow rates can be determined based on upstream conditions as (20) yields As mentioned in, 26,27 the power, estimated from the nonuniform pressure, is not directly involved in the energy balance equation so the power and the heat transfer rates do not satisfy the energy balance. In this study, the existence of the nonuniform pressure in Equation (34) guarantees the energy balance between the indicated power and the heat transfer rates for the final stable oscillatory regime. The n simultaneous equations described by Equation (33) can be solved by the Gauss-Seidel method 28 for temperature solutions.
The indicated power and heat transfer rates are updated with The time-marching solution procedure with the thermodynamic model is terminated when changes in indicated power, rate of heat input, and rate of heat output between two consecutive cycles are all below 10 −3 W.
It is noted that in the present model, there coexist two absolute pressures with different roles in the current thermodynamic model; one is the uniform pressure P and the other is the nonuniform pressure P k . The assumption about uniform pressure P in the ideal-gas equation of state is exploited to evaluate the mass flow rates and masses of subchambers. These quantities are, in turn, used for determined the relative nonuniform pressure P rel,k and the absolute nonuniform pressure P k with the help of the global mass conservation equation. The absolute pressure P k is then directly involved in the energy equation for determining temperatures T k and the indicated power • ��� ⃗ W . In summary, the procedure for calculating P and P k is summarized as follows: It is interesting to investigate the subtle influence of the total number of the subchambers on the predictions of the engine performance. Figure 3 shows the variations in numerical predictions by the thermodynamic model as the total number of the subchambers is changed. The relative difference in indicated power and thermal efficiency between two last cases are 0.05% and 2.5%, respectively. It is found that the case having 120 subchambers in total leads to satisfactory accuracy and acceptable computation time. Therefore, it is chosen in the numerical analysis. For this case, the number of subchambers for heater tubes, regenerator, and cooler fins in the thermodynamic model is 48, 45, and 25, respectively.
It should be noted that time step sizes of both the CFD and thermodynamic models are identical with value of 2 × 10 − 4 s for the baseline case. To match between two models, it should be guaranteed that the varying volumes of expansion and compression chambers between are closely equal, as shown in Figure 4.

BET WEEN CFD AND THERMODYNAMIC MODELS
In this section, the approach of data exchanging between the CFD and the thermodynamic models is described clearly for the particular test case shown in Table 2. Note that this study is based on an assumption that the rotation speed is known and constant. As shown in Table 2, the rotation speed of the test case is set to be 1500 rpm for both the CFD model and thermodynamic model so that the displacement and velocity of piston and displacer are confirmed.
The thermodynamic model contains several coefficients that should be given beforehand, such as the thermal resistances R k in Equation (20) on the walls of the expansion and compression chamber, heater tubes, and cooler fins, and the friction coefficient C in Equation (27) for the regenerator. In the computation, the CFD model predicts the required coefficients for the thermodynamic model, and then, the latter can be used to run over 1000 cycles to readily reach the stable oscillatory regime and obtain the stable temperatures and pressures in all the subchambers. The information of the stable oscillatory regime predicted by the thermodynamic model is introduced back to the CFD model to renew its initial guesses. Since the predicted initial guesses of the CFD analysis by the thermodynamic model are actually very close to the solutions in the stable oscillatory regime, the CFD model can achieve detailed solutions for the stable oscillatory regime in just a finite number of cycles. Then, the CFD model corrects the required coefficients again for the thermodynamic model. The iteration of prediction-and-correction continues until the exchanged information converges. In general, it needs only a couple of exchanges to achieve convergence. As a result, the CFD analysis of the stable oscillatory regime can be accelerated so that overall computation time can be dramatically reduced. Besides, accuracy of the numerical analysis is improved through the matching between the two models.
The exchanging procedure between the two models is illustrated in Figure 5. First, the CFD simulation is started with initially guessed conditions and it calculates the temporal variations in flow and thermal fields and heat transfer rates in each chamber. These data with the last CFD cycle (typically, the 5th cycle) are recorded and used for matching with the thermodynamic model. Then, the unknown coefficients of the thermodynamic model are calculated based on the results of the CFD analysis in the data exchanging process. For example, the heat transfer rates can be determined from the thermal fields and the pressure drops in the workspace from the flow fields. The heat transfer rate on the walls and the volume-averaged pressures in the chambers obtained from the CFD analysis are further used to determine the thermal resistances, pressure drops, and the friction factors of the thermodynamic model. If relative differences in thermal efficiency and indicated power between the two models are within 3%, then the procedure is terminated. Otherwise, temperatures of subchambers and uniform pressure at the crank angle of 0 degree from the thermodynamic model are used to renew the initial guesses of the CFD model via a User-Defined Function (UDF) of the commercial software. This UDF setting must guarantee that the initial guessed temperatures of all the subchambers, which were extracted from the thermodynamic model, are precisely assigned to their corresponding locations in the solution domain of the CFD model such that total masses of the two models are equal. It should be noted that the equality in total mass is an essential factor to the success of the exchanging process. It is noted that the thermodynamic

| RESULTS AND DISCUSSION
In this study, the commercial CFD software, ANSYS Fluent ver. 16.2, and the in-house code of the thermodynamic model are executed on a personal computer with Intel® Xeon ® CPU E5-2620 v3@2.40 GHz. From a set of four meshes, 1 256 075, 1 750 279, 2 397 036, and 3 647 349 cells, the mesh of 2 397 036 cells is selected as a result of the gridindependence check procedure. The relative differences of indicated power and thermal efficiency between two later meshes are both below 0.75%. For the selected mesh, the minimum orthogonality is equal to 0.21 while the maximum skewness is 0.89. Additionally, the time step size of 2 × 10 − 4 s is selected because the relative differences in indicated power and thermal efficiency between two time steps, 2 × 10 − 4 and 1.667 × 10 − 4 s, are 0.06% and 0.13%. The firstorder implicit scheme is used for the temporal discretization. The coupled scheme is selected for the pressure-velocity coupling. Furthermore, the second-order upwind scheme is used for the momentum and energy transport equations. Details about the CFD model are available in Cheng and Phung. 18 Typically, the matching procedure requires only three iterations. Results of required coefficients of thermodynamic model are tabulated in Table 3 for the test case using the parameters given in Tables 1 and 2 after the data exchanging process. For this particular case, the CFD simulation gives the indicated power of 85.9 W and the thermal efficiency of 35.0% whereas the thermodynamic model gives the indicated power of 88.2 W and the thermal efficiency of 34.8% once the data exchanging process is completed. The relative differences between two models in terms of the indicated power and the thermal efficiency are within 3%. Results of the matching process between the CFD simulation and the thermodynamic model in the final stable Compression chamber 0.892

Friction factor
Friction coefficient C in Equation (27)  700 Initial guesses of CFD model (at zero crank angle) provided by the thermodynamic model Initial CFD temperatures Changing during data exchanging Initial CFD pressures Changing during data exchanging oscillatory regime are displayed in Figure 6. It is found that the P-V diagrams for the expansion and compression chambers, input and output heat transfer rates, and the volumeaveraged regenerator temperatures yielded by the two models are in good agreement through the data exchanging process. In the thermodynamic model, nonuniform absolute pressure distribution in all the subchambers can be predicted using Equation (31). Based on the pressure data of the subchambers, the pressure difference between the expansion and the compression chambers can be further calculated. The variation in the pressure difference between the expansion and the compression chambers against the crank angle, predicted by the thermodynamic model, is illustrated in Figure 7. It is noted that the maximum pressure difference between these two chambers reaches 47.66 kPa at the crank angle of 316.8°, whereas the minimum pressure difference is −16.09 kPa at the crank angle of 176.4°. The pressure difference between the two chambers is not large.
In order to demonstrate the relative performance of the present approach, Table 4 shows the comparison of the time consumptions of three different modes, namely, thermodynamic model only, CFD model only, and data exchanging modes. All three models are run on the same personal computer. It is seen that running over 1030 cycles to reach the stable oscillatory regime, the CFD-model-only mode requires about 227 days, whereas the data exchanging mode needs merely 80 hours. The CFD model with the help of the thermodynamic model can greatly reduce the computational time. This proves that the present approach indeed is superior to the CFD-model-only mode in terms of the computation time. One can also see that  Figure 8 are the profiles for the 10th, 25th, 50th, 100th, 250th, and 500th cycles. The end position at x∕l r = 0 corresponds to the interface between the regenerator and the heater tubes, and x∕l r = 1 corresponds to the interface between the regenerator and the cooler fins. It is shown that the temperature profile changes rapidly from the first cycle to the 100th cycle, and then gradually reaches a stable linear profile at the 500th cycle. Actually, if the temperature profile in the regenerator is not stable yet, the flow and thermal fields in the cylinder might not be stable. According to Table 4, if it is the CFD model used for simulation, the estimated time would be over 110 days to complete the computation of 500 cycles. This is why it is not practical and requires improvement.

METHODS
In this section, the experimental data, provided by Cheng and Phung, 18 are used here for validating the data exchanging approach. The specifications of the β-type Stirling prototype engine are detailed in Table 5. Meanwhile, final values of the coefficients required by the thermodynamic model, which are provided by the CFD model as described in the precedent sections, are illustrated in Table 6. Detailed information regarding the experimental and numerical methods of determination of the shaft power of the engine is available in. 18 To save space of this report, the information is not provided herein.
A comparison in the shaft power versus rotation speed between the numerical predictions and experimental data is exhibited in Figure 9. The numerical predictions are carried out for different heating temperatures by the present CFD model with the help of the thermodynamic model. It is seen that the curves of the CFD model closely agree with the experimental data. This proves that the present approach can provide accurate predictions. For example, an increase in the heating temperature raises the maximum shaft powers and also shifts the optimal rotation speed rightwards. It is also found that the curves of the thermodynamic model are matched with the experimental data to a certain extent when the rotation speed is below the maximum power point. However, the thermodynamic model underestimates the shaft power as the rotation speed exceeds the maximum power point. This is expected even though the reasons may be complicated due to the ideal assumptions made with the model.

| CONCLUSIONS
In this study, a method of exchanging data between the 3rdlevel CFD and 2nd-level thermodynamic models is proposed to facilitate the numerical simulation of Stirling engines. In this method, the CFD model provides some coefficients required by the thermodynamic model, and the latter is used to advance the simulation to the stable oscillatory regime within short time. The information of the stable oscillatory regime predicted by the thermodynamic model is introduced back to the CFD model to renew its initial guesses. In general, the CFD models need huge computation time to reach the stable oscillatory regime while the thermodynamic model requires only seconds. The overall time consumption can be dramatically reduced, and accuracy of the CFD analysis can be improved as well. It is found that the cost in CFD simulation is indeed remarkably reduced with the help of the thermodynamic model. To achieve the goal, the accuracy of the existing thermodynamic model needs to be improved so that its predictions match the numerical predictions from the CFD model to a certain extent. Therefore, a modified thermodynamic model is proposed by eliminating some assumptions about the thermal boundary conditions, and in the meantime, the pressure losses effects are taken into account in the energy equations. Meanwhile, to capture spatial variations in temperature and pressure, the engine is divided into five main chambers, and among the five main chambers, the cooler fins, regenerator, and heater tubes are further divided into a number of subchambers.
For a particular case considered in this study, it is seen that running over 1,030 cycles to reach the stable oscillatory regime, the CFD-model-only mode requires about 227 days, whereas the data exchanging mode needs merely 80 hours. The CFD model with the help of the thermodynamic model can greatly reduce the computational time.
In a comparison in the shaft power versus rotation speed between the numerical predictions and existing experimental data for different heating temperatures, the present approach leads to numerical solutions which closely agree with the experimental data.