Today's energy production is still comprised of large-scale cogeneration to fulfil the demand for heat and power in the industry. However, the power generation is rapidly starting to transform. Centeralised systems, where a single power plant distributes the electricity to household buildings, will be replaced by more flexible decentralised systems or smart grids. Such systems will produce electricity from renewable sources and also thermal power units that are independent of the main distributor. Indeed, small-scale highly efficient and flexible thermal power units are still an absolute necessity to ensure grid stability and flexibility in case of a lack of renewable energy. Therefore, there is a growing interest for more efficient and cost-effective energy systems to meet today's flexibility requirements and emission reduction targets. In this framework, small-scale power generation machines and more specifically Micro Gas Turbines (mGT) seem to show significant potential (Ismail et al., 2013). Their multi-fuel capabilities coupled with the low emissions and maintenance costs make them an ever-growing competitor against engines with similar electric output (Pepermans et al., 2005).

In this context, mGT operation characteristics should be able to change dynamically to respond to the fluctuating demand of a flexible power grid. The engine has to effectively transit from one operating point to another and also work at part load. Moreover, during transient periods, the compressor could be facing potentially dangerous operation due to degradation (Spakovszky et al., 2000). Thus, it is crucial to have at our disposal a reliable transient model that can predict accurately how the cycle parameters evolve in dynamic conditions.

In principle, mGT or generally gas turbine performance simulation is firmly relying on a comprehensive understanding of the engine components' behaviour. More specifically, the accurate performance prediction of each turbomachinery part at off-design states is presented through the component maps. In mGT cycles not only the compressor but also the turbine maps provide an adequate insight into the operating range as the turbine is usually operating below choking conditions contrary to larger GTs.

Commonly, the performance maps are provided by the original equipment manufacturer. Additionally, maps could be calculated using streamline curvature and computational fluid dynamics (CFD) or performing time consuming experimental studies (Zanger et al., 2011; Kislat et al., 2019). On the other hand, detailed experimental tests can be potentially hazardous to the engine performance as the compressor is forced to work near the surge line.

The look-up table is one of the oldest methods to display the component's performance characteristics (Walsh and Fletcher, 2004). If a point is not located at the iso-speed line, it is interpolated linearly or with a spline. As the performance maps present a significant nonlinearity between the parameters of reduced speed, mass-flow and pressure ratio, the look-up tables' precision is strongly linked with the amount of experimental data. Thus, the interpolation of limited experimental data jeopardises the quality of prediction.

This precise representation of component maps is a subject that attracted many researchers throughout the years. For this reason, Kong et al. proposed the method of shifting and scaling the shape of the performance map by matching the actual component's features (Kong et al. 2003). In addition, Tsoutsanis et al. suggested nonlinear component map fitting methods to build analytical relations between the characteristic parameters and replaced the look-up tables (Tsoutsanis et al., 2015). Their method is proven to display satisfying results in transient gas turbine simulations.

Apart from the correct simulation of the turbomachinery part, mGT dynamic models should be able to calculate the key parameters that characterise the recuperator in real time. Many studies have been performed to capture the energy exchange dynamic phenomena in heat exchangers integrated in Brayton cycles. Camporeale et al. presented an analytical mathematic model by applying the partial differential equations expressing unsteady heat transfer (Camporeale et al., 2000). Ghigliazza et al. developed and validated a numerical model for the recuperator of the Turbec T100 engine within the framework of the TRANSEO software package (Traverso, 2005). For the dynamic computation of the energy equation, he discretised the heat exchanger into N cells. He also assumed a simplified “lumped volume” approach for the continuity and momentum equations (Ghigliazza et al., 2009). Furthermore, Henke et al. pursued an approach quite similar to the one of Ghigliazza et al. when he simulated the performance of Turbec T100 recuperator with a different treatment on the calculation of heat transfer coefficients (Henke et al., 2016).

The above-mentioned modelling techniques present without a doubt several advantages in favour of the accurate performance prediction. However, the link between precision, computational time and overall complexity remain to be examined especially in the field of mGTs. Therefore, we plan to use them in a control application to account for a real time performance.

In this paper, we present the development of an accurate and efficient dynamic mGT model. First, we discuss a component map modelling method for the compressor and the turbine to enhance the prediction of the mGT performance in a transient simulation. Various analytical relations are tested to fit in the iso-speed lines and the curve with the lower deviation is used in the model. Furthermore, a modelling technique for the heat exchanger is developed and compared regarding its adjustment with the measurements in the T100 test rig. Finally, experimental data of the recuperator are effectively used to tune the heat transfer coefficient of the recuperator model with the intention to increase even more the accuracy of the model.


