**1. Introduction**

Nowadays internal combustion engines play a crucial role in worldwide mobility. A general concern about climate change imposes the adoption of an environmentally sustainable policy, as well relevant efforts towards a reduction of fuel consumption and compliance with strict emission regulations [1,2]. In the past decades, analytical and numerical models provided fundamental support to the design of IC engines, in addition to the necessary experimental campaign on real engine prototypes. Today reliable 1D CFD tools are commonly used in the community of engineers and researchers to investigate the behavior and performance of new engines [3,4]. In general, state of the art 1D simulation models can calculate the steady-state engine maps (in terms of air mass flow rate, torque and power, max. cylinder pressure, fuel consumption and emissions), being able to represent every operating condition just starting from a few operating points analyzed in

**Citation:** Marinoni, A.; Tamborski, M.; Cerri, T.; Montenegro, G.; D'Errico, G.; Onorati, A.; Piatti, E.; Pisoni, E.E. 0D/1D Thermo-Fluid Dynamic Modeling Tools for the Simulation of Driving Cycles and the Optimization of IC Engine Performances and Emissions. *Appl. Sci.* **2021**, *11*, 8125. https://doi.org/ 10.3390/app11178125

Academic Editors: Georgios Karavalakis and Ricardo Novella Rosa

Received: 3 May 2021 Accepted: 28 August 2021 Published: 1 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the experimental laboratory (map-based approach) [5,6]. This is a very powerful and fast simulation approach, but the hypothesis of steady state working condition represents a major limitation. In fact, the IC engine is an unsteady volumetric machine, and its behavior is deeply influenced by transient conditions due to mechanical, thermal and fluid-dynamics inertias, not considered in a map-based approach. A robust simulation method is required to account for the transient conditions of the IC engine coupled to the vehicle and its components, to obtain more realistic results in terms of fuel consumption and emissions [7].

The present work wants to highlight the application of the new simulation platform based on Gasdyn, for the 1D thermo-fluid dynamic engine simulation, and Velodyn, for the simulation of the vehicle longitudinal dynamics and powertrain composition, to cope with a typical transient condition simulation. The two simulation tools are integrated in a "co-simulation" framework, which allows the two software to interact in real-time within the MATLAB® Simulink® environment. The vehicle simulation tool interacts with the virtual engine to proceed forward in time and follow the prescribed test cycle.

In particular, the adoption of a specific numerical scheme, based on a staggered leapfrog, finite volume method applied to coarse meshes, allows to preserve the mass conservation and accuracy of the 1D fluid dynamic modeling of the IC engine, resulting in an optimal solution for the simulation of driving cycles in time frames comparable to real-time. In this way, the 1D simulation tool allows the investigation and comparison of different solutions of hybrid architectures in terms of performance, fuel consumption and emissions.

A six-cylinder, 210 kW diesel engine has been modeled in detail and then validated in steady-state conditions, on the basis of experimental engine map data. Then the 1D engine model was made compatible for the two-way coupling with the vehicle simulation tool, via an S-function block in the Simulink® framework.

Finally, the Co-Simulation environment has been applied and validated against the traditional map-based approach, which is based on the available experimental engine maps of the IC engine. The test case considered is a hybrid powertrain of a commercial bus. An analysis of fuel consumption in real driving emission (RDE) cycles is discussed, considering both the conventional and hybrid architectures of the bus, mainly focusing on the hybrid components design and hybrid control unit (HCU) logic settings.

#### **2. 0D/1D Thermo-Fluid Dynamic Modeling of IC Engine**

The 1D thermo-fluid dynamic tool used to carry out the simulation of the Diesel engine is based on the solution of the conservation equations of mass, momentum and energy along the IC engine domain, modelled by means of one-dimensional pipes, to represent the complete intake and exhaust duct systems, and zero-dimensional elements, to represent junctions of pipes, abrupt area changes, valves, turbocharger and cylinders. For what concerns the fluid dynamics, the gas flow is assumed to be one-dimensional, unsteady, compressible with friction and heat transfer at the pipe walls. Moreover, the transport equations for the tracking of chemical species along the pipes could have been considered as well, but in this study, this aspect was not considered. This comprehensive approach is extensively described in literature, see for example [5,6,8–10]. In this work, an innovative contribution is represented by the development and application of a fast simulation method (FSM) to reduce the computational time and reduce the CPU/real-time ratio. The FSM numerical solver developed can achieve satisfactory conservation of mass with coarse meshes adopted to discretize the ducts, up to 10 cm, to significantly reduce the computational effort, because of the decreased number of computational nodes and increased time step. Traditional numerical methods, based on finite difference techniques with 2nd order accuracy [11,12], behave as non-conservative when the calculation domain is so coarsely discretized [13]. Hence, an alternative numerical solver which ensures a good conservation with such large meshes has been applied, as briefly described below. The numerical method is a one-dimensional finite volume technique based on cells and connectors, named 1Dcell, explicit in time and 2nd order accurately derived from a quasi3D formulation [14]. The transported quantities are defined in the cell centers, while their fluxes are established on the connectors. The space-time grid is staggered, as shown in Figure 1. The centers of the computational control volumes for the solution of the momentum equation are represented as small squares, while the cell control volume centers for the solution of mass and energy equations are represented by circles. They are located at different positions in both space and time.

