Proposal of a Predictive Mixed Experimental-Numerical Approach for Assessing the Performance of Farm Tractor Engines Fuelled with Diesel-Biodiesel-Bioethanol Blends

: The e ﬀ ect of biofuel blends on the engine performance and emissions of agricultural machines can be extremely complex to predict even if the properties and the e ﬀ ects of the pure substances in the blends can be sourced from the literature. Indeed, on the one hand, internal combustion engines (ICEs) have a high intrinsic operational complexity; on the other hand, biofuels show antithetic e ﬀ ects on engine performance and present positive or negative interactions that are di ﬃ cult to determine a priori. This study applies the Response Surface Methodology (RSM), a numerical method typically applied in other disciplines (e.g., industrial engineering) and for other purposes (e.g., set-up of production machines), to analyse a large set of experimental data regarding the mechanical and environmental performances of an ICE used to power a farm tractor. The aim is twofold: i) to demonstrate the e ﬀ ectiveness of RSM in quantitatively assessing the e ﬀ ects of biofuels on a complex system like an ICE; ii) to supply easy-to-use correlations for the users to predict the e ﬀ ect of biofuel blends on performance and emissions of tractor engines. The methodology showed good prediction capabilities and yielded interesting outcomes. The e ﬀ ects of biofuel blends and physical fuel parameters were adopted to study the engine performance. Among all possible parameters depending on the fuel mixture, the viscosity of a fuel blend demonstrated a high statistical signiﬁcance on some system responses directly related to the engine mechanical performances. This parameter can constitute an interesting indirect estimator of the mechanical performances of an engine fuelled with such blend, while it showed poor accuracy in predicting the emissions of the ICE (NO x , CO concentration and opacity of the exhaust gases) due to a higher inﬂuence of the chemical composition of the fuel blend on these parameters; rather, the blend composition showed a much higher accuracy in the assessment of the mechanical performance of the ICE. analysis, investigation, M.B., C.C.; curation, draft preparation, visualization, M.B.; M.B.


Biodiesel and Bioethanol as Alternative Fuels
Biofuels are commonly considered one of the main paths towards a new sustainable paradigm in transportation technologies. As stated by the World Energy Council, the definition of a sustainable large-scale model in the sectors of energy production and transportation is based on three core Table 1. Diesel oil/biodiesel/bioethanol main physical and chemical properties and typical values [4,12,13].