In this section, the calculation strategies for the compressor, turbine and recuperator of the mGT dynamic model are analysed. First of all, the dynamic model, that is employed in this study, is described along with the thermodynamic cycle and the nominal conditions of the machine. Consequently, the map fitting method is presented for the compressor and turbine and finally, the heat exchanger model is explained in detail.

mGT dynamic model

The mGT dynamic model developed for and described in this paper, aims at simulating the behaviour of the Turbec T100 engine which is located at Vrije Universiteit Brussel (VUB). This mGT, exploiting the recuperated Brayton cycle can operate in dry and wet modes as it was converted into a micro humid air turbine (De Paepe et al., 2014). For the purposes of this paper, we are only concerned with the dry mode simulation. The T100 has a nominal power output of 100 kWe, 165 kWth with a maximal rotational speed close to 70,000 rpm at nominal operating conditions. It is equipped with a control system that keeps the produced electric power constant by adjusting the rotational speed. Moreover, Turbine Outlet Temperature (TOT) is also controlled and fixed at around 645°C by managing the fuel flow in the combustion chamber to maintain high thermal efficiency.

The simulation tool is coded in MATLAB/SIMULINK environment. It offers the ability to combine different blocks and user-defined functions into a subsystem as it is depicted in Figure 1. In all cycle components, we pursued a 0-D modelling with a lumped parameter approach to account for the energy and mass conservation dynamics (Traverso, 2005). The control system consists of the power control which changes the produced power to achieve the desired speed and the fuel control that keeps a constant TOT. The scheme and parameters of the control tool have been received from the manufacturer and are based on a similar approach as described in Montero Carrero et al. (Montero Montero Carrero et al., 2015). For the combustion chamber, we assume a constant volume inside it and a homogeneous air/natural gas mixture. Moreover, a plenum is placed after the instant combustion to account for the pressure and mass flow dynamics.

Figure 1.

Model layout of the mGT cycle.

Turbomachinery map fitting

The purpose of this process is to fit the accessible map data with a specific mathematical equation for both turbomachinery (compressor and turbine) component performance maps. For instance, the fitting equation for each compressor map iso-speed line should have a single form that is controlled by the coefficients of the equation. The map data of the T100 mGT, used for this study, is provided by the manufacturer. Therefore, the precision of our fitting techniques relies on the deviation from the data and the form of the mathematical equation. For this reason, 3 fitting cases will be analysed for both compressor and turbine maps and tested and the technique that minimizes the error from the data will be integrated into the dynamic model.

Compressor map

For the representation of pressure ratio πc in relation to the reduced mass flow m˙c rate and shaft speed N, three fitting curves are tested with different complexities. The equations of the three cases are presented below.


Equation 1 is an ellipse curve with a center at (0,0) and axes rotation. The Superellipse is a curve with the same parameters as the Ellipse and a coefficient n which alters the shape and adds an extra degree of freedom. For equation 3, a Supershape curve is used with a center at (x0,0), 5 degrees of freedom and no axes rotation. The fitting of the several speed lines for the three cases is presented in Figure 2. It is easily observed that Ellipse and Superellipse are not able to simulate the steep decrease in pressure ratio (choking conditions). Hence, the Supershape is employed to test if we can fit the data points in that region with a single equation. It seems that this equation can fit all the data points but with questionable accuracy. Therefore, we calculated the relative Root-Mean Squared Error (RMSE) for each case in every speed line of the data points that do not belong to the choking region and compared the different methods. Figure (a) shows that Ellipse has the most effective fitting behaviour in low as well as high rotational speeds with a relative error below 0.9%. Moreover, Supershape presented increased error (above 1.5%) despite the fitting of all data points. Therefore, Case 1 is chosen to model the speed lines due to its satisfactory results among the proposed fitting methods. This technique is also coupled with a linear interpolation between the points of maximum mass flow and minimum pressure ratio at the chocking region to account for the accurate representation of the steep decrease in the πc value. Then, the parameters of the adopted equation (α,b,ϑ) are expressed as a function of rotational speed with the use of 1-D interpolation function provided by Python programming language.

Figure 2.

Compressor maps fitting approaches: Ellipse, Superellipse and Supershape.

Turbine map

A comparable approach is followed for the turbine map fitting. Although for this occasion a different group of fitting curves is applied. Due to the distinct form of the iso-speed lines in a turbine map, an exponential curve with a rotating axis is used (equation 4). The fitting behaviour for his equation is depicted in Figure 3 and is compared with a 4th and 5th degree polynomials that are the equations 5 and 6 consequently.

Figure 3.

The Exponential fitting approach shows good fitting characteristics on the data points of the turbine map.