**Figure 1.** (**a**) One-dimensional 1Dcell approach, cells and connectors; (**b**) time marching procedure adopted for the solution of the governing equation in the 1Dcell method: (A) represents the cell control volume and (B) represents the port control volume.

With the above-mentioned reference frame, the continuity and energy conservation equations are solved with respect to the cell element, whereas the momentum equation is solved with respect to the connector elements between cells. The grid is also staggered in time as depicted by Figure 1. The closing equation is the perfect gas equation of state. The conservation equations of mass, energy and momentum are reported below:

$$
\rho\_{cell}^{t+1} = \rho\_{cell}^t + \frac{\Delta t}{V} \sum\_{P}^{N\_{Prot}} \left(\rho v A\right)\_P^{t+\frac{1}{2}} \tag{1}
$$

$$\left(\rho \epsilon^{0}\right)\_{cell}^{t+1} = \left(\rho \epsilon^{0}\right)\_{cell}^{t} + \frac{\Delta t}{V} \sum\_{P}^{N\_{\text{Prot}}} \left(\rho \upsilon h^{0} A\right)\_{P}^{t+\frac{1}{2}} + S\_{t} \tag{2}$$

$$\left(\rho vA\right)\_{cell}^{t+\frac{1}{2}} = \left(\rho \stackrel{\rightarrow}{v} A\right)\_{cell}^{t-\frac{1}{2}} + \frac{\Delta t}{\Delta L} \sum\_{P}^{Nports} \left[\left(p + \rho v^{2}\right)\_{P} A\_{P}\right] + S\_{M\nu} \tag{3}$$

where *ρ* is the gas density, *p* is the gas pressure, *A* is the pipe section and *V* is the cell volume; Δ*t* and Δ*x* the time step size and the distance between cell centers, respectively; *h*0 is the total enthalpy of the cell; *e*0 is the total internal energy of the cell; *Se* and *Sm* terms are the sources of the corresponding equations which consider heat transfer and friction with the duct wall.

This 1Dcell method has been applied for the solution with coarse meshes, to achieve a robust fast simulation method (FSM). Regarding the accurate solution in the case of ordinary mesh sizes (1 cm), the Corberán–Gascon TVD, finite difference, symmetric, explicit method was applied [15]. The TVD feature of this explicit method allows the reduction of the intrinsic instabilities or spurious oscillations which may appear using the explicit formulation. For both numerical methods, the time step is computed according to the Courant–Friedrichs–Lewy condition described by [16].

Intra-pipe boundary conditions are treated resorting to a characteristic-based approach [11,16] to handle junctions, valves, sudden area changes, intercoolers and turbochargers. In particular, for pressure loss junctions experimentally derived loss coefficients are used. These data are available from the literature and are obtained mainly in steady-state flow conditions. Similarly, for poppet valves, measured flow coefficients have been used. Their applicability can be argued but, in any case, they allow a good representation of the physical phenomena with a reasonable approximation. The coupling with 3D models can be realized, improving the accuracy of the simulation, but this brings, as expected, an increase of computational burden that cannot be accepted when the calculation of a driving cycle is addressed [9]. A major concern in this type of simulation is also the approximation introduced by the modeling of the turbocharger. Usually, the compressor and the turbine are treated relying on performance maps, obtained experimentally in state-state flow conditions. The maps are provided by the supplier of the turbocharger or, alternatively, they can be obtained by measuring the turbocharger performances on the test bench [17]. In both cases, the result is a discrete map where information is available for specific revolution speeds. To extend the validity of the maps in regions outside of the experimental ones and between the single iso-velocity lines, interpolation techniques are used. The reliability of the procedure depends mainly on the interpolation technique adopted, but also on the pre-processing of the raw data of the map, which needs to be properly smoothed [18]. Alternatively, a few 1D models can solve the 1D conservation equations for the stator and the rotor, considering the effect of the rotation of the turbine and compressor wheel [19,20]. These are known as map-less approaches, which have a reduced level of approximation, but they need in any case a bit of calibration and detailed data about the turbo geometry. This approach, going in the direction of a better approximation, brings the increase of the computational burden, which, once again, cannot be accepted when the goal is to reach a real-time simulation of a driving cycle.

