Study of the process of combustion in a high swirl engine using computational fluid dynamics

Resumen

La simulación numérica de motores de encendido por chispa presenta incertidumbres que comprometen la fiabilidad de los resultados debido a la interrelación existente de la cinética química y la cinemática del flujo dentro de la cámara de combustión. En este trabajo se evalúa el efecto que tiene el uso de distintos mecanismos de reacción aplicados a un modelo de combustión tipo ecuación G y su interacción con el campo de velocidades dentro de la cámara de combustión de un motor CFR que posee características técnicas especiales (alto swirl y bajas rpm) que permiten estudiar dicho fenómeno. La investigación se realizó numéricamente empleando mecánica de fluidos computacional (CFD), utilizando la ecuación G(con y sin cinética de reacción) como modelo de combustión y el modelo de turbulencia RANS -k − ϵ- RNG. Se compara y valida la simulación contra una prueba experimental que reporta la curva del diagrama indicador. Se estudió la influencia que tienen los mecanismos de reacción y el swirl en el desarrollo del frente de llama y sus variables asociadas dentro de la cámara de combustión. Finalmente se encontraron ligeras discrepancias entre los mecanismos de reacción empleados que podrían ser importantes en simulaciones de motores encendidos por chispa distintos del aquí estudiado. Así mismo se apreció como un aumento en la intensidad del swirl desvía la llama notablemente de los planteamientos teóricos. Las técnicas y códigos actuales requieren ajustar distintos parámetros y por tanto se hace necesaria la validación contra datos experimentales antes de sacar conclusiones de los casos numéricos.

Abstract

The numerical simulation of spark ignition engines has uncertainties that compromise the reliability of the results due to the interrelationship between chemical dynamics and flow kinematics within the combustion chamber. This work evaluated the effect of using different reaction mechanisms applied to a G-equation combustion model and its interaction with the velocity field inside the combustion chamber of a CFR engine, which has special technical characteristics (high swirl and low rpm) that allowed the study of this phenomenon. The research was carried out numerically by computational fluid dynamics (CFD) techniques: using a G-equation combustion model (with and without reaction kinetics) and a RANS -k − ϵ- RNG turbulence model. The numerical simulation was compared and validated against reported experimental pressure indicator diagram. The influence of reaction mechanisms and swirl in the development of the flame front and associated variables within the combustion chamber was studied. Finally slight discrepancies were found between reaction mechanisms used in this study that could be important in simulations of spark ignited engines that are different from the one studied. It was noted that an increase in the intensity swirl yields a flame development deviation from the theoretical models. The current codes and simulation techniques require tuning of various parameters, and thus, the validation against experimental data is needed before drawing conclusions from the numerical cases.

Palabras clave

Motores de Combustión Interna;Simulación Numérica;Ecuación G

Keywords

Internal Combustion Engines;Numerical Simulation;G-Equation

1. Introducción

La combustión es considerada la mayor fuente generadora de energía del planeta; representando más del 90% del consumo global, entre transporte, calefacción y centrales termoeléctricas [1]. El estudio de la combustión es importante debido al requerimiento de producir emisiones mucho más bajas manteniendo una alta eficiencia térmica. Entre los sistemas de combustión existentes, los motores a combustión interna encendidos por chispa representan una gran proporción y ligeras mejoras en su eficiencia dan a lugar a grandes impactos en términos ambientales y económicos [2]. Con el fin de lograr esto, se ha buscado entender los fenómenos físicos involucrados en la combustión a través de diversos experimentos y simulaciones numéricas. Siendo esta última la herramienta de investigación dominante en épocas recientes. No obstante, el acople de las ecuaciones de transporte a las de cinética química junto a la modelación del régimen del flujo turbulento que son requeridos para realizar las simulaciones numéricas presentan obstáculos que comprometen la fiabilidad de los resultados [3].

La investigación actual de motores encendidos por chispa se centra en reducción de las emisiones y aumento de la eficiencia de la combustión del motor [2]. Estas investigaciones cobijan campos experimentales y numéricos que buscan extender el entendimiento y modelación del fenómeno y sus variables asocidas como: acople turbulento-cinética química, ignición de mezcla aire-combustible, mecanismos cinéticos de reacción, efectos de evaporación de combustible, entre otros.

El efecto que tiene swirl en la combustión de los motores de encendido por chispa ha sido estudiado desde de hace más 50 años [4], encontrándose que la condición dinámica que tiene la mezcla al momento de entrar a la cámara de combustión tiene una incidencia preponderante en el desarrollo y propagación de la llama. Sin embargo, a pesar de los esfuerzos experimentales muchos elementos están aún sin resolver [5]; [6] ;  [7]. Con el advenimiento del CFD se ha deseado continuar indagando en este fenómeno, sin embargo, el acople entre las ecuaciones de momento, el modelo de turbulencia y cinética química presentan incertidumbres que no han sido abordadas en estudios previos. La ampliamente validada ecuación G [8] como modelo de combustión permite evaluar las implicaciones que tiene el uso de distintos modelos cinéticos y su interacción con la turbulencia inducida por el swirl dentro de la cámara de combustión del motor y así determinar sí estas interacciones afectan los resultados numéricos predichos.

Este trabajo se enfoca en las consecuencias del uso de distintos mecanismos de reacción y su interacción con el campo de velocidades de la mezcla a ser encendida. La metodología empleada involucra una simulación tridimensional completa que es validada contra una prueba experimental realizada en la cámara de combustión de un motor CFR (Cooperative Fuel Research). Contrario a los motores modernos [9] la morfología de este motor no busca promover tumble y es el swirl el principal agente que modifica el campo de velocidad inicial de entrada de la mezcla. Esto hace al motor interesante desde el punto de vista fenomenológico de la combustión ya que sus características ténicas: baja velocidad, alto swirl y turbulencia permiten hacer simplificaciones de cinética química que ayudan a indagar acerca de los efectos de estos dos factores en el desarrollo y propagación de la llama dentro de la cámara de combustión. Además la ubicación longitudinal de la bujía [10] provoca un cambio en la forma de propagación de la llama lo cual es atractivo para este estudio.

