Issue 
Mechanics & Industry
Volume 18, Number 5, 2017



Article Number  509  
Number of page(s)  7  
DOI  https://doi.org/10.1051/meca/2017029  
Published online  01 November 2017 
Regular Article
βNTF reduction and fast kriging simulation of optimal engine configurations
^{1}
ESTIA Recherche,
64210
Bidart, France
^{2}
Université de Bordeaux, I2M CNRS, UMR 5295,
33607
Bordeaux, France
^{3}
UVHC, LAMIH CNRS UMR 8201, Le Mont Houy,
59313
Valenciennes cedex 9, France
^{4}
UVHC, ENSIAME, Le Mont Houy,
59313
Valenciennes cedex 9, France
^{5}
AKIRA Technologies, ZA Saint Frédéric rue de la Galupe,
64100
Bayonne, France
^{*} email: s.cagin@estia.fr
Received:
16
February
2017
Accepted:
18
August
2017
In an optimization process, models are applied to simulate different design behaviors in order to determine the most suitable one. However, this requires the use of a structured methodology to correctly explore the design space and truly converge to the best solution. It is therefore necessary to test and validate the optimal design. For engines, two ways are essentially used: building and testing a real cylinder, or simulating the new design with ComputationalFluidDynamics (CFD) models. These two techniques are both expensive and time consuming. An alternative way is proposed to test new designs with a fast simulation based on a kriging method. The exploration of the design space is based on 27 cylinder configurations and the results of their CFD models. It converged to an optimal design depending on the objective function. A kriging method was used to interpolate the behavior of the optimal design just found. In this paper we present the βNTF model reduction (to define the data set used by the kriging method) and the principle of the kriging technique. We then briefly discuss the results. The results underline the method's advantages despite the small gap between the expected results and those for kriging.
Key words: kriging / fast simulation / βNTF reduction / design space / 2stroke engine optimization
© AFM, EDP Sciences 2017
1 Introduction
As engineering systems become increasingly sophisticated, the complexity of the problems faced by engineers also increases. The lack of efficient physical models is one consequence of this complexity. Empirical studies must be developed to understand the behavior of systems, to enhance models and allow for design optimization. This is particularly true for engine design, where typically, the engine is controlled by a large number of parameters in order to meet multiple performance objectives [1]. It is complicated to take all of them into account in an optimization process.
The lack of good models also has an impact on optimization processes because new designs need to be tested in order to define their performances. In the field of engines, engineers usually have two ways to validate the optimal design found: to run tests with a prototype or to analyze the results of CFD simulations. Both solutions are expensive so a third option is proposed in this paper based on a kriging method.
The kriging method is already used in various fields, including the field of engines, and for multiobjective optimization issues. In particular, kriging is fully integrated into the design space exploration phase using Response Surface Models. Multiple cylinder designs have already been improved using this method, such as Jeong et al. [2], Castric [3] or Brahmi et al. [1]. In the design optimization process presented in this paper, kriging is used in a different way: it is not part of the exploration of the design space, but it used after this phase. Indeed, after finding the optimal cylinder design, the behavior of gas flows inside the optimal chamber is interpolated from previous results thanks to the kriging technique and the gas flow distribution can be visualized. This visualization of the interpolated flow behavior is called the “Fast Kriging Simulation”.
To validate the method, it has been applied to the optimization of a 2stroke Diesel engine with ports. The complexity of the engine cycle led us to focus on scavenging. Less studied than combustion, the scavenging process also plays a major role in the formation of pollutants. Indeed, it is during scavenging that fresh air enters the cylinder pushing out the burnt gases. As with all fluid flows, any change in the cylinder design has consequences on fluid dynamics. Thus, improving cylinder design cannot be done without studying scavenging behavior.
After defining the cylinder variables influencing the scavenging process [4], a new model was developed: it explicitly integrates these variables. Based on this model, the design space was explored converging to different design solutions depending on the form of the fitness function. For the designer, the last step in design optimization is to finally test the different designs which have been defined.
First, the methodology we developed is introduced. Then, we present the βNTF and kriging methods. βNTF is used to establish the data set needed for kriging interpolation. Finally, the numerical results are given followed by a discussion underlining the interest of the method.
2 Overall methodology and case of application
Kriging visualization is a part of the overall design optimization method we developed. It is used as the final step in our methodology. The method was applied to enhance the cylinder design of a 2stroke Diesel engine. The particularity of our engine is the use of ports for both the intake and emission of gases. This creates a singular flow in the cylinder during the scavenging process which justifies us focusing on this part of the engine cycle.
The whole methodology is represented in Figure 1 and is briefly summarized in this section.
The first step is to define the design parameters influencing the flow and their range of values [4]. Then 27 cylinder configurations were simulated thanks to CFD models. The evolution of distribution of gas in the chamber during the scavenging process for each configuration was extracted [5]: these data will be used for kriging approximation as the known points. Based on CFD results, certain scavenging models were developed and used to determine the most suitable cylinder configuration.
After finding an optimal cylinder configuration, the designer needs to evaluate its new performances in order to validate the new design. Experimental tests on prototypes or numerical simulations are the most common ways to do so. However, they are both expensive: CFD is generally timeconsuming; the experiments are expensive unique models and require special instruments (e.g. transparent cylinder). Before investing time and money in simulations or experiments, the designer should be certain to have the best design. In this approach, we have developed a transitional solution based on the kriging method. It enables in a very short time for the distribution of gases to be visualized during the scavenging process whatever the cylinder design. Before going further, the designer can do a qualitative analysis and evaluate the relevance of the optimal solution. This is the purpose of the last step of the methodology presented in this article.
Fig. 1 Method of model reduction and optimization. 
3 NTF reduced model and Kriging
The kriging method implies the use of a data set. Twentyseven engine configurations have been simulated, the simulations last 174 crankshaft degrees. Based on CFD simulations and results, several quantities (fields of pressure, temperature, speed, etc.) have been extracted at each crankshaft angle degree. This represents over 30 million data items. Not all of them can be used with the kriging method; it is too much data to deal with. A model reduction was carried out using the Nonnegative Tensor Factorization (NTF) algorithm to prepare the data for kriging exploitation.
3.1 The βNTF reduction
NTF was first proposed by Shashua et al. [7]. This algorithm is a generalization of Nonnegative Matrix factorization (NMF) presented by Lee et al. [8]. NTF is based on a Parallel Factor Analysis (PARAFAC) analysis. Its particularity is that it imposes nonnegative constraints on tensor and factor matrices. PARAFAC is a tensor (multiway array) factorization method which allows us to find hidden factors (component matrices) from multidimensional data.
We consider as a 3D tensor, . We have no negative value y_{ijk} ≥ 0 ∀ i, j, k. The NTF concept is illustrated in Figure 2. (1)
With the complete model, the reduced model and the error or the noise, the decomposition generates a set of three unknown matrices U^{(n)} such as , with n the number of modes and l = (I, J, K). The NTF algorithm is very attractive because of its ability to take into account spatial and temporal correlations between variables more accurately than 2D NMF. As is the case with NMF, NTF also provides greater stability and a unique solution, as well as meaningful latent (hidden) components or features with physical or physiological meaning and interpretation [6]. Finally, the NTF algorithm is very simple to implement in python and provides powerful implementation with multiarray data. The NTF method also usually provides sparse common factors or hidden (latent) components.
In this study, in addition to its efficiency and speed, the capability of the NTF algorithm to deal with spatial and temporal correlations between variables is highly useful. Indeed, the correlations between variables have to be identified depending on numerous parameters, such as spatial localization in the cylinder, temporal moment during the engine cycle and the cylinder's design configuration. Moreover, the spatial, temporal and parametric dimensions are totally independent of each other.
Fig. 2 NTF concept. 
3.2 Application to the data set
The data are collected from the CFD results of the 27 cylinders tested (defined by an experiment design). The data extracted from the CFD results characterize gas distribution, temperature, pressure, etc. Only gas distribution will be presented here, but the method has been applied to all the fields extracted from the CFD simulations.
Because of the cylinder's piston motion during scavenging, some cells are added or deleted between two time steps. That is why all the results have been projected onto the same regular mesh (with 6561 square cells). Each cell is associated with its spatial location in the cylinder (from the regular mesh), the configuration of the cylinder (i.e. the value of each design parameter), the crankshaft angle at which the distribution is represented and the mass fraction of fresh gases (the output we want to forecast). The scavenging lasts 171 crankshaft angles: the results are collected from θ = 95 crank degrees to θ = 266 crank degrees after Top Dead Center (TDC). The piston motion in the cylinder is considered between 2 extreme positions: the TDC represents the highest position, whereas the Bottom Dead Center (BDC) is the lowest, as illustrated in Figure 3.
Due to the number of cylinder configurations, the duration of the scavenging process and the number of mesh cells, the number of known points is equal to 30 823 578. The data set therefore had to be reduced: this was done thanks to the βNTF method.
To use the NTF algorithm, the data should be organized so as to obtain the model. The three dimensions were respectively associated with space, time and input parameter combinations (cf. Fig. 4)
The NTF decomposition was done with 160 modes, the best compromise between efficiency and accuracy. The average relative error is 2.5% considering all the data. The average error reaches 8% at 180 crankshaft degrees (when the ports are fully opened and the gas exchanges are greatest).
The maximum error of 19.6% is represented in Figure 5. In this figure we can observe the gas distribution in a cylinder before and after βNTF reduction.
The data set now consists of 1 081 920 points, which represents a reduction of 96.5% compared to the CFD data thanks to the NTF method.
Now the NTF model is established. The 1 081 920 points included in the model are associated with the two spatial coordinates, the temporal moment during the engine cycle and the seven parameters of the cylinder design. All of this constitutes the data set of known points necessary to apply the kriging method. Without the reduction, the number of points is too high to apply the kriging method: the memory is full and cannot deal with all the data points, and solve the kriging algorithm at the same time.
Fig. 3 TDC and BDC piston positions. 
Fig. 4 NTF model of engine quantity. 
Fig. 5 Example of fresh gas distribution from (a) CFD results (b) βNTF reduced model. 
3.3 Kriging principles
Kriging is an optimal interpolation technique based on regression against observed Z values of surrounding data points, weighted according to spatial covariance values. The method was first developed by Krige [9] for the spatial interpolation of soil properties, but Matheron [10] was the first person to formalize the method. The kriging principle is to interpolate the value of an unknown function Z at a located point, given the values of function Z at other points. The interest of the method lies in its capability to take into account not only the distance from the evaluated point to the known points, but also the distance between two known points.
Many kriging algorithms have been developed since 1973 [10], but they are all based on the same basic linear regression estimator Z*(u) defined as [11]: (2)
With u the point where the Z function is estimated and u_{α} the neighborhood points; n(u) is the number of local points used to estimate Z*(u); m(u) and m(u) are the expected values of Z respectively for points u and u_{α} and λ_{α}(u) the kriging weights associated to u_{α} estimate the point u.
m(u) is the trend component associated with the random field Z(u). The aim of the kriging method is to evaluate and to lower the residual R(u) = Z(u) – m(u) as the weighted sum of the surrounding residuals calculated for the u_{α} points. Therefore, the key point of this technique is the way the kriging weights are evaluated: they are derived from covariance function or semivariogram.
3.4 The semivariogram
The semivariogram is a function representing the spatial dependency of the known points and is obtained from the stationarity definition. Indeed, the stationarity hypothesis is an essential condition for using kriging method. So the variation of a data set is only dependent on the distance r between two locations where variables values are Z(u_{i} + h) and Z(u_{i}) with r = h, can be given by the following semivariogram: (3)
With N(r) ={ (i, j)such asu_{i} − u_{j} = r } is the pair number of Z(u_{i} + h) and Z(u_{i}) and is the experimental semivariogram.
The experimental semivariogram presented in equation (3) estimates the theoretical semivariogram only for a finite number of distances. It is assumed that two spatially close points should have similar values and the similarity decreases as the distance increases.
Then, the variogram is established as a model adjustment based on the semivariogram points. Various forms of variogram model are available. If the semivariogram represents the spatial autocorrelation between the known points, it does not indicate the possible directions or distances which justify the need for the model's adjustment. This is usually done by linear regression or by the least squares method.
Finally the weights λ_{α} in (2) are computed solving the following system: (4)
With (5) (6) (7) λ is the Lagrange multiplier.
The variance of the square of the standard error each point is obtained thanks to the relationship: (8) indicates the incertainity of the estimated value. The accuracy of the prediction value depends directly on the distance between the evaluated point and the known points. If we consider the estimation errors normally distributed around the true value, the probability that the true value will be in Z(u_{i}) ± s_{i} is 68%, while the probability that the true value will be in Z(u_{i}) ± 2s_{i} is 95% (Davis [12]).
4 Application to the cylinder design
4.1 The data set
The use of the kriging method is motivated by its capability to take into account all the information relating to the data points in order to determine the weights λ. Not only the spatial coordinates are compared but also the temporal parameter (crankshaft angle at which the quantity is evaluated) and the 7 variables characterizing the cylinder design [4]:

the height of the exhaust port θ_{end_exh_port} which is the crankshaft angle (before the BDC) the exhaust port is fully uncovered;

the advance of opening of intake and exhaust ports, θ_{in_advance} and θ_{exh_advance}, are respectively the crankshaft angles at which the piston begins to uncover the intake and exhaust ports;

the boost pressure P_{boost}, pressure of intake gases;

the difference between intake and exhaust pressures ΔP;

the angles β_{in} and β_{exh} are respectively the angle between the piston plane at BDC and the intake or exhaust ports (Fig. 6).
Other parameters can be chosen or can be removed depending on the designer's framework.
The study is focused on the scavenging process. During the phase, the fresh gases enter the cylinder pushing out the burnt gases. During the gas exchanges, backflow and shortcutting should be avoided, which means no burnt gases go through the intake ports (this happens when the pressure in the chamber is higher than the intake pressure) and no fresh gases go through the exhaust ports.
Fig. 6 β angle. 
4.2 Fast kriging simulation
From the database, the designer is able to evaluate any quantity in the cylinder for any configuration at any moment during scavenging: the kriging method defines the quantity depending on the spatial location in the cylinder and the value of the chamber variables at any moment. Indeed, the “distance” between the evaluated point and the known points is also linked to the cylinder design in this approach.
Coded in Matlab^{©}, the kriging method is put into practice with the optimal cylinder found in [13] thanks to genetic algorithms. Kriging is used to determine the distribution of gas in the optimal design to illustrate the method's results.
4.3 Results of kriging fast simulations
The optimal cylinder is defined by the following variable values [13]:

β_{in} = 60°;

β_{exh} = 5.88°;

θ_{end_exh_port} = 9.81 crank degrees before BDC;

θ_{in_advance} = 49.02 crank degrees before BDC;

θ_{exh_advance} = 76.43 crank degrees before BDC;

P_{boost} = 2.16 bar;

ΔP = 0.29 bar.
Concerning the execution time, kriging interpolation requires around 20 min to map the gas distribution for a given crankshaft angle. Hence, 2.5 days are needed to map the whole scavenging process for one cylinder configuration. To compare, 10 h are needed to simulate the whole process of scavenging with a 2D numerical model.
Prima facie, kriging appears quite long to model the process. Because of the previous model reductions, we would expect less time. Indeed, it is just a way to globally observe the flow in the optimal geometry in order to validate it, so the execution time should be as short as possible. However, the interpolation done with the kriging method is a loop process which requires time. For each point, the kriging algorithm determines the weights comparing spatial coordinates and cylinder design of the data set points, and then approximates the expected value for the new design. Moreover, it is still faster than to get results from CFD simulations or any experiment.
In addition, there is no need to krige the whole process: we can focus on certain significant crankshaft angles (intake opening angle, BDC position, etc.) and, with the kriging method, there is no need to solve the previous time step. Therefore, it only takes 20 min to obtain the gas distribution in the cylinder after closing the ports, for example. With a numerical model, the whole process needs to be simulated to get the same information. Considering a global quantitative validation, only few crankshaft angles are needed to validate (or not) the design and to choose if the optimal design should be fully tested.
Concerning the error, Figure 7 shows the distribution of the fresh gas for the optimal configuration previously defined. Three crankshaft angles were selected:

150 crank angles after TDC, ten crankshaft angles after we see the first fresh gases entering the cylinder;

180 crank angles after TDC (the piston is at the TDC position);

256 crank angles after TDC, all ports are closed.
These 3 moments are the most significative piston positions: opening of the ports, ports fully opened and closing of the ports. To get the 3 pictures of the gas distribution, a few minutes are needed with the kriging method, only these moments are interpolated. With the CFD simulation, it is necessary to run the whole process in order to get the data about the 3 moments, around 10 h of calculation are needed.
To evaluate the error of the kriged interpolation, the CFD model for the optimal configuration was used. Figure 7 illustrates the difference between the kriged (up) and numerical results (down).
Some differences can be observed. One striking difference is the gap between the expected opening of intake ports and the kriged opening. With the fast simulation, the fresh gas enters the cylinder 139 crank degrees after TDC whereas it is supposed to open at 131 crank degrees after TDC. There is a difference of 8 crank degrees which cannot be justified by the blackflow of burnt gases in the intake duck blocking the entry of fresh gases. Indeed, all simulations show that the backflow lasts less than 5° whatever the design. More realistic, fast kriging is an interpolation of values which implies an error in terms of estimation. In addition, the data set used to interpolate the value comes from a reduced model which also implies some reduction errors as illustrated by Figure 5. The data set points should be carefully chosen and used.
The first observation is that the global path of fluids is well represented. In both simulations, the entering fresh gases follow the cylinder wall in a loop path: gases enter on the bottom left side of the cylinder, touch the cylinder head and go out on the bottom right side. This motion is due to the inclination of the intake duct. So, the kriging method well represents the general expected flow path of gases.
In addition, the evolution of the average mass fractions of fresh gases in the cylinder was drawn comparing the CFD and kriging results. Figure 8 shows that the two curves are close except at the beginning. As explained before, there is a gap between the opening of the intake ports and the entry of the first fresh gases with the kriging method. Nevertheless, it does not seem to influence the overall proportion of fresh gases (except for the beginning). At the end of the two simulations, the averaged mass fraction is almost the same with both methods.
The fast kriging simulation provides the flow trend, as we can see in Figure 7 (its path, the mixing of gases, etc.) with a reduced cost compared to CFD simulations or prototype testing. The fast simulation therefore allows a first validation of the optimal configuration found before going further or looking for another configuration.
Fig. 7 Comparison of fresh gases distribution in the cylinder between the kriged results (up) and the CFD results (down). 
Fig. 8 Comparison of average mass fraction of fresh gases between CFD and kriging results. 
5 Conclusion
The kriging method is a good alternative option to globally reproduce the flow trend. It offers initial results at a very low cost (in terms of time and resources) compared to CFD simulations or experiments. If the data are interpolated at several well selected moments, the fast kriging simulation provides a good observation of optimal design behavior. Comparing CFD and kriging results, the flow path and the average mass fraction of fresh gases in the cylinder are close enough to validate the method.
However, use of the method is mostly restricted by the data set needed for the interpolation. Indeed, the number of known points impacts the time and the precision of the interpolation. Thanks to CFD simulations, a lot of data items were extracted (more than 30 million) which motivates model reduction with the βNTF algorithm. However, model reduction induces an error which is added to the kriging interpolation error. That is why the kriging technique is a very interesting alternative way to validate a design, but does not replace a CFD simulation or tests on real prototypes.
Nevertheless, the approach presented here is also a new way to apply the powerful kriging method in a design optimization process. Instead of using it to explore the design space and search for optimal solutions, it is used to interpolate the behavior of the system during its use and validate the optimal design early on.
References
 E.H. Brahmi, L. DenisVidal, Z. Cherfi, Simulation and optimization of engine performance using kriging model and genetic algorithm, Mathematical Methods in Reliability, Moscow, 2009 [Google Scholar]
 S. Jeong, Y. Minemura, S. Obayashi, Optimization of combustion chamber for diesel engine using kriging model, J. Fluid Sci. Technol. 1 (2006) 138–146 [CrossRef] [Google Scholar]
 S. Castric, Méthodes de recalage de modeles et application aux emissions des moteurs diesel, Université de Technologie de Compiègne, 2007 [Google Scholar]
 S. Cagin, X. Fischer, N. Bourabaa, E. Delacourt, C. Morin, D. Coutellier, A methodology for a new qualified numerical model of a 2stroke diesel engine design, The International Conference On Advances in Civil, Structural and Mechanical Engineering − CSME 2014 HongKong, 2014 [Google Scholar]
 S. Cagin, et al. Scavenging process analysis in a 2stroke engine by CFD approach for a parametric 0D model development, in: 7th International Energy, Energy and Environment Symposium, Valenciennes, 2015 [Google Scholar]
 A. Cichocki, R. Zdunek, S. Choi, R. Plemmons, S.I. Amari, Nonnegative tensor factorization using alpha and beta divergences, in: IEEE International Conference on Acoustics, Speech and Signal Processing, Honolulu, Hawaii, 2007 [Google Scholar]
 A. Shashua, T. Hazan, Nonnegative tensor factorization with applications to statistics and computer vision. ICML '05, in: Proceedings of the 22nd International Conference on Machine learning, New York, 2005, pp. 792–799 [Google Scholar]
 D.D. Lee, H.S. Seung, Learning the parts of objects by nonnegative matrix factorization, Nature 401 (1999) 788–791 [NASA ADS] [CrossRef] [PubMed] [Google Scholar]
 D. Krige, A statistical approach to some basic mine valuation problems on the witwatersrand, J. Chem. Metal. Min. Soc. South Africa 52 (1964) 119–139 [Google Scholar]
 G. Matheron, The intrinsic random functions and their applications, Adv. Appl. Probab. 5 (1973) 439–468 [CrossRef] [Google Scholar]
 P. Goovaert, Geostatistics for natural resources evaluation, Oxford Uni., 1997 [Google Scholar]
 J.C. Davis, Statistics and data analysis in geology, second ed., New York, 1986 [Google Scholar]
 S. Cagin, et al. A new reduced model of scavenging to optimize cylinder design, Simulation − Transactions of the Society for Modeling and Simulation International. 92 (2016) 507–520 [CrossRef] [Google Scholar]
Cite this article as: S. Cagin, X. Fischer, É. Delacourt, N. Bourabaa, C. Morin, D. Coutelier, B. Carré, S. Loumé, βNTF reduction and fast kriging simulation of optimal engine configurations, Mechanics & Industry 18, 509 (2017)
All Figures
Fig. 1 Method of model reduction and optimization. 

In the text 
Fig. 2 NTF concept. 

In the text 
Fig. 3 TDC and BDC piston positions. 

In the text 
Fig. 4 NTF model of engine quantity. 

In the text 
Fig. 5 Example of fresh gas distribution from (a) CFD results (b) βNTF reduced model. 

In the text 
Fig. 6 β angle. 

In the text 
Fig. 7 Comparison of fresh gases distribution in the cylinder between the kriged results (up) and the CFD results (down). 

In the text 
Fig. 8 Comparison of average mass fraction of fresh gases between CFD and kriging results. 

In the text 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.