## **3. Engine Description**

The engine considered is a 6.7 L, turbocharged Diesel engine for heavy-duty applications such as buses and trucks, developed by FPT Industrial. This engine belongs to a commercial family, with a maximum power ranging from 184 kW to 235 kW. In this work, the medium-size 210 kW engine has been investigated.

Specific information about the engines studied is summarized in Tables 1 and 2.


**Table 1.** Engine geometrical data.

**Table 2.** Engine performance data.


The engine is equipped with a twin-entry turbocharger characterized by a fixed geometry and a wastegate valve.

A wide set of experimental data was available to validate the engine thermo-fluid dynamic model. All these data refer to different steady-state operating points, defined by engine speed and brake mean effective pressure (BMEP). The operating map of the 210 kW engine, investigated in this study, is represented by the red points in Figure 2.

**Figure 2.** 210 kW (red points) and 235 kW (blue circles) engine operating maps (BMEP as function of engine speed).

## **4. Calculation Methodology**

As anticipated, the final aim of this work is to carry out the transient simulation of the IC engine integrated with a complete vehicle model, considering the real driving conditions. The analysis conducted is summarized by the following steps:


During the first two steps of the activity, the 1D engine model is refined and then validated under steady-state conditions throughout the whole engine map. The simulation code iterates the calculation by studying all the operating points in sequence. Once the results obtained match the experimental data with minor discrepancy, the simulation of transient operation of the engine can be investigated through step 3 and step 4. The 1D model is reported in Figure 3.

**Figure 3.** Gasdyn 1D schematic of the 6-cylinder, turbocharged Diesel engine investigated in this study. The "CoSim" input/output objects are clearly visible; "CoSim" input ports are represented in blue circular elements, whereas the "CoSim" output ports are highlighted in green.

Relying on the coupling functionality of the simulation tool, the 1D engine model can be exported to the Simulink® environment by means of an S-Function block. This possibility is exploited for the vehicle simulation, which is performed by means of the vehicle tool: this software is a Simulink® add-on that allows the simulation of the longitudinal dynamics of vehicles with a large variety of architectures [21]. The vehicle model consists of different sub-system blocks, each one in charge of a specific task, to represent the specific devices along the powertrain such as gearbox, electric motor, and battery as described with more details in [4].

The IC engine model operating in transient conditions is represented by the simplified scheme reported in Figure 4.

**Figure 4.** Schematic representation of the IC engine transient model.

In Figure 4 the S-Function block in light yellow contains the engine 0D/1D thermofluid dynamic model, which acts as virtual engine solver. All the parameters to be set or monitored and controlled while the engine is running are highlighted in the green and blue box.

The complete transient model receives as input the target torque and the engine speed as function of time: these are used to define the required engine parameters for the current engine speed and load by means of interpolation tables. The controllers (green blocks)

determine the actuated quantities such as injected fuel mass, wastegate valve opening and inter-cooler wall temperature, which enter the virtual engine S-Function. Finally, the S-Function block provides the desired outputs.

This complete engine model, operating in transient conditions, can be coupled to the conventional vehicle model, improving the basic representation of the IC engine by means of measured steady-state operating maps. The standard vehicle simulation tool reproduces the behavior of the different vehicle components, including the engine, by pre-determined tables and via interpolation, every output quantity requested can be instantly calculated without considering the dynamics of the device.

To introduce the engine model by an S-Function block, the pre-existent vehicle velocity control is modified as reported in Figure 5, which shows the simulation architecture.

**Figure 5.** Updated vehicle velocity control; ω: angular speed, V: vehicle speed, T: torque.

While in the conventional approach, the output torque of the engine is simply determined via interpolation of the engine maps. The proposed methodology instead asks the running engine model to produce the required torque. This implies that the output torque of the engine might not be equal to the target but reflects the fluid dynamic solution obtained, hence it is affected by the intrinsic inertia and unsteadiness of the engine. It is the actual torque produced by the engine that is fed to the vehicle dynamics.

#### **5. Calibration of Engine Model**