. Fuel Properties Effects in Compression-Ignition Engines
CI engines are the most widespread power-generation units in agricultural machines and tractors [14]. The effects of the fuel properties on the combustion process are crucially relevant and must be taken into great consideration when alternative fuels are used. The viscosity has a direct effect on the fuel injection process. Since the fuel viscosity increases as the ambient temperature decreases, when performing a cold start at low ambient temperatures the resulting higher viscosity of the fuel causes poorer fuel atomization and vaporization, larger droplets, and a greater in-cylinder penetration of the fuel spray [4,[15][16][17]. High viscosity may also reduce fuel flow rates, thus affecting the optimal functioning of auxiliary equipment such as the fuel feed pump [18,19]. The fuel density increases as temperature decreases. The fuel injection process is also affected by the higher density when performing a cold start at low ambient temperatures. The cetane number (CN) describes the auto-ignition capabilities of the fuels used in a CI engine. A high CN indicates a quick ignition of the fuel which in turn leads to a more complete fuel combustion and to a quicker start-up phase of the engine. On the other hand, harmful NO x emissions are reduced when the CN decreases because of the consequent ignition delay and the lower combustion temperature [20]. Water injection into the combustion chambers of diesel engines [21] is a solution that has been proposed with the aim of curbing the combustion temperature and, therefore, the NOx emissions. The evaporation characteristics (boiling point, latent heat of evaporation) influence the spray structure and bear consequences on the air-fuel mixture preparation [4]. The flash point (FP) temperature is the minimum temperature (at a pressure of 101.3 kPa) at which the fuel ignites. Low FP values might represent issues for the safety and a proper handling of fuels [4]. High oxygen fractions are related to higher fuel-combustion temperatures and consequently to higher NO x -formation rates for diffusion flames [20].

Biodiesel and Bioethanol Use in Compression Ignition Engines
Several studies [22,23] state that biodiesel has positive effects on exhaust emissions of CI engines: carbon monoxide (CO) is lowered up to 30-50 % depending on the share of biodiesel in the blend [24], as well as hydrocarbons (HC) and particulate matter (PM). The main responsible for these positive effects is the higher oxygen content and the lower hydrogen and carbon content compared to petro-diesel [6,24]. On the other hand, a variety of studies recognize [20] that biodiesel properties bring in negative effects with respect to nitrogen oxides (NO x ) emissions, whose substantial relative increase represents one of the main cons of using biodiesel as a petro-diesel substitute [25]. Controlling NO x emissions is crucial Energies 2019, 12, 2287 4 of 45 in order to meet the EURO and TIER emission standards. Improvements in terms of NO x emissions can be obtained by varying the biodiesel substitution rate in diesel-biodiesel blends [26], as well as introducing small percentages of a third alcohol-based fuel (ethanol, methanol or fusel-oil [27]) to produce ternary diesel-blends. The use of these ternary blends in diesel engines increases the brake specific fuel consumption (BSFC) but it also decreases PM concentrations in the exhaust gases, which is the greatest concern of traditional diesel engines [28]. Moreover, the effect of biodiesel and bioethanol on the power output must be considered. It has been frequently observed that the engine power and torque increase with the increase in biodiesel percentage in the blends [29]. Such behaviours should be credited to the higher oxygen content, the higher biodiesel consumption, an advance of injection timing and a shorter ignition delay time [29]. On the other hand, when bioethanol is introduced in the fuel blend, the consequent decrease of the Lower Heating Value (LHV) is responsible for a drop in the maximum power output levels [28]. However, different indications can be found in literature showing that the benefits in terms of cold flow properties caused by the addition of bioethanol could boost the power output levels [30]. Hence, when using ternary blends, the substitution of fuels in the mixture should respect a trade-off between the de-rating of performance and the reduction of polluting emissions. For example, in Shahir et al. [28] the best results have been obtained with 5%-ethanol blends: a small increase in BSFC is measured together with sensible reductions of PM, HC, CO and CO 2 emissions [28]. However, due to the complex fuel-interactions occurring in the combustion of biodiesel-bioethanol mixtures, the effects on NO x emissions are still not clear. Some studies register a NO x reduction [31,32], others report an increase in pollutant concentration [28,33]. These controversial results highlight the need for further investigation in terms of assessment of emissions and performances by developing or applying numerical optimization methods. These tools should allow the users to predict a priori the characteristics of a fuel blend and, above all, its effect on the motor outputs (mechanical and environmental performances). In the explained terms, these methods can be regarded to as predictive methods.

Emissions and Performances Prediction: Application of the Response Surfaces Methodology
The Response Surface Methodology (RSM), as extensively described by Bezerra et al. [34], is a collection of mathematical and statistical techniques based on the fit of a polynomial equation to the experimental data. It was first developed by George E. P. Box and K. B. Wilson in 1951 [35]. The final goal of the application of RSM is to describe in a mathematically simple and intuitive way (with multilinear models only) the behaviour of a dataset in order to carry out statistical previsions. A response of interest relies on several significant variables, and the aim of the method is to model and optimise such response [36]. In the last decade, several researchers in the field of alternative fuels began applying RSM to the modelling and optimization issues. A typical engine-related RSM analysis involves engine parameters such as load, speed, fuel injection timing to reduce fuel consumption, exhaust emissions and noise, and aims to improve engine performance [37][38][39] or to study the behaviour of auxiliary organs [19]. A selection of works from the extensive review paper by Yusri et al. [36], and how to proceed with RSM application in the field of biodiesel and bioethanol, is presented hereinafter.
Bose et al. [40] have applied a D-optimal RSM to a single-cylinder diesel engine (with indirect-injection system) fuelled with diesel-ethanol mixtures. Engine load and ethanol concentrations have been used as input parameters, while BSFC and NO x , CO and HC concentrations have been set as output variables. RSM application has shown a high precision, with a R 2 in the range of 0.958 (HC) and 0.983 (BSFC). Fang et al. [41] have applied a central composite design RSM to a Isuzu 4HK1-TC direct-injection medium-duty diesel engine fuelled with diesel-ethanol mixtures. Gross indicated mean effective pressure and engine speed have been used as input parameters, while CO 2 , HC, CO and NO x concentrations have been set as output variables. Pandian et al. [42] have applied RSM to a twin-cylinder water-cooled diesel engine (with a direct-injection system) fuelled with diesel-biodiesel mixtures. Injection timing, injection pressure and nozzle tip protrusion have been

Aims and Novelty of this Study
The first aim of this study is to demonstrate the effectiveness of a mixed experimental-numerical approach, based on the RSM for data processing, in the investigation of the operation of a particularly complex system such as an internal combustion engine fuelled with different biofuels blends. Indeed, an ICE consists of several devices whose behaviour can change according to the properties of the fuel. The interactions between the many possible influencing parameters are rather hard to grasp, and it is difficult to define a single or simplified mathematical relationship between these parameters and the engine performance and emissions. A first analysis is carried out to assess the effect of the different shares of biofuel on the engine performance (torque and BSFC).
Secondly, the present study intends to investigate whether the kinematic viscosity at the temperature stated by international norms (40 • C in this case, hence not necessarily linked to the combustion-chamber conditions yet easily and quickly measured with proper equipment) could be Energies 2019, 12, 2287 6 of 45 a useful parameter for indirectly estimating some performance parameters related to the engines of agricultural machines. Overall in this study, specific attention has been devoted to the engine toque, the hourly fuel consumption, the NO x and CO concentration, and the opacity of the exhaust gases.
The relevance of this study from a practical point of view is demonstrated by the easy application of the models and the simple applicability of the methodology; indeed, few and easily-accessible input parameters are required to handle the mathematical relationships presented in this work. Also unexperienced users will be able to obtain reliable forecasts of engine performance and emissions when an ICE is fuelled with alternative fuel blends without the need to assess the detailed and complex behaviour of each single sub-component. By doing so, it will be possible to apply this method also to other fuel blends and/or to other CI motors, regardless their type (e.g., direct/indirect injection, atmospheric/turbocharged) and size (i.e., light/heavy-duty).

Experimental Setup
The experimental procedure consisted of several tests on an agricultural tractor engine and aimed to record the engine performance as well as the main polluting characteristics of the exhaust gas within the full operating-range of the engine. The tests were performed on a "T4020V" farm tractor by New Holland Agriculture (Torino, Italy). The machine was equipped with a 4-cylinder 3200-cm 3 diesel engine with direct fuel injection system (Table 2). The tractor was fuelled with five fuels: pure diesel oil, two binary and two ternary blends of diesel oil (D), biodiesel (B; Table 3) and bioethanol (E; Table 3).  Specifically, these blends were: B15 (85% D, 15% B in volume), B25 (75% D, 25% B), B15E3 (82% D, 15% B, 3% E) and B25E3 (72% D, 25% B, 3% E). An external tank was used during these tests for rapid fuel-substitution and for the measurement of the instant consumptions; a description of such this apparatus is available in [48]. The kinematic viscosity of the fuels involved were measured with a "Simple VIS" portable automated viscometer by Cannon Instrument Company (State College, PA, USA) [49]. A total of 73 measurements have been performed on the fuel samples (from 9 to 19 samples per fuel mix at the standard temperature of 40 • C, see Table 4). During the tests at the Power Take-Off-dynamometer (PTO-dyno), the following variables were measured: • the rotational speed, the torque at the rear PTO of the tractor and the hourly fuel consumption; • the NO x concentration (ppm), the CO concentration (ppm) and the opacity (m −1 ) of the exhaust gases.
Each test at the dyno started only after engine warm-up and included a classic test as prescribed by the Organization for Economic Co-operation and Development (OECD-like bench test) with the fuel pump rack fully opened [50], i.e. a test starting from the maximum engine speed with a rising brake-force controlled by the experimenters. Although the wide-open throttle prescription does not correspond to the most common operation mode of an engine, this test condition is a usual prescription also in international norms (ECE R24 [45], ECE R120 [51], SAE J1349 [52]), aimed at obtaining experimental data that are easily reproducible and fully comparable. Thanks to an instantaneous reading of the PTO speed, it was possible to modulate the braking force applied by the dyno so as to stabilize the tractor at six speeds of rotation (300-400-500-600-700-800 rpm at the PTO, corresponding to 822-1096-1370-1644-1918-2192 rpm at the engine), controlled by the dyno. Each operating condition was kept constant for at least 2 minutes to allow the gas analyser instruments to collect a sufficient number of measurements in steady-state conditions. A total of 30 experimental cases were studied, corresponding to six different speeds per each of the five above-listed fuel blends. The trials were all performed on the same day to have minimal influencing effects of the external environmental conditions, i.e.: an atmospheric pressure of 99 kPa, an ambient temperature of 301 K, a relative humidity of 42% (hence the air-intake correction factor according to [53] is 0.9788). All the measurements were performed thanks to the mobile testing equipment for agricultural tractors that is part of the "Laboratory for Agroforestry Innovation" of the Free University of Bozen-Bolzano. In particular: a "SIGMA 50 Mobile" electric-generator PTO-dyno by N. J. Froment & Co. Ltd. (Easton-on-the-Hill, East Northamptonshire, UK) [54] was used to measure the engine torque and power at the above-mentioned rotational speeds ( Figure 1);

•
the NO x and CO concentration in the exhaust gases [55] at full-load and at six different speeds was measured with a "Vario Plus Industrial" gas analyser by MRU Instruments (Houston, TX, USA) [56]; the gas samples (360 in total; 12 per each of the 30 fuel and speed combinations of values) were taken only after the engine speed was stable at each set value; • the smoke opacity of the exhaust gases (related to the soot measured by the filter smoke number [43,44]) was measured through a "OpTrans 1600" optical Diesel smoke-meter by MRU Instruments [45,56], conforming to NF R 10-025, PTB EO 18-09, ISO 11614 draft and ECE R24, Annex 8/9 specifications; also in this case, the gas samples (90 in total; three per each of the 30 fuel-and-speed combinations of values) were taken only after the engine speed was stable at each set value; the output of this instrument, based on the photometric absorbance principle, is a measurement of the smoke density expressed as absorption coefficient from 0.00 to 9.99 m −1 .  [54] was used to measure the engine torque and power at the above-mentioned rotational speeds ( Figure 1); • the NOx and CO concentration in the exhaust gases [55] at full-load and at six different speeds was measured with a "Vario Plus Industrial" gas analyser by MRU Instruments (Houston, TX, USA) [56]; the gas samples (360 in total; 12 per each of the 30 fuel and speed combinations of values) were taken only after the engine speed was stable at each set value; • the smoke opacity of the exhaust gases (related to the soot measured by the filter smoke number [43,44]) was measured through a "OpTrans 1600" optical Diesel smoke-meter by MRU Instruments [45,56], conforming to NF R 10-025, PTB EO 18-09, ISO 11614 draft and ECE R24, Annex 8/9 specifications; also in this case, the gas samples (90 in total; three per each of the 30 fuel-and-speed combinations of values) were taken only after the engine speed was stable at each set value; the output of this instrument, based on the photometric absorbance principle, is a measurement of the smoke density expressed as absorption coefficient from 0.00 to 9.99 m −1 .
(a)  It is worth mentioning that the emissions of NO x and CO are supplied by the measurement equipment in terms of ppm on a volume basis [12,24,46,[57][58][59][60]. The data elaboration and the following RSM have been carried out using the data with this measurement unit rather than g (kW h) −1 , which is a quantity that can supply an indication of the emissions in compliance also with the indications of the international norms on pollutants (EURO, TIER). If not disposing of a mass flowmeter for the exhaust gases, it is however possible to have an indirect indication of the emissions in g (kW h) −1 . The conversion procedure between ppm and g (kW h) −1 is rather straightforward once the operating conditions and the power output of the engine are known, as reported in the results section of this work. Specifically, a paper by Pilusa et al. [61] supplies the correlations required for the conversion and indicates some practical conversion factors for heavy-duty diesel engines such as the engines equipping agricultural machines.

The Response Surface Methodology
One of the key points of this experimental procedure is to avoid any specific modification to the power unit when the fuelling is changed, hence no re-calibration of the mechanical and of the exhaust gas cleaning after-treatment systems have been performed in the tested tractor. Since the engine responses are not directly proportional to the fuel feed change, it is not possible to define direct analytical correlations based on the physics of the system. Therefore, having simple mathematical relationships at the users' disposal to assess the machine outputs in case of use of alternative fuel blends can be very interesting from a practical point of view. In order to assess and overcome these limitations, a very practical but effective methodology based on a statistical analysis was adopted. Specifically, for each of the investigated variables, an interpolating function was derived to estimate the value of the independent variables (here: the engine torque, the hourly consumptions, the NO x and CO concentration and the opacity of the exhaust gases) by using the values of the dependent variables. From a formal point of view, the final formulation, and hence the mathematical relationships among the quantities, is completely independent from the physics of the represented process, being mainly oriented to obtain numerically-acceptable results within the validity domain and therefore a reliable prediction ability.
All the collected data concerning the mechanical and environmental performances of the engine have been elaborated statistically. After an Analysis of Variance (ANOVA) test, aimed at highlighting the parameters of influence on the above-indicated variable (response), the subsequent application of the Response Surface Methodology (RSM) allowed the calculation of the regression function coefficients that describe the effects of the statistically-significant independent variables on the response [62] ( Figure 2). It is worth mentioning that the emissions of NOx and CO are supplied by the measurement equipment in terms of ppm on a volume basis [12,24,46,[57][58][59][60]. The data elaboration and the following RSM have been carried out using the data with this measurement unit rather than g (kW h) −1 , which is a quantity that can supply an indication of the emissions in compliance also with the indications of the international norms on pollutants (EURO, TIER). If not disposing of a mass flowmeter for the exhaust gases, it is however possible to have an indirect indication of the emissions in g (kW h) −1 . The conversion procedure between ppm and g (kW h) −1 is rather straightforward once the operating conditions and the power output of the engine are known, as reported in the results section of this work. Specifically, a paper by Pilusa et al. [61] supplies the correlations required for the conversion and indicates some practical conversion factors for heavy-duty diesel engines such as the engines equipping agricultural machines.

The Response Surface Methodology
One of the key points of this experimental procedure is to avoid any specific modification to the power unit when the fuelling is changed, hence no re-calibration of the mechanical and of the exhaust gas cleaning after-treatment systems have been performed in the tested tractor. Since the engine responses are not directly proportional to the fuel feed change, it is not possible to define direct analytical correlations based on the physics of the system. Therefore, having simple mathematical relationships at the users' disposal to assess the machine outputs in case of use of alternative fuel blends can be very interesting from a practical point of view. In order to assess and overcome these limitations, a very practical but effective methodology based on a statistical analysis was adopted. Specifically, for each of the investigated variables, an interpolating function was derived to estimate the value of the independent variables (here: the engine torque, the hourly consumptions, the NOx and CO concentration and the opacity of the exhaust gases) by using the values of the dependent variables. From a formal point of view, the final formulation, and hence the mathematical relationships among the quantities, is completely independent from the physics of the represented process, being mainly oriented to obtain numerically-acceptable results within the validity domain and therefore a reliable prediction ability.
All the collected data concerning the mechanical and environmental performances of the engine have been elaborated statistically. After an Analysis of Variance (ANOVA) test, aimed at highlighting the parameters of influence on the above-indicated variable (response), the subsequent application of the Response Surface Methodology (RSM) allowed the calculation of the regression function coefficients that describe the effects of the statistically-significant independent variables on the response [62] (Figure 2   The RSM is a very effective numerical tool that allows the calculation of an explicit polynomial regression-function from a set of input data that is the best approximation, in a limited validity domain, of the real function governing the phenomenon under study [63][64][65][66][67]. As opposed to other mathematical tools such as artificial neural networks [68], the RSM gives as a result an explicit polynomial function that is the first part of the Taylor series of an unknown real function f(x i ). This polynomial function can be used to study and/or optimize a system by making some quantitative predictions about the involved quantities. Design-Expert 7.0.0 software (Stat-Ease, Minneapolis, MN, USA) was used for these calculations.
If y and x i are, respectively, a generic response (i.e., an independent variable) and a generic numerical factor (i.e., a dependent variable), non-coded (i = 1 to m, with m the number of investigated variables; m ≥ 1), a 0 is the interception coefficient, and a i , a ii , a iii , a ij and a ijk (i j k) are respectively the coefficients of the linear, quadratic, cubic, 2 nd -order and 3 rd -order interaction terms, the most generic regression model used in RSM (i.e., a "full-cubic model") is: For the present study, after analysing the data, the software automatically selected for most of the models a (multi-) linear 2 nd -order model with interactions (more precisely, a "full-quadratic model") as starting point to build the mathematical model for each of the investigated responses: The validity domain of each polynomial function f is given by the lower and upper values of each independent variable x i (i = 1 to m) and, therefore, is the following hyperspace: The ANOVA, which is part of this methodology, allows the analyst to identify the most significant factors and polynomial terms, thus operating a partial simplification of the resulting polynomial models based on the p-values (threshold value: 0.05).

System Modelling
The need to introduce here an approach based on the RSM rather than on other tools (e.g., numerical simulations) follows the awareness that the determination of the operative point of a generic machine equipped with an internal combustion engine is not simple since the power conversion occurring inside the combustion chamber is very complex (Figure 3). Indeed, numerical simulations require more technical and physical details to model the involved phenomena. Indicatively, it is necessary to know and to model all the parameters indicated in Figure 3, so this approach requires a complete knowledge of the machine, including several parameters typically known only by the manufacturer or by disassembling the machine. Instead, RSM allows to easily obtain a precise prediction of the value of some quantities that represent the global expression of the operation of the machine (e.g., the torque or the NO x emissions). Indeed, by writing a sort of (polynomial) "transfer function" that characterizes the system "tractor" under consideration (i.e., with its specific engine, auxiliaries, settings, maintenance history, age), it is possible to predict the outputs even in correspondence to input values different from the values used to set up the transfer functions coefficients (i.e., the test conditions) and plot complete output curves/surfaces. While the approach and its application procedure are absolutely general and applicable to every farm tractor, the so-developed polynomial models, characterized by a series of numerical features (i.e., the polynomial order, the selected and remaining terms after a backward elimination leaded by an ANOVA, the values of each coefficient) are specific for the considered tractor. polynomial order, the selected and remaining terms after a backward elimination leaded by an ANOVA, the values of each coefficient) are specific for the considered tractor. Much of the variability of outputs that can be recorded during the experiments on real machines is due to a variety of structural characteristics of the machines under test. By carrying out more tests on the same day and on the same machine as in this study, many parameters of influence (specifically: the engine structural characteristics, the fuel system settings, the environmental conditions) remain constant. Furthermore, the test protocol of OECD standards reduces the variability of the operating settings as it requires the fuel pump rack fully-opened (i.e., the accelerator position at 100%) to be able to detect the properly-said "engine characteristic curves", i.e. the limit curves of the motor in its operating range.   Much of the variability of outputs that can be recorded during the experiments on real machines is due to a variety of structural characteristics of the machines under test. By carrying out more tests on the same day and on the same machine as in this study, many parameters of influence (specifically: the engine structural characteristics, the fuel system settings, the environmental conditions) remain constant. Furthermore, the test protocol of OECD standards reduces the variability of the operating settings as it requires the fuel pump rack fully-opened (i.e., the accelerator position at 100%) to be able to detect the properly-said "engine characteristic curves", i.e. the limit curves of the motor in its operating range.
Among all possible factors of influence on the operating point of a tractor, this study aimed at evidencing the effect of the fuel mix composition by varying the operative conditions in terms of applied load/engine speed (numerical factor controllable by the experimenters through the PTO-dyno). Specifically, the here-reported analyses considered the following characteristics to characterize the fuel blends: 1.
The percentage parts of the three base fuel-components (numerical factors controllable by the experimenters; biodiesel ranges from 0% to 25%, bioethanol can assume only two possible values: 0% or 3% to avoid damages to the fuel system and the engine); please notice that the percentages of the three components composing the fuel mixes are not fully independent from each other as their sum must be always 100; therefore, one of these was excluded by the analyses as it is not properly an "independent variable"; the proposed numerical models have therefore only two of three fuel percentages; 2.
It is worth mentioning that the reference temperature for viscosity, although quite far from the injection temperature of the fuel in the combustion chamber, allows the comparison among different mixes and the formulation of hypothesis on its influence in the spray formation, hence on the performances. Indeed, viscosity at the reference temperature of 40 • C is used in the characterization of all the fuel samplings, thus granting a metrological uniformity. Then again, the fuel viscosity in the injection conditions is not known; it is not the object of any measurement in the literature and it is hardly measurable. Even if available, it would not be significant in the construction of an empirical predictive model. Vice versa, the temperature at 40 • C can be measured with commonly-available instruments. A total of 30 experimental cases were investigated: six engine speeds with five different fuels/viscosity values.

Numeric Models of the System Making Use of the Blend Composition
During the PTO tests, the system tractor, fuelled with one of the presented fuels, stabilizes on a different engine speed every time a new value of the load (expressed as braking torque) is applied by the dyno; therefore, the following equation holds: Thanks to the subsequent elaborations of the quadruple [n(rpm);T(Nm);B(%);E(%)] made with RSM, it is possible to obtain a first model of the torque as a function of the engine speed and the blend composition (expressed through its percentages of biodiesel and ethanol), to be used to make reliable forecasts: During the same experiments at the PTO-dyno, the hourly consumption, the NO x concentration, the CO concentration and the opacity of the exhaust gases by varying the value of the load (i.e., braking torque) applied by the dyno have also been detected, formally stating the following correlations: Opacity Considering the expression for the torque (Equation 5) above, it is possible to correlate each of the above-listed variables directly with the engine speed, thus raising the predictive capabilities of the so-created models by using a parameter very easy to detect: At the end of these elaborations, this approach allows to use a simple tool to infer the desired outputs, as represented in the following Figure 4.
At the end of these elaborations, this approach allows to use a simple tool to infer the desired 439 outputs, as represented in the following One of the key aspects of this study is to propose a correlation between the four explicated

Numeric Models of the System Making Use of the Viscosity
The blend viscosity at 40 • C depends on the blend composition and can be generically expressed by a further equation by using, as independent variables, the percentages of biodiesel and bioethanol in the blend, as follows: One of the key aspects of this study is to propose a correlation between the four explicated responses and the viscosity ( Figure 5). Four new formulations have been obtained, as follows:   (19)   The validity domain of every model presented in this study is the hyperspace defined by the Cartesian product of the following intervals:

Validation of Numerical Models
Analogously to what is commonly done for neural networks [68,74], the acquired data were randomly subdivided into two subsets: 80% of the experimental cases (i.e., 24 cases) were used for building a first set of models (one per response), the remaining 20% of the experimental cases (i.e., 6 cases) were used to verify the predictive capabilities of the different models built from the other 24 cases. By doing so, it was possible to validate the proposed approach also for the present study (i.e., on data of an engine fuelled with biofuel blends). The different models were evaluated in terms of: fitting of the initial set of data (through the coefficient of determination, the adjusted coefficient of determination, the average relative fitting error, the mean squared error) and 2.
predictive capabilities on the values of the several responses for the six verification cases (through the average relative prediction error and the mean squared error).
Below there are the six experimental cases randomly chosen for the validation of the models (control cases). Each is identified by a 3-tuple listing the engine speed (rpm), the biodiesel content (%), the bioethanol content (%), or by the equivalent 2-tuple listing the engine speed (rpm), the kinematic viscosity at the reference temperature of 40 • C (mm 2 s −1 ): The relative error of fitting or prediction (for the six control cases) for the i-th case (identified by the 4-tuple [n(rpm);T(Nm);B(%);E(%)] or 3-tuple [n(rpm);T(Nm);v(mm 2 s -1 )] has been calculated as follows: where: • p i is the prediction of the model for the selected quantity in the i-th case; • m i is the experimental measurement of the selected quantity in the i-th case.
After validating the method, the RSM was applied once again to the whole dataset (30 cases) to give the reader an extended, more complete, model for each of the investigated responses ( Figure 6). Then these models were also evaluated using the same metrics (e.g., the coefficient of determination and the average relative fitting error) to have an indication of the enhanced fitting capabilities given by the addition of six more points to the initial set of 24.

Results and Discussion
The following paragraphs show the final correlations found with the proposed numerical model. The mathematical relationships are obtained using the set of selected terms resulting from the systematic backward-elimination of not-significant terms (one by one). This process was performed as follows: after each single elimination an ANOVA was performed, and the new p-values were calculated. On such basis it was possible to decide whether to proceed with the selection of factors until a set of solely significant factors which would then appear in the final regression model.

Engine Torque as a Function of the Fuel Blend Composition
The following Table 5 shows the numerical regression models, based on 24 and 30 experimental cases respectively, as resulted from the backward selection of the terms forming the most complete model handled by RSM. The equations are presented in terms of actual factors, i.e. the terms are used without any conversion, scaling or non-dimensionalization. The following Figures 7 and 8 give a graphical (2D/3D) representation of these equations, allowing the visualization of the trends and absolute values. Table 5. Main data concerning the two numerical models (respectively based on 24 and 30 cases) and the predictions of the models in six selected cases for the engine torque.

nr. 1 (on 24 cases) nr. 2 (on 30 cases)
Regression equation   Both the models resulting from the backward selection of terms are quadratic even if there is only one quadratic term related to the engine speed. For this reason, by observing Figures 7 and 8, the resulting response surface is a single-curvature surface.
With regards to only the cases used for the model building phase, the predictive capability of both the models is very high: the R 2 is higher than 0.9918 and the adjusted R 2 is higher than 0.9930. When considering the six control cases excluded in the building of the first model, it is possible to confirm the very high predictive capability of this model: for such cases the average error (0.38%) and the mean square error (0.01%) are very low.
The prediction of the engine power at the different rotational speeds and the variation of the fuel mix composition is simple, starting from the same equation of the torque. The predicted value of the torque should be multiplied by the engine speed expressed in rad s −1 . For the above-mentioned reasons, no model for the power has been calculated.

535
With regards to only the cases used for the model building phase, the predictive capability of 536 both the models is very high: the R 2 is higher than 0.9918 and the adjusted R 2 is higher than 0.9930.  The following Table 6 shows the numerical regression models, respectively based on 24

Hourly Fuel Consumption at Full Load as a Function of the Composition of the Fuel Blends
The following Table 6 shows the numerical regression models, respectively based on 24 and 30 experimental cases, as resulted from the backward selection of the terms forming the most complete model handled by RSM. The equations are presented in terms of actual factors, i.e. the terms are used without any conversion, scaling or nondimensionalization. The following Figures 9 and 10 give a graphical (2D/3D) representation of these equations, allowing the visualization of the trends and absolute values. Table 6. Main data concerning the two numerical models (respectively based on 24 and 30 cases) and the predictions of the models in six selected cases for the hourly fuel consumption at full load.

nr. 1 (on 24 cases) nr. 2 (on 30 cases)
Regression equation   Both the models resulting from the backward selection of terms are quadratic, even with a different number of terms. In particular, the 24-case model has a term more, related to the interaction between biodiesel and bioethanol and this is also responsible for the slightly high value of R 2 and adjusted R 2 (the differences are however very low: less than 0.02).
With regards to only the cases used for the model building phase, the predictive capability of both the models is high: the R 2 is higher than 0.9040 and the adjusted R 2 is higher than 0.8790. By looking at the six control cases excluded in the building of the first model, it is possible to confirm a good predictive capability of this model: for those cases, the mean square error is 0.92% and the average error is 2.48%. Even if it is not as low as for other investigated responses, it should be considered that it is however lower than 5%, which is commonly considered the threshold for the acceptability of every measurement and predictive system.
Analogously to the engine power, which can be calculated from the torque prediction, the prediction of the brake specific fuel consumption at the different rotational speeds and varying the fuel mix composition can be obtained from the same equation of the hourly consumption. For the explained reason, no model for the BSFC has been calculated.

570
Both the models resulting from the backward selection of terms are quadratic, even with a 571 different number of terms. In particular, the 24-case model has a term more, related to the interaction 572 between biodiesel and bioethanol and this is also responsible for the slightly high value of R 2 and 573 adjusted R 2 (the differences are however very low: less than 0.02).

574
With regards to only the cases used for the model building phase, the predictive capability of 575 both the models is high: the R 2 is higher than 0.9040 and the adjusted R 2 is higher than 0.8790. By 576 looking at the six control cases excluded in the building of the first model, it is possible to confirm a 577 good predictive capability of this model: for those cases, the mean square error is 0.92% and the 578 average error is 2.48%. Even if it is not as low as for other investigated responses, it should be 579 considered that it is however lower than 5%, which is commonly considered the threshold for the 580 acceptability of every measurement and predictive system.

581
Analogously to the engine power, which can be calculated from the torque prediction, the  The following Table 7 shows the numerical regression models, respectively based on 24 and 30

NO x Concentration in the Exhaust Gases as a Function of the Composition of the Fuel Blends
The following Table 7 shows the numerical regression models, respectively based on 24 and 30 experimental cases, as resulted from the backward selection of the terms forming the most complete model handled by RSM. The equations are presented in terms of actual factors, i.e. the terms are used without any conversion, scaling or nondimensionalization. The following Figures 11 and 12 give a graphical (2D/3D) representation of these equations, allowing to visualize the trends and absolute values. Table 7. Main data concerning the two numerical models (respectively based on 24 and 30 cases) and the predictions of the models in six selected cases for the NO x concentration in the exhaust gases.

nr. 1 (on 24 cases) nr. 2 (on 30 cases)
Regression equation    Within the validity domains, the two resulting numerical models can be used to predict the engine NO x concentration at the different rotational speeds and the variation of the fuel mix composition by simply substituting the values of the input quantities. This is what has been done with the first model and the six control cases.
Both the models resulting from the backward selection of terms have a quadratic form with two of the factors (biodiesel content, engine speed) raised to the second power. The resulting response surface presents double curvature in a shape resembling a roofing-tile. The values of R 2 and adjusted R 2 are very high for both the models (24 and 30 cases) and approximately at the same values (about 0.982).
Firstly, this means that the predictive capability of both the models is high with respect to only the cases used for the model-building phase. Secondly, concerning the very small difference of R 2 values, this means that the six control cases are well aligned with the calculated response surface and that the trend has been completely respected when including also those cases in the second model.
By looking at the same six control cases, it is possible to confirm the very high predictive capability of the first model: indeed, for those cases the mean square error is 0.17% and the average error is 0.36%, both essentially negligible. Within the validity domains, the two resulting numerical models can be used to predict the engine NO x concentration at the different rotational speeds and the variation of the fuel mix composition by simply substituting the values of the input quantities. This is what has been done with the first model and the six control cases.
Both the models resulting from the backward selection of terms have a quadratic form with two of the factors (biodiesel content, engine speed) raised to the second power. The resulting response surface presents double curvature in a shape resembling a roofing-tile. The values of R 2 and adjusted R 2 are very high for both the models (24 and 30 cases) and approximately at the same values (about 0.982).
Firstly, this means that the predictive capability of both the models is high with respect to only the cases used for the model-building phase. Secondly, concerning the very small difference of R 2 values, this means that the six control cases are well aligned with the calculated response surface and that the trend has been completely respected when including also those cases in the second model.
By looking at the same six control cases, it is possible to confirm the very high predictive capability of the first model: indeed, for those cases the mean square error is 0.17% and the average error is 0.36%, both essentially negligible.

CO of the Exhaust Gases as a Function of the Composition of the Fuel Blends
The following Table 8 shows the numerical regression models, respectively based on 24 and 30 experimental cases, as resulted from the backward selection of the terms forming the most complete model handled by RSM. The equations are presented in terms of actual factors, i.e. the terms are used without any conversion, scaling or non-dimensionalization. In this case the multilinear model is given as a function of the natural logarithm of the CO concentration, to better handle the high ratio between the minimum and the maximum experimental values. The following Figures 13 and 14 give a graphical (2D/3D) representation of these equations, allowing the visualization of the trends and absolute values. Within the validity domains, the two resulting numerical models can be used to predict the engine CO concentration at the different rotational speeds and the variation of the fuel mix composition by simply substituting the values of the input quantities. This is what has been done with the first model and the six control cases. Within the validity domains, the two resulting numerical models can be used to predict the engine CO concentration at the different rotational speeds and the variation of the fuel mix composition by simply substituting the values of the input quantities. This is what has been done with the first model and the six control cases.   Due to the very high ratio between the minimum and the maximum experimental values (respectively: 90 and 10,597 ppm), the method and the software suggest operating a transformation of the response data before performing the ANOVA and to build the model in the same form (polynomial). Consequently, the models give an estimation of the natural logarithm of the CO concentration in correspondence at an n-tuple of values for the involved factors. The subsequent calculation of the CO concentration is simple, involving a rising to the power of the Napier's constant. Apart from this initial transformation, both the models resulting from the backward selection of terms have a quadratic form with two of the factors (biodiesel content and engine speed) raised to the second power. The effect of the biodiesel content on the global shape of the surface is however inferior with respect to the engine speed, so the curvature due to the biodiesel content is very small.
The values of R 2 and adjusted R 2 are very high for both the models (24 and 30 cases) and substantially at the same levels (about 0.987). This means, firstly, that the predictive capability of both models is high with regards to only the cases used for the model building phase. Secondly, concerning the very small difference of R 2 values, this means that the six control cases are on average well aligned with the calculated response surface, and the general trend has been completely followed when including also those cases in the second model. By looking at the same six control cases, it is possible to Due to the very high ratio between the minimum and the maximum experimental values (respectively: 90 and 10,597 ppm), the method and the software suggest operating a transformation of the response data before performing the ANOVA and to build the model in the same form (polynomial). Consequently, the models give an estimation of the natural logarithm of the CO concentration in correspondence at an n-tuple of values for the involved factors. The subsequent calculation of the CO concentration is simple, involving a rising to the power of the Napier's constant. Apart from this initial transformation, both the models resulting from the backward selection of terms have a quadratic form with two of the factors (biodiesel content and engine speed) raised to the second power. The effect of the biodiesel content on the global shape of the surface is however inferior with respect to the engine speed, so the curvature due to the biodiesel content is very small.
The values of R 2 and adjusted R 2 are very high for both the models (24 and 30 cases) and substantially at the same levels (about 0.987). This means, firstly, that the predictive capability of both models is high with regards to only the cases used for the model building phase. Secondly, concerning the very small difference of R 2 values, this means that the six control cases are on average well aligned with the calculated response surface, and the general trend has been completely followed when including also those cases in the second model. By looking at the same six control cases, it is possible to notice an acceptable predictive capability of the first model: indeed for those cases, the average error is −0.91%, essentially negligible, but the mean squared error has a bit higher value (4.47%), which is still acceptable but close to the threshold of 5% commonly considered for predictive models. This means that the deviations of the predicted values from the real values are equally distributed above and below the zero, and nevertheless quite significant in absolute terms. It is possible to therefore conclude that the model gives a good approximation of the trend of the plotted quantity (the CO concentration, in this case), but should be used with caution in formulating specific predictions. Nevertheless, this does not constitute a problem in the present case since it means that the model is fully reliable to draw conclusions concerning the influence of the factors on this response.

Opacity of the Exhaust Gases as a Function of the Composition of the Fuel Blends
The following Table 9 shows the numerical regression models, respectively based on 24 and 30 experimental cases, as resulted from the backward selection of the terms forming the most complete model handled by RSM. The equations are presented in terms of actual factors i.e., the terms are used without any conversion, scaling or non-dimensionalization. In this case the multilinear model is given as a function of the natural logarithm of the opacity, to better handle the high ratio between the minimum and the maximum experimental values. Table 9. Main data concerning the two numerical models (respectively based on 24 and 30 cases) and the predictions of the models in six selected cases for the opacity of the exhaust gases. Within the validity domains, the two resulting numerical models can be used to predict the exhaust gas opacity at the different rotational speeds and the variation of the fuel mix composition by simply substituting the values of the input quantities. This is what has been done with the first model and the six control cases. The following Figures 15-17 give a graphical (2D/3D) representation of these equations, allowing the visualization of the trends and absolute values.    Due to the very high ratio between the minimum and the maximum experimental values (respectively: 0.002 and 3.081 m −1 ), the method and software suggest operating a transformation of the response data before performing the ANOVA and to build the model in the same form (polynomial). Consequently, the models give an estimation of the natural logarithm of the opacity of exhaust gases in correspondence at an n-tuple of values for the involved factors. The subsequent calculation of the opacity is simple, involving a rising to the power of the Napier's constant. Apart from this initial transformation, both the models resulting from the backward selection of terms have a quadratic form with two of the factors (biodiesel content, engine speed) raised to the second power. The resulting shape of the response surface is a basin, hence with a sensible contribution of the biodiesel content and engine speed to the curvature of the surface.
The values of R 2 and adjusted R 2 are quite high for both the models (24 and 30 cases) and substantially at the same levels (beyond 0.899). This means, firstly, that the predictive capability of both the models is high with regards to only the cases used for the model building phase. Secondly, concerning the very small difference of R 2 values, this means that the six control cases are, on average, well aligned with the calculated response surface, and the general trend has been completely respected when including also those cases in the second model. The slight differences of R 2 and adjusted R 2 are essentially due to the simplification of the model form (exclusion of the term with the interaction between biodiesel and engine speed) operated by the software.
By looking at the same six control cases, it is possible to notice a not-acceptable predictive capability of the first model: indeed, for those cases, the average error is 5.62% and the mean squared error has a Due to the very high ratio between the minimum and the maximum experimental values (respectively: 0.002 and 3.081 m −1 ), the method and software suggest operating a transformation of the response data before performing the ANOVA and to build the model in the same form (polynomial). Consequently, the models give an estimation of the natural logarithm of the opacity of exhaust gases in correspondence at an n-tuple of values for the involved factors. The subsequent calculation of the opacity is simple, involving a rising to the power of the Napier's constant. Apart from this initial transformation, both the models resulting from the backward selection of terms have a quadratic form with two of the factors (biodiesel content, engine speed) raised to the second power. The resulting shape of the response surface is a basin, hence with a sensible contribution of the biodiesel content and engine speed to the curvature of the surface.
The values of R 2 and adjusted R 2 are quite high for both the models (24 and 30 cases) and substantially at the same levels (beyond 0.899). This means, firstly, that the predictive capability of both the models is high with regards to only the cases used for the model building phase. Secondly, concerning the very small difference of R 2 values, this means that the six control cases are, on average, well aligned with the calculated response surface, and the general trend has been completely respected when including also those cases in the second model. The slight differences of R 2 and adjusted R 2 are essentially due to the simplification of the model form (exclusion of the term with the interaction between biodiesel and engine speed) operated by the software.
By looking at the same six control cases, it is possible to notice a not-acceptable predictive capability of the first model: indeed, for those cases, the average error is 5.62% and the mean squared error has a high value (40.07%), both beyond the threshold of 5% commonly considered for predictive models. This means that the deviations of the predicted values from the real values are significant and unevenly distributed with respect to the zero (the positive sign means that model generally overestimates the experimental situation). Although the model can give a good approximation of the trend of the plotted quantity (in this case, the opacity) and it can be used to draw conclusions concerning the influence of the factors on this response, it should be used with extreme caution in formulating specific predictions.

Mechanical and Environmental Outputs as a Function of the Viscosity of the Fuel Blends
In the following paragraphs the presented methodology has been applied to outline a correlation between the viscosity of the fuel blend and the motor outputs. This attempt was made despite it being known that the viscosity is certainly not the only parameter influencing the outputs of an engine (other such parameters surely being the LHV and the density), as confirmed by some of the results presented below. However, this attempt was carried out since the supply of bioethanol or biodiesel may not be constant during the year if considering the production process of these fuels from the perspective of a fully self-sufficient farm. Therefore, the most practical way to characterize a blend contained in a tank without using complex and expensive procedures is a quick viscosity test.

Engine Torque as a Function of the Viscosity of the Fuel Blends
The following Table 10 shows the numerical regression models, respectively based on 24 and 30 experimental cases, as resulted from the backward selection of the terms forming the most complete model handled by RSM. The equations are presented in terms of actual factors, i.e. the terms are used without any conversion or scaling. The following Figures 18 and 19 give a graphical (2D/3D) representation of these equations, allowing the visualization of the trends and absolute values.
Within the validity domains, the two resulting numerical models can be used to predict the engine torque at the different rotational speeds and the variation of the fuel mix kinematic viscosity by simply substituting the values of the input quantities. This is what has been done with the first model and the six control cases. Table 10. Main data concerning the two numerical models (respectively based on 24 and 30 cases) and the predictions of the models in six selected cases for the engine torque.

777
The values of R 2 and adjusted R 2 are high for both the models (24 and 30 cases) and 778 substantially at the same levels (beyond 0.9806). This means, firstly, that the predictive capability of

777
The values of R 2 and adjusted R 2 are high for both the models (24 and 30 cases) and  Table. are reported as labels directly on the contours). The same pictures show also the positions of the experimental measurements (black dots positioned within the investigated hyperspace; notice that many points are positioned under the represented surface and hence are not visible).
Both the models resulting from the backward selection of terms have a quadratic form, even if the 30-case model is more complete, having two of the factors (engine speed, kinematic viscosity) raised to the second power. The resulting shape of the response surface is somehow similar, due to the low influence of the kinematic viscosity on the curvature.
The values of R 2 and adjusted R 2 are high for both the models (24 and 30 cases) and substantially at the same levels (beyond 0.9806). This means, firstly, that the predictive capability of both the models is high with regards to only the cases used for the model building phase. Secondly, concerning the very small difference of R 2 values, this means that the six control cases are on average well aligned with the calculated response surface, and the general trend has been completely followed when including also those cases in the second model. The slight differences of R 2 and adjusted R 2 of the two models are essentially due to the completion of the model form (inclusion of the term with the kinematic viscosity raised to the 2 nd power) operated by the software.
By looking at the same six control cases, it is possible to notice a high predictive ability of the first model: indeed, for those cases the average error is 0.40% and the mean squared error has an even lower value (0.02%), both far below the common threshold of 5% usually considered for predictive models. This means that the deviations of the predicted values from the real values are equally distributed above and below the zero, and the predicted values are quite close to real values. The model can give a very good approximation of the trend of the plotted quantity (the torque, in this case), so it can be used to draw conclusions concerning the influence of the factors on this response. Moreover, it can also be used to make precise specific predictions.
Here, the prediction of the engine power at the different rotational speeds and the variation of the fuel mix viscosity is obtained starting from the same equation of the torque; the predicted value of the torque should be multiplied by the engine speed expressed in rad s −1 . For the explained reason, no model for the power has been calculated.

Hourly Fuel Consumption at Full Load as a Function of the Viscosity of the Fuel Blends
The following Table 11 shows the numerical regression models, respectively based on 24 and 30 experimental cases, as resulted from the backward selection of the terms forming the most complete model handled by RSM. The equations are presented in terms of actual factors, i.e. the terms are used without any conversion, scaling or non-dimensionalization. The following Figures 20 and 21 give a graphical (2D/3D) representation of these equations, allowing the visualization of the trends and absolute values. Table 11. Main data concerning the two numerical models (respectively based on 24 and 30 cases) and the predictions of the models in six selected cases for the hourly fuel consumption at full load.

819
Within the validity domains, the two resulting numerical models can be used to predict the 820 engine hourly fuel consumption at full load at the different rotational speeds and the variation of the

828
The values of R 2 and adjusted R 2 are quite high for both the models (24 and 30 cases) and 829 substantially at the same levels (beyond 0.857). This means, firstly, that the predictive capability of

819
Within the validity domains, the two resulting numerical models can be used to predict the 820 engine hourly fuel consumption at full load at the different rotational speeds and the variation of the

828
The values of R 2 and adjusted R 2 are quite high for both the models (24 and 30 cases) and 829 substantially at the same levels (beyond 0.857). This means, firstly, that the predictive capability of 830 both the models is quite high with regards to only the cases used for the model building phase. The 831 very small difference of R 2 values between the two models means that the six control cases are, on Within the validity domains, the two resulting numerical models can be used to predict the engine hourly fuel consumption at full load at the different rotational speeds and the variation of the fuel mix composition by simply substituting the values of the input quantities. This is what has been done with the first model and the six control cases.
Both the models resulting from the backward selection of terms have a quadratic form with the same number of terms, having both the two considered factors (engine speed and kinematic viscosity) raised to the second power. The resulting shape of the response surface is substantially equal (a dome), with the kinematic viscosity having a little influence on the curvature in the considered domain.
The values of R 2 and adjusted R 2 are quite high for both the models (24 and 30 cases) and substantially at the same levels (beyond 0.857). This means, firstly, that the predictive capability of both the models is quite high with regards to only the cases used for the model building phase. The very small difference of R 2 values between the two models means that the six control cases are, on average, well aligned with the calculated response surface, and the general trend has been completely followed also including those cases in the second model. By looking at the same six control cases, it is possible to notice a high predictive capability of the first model: indeed, for those cases, the average error is −0.86% and the mean squared error is 1.05%, both far below the common threshold of 5% usually considered for predictive models. This means that the deviations of the predicted values from the real values are equally distributed above and below the zero and the predicted values are quite close to real values. The model can give a good approximation of the trend of the plotted quantity (the hourly consumption, in this case), so it can be used to draw conclusions concerning the influence of the factors on this response. Moreover, it can also be used to make quite precise specific predictions.
Also in this case, as presented above the prediction of the brake specific fuel consumption at the different rotational speeds and the variation of the fuel mix viscosity is simple starting from the same equation of the hourly consumption. For the explained reason, no model for the BSFC has been calculated.

NO x Concentration in the Exhaust Gases as a Function of the Viscosity of the Fuel Blends
The following Table 12 shows the numerical regression models, respectively based on 24 and 30 experimental cases, as resulted from the backward selection of the terms forming the most complete model handled by RSM. The equations are presented in terms of actual factors, i.e. the terms are used without any conversion or scaling or non-dimensionalization. The following Figures 22 and 23 give a graphical (2D/3D) representation of these equations, allowing the visualization of the trends and absolute values.
Within the validity domains, the two resulting numerical models can be used to predict the NO x concentration at the different rotational speeds and the variation of the fuel mix kinematic viscosity by simply substituting the values of the input quantities. This is what has been done with the first model and the six control cases.
Both the models resulting from the backward selection of terms have a cubic form with only the kinematic viscosity) raised to the third power. The resulting shape of the response surface is substantially equal (a sort of wave), with the kinematic viscosity having a significant influence on the curvature in the considered domain. Table 12. Main data concerning the two numerical models (respectively based on 24 and 30 cases) and the predictions of the models in six selected cases for the NOx.

876
The values of R 2 and adjusted R 2 are very high for both the models (24 and 30 cases) and 877 substantially at the same levels (beyond 0.939). This means, firstly, that the predictive capability of response. Moreover, it can also be used to make quite precise punctual predictions.
By looking at the same six control cases, it is possible to notice a generally low predictive 883 capability of the first model: indeed, for those cases, the average error is 7.56% and the mean squared 884 error is 2.67%, the latter far below the common threshold of 5% usually considered for predictive  response. Moreover, it can also be used to make quite precise punctual predictions.
By looking at the same six control cases, it is possible to notice a generally low predictive 883 capability of the first model: indeed, for those cases, the average error is 7.56% and the mean squared 884 error is 2.67%, the latter far below the common threshold of 5% usually considered for predictive  The values of R 2 and adjusted R 2 are very high for both the models (24 and 30 cases) and substantially at the same levels (beyond 0.939). This means, firstly, that the predictive capability of both the models is very high with regards to only the cases used for the model building phase. The model can give a good approximation of the trend of the plotted quantity (the hourly consumption, in this case), so it can be used to draw conclusions concerning the influence of the factors on this response. Moreover, it can also be used to make quite precise punctual predictions.
By looking at the same six control cases, it is possible to notice a generally low predictive capability of the first model: indeed, for those cases, the average error is 7.56% and the mean squared error is 2.67%, the latter far below the common threshold of 5% usually considered for predictive models. Looking at the deviations of the predicted values with respect to the real values, it is possible to notice that the high average error is substantially due to the previsions on the two cases using B25E3, which have the same kinematic viscosity of B25 (3.264 vs. 3.265 mm 2 s −1 ) but surely not the same NO x emissions: at the same engine speed the emissions with B25 are 100 ppm higher than B25E3. Excluding those two cases, the average error drops down at −1.10%.

CO of the Exhaust Gases as a Function of the Viscosity of the Fuel Blends
The following Table 13 shows the numerical regression models, respectively based on 24 and 30 experimental cases, as resulted from the backward selection of the terms forming the most complete model handled by RSM. The equations are presented in terms of actual factors, i.e. the terms are used without any conversion, scaling or non-dimensionalization. In this case the multilinear model is given as a function of the natural logarithm of the CO concentration, to better handle the high ratio between the minimum and the maximum experimental values. The following Figures 24 and 25 give a graphical (2D/3D) representation of these equations, allowing the visualization of the trends and absolute values.
Within the validity domains, the two resulting numerical models can be used to predict the CO concentration at the different rotational speeds and the variation of the fuel mix kinematic viscosity by simply substituting the values of the input quantities. This is what has been done with the first model and the six control cases. Table 13. Main data concerning the two numerical models (respectively based on 24 and 30 cases) and the predictions of the models in six selected cases for the CO concentration in the exhaust gases.

891
The following Table 13 shows the numerical regression models, respectively based on 24 and 30 892 experimental cases, as resulted from the backward selection of the terms forming the most complete    response. Moreover, it can also be used to make quite precise specific predictions.

929
By looking at the same six control cases, a generally low predictive ability of the first model can 930 be observed: indeed, for those cases the average error is −5.34% and the mean squared error is 2.88%.

931
Notwithstanding the good fitting operated on the 24 cases used to set up the model coefficients, the 932 use of the fuel viscosity to predict the value of the CO concentration at several engine speeds is not 933 recommended. The following Table 14    Although the two models resulting from the backward selection of terms have different forms (the 24-case model is cubic, the 30-case model is quadratic), the resulting shape of the response surface is quite similar. The contribution of the engine speed on the response surface shape is prevailing over the viscosity, and the values of the natural logarithm of the CO concentration sharply decrease with the rising of the speed, reaching a sort of plateau above 1370 rpm.
The values of R 2 and adjusted R 2 are very high for both the models (24 and 30 cases) and substantially at the same values (beyond 0.985). This means, firstly, that the predictive capability of both the models is very high with regards to only the cases used for the model building phase. The model can give a good approximation of the trend of the plotted quantity (the CO concentration, in this case), so it can be used to draw conclusions concerning the influence of the factors on this response. Moreover, it can also be used to make quite precise specific predictions.
By looking at the same six control cases, a generally low predictive ability of the first model can be observed: indeed, for those cases the average error is −5.34% and the mean squared error is 2.88%. Notwithstanding the good fitting operated on the 24 cases used to set up the model coefficients, the use of the fuel viscosity to predict the value of the CO concentration at several engine speeds is not recommended.

Opacity of the Exhaust Gases as a Function of the Viscosity of the Fuel Blends
The following Table 14 shows the numerical regression models, respectively based on 24 and 30 experimental cases, as resulted from the backward selection of the terms forming the most complete model handled by RSM. The equations are presented in terms of actual factors, i.e. the terms are used without any conversion, scaling or non-dimensionalization. In this case the multilinear model is given as a function of the natural logarithm of the opacity, to better handle the high ratio between the minimum and the maximum experimental values. The following Figures 26 and 27 give a graphical (2D/3D) representation of these equations, allowing the visualization of the trends and absolute values.     Within the validity domains, the two resulting numerical models can be used to predict the exhaust gases opacity at the different rotational speeds and the variation of the fuel mix kinematic viscosity by simply substituting the values of the input quantities. This is what has been done with the first model and the six control cases.
The two models resulting from the backward selection of terms have the same forms: both have a quadratic form with the viscosity term reaching the second power, even if with a different number of terms. The resulting shape of the response surface is quite similar and complex for the logarithm of the CO concentration by varying the viscosity of the blend, and simpler when considering the CO concentration directly: there is a wide plateau involving all cases having the engine speed greater than about 1000 rpm. The contribution of the engine speed on the response surface shape is well outlined in both the models (i.e., by increasing the engine speed the value of the logarithm decreases). The viscosity seems to have little influence on the figure of the opacity; only in the graph of the logarithm of the opacity this (however little) influence is evidenced by a slight waviness of the surface.
The values of R 2 and adjusted R 2 are quite high for both the models (0.8781 and 0.8616 respectively for the 24-and 30-case model). This means, firstly, that the predictive capability of both the models is quite high with regards to only the cases used for the model-building phase. Both the models can give a quite good approximation of the trend of the plotted quantity (the opacity, in this case), so they can be used to draw conclusions concerning the influence of the factors on this response. However, it is not recommended to use them to make precise specific predictions due to the high prediction errors (using the first model to make predictions, the average error is 18.95% and the mean squared error is  20.84%, both far above the threshold of 5%) shown for the six control cases. Therefore, it is possible to conclude that the use of the fuel viscosity to predict the opacity value at the several engine speeds is not recommended. Clearly, there are different fuel blends having the same viscosity but producing very different effects of the exhaust gases in terms of opacity, in this case.

Comments
As can be seen in the following Table 15, the use of the RSM allows the highly precise expression of all the investigated variables concerning the main mechanical and environmental performances of an engine as a function of the fuel blend composition (in terms of percentages of biodiesel and bioethanol in the blends) and the engine speed. Observing the reported values of R 2 and adjusted R 2 , it is possible to conclude that the chosen factors are very good descriptors of the complex phenomena occurring in the whole energy generation system such as the engine of an agricultural machine. In particular, the models for the torque, the hourly consumption, the NO x concentration, the CO concentration and the opacity all present very high values of R 2 (above 0.90).
Moreover, within their validity domain, the proposed polynomial numerical models for all these factors are very accurate approximations of the real response functions, and their trends (expressed as first and second derivative) can be used instead of the trends of the real (unknown) response functions to investigate the machine response to a change of fuelling and formulate interesting considerations.
Finally, excluding the models for hourly consumption and opacity, the same models are also applicable in formulating numerical predictions of the investigated factors. It is therefore possible to conclude that RSM can express its full potential also for ICEs when it is applied to develop correlations involving the fuel blend composition and the engine speed. Hence, the use of RSM is highly recommended in such a scenario as the one presented in the current paper. In the following Table 16 it is possible to see that the use of RSM expresses the principal mechanical performances of an engine (i.e., the torque) as a function of the fuel blend kinematic viscosity and the engine speed with a high precision. This means that the fuel blend viscosity is a good descriptor for the Energies 2019, 12, 2287 38 of 45 mechanical behaviour of the tractor power unit (both R 2 and adjusted R 2 are above 0.98), probably due to the operative flexibility of the injection system, which is able to compensate small differences in fuels LHV by varying the injected fuel quantity. Within their validity domain, the proposed polynomial numerical models for this variable are good approximations of the real response functions, and their trends (which can be expressed as first and second derivative) can be used instead of the trends of the real (unknown) response functions to investigate the machine response to a change of fuelling and to formulate interesting considerations. Also, considering that some blends present a different composition (hence also LHV) but a very similar viscosity (D100 and B15, B25 and B25E; please see Table 3), this entails that the fuel blend viscosity is the principal driver for a correct combustion of the fuel and in determining machine performance more than the different LHVs of the blends. This is displayed also by the high predictive capability of the torque model based on the kinematic viscosity and the engine speed. Nevertheless, the same conclusion cannot be drawn for the variable related to the fuel consumption (hourly fuel consumption) where, instead, the influence of the blend composition (through the values of LHV) is evident and reasonable. Finally, although the models concerning the three factors related to environmental performances of an engine (NO x and CO concentration, the opacity of exhaust gases) have a quite high R 2 and adjusted R 2 (both above 0.87), they cannot be considered reliable in the point predictions of the system response. Indeed, the high R 2 values are due to the very complex models obtained by RSM (i.e., with a high polynomial degree and perhaps with a logarithmic transformation of initial data), but the high non-linearity of the response as well as the poor reliability of viscosity as indirect estimator for these variables make it inconvenient to use these models to predict exclusively the above-mentioned emissions. The very high values of the average error and mean squared error obtained by the first model for the six control cases confirm this, especially for the NO x concentration and the opacity of the exhaust gases.
Therefore, it is possible to conclude that RSM applied to the performance of ICEs can be applied also in combination with the kinematic viscosity, but only for the prediction of the mechanical performances of the ICE. Rather, the use of RSM applied on viscosity data is not recommended for pollution-related factors such as NO x concentration, CO concentration and the opacity of exhaust gases, even though literature demonstrates that a physical correlation exists between these parameters (viscosity and emissions).
In the following Table 17, the positive or negative correlations of each factor (column headings) with the indicated response (row headings) are reported (30-case models only). A positive correlation between a factor and a response indicates that, when the response value increases/decreases, an increase/decrease of the factor value is detected, meaning that the first partial derivative of the response with respect to the investigated factor has a positive value. Vice versa, there is a negative correlation of a response with a factor when the increasing/decreasing of that factor causes a decrease/increase of the response (being all other factors constant); in this case the value of the partial derivative of the response with respect to the investigated factor is negative. In some cases, the trend of the response is not fully determined: for example, the response can first increase and then decrease with the increasing of a factor, meaning that the first partial derivative of the response, made with respect to the investigated factor, is not constant. There is also another scenario to be considered only for the models using the blends composition as factors: the response can increase with the increasing of a factor only for a set of values of the other factor and decrease in other parts of the domain. In this case, the investigated first derivative has a dependence also with the other factor. For both of these two cases a complex correlation is reported in Table 17, thus directing the reader towards the graphical representations of the response surface.
As can be noticed, the dependence on the viscosity is complex in almost all the considered responses. This is to be expected, since a single influence factor has been voluntarily considered over a set of influence parameters, including e.g. the LHV and the density. In practice, when looking at responses related to the mechanical outputs, the kinematic viscosity is a good descriptor for the system response while the same parameter, even if not so influent, is not recommended for the pollution-related responses.

Conclusions
This article presents the application of the RSM to the prediction of the performance and the emissions of a tractor engine fuelled with different biofuels blends. The RSM is a numerical methodology based on finding a local polynomial approximation of the response of a system while varying many factors within a set validity domain. Its use is particularly indicated for very complex systems that can hardly be described analytically or with a simulation, but require a quite wide set of cases with many combinations of the values for the possible influencing factors. In particular, the numerical models for torque, fuel consumption and emissions of CO, NO x and exhaust opacity have been assessed for a New Holland T4020V farm tractor fuelled with five fuel blends of diesel oil, biodiesel and bioethanol in different percentages. The application of the RSM indicated that the fuel composition (expressed as percentages of alternative fuels in the blend) has clear effects on the above-listed engine performances (mechanical and environmental outputs) and the related mathematical models show a high fitting (R 2 > 0.90 for both mechanical and environmental outputs) and predictive capabilities (absolute average errors spanning from 0.38% to 5.62% on six check cases). Therefore, it can be concluded that RSM applied to fuel blend composition is very effective in predicting the performance of an agricultural engine.
The kinematic viscosity, a physical parameter very easily obtained also for a fuel mix of unknown composition, is a particularly powerful indirect estimator only for the mechanical outputs of the motor (torque, hourly consumption) of this tractor. The average errors in the prediction of the outputs are all below 0.86% in this case. Conversely, the study also demonstrated that the motor environmental outputs (CO, NO x and opacity) cannot be reliably predicted by the kinematic viscosity and the engine speed alone, indicating that those responses are influenced also by other not-negligible parameters. Indeed, the RSM showed an acceptable fitting on the initial data (R 2 > 0.86), but poor predictive capabilities (absolute average errors above 5.34% on six check cases).
In all the illustrated cases (models using the percentages of alternative fuels in the blend or using the kinematic viscosity), even if the predictive capability of the model is not high, the RSM proves to be a very effective tool for experimenters and machine manufacturers. Indeed, the high fitting obtained in all the RSM models indicates that the proposed numerical models adequately describe the system response to a variation of the considered factors. Hence, the proposed polynomial models can be used to investigate the effect of a change of fuelling on the machine response in term of increasing/decreasing trends: the models can be used to formulate interesting hypothesis on the characteristics which the fuel blends should exhibit before starting new experimental surveys. Funding: This work was supported by the Open Access Publishing Fund of the Free University of Bozen/Bolzano. Also, this research was developed within the "VISCOMOTOR" ("Effects of biofuels on lubricants performance in internal combustion for agricultural purposes") internal research project, grant number TN200D, and the "BIO-TRACT-EFFICIENCY" ("Experimental investigation on the efficiency of agricultural machines powered with different fuels") internal research project, grant number TN200F.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.