TOWARDS REAL TIME TRANSIENT MGT PERFORMANCE ASSESSMENT: EFFECTIVE PREDICTION USING ACCURATE COMPONENT MODELLING TECHNIQUES

In an energy mix driven by renewables there is a need for small-scale highly efficient and flexible cogeneration units, such as Micro Gas Turbines (mGTs). These mGTs should perform transient operations and work at part-load to meet the power grid requirements. Therefore, full transient characterisation is necessary. One of the most crucial factors is to accurately incorporate the dynamic behaviour of each component. Compressor and Turbine performance maps, although essential, are usually obtained in costly test rigs or CFD simulations. Also, the accurate modelling of heat exchanger affects the efficiency of the whole cycle. The aim of this work is the development of real time transient mGT model, where we focus in the first step on accurate component modelling. Hence, an effective performance map representation method for mGT's compressor and turbine was developed. Moreover, a Recuperator 1-D numerical model is developed. Those modelling techniques were tested in a MATLAB/SIMULINK model in transient conditions. The fundamental target of this study is to enhance the fidelity of dynamic simulations for small scale gas turbines. Key parameters like shaft speed and combustor inlet temperature, show deviation from experiments less


INTRODUCTION
Today's energy production is still comprised with large-scale cogeneration to fulfil the demand of heat and power in the industry. However, the power generation is rapidly starting to tranform. 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 from a 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 the 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 evergrowing 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 effectivelly 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 accuratelly how the cycle parameters evolve in dynamic conditions. In principle, mGT or generally gas turbine performance simulation is firmly relying on the comprehensive understanding of the engine components behavior. More specifically, the accurate perfomance 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 of 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 stream line 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, 1998). 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 preformance 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 appoach quite similar to the one of Ghigliazza et al. when he simulated the performance of Turbec T100 recuperator with different treatment on the calculation of heat transfer coefficients (Henke et al., 2017).
The above mentioned modeling techniques present without a doubt a number of advantages in favor of the accurate performance prediction. However, the link between precision, computational time and overall complexity still 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 modeling method for the compressor and the turbine to enhance the prediction of the mGT performance in 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 modeling 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.

METHODOLOGY
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 70000 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 o C by managing the fuel flow in the combustion chamber to maintain a 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 is 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.

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 in the dynamic model.

Compressor Map
For the representation of pressure ratio in relation to the reduced mass flow ̇ rate and shaft speed N, three fitting curves are tested with different complexities. The equations of the three cases are presented below.
: |̇c os − sin | 2 + |̇s in + cos | 2 = 1 (1) : |̇c os − | + |̇s in + cos | = 1 (2) : 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 the equation 3, a Supershape curve is used with a center at ( 0 , 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 a questionable accuracy. Therefore, we calculated the relative Root-Mean Squared Error (RMSE) for each case in every speed line of the Figure 1 Model layout of the mGT cycle 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 value. Then, the parameters of the adopted equation ( , , ) are expressed as a function of rotational speed with the use of 1-D interpolation function provided by Python programming language.

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 rotated axis is used (equation 4). The fitting behaviour for his equation is depicted in Figure 3 and is compared with a 4 th and 5 th degree polynomials that are the equation 5 and 6 consequently.
: ̇= • 4 + • 3 + • 2 + • + : ̇= • 5 + • 4 + • 3 + • 2 + • + Exponential has 4 coefficients ( , , , ) 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, Exponential equation is employed for the low-speed region ( < 1750 0.5 ) and Poly4 for the rest in order 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 method and fitting equations were used for the calculation of the efficiency curves for both the compressor and turbine (not included in this work). Figure 3 The Exponential fitting approach shows good fitting characteristics on the data points of the turbine map.

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: ̇, +1 =̇, ( +1 − ) = [ , − 1 2 ⁄ ( +1 + )] (9) ̇ℎ , +1 =̇ℎ ,ℎ ( +1 ℎ − ℎ ) = ℎ [ 1 2 ⁄ ( +1 ℎ + ℎ ) − , ] (10) 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 is kept constant and is taken by the work of Lagerström at nominal conditions. In the second case, the is expressed as a function of ̇. This function is constructed by tuning the in order for the component to match the effectiveness ( ) of experimental data of the T100 in 4 different operating points (100, 90, 80, 70kW). Then a linear interpolation is applied between the tuned points which have a range between 5117 W/K and 5279 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: -10kWe, -20kWe, +10kWe, +20kWe. 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 ( , ) and their difference from the experimental data at a load step -20kWe from the nominal power (100kWe). 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 100kWe and 80kWe (in 0 sec and 200 sec 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, , 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, , present a maximum difference in nominal power no more than 1 o C. Also, the results tend to match those from the experiment at 200 sec. The "Polynomials" model in this case shows very close values with our method. This is explained as the temperature is more depended on the connected "lumped volume" in dynamic conditions than the pressure ratio. Moreover, the variation in the steady-state value of , 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 100kWe to 80kWe.

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 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,100kWe). Then the 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 ( , ) with negative steps in the demanded electrical power. The transition from 100kWe to 90kWe (Fig. 7 (a)) presents a small deviation in the steady state values 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) 100kWe to 90kWe and b) 100kWe to 80kWe.
from the experimental data. This absolute difference is approximately 1 o C at 500 sec and is a product of the constant . In dynamic operation the specific heat flux (̇) along with the temperature difference between the wall and the fluid ( ) change. As a result, does not remain constant from a physical point of view. Thus, the 1-D model slightly miscalculates the , during the transition to 90kWe. We easily overcome this problem by using a variable being a function of the mass flow of the cold medium and tune it to simulate the that was observed from experiments in each operating point. After the introduction of a variable the , follows satisfactorily the experimental data in the two load steps of Figure  7. Moreover, the transition curve in 50-100 sec is faithfully followed in both Figure 7 (a) and (b) by the model. Figure 8 (a) 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 8 (b). Although, in this case the model cannot simulate efficiently the large drop in temperature at the time period 50-150 sec. This problem is not truly related with the features of the 1-D numerical model. This steep decrease in the temperature of the experimental data is caused by the drop in 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 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 ( ) decreases along with the . This drop in cannot be properly calculated by the mGT dynamic model as the model's control system tries to immediately stabilise the at 645 o C. Therefore, this behaviour can be fixed by adjusting the characteristics of the control system model. As it is depicted in Figure 8 (c) the of the measurements drops substantially compared with the simulated . The effect of this phenomenon in the , is observed when we use the from the measurements as an input in the model ( Fig. 8 (b) red lines). As the drops, the simulated , 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 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 a future work.

CONCLUSIONS
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 in the MATLAB/SIMULINK dynamic model of the Turbec T100 in order 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 maps modelling. The calculated rotational speed followed more accurately the experimental data after the implementations 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 change. The variable heat transfer coefficient increased the accuracy of the model regarding , . This demonstrated the impact of the in 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.