Exponential has 4 coefficients (α,b,c,ϑ) and its fitting behaviour on the datapoints is shown in Figure 4. This equation can only be solved indirectly, contrary to the other two cases. This characteristic makes it computational heavier solution against polynomial curves. However, in Figure 4(b) we can see that equation 4 presents exceptional fitting accuracy, especially in the low-speed region. For this reason, again a hybrid scheme was adopted for the turbine map fitting. Thus, the Exponential equation is employed for the low-speed region (N<1750rpmK0.5) and Poly4 for the rest to ensure a balance between accuracy and fast calculations. Additionally, the coefficients of the adopted curves are demonstrated as a function of rotational speed with 1-D interpolation. Similar methods and fitting equations were used for the calculation of the efficiency curves for both the compressor and turbine (not included in this work).

Figure 4.

The relative RMSE of the 3 cases as a function of the reduced speed for (a) the compressor and (b) the turbine map fitting. The relative RMSE of the fitting curves for the compressor and turbine highlight the superiority of the first method (red lines).

1-D heat exchanger model

The current work focuses on the recuperator developed by RSAB and specifically designed for mGT applications (Lagerström and Xie, 2002). Such a 1-D model is necessary as a simpler 0-D approach assumes a constant wall temperature inside the heat exchanger. This creates numerical problems in dynamics conditions and especially in compact recuperators used for mGTs where the temperature wall gradient is high (McDonald, 2003). The heat exchanger is discretised into a N cells (5 cells used for good accuracy) according to Figure, where the energy equation is applied using the following equation:


The heat capacity and the mass of the casing are expressed with cm and m respectively. The remaining values of the equation are presented in Figure 5. For the heat balance of the air and gas stream, we applied a steady-state approximation. This is fairly justified as the temperature rate of change in each “stream” cell is greater than the temperature rate of change in each “wall” cell. So, the heats that are absorbed and released in each cell are presented below:

Figure 5.

Modelling scheme of the heat exchanger.

We also assumed that the thermal resistance of the recuperator plates is fairly small and could be ignored.

Moreover, two different cases for the heat transfer coefficient will be considered. In the first one, the UA is kept constant and is taken by the work of Lagerström at nominal conditions. In the second case, the UA is expressed as a function of m˙c. This function is constructed by tuning the UA in order for the component to match the effectiveness (ε) of experimental data of the T100 in 4 different operating points (100, 90, 80, 70 kW). Then a linear interpolation is applied between the tuned UA points which have a range between 5,117 and 5,279 W/K.

Results and discussion

The map fitting method and the recuperator model are tested in respect to their fidelity to transient experimental data taken from the T100. We considered four different steps in the demanded power output using the control system of the model in order to ascertain the accuracy of the model in various load transitions. The different steps are the following: −10 kWe, −20 kWe, +10 kWe, +20 kWe. A steady state validation was also performed with a relative mean error from experiments below 2% but due to space limitations we only opted to focus on the dynamic part.

Map fitting

The map fitting technique that is implemented and described above is compared with the approach that is proposed by Sieros et al. (called “Polynomials”) for the compressor and assumes that the turbine is chocked (Sieros et al., 1997). Both methods are considered common in turbomachinery modeling. The rotational speed is closely related to the pressure of the system by the compressor and turbine map. Moreover, the air temperature at the outlet of the compressor is linked with the compressor map through the isentropic law. Hence, Figure 6 shows the response of shaft speed (N) and compressor outlet temperature (Tc,out) and their difference from the experimental data at a load step −20 kWe from the nominal power (100 kWe). The rotational speed behaviour shows a satisfying agreement with the experimental data. The transition to lower rotational speed is modeled effectively as the simulation follows the slope of the experimental data rather well. Furthermore, the deviation from the real shaft speed at 100 and 80 kWe (in 0 and 200 s respectively) is reduced by the implementation of the “Map fitting” method as the “Polynomials” approach has a larger relative mean error (1.15% and 0.38% respectively). Similarly, Tc,out satisfactorily follows the experimental data for the demanded power step. The accurate simulation of this value is crucial as a deviation from the measured data can propagate to the recuperator component. This will prevent us from accurately validating the heat exchanger 1-D model. Therefore, as shown in Figure 6, Tc,out presents a maximum difference in nominal power of no more than 1°C. Also, the results tend to match those from the experiment at 200 s. The “Polynomials” model, in this case, shows very close values with our method. This is explained as the temperature is more dependent on the connected “lumped volume” in dynamic conditions than the pressure ratio. Moreover, the variation in the steady-state value of Tc,out at nominal power is caused by the assumption of constant pressure losses inside the compressor.

Figure 6.

The behaviour of Rotational speed and Temperature in the compressor outlet for a demanded electrical power output step from 100 to 80 kWe.

Recuperator model

