Microsoft Word - 66iervolino.docx CHEMICAL ENGINEERING TRANSACTIONS VOL. 65, 2018 A publication of The Italian Association of Chemical Engineering Online at www.aidic.it/cet Guest Editors: Eliseo Ranzi, Mario Costa Copyright © 2018, AIDIC Servizi S.r.l. ISBN 978-88-95608-62-4; ISSN 2283-9216 CFD Modelling of a Spark Ignition Internal Combustion Engine Fuelled with Syngas for a mCHP System Daniele Piazzulloa Michela Costaa,*, Zvonimir Petranovicb, Milan Vujanovicb, Maurizio La Villettac, Carmine Caputod, Domenico Cirilloc a Istituto Motori – CNR, Viale Marconi, 8, 80125 Naples - Italy b University of Zagreb, Faculty of Mechanical Engineering and Naval Architecture, Ivana Lučića, 5, 10000 Zagreb – Croatia c C.M.D. Costruzioni Motori Diesel S.p.A, Via A. Pacinotti, 2 81020 San Nicola La Strada (CE) - Italy d University of Rome “Tor Vergata”, Via del Politecnico 1, Rome - Italy m.costa@im.cnr.it Micro Combined Heat and Power (mCHP) powered with biomass is nowadays a technology attracting increasing interest to develop a local supply chain to produce, process and valorise the available material in territorial areas as much as possible circumscribed, with a considerable reduction also of the CO2 related to transportation. Application for biomass powered mCHP produces environmental benefits by reducing primary energy consumption and associated greenhouse gas emissions and complies with the need for increased decentralization of energy supply. Of particular relevance is mCHP based on biomass gasification due to the negligible particulate matter release with respect to combustion. The present work describes a 3D CFD model of the spark ignition (SI) internal combustion engine (ICE) fuelled with syngas installed in the mCHP pilot system ECO20 manufactured by the Italian company Costruzioni Motori Diesel S.p.A. (CMD). The considered system is made of a gasifier combined with proper syngas cleaning devices, an ICE and a generator to deliver a maximum electrical and thermal power of 20 kW and 40 kW, respectively. For the proper initialisation of the 3D CFD model, the syngas composition is experimentally characterised using a gas-chromatograph on samples collected under real operation. The calculated pressure cycle is verified by comparison with the one calculated through a properly developed 1D ICE model. Main goals of the performed numerical analysis are to study into detail the combustion process and to assess the engine performance characteristics related to the use of syngas. 1. Introduction The process of thermochemical conversion of biomasses through gasification is known since a long time, but it is currently receiving renewed interest as a valuable alternative to incineration, especially for its lower environmental impact. A unique composition of syngas does not exist but its quality in terms of lower calorific value (LHV) and presence of undesirable compounds as TAR and particulate matter strongly depends upon the treated biomass, the gasification equivalence ratio, as well as on a number of design features of the reactor, such as the locations of air inlets and the volume of the gasification zone. According to Banapurmath et al. (2009), the power reduction of an engine operated with syngas is mainly attributed to the lower LHV of the used fuel. The volume of gas/air mixture that enters the engine cylinders limits the power output because it reduces the volumetric efficiency (Lapuerta et al., 2013). Another important parameter is the energy density of the producer gas/air mixture that mainly depends upon the concentration of the combustible components in the gas. The stoichiometric ratio of the air/producer gas mixture is between 1.0 and 1.3, compared to 17 for methane. The LHV of the producer gas obtained from biomass gasification in a downdraft gasifier using air as an oxidizing agent is about 5000 kJ/Nm 3 with the following average composition of the combustible components: 20% H2, 20% CO and 1% CH4. With such concentrations of combustible substances, the energy density of the producer gas/air mixture is about 2400 kJ/Nm3. This value is lower than the energy density of a natural gas/air mixture that is equal to 3400 kJ/Nm3 (natural gas is assumed to be 13 DOI: 10.3303/CET1865003 Please cite this article as: Piazzullo D., Costa M., Petranovic Z., Vujanovic M. , La Villetta M., Caputo C. , Cirillo D., 2018, Cfd modelling of a spark ignition internal combustion engine fuelled with syngas for a mchp system, Chemical Engineering Transactions, 65, 13-18 DOI: 10.3303/CET1865003 100% CH4). Therefore, the theoretical value of power de-rating when a natural gas engine is switched to operate on producer gas is of about the 30% (Costa et al., 2017). In the scientific literature there are few works analysing the operative characteristics of ICEs under syngas powering. From a CFD perspective, Monteiro et al. (2013) developed a multi-zone thermodynamic model for the simulation of the SI combustion of syngas at elevated rotational speed, while Gamiño et al. (2010) simulated the syngas combustion process in a multi-spark engine using a 2D model. Zhang et al. (2012) also studied syngas/air combustion with a 2D model in an engine under partially premixed conditions. The study proposed in the present work is focused on the development of a 3D CFD model of a SI engine fuelled with syngas deriving from woodchip gasification within the commercial tool AVL FIRETM. The considered ICE is actually mounted in an integrated layout aimed at CHP production, coupling a gasifier, a syngas cleaning system and the studied ICE. The computational domain is built according to measurements made on the combustion chamber. Due to the lack of boundary conditions, the CFD model is initialized according to a previously developed 1D ICE model developed within the GT-Power environment. The 3D ICE model reproduces the real operation during the closed valve period. 2. Methodology The methodology of the conducted research, consisting of experimental work, 1D calculation, and 3D CFD modelling is presented. Firstly, the experimental system configuration is explained and the engine main characteristics are shown. The description of the 1D and 3D CFD simulations setup is then given. 2.1 Experimental research of GM Vortec engine The considered CHP unit, as already mentioned, is an integrated system combining a downdraft gasifier, syngas cleaning devices, a SI ICE and an electrical generator. Syngas pureness is indeed guaranteed through a process of cleaning, cooling and filtering, before secondary conversion in a 3.0L GM Vortec I-4 engine. The ICE is 4-cylinder, liquid-cooled, usually fuelled by natural gas and properly modified to work with syngas. It is equipped with an oil tank-refilling pump and oil filter for continuous operation. Further characteristics are reported in Table 1. Table 1. Vortec 3.0 L engine characteristics. Displacement 3000 cc Compression Ratio 10.5 Bore x Stroke [mm] 101.60 x 91.44 Connecting Rod [mm] 203.2 Intake Valve Open (IVO) 101° Intake Valve Close (IVC) 256° Exhaust Valve Open (EVO) 470° Exhaust Valve Close (EVC) 659° Maximum intake/exhaust valve open 11.25 mm RPM 1500 2.2 3D CFD Grid Generation and Grid Sensitivity Analysis The computational domain of the ICE combustion chamber is generated as a moving numerical grid using the Fame Engine Plus tool available in the AVL FireTM software. As the simulations are performed only during the closed valve period, and due to the lack of data regarding the geometry of the intake and exhaust ducts, only the engine combustion chamber is modelled. For a more reliable initialization of the flow-field inside the combustion chamber, in a future work the model of the charge exchange parts will be considered. Here, the flow field initialization at IVC is performed by exploiting a 1D numerical schematization of the engine layout developed in the GT-Power environment. To study the effects of grid size on the model results, three computational grids with different resolutions are generated, with main features shown in Table 2. A coarse, an intermediate and a refined grid are built with a maximum cell dimension chosen equal to 2 mm, 1.5 mm and 1 mm, respectively. In the spark plug zone, to obtain the plug contour details, the cell dimension is set equal to 1/16 of the maximum cell size. Numerical results of this analysis are obtained considering a homogeneous mixture trapped in the combustion chamber at IVC, composed of a stoichiometric mixture of air and syngas. Syngas is obtained from woodchip gasification at an equivalence ratio equal to 0.3. The composition is considered according to samples collection under real operation and analysis through a gas-chromatograph. 14 The turbulent flow field is reproduced within a Reynolds Averaged Navier Stokes (RANS) approach, basing on the k-ζ-f model (Durbin et al., 1991), while the combustion process is modelled using the general gas phase reaction approach through the implementation of the GRI-Mech 3.0 detailed kinetic mechanism (Berkley University, 2017). GRI-Mech 3.0 is an optimized mechanism designed to model natural gas combustion, including NO formation and re-burn chemistry. It is based on 325 elementary reactions involving 53 species, and its prediction capability has long been evaluated against more than 60 experimental studies that distinguish this mechanism from others proposed in the literature. Representative initial conditions resembling those of real operation are assumed, according to the calculations described in the following paragraph. Table 2. CFD numerical grids characteristics at the bottom and top dead centre Numerical Grid Coarse Intermediate Refined Total cells (BDC) 76839 138796 252035 Total cells (TDC) 38451 68387 125018 The results of the grid sensitivity analysis performed are shown in Figure 1. It is possible to notice how the results between refined and intermediate grids are very similar. On the other hand, the coarse grid gives a lower peak of pressure (of about 4 bar). Therefore, for all simulations, due to the calculated results, and due to the computational time, the intermediate grid was chosen. It should be reported, that for a full closed-valve cycle simulation calculation time was about 32 hours for the refined grid, 14 hours for the intermediate, and 6 hours for the coarse grid on the workstation with 4 Intel Xeon CPUs. 2.3 Engine characterization and 1D model The GT-Power flow model involves the solution of the Eulerian equations, namely the solution of the equations of conservation of mass, momentum and energy in one dimension in the absence of viscosity (ideal case). The whole system is discretized into discrete volumes and these volumes are connected each to each other by boundaries. Scalar variables like pressure, temperature, density, internal energy, enthalpy and species concentrations are assumed to be uniform inside each volume, namely they are a function of the longitudinal coordinate, hence the specific position varying from the intake mouth to the exhaust section. The characterization of the considered engine is performed by manually measuring the length and size of each component and by taking information about the adopted operative conditions, as the valves lift law and the engine speed and syngas composition fuelling the engine. Moreover, the predictive combustion model EngCylCombSITurb is chosen due to the possibility of evaluating the influence of variation in the composition of the producer gas. In particular, the parameters reporting the dependency of the laminar flame speed upon the syngas composition and producer gas/air equivalence ratio are tuned following the approach of Hernandez et al. (2005). The sub-model chosen to describe the wall heat transfer is the classical Woschni model (Woschni, 1967) that in the present simulation is tuned by choosing the default parameters. This model allows an estimation of the wall heat transfer by defining the heat convection coefficient as a function of the cylinder pressure, temperature and characteristic length and velocity: ℎ = 3.26 ∙ . ∙ . ∙ . ∙ . (1) where b is the bore length in m, P is the pressure in kPa, T is the gas temperature in K and w is the average cylinder gas velocity in m/s. On the other hand, the 3D CFD simulation are conducted by considering the hypothesis of Reynolds' analogy between the viscous and thermal boundary layers in the momentum and energy transport equations (Viegas et al., 1985) to evaluate the heat flux at the wall. The temperature value is instead imposed. Woschni has converted the equation of the boundary layer theory, so that it is only a function of pressure and temperature (besides the characteristic length and velocity). The presented 1D model is initialized according to the operative conditions reported in Table 3. In particular, the syngas composition is obtained from the measurement campaign through sampling (La Villetta, 2017). In the present simulation, the producer gas at a temperature of 333.15 K is supposed to homogeneously mix with ambient air at 279.15 K in stoichiometric ratio before fuelling the ICE. The pressure of the syngas produced is slightly below the ambient condition, as it takes into account the pressure drops that occurs between the exit of the gasifier and the intake of the ICE due to the presence of cleaning modules mounted on the real layout. The results obtained from the 1D simulation are used to define the initial condition for the 3D CFD model at IVC. In particular, the composition of the mixture trapped inside the combustion chamber at IVC taken from the 1D simulation also considers the presence of a 6.7% of burnt residues composed by CO2, H2O and N2. 15 The mixture composition of syngas, air and residues at IVC and the initial conditions for the 3D CFD simulations are reported in Table 4. Figure 1: Calculated mean pressure, coarse (black), intermediate (blue), and refined (red) grids. Table 3. Initial conditions considered for the 1D simulation. Raw Syngas P = 0.95 bar; T = 333.15 K Ambient Air P = 1.00 bar; T = 279.15 K RPM 1500 Spark Timing -33° BTDC Wall Temperature Head=400 K; Piston=400K; Cylinder=400 K Table 4. Initial conditions at IVC. Species mass fractions [-] CO 0.075 Pmixture @IVC 1.56 bar CO2 0.097 Tmixture @IVC 419.93 K H2 0.0052 TKE @IVC 33.84 m2/s2 O2 0.099 Turbulent length scale @IVC 3.82 mm H2O 0.004 CH4 0.004 N2 0.7146 3. Results Figure 2 reports the comparison between the pressure, rate of heat release (ROHR) and temperature calculated with the 1D simulation and 3D CFD models. These quantities are averaged over the entire volume. The curves show a well captured agreement between the simulations despite the differences in the combustion process and wall heat-flux estimation between the two different numerical approaches, noticed also during the expansion phase after the end of combustion. In Figure 3, the evolution of the mass fractions of the chemical species inside the cylinder is reported. The general species transport model available in the AVL Fire software is activated in the 3D CFD simulation, solving the mass balance equations relative to the species belonging to the GRI-Mech 3.0 kinetic mechanism. A complete depletion of all the species characterizing the producer gas is noticed, as methane, hydrogen and carbon monoxide are totally burned and converted in combustion product as vapour and carbon dioxide. Finally, Figure 4 reports a comparison between 3D visualizations of the temperature and O2 mass fraction evolution for different crank angles in a plane passing through the spark plug symmetry plane. The evolution of these two quantities is strictly related, and allows a clear evaluation of the flame speed propagation inside the cylinder. Finally, it is possible to perform a rough estimation of the electric efficiency of the engine by considering the ratio between the net power produced by the ICE and the syngas primary energy. In particular, Costa et al. (2017) performed an energetic analysis based on experimental measurements on the whole CMD ECO20 330 340 350 360 370 380 390 CAD (°) 5.0x105 106 1.5x106 2.0x106 2.5x106 3.0x106 3.5x106 4.0x106 P re ss u re ( P a ) Coarse Grid Refined Grid Intermediate Grid 16 system as woodchip is used as feedstock. For a value of return user temperature of 20.50°C, using the syngas primary energy equal to 62.25 kW, the relative electric efficiency of the engine ηe_is equal to 22.5%. a) b) Figure 2: Comparison between a) pressure cycle and ROHR, and b) temperature evolution of the 1D and 3D CFD simulation. Figure 3: Evolution of the chemical species mass fractions inside the cylinder in the 3D CFD simulation. Crank Angle O2 mass fraction Temperature (K) 345° 350° 355° 360° 365° Figure 4: Comparison between 3D evolution of the O2 mass fraction (left) and gas temperature (right) in a plane passing through the symmetry axis of the spark plug. Considering the same syngas primary energy, the net power produced by the engine is evaluated as the indicated power over the entire cycle P deriving from the 1D simulation. Therefore, the indicated efficiency is equal to: = , = .. = 32.28% (2) P re s s u re ( P a ) 300 320 340 360 380 400 420 CAD (°) 400 800 1200 1600 2000 GT - Power 3D CFD 300 320 340 360 380 400 420 CAD (°) 0 0.05 0.1 0.15 0.2 0.25 CH4 CO CO2 H2 H2O O2 17 The engine electric efficiency has to consider the mechanical losses η and the generator efficiency η : _ = ∙ ∙ (3) Assuming the mechanical efficiency equal to 0.75 and the generator efficiency equal to 0.95, the electric efficiency of the engine results equal to 0.23. This value is in line with the work of Costa et al. (2017) and the literature data concerning cogeneration systems powered with syngas. 4. Conclusion A 3D CFD numerical model of a spark ignition engine fuelled with syngas deriving from woodchip gasification is presented. The considered engine is actually mounted in an integrated layout aimed at the combined heat and power production, coupling a gasifier, a syngas cleaning system, the spark ignition engine and a generator. The numerical grid is modelled thanks to measurements taken to define the geometrical characteristics of the combustion chamber, and the 3D CFD model is initialized thanks to a proper 1D model developed in the GT-Power environment. The presented numerical tool reproduce the characteristics of the real engine operation during the closed valve period, in good agreement with the 1D simulation, but also gives a greater detail of the chemical species concentration rate of change during the combustion process. It is noticed that under syngas fuelling the electric efficiency of the engine is equal to 23%, thus revealing a depowering with respect to the powering mode by natural gas. Reference Banapurmath N.R., Tewari P.G., 2009, Comparative performance studies of a 4-stroke CI engine operated on dual fuel mode with producer gas and Honge oil and its methyl ester (HOME) with and without carburetor, Renewable Energy,; 24:1009-15. Berkley University, 2017, , accessed 14.01.2017. CMD S.p.A., 2017, , accessed 18.01.2017. Costa M., La Villetta M., Massarotti N. et al, 2017, Numerical analysis of a compression ignition engine powered in the dual-fuel mode with syngas and biodiesel, Energy,137:969-979. Costa M., La Villetta M., Cirillo D. et al, 2017, Performance analysis of a small-scale combined heat and power system powered by woodchips, Proceedings of the 11th Conference on Sustainable Development of Energy, Water, and Environment Systems, SDEWES. Durbin P.A., 1991, Near-wall turbulence closure modelling without damping functions. Journal of Theoretical and Computational Fluid Dynamics, 3:1-13. Gamiño B., Aguillón J., 2010, Numerical simulation of syngas combustion with a multi-spark ignition system in a diesel engine adapted to work at the Otto cycle, Fuel, 89(3):581-591. General Motors, 2017, , accessed 14.01.2017. Hernández J.J., Tinaut F.V., Horrillo A., 2001, Thermochemical behaviour of producer gas from gasification of lignocellulosic biomass in SI engines, SAE Technical Paper No. 2001-01-3586. Hernandez, J. J., Lapuerta, M., Serrano, C. et al, 2005, Estimation of the laminar flame speed of producer gas from biomass gasification, Energy & fuels, 19(5), 2172-2178. Lapuerta M., Monteiro E., et al., 2013, Performance of a syngas-fuelled engine using a multi-zone model. 1st International Congress on Bioenergy, 23-25. La Villetta M., 2017, Biomass-to-Energy systems: techno-economic aspects and modelling approaches for characterization and improvement of performance, Ph.D. Thesis, University of Naples “Parthenope”, Naples, Italy. Monteiro E., Rouboa A., Sotton J., Bellenoue, M., 2013, Performance of a syngas-fuelled engine using a multi- zone model. 1st International Congress on Bioenergy, 23-25. Thermoflex software, 2017, , accessed 14.01.2017. Viegas J.R., M. W. Rubesin, C. C. Horstman, 1985, On the Use of Wall Functions as Boundary Conditions for Two-Dimensional Separated Compressible Flows. Technical Report AIAA-85-0180, AIAA 23rd Aerospace Sciences Meeting, Reno, Nevada. Woschni G., 1967, A Universally Applicable Equation for the Instantaneous Heat transfer Coefficient in the Internal Combustion Engine, SAE Technical Paper 670931. Zhang F., Yu R., Bai X. S., 2012, Detailed numerical simulation of syngas combustion under partially premixed combustion engine conditions, International Journal of hydrogen energy, 37(22):17285-17293. 18