This is a self-archived – parallel published version of this article in the publication archive of the University of Vaasa. It might differ from the original. Toward a digital twin of a mid-speed marine engine: From detailed 1D engine model to real-time implementation on a target platform Author(s): Hautala, Saana; Mikulski, Maciej; Söderäng, Emma; Storm, Xiaoguo; Niemi, Seppo Title: Toward a digital twin of a mid-speed marine engine: From detailed 1D engine model to real-time implementation on a target platform Year: 2023 Version: Accepted Manuscript Copyright ©2023 IMechE. Published by Sage Publications. The article is protected by copyright and reuse is restricted to non-commercial and no derivative uses. Users may also download and save a local copy of an article accessed in an institutional repository for the user's personal reference. Article reuse guidelines: sagepub.com/journals-permissions, DOI: 10.1177/14680874221106168 journals.sagepub.com/home/jer Please cite the original version: Hautala, S., Mikulski, M., Söderäng, E., Storm, X. & Niemi, S. (2023). Toward a digital twin of a mid-speed marine engine: From detailed 1D engine model to real-time implementation on a target platform. International Journal of Engine Research, 24(12), 4553-4571. https://doi.org/10.1177/14680874221106168 Towards a digital twin of a mid-speed marine engine: from detailed 1D engine model to real-time implementation on a target platform Saana Hautala 1)*, Maciej Mikulski 1)*, Emma Söderäng 1), Xiaoguo Storm 1), Seppo Niemi 1) 1) University of Vaasa, School of Technology and Innovation, Wolffintie 34, FI-65200 Vaasa, Finland * Corresponding author: saana.hautala@uwasa.fi, maciej.mikulski@uwasa.fi Abstract System complexity is challenging for development of marine mid-speed engines when striving to meet increasingly stringent emission targets. Control-oriented modeling offers a solution, cutting calibration time and enabling robust control strategies. Simultaneously, real-time, physics-based engine models (digital twins) are emerging as they offer better predictive capability and scalability than typical mean-value, data-driven approaches. This study explored development of a control-oriented digital twin of a Wärtsilä 4L20 marine engine. Starting from a detailed one-dimensional model (GT-Suite), it explored reduction strategies towards a fast-running engine model (FRM), balancing the calculation speed and accuracy trade-off. Finally, the FRM was tested for real-time implementation on a target machine. Comprehensive experimental data from the 4L20 platform in the VEBIC Engine Laboratory provided the baseline for model calibration. Model calibration and validation covered four representative operating points and involved correlation of crank-angle, resolved in- cylinder pressures, thermal state at several locations of the engine air-path and relevant performance indicators. The results shed new light on the feasibility of digital twins in the marine engine domain. The obtained FRM was three times faster than real-time, while the accuracy loss was comfortably within the 5% tolerance levels for the governing outputs, including crank angle resolved in-cylinder pressure. The grid-resolved simulation was obtained with 4 times fewer flow components and internal discretization length of 100% and 150% of the cylinder bore for intake and exhaust components respectively. The balance between predictivity, accuracy and real-time surplus, was ultimately more favourable than in state of the art automotive applications and enables exploring further coupling with semi-predictive emission submodels. The real-time capable FRM is considered applicable in hardware-in-the-loop simulation, and this application is scheduled in a follow-up project. mailto:saana.hautala@uwasa.fi mailto:maciej.mikulski@uwasa.fi Keywords: marine diesel engine, mid-speed engine, 1D-simulation, fast-running engine model, digital twin, engine control Nomenclature 1D one-dimensional BMEP brake mean effective pressure BSFC brake specific fuel consumption BTE brake thermal efficiency CA50 crank angle at 50% fuel mass fraction burned CAC charge-air-cooler Cd discharge coefficients CA deg crank angle degree CO carbon monoxide CO2 carbon dioxide CPU central processing unit DF dual-fuel ECS engine control system EEDI Energy Efficiency Design Index EEOI Energy Efficiency Operational Index FRM fast-running engine model GHG greenhouse gas HiL hardware-in-the-loop HRR heat release rate ICE internal combustion engine IMEP indicated mean effective pressure IMO International Maritime Organization MVM mean value model NOX nitrogen oxides pX pressure in location X RCCI reactivity controlled compression ignition RT real-time SOI start of injection SOX sulfur oxides TX temperature in location X TPA three pressure analysis UHC unburned hydrocarbon VEBIC Vaasa Energy Business Innovation Centre 1. Introduction The global environmental crisis and increasingly stringent emission legislation are driving internal combustion engine (ICE) research towards ultra-clean and efficient powertrains. Electrification is making rapid progress in the light-duty transport sector, but applications requiring high energy density, like off-shore shipping, will continue to rely on the ICE. According to the third International Maritime Organization’s (IMO) Greenhouse Gas Study [1], shipping is responsible for about 2.5% of global greenhouse gas (GHG) emissions and significantly contributes to local air pollution in coastal and harbor areas. There are several national and international regulations for lowering waterborne transport emissions including nitrogen oxides (NOX), sulfur oxides (SOX), carbon monoxide (CO) and unburned hydrocarbons (UHC). Propulsion and ship-level efficiency indexes, including the Energy Efficiency Design Index (EEDI) and the Energy Efficiency Operational Indicator (EEOI) also serve to control GHG emissions. A good overview of current waterborne emission legislation can be found in the IMO Marine Engine Regulations [2]. System complexity becomes a limiting factor for the development of marine mid-speed engines needing to meet tough emission targets. Physics-based modeling and simulation addresses this issue of complexity, drastically cutting the time and effort needed to develop new hardware and robust control strategies. Among various modeling approaches, one-dimensional (1D) engine simulation toolchains [3‒6] are gaining traction for many applications because of their favourable trade-off between accuracy and calculation time. This trend is recognized in the marine domain. Table 1 reviews the most recent applications of control-oriented models in marine domain, putting them in the context of predictive capabilities, accuracy and simulation speed. Automotive state of the art in this field is included for reference along with the characterisation of the present work. Table 1. Control-oriented models in marine engine development. Color coding: green – strongest in category; yellow – middle in category; red – weakest in category; grey – out of scope/not feasible Reference Engine Type Fidelity/predictivity Accuracy Simualtion Speed HiL application airpath combustion Mikulski et al. (2019) [7] MAN 6L35/44DF 1D-detail; fully predictive 0D - multizone fully predictive; Validated on single point; 5% error ~1000 times slower than real-time Not feasible Stoumpos et al. (2020) [10] Wärtsilä 9L50DF engine 1D-detail; fully predictive 0D - two-zone; semi-predictive Not validated 20-25 times slower than real-time Not feasible Perabo & Zadeh (2020) [11] Generic marine diesel engine Data driven; non predictive Data driven; non predictive Not validated ~100 times faster than real-time HiL simulation complete ship propuslion Llamas & Eriksson (2019) [13] MAN 6S80ME-C9.2 0D-MVM; semi-predictive Anatlytical model; semi-predictive Full factorial validation; 4% error 50 times faster than real-time Not demonstrated Wu & Li (2016) [22] Automotive state of the art 1D-FRM; fully predictive 0D; single zone; semi-predictive Full factorial validation Close to real time Full HiL tests; engine control Hautala et al. (2022) [This work] Wärtsilä 4L20 1D-FRM; fully predictive 0D - two-zone; semi-predictive Full factorial validation; 5% error 3 times faster than real-time Demonstarted on a target machine According to Table 1, on the high-end of control-oriented model applications in marine domain is the approach by Mikulski et al. [7]. The authors studied reactivity controlled compression ignition (RCCI) in a medium-speed MAN 6L35/44DF marine engine. In this work, a GT-Power model, was coupled with a predictive, chemical kinetics-based multi-zone combustion model. This coupling allowed extrapolation of engine performance and emissions to RCCI combustion conditions. Large-scale multi-processor cluster simulations were performed to seek optimum efficiency within certain constraints (NOX, CO, fuel blend ratio). The results confirmed the operational advantages and challenges of the RCCI concept compared to conventional dual-fuel (DF) operation in large-bore engines. Based on the simulations, the authors concluded that RCCI implementation in existing medium-speed ICE systems does not require large-scale hardware change. Mavrelos et al. [8] and Stoumpos et al. [9] developed a detailed air-path model of a large marine dual- fuel engine (Wärtsilä 9L50DF), including semi-predictive combustion and emission models. The authors used it to investigate steady-state performance of both diesel and gas operation with respect to efficiency and emissions. Additionally, parametric studies were carried out to determine pathways for simultaneous NOX and CO2 reduction during DF operation. More recently, Stoumpos et al. [10] combined their modeling approach with an in-house engine control system (ECS) model, aiming to achieve a digital twin. Although the engine and ECS were modeled to a detailed level, based on manufacturer´s specifications, they ultimately lacked experimental validation. Furthermore, actual digital twin implementation was hindered by relatively long simulation times, attributed to the model´s level of detail. The authors reported computational times approximately 20‒25 times greater than real-time requirements. As the above studies show, 1D models are efficient tools for studying novel engine technologies and can be used for control development, assuming real-time simulation is not required. For the second application, usually non-predictive data-driven approaches were used. In the marine domain, for instance, Perabo and Zadeh [11] represented the main diesel propulsion engine as a transfer function for hardware-in-the-loop (HiL) simulation of a complete shipboard power system. This type of simple engine model couples well with electric component submodels, including local controllers like voltage regulators and speed governors. Models of this type are far faster than real-time and provide good results in system-level simulation, as described above. However, they ultimately are too simplistic and lack predictivity for dedicated engine control applications, or when considering coordinated power-management–engine control. Dedicated model-based engine-control applications are currently more prevalent in the automotive sector [12]. Traditionally, many studies focused on engine control used zero-dimensional, mean-value engine models to develop controls. Such models are, in principle, physics-based but, being combustion cycle oriented, completely disregard high-frequency phenomena like system pulsation or combustion dynamics. Llamas and Eriksson [13] have recently used such an approach to model a large, two-stroke marine engine with exhaust gas recirculation, applying it in vessel-level simulation with promising results in terms of accuracy and speed. Nowadays, however, a more tailored modeling approach is necessary. Combustion engines are becoming more complex phenomenologically and advanced combustion concepts are being introduced [14‒ 16], so coordinated system controls must extend beyond conventional applications. By means of selective model reduction strategies, one can maintain the benefits of 1D framework and still make the model real-time- capable. For example, a fast-running engine model (FRM) allows focus on details in the areas of the model that are restrictive for a particular application, while gradually introducing simplifying assumptions in parts of the model that are less critical for the given simulation approach. The model reduction process lasts until the target simulation speed is obtained. An FRM can be used in system-level simulations where long transient events need to be simulated, such as those involving waste-heat-recovery [17] or driving-cycle simulation comparing multiple engine / vehicle configurations [18]. Successful FRM implementation in HiL simulation, where fast execution time is critical, were recently demonstrated in the automotive sector by amongst others: Wei et al. [19] and Andric et al. [20, 21]. The state of the art in the achieved trade-off between simulation time and fidelity in this sector can be well illustrated based on the work of Wu and Li [22], referenced in Table 1. The authors developed a real-time-capable (with narrow margins) FRM of a high-speed diesel engine, by iterative reduction of the detailed 1D model. The effects on engine parameters were compared between rapid control prototyping (FRM on HiL) and GT-Power controllers used in the detailed engine model. It was verified that the HiL simulations correlated well with the detailed model results. The simplifications caused some differences in the predicted fuel value but this was not a critical concern in the context of the HiL application, so the model was considered fit for purpose. Such HiL-implemented, physics-based, predictive models are commonly referred to as digital twins since they can replicate the engine operation (including non-nominal conditions) in real-time. In future, digital twin technology will enable new functionalities, like extensive condition-based monitoring, and provide a route towards self-learning, integrated powertrains. This is expected ultimately to defeat the system complexity issue, which is the main bottleneck for advanced engine development. As evident from Table 1, engine digital twin technology is already quite mature for high-speed automotive engines, but this is not so in the marine sector. The original equipment manufacturers in marine propulsion are making progress in this technology, but digital twin functionality for large-bore engines has not been openly demonstrated. One should note that the techniques for model reduction towards real-time-capable surrogates are not directly scalable from smaller engines because typical control applications differ. Furthermore, the scale of elements influences numerical capabilities of the solvers (refer to the so-called Courant condition [23]). The present work builds on the above-mentioned knowledge gaps in the application of digital twins for marine engines. Specifically, it aims to answer three research questions: (i) How can large volumes specific to marine engines be discretized to capture system dynamics (including pulsations) with high calculation speed? (ii) What mechanisms are responsible for accuracy loss when transferring the code to FRM and further towards a real-time license? (iii) What are the enablers and limiting factors of the developed model in typical engine control applications? Finally, note that the majority of the published works that consider marine engine modeling lack validation via real test data. This shortcoming is a direct result of the difficulty in accessing large-bore engine facilities and of the highly competitive environment. The present work excels in this respect, for the first time providing a complete and thoroughly validated modeling framework, including a detailed 1D marine diesel engine model and its FRM surrogate tested for real-time capabilities on a target platform. The main merits of our modeling approach are included in Table 1 for reference. 2. Materials and methods 2.1 The research object The object of this research was a Wärtsilä 4L20 diesel engine. The detailed technical parameters of the engine were used to identify the model and the engine test results were used for model calibration and validation. Table 2 provides the engine´s basic specification: a detailed specification is in the Wärtsilä 4L20 Project Guide [24]. The engine was turbocharged and equipped with a water-cooled, charge-air-cooler (CAC). The fuel used was EN ISO 590 standard, sulfur-free diesel, injected via centrally mounted nine-hole solenoid injectors into the engine´s “Mexican-hat” combustion chambers. The engine was connected to a squirrel-cage induction motor, which acted as a generator and was used to load the engine in the laboratory setup. The full load for the engine was 848 kW, providing 800 kW on the generator. Table 2. Specification of the Wärtsilä 4L20 research engine Cylinder configuration four, in line Bore 200 mm Stroke 280 mm Swept volume / cylinder 8.8 dm3 Compression ratio 16:1 Rated Speed 1000 rpm Brake power 848 kW Fuel system Common rail Injector Solenoid / nine-hole axisymmetric / 153º umbrella angle Turbocharger ABB TPS48E01 Valve system four valves/cylinder, Miller timing- capable Exhaust Valve Opening / Closing (EVO/EVC) 139 / 386 CA deg * Intake Valve Opening / Closing (IVO/IVC) 334 / 507 CA deg * * valve timings are provided in CA scale from -180 to 540 CA (0 CA is TDC firing) as used in Fig. 8 The engine was calibrated for constant-speed operation at 1000 rpm. For model validation, measurements were performed at four representative load points of 100, 75, 50 and 25% of the full load. The measurements included temperature and pressure readings taken at several locations in the engine air-path, indicated in Fig 1. Note that care has been taken to provide the data necessary to calibrate and validate different model components. Thus, the measurement locations included the inlet/outlets of CAC, turbine and compressor respectively. Air-path parameters were determined as time-averaged values (60 s recording window) for a given steady-state operating point. Engine air mass flow rate was not measured directly, but calculated from the thermodynamic state at the charge air inlet, according to Eq. (1): 𝑚̇𝑎𝑖𝑟[𝑘𝑔/𝑠] = 𝑘 ∙ √(𝜌𝑎𝑖𝑟 ∙ ∆𝜌𝑛𝑜𝑧𝑧𝑙𝑒 ∙ 10.19) ∙ 0.985 (Eq .1) where k is coefficient for air measurement nozzle, ρair is air density (for ambient conditions), and ∆pnozzle is pressure drop over the air measurement nozzle. Since k was calculated on mmH2O, the conversion factor to mbar for ∆pnozzle was 10.19. Hypothetical sealing air losses from the turbocharger were 1.5%, thus 0.985 factor was used. To this end, the maximum uncertainty of ṁair, and other indirectly measured parameters, were calculated using the partial divertive method, following the original work by Kline and McClintock [25]. For directly measured quantities, either device accuracy or standard deviation, whichever was higher, was used for measurement uncertainty. For the sake of brevity, the detailed uncertainty analysis is not discussed in the paper. However, measurement uncertainties are mentioned in the model validation section wherever differences between model and measurement results were close to the level of significance, thus influencing the interpretation. Fig. 1. High-level diagram of the Wärtsilä 4L20 research engine and measurement points. Subscripts at corresponding pressures (p) and temperatures (T) denote respective locations: 0 – ambient; 1 – compressor inlet ; 2 – compressor outlet / CAC inlet; 3 – CAC outlet / intake manifold; 4 – in-cylinder; 5 – turbine inlet; 6 – turbine outlet. Other measurements for model identification/validation included the fuel consumption (gravimetrically) and instantaneous in-cylinder pressure at all cylinders. Kistler 6125C piezoelectric transducers, connected via charge amplifiers to the Kistler KiBox acquisition module, measured cylinder pressures, and the high- frequency signal was sampled with an angular resolution of 0.1 crank angle degree (CA deg). The pressure signal was automatically pegged, filtered and cycle-averaged (100 cycles) by the KiBox module. The real-time pegging for all in-cylinder pressures was performed using the intake manifold pressure (p3) signal values at respective cylinder bottom dead centres / intake valve opening (IVO) coincidences. Such conditioned cycle- averaged in-cylinder pressures were further post-processed to obtain burn rates and related combustion indicators, like the crank angle at 50% fuel mass fraction burned (CA50). To provide the highest degree of accuracy for those essential model inputs, the detail three pressure analysis (TPA) method was used for this purpose. The method utilizes a 1D gas exchange model of the representative cylinder, and as such excels over the standard second-law thermodynamic analysis, as it enables accurate estimation of cylinder-specific mixture conditions and disregards errors related to heat loss estimation [26]. Note that the test bench instrumentation included many other measurement devices, including several analyzers for emissions and detailed monitoring of the engine media and the generator circuit. These are not discussed here, as they are not pertinent to the scope of the present study. Fig. 2 shows the complete engine test rig. Fig. 2. Laboratory setup for the Wärtsilä 4L20 engine in the VEBIC engine laboratory 2.2 The modeling workflow The route towards the real-time capable digital twin explored in this research followed the workflow depicted in Fig. 3. The detailed 1D engine model, forming the basis for this development, is introduced in subsection 2.3. The reduced-order FRM and its simplifying assumptions are then discussed in subsection 2.4. The methods applied to calibrate and validate the models are commonly discussed for detailed 1D model and FRM in subsection 2.5. Finally, the real-time methods used in this research are outlined in subsection 2.6. Fig. 3. Process flow chart presenting the development of a real-time capable engine model. 2.3 The 1D model The 1D model was originally created by Wärtsilä for a stock 200 mm-bore engine, and further modified in this research to fit the particularities of the 4L20 research engine in Vaasa Energy Business Innovation Centre (VEBIC). The model block diagram is provided in Fig. 4 and assumptions regarding the discretization and solution methodology are provided in Table 3. Fig. 4. Overview of detailed 1D model representing the Wärtsilä 4L20 engine The model shown in Fig. 4 closely reproduces the geometry of the actual air-path. To this end, the system has been discretized by creating an individual flow element at every distinctive bend, contraction and flow split, according to the GT-Suite Flow Theory Manual [27]. The flow elements have been further identified using exact 3D-CAD drawings of the elements. Aiming for high accuracy at sensitive locations, (high flow velocities) the default GT-Suite calculation of discharge coefficients (Cd) was exchanged for dedicated identification measurements on the flow bench. This particularly concerned the CAC and cylinder ports, along with valves. For the latter, the Cd was identified at relevant points of the valve-lift profile, independently for forward and reverse flows. Achieving high accuracy was supported by proper model discretization and adequate selection of solver (see Tab. 3 for details). Further distinctive assumptions on the air-path side included a predictive wall heat transfer model on both the intake and exhaust side, and a semi-predictive model for the CAC, where output temperatures were correlated to individual load conditions. The turbine and compressor submodels were standard GT-Power map-based approaches, using maps provided by the turbocharger manufacturer. Table 3. Assumptions for the detailed 1D model and its FRM surrogate Feature 1D detailed FRM Minimum discretization length intake 50 mm 200 mm Minimum discretization length exhaust 50 mm 300 mm Number of flow components 181 47 Solver explicit, forward Runge-Kutta method explicit, forward Runge-Kutta method / explicit Euler real- time (in the RT license) Maximum simulation time step * 0.00017 s 0.12 s Average simulation time ** 77 s per steady-state case 19 s per steady-state case (4 s with a real-time license) * user-imposed limit ** on single CPU (Intel i7-8750H 2.20 GHz) with 16 GB RAM To prevent the uncertainties of reproducing combustion from affecting overall model accuracy, burn rates were determined from experimental data (mentioned TPA analysis) and imposed directly in the model. The model is semi-predictive in terms of combustion since the cumulative burn rate at imposed distinctive operating points (100, 75, 50 and 25% load) is corralled to the fuel value (model input). This allows simulating intermediate load conditions and finally provides load-transient simulation capability. Note that, to avoid additional computational load to the targeted FRM, at this stage, the in-cylinder model does not involve any emission submodels. 2.4 Fast-running engine model For real-time applications, the execution time of the detailed Wärtsilä 4L20 GT-Power model needed to be reduced. The calculation time of the detailed 1D model was approximately six times longer than real-time (refer to Tab. 3 for details on simulation time). The approach of building a real-time surrogate relied on incremental and target-oriented simplification of the detailed model, leading to a so-called fast-running engine model (FRM). This methodology prevails over a commonly used map-based mean value model (MVM) approach because it enables more physics, hence resulting in a more predictive model. The GT-Suite FRM converter tool was used to support the model reduction process. In the FRM conversion, the predefined subsystems were simplified in terms of resolution and/or submodel accuracy. The macro discretization level was reduced by aggregating individual flow elements into subsystems. The selection of elements to be modeled independently was determined by simulation speed and the relevance of pressure pulsation reproduction. For instance, the intake manifold, where pulsations largely shape individual cylinder aspiration efficiencies, was still modeled in the FRM with individual runners. The corresponding elements on the exhaust side were aggregated into a single volume, yet twin-turbine entry needed for capturing turbocharger efficiency was maintained. The overall process reduced the number of individual flow components from 181 in the detailed model to 47 (Tab. 3). Fig. 5 depicts the resultant FRM model assembly: comparison with Fig. 4 underscores the degree of simplification made during the FRM conversion. The individual steps 1-3 highlighted in the figure, were chosen arbitrarily, to allow analysing the effect of incremental reduction on model accuracy and performance. Those are referred in section 3.3. Fig. 5. Overview of the 1D FRM representing the Wärtsilä 4L20 engine including main FRM simplification steps. Reduced amount of flow components in each step together with improvement in average simulation speed in comparison to detailed 1D engine model are presented. All created subsystems, i.e., exhaust manifold, exhaust pipes, intake manifold, compressor outlet pipes (including CAC), and intake pipes were simplified in terms of accuracy, by moving from predictive, very detailed submodels of friction and heat loss to simpler approaches. This strategy was used only if it was supported by sensitivity analysis. For example, aggregated elements of the exhaust manifold were left modeled with a predictive wall thermal solver, as it was by far the most sensitive to the accuracy of wall temperature estimation. Other elements were modeled with wall temperature imposed as an average value from detailed 1D simulations. Since the imposed combustion profile was already used in the detailed model, simplifications https://www.sanakirja.fi/english-finnish/arbitrarily for combustion or cylinder elements in general were not required. The same applied for the turbocharger model. In the FRM, the trade-off between model accuracy and calculation time can be tailored by appropriate selection of the internal discretization lengths. This was initially set for the exhaust and intake sides to 300 mm and 200 mm respectively. The values were, thus, roughly 4‒6 times greater than the detailed model, altogether with macro-discretization measures, resulting in a 74% reduction in the number of flow components to be solved for the FRM (see Tab. 3 for absolute quantities). The effect of discretization length on accuracy/calculation speed is further analyzed in the Results section. The FRM shared the same explicit solver as the detailed model, meaning that the calculation step size was finally dependent on the discretization length through the Courant condition [23]. This was additionally constrained in the detailed model by setting the maximum time step to the equivalent of 1 CA deg. In the FRM this accuracy-enhancing constraint was waived to support fast run times. Instead, the equivalent of a full engine cycle was imposed as maximum time step. Table 3 summarizes the main assumptions and resulting performance indicators of both the detailed model and its FRM surrogate. 2.5 Model calibration and validation Section 3, Results, is the focal point of this paper, demonstrating the final validity of the models. In order to achieve the accuracy level discussed there, the models had to be thoroughly calibrated. At this point, the reader should bear in mind that model calibration is a complex process, so the ensuing description merely highlights the considerations and the principle steps taken. The short discussion aims to support the credibility of the approach. To decouple the complex calibration problem for the detailed 1D model, component calibration was performed before tackling full model simulations. The CAC, turbocharger and cylinder components were calibrated separately against measurement data. The components were decoupled from the engine model at the points where temperatures and pressures had been measured in the laboratory (see Fig. 1 for reference). The measured values were then used as boundary conditions for the component models. The calibration of the CAC was straightforward, entailing manual tuning of the intercooler efficiency map to match the output temperatures for all measured load points. For turbocharger elements to match the experimental output, the considered calibration factors were pressure, temperature and mass flow rates, efficiency and mass flow multipliers on turbine and compressor. This process involved knowledge-driven sensitivity analysis: in fact, only the efficiency multiplier (affects mainly respective output temperature [28]) on the turbine side required tuning. Other outputs were considered within the adopted tolerances. Cylinder calibration is a complex and time-consuming process. Briefly, it entailed first matching the simulated (TPA) and experimental in-cylinder pressures. Cylinder-to-cylinder variations (standard deviations at selected points of the pressure trace) of the experimental data were considered as calibration tolerances. This was used to obtain the accurate burn rate for individual cylinders. With the burn rate imposed on the actual cylinder model, the second stage matched synthetic cylinder outputs with experimental data, namely; brake power, brake mean effective pressure (BMEP), brake thermal efficiency (BTE), maximum cylinder pressure, air and fuel mass flow rates, brake specific fuel consumption (BSFC) and volumetric efficiency. In both stages, intake valve lash, port wall temperature (heat flow in ports) and intake valves’ flow area multiplier were considered for tuning as the parameters that were not specifically identified in the model, thus subjected to biggest uncertainty. The tuning process was conducted repetitively until convergence between the reverse run model (TPA) and the forward run model was reached. This was done for all operating points to find a global (case independent) optimum tuning. For the full engine model, only the injection mass and timing, together with ambient boundary conditions, were used as inputs. The flow characteristics of the complete model differed from the decoupled submodels due to the cross-transfer of their respective tolerances. Thus, in the final stage, the effects of the tuning parameters were reviewed and fine-tuning on the whole model was performed. In addition to the tuning parameters already considered, the exhaust outlet pressure (p6) was optimized (via the inbuilt GT-Suite optimizer) to calibrate the exhaust temperatures before (T5) and after the turbine (T6). Calibration of the FRM was embedded in the reduction process and performed after individual reduction steps. Calibration was against the validated 1D model simulation results. Since the FRM reduction for larger volumes affects mainly subsystem capability to predict pressure drops and output temperatures, respective orifice diameters or heat transfer multipliers were used to calibrate the model on the subsystem level. Calibration was performed for the highest load point only since the simplifications manifest most strongly in this regime [27]. The optimized tuning parameter value from the calibration was used for validating the full model at all operating points. The key variables considered under final FRM validation are presented in Table 4, together with the targeted tolerance level. Additionally, the factor of real-time, which presents the model’s capability to perform in real-time, was under investigation throughout the process. Table 4. FRM key variables and tolerances Variable Units Tolerance type Tolerance BMEP bar absolute 0.5 BTE % percentage 5.0 Volumetric efficiency fraction percentage 5.0 Average max. cylinder pressure bar percentage 5.0 Intake manifold pressure bar percentage 5.0 Exhaust temperature at turbine inlet K percentage 5.0 2.6 Real-time implementation During the FRM conversion process, the calculation time was decreased from roughly six times real- time to one to two times real-time (air-path response only; timescales were much more demanding for in- cylinder processes). However, additional steps were needed to achieve capability for real-time simulations. Therefore, the license type was changed into GT-Suite-RT, which uses an real-time (RT) solver optimized for speed and is specifically intended for real-time simulations. Simultaneously, RLT Creator was applied to recreate the RLT dependencies for the RT solver as look- up tables and enable the model to work fully stand-alone. In addition, an FRM accelerator tool in the software was used to facilitate settings to further decrease the model's calculation time. The critical change was that the RT solver did not store results that were not required by the user. Therefore, the inputs for the Simulink model needed to be determined through Simulink Harness to pass the results. These actions were performed to smooth the coupling with Simulink, which governed the real-time simulations and later enables the connection with the physical laboratory system (Speedgoat). The conversion of the model to Simulink followed standard GT- Suite guidelines [29]. Note that both the detailed model and the off-line FRM used explicit, continuous state solvers. The RT simulation in Simulink employed a discrete solver, giving direct control over model execution by imposing a fixed time step. The user could therefore adjust the simulation time step to match the application demands and available computational power, while avoiding stability issues imposed by the Courant condition. Solver choice was therefore the most fundamental change compared to off-line FRM discussed earlier and ultimately the one that had greatest effect on accuracy. Matlab R2017a with GT Simulink real-time library was used. GT-Suite provided all the files needed for importing the engine model to Simulink and for converting the model into a real-time application. These files included the GT Simulink interface library, C source code and header files, and object file library files. The GT-Suite-RT model block represented the engine model by using an s-function and the mentioned C code. A Dell Optiplex 760 PC (2008 model year) was used for running real-time simulations. This target PC ran the Simulink Real-Time kernel and had a Modbus card to communicate with the Speedgoat Performance Machine. The Speedgoat target machine runs the engine control system in the laboratory and enables communication between the model and the real engine in complete digital twin implementation. 3. Results 3.1 Component-level validation The results of component-level validation are discussed here briefly. As described earlier, validation was done after calibrating the CAC, turbocharger and in-cylinder submodels separately. Fig. 6 depicts the model´s CAC accuracy by comparing measured and simulated outlet temperatures (T3) for all operating points. Fig. 6. Measured and simulated CAC outlet temperature; errorbars refer to the tolerance of temperature sensor of ±0.5 K The model used the imposed mass flow rate (calculated from imposed inlet conditions) and imposed coolant temperature to calculate the airflow outlet temperature. Intercooler efficiency was the calibrated parameter. According to Fig. 6, the semi-predictive, efficiency-based CAC model was able to capture the T3 within the assumed tolerances for all operating conditions. Note that pressure losses were calibrated directly using the friction loss multiplier and the resulting CAC outlet pressure matched the measured pressure drop almost 1:1 (figure omitted for brevity). The turbine and compressor models were calibrated only at the 75% load point and Fig. 7 shows how well the calibration holds for all other operating points. In the case of both component models, they can be verified as predictive since the trends in their respective outputs were captured correctly across the whole load sweep. A calibration dilemma was unavoidable in both cases. Because the compressor maps (basic model input) were obtained in static conditions, they did not fully replicate the dynamics in a real engine. Fig. 7. Simulated turbocharger results together with measured values and estimated error bars This was clearly seen in the compressor validation results (upper plots in Fig. 7), where it was necessary to prioritize either the accuracy of outlet temperature (T2) or the outlet pressure (p2) calculation accuracy. For the compressor, the tuning targeted p2 as the determinant for the subsequent cylinder submodel accuracy in the combined engine model. The p2 value was predicted to the point, while a constant error (from the calibration at 75% load) in T2 was transferred across the whole load sweep. Nevertheless, the error was close to the sensor accuracy. At this point one should note, that flow temperature measurement in the engine’s air-path systems is generally challenging. The thermocouple´s location, shielding and associated exposure to flow field and heat transfer from the pipe wall affect the sensor’s reading. This is especially important in exhaust flow measurement where the gradients between the wall and flow temperature are particularly large. In this case, nether the device accuracy nor standard deviation gave good estimation of measurement accuracy since it was biased with an in-principle unknown systematic error. Bearing in mind the above consideration, the turbine model was considered valid in terms of outlet temperature (T6) calculation. The trends were predicted correctly. The detail results are shown in the lower plots of Fig. 7. The error bars for T6 measurement were omitted for the reason discussed in the previous paragraph, and the values should be treated as indicative. A deeper insight into the issue of exhaust temperature measurement uncertainty is provided in Appendix A (Fig. A1), depicting exhaust temperature sensor mount sensitivity. The TPA cylinder model was calibrated in order to produce the calculated burn rates for the combustion profile. Afterwards, the cylinder component validation was carried out. Fig. 8 shows the validation plots for the decoupled cylinder model after the combustion profile was applied. The in-cylinder pressures were compared to measurements at all operating points. In addition, Fig. 8 also presents the simulation results of other in-cylinder parameters at the operating point of 75%, highlighting the particularities of the modeling approach and realized combustion strategy. From Fig. 8´s cylinder component validation figures it can be seen that the calculated cylinder pressures remained slightly below the averaged measured pressures (top left plot in Fig. 8). The error bar in the peak pressures illustrates the cylinder-to-cylinder deviation. It appeared that cylinder pressure curves were within this tolerance in other respects but the peak pressure remained lower at high loads. The tuning parameters calibrated for the TPA model were applied to the decoupled cylinder model. These parameters directly affected the peak pressures and flow characteristics on the intake side. Thus, more precise optimization of the tuning parameters could provide a solution to the low peak pressure values. Fig. 8. Measured and simulated cylinder pressures of all four load points with estimated error bars of maximum pressures and simulated cylinder parameters for 75% load point. The x-axis denotes CA deg. Different parts of the cycle are the focus of subsequent subplots. Instantaneous measurements were not available for the other in-cylinder parameters presented in Fig. 8. Those are discussed briefly, however, to highlight the predictive capabilities of the model. The simulated results of the in-cylinder trapped mass, bulk temperature and heat release rate (HRR) followed the expected trends. For instance, the second peak in the in-cylinder temperature (bottom left plot in Fig. 8) that appeared around 230‒340 CA deg, was associated with the exhaust valve lift profile (consult Table 2 for respective valve timings). The valve stopped lifting at 230 CA deg and remained fully open until IVO at 334 CA deg. The outflow became chocked in this period, which together with increasing piston speed, manifested in some in-cylinder charge compression and hence momentarily elevated temperature. This phenomenon was also clearly seen in the in-cylinder trapped mass (top right plot in Fig. 8). Further, observe that the `saddle´ around the peak HRR (bottom right plot in Fig. 8) indicates a superposition of premixed and diffusion flame combustion regimes – characteristic for the state of the art compression ignition engines. The model did not phenomenological differentiate premixed and diffusion combustion. However, as it exactly reproduced the imposed experimental burn rates (scaled by the fuel value for different load conditions), the manifestation of different combustion phases in a real engine was replicated in the simulated heat release. Note, that HRR in Fig. 8 was given as apparent (gross) and normalized by total fuel energy input. 3.2 Validation of the detailed 1D engine model The full model was validated in terms of multiple model parameters. For brevity, the ones selected for discussion were intake manifold pressure (p3), exhaust temperature at turbine inlet (T5), maximum cylinder pressure (cycle and cylinder averaged) and BTE. Fig. 9 shows these parameters´ simulated values and the experimental data, together with their respective uncertainties. This is presented for the engine´s complete load range. The first two parameters, intake manifold pressure and exhaust temperature at the turbine inlet, were pertinent for model validation for two reasons. First, they accumulate all intermediate errors of the detailed air-path model, and second, they shape the overall model performance through feedback/forward mechanisms. Namely, T5 directly affects turbine power and so shapes p3. The p3 affects volumetric efficiency and hence the indicated and peak pressure for all cylinders. Cycle and cylinder-averaged peak pressure was selected for model validation to reflect the in-cylinder model accuracy, providing a clear and consistent measure to assess experimental data uncertainty (standard deviation from cylinder-to-cylinder). Finally, the mentioned accuracy of in-cylinder pressure reproduction shaped the BMEP of the whole engine model. BMEP already accumulated friction loss submodel uncertainty. However, in order to capture the model performance holistically this was further post-processed to BTE, adding the effect of simulated fuel value to the equation. Note that indicated mean effective pressure (IMEP)/BMEP and exhaust temperature are relevant control variables for future uses of the model (see Discussion and outlook section). Fig. 9. Measured and simulated 1D engine model validation parameters of all four load points together with estimated error bars. The complete model´s simulated p3 was below the measured values. This was consistent for all operating points, though at low loads the simulation results were within the estimated measurement accuracy. To this end, a standard deviation from the time-averaged p3 value at given operating point was used as a measure of maximum error caused by the system´s dynamic behaviour (natural cyclic instability of the engine). The inconstancy between simulated and measured results came from an attempt to balance the calibration effort of the model. Namely, the complete engine model was not re-calibrated for p3 and the simulation relied solely on the tuning of the compressor for outlet pressure (p2 – discussed in the previous section for component level validation) and accurate identification of the remaining air-path elements. Considering this, the inconstancy was simply acknowledged at this stage of the research. At this point, one should note rather inclusive tolerance levels of the targeted FRM implementation (Tab. 4). Nevertheless, the validation results for p3 in Fig. 9 indicate that more fine-tuning could still improve performance of the detailed model for high-fidelity simulations. In the complete engine simulation, the underestimation of p3 directly transferred to the error in cylinder model accuracy, illustrated in Fig. 9 (cycle and cylinder-averaged maximum cylinder pressure). Still the inaccuracies were negligible when the model predictions were narrowed to the synthetic control parameters. The simulated BMEP and BTE values in Fig. 9 are within their respective measurement accuracies, even though they incorporated additional friction and fuel value submodel uncertainties. The good accuracy of the model in terms of BMEP and BTE was complemented by good fit of other control parameters like CA50 and the T5. The CA50 results are not shown explicitly but it is suffice to say that the model predictions were well within both the cylinder-to-cylinder and cycle-by-cycle variations in the experimental measurements. Further insight on the actual values is provided in the real-time model simulation results in subsection 3.4. The T5 simulation and experimental results matched well, bearing in mind the earlier- mentioned issue with exhaust flow temperature measurement (subsection 3.1. and Appendix A). Taking account of the above results, the detailed model was considered sufficiently accurate to serve as the reference for the FRM. 3.3 Validation of the FRM The simulation results of the FRM were validated against the results of the detailed 1D engine model. The same parameters discussed in section 3.2 support the discussion on FRM model validity, but with the addition of volumetric efficiency, as this was an important calibration factor when converting the detailed model to FRM. These parameters are presented in Fig. 10, along with the so-called real-time factor for both models. This was the ratio of simulation time per cycle compared to the time the single engine cycle actually takes. The results in Fig. 10 are discussed assuming the FRM tolerance levels in Table 4. As the FRM conversion was an incremental process, steps 1‒3 refer to the governing simplification steps of the FRM. This illustrates how different model reduction measures influenced the accuracy/simulation speed trade-off. Fig. 10. Simulated engine performance parameters and real-time factors of the key steps of the FRM conversion. Step 1 entailed the first simplification of the exhaust manifold and was expected to show some effects on pressure and temperature behaviour. This is seen in Fig. 10 as reduced simulated p3 and higher simulated T5. Step 1´s simplifications were also causing the underprediction of the simulated volumetric efficiencies. Over the operating points of 75‒25%, the volumetric efficiencies differed from the detailed model simulations by more than the set tolerances. Step 2 included the simplification of all other subvolumes in the model´s intake side. This manifested mostly as reduced p3 and lower T5, making the predicted volumetric efficiency lower than the detailed model´s simulations and beyond the tolerances. The results of step 2 have been intentionally omitted in Fig. 10 to support clarity. Note the significant reduction in simulation time when moving from step 1 to step 2. Step 3 achieved another reduction in simulation time by further simplifying the exhaust manifold subsystem into a single element. In this step, the FRM was finally calibrated with friction multipliers. The FRM had a real-time factor of one to two, depending on the operating point, and so was slightly above the real-time target. The simulations mostly remained inside the target tolerances, with the exception of the T5. In the light of the previously discussed challenges and inaccuracies in exhaust temperature measurement, the slightly larger uncertainty in the turbine inlet temperature estimator was considered acceptable. Volumetric efficiency at the lowest operating point of 25% was on the cusp of acceptable tolerance. However, bear in mind that the calibrations during the conversion process were performed only for the highest load point to avoid case-dependent values at a later stage of real-time simulations. As acknowledged, discretization directly impacts simulation accuracy and simulation speed. Therefore, a study was carried out for the FRM model, where all subsystems were simplified in terms of accuracy (step 2). Both shorter and longer discretization lengths were assessed against the lengths originally selected during the FRM process. Only volumetric efficiency and average maximum cylinder pressure were reviewed, together with the factor of real-time to indicate simulation speed. Fig. 11 depicts the results for the different discretization lengths on the intake side: the exhaust side results are in Fig. 12. Both figures also include the results for the detailed model for reference. The model marked `Default´ is the optimum accuracy/speed trade- off, using discretization lengths of 200 mm and 300 mm for the intake and exhaust side respectively. Fig. 11. Effects of reduced and increased intake discretization lengths (accuracy tolerances marked with dash lines). Coarser discretization in general should reduce calculation time because fewer volumes are to be solved. Finer discretization increases accuracy at the expense of calculation time. Movement in either direction is governed by the point at which the trade-offs become unfavourable. This is evident from Fig. 12. On the coarse side, increasing the exhaust discretization length beyond the optimal trade-off point did not bring substantial calculation time reduction, nor did the FRM´s accuracy deteriorate significantly. However, Fig. 11 shows that the intake was far more sensitive to changes in discretization lengths in respect of volumetric efficiency and hence in-cylinder pressure. Here, increasing discretization length (i.e., coarser discretization) significantly reduced the volumetric efficiency due to inability fully resolve intake pulsations. Importantly, the higher discretization lengths also slowed the model´s simulation time, despite fewer volumes to be resolved. The system of differential equations became stiff and the solver drastically increased the simulation time step to keep the solution stable. Fig. 11 and 12 clearly show that finer discretization (i.e., shorter discretization length) both on the intake and exhaust side was unjustifiable because brought only a negligible increase in accuracy but increased calculation time by over 65% compared with the default optimum. Accuracy did not increase substantially because the FRM was grid-resolved in the optimum set point. In other words, the resolution of the model was well-matched to the simulation problem timescale. Fig. 12. Effects of reduced and increased exhaust discretization lengths (accuracy tolerances marked with dash lines). 3.4 Towards real-time implementation Moving the FRM from the standard GT-Power license towards real-time implementation by converting it to GT-Suite-RT, imposed some changes on the model (refer to the Methods section for details). Fig. 13 summarizes the effects of those changes on relevant model output parameters. The results are partially made redundant by Fig. 10 from the previous subsection. At this point, this redundancy was considered beneficial to support the summarizing aspect of following discussion. Fig. 13. Real-time factors and simulated engine performance parameters of FRM in GT-Power and FRM in GT-Suite-RT after required modifications for Simulink coupling. The primary outcome of implementing the real-time license was a reduction in calculation time of over 70% compared to the standalone FRM. This was primarily attributed to the solver change, with the explicit Euler Real-Time used by GT-Suite-RT being particularly optimized for speed. The solver cut down the calculation time to around 4 s (0.35 real-time factor), largely independently of the operating point. Note, that the conventional explicit solver used by GT exhibited large case-by-case variations in this respect. Interestingly, the acceleration of the simulation did not bring substantial change in model accuracy. This is evident from Fig. 13 and so does not require further comments. However, the evident conclusion was that both the stand-alone desktop FRM and its real-time-executable compilation were well within the targeted accuracy limits. The RT-executable model (desktop simulation) was further embedded in Simulink and tested on a target platform. Note, that the target platform had considerably lower central processing unit (CPU) performance than the desktop PC used for prior simulations (Intel Core Q9400@2.66 GHz CPU/8 GB RAM vs Intel Core i7-8750H@2.2 GHz CPU/16GB RAM respectively). Thus, the Simulink RT simulation had to be run with the imposed integration step size of 1 s to still meet the real-time requirements with limited computational performance. This, of course, was insufficient to capture the detailed dynamics of in-cylinder processes, so the results presented in Fig. 14 should be treated as a functionality check of implementation on the target platform model, rather than an in-depth validation. Fig. 14 shows the sample load transient simulation results under the above assumptions. The load and start of injection (SOI) were provided as inputs to the model: BMEP, CA50 and exhaust temperature before turbine were considered as model outputs (responses) in this test case. Additionally, BMEP recalculated directly from model inputs is shown as reference. Fig. 14. Simulated engine inputs and outputs in Simulink RT when time step of 1 s was used. It is evident from Fig. 14 that the more stable BMEP and CA50 mean values were not affected by the large sample time, other than during load-changing where some aliasing effects manifested. Contrarily, the model failed to reproduce the in-cylinder pressure because the sample time was too large to capture the combustion timescale. This was not shown explicitly but one of its results is the unrealistic (oscillating) behaviour in exhaust temperature. 4. Discussion and outlook Despite the evident shortcomings discussed in the previous subsection, Fig. 14 demonstrates a stable, real-time implementation of the FRM on the target platform. Bearing in mind that the capabilities of cutting- edge rapid prototyping machines match the current desktop specification, it can be hypothesized that the control-oriented, physics-based, digital twin of a marine engine is achievable by means of an FRM developed through reduction of a detailed model. Importantly, the results implied that real-time simulation of a marine engine with detailed handling of in-cylinder pressure is possible. This ultimately takes the foreseen application of this model reduction concept far beyond typical digital twin engine applications based on black-box or mean value engine models, as stated in the introduction, mainly limited to on-board monitoring. The detailed CA- based resolving of combustion, demonstrated here, enables coupling with semi-predictive emission models [30, 31] to be captured with the same real-time FRM implementation. According to the offline test in Fig. 13, there would be a sufficient real-time surplus available for this purpose. This opens up a route towards fully integrated model-based control approaches currently not openly demonstrated on a real target platform [32]. Integrated combustion/air-path control is currently the cutting edge and considered an enabler for advanced combustion concepts like RCCI to be developed beyond laboratory-stage tests [7]. As the above developments are still to be followed-up with the baseline created in this research, the developed FRM already can be used as a plant model to design and calibrate an engine performance controller (closed-loop IMEP or combustion onset control, cylinder balancing etc.). All the mentioned prospects are already planned within the next project of the University of Vaasa Combustion Engine Research Group. The Clean Propulsion Technologies project [33] is currently one of the biggest consolidated research endeavours for development of the next generation powertrains for the marine and off-road sectors. The project re-validates the present modeling toolchain to cover extended hardware functionalities related to the engine’s conversion towards RCCI. It further couples the improved engine model with a newly developed fully predictive (performance and emissions) RCCI combustion model, explores the addition of aftertreatment toolchain coupling, aiming for complete powertrain simulation. The process follows the methodology validated in the present research, aiming to develop both the detailed predictive models based on 1D framework and their reduced-order control surrogates. 5. Conclusions A fully physics-based digital twin of a large-bore marine engine, obtained through systematic reduction of a detailed 1D model has been demonstrated for the first time in this work. The solution was tested on the target platform and thoroughly validated against experimental data. The work yielded the following conclusions: • For state of the art mid-speed engines, it was feasible to create a Fast Running Model (FRM) that was real-time capable while maintaining a high level of predictivity, including crank-angle resolved two-zone combustion and detailed thermal management. • This FRM reached grid-resolved simulation with 47 flow components (4 times fewer elements than the detailed model) and internal discretization length of 100% and 150% of the cylinder bore for intake and exhaust components respectively (4 and 6 times coarser discretization). • At such settings, the FRM was three times faster than real-time, while the accuracy loss was comfortably within the 5% tolerance levels for the governing outputs. • The obtained level of predictivity and real-time surplus was more favourable than in the state of the art automotive applications and enables exploring further coupling with semi-predictive emission submodels. For final hardware in the loop demonstration of this modeling approach and its functional extension towards fully predictive advanced low-temperature combustion, the reader is referred to the Clean Propulsion Technologies project repository [33]. Acknowledgement: This research was funded by Business Finland partly under the project ‘Integrated Energy Solutions to Smart and Green Shipping’ and partly under the project ‘Clean Propulsion Technologies’. Appendix A. Fig. A1. Impact of thermocouple location on measured exhaust temperatures (T5, see Fig. 1). The first measurement was conducted from a single exhaust pipe; the second measurement was taken afterwards when the sensors had been replaced and mounted in both the upper pipe (from cylinders 2 and 3) and the lower pipe (from cylinders 1 and 4) before the turbine. Note that the distance between the locations is no more than 5 cm from each other and within a single flow component when the FRM is considered. References [1] European Commission. 3rd IMO GHG study, https://ec.europa.eu/clima/policies/transport/shipping_en (2014, accessed 26 April 2021). [2] DieselNet. IMO Marine Engine Regulations, https://dieselnet.com/standards/inter/imo.php (2020, accessed 7 July 2021). [3] Gamma Technologies. GT-Suite Overview, https://www.gtisoft.com/gt-suite/gt-suite-overview/ (accessed 1 July 2021). [4] Ricardo PLC Software. Wave, https://software.ricardo.com/products/wave (accessed 1 July 2021). [5] Polilink, Innovation channel. Gasdyn - software, https://www.polilink.polimi.it/en/casi-di-successo/gasdyn-software-2/ (accessed 1 July 2021). [6] AVL. AVL Advanced Simulation Technologies, https://www.avl.com/documents/10138/885977/AVL+Advanced+Simulation+Technologies+Catalog.pdf (2017, accessed 1 July 2021). [7] Mikulski M, Ramesh S and Bekdemir C. Reactivity Controlled Compression Ignition for clean and efficient ship propulsion. Energy 2019; 182: 1173-1192. [8] Mavrelos C and Theotokatos G. Numerical investigation of a premixed combustion large marine two-stroke dual fuel engine for optimising engine settings via parametric runs. Energy Convers. Manag 2018; 160: 48-59. [9] Stoumpos S, Theotokatos G, Boulougouris E, et al. Marine dual fuel engine modelling and parametric investigation of engine settings effect on performance-emissions trade-offs. Ocean Eng. 2018; 157: 376-386. [10] Stoumpos S, Theotokatos G, Mavrelos C, et al. Towards Marine Dual Fuel Engines Digital Twins—Integrated Modelling of Thermodynamic Processes and Control System Functions. J. Mar. Sci. Eng. 2020; 8(3): 200. [11] Perabo F and Zadeh M. K. Modelling of a Shipboard Electric Power Systemfor Hardware-in-the-Loop Testing. In: 2020 IEEE Transportation Electrification Conference & Expo (ITEC), June 2020, pp: 69-74. [12] Eriksson L. and Nielsen L. Mean Value Engine Modeling. In: Modeling and Control of Engines and Drivelines. West Sussex UK: John Wiley & Sons Ltd, 2014, pp. 145–210. [13] Llamas X and Eriksson L. Control-oriented modeling of two-stroke diesel engines with exhaust gas recirculation for marine applications. Proc IMechE, Part M: J Engineering for the Maritime Environment 2019; 233(2): 551-574. [14] Salminen HJ, Ciccateri F, Ijäs T, et al. Experimental Demonstration of a Novel Deflagration-based Pressure Gain Combustion Tehnology. In: AIAA Propulsion and Energy 2020 Forum, Virtual event, 24-28 August 2020, AIAA 2020-3871. [15] Mikulski M, Balakrishnan PR and Hunicz J. Natural gas-diesel reactivity controlled compression ignition with negative valve overlap and in-cylinder fuel reforming. Appl. Energy 2019: 254. [16] Hunicz J, Matijošius J, Rimkus A, et al. Efficient hydrotreated vegetable oil combustion under partially premixed conditions with heavy exhaust gas recirculation. Fuel 2020: 268. [17] Li X, Tian H, Shu G, et al. Potential of carbon dioxide transcritical power cycle waste-heat recovery systems for heavy-duty truck engines. Appl. Energy 2019; 250: 1581-1599. [18] Millo F, Di Lorenzo G, Servetto E, et al. Analysis of the Performance of a Turbocharged S.I. Engine under Transient Operating Conditions by Means of Fast Running Models. SAE Int. J. Engines 2013; 6: 968-978. [19] Wei Y, Uppalapati LR and Vernham B. Use of Predictive Engine and Emission Model for Diesel Engine Model Based Calibration. SAE technical paper 2019-01-2227, 2019. [20] Andric J, Schimmel D, Sediako AD, et al. Development and Calibration of One Dimensional Engine Model for Hardware-In-The-Loop Applications. SAE technical paper 2018-01-0874, 2018. [21] Andric J, Schimmel D and Heide J. Calibration Procedure for Measurement-Based Fast Running Model for Hardware-in-the-Loop Powertrain Systems. SAE technical paper 2020-01-0254, 2020. [22] Wu H and Li M. A Hardware-in-the-Loop (HIL) Bench Test of a GT-Power Fast Running Model for Rapid Control Prototyping (RCP) Verification. SAE technical paper 2016-01-0549, 2016. [23] Courant R, Friedrichs K and Lewy H. On the partial difference equations of mathematical physics. IBM J Res Dev 1967; 11(2): 215–234. [24] Wärtsilä. Wärtsilä 20 Product Guide, https://cdn.wartsila.com/docs/default-source/product-files/engines/ms-engine/product-guide-o-e-w20.pdf?sfvrsn=6 (2020, accessed 14 June 2021). [25] Kline S and Mcclintock F. Describing uncertainties in single-sample experiments. Mech Eng. 1953; 75: 3-8. [26] GT-Suite. Engine Performance Application Manual, 2019. [27] GT-Suite. Flow Theory Manual, 2019. https://dieselnet.com/standards/inter/imo.php https://www.gtisoft.com/gt-suite/gt-suite-overview/ https://software.ricardo.com/products/wave https://www.polilink.polimi.it/en/casi-di-successo/gasdyn-software-2/ https://www.polilink.polimi.it/en/casi-di-successo/gasdyn-software-2/ https://www.avl.com/documents/10138/885977/AVL+Advanced+Simulation+Technologies+Catalog.pdf https://cdn.wartsila.com/docs/default-source/product-files/engines/ms-engine/product-guide-o-e-w20.pdf?sfvrsn=6 https://cdn.wartsila.com/docs/default-source/product-files/engines/ms-engine/product-guide-o-e-w20.pdf?sfvrsn=6 [28] Eriksson L. and Nielsen L. Turbocharging Basics and Models. In: Modeling and Control of Engines and Drivelines. West Sussex UK: John Wiley & Sons Ltd, 2014, pp. 211–262. [29] GT-Suite. Engine Performance Tutorials, 2019. [30] Sharma S, Sun Y and Vernham B. Predictive Semi-Empirical NOx Model for Diesel Engine. IJEPE 2019; 13(5): 365–376. [31] Karaky H, Marty P, Tauzia X, et al. Semi-physical NOx and soot model for CI engines: Study of its calibration procedure and portability. IMechE 2020; 234(14): 3414-3428. [32] Indrajuana A, Bekdemir C, Feru E, et al. Towards Model-Based Control of RCCI-CDF Mode-Switching in Dual Fuel Engines. SAE technical paper 2018-01-0263, 2018. [33] Clean Propulsion Technologies. Clean Propulsion Technologies, https://cleanpropulsion.org/ (2021, accessed 4 July 2021). https://cleanpropulsion.org/