As it was described in the previous section, the heat exchanger model is initially developed with a constant heat transfer coefficient taken at nominal conditions. This approach was tested and it seemed more appropriate to use a variable UA in order to capture the temperature dynamics inside the heat exchanger with precision. Therefore, experimental data from the T100 were utilized and we calculated the effectiveness (ε) of the component in four different operating modes (70, 80, 90 and 100 kWe). Then the UA is tuned until the model simulates the effectiveness that was observed from experiments in those operating modes.

Figure 7 shows the behaviour of the recuperator outlet temperature (Tcc,in) with negative steps in the demanded electrical power. The transition from 100 to 90 kWe (Figure 7a) presents a small deviation in the steady-state values from the experimental data. This absolute difference is approximately 1°C at 500 s and is a product of the constant UA. In dynamic operation the specific heat flux (q˙) along with the temperature difference between the wall and the fluid (ΔT) change. As a result, UA does not remain constant from a physical point of view. Thus, the 1-D model slightly miscalculates the Tcc,in during the transition to 90 kWe. We easily overcome this problem by using a variable UA being a function of the mass flow of the cold medium and tuning it to simulate the ε that was observed from experiments in each operating point. After the introduction of a variable UA the Tcc,in follows satisfactorily the experimental data in the two load steps of Figure 7. Moreover, the transition curve in 50–100 s is faithfully followed in both Figure 7a and b by the model.

Figure 7.

The calculated recuperator outlet temperature (cold side) follows well the experimental data for a step in demanded electrical power output after tuning from (a) 100 to 90 kWe and (b) 100 to 80 kWe.

Figure 8a shows that after the 1-D model is tuned, the deviation for the experiments is minimized in a positive load step-change as well. This is observed in Figure 8b. Although, in this case, the model cannot simulate efficiently the large drop in temperature at the time period of 50–150 s. This problem is not truly related to the features of the 1-D numerical model. This steep decrease in the temperature of the experimental data is caused by the drop in TOT and it is the product of the control system of the test rig. During the experimental tests, the demanded power is stabilised first and afterwards the TOT is manipulated to remain constant. A positive step-change in the demanded power leads to an increase in the mass flow rate of the engine. As a result, the mixture in the combustion chamber becomes leaner and the Turbine Inlet Temperature (TIT) decreases along with the TOT. This drop in TOT cannot be properly calculated by the mGT dynamic model as the model's control system tries to immediately stabilize the TOT at 645°C. Therefore, this behaviour can be fixed by adjusting the characteristics of the control system model. As it is depicted in Figure 8c the TOT of the measurements drops substantially compared with the simulated TOT. The effect of this phenomenon in the Tcc,in is observed when we use the TOT from the measurements as an input in the model (Figure 8b red lines). As the TOT drops, the simulated Tcc,in is affected as well and gets closer to the experimental data. Thus, it is clear that the model simulates the recuperator performance with good precision. The effect of the control system on the design parameters is not part of the current work due to space limitations. However, it seems crucial to focus on the behaviour of the control in future work.

Figure 8.

Recuperator outlet temperature (cold side) for a step in demanded electrical power output from (a) 90 to 100 kWe and (b) 80 to 100 kWe.


Simulation strategies for the correct modelling of mGT components dynamics are developed. The map fitting method for the correct modelling of the compressor and turbine is thoroughly described. Additionally, the discretisation scheme together with the implemented equations of the 1-D numerical model for the heat exchanger have been discussed in detail. All the presented modelling methods have been integrated into the MATLAB/SIMULINK dynamic model of the Turbec T100 to increase the accuracy of the calculations.

The results of the T100 mGT model, which are focused on the compressor and turbine maps as well as the recuperator, were compared with measurement data for the VUB T100 test rig. The map fitting method increased the overall accuracy of the model and provided an effective technique for mGT turbomachinery map modelling. The calculated rotational speed followed more accurately the experimental data after the implementation of the map fitting approach. In addition, the 1-D heat exchanger model shows promising results and accurate agreement with the experimental data at both positive and negative load step changes. The variable heat transfer coefficient increased the accuracy of the model regarding Tcc,in. This demonstrated the impact of the UA on the overall performance of compact recuperators utilized for mGT operations and should be always considered in dynamic modelling. An interesting question is how the key performance parameters of the mGT affect the surge margin of the compressor and move the operating point of the turbine towards choking. This aspect needs further research in future work and comparison with CFD results. Also, the role of the control system in such machines should be thoroughly analysed.



specific isobar heat capacity (J/kg/k)


casing heat capacity (J/kg/k)


mass (kg)


mass flow (kg/s)


shaft speed (RPM)


pressure (Pa)


heat flow or work (W)


time (s)


temperature (K)


heat transfer coefficient (W/K)

Greek Symbols


angle (radians)


pressure ratio (−)


effectiveness (−)



Computational Fluid Dynamics


Micro Gas Turbine


Rekuperator Svenska AB


Turbine Inlet Temperature


Turbine Outlet Temperature