An important step is related to the calibration of the 1D thermo-fluid dynamic model of the engine; in particular, the set up and consequent tuning of the combustion model is fundamental to achieve reliable results.

#### *5.1. Description of Complete Engine Model*

The 6-cylinder turbocharged Diesel engine schematic is reported in Figure 6, highlighting the typical components introduced for the simulation.

Figure 6 shows the air and exhaust gas paths, including the compressor and the fixed geometry twin-entry turbine with a waste-gate valve. Six fuel injectors are present, together with dedicated controllers for intercooler outlet temperature and engine BMEP, to make the engine work in the same operating points of the experimental campaign. In the case of steady-state operating points, the engine cycles calculations are performed until convergence to the targets is reached.

It is important to underline that the model does not feature any after-treatment system (ATS) and that the transport of chemical species along the ducts was not activated: this means that the emissions could be evaluated only as cylinder-out.

The calibration procedure of the model was carried out with a detailed 1D schematic, with a very refined mesh. After that, the subsequent simulations with the complete enginevehicle system were based on the FSM approach with a coarse mesh, to achieve a significant reduction of the computational burden preserving the accuracy.

**Figure 6.** Six-cylinder turbocharged Diesel engine, 0D/1D schematic.

#### *5.2. Engine Model Calibration*

The CI engine combustion process is modeled with this approach: the cylinder volume is considered as a 0D system, in which the different chemical species are accounted for by a mixture of ideal gases. The process takes place with closed valves and is governed by the following energy equation (1st law of thermodynamics) [22]:

$$\frac{\text{V}}{\overline{\gamma} - 1} \cdot \frac{\text{dP}}{\text{d}\theta} + \frac{\overline{\gamma}}{\overline{\gamma} - 1} \cdot \text{P} \cdot \frac{\text{dV}}{\text{d}\theta} = \text{ } \text{m}\_{\text{F}} \cdot \text{LHV} \cdot \frac{\text{d}\mathbf{x}\_{\text{b}}}{\text{d}\theta} \cdot \eta\_{\text{c}} - \frac{\dot{\mathcal{Q}}\_{\text{w}}}{\omega} \tag{4}$$

where V is the instantaneous in-cylinder volume, P the in-cylinder pressure (known from the experimental pressure traces), γ the average adiabatic constant of the equivalent gas mixture, θ the crank angle, mF the injected fuel mass (per cylinder per cycle), LHV the lower heating value of the fuel, xb the "fuel burnt mass fraction" (i.e., the mass fraction of fuel that has already taken part to the combustion, at a specific crank angle position), ηc the efficiency of the combustion process, . Qw the heat loss through the cylinder walls and ω the crankshaft angular speed. In the energy balance, the unknown to be computed is xb.

In Equation (4), the left side member is the so-called "actual heat release rate" (AHRR), an energy contribution due to pressure and volume variation representing the energy effectively used to produce work. This term equals the difference between the "heat release rate" (HRR) by fuel combustion and the heat losses through the cylinder walls.

While the quantities AHRR and HRR have a known expression, the heat losses typically do not. The simple application of predictive formulas such as the Woschni model [23,24] ([6,7]) does not guarantee a reliable calculation of this term, since those models need to be recalibrated for each test case (as documented in [25,26]). Hence, by integrating Equation (4) during the combustion phase, it is possible to write that:

$$\mathbf{Q}\_{\rm w,tot} = \int\_{\theta\_i}^{\theta\_l} \frac{\dot{\mathbf{Q}}\_{\rm w}}{\omega} \mathbf{d}\theta = \mathbf{m}\_{\rm F} \cdot \text{LHV} - \int\_{\theta\_i}^{\theta\_l} \left( \frac{\mathbf{V}}{\overline{\mathbf{}} - 1} \cdot \frac{\mathbf{dP}}{\mathbf{d}\theta} + \frac{\overline{\mathbf{y}}}{\overline{\mathbf{y}} - 1} \cdot \mathbf{P} \cdot \frac{\mathbf{dV}}{\mathbf{d}\theta} \right) \mathbf{d}\theta \tag{5}$$

where Qw,tott is the total heat exchanged during combustion, by assuming the combustion efficiency ηc equal to 1. This quantity can be computed for all the engine map operating points, once the in-cylinder pressure is known by measurements and used to calibrate the Woschni model point by point. In this work the tuning has been carried out by means of the two coefficients and (for more details please refer to [23,24]) and by minimizing the following error function for each operating point:

$$z(\mathbb{C}\_1, \mathbb{C}\_2) = \left| \frac{Q\_{w, \text{tot}} - Q\_w(\mathbb{C}\_1, \mathbb{C}\_2)}{\text{CHR}} \right| \tag{6}$$

where CHR is the cumulative heat release, and *Qw* (*C*1, *C*2) is the computed wall heat exchanged for a given set of values (*C*1, *C*2). The distribution of *C*1 and *C*2 values on the engine map is shown in Figure 7. To retain a simpler set of *C*1 and *C*2 coefficients to be inserted in the 1D simulation model, their distribution on the engine map has been averaged along the BMEP axis, achieving a dependence only on engine speed.

**Figure 7.** (**a**) *C*1 distribution across the engine map; (**b**) *C*2 distribution across the engine map. *C*1 and *C*2 are the coefficients of the classical Woschni heat transfer model.

A little distortion in the distribution can be seen in some full load points: hence, these values have been reconstructed based on the adjacent points.

The total heat losses during the combustion process can be computed by applying the Woschni model, as reported in Figure 8, where the heat loss is expressed as a percentage of the total heat released by fuel combustion (i.e., simply, assuming ηc = 1).

As shown by the contour plots of Figure 8, the recalibration of Woschni's formula is fundamental to correctly evaluate the heat loss term, especially for medium-large displacements such as the one considered: without this tuning process, the heat losses would have been underestimated almost by a factor of two.

**Figure 8.** (**a**) Percentage heat losses with recalibrated Woschni model; (**b**) Percentage heat losses with original Woschni model (i.e., using the default coefficients of the model).

The fuel burnt mass fraction xb can be computed once the energy balance has been verified. The calculation is performed by integrating Equation (5), which is implicit and nonlinear, since γ depends on the chemical composition, which varies during the combustion. Once the xb profiles are computed for each sampled operating point, they must be extended across the entire engine map, since the engine could run in whichever working condition. Instead of interpolating among the calculated profiles, which are functions defined over different domains, it is advisable to approximate them by an analytical function (in this case the "double Wiebe's" one) and then interpolate the fitting coefficients of these functions. The double Wiebe's function can be written as follows (Equation (7)) [27]:

$$\alpha\_{\mathsf{b}}(\theta) = \beta \cdot \left(1 - \mathbf{e}^{-\mathsf{a}\_{1} \cdot \left(\frac{\theta - \theta\_{\mathsf{i},1}}{\theta\_{\mathsf{i},1} - \theta\_{\mathsf{i},1}}\right)^{m\_{1} + 1}}\right) + (1 - \beta) \cdot \left(1 - \mathbf{e}^{-\mathsf{a}\_{2} \cdot \left(\frac{\theta - \theta\_{\mathsf{i},2}}{\theta\_{\mathsf{i},2} - \theta\_{\mathsf{i},2}}\right)^{m\_{2} + 1}}\right) \tag{7}$$

where β is the "weight factor", a1 and a2 are the "efficiency parameters", m1 and m2 the "shape factors", θi,1 and θi,2 the "starts of combustion", <sup>θ</sup>f,1 and <sup>θ</sup>f,2 the "ends of combustion". The fitting is performed by means of MATLAB® built-in algorithms and the only parameters used are β, <sup>θ</sup>f,1, m1 and m2; moreover, the following hypotheses are done:


Once the fitting coefficients are defined for all the sampled points, their distribution across the engine map is further approximated with polynomial surfaces as function of engine speed and BMEP, to remove local distortions and have an easy computation of the coefficients for any combination of engine speeds and loads.

A comparison for the computed trends of the fuel burnt mass fraction xb is represented in Figure 9, for two operating points as examples.

**Figure 9.** (**a**) Fuel burnt mass fraction profile for the operating point at n = 800 rpm, BMEP = 12 bar. This is characterized by a double injection (a pilot and a main one), however the assumption of neglecting the pilot injection (for the fitting of the double-Wiebe function) is accurate; (**b**) Fuel burnt mass fraction profile for the operating point at n = 1700 rpm, BMEP = 8 bar (with single injection, in this case).

The profiles calculated by polynomial approximation are then considered for the remaining part of this study, given the satisfactory approximation they can provide.

#### **6. Steady-State Map Results**

After simulating all the map operating points, the predicted results can be collected and represented by means of contour plots, indexed by engine speed and load, which are used to identify the engine working conditions. These contour plots can be addressed as the steady-state maps of the IC engine, each one detailing a particular property. This kind of representation offers a quick, graphical way to verify the satisfactory correspondence between simulation results and experimental data.

As an example, Figure 10 shows a group of steady-state maps referring to the most important engine performance and emission characteristics: in particular, Figure 10a–f highlight the discrepancy between measured and calculated results in terms of percentage error, whereas Figure 10g,h directly show the measured and simulated (non-dimensional) cylinder-out NOx emissions.

**Figure 10.** *Cont.*

**Figure 10.** Steady-state maps of the six-cylinder Diesel engine predicted by *Gasdyn*: (**a**) % error—BSFC; (**b**) % error—air mass flow rate; (**c**) % error—turbocharger Speed; (**d**) % error—turbine inlet temperature; (**e**) % error—turbine inlet pressure; (**f**) % error—cylinder peak pressure; (**g**) cylinder-out NOx concentration (non-dimensional); (**h**) cylinder-out experimental NOx concentration (non-dimensional).

The prediction of specific fuel consumption by means of the 0D/1D code is very accurate at low and medium engine speeds, as reported in Figure 10a, whereas a larger discrepancy can be found at high regimes. However, considering the typical applications

of the engine investigated, that part of the engine map is not frequently exploited, so that a reliable evaluation of fuel consumption is expected during an RDE test cycle.

Considering further important engine parameters (Figure 10b–f), the percentage errors are very small and local distortions can be acceptable at the edge of the maps (typically due to interpolation issues). Finally, the calculated NOx emission map (Figure 10g) agrees with the corresponding measured trends and values (Figure 10h), represented in normalized form with respect to the maximum experimental value. The satisfactory comparison confirms the good set-up of the combustion model, being NOx formation dependent on the correct prediction of in-cylinder temperature and local air-to-fuel ratio. The most important result of this activity is that the engine model can be run across non-mapped conditions reflecting the behavior of the real engine. When the model runs in transient, interpolations are performed to set up the engine model exploiting the know values of the nearest mapped operating conditions parameters.

#### **7. Description of Vehicle and Test Cycles**

The vehicle used for the simulations is an interurban hybrid bus, which is coupled to the six-cylinder Diesel engine analyzed for this investigation. The main characteristics of this vehicle (mass, size, gearbox ratios, target mission) with a classical architecture were taken from the literature. Moreover, it is assumed that the vehicle runs on medium payload, with an overall vehicle mass equal to 12,600 kg. In the simulations, the vehicle is equipped with a six-speed manual gearbox and, according to the software gear shifting strategy, the gear shift lines are placed to allow the engine to work in its high efficiency region and to prevent it from running at excessive speeds. These lines are represented by engine torque and speed curves in Figure 11.

**Figure 11.** Gear shift lines corresponding to a six-speed manual gearbox mounted on the studied inter-urban bus, reported on the BSFC map of the six-cylinder Diesel engine.

To calculate the vehicle longitudinal dynamics, all the vehicle model blocks are provided with the requested architecture-layout data. Firstly, the simulations are carried out on the Diesel bus version with manual transmission. Then, a full hybrid parallel P2 configuration of the bus has been set up and investigated. The electric motor, placed between the engine flywheel and the gear box is a permanent magne<sup>t</sup> synchronous one. A "torque split" algorithm, together with a first-order thermal model, is used to test and validate the motor, which is simulated through steady-state maps. To match the different engine and motor

speed ranges, a reduction gear connects the latter to the drive shaft between the clutches. The system is powered by a lithium-ion battery pack, characterized by a 1p110s (1 parallel, 110 series) arrangemen<sup>t</sup> working at 400 V, which has been designed taking into account the peak-power request, calculated by means of the "torque split" algorithm adopted, and of the overall energy storage requirements.

Once the electric components are chosen, a proper tuning of the hybrid control unit (HCU) is performed in the vehicle simulation tool, acting on the specific torque-speed curves that define the operating mode of the hybrid system, and the battery SoC (state of charge) threshold levels [28]. As mentioned, a first set of parameters to be defined is represented by the torque curves, which are indexed by the revolution speed, generalized to the drive shaft between the clutches. These curves characterize different regions on a torque-speed plane: depending on the total torque required by the vehicle and its actual speed, the position of such working condition on the torque-speed plane determines the hybrid system operating mode. For the hybrid architecture considered, the torque-speed curves are referred to the engine BSFC map: this allows their correct positioning, aimed at reaching the best overall efficiency.

In Figure 12, it is possible to distinguish four different curves, cited with their actual name in the vehicle simulation tool. To correctly place these curves, the underlying rationale of the hybrid powertrain must be pursued: letting the engine work in its highefficiency regions.

**Figure 12.** Position of the HCU-Parallel block curves on the engine BSFC map.

The lowest purple line defines a map region below which the engine is generally prevented from working: when the requested engine torque lies below this curve, the HCU determines an increase of torque request, so as to bring the operating point above the purple line. Then, the torque surplus is used to charge the battery through the motor generator mode. This curve has been placed to achieve a compromise: it prevents the engine from working in low-efficiency points, such as the bottom part of the map, so that above this curve an optimal operation of the engine is achieved.

The region between the purple line and the yellow one defines the area in which the engine can work. If the required torque is located between the yellow line and the red one the engine is supported by the electric motor to provide the requested torque. The blue curve represents the maximum efficiency curve for the engine, to achieve the minimum fuel consumption and CO2 production. It is important to remark that, due to the characteristics of this engine, the maximum efficiency is often located around the full load

in the engine map. If the overall requested torque exceeds the full load torque of the engine, the additional torque requested is provided by the electric motor. The vehicle control unit also decides when to turn off and on the engine. It is important to highlight that the engine model has the capability to be turned off and on, with the corresponding effects on the fluid dynamic solution.

To evaluate the bus performance under various scenarios, two different missions were chosen. The first is the orange county (OC) cycle, a chassis test for heavy-duty vehicles such as buses. Its speed profile is depicted in Figure 13a: with relatively low speeds, but frequent and steep accelerations and decelerations. This cycle is suitable to test the bus in urban driving conditions. The second one is the world harmonized vehicle cycle (WHVC), a generic test for heavy-duty vehicles, usually used for research purposes to compare the vehicle and engine emissions, derived from its corresponding engine dynamometer cycle: the WHTC. By looking at Figure 13b, this driving cycle consists of three different phases: an urban one, from the start of the cycle until the red dotted line; a rural phase, limited by the red and green dashed lines; a motorway phase, in the last part of the cycle. The composition of the cycle makes it suitable to test the bus in various conditions, including travelling on extra-urban roads, during the rural and motorway phases. Further information about these cycles can be found in [29,30].

**Figure 13.** (**a**) Orange County cycle—speed profile; (**b**) WHVC—speed profile.

#### **8. Vehicle Dynamics Simulations**

The application and validation of this new co-simulation procedure are based on a direct comparison between two different approaches. In this section, for both bus versions (conventional Diesel and hybrid), the main outcomes of the map-based approach and the transient coupled model are reported and discussed. Only the results related to the WHVC cycle are shown since it represents a mixed mission able to test the new methodology into a wide range of accelerations and speeds, even if the full exploitation of a hybrid system can be reached in urban driving.

The fuel economy series for the different Diesel simulation models are collected in Table 3. The increase of fuel consumption for the transient model is probably due to the dynamics of the engine (inertia is considered with respect to the map-based steady-state model) as well as to a little overestimation of fuel consumption, as shown by Figure 10a and by the midway model (based on the map of fuel consumption, obtained by the analysis of steady-state operating conditions).

**Table 3.** WHVC: Diesel bus fuel economy comparison.


1 Map-Based approach with fuel consumption map from experimental data. 2 Map-Based approach with fuel consumption map from engine model steady-state analysis. 3 Vehicle-engine transient model.

To carry out a fair comparison between the fuel consumption of the Diesel and hybrid versions, a normalization technique must be applied. The procedure devised by the "low carbon vehicle partnership" (LowCVP) and adopted in the United Kingdom is described in detail in [29]: intended just for buses, it perfectly fits the case considered in this work. According to this procedure, it is necessary to calculate a quantity called "net energy change" (NEC), expressed by the following formula:

$$\text{NEC} = (\text{SoC}\_{\text{f}} - \text{SoC}\_{\text{i}}) \cdot \text{V}\_{\text{sys}} \cdot \text{K}\_{\text{1}} \tag{8}$$

where SoCf and SoCi are respectively the final and initial battery state of charge expressed in [Ah] (calculated multiplying the battery capacity by the percentage state of charge), Vsys is the nominal DC operating voltage of the battery, and K1 is a constant equal to 3600 s/h, necessary to pass from hours to seconds. Besides the NEC, two further quantities must be calculated. These are the "total fuel energy" (TFE) and the "total cycle energy" (TCE):

$$\text{TFE} = \mathbf{m}\_{\text{F,tot}} \cdot \text{LFV} \tag{9}$$

$$\text{TCE} = \text{TFE} - \text{NEC} \tag{10}$$

where mF,tot is the total fuel mass consumed. The last formula (10) computes the total energy used to complete the driving cycle, provided by the fuel combustion and the battery. The final step is to calculate the "NEC variance" according to Equation (11) and check which range of Table 4 it falls into.

$$\left| \text{NEC Variance} \right| = \left| \frac{\text{NEC}}{\text{TCE}} \cdot 100 \right| \tag{11}$$

**Table 4.** Possible |NEC Variance| ranges.


The correction requires to simulate the driving cycle multiple times, until getting three valid runs at least (i.e., each one with a |NEC Variance| ≤ 5%). Once this is achieved, each test run must be identified by an NEC Variance value and a fuel consumption; these points are then represented in a plane and processed with a linear regression (see Figure 14). If the "coefficient of linear correlation" of the regression is R<sup>2</sup> ≥ 0.8, then the set of points constitutes a family of valid test runs and the corrected fuel consumption is given by the linear regression straight line, evaluated at NEC Variance = 0%. Hence, as indicated by the normalization procedure, three different SoCi are chosen: 35%, 50% and 65% for both the map-based experimental procedure and the transient simulation. Note that all six runs are valid (i.e., |NEC Variance| ≤ 5%).

**Figure 14.** Calculation of WHVC fuel consumption by different models.

The normalized fuel consumption, calculated for the hybrid vehicle by both simulation approaches, is collected in Table 5 and directly compared to that for the conventional Diesel vehicle. As already discussed, the transient simulation model reveals slightly higher fuel consumption.



1 Map-Based approach with fuel consumption from experimental data. 2 Vehicle with engine transient model.

In Figure 15a, the resulting cloud of engine operating points during the WHVC test cycle with a SoCi of 50%, simulated by the transient approach, are shown together with the HCU curves of Figure 12. The engine operation follows the rules of the logic previously described: it works very rarely where the engine is not efficient, except when the engine is off or it is turned on trying to the target operating condition. Figure 15b shows the SoC of the battery during the same run, calculated by both simulation methods. The two curves are similar, but not identical: this depends on how the HCU manages and uses the two sources of torque (i.e., the e-motor and the IC engine). It bases the splitting on the required torque at that time instant: since the requested torque can be a little bit different between the two modeling approaches (due to the influence of engine dynamics), a particular engine condition can fall into different HCU regions on the map. Therefore, the managing of the required torque can be different, as well as the SoC profile.

**Figure 15.** WHVC with SoCi = 50%, simulation results: (**a**) cloud of engine operating points calculated by the transient *Gasdyn* approach; (**b**) battery SoC.

Moreover, the SoC at the end of simulation (SoCi) reaches half of the battery storage energy. This convergence condition can be noticed for both the driving cycles and the three SoCi values chosen, emphasizing the good design of the HCU logic.

Even if the WHVC is a mixed driving cycle, significant fuel savings can be achieved, probably mostly during its first phase: the urban drive. In fact, by comparing the Diesel and hybrid map-based calculations, the fuel economy improves by 15.96%, while for the hybrid transient simulation approach the reduction of fuel economy is equal to 16.16% (with respect to the Diesel transient model). Overall, a satisfactory agreemen<sup>t</sup> can be found, when comparing the two different powertrain architectures with the same modeling approach: independently of how the engine operation is calculated (through the map-based or transient method), the HCU is able to manage both, highlighting the advantages of the bus powertrain investigated in this work: the hybrid architecture.

Dealing with the simulation of RDE and WHTC test cycles, it is important to discuss the computational effort requested and the prediction accuracy achieved by the FSM (fast simulation method) methodology, compared to the standard solver applied to a 1D schematic with a refined mesh. This can be analyzed in terms of CPU time and of variation of engine output quantities. To this purpose, the simulation of the same WHTC cycle has been performed with both a refined mesh (1 cm) and a coarse one (10 cm with the FSM methodology). As shown in Figure 16a,b the reduction of computational effort is relevant, around 80% with respect to the same simulation in the refined mesh solution. To estimate if the accuracy has changed between the two approached, the cumulative fuel consumption and cylinder out emission are compared. As it can be seen, the two approaches differ by less than 1%, revealing that the overall accuracy has been preserved and that long-lasting driving cycle can be simulated by resorting to fast simulation strategies without introducing excessive approximations.

**Figure 16.** (**a**) Comparison between CPU time/real time of the driving cycle simulated by means of both refined mesh (1 cm) and coarse mesh (10 cm, FSM); (**b**) difference in the prediction of cumulative emissions of the engine emissions between the refined mesh case and FSM.