Estas características fueron ya aprovechadas en estudios anteriores referentes a encendido de mezclas pobres [11], específicamente, en el uso de bujías de plasma (Plasma Jet Ignition) para el estudio del uso de mezclas de hidrógeno y metano como combustible para vehículos comerciales [12] y en la modelación numérica de motores HCCI [13].

La investigación se desarrolló en dos partes fundamentales. Primero, se aproximaron los valores de la variable presión de las simulaciónes numéricas por medio de una búsqueda sistemática a la curva de presión experimental reportada. Segundo, se analizó el campo de datos tridimensional extraido de las simulaciones ya validadas. Se compararon variables de combustión (especies, número de Damkölher) entre dos modelos de ecuación G utilizados.

Las pruebas experimentales emplearon gasolina y gasolina-etanol como combustible en un CFR reportando la curva de presión, la fracción de masa quemada, y la energía liberada [14]. Este motor no es transparente por lo cual no fue posible realizar mediciones de campos de velocidad dentro de la cámara de combustión.

Tras la validación de la curva de presión dentro de un rango de error aceptable, se realizó un segundo estudio modificando el valor del swirl -el componente cinématico/turbulento- con el fin de apreciar las consecuencias entre este cambio y los modelos de combustión, encontrándose diferencias en las variables globales de presión, temperatura media y energía liberada.

Este trabajo modeló numéricamente la combustión en la cámara del CFR empleando modelos del tipo ecuación G, que ha sido reportada como la metodología mas adecuada para este tipo de simulaciones [2]. Se usó un modelo de ecuación G más cinética en equilibrio y otro empleando cinética química. Además, se analizó el efecto del swirl con respecto al avance y forma del frente de llama contra los modelos teóricos [15].

Los resultados permiten apreciar ligeras discrepancias entre los modelos de combustión empleados y su interacción con el campo de velocidades 3D presente en la cámara de combustión. Dichas diferencias no podrían ser apreciadas empleando los modelos cero o unidimensionales tradicionales para motores de combustión interna [16].

2. Descripción del estudio experimental

El estudio experimental realizado por Mantilla [14] requirió el uso de una tarjeta de control marca MoTeC 3.41G2 Engine Control Unit (ECU) para manejar los eventos de inyección de combustible y avance al encendido. La medición de presión en la cámara de combustión se realizó con un transductor Kistler 6125B, a 1/10 de ángulo del cigüeñal. Los flujos de aire y combustible fueron medidos. La temperatura de escape se obtuvo empleando una termocupla tipo Omega K con +/-0.01C de exactitud. La figura 1 ilustra el esquema del montaje experimental.

Figura 1. Montaje Experimental [14].

Los gases de escape fueron medidos empleando una probeta de acero inoxidable ubicada entre el múltiple de escape y el muffler, seguido a esto se encontraba una trampa de vapor del tipo membrana. El gas en el escape se midió empleando los siguientes monitores de emisión continua:

• Rosemount Analytical Model 755R (± 0.1%)
• Rosemount Analytical Model 880 (± 0.05%)
• Rosemount Analytical Model 880 (± 0.01%)
• Rosemount Analytical Model 955 (± 1 % Full Scale)
• Heated FID J.U.M. Engineering (± 10ppm)
• Beckman Industrial Model 400A Flame Ionization Detector (± 1PPM)

2.1. Datos experimentales

El estudio experimental recolectó las condiciones termodinámicas de la mezcla al momento de ser introducida en la cámara de combustión. La tabla 1 registra los valores de presión atmosférica y temperatura en distintas partes del cilindro y la camisa, así como la relación de equivalencia y la fracción de masa residual. Las dimensiones del motor así como los ángulos de apertura y cierre de las válvulas de escape son resumidas en la tabla 2. Algunos de los datos fueron extrapolados y escalados de los planos del manual del fabricante [10]. El eje coordenado Z está sobrepuesto al eje de simetría del motor, con origen a la altura de la culata del motor.

Tabla 1. Condiciones de entrada experimento [14]

Parámetros	Magnitud
Presión	99 KPa
Temperatura de gases de entrada (Puerto)	303 K
Relación de equivalencia ϕ gases de entrada	1
Fracción de masa residual	0.13
Temperatura de la cara del pistón	353 K
Temperatura de la camisa	373 K
Temperatura culata	353 K

Tabla 2. Dimensiones de la cámara y el cilindro [14]

Parámetros	Magnitud
Diámetro del Pistón	0.0889 m
Longitud de la biela	0.2000 m
Radio del cigüeñal	0.05715 m
Relación de compresión	7
Ángulo de apertura válvula de escape	500°
Ángulo de apertura válvula de admisión	710°
Ángulo de cierre válvula de escape	15°
Ángulo de cierre válvula de admisión	210°
Ángulo de Chispa	345°
Duración de la chispa	8°
Distancia Horizontal a la bujía	0.0439 m
Distancia Vertical a la bujía	-0.00925 m
Revoluciones del Motor	900 rpm

3. Descripción de los casos de estudio numéricos

Los casos de estudio ilustrados en este artículo fueron simulados durante ciclo cerrado, comprendido entre el cierre de la válvula de admisión y la apertura de la válvula de escape.

3.1. Elección del software

Entre las diversas opciones de software comercial existente para simulación de motores se opta por el paquete comercial CFD CONVERGE™.

CFD CONVERGE™ es un código computacional que permite la simulación tridimensional de flujo compresible e incompresible, flujo reactivo y transitorio en geometrías complejas con superficies móviles e inmóviles. Las ecuaciones gobernantes son discretizadas empleando un esquema de volúmenes finitos y los modelos de turbulencia disponibles emplean la metodología Reynolds Averaged Navier Stokes (RANS) ó Large Eddy Simulation(LES).

CFD CONVERGE™ genera la malla para el dominio del caso, siempre que el usuario ingrese la geometría y esta pueda ser triangulada empleando el formato STL [17]. Este método permite que la malla pueda ser modificada (comprimida, refinada, adaptada) durante la simulación. La malla adaptativa de CFD CONVERGE™ usa el denominado campo de subgrilla que se define como la diferencia entre el campo resuelto y el campo actual. Esta nomenclatura es proveniente de la modelación de los esfuerzos de sub-grilla de LES [18]. La idea principal es emplear series infinitas para modelar los campos de sub-grilla que no son capturados debido al filtro LES, sin embargo, al ser imposible tomar todos los términos se aproxima al primer término de la serie que está en función del campo resuelto. Esta técnica es adaptada a la generación de malla adaptativa y nuevas celdas son creadas acorde a una tolerancia entre la variable de interés (velocidad, especies, o temperatura) y el campo resuelto.

CFD CONVERGE™ incluye además un solucionador de cinética de reacción denominado SAGE™que admite mecanismos de reacción en formato CHEMKIN™. Este solucionador resuelve para cada paso tiempo la cinética química en cada celda del dominio. El solucionador toma los valores promedio de la celda, considerándola un reactor homogéneo abierto.

3.2. Modelos Termodinámicos

Los modelos termodinámicos empleados hacen uso de la teoría de gas ideal, junto a la ponderación de los coeficientes JANAF [19] para diversas propiedades. Para los casos simulados las propiedades de la mezcla aire-combustible fueron ingresadas al programa. El combustible (Gasolina) se modela como isooctano iC₈H₁₈. Los modelos de las propiedades de la mezcla se incluyen en el Apéndice A.

3.3. Modelo de Turbulencia

El modelo de turbulencia empleado fue del tipo RANS -k − ϵ- RNG [20]. Este modelo RANS ha sido reportado como exitoso por diversos autores [2] ;  [21] y es el más empleado en la simulación estándar de motores a combustión interna, debido que captura el efecto del swirl en la turbulencia y mejora la predicción de la rata de deformación. Variables críticas que mejoran la exactitud de la modelación numérica.

3.4. Modelo de combustión

Se seleccionan dos modelos:

• Empleando la ecuación G [22] y equilibrio:

${\frac {\partial {\tilde {G}}}{\partial t}}+{\tilde {\overset {\rightarrow }{v_{f}}}}\cdot \triangledown {\tilde {G}}=$${\frac {\overline {{\rho }_{u}}}{\overline {\rho }}}s_{t}^{0}\mid \triangledown {\tilde {G}}\mid -$$D_{t}{\tilde {k_{M}}}\mid \triangledown {\tilde {G}}\mid$

(1)

donde ${\textstyle v^{'}}$ es la velocidad turbulenta, a₄ = 0.78, b₁ = 2, b₃ = 1.0, y l es la escala integral del flujo.

Esta metodología propuesta por Peters [22] modela el Corrugated Flametes Regime. Entre sus ventajas está el evitar los inconvenientes referentes al counter-gradient diffusion presentes en los modelos del tipo variable de progreso. Acorde con Poinsot [23] esta forma de la ecuación G propuesta por Peters guarda mucha similitud con los modelos Flame Density Function o modelos Sigma. La figura 15 y el gráfico 9-18 presente en [15] illustran que este motor opera en la zona donde u′ > s_l lo que implica acorde con la bibliografía [23] que el régimen de operación se encuentra en la zona corrugated validando así el uso de esta metodología en este estudio.

• Empleando el solucionador de cinética química -SAGE- incluido en el software [17] Se configuraron los siguientes casos:
• Modelo cinético reducido de 56 especies y 168 reacciones. Desarrollado por Liu, Ming, Zhao y Pang. [24].
• Mecanismo cinético de 140 especies y 663 reacciones. Desarrollado por Yoo, Luo, Lu, Kim, y Chen [25].

3.4.1. Modelo de ignición

La ubicación angular de la bujía es insensible en su posición siempre que se encuentre a la distancia determinada por las condiciones geométricas -ver tabla 2. Se optó por no simular la geometría de la bujía, aunque se reconoce que tiene un efecto sobre el campo de velocidades [26]. Para modelar las distintas etapas de la bujía, (Glow Discharge, Plasma Arc, entre otros[15]) se ubica una fuente esférica de diámetro igual a 1 mm. Se considera que la duración del término fuente es equivalente a 1.4 milisegundos considerando experimentos previos [27] lo que es equivalente a 8 grados de giro del cigüeñal.

3.4.2. Velocidad laminar de llama

El modelo de velocidad laminar de llama elegido fue el desarrollado por Gülder [28], con las constantes para isooctano.

3.4.3. Velocidad turbulenta de Llama

Determinar la velocidad turbulenta de llama sigue siendo un problema abierto, que continua bajo investigación [22] ;  [23]. Tradicionalmente se utilizan correlaciones entre variables del fluido. Estos modelos tienen constantes que deben ajustarse a las condiciones de turbulencia presentes en la cámara de combustión. Debido a la baja relación de compresión y alto grado de swirl la constante b1 de la ecuación (3), correspondiente a la correlación de llama turbulenta de Peters [22] fue modificada considerablemente siendo 31 su valor final (el valor teórico es igual 2). El desarrollador del software indica que esta constante debe ser ajustada de acuerdo con la discretización del dominio [17].

3.5. Escalas de tiempo y fenómeno transitorio

Al ser el ciclo de un motor un fenómeno transitorio y con diversas escalas de tiempo (chispa, ignición, propagación de llama) se debe emplear un paso de tiempo adaptativo. La duración total de la simulación se encuentra dada en grados de giro del cigüeñal. La escala temporal más pequeña se tomó como 1e-8 s y la más grade como 1e-5 s.

3.6. Configuración de malla

La malla es autogenerada por el software acorde con las condiciones geométricas del modelo a simular. Dicha malla es estructurada y el software internamente trunca las celdas (paralelepípedos) de manera que se ajusten a la geometría del dominio. Las características de malla tales como skewness, aspect ratio y non-orthogonality cumplen con las condiciones de una malla estructura perfecta con contadas excepciones en la frontera. Esto puede ser verificado empleando un posprocesador como Paraview™ y su filtro mesh quality[29]. Respecto a su relación con la discretización espacial, para la iteración numérica se se ingresó al programa un valor de deferred correction[21] igual a 0.5 empleando esquemas de primer y segundo orden que implica que todo el dominio se discretizó empleando CDS. Esta simulación se llevó a cabo considerando la cámara de combustión completa, es decir, no se modelaron sectores circulares, debido a la posición de la bujía que impide el uso de la simetría axial. La figura 2 ilustra el dominio 3d cilíndrico empleado y un corte transversal del mismo.

Figura 2. Malla en −150° BTC.

Las cuatro mallas base al inicio de la simulación usadas en este estudio estaban compuestas de celdas que tenian 4, 3, 2 y 1 mm de arista que corresponden a 14000,33000,100000 y 800000 elementos respectivamente. La malla elegida para realizar el reporte (ver sección 4.3) fue la de 2 mm, debido a que su número de elementos se encuentra entre los valores típicos encontrados en la literatura de simulación de motores [30]; [31]; [32] ;  [33].

3.6.1. Adaptación de malla

Debido al alto swirl presente en el caso de estudio sumado a los altos gradientes de temperatura y velocidad a través de la zona de llama, se debe realizar adaptación de malla local. Para el caso de temperatura y velocidad esta refinación fue de tres y dos niveles respectivamente. Estos valores fueron asignados directamente al software y están relacionados en forma de potencias de 2 [17]. La opción de refinamiento local estuvo activa para la duración total de los casos de estudio. Además, para el caso de la fuente de ignición se realizó una refinación esférica de 4 y 5 niveles. Sin embargo, esta última solo estuvo presente desde −15° ATC hasta el TDC.

3.7. Condiciones iniciales y de contorno

Las condiciones iniciales de las simulaciones fueron calculadas siguiendo los siguientes supuestos:

• Los productos residuales en la cámara de combustión, son especies en equilibrio equivalentes a una combustión completa.
• Se hace un ponderado másico entre el porcentaje de masa residual y la relación estequiometrica aire-combustible (ver tabla 3).

Tabla 3. Especies Iniciales en Fracción Másica

Especies	Fracción másica
IC8H18	0.05414
O2	0.18960
N2	0.71980
CO2	0.02498
H2O	0.01148

• Otros estudios empleando motores CFR han reportado que el rango de swirl está entre (2 y 8) según sea el levantamiento de la válvula de admisión que se tenga [34]; [35] ;  [36]. Acorde con el ensayo experimental y levantamiento de válvula dado [14], se considera que el motor está trabajando con un swirl igual a 8 (ver tabla 4) y su perfil es equivalente a la constante estándar de Amsden [37].

Tabla 4. Parámetros generales de la simulación

Parámetro	Valor
Tiempo de inicio simulación	-150 BTC
Tiempo de fin simulación	140 ATC
Swirl	8.0
Perfil de Swirl (Constante de Amsden [37])	3.11
Número de Schmidt Global	0.78
Número de Prandtl Global	0.9
Solucionador de Navier-Stokes	PISO COMPRESIBLE
Energía Cinética turbulenta	98 m2/s2
Disipación turbulenta	2E06 m2/s3
CFL Convectivo	⩽ 1
Especies Involucradas	IC8H18 O2 N2 CO2 H2O CO H2 N2O NO NO2
Constantes velocidad Turbulenta de llama	a4=0.78 b1=31/b1=14 b3=1.0
Valor Término Fuente(Ignición)	8e+11 m−3s−1

• La simulación numérica inicia, justo al cierre de la válvula de admisión.
• La disipación turbulenta y la energía cinética turbulenta fueron ajustadas de forma iterativa. Mantilla [14] indica que estas condiciones parecen ser insensibles en los rangos empleados en este caso, y que las diversas correlaciones existentes para deducirlas dan origen a resultados muy parecidos.

Estas condiciones iniciales fueron tomadas considerando estudios previos 0D y 1D hechos por Mantilla [14].

Condiciones de Contorno

Las condiciones de contorno de temperatura (ver tabla 1) y momento se fijaron considerando la ley de la pared y la ley logarítmica siguiendo el esquema propuesto por Launder [38] referente al uso de una u otra ley según sea el valor de y⁺. El uso de estas condiciones de contorno se justifica en que estos perfiles y⁺ no varían en la misma medida que el flujo dentro del dominio durante la simulación.

3.8. Solucionador de las ecuaciones de Navier-Stokes

Las ecuaciones de momento, continuidad y energía junto con las de especies químicas son resueltas empleado el solucionador PISO [21] para flujo compresible. La tabla 4 lista los parámetros generales más importantes empleados en las simulaciones.

Los residuales empleados en la solución numérica se listan a continuación:

• Momento, energía, especies y escalares pasivos: 1e-4.
• Presión: 1e-8.
• TKE y ϵ: 1e-3.

4. Resultados

A continuación se presentan los resultados de ambas simulaciones, en el primer caso se aborda la ecuación G y posteriormente la simulación empleando los mecanismos de isooctano reducido.

4.1. Ecuación G

Se realizaron 30 experimentos numéricos en los cuales se buscaba replicar la curva de presión indicada de la forma más exacta posible. Esto fue llevado a cabo de forma iterativa hasta alcanzar un nivel de error aceptable. Tras diversos análisis de sensibilidad, los parámetros modificados en este caso fueron:

• La constante b1 (ver tabla 4)
• El valor del término fuente (ver tabla 4)

El primer parámetro modifica la velocidad turbulenta de llama y por tanto la propagación de la misma dentro de la cámara, el segundo se relaciona con la ignición de la mezcla.

4.2. Ecuación G + Modelo cinético

En cuanto a la simulación empleando el solucionador de cinética química incluido en el software SAGE y el modelo de cinética reducido de iso-octano propuesto por Liu et al. [24]. Se encontraron diversas dificultades que debieron ser superadas de forma iterativa.

La primera dificultad se debe al mecanismo de reacción empleado. Este no fue planteado para todo el rango de presión descrito por la simulación. El problema principal era la baja presión al momento de la ignición (debido a la baja relación de comprensión). Esta situación impide al mecanismo detectar o iniciar el encendido de la mezcla homogénea. Una prueba en un reactor a volumen constante (CVR) hecha en CHEMKIN™ reveló este comportamiento. Debido a que al utilizar un mecanismo completo para el isooctano como el propuesto por Curran [39] la demanda de recursos crece de forma exponencial. Se optó por emplear una metodología híbrida que consiste en lo siguiente:

• Tomar el experimento con el modelo anterior (ya validado). Emplear la ecuación G y el modelo de ignición que garantiza el encendido de la mezcla.
• Una vez encendida la mezcla y alcanzadas las presiones de diseño del mecanismo reducido, se emplea el solucionador cinético. El software permitió realizar este algoritmo de forma sencilla, configurando la opción SAGE dentro y fuera del frente de llama representado por la ecuación G [17].

Los parámetros de la simulación son los mismos de la tabla 4. El solucionador cinético se usó desde -15 BTC hasta +40 ATC. Las especies, reacciones del mecanismo y su base de datos termodinámica pueden encontrarse en el artículo de Liu et. al [24] en formato CHEMKIN™. Adicionalmente, la constante b1 de la ecuación de Peters (3) fue ajustada iterativamente a un valor de 14. Para esto, se requirió realizar 10 simulaciones con los parámetros ya mencionados.

4.3. Curva de presión

Los resultados apropiados se encontraron al realizar la simulación con una malla típica de 100000 al inicio de la simulación, utilizando los parámetros descritos en la tabla 4. Para validar dichos resultados y evaluar su independencia de malla se realizaron pruebas con mallas de 14000, 33000 y 800000 elementos utilizando los mismos paramétros. Se encontró que al hacerse la malla más burda la curva de presión se corría a la derecha de la curva experimental producto de un avance de llama más lento y al hacerse más fina se corría hacia la izquierda producto de un avance de llama más veloz. Esto es consecuencia directa de la velocidad turbulenta de llama y la constante b₁ presente en las ecuaciones 2 y 3 que dependen del tamaño de la grilla y por tanto al usar la misma constante para distintas mallas se presentan discrepancias. Esta técnica implica que la constante b₁ debe ajustarse caso a caso y al usar o no mecanismos cinéticos.

La curva de presión obtenida se reporta (fig. 3) junto a la de los resultados experimentales.

Figura 3. Curvas de presión de este trabajo.

El análisis de esta figura indica lo siguiente:

• Para la ecuación G:
• Error en la ubicación angular del pico de presión es menor al uno por ciento (1%).
• Error en la magnitud del pico de presión es igual al ocho por ciento (8%).
• Para la ecuación G + SAGE empleando el mecanismo de Liu et al.:
• Error en la ubicación angular del pico de presión es menor al cero uno por ciento (0.1%).
• Error en la magnitud del pico de presión es igual al seis punto cero cinco por ciento (6.05%)

Se evidencian dos fuentes de error para la discrepancia en la magnitud del pico de presión:

• La constante b1 (ver tabla 4). Afecta la velocidad turbulenta de llama, [17] ajusta una mayor velocidad de llama a un valor mayor de dicha constante.
• El hecho de usar isooctano, como combustible sustituto (surrogate) de la gasolina.

Se evidencia que el uso del mecanismo cinético mejora ligeramente la predicción de la magnitud pico presión y de la ubicación del mismo, esto debido al mayor número de especies consideradas que capturan mejor la liberación de energía. No obstante que ambas metodologías arrojen resultados tan cercanos, ilustran que debido a las bajas RPM del motor, la combustión se encuentra cercana al equilibrio químico. De la misma forma se aprecia como el modelo cinético reducido propuesto por Liu [24] (Figura 4 ;  Figura 5 y 7) modela satisfactoriamente la combustión en el motor, a un menor costo computacional al momento de ser comparado con mecanismos más robustos como el de Yoo [25].

Figura 4. a)Efecto del swirl ecuación G b) Efecto del swirl ecuación G + SAGE.

Figura 5. Curvas de temperatura.

4.4. Efecto del swirl en la combustión

La figura 3 fue obtenida de forma iterativa, procurando el mejor acercamiento a la curva experimental; variando tanto el swirl como la constante b1. Un estudio adicional buscó comparar el efecto del swirl en la forma de la curva de presión manteniendo constantes los demás parámetros.

El swirl es un parámetro de entrada, relacionado con el movimiento del aire dentro del cilindro cuando éste pasa a través de los puertos de admisión. La figura 4 ilustra para ambas metodologías que para un mayor swirl (y por ende más turbulencia) aumenta la velocidad a la cual avanza la llama y corre el punto de presión máxima a la izquierda de la curva experimental. Por el contrario, un menor swirl retarda la combustión y corre el punto de presión máxima a la derecha de la curva experimental. Existe entonces una fuerte dependencia y sensibilidad entre el valor del swirl y la ubicación del pico de presión. La curva correspondiente a la reportada en la figura 3 es la de swirl igual a 8, que es la que mejor se ajusta a la curva experimental.

4.5. Temperatura promedio de la cámara de combustión

En la figura 5, se aprecia como los resultados presentan una buena correlación, en cuanto al pico de temperatura y su magnitud. Adicionalmente se hace una comparación entre las temperaturas promedio y el nivel de swirl -ver figura 6. A mayor swirl el pico de temperatura se corre hacia la izquierda, en congruencia con el pico de presión de la figura 4a.

Figura 6. Curvas de temperatura vs swirl.

El análisis de la figura 5 indica lo siguiente:

• Para la ecuación G:
• Error en la magnitud del pico de temperatura es igual al cuatro por ciento (4%).
• Para la ecuación G + SAGE empleando el mecanismo de Liu et al.:
• Error en la magnitud del pico de temperatura es igual al tres por ciento (3%).

4.6. Energía liberada

La energía liberada calculada se comparó con su contraparte experimental. La figura 7 ilustra la liberación de energía para ambos casos (numérico y experimental).

Figura 7. Liberación de calor entre simulaciones.

La discrepancia en el calor liberado puede ser explicada con ayuda de la figura 6. Sí se considera la diferencia de picos de temperatura, la masa dentro del cilindro y que el calor específico de la mezcla es cercano a 1, un cálculo rápido arroja un resultado que se ajusta al rango de la diferencia en la liberación de energía entre las simulaciones y la parte experimental que se aprecia la figura 7. Además acorde con la figura 21 del artículo original de validación del mecanismo de Liu et. al. [24], los autores hacen una clara referencia que esta sobrepredicción se debe a incertidumbres referentes a las condiciones iniciales, la simplificación de la geometría interna y al uso de isooctano como surrogate para simular gasolina.

Curvas como la fracción de masa quemada (X_b) permiten caracterizar la combustión dentro de motores de ignición por chispa [15]. Heywood [15] propone dos formas de calcular la curva de masa quemada, una basada en resultados experimentales (6) y otra basada en la función de Wiebe(7).

$X_{b}={\frac {P^{1/n}V-P_{0}^{1/n}V_{0}}{P_{f}^{1/n}V_{f}-P_{0}^{1/n}V_{0}}}$

(6)

Δθ representa la duración de la combustión y θ₀ el momento en el cual salta la chispa. Los valores de: a, y m hacen parte de regresiones ya establecidas [15]. El valor de n, se toma como 1.26, siendo este un valor de ajuste empírico [15]. La duración de la combustión se considera al punto donde la relación PV^1.15 es máxima [40]. Esta metodología se aplicó tanto a los resultados númericos como a los experimentales. La figura 8 ilustra la fracción de masa proveniente de los datos experimentales y la calculada empleando la fórmula (6).

Figura 8. Curva de fracción de masa quemada.

Existe una ligera discrepancia al inicio de la curva y entre las pendientes de los casos de estudio contra la parte experimental. Esto se encuentra íntimamente relacionado a la ignición y propagación del frente de llama y que sugiere que existe una mayor turbulencia en los resultados experimentales con respecto a las condiciones consideradas en los estudio numéricos. Que la fracción de masa quemada sea una cantidad normalizada permite aseverar lo anterior, a pesar de la influencia de la discrepancia en la liberación de energía expuesta previamente.

4.7. Especies

Para evaluar el desempeño de las especies, se toma un punto del diagrama indicador donde la curva de liberación de energía ha alcanzado un valor constante. Estos valores no se pueden comparar con datos experimentales (ya que existen aún reacciones mientras el gas se mueve hacia el tubo de escape donde se midieron las emisiones). Sin embargo permiten comparar las diferencias en términos de especies químicas de ambos modelos empleados. Cabe recordar que los modelos reducidos son de 56 especies [24] y 140 especies [25] respectivamente. El punto tomado se encuentra en 85° ATC. Las fracciones másicas de especies comunes se ilustran en la figura 9.

Figura 9. Fracciones Másicas Comunes.

Los resultados muestran que ambas metodologías predicen de forma similar las especies estables (CO2, H2O), justificando la suposición de equilibrio debida a la baja velocidad del motor. Sin embargo, la diferencia de las fracciones de las especies intermedias (CO,NO) entre modelos es considerable. El CO y NO poseen una discrepancia de más de casi el doble al comparar el modelo de Liu contra la ecuación G. Esto puede ser explicado con las reacciones adicionales que involucran CO, NO presentes en el mecanismo cinético con respecto a la cinética en equilibrio. Adicionalmente, la fracción de O2 en el caso de Liu es inferior en dos ordenes de magnitud a los otros modelos, sustentando que la cinética química en este caso, juega un papel preponderante y explica el aumento de especies intermedias cuando se usa el mecanismo de Liu. Además, el comportamiento de estas especies afecta según Liu et. al. [24] la predicción de la energía liberada. El modelo de Yoo et. al. [25], carece de las especies N2O, NO y NO2 siendo esto la fuente de error en la predición del O2 para este caso.

4.8. Desarrollo de la combustión dentro de la cámara

Hecha la validación con las curvas anteriores. Se procede a analizar la simulación tridimensional (3D). Se examina el perfil y el desarrollo de la llama dentro de la cámara de combustión obtenidos desde la simulación. Las figuras fueron realizadas haciendo un corte por el plano medio del término fuente, su ubicación se aprecia en la tabla 2.

La figura 10a ilustra la fracción de isooctano en un plano medio del cilindro empleando el modelo de la ecuación G. La gráfica ilustran las zonas correspondientes a: la zona de precalentamiento, la zona de reacción y la zona de mezcla quemada. Además, se aprecia el crecimiento de la zona quemada, tras la ignición; cuya forma es asimétrica, globular y bastante alejada de planteamientos teóricos que postulan el frente de llama como esférico [15].

Figura 10. Consumo isooctano Ecuación G a) −10° BTC, Fracción CO en TDC b) ecuación G c) ecuación G +SAGE Liu et al..

4.8.1. Efectos del modelo de combustión

Ambas metodologías tienen un buen acierto -con ligeras discrepancias- tanto en la curva de presión como en la curva de fracción de masa quemada. Ahora, se indaga sobre los posibles cambios en las variables dentro de la cámara que conlleva el uso del modelo de la ecuación G contra el modelo de la ecuación G empleando el SAGE. En primera instancia las figuras 10b y 10c indican como el uso de la metodología SAGE; cinética química posterior y anterior a la llama provoca la generación de zonas de alta concentración de la especie CO con respecto a la otra metodología que emplea equilibrio y explicando la discrepancia acumulada de la figura 9b.

El número de Damköhler relaciona el tiempo característico turbulento con el tiempo característico químico. Las figuras 11 y 12 muestran el número de Damköhler para ambas metodologías. La diferencia entre los valores pico de dicho número de Damköhler es menor al 3%, adicionalmente ambas gráficas guardan una alta similitud geométrica. La discrepancia en el modelo híbrido parece provenir del uso del mecanismo cinético, cuyas reacciones se apagan debido al agotamiento de los reactivos aumentando así el tiempo químico. Es notable apreciar que a medida que el frente de llama avanza, la zona aún sin quemar posee un número de Damköhler mucho más alto.

Figura 11. Comparación Número de Damköhler en −4° BTC a) ecuación G b) ecuación G +SAGE Liu et al..

Figura 12. Comparación Número de Damköhler en TDC a) ecuación G b) ecuación G +SAGE Liu et al..

El análisis que se puede apreciar en las figuras 13 y 14 conduce a la conclusión en la cual la zona delante del frente de llama es comprimida y presenta mayor turbulencia aumentando la magnitud de la vorticidad e incrementando localmente el número de Damköhler. Las lineas de corriente superpuestas a los contornos de vorticidad soportan esta afirmación.

Figura 13. Comparación Vorticidad en −4° BTC a) ecuación G b) ecuación G +SAGE Liu et al..

Figura 14. Comparación Vorticidad en TDC a) ecuación G b) ecuación G +SAGE Liu et al..

Finalmente, se grafican los resultados para toda la cámara de Reynolds turbulento Re_t vs número de Damköhler Da. Se nota que este motor se encuentra en la zona de operación (recuadro azul) para motores de chispa [15] (ver fig. 15). Este diagrama claramente ilustra a que a medida que avanza el frente de llama dentro de la cámara de combustión existe una reducción, tanto del número de Damköhler como del número de Reynolds turbulento y ambos casos (Ecuación G y Ecuación G + SAGE) reportan tendencias similares en cuanto a este comportamiento. Dicha reducción se soporta en las figuras Figura 11 ;  Figura 12, 13 y 14. En el primer caso ya que a medida que la llama avanza el valor global de Da se va a reduciendo al igual que el de la vorticidad, disipándose la turbulencia y reduciendo el valor de Re_t.

Figura 15. Diagrama Comparativo Ret vs Da, adaptado de Turns [1].

4.8.2. Efecto del swirl dentro de la cámara de combustión

El uso de un modelo 3D permite examinar los vectores de velocidad en la zona de llama (G=0) y apreciar que esta asimetría, se debe principalmente al intenso swirl y turbulencia dentro de la cámara de combustión que deforma el frente de llama (ver fig. 16). Estas interacciones entre el campo de velocidades y la cinética química y sus consecuencias no pueden ser adquiridas empleando modelos cero o unidimensionales.

Figura 16. Efecto del swirl = 8 en TDC ecuación G. (a) Avance de llama. (b) Energía Cinética Turbulenta.

Se analiza el efecto en el campo de velocidades del uso de un mayor o menor swirl al momento de introducir la mezcla dentro de la cámara. Desde la curva de presión se dedujo que un mayor swirl implicaba una mayor turbulencia y una combustión más rápida. Las figuras 17 y 18 soportan esta deducción e ilustran la estrecha relación entre el avance del frente de llama y el campo de velocidades inicial, siendo este efecto advectivo del swirl el encargado de deformar la llama e inducir una mayor energía cinética turbulenta (tke) que promueven un avance de llama mucho más veloz. Tomando el área del círculo del corte de la cámara de combustión como referencia, la llama ha recorrido aproximadamente el 68%,80% y 88% del área del corte para un swirl de 6,8 y 10 respectivamente.

Figura 17. Efecto del swirl = 6 en TDC ecuación G. (a) Avance de llama. (b) Energía Cinética Turbulenta.

Figura 18. Efecto del swirl = 10 en TDC ecuación G. (a) Avance de llama. (b) Energía Cinética Turbulenta.

4.8.3. Evolución de la energa cinética turbulenta

El nivel de turbulencia es máximo al momento de la carrera de admisión y decrece hasta llegar cerca del punto muerto superior donde tiene un ligero incremento para continuar descendiendo [41]. Este comportamiento es bien conocido [42] y el correspondiente a toda la cámara se presenta en la figura 19 para el caso de la ecuación G. Se aprecia el descenso al inicio y el pico de turbulencia correspondiente a la propagación de la llama.

Figura 19. Evolución de la energa cinética turbulenta.

5. Conclusiones

En este estudio se evaluaron dos metodologías de modelos de combustión empleando la ecuación G y su interacción con el campo de velocidades dentro de la cámara de combustión de un motor CFR durante el ciclo cerrado. La comparación de los resultados de ambos casos indican que aunque ambas metodologías representan con una tolerancia aceptable la curva de presión (variable global) existen discrepancias en cuanto al comportamiento de las especies y el campo flujo. Estas ligeras discrepancias podrían volverse notables en otros casos de simulación de motores a gasolina donde simplificaciones como la de equilibrio no fuesen válidas y por tanto la interacción del campo de velocidades y la combustión fuese mucho más crítico y dependiente.

Además se nota el importante efecto que tiene el swirl inicial en el desarrollo de la llama y su relación con el modelo de combustión. Se evidencia como a medida que el swirl aumenta la deformación de la llama es mucho más pronunciada y alejada del modelo esférico teórico cuando se considera al swirl como el factor preponderante en dicha deformación. Sin embargo, el modelo de combustión influye ligeramente en la curva de presión y es evidente desde la gráfica que para para la metodología SAGE, un aumento substancial en el swirl no representa un incremento considerable en la curva de presión.

A la a luz de estos resultados más investigaciones son necesarias. El efecto entre la turbulencia y la combustión dentro de la cámara de combustión del motor y su interacción con su campo de velocidades aún no puede ser claramente establecido, más cuando el modelo de simulación considerado el más eficiente, debe ajustarse en sus constantes radicalmente para los cambios en el modelo de combustión empleado.

Agradecimientos

Este estudio fue posible gracias a la empresa Convergent Science que fue la proveedora del software y las licencias académicas que permitieron obtener los resultados aquí expuestos.

Anexo A. Propiedades de la mezcla gaseosa

El archivo anexo PropiedadesGas.txt incluye las propiedades de la mezcla de gases empleada en el estudio. Las tres columnas están organizadas acorde con la tabla A.1.

Tabla A.1. Propiedades del Gas

Temperatura	Viscosidad	Conductividad
[K]	[N * S/m2]	[W/(m * k)]

References

1. [1] S. Turns; An Introduction to combustion; McGrawHill (2011)
2. [2] R. Reitz; Directions in internal combustion engine research; Combustion and Flame, 160 (2013), pp. 1–8
3. [3] D. Haworth; A review of turbulent combustion modelling for multidimensional in-cylinder cfd; Combustion and Flow Diagnostics, and Fundamental Advances in Thermal Fluid Sciences (2005)
4. [4] J.H. Johnson, Effect of swirl on flame propagation in a spark ignition engine, Tech. rep., SAE Technical Paper (1962).
5. [5] P. Hill, D. Zhang; The effects of swirl and tumble on combustion in spark-ignition engines; Progress in energy and combustion science, 20 (5) (1994), pp. 373–429
6. [6] K. Lee, C. Bae, K. Kang; The effects of tumble and swirl flows on flame propagation in a four-valve si engine; Applied thermal engineering, 27 (11) (2007), pp. 2122–2130
7. [7] I. Khan, G. Greeves, C. Wang, Factors affecting smoke and gaseous emissions from direct injection engines and a method of calculation, Tech. rep., SAE Technical Paper (1973).
8. [8] L. Liang, R.D. Reitz, Spark ignition engine combustion modeling using a level set method with detailed chemistry, Tech. rep., SAE Technical Paper (2006).
9. [9] W.W. Pulkrabek; Engineering fundamentals of the internal combustion engine; Prentice Hall (2004)
10. [10] Waukesha, The cfr engine: A standard of excellence for a half of century (June 1980).
11. [11] L. Tozzi, E. Dabora; Plasma jet ignition in a lean-burn cfr engine; T. combustion institute (Ed.), Nineteenth Symposium on Combustion (1982), pp. 1467–1474
12. [12] G. Karim, I. W. I, Y. Al-Alousi, Methane-hydrogen mixtures fuels, Hydrogen Energy (1996) 625-631.
13. [13] D. Flowers, S. Aceves, R. Smith, J. Torres, J. Girard, R. Dibble; Hcci in a cfr engine: Experiments and detailed kinectic modeling; SAE (Ed.), World Congress (2000)
14. [14] J. Mantilla; Modelado de la combustion en mezclas gasolina-etanol en motores de combustion interna; Universidad Nacional de Colombia, sede Medellin (2010) Ph.D. thesis
15. [15] J. Heywood; Internal Combustion Engine Fundamentals; McGrawHill (1989)
16. [16] E. Pariotis, G. Kosmadakis, C. Rakopoulus; Comparative analysis of three simulation mmodel applied on a motored internal combustion engine; Energy conversion and management, 60 (2012), pp. 45–55
17. [17] C. Science; Converge CFD User Manual; Convergent Science (2014)
18. [18] E. Pomraning; Development of large eddy simulation turbulence models; University of Wisconsin-Madison (2000) Ph.D. thesis
19. [19] Janaf thermochemical tables, Tech. rep., National Bureau of Standards (1971).
20. [20] V. Yakhot, S. Orzag, S. Thangam, T. Gatski, C. Speziale; Development of turbulence model s for shear flows by a double expansion technique; Physics of Fluids, 4 (1992), pp. 1510–1520
21. [21] H. Versteeg, W. Malalasekera; An Introduction to Computational Fluid Dynamics: The Finite Volume Method; Pearson (2007)
22. [22] N. Peters; Turbulent Combustion; Cambridge University Press (2000)
23. [23] T. Poinsot, D. Veynante; Theoretical and Numerical Combustion; Edwards (2005)
24. [24] Y. Liu, M. Jia, M. Xie, B. Pang; Enhancement on a skeletal kinectic model for primary reference fuel oxidation by using a semi-decoupling methodology; Energy & Fuels (2012), pp. 7069–7083
25. [25] C. Yoo, Z. Luo, T. Lu, H. Kim, J. Chen; A dns study of ignition characteristics of a lean iso-octane/air mixture; Proceedings of the Combustion Institute, 34 (2) (2013), pp. 2985–2993
26. [26] Z. Tan, R. Reitz; An ignition and combustion model based on the level-set method for sparking ignition engine multidimensional modeling; Combustion and Flame, 145 (2006), pp. 1–15
27. [27] R. Maly, M. Vogel; Initiation and propagation of flame fronts in lean ch 4-air mixtures by the three modes of the ignition spark; Symposium (International) on Combustion, Vol. 17, Elsevier (1979), pp. 821–831
28. [28] O. Gulder; Sae paper 841000; SAE (Ed.), West Coast Internacional Meeting & Exposition (1984)
29. [29] U. Ayachit, The paraview guide: A parallel visualization application.
30. [30] F. Perini, D. Sahoo, P.C. Miles, R.D. Reitz; Modeling the ignitability of a pilot injection for a diesel primary reference fuel: impact of injection pressure, ambient temperature and injected mass; SAE International Journal of Fuels and Lubricants, 7 (2014-01-1258) (2014), pp. 48–64
31. [31] F. Perini, E. Galligani, G. Cantore, R.D. Reitz, Validation of a sparse analytical jacobian chemistry solver for heavy-duty diesel engine simulations with comprehensive reaction mechanisms, Tech. rep., SAE Technical Paper (2012).
32. [32] C.F.F. Rodriguez, J. Mantilla, Modeling hcci engine combustion coupling cantera to kiva 4, Tech. rep., SAE Technical Paper (2015).
33. [33] F. Perini, G. Cantore, R.D. Reitz, An analysis on time scale separation for engine simulations with detailed chemistry, Tech. rep., SAE Technical Paper (2011).
34. [34] Y. Zhang, E. Kung, D. Haworth; A pdf method for multidimensional modeling of hcci engine combustion: effects of turbulence/chemistry interactions on ignition timing and emissions; Proceedings of the combustion institute, 30 (2) (2005), pp. 2763–2771
35. [35] W. Moore, M. Foster, K. Hoyer, Engine efficiency improvements enabled by ethanol fuel blends in a gdi vva flex fuel engine, Tech. rep., SAE Technical Paper (2011).
36. [36] S.J. Fygueroa, J.O. Araque, C.G. Villamar; Perfiles de velocidad en el cilindro de un motor alternativo; Información tecnológica, 18 (1) (2007), pp. 73–78
37. [37] B.T. Amsden A.A., O’ Rourle P.P, Kiva ii: A computer program for chemically reactive flows with sprays, Tech. rep., Los Alamos National Laboratory (1989).
38. [38] B.E. Launder, D. Spalding; The numerical computation of turbulent flows; Computer methods in applied mechanics and engineering, 3 (2) (1974), pp. 269–289
39. [39] H. Curran, M. Mehl, W. Pitz, C. Westbrook; Modeling of component mixtures revelant to gasoline; European Combustion Meeting (2009)
40. [40] R. Brown, Combustion data acquistion and analysis, Tech. rep., Loughborough University Department of Aeronautical and Automotive Engineering (2001).
41. [41] G. Davis, R. Tabaczynski, R. Belaire, The effect of intake valve lift on turbulence intensity and burnrate in si engines-model versus experiment, Tech. rep., SAE Technical Paper (1984).
42. [42] G. Davis, A. Mikulec, J. Kent, R. Tabaczynski, Modeling the effect of swirl on turbulence intensity and burn rate in si engines and comparison with experiment, Tech. rep., SAE Technical Paper (1986).