Application of the level set method to model the premixed combustion process in an otto two stroke engine

Resumen

En este artículo se ha utilizado el método numérico conocido como level set para modelar el proceso de combustión en un motor Otto de dos tiempos. Se ha implementado el level set en un software libre de mecánica de fluidos computacional (CFD) basado en volúmenes finitos, el cual ha proporcionado los campos de presión y temperatura, así como la propagación del frente de llama. Con el fin de validar el modelo, los resultados obtenidos numéricamente se han comparado con datos experimentales, obteniendo una concordancia bastante satisfactoria. Asimismo, se ha comparado el método de level set con otro procedimiento numérico muy utilizado, mostrando la diferencia entre ambos resultados.

Abstract

In this paper, the numerical method level set has been used to model the combustion process in an Otto two-stroke engine. The level set has been implemented in a CFD (Computational Fluid Dynamics) software based in finite volumes. The pressure and temperature fields have been obtained, such as the propagation of the flame front. In order to validate this model, the numerically obtained results have been compared with experimental data, verifying a satisfactory concordance between both of them. Besides, the level set method has been compared with other numerical procedure, showing the difference between both results.

Palabras clave

Mecánica de Fluidos Computacional ; Ecuación G ; Combustión ; Motores

Keywords

Computational Fluid Dynamics ; G equation ; Combustion ; Engines

1. Introducción

Hay un gran número de procesos en la naturaleza que envuelven varias fases, o varias sustancias separadas por una interfaz. Ejemplos son las olas del mar, una gota de aceite flotando en agua, la fusión de un material, etc. Una técnica numérica para tratar el seguimiento de estas interfaces móviles es el método de level set. Este fue introducido por los matemáticos Osher y Sethian [1] en 1988, y la principal ventaja es la eficacia en cuanto al seguimiento de formas que cambian bruscamente su geometría y cuya velocidad de interfaz depende de la curvatura o la dirección normal. Desde su implementación, el procedimiento de level set ha sido utilizado para diversas aplicaciones que envuelven interfaces móviles, tales como crecimiento de burbujas en ebullición [2]  and [3] , grietas en materiales [4] , movimiento de olas debido al viento [5] , solidificación [6]  and [7] , etc. En cuanto a combustión, el método level set se utiliza para caracterizar el movimiento del frente de llama que separa los productos quemados de los no quemados. Ejemplos de trabajos de combustión utilizando el método de level-set son los de [8]  and [9] , y en cuanto a combustión en motores, se pueden encontrar en la literatura los trabajos de [10] , [11] , [12] , [13] , [14]  and [15] .

En este trabajo se ha utilizado el método level set para modelar el proceso de combustión en un motor Otto de dos tiempos. El proceso de combustión influye notablemente en el rendimiento del motor y la emisión de contaminantes. Además, el calor generado afecta a la vida de los componentes del motor, sobre todo el pistón y aros. El aporte más relevante del presente artículo es haber hecho una comparación entre el método de variable de progreso de reacción, implementado en un gran número de softwares comerciales, y el level set. A partir de dicho análisis se han citado las ventajas del level set empleando resultados experimentales. Asimismo, se ha discutido el efecto del proceso de reinicialización (que forma parte del método de level set) en los resultados. La estructura del artículo es la siguiente: en primer lugar, se resumen los conceptos básicos del funcionamiento del motor estudiado y del proceso de combustión. A continuación se explica la ecuación de level set empleada para caracterizar el frente de llama, así como su implementación en un código de CFD y, finalmente, se exponen los resultados y las conclusiones.

2. Especificaciones técnicas y funcionamiento del motor

El motor estudiado en el presente trabajo ha sido analizado en un artículo previo [16] , en el cual se simuló numéricamente el proceso de barrido mediante un software comercial de CFD. Las especificaciones técnicas se resumen en la tabla 1 .

Tabla 1. Especificaciones técnicas

Parámetro	Valor
Tipo de motor	Dos tiempos, Otto
Potencia (kW)	7,5
Velocidad (rpm)	6.500
Cilindrada (cm3 )	127,3
Diámetro (mm)	53,8
Carrera (mm)	56
Relación de compresión	9,86:1
Sistema de inyección	Inyección directa
Instante de ignición	20° antes del PMS*

• . PMS: punto muerto superior.

Aunque el objetivo primordial de este trabajo es la descripción del método numérico, es conveniente explicar primero el funcionamiento del motor, que es prácticamente similar en todos los motores de dos tiempos de inyección directa. Al comienzo del ciclo, tras producirse la combustión, el pistón desciende desde el punto muerto superior (PMS), fig. 1 (a) (nótese que las flechas verticales de la figura indican la dirección del pistón). A medida que el pistón desciende, la lumbrera de escape es descubierta (y consecuentemente abierta) y con ello los gases quemados empiezan a ser expulsados, fig. 1 (b). El pistón continúa descendiendo. Cuando el pistón descubre (y consecuentemente abre) las lumbreras de transferencia, el aire a presión entra desde el cárter y ayuda a expulsar el resto de gases de escape al exterior, fig. 1 (c). Después de alcanzar el punto muerto inferior, el pistón asciende. Las lumbreras de transferencia y luego las de escape se cierran, fig. 1 (d), y el aire es comprimido a medida que el pistón asciende. En esta etapa es cuando tiene lugar la inyección de combustible, el cual se mezcla homogéneamente con el aire. Un poco antes de que el pistón alcance el PMS, salta la chispa en la bujía y tiene lugar la combustión, comenzando el ciclo de nuevo.

Figura 1. Funcionamiento del motor.

Este motor es de encendido por chispa (por bujía). El proceso mediante el cual se produce la ignición se explica a continuación. La inyección de combustible comienza tras el cierre de la válvula de escape (en este momento el cilindro está lleno de aire, y el pistón se encuentra ascendiendo, figura 1 (d)). Como este motor es de inyección directa, el combustible entra directamente en la cámara de combustión. Tal y como se indica en la figura 2 , la cámara de combustión es hemisférica y se encuentra situada en la parte superior del cilindro. El inyector se encuentra centrado y la bujía se encuentra en un lateral de la cámara. El combustible se introduce en forma de spray y se mezcla homogéneamente y en proporción estequiométrica con el aire. Tras saltar la chispa en la bujía, se forma un frente de llama que se propaga a toda la carga.

Figura 2. Inicio de la combustión en un motor de encendido por chispa.

En el presente trabajo se ha utilizado un sensor piezoeléctrico acoplado a la bujía para medir experimentalmente la presión en el interior del cilindro a lo largo de cada ciclo de funcionamiento, obteniendo el diagrama representado en la figura 3 .

Figura 3. Presión en el interior del cilindro medida experimentalmente.

3. Método de level set

Una vez explicados los fundamentos del funcionamiento del motor y el proceso de combustión, se procederá a describir el método numérico que se ha utilizado para modelar el avance del frente de llama.

3.1. Ecuación de level set aplicada a combustión

Para caracterizar el movimiento del frente de llama que divide el campo fluido en una región quemada y una región no quemada se ha utilizado lo que se conoce como ecuación de level set. El frente de llama se caracteriza mediante el valor cero de la función level set, la cual en el presente trabajo se ha definido como la distancia al frente de llama, con su signo, y le ha asignado como G . En la figura 4 se indican los isocontornos de la función level set, donde cada línea representa una determinada distancia al frente de llama. En este caso, se ha elegido G > 0 para la región de gases quemados y G < 0 para la mezcla no quemada de aire-combustible, aunque también podría haberse elegido lo contrario.

Figura 4. Frente de llama e isocontornos de level set alrededor del mismo.

Supóngase que x_f(t) es el trayecto de un punto en el frente de llama, es decir, donde la distancia es cero:

$G(x_{f}(t),t)=0$

( 1)

Derivando respecto al tiempo resulta la siguiente expresión:

${\frac {\partial G}{\partial t}}+\nabla G\cdot {\frac {d{\overset {\rightarrow }{x}}_{f}}{dt}}=$$0$

( 2)

Esta última ecuación fue la introducida por Osher y Sethian [1] cuando desarrollaron el método de level set, aunque en el método inicial no se refería a combustión sino a flujo bifásico en general. En combustión, es muy común denominar la ecuación anterior como ecuación G .

A modo de ejemplo, se muestra la utilización de la ecuación (2) para modelar el crecimiento de una esfera. Particularmente, se ha tomado una esfera de radio inicial 1 m que crece radialmente a una velocidad de 1 m/s. Se ha elegido una malla uniforme constituida por hexahedros. Con estos parámetros, a medida que transcurre el tiempo la superficie mantiene siempre la forma esférica, pero con el radio incrementándose progresivamente. El resultado de aplicar la ecuación (2) se muestra en las figura 5 (a) y (b) para los instantes de tiempo 1 s y 5 s respectivamente.

Figura 5. Crecimiento de una esfera de radio inicial 1 m y velocidad de crecimiento radial 1 m/s. (a) 1 s; (b) 5 s.

La ecuación (2) «mueve» el cero de la función level set, G , exactamente a la velocidad del frente de llama de una manera muy precisa. El problema es que, después de cada iteración, G deja de ser la distancia al frente de llama. Para resolver esto, Sussmann et al. [17] introdujeron un procedimiento iterativo para reinicializar G , el cual se basa en la siguiente ecuación:

${\frac {\partial G}{\partial \tau }}=sign(G_{0})\left(1-\right.$$\left.\left|\nabla G\right|\right)$

( 3)

siendo τ un tiempo artificial que se introduce porque la ecuación (3) es una ecuación temporal pero que debe ser resuelta en cada paso de tiempo del problema genérico, G₀ el campo de la función level set al principio de cada iteración en este tiempo artificial introducido y sign(G₀ ) el signo de G₀ , dado por:

$sign(g_{0})=\left\{{\begin{array}{ccc}-1&si&G_{0}<0\\0&si&G_{0}=0\\+1&si&G_{0}<0\end{array}}\right\}$

( 4)

Como el signo vale cero cuando la función level set es cero, G tiene el mismo cero que G₀ . Este proceso de reinicialización tiene la importante peculiaridad de que la función level set comienza a ser reinicializada desde el frente de llama. Para ver esto más claramente, la ec. (3) se puede reescribir como:

${\frac {\partial G}{\partial \tau }}+{\overset {\rightarrow }{w}}\cdot \nabla G=$$sign(G_{0})$

( 5)

siendo:

${\overset {\rightarrow }{w}}=sign(G_{0}){\frac {\nabla G}{\left|\nabla G\right|}}$

( 6)

El vector ${\textstyle {\overset {\rightarrow }{w}}}$ tiene módulo unitario, y su dirección es normal al frente de llama y apuntando hacia uno u otro lado dependiendo del signo, por lo cual la reinicialización se propaga hacia el exterior del frente de llama con velocidad unitaria. Esto significa que G se reinicializa cerca del frente de llama primeramente, y se va propagando a ambos lados del mismo. A modo de ejemplo, la figura 6 muestra el procedimiento de reinicialización para una esfera de radio unitario. Se ha representado la función level set en un plano que pasa por el centro de la esfera. El valor donde la level set es nula se indica mediante una línea más ancha. La ecuación de reinicialización es tan eficaz que funciona aunque las condiciones iniciales para la función level set estén lejanas de la función distancia. Concretamente, para este caso se ha comenzado con G =+1 en la región externa a la esfera y G =−1 en la región interna a la esfera (ver instante τ = 0 en la fig. 6 ). Si lo que se pretende es reconstruir la función level set a lo largo de una franja nΔτ a ambos lados del frente de llama, se necesita resolver la ecuación (3) para un tiempo artificial τ=0,...,nΔτ. Para este caso, se ha considerado un intervalo de tiempo artificial Δτ=0,01. Cuando se resuelve la ecuación (3) para 10 intervalos de tiempo artificial, es decir, τ=10Δτ, la level set se actualiza una distancia de 10 × 0,01 = 0,1 a cada lado (ver instante τ=10Δτ en la figura 6 ), por tanto un grosor total de 0,2. Para un tiempo de τ=20Δτ, la level set se actualiza a cada lado una distancia de 20 × 0,01 = 0,2 (ver instante τ=20Δτ en la figura 6 ) y por último para un tiempo de τ=50Δτ, la level set se actualiza una distancia de 50 x 0,01 = 0,5 (ver instante τ=50Δτ en la figura 6 ) a cada lado del frente de llama.

Figura 6. Procedimiento de reinicialización de la level set para una esfera.

Este proceso de reinicialización tiene la ventaja de que, tanto el vector unitario normal, ecuación (7), como la curvatura, ecuación (8), se determinan con gran precisión. Además, para caracterizar estas magnitudes no es necesario resolver la ecuación (3) hasta el estado estacionario, sino en una fina franja alrededor del frente de llama, lo cual ahorra mucho tiempo computacional. Se verá en el apartado siguiente que el vector unitario normal y la curvatura intervienen en la velocidad de propagación del frente de llama. Por este motivo, el método de level set proporciona una gran precisión en la simulación del proceso de combustión¹ .

${\overset {\mbox{ˆ}}{n}}={\frac {\nabla G}{\left|\nabla G\right|}}$

( 7)

$\kappa =-\nabla .{\overset {\mbox{ˆ}}{n}}=-\nabla .\left({\frac {\nabla G}{\left|\nabla G\right|}}\right)$

( 8)

3.2. Incorporación de la velocidad del frente de llama en la ecuación de level set

Muchos casos de combustión, incluido el del presente trabajo, operan en lo que se conoce como régimen flamelet, que consiste en que la química ocurre tan rápidamente que tiene lugar en finas capas del frente de llama llamadas flamelets. En este caso, la química puede desacoplarse del problema [8] , [9] , [10] , [11] , [12] , [13] , [14] , [15]  and [18] . En cuanto a la velocidad de propagación del frente de llama, para el caso de flujo laminar esta es la suma de la velocidad del fluido en el frente de llama y la velocidad de quemado en la dirección normal:

${\frac {d{\overset {\rightarrow }{x}}f}{dt}}={\overset {\rightarrow }{u}}f+$${\overset {\rightarrow }{n}}u_{L}$

( 9)

Sustituyendo la ecuación anterior en la ecuación (2), y teniendo en cuenta que la normal viene dada por la ecuación (7), se obtiene:

${\frac {\partial G}{\partial t}}+{\overset {\rightarrow }{u}}f.\nabla G=$$uL\left|\nabla G\right|$

( 10)

Markstein [19] indicó que la velocidad de llama viene dada por la siguiente expresión:

$u_{L}=u_{L}^{0}-D\kappa$

( 11)

siendo ${\textstyle {\mu }_{L}^{0}}$ la velocidad que tendría una llama plana (sin corrugar), obtenida mediante el modelo de Metghalchi y Keck [20] , D la difusividad de Markstein y κ la curvatura, dada por la ecuación (8). Sustituyendo la expresión de la velocidad de llama en la ecuación (10), se obtiene:

${\frac {\partial G}{\partial t}}+{\overset {\rightarrow }{u}}f.\nabla G=$$u_{L}^{0}\left|\nabla G\right|-D\kappa \left|\nabla G\right|$

( 12)

La ecuación anterior fue propuesta por Williams [21] para flujo laminar. Sin embargo, en un motor tiene lugar flujo turbulento, régimen para el cual Peters [22] propuso la siguiente ecuación:

${\frac {\partial G}{\partial t}}+{\overset {\rightarrow }{u}}f.\nabla G{\frac {{\rho }_{n}}{\rho }}u_{T}\left|\nabla G\right|-$$D_{T}\kappa \left|\nabla G\right|$

( 13)

siendo ρ_n la densidad de los gases no quemados, ρ   la densidad del gas en la llama definida por ${\textstyle G({\overset {\rightarrow }{x}},t)=}$$0$ y u_T la velocidad de llama turbulenta. La relación entre la velocidad de llama laminar y turbulenta que propuso Peters [22] viene dada por la siguiente expresión:

$u_{T}=u_{L}+{\upsilon }^{'}\left\{-{\frac {a_{4}b_{3}^{2}}{2b_{1}}}Da+\right.$$\left.{\left[{\left({\frac {a_{4}b_{3}^{2}}{2b_{1}}}Da\right)}^{2}+a_{4}b_{3}^{2}Da\right]}^{1/2}\right\}$

( 14)

siendo a₄ , b₁ , y b₃ constantes de valor 0,78, 2 y 1 respectivamente, υ^' la intensidad turbulenta y Da el número de DamKohler, dado por:

$Da={\frac {u_{L}l}{{\upsilon }^{'}l_{F}}}$

( 15)

siendo l la longitud de escala turbulenta y l_F el espesor de llama, dado por:

$l_{F}={\frac {\lambda /c_{p}}{\rho u_{L}}}$

( 16)

siendo λ la conductividad térmica y c_p el calor específico.

4. Implementación de la ecuación de level set en un código de mecánica de fluidos computacional

El procedimiento de level set descrito en el apartado anterior se ha implementado en un software de CFD. La simulación se ha llevado a cabo desde 30° antes del punto superior (con el pistón ascendiendo) hasta 30° después del mismo (con el pistón descendiendo), intervalo de tiempo para el cual las lumbreras de transferencia y escape permanecen cerradas. A continuación se describen las ecuaciones gobernantes y posteriormente el procedimiento numérico empleado.

4.1. Ecuaciones gobernantes

Además de la ecuación de level set, ecuación (13), el proceso de combustión turbulenta en el motor se ha modelado mediante las ecuaciones Reynolds averaged Navier Stokes (RANS), junto con la ecuación de conservación de la energía. El flujo se ha tratado como compresible, concretamente como gas ideal.

En cuanto a la ecuación de conservación de la masa, esta es igual que para flujos no reaccionantes (la combustión no genera masa):

${\frac {\partial \rho }{\partial t}}+{\frac {\partial }{\partial x_{i}}}(\rho u_{i})=$$0$

( 17)

La ecuación de conservación de cantidad de movimiento es también igual que para flujos no reaccionantes:

${\frac {\partial }{\partial t}}(\rho u_{i})+{\frac {\partial }{\partial x_{j}}}(\rho u_{i}u_{j})=$$-{\frac {\partial p}{\partial x_{i}}}+{\frac {\partial {\tau }_{ij}}{\partial x_{j}}}{\frac {\partial }{\partial x_{j}}}(-$$\rho {\overline {{u^{'}}_{i}{u^{'}}_{j}}})$

( 18)

siendo τ_ij el tensor de tensiones y el término ${\textstyle \left(-\right.}$$\left.\rho {\overline {{u^{'}}_{i}{u^{'}}_{j}}}\right)$ una representación de las tensiones de Reynolds, que en el presente trabajo se han calculado utilizando el modelo de turbulencia k-ɛ debido a su robustez, economía computacional y razonable precisión para un gran rango de flujos turbulentos.

La ecuación de conservación de la energía viene dada por:

${\frac {\partial }{\partial t}}(\rho H)+{\frac {\partial }{\partial x_{i}}}(\rho u_{i}H)=$${\frac {\partial }{\partial x_{i}}}\left[\left({\frac {\mu }{\sigma }}+\right.\right.$$\left.\left.{\frac {{\mu }_{t}}{{\sigma }_{t}}}\right){\frac {\partial H}{\partial x_{i}}}\right]+$$S_{rad}$

( 19)

siendo H la entalpía, μ la viscosidad, μ_t la viscosidad turbulenta, σ el número de Prandtl, σ_t el número de Prandtl turbulento y S_rad un término fuente de entalpía debido a la transferencia de calor por radiación. En el presente trabajo se ha utilizado el modelo de radiación P1, empleando un valor del coeficiente de absorción de 0,2 m^-1 , Bose [23] . Otros trabajos que han utilizado el modelo P1 aplicado a combustión son los de Kontogeorgos et al. [24] , Estrada [25] , Copete et al. [26] y Liu et al. [27] .

La temperatura se ha obtenido a partir de la entalpía mediante la definición de esta última:

$H=H_{ref}+{\int }_{T_{ref}}^{T}c_{p}dT+{\alpha }_{p}H_{c}$

( 20)

siendo T la temperatura, T_ref la temperatura de referencia, H_ref la entalpía de referencia, c_p el calor específico, α_p la fracción volumétrica de los productos y H_c la entalpía de combustión.

4.2. Condiciones de contorno e iniciales

Todo modelo de CFD requiere condiciones iniciales y de contorno. Respecto a las condiciones de contorno, se ha asumido convección desde el cilindro hacia el agua de refrigeración:

$q=h(T_{gas}-T_{agua})$

( 21)

siendo q la tasa de transferencia de calor, T_gas la temperatura en el interior del cilindro, T_agua la temperatura del agua de refrigeración (para este caso se ha medido experimentalmente obteniendo un valor de 345,6 K) y h el coeficiente de transferencia de calor, dado por la siguiente expresión [28] :

$h=10,4kb^{-1/4}{\left(u_{piston}/\nu \right)}^{3/4}$

( 22)

siendo b la carrera, k la conductividad térmica, u_piston la velocidad media del pistón y ν la viscosidad cinemática. Sustituyendo valores en la ecuación anterior resulta un valor de h = 7.853 W/m² K.

Respecto a las condiciones iniciales, la presión en el instante correspondiente a 30° antes del PMS es de 3,85 bar, figura 3 . Por desgracia, actualmente no se dispone de medios para medir con la suficiente precisión la temperatura en el interior del cilindro debido a que el tiempo de respuesta de un sensor de temperatura no es lo suficientemente rápido como para seguir la velocidad del motor. Por este motivo, utilizando los datos experimentales de la presión se ha determinado el ciclo termodinámico ideal, a partir del cual se ha obtenido que la temperatura inicial es de 511 K. Detalles del procedimiento se pueden consultar en bibliografía de termodinámica o en [16] .

4.3. Procedimiento numérico

Las ecuaciones gobernantes se han resuelto mediante el software de CFD Open Field Operation and Manipulation (OpenFOAM). Se indican en el apéndice detalles de la implementación y una síntesis de los principales fragmentos del código. Se ha acudido al OpenFOAM porque, al tratarse de un software libre, permite una total manipulación del mismo, lo cual ha facilitado en gran medida la programación del problema.

La malla se muestra en la figura 7 . Se ha utilizado una malla deformable correspondiente al movimiento ascendente/descendente del pistón. Concretamente, en la figura se muestra la malla al comienzo y final de la simulación (−30° y 30° de ángulo de cigüeñal respectivamente), y la posición en el PMS. Se han empleado elementos hexaédricos, cuyo número ha variado desde aproximadamente 100.000 en la posición del PMS hasta 130.000 en la posición correspondiente a −30 y 30° de ángulo de cigüeñal.

Figura 7. Movimiento y adaptación de la malla para el comienzo y final de la simulación y el PMS.

El paso de tiempo ha sido equivalente a 0,1° de ángulo de cigüeñal, por lo cual han sido necesarios 600 pasos de tiempo para simular el recorrido desde −30° hasta 30°. Las ecuaciones gobernantes han sido discretizadas mediante un esquema de segundo orden y el tratamiento temporal ha sido de primer orden. El acoplamiento entre presión-velocidad se ha tratado mediante el algoritmo PISO. En cuanto al tamaño de malla y de paso de tiempo, se ha hecho un análisis de sensibilidad para verificar que los resultados obtenidos son prácticamente insensibles ante refinamientos de malla o tamaño de paso.

5. Resultados

Para el intervalo de tiempo considerado, es decir, desde 30° de ángulo de cigüeñal antes del PMS hasta 30° después del PMS, la presión media dentro del cilindro obtenida numérica y experimentalmente se muestra en la figura 8 . La ignición se produce 20° antes del PMS, momento en el cual el pistón se encuentra ascendiendo. La presión en el interior del cilindro asciende drásticamente a partir de la ignición. Tras alcanzar el PMS, la presión sigue ascendiendo durante unos instantes, pero pronto desciende debido a la expansión propiciada por el pistón en su carrera descendente. Esta caracterización de la distribución de presión proporciona información muy útil puesto que el rendimiento del motor es óptimo si el momento de máxima presión coincide con el PMS (es decir, cuando el volumen del cilindro se reduce al mínimo). Este hecho se consigue regulando el instante de salto de la chispa, para lo cual el modelo desarrollado en el presente trabajo es muy útil a la hora de predecir el comportamiento.

Figura 8. Presión en el interior del cilindro obtenida numérica y experimentalmente.

En cuanto a la precisión obtenida, se observa que el método de level set proporciona resultados bastante satisfactorios, siendo el error promedio 4,2%. Con el objetivo de cuantificar la mejoría obtenida al aplicar el método de level set, se ha repetido la simulación sin hacer el procedimiento de reinicialización. Como se puede observar en la figura 8 , los resultados para este caso proporcionan un error mayor, concretamente un 8,6%. Asimismo, también se ha repetido la simulación utilizando un método basado en lo que se conoce como una variable de progreso de reacción (reaction progress variable en inglés), implementado en un gran número de softwares comerciales como Fluent, StarCD, CFX, etc. En este método, para caracterización de la propagación del frente de llama también se define una magnitud escalar, la variable de progreso de reacción, pero en lugar de ser la distancia al frente de llama es la fracción volumétrica de uno de los componentes. Para este caso concreto se ha elegido como c = 0 en el componente no quemado y c = 1 en el quemado. En lugar de resolver la ecuación G , ecuación (13), la advección de c viene dada por la siguiente ecuación:

$\rho {\frac {\partial c}{\partial t}}+\rho {\overset {\rightarrow }{u}}_{f}\cdot \nabla c=$$\nabla \cdot \left({\frac {{\mu }_{t}}{{\sigma }_{t}}}\nabla c\right)+$${\rho }_{u}u_{T}\left|\nabla c\right|$

( 23)

siendo μ_t la viscosidad turbulenta y σ_t el número de Prandtl turbulento. Lógicamente, este procedimiento no da lugar a proceso de reinicialización. Los resultados obtenidos con este método también se indican en la figura 8 . Como en la escala de la gráfica no se puede apreciar el error obtenido con demasiada claridad, en la tabla 2 se indica el error promedio, en el valor máximo y en las pendientes de las etapas de compresión expansión. Se puede observar que el error cometido al utilizar el método de variable de progreso de reacción es mayor que el obtenido al utilizar el procedimiento de level set, y los resultados son muy parecidos a los obtenidos con el método de level set sin reinicialización. El motivo de que el método de level set sin reinicialización y el método de variable de progreso de reacción se alejen más de los resultados numéricos es debido a una menor imprecisión en la determinación de la curvatura y el gradiente de c (en el método de variable de progreso de reacción) o de G (en el método de level set sin reiniliciazión).

Tabla 2. Error en la presión promedio, en el valor máximo y en las pendientes de las etapas de compresión y expansión utilizando el método de level set con y sin reinicialización y el método de variable de progreso de reacción

Método	Level set con reinicialización	Variable de progreso de reacción	Level set sin reinicialización
Error (%)
Promedio	4,2	7,7	8,6
Valor máximo	2,5	7,6	9,5
Pendiente compresión	3,6	11,2	13,0
Pendiente expansión	1,25	1,5	1,9

El campo de temperaturas se muestra en la figura 9 . Asimismo, se ha representado el frente de llama mediante una línea que indica la isosuperficie donde G = 0. Como se puede observar, el frente de llama se propaga desde la bujía. La temperatura es muy elevada, llegando a superar los 2.500 K cuando el pistón está cerca del PMS. Sin embargo, a medida que el pistón desciende, los productos se enfrían debido a la expansión propiciada por el descenso del pistón. Desgraciadamente, los resultados del campo de temperaturas no se han podido validar con resultados experimentales debido a que los sensores de temperatura que existen actualmente no tienen un tiempo de respuesta lo suficientemente rápido como para registrar datos a lo largo del ciclo.

Figura 9. Campo de temperaturas (K) y frente de llama para varios valores de ángulo de cigüeñal.

Esta caracterización del campo de temperaturas es muy importante para cuantificar las emisiones de NO_x , las cuales se intensifican cuanto mayor sea la temperatura alcanzada en la combustión.

6. Conclusiones

En este trabajo se ha implementado el método numérico level set en un código de CFD para simular el proceso de combustión en un motor Otto de dos tiempos. El procedimiento de level set que se ha programado consiste en resolver la ecuación G junto a su reinicialización. Se han obtenido resultados de la propagación del frente de llama, campo de temperaturas y campo de presiones. El campo de presiones obtenido numéricamente se ha comparado con valores experimentales y con otro método numérico, obteniendo mejores resultados mediante el procedimiento de level set. Asimismo, se han mostrado los resultados que proporciona el método de level set sin el proceso de reinicialización.

La información obtenida a partir de este trabajo es muy útil para diseñar nuevos modelos de motores. El campo de presiones está relacionado con el rendimiento del motor y el campo de temperaturas con las emisiones de NO_x . En cuanto a trabajos futuros, es importante ampliar el código programado para cuantificar las emisiones. Otro trabajo futuro muy interesante es el estudio de la optimización de la combustión variando parámetros como la geometría de la cámara de combustión, presiones de trabajo, relación de compresión, etc. Un aspecto también muy interesante sería el estudio de la reducción de las emisiones de NO_x por medio de una disminución de la temperatura mediante recirculación de los gases de escape.

Apéndice 1. Programación de las ecuaciones gobernantes

El software de CFD empleado en el presente trabajo, OpenFOAM, es básicamente una biblioteca de C++ para resolver numéricamente problemas de mecánica de fluidos. Como el método de level set no viene implementado por defecto en el software, se ha programado a partir de otros solvers, sobretodo el «xiFoam», (utilizado para simular el proceso de combustión premezclada), creando con ello un código propio. Se ha utilizado la ecuación G y su reiniliciación, ecuaciones (13) y (5) para obtener el campo de level set. La densidad se ha obtenido a partir de la ecuación de conservación de la masa, ecuación (17), la temperatura a partir de la ecuación de conservación de la energía y definición de entalpía, ecuaciones (19) y (20), la presión a partir de la ecuación de estado (P = ρRT ) y la velocidad a partir de la ecuación de conservación de la cantidad de movimiento, ecuación (18).

La ecuación G , dada por la Ec. (13), se ha programado de la siguiente manera:

• solve
• (
• fvm::ddt(G)
• + fvm::SuSp(fvm::grad(G),phi)
• ==
• (rhouT-Dk)*fvm::grad(G)
• );

La ecuación de reinicialización, dada por la ecuación (5), se ha programado desarrollando la derivada en el tiempo artificial mediante el método de Euler:

• GCorrNew = G;
• for (label i = 0; i < imaximum; i++)
• {
• GCorr = GCorrNew;
• surfaceVectorField gradGCorrf = fvc::interpolate(fvc::grad(GCorr));
• surfaceVectorField w = sign(G)* gradGCorrf/(mag(gradGCorrf) + SMALL);
• GCorrNew = GCorr + step*(-fvc::grad(GCorr), w) + sign(G));
• }
• G = GCorrNew;

La ecuación de conservación de la masa, ecuación (17) , se ha programado de la siguiente manera:

• solve
• (
• fvm::ddt(rho)
• + fvm::div(rho,phi)
• );

La ecuación de conservación de la energía, ecuación (19), se ha programado de la siguiente manera, en la cual se ha llamado ${\textstyle \Gamma =}$${\frac {\mu }{\sigma }}+{\frac {{\mu }_{t}}{{\sigma }_{h}}}$ :

• volScalarField gamma
• (
• «gamma» ,
• turbulence->nu()/Pr + turbulence->nut()/Prt
• );
• solve
• (
• fvm::ddt(rho,h)
• + fvm::div(rhophi,h)
• fvm::laplacian(gamma,h)
• ==
• radiation->Sh(thermo)
• );

La ecuación de conservación de cantidad de movimiento, ecuación (18) , se ha reordenado como sigue:

${\frac {\partial }{\partial t}}(\rho u_{i})+{\frac {\partial }{\partial x_{j}}}(\rho u_{i}u_{j})+$${\frac {\partial }{\partial x_{j}}}(\rho {\overline {{u^{'}}_{i}{u^{'}}_{j}}}-$${\tau }_{ij})=-{\frac {\partial p}{\partial x_{i}}}$

( A.1)

la cual se ha programado de la siguiente manera:

• solve
• (
• fvm::ddt(rho,U)
• + fvm::div(rhophi, U)
• + turbulence->divDevRhoReff(U)
• ==
• fvc::grad(p)
• );

Bibliografía

1. [1] S. Osher, J.A. Sethian; Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations; J Comput Phys, 79 (1988), pp. 12–49
2. [2] G. Son, V.K. Dhir; Numerical simulation of film boiling near critical pressures with a level set method; J Heat Trans, 120 (1998), pp. 183–192
3. [3] S. Tanguy, T. Ménard, A. Berlemont; A level set method for vaporizing two-phase flows; J Comput Phys, 221 (2007), pp. 837–853
4. [4] N. Sukumar, D. Chopp, B. Moran; Extended finite element method and fast marching method for three-dimensional fatigue crack propagation; J Mechanics Phys Solids, 7 (2003), pp. 29–48
5. [5] D. Chambers, D. Marcus, M. Sussman; Relaxation spectra of surface waves; En Proceedings of the 1995 International Mechanical Engineering Congress and Exposition (1995), p. 1995
6. [6] L. Tan, N. Zabaras; A level set simulation of dentritic solidification of multi-component alloys; J Comput Phys, 221 (2007), pp. 9–40
7. [7] B. Li, J. Shopple; An interface-fitted finite element level set method with application to solidification and solvation; Commun Comput Phys, 10 (2011), pp. 32–56
8. [8] V. Smiljanovski, V. Moser, R. Klein; A capturing-tracking hybrid scheme for deflagration discontinuities; J Combust Theory Model, 2 (1997), pp. 183–215
9. [9] R. Fedkiw, T. Islam, S. Xu; The ghost fluid method for deflagration and detonation discontinuities; J Comput Physics, 154 (1999), pp. 393–427
10. [10] R. Dahms, N. Peters, D.W. Stanton, Z. Tan, J. Ewald; Pollutant formation modelling in natural gas SI engines using a level set based flamelet model; Int J Eng Res, 9 (2008), pp. 1–14
11. [11] S. Yang, R.D. Reitz; Improved combustion submodels for modelling gasoline engines with the level set G equation and detailed chemical kinetics; P I Mech Eng D-J Aut, 223 (2009), pp. 702–726
12. [12] S. Singh, L. Liang, S.C. Kong, R.D. Reitz; Development of a flame propagation model for dual-fuel partially premixed compression ignition engines; Int J Eng Res, 7 (2006), pp. 64–73
13. [13] L. Liang, R.D. Deitz; Spark ignition engine combustion modelling using a level set method with detailed chemistry; SAE paper 2006 (2006), pp. 01–2043
14. [14] Z. Tan, R.D. Reitz; An ignition and combustion model based on the level-set method for spark ignition engine multidimensional modelling; Combust Flame, 145 (2006), pp. 1–15
15. [15] J. Ewald; A level set based flamelet model for the prediction of combustion in homogenous charge and direct injection spark ignition engines; (1^st Edition)Curvillier Verlag, Göttingen (2006)
16. [16] M. Lamas-Galdo, C. Rodríguez-Vidal, J. Rodríguez-García, M. Fernández-Quintás; Modelo de Mecánica de Fluidos Computacional para el proceso de barrido en un motor Otto de dos tiempos; DYNA Ingeniería e Industria, 86-2 (2011), pp. 165–172
17. [17] M. Sussmann, P. Smereka, S. Osher; A level set approach for computing solutions to incompressible two-phase flows; J Comput Phys, 114 (1994), pp. 146–159
18. [18] H.K. Versteeg, W. Malalasekera; An introduction to computational fluid dynamics the finite volume method; (2^nd Edition)Pearson Prentice Hall (2007)
19. [19] G. Markstein; Nonsteady Flame Propagation; Pergamon Press, Oxford (1964)
20. [20] M. Metghalchi, J.C. Keck; Burning velocities of mixtures of air with mechanol, isooctane, and indolene at high pressure and temperature; Combust flame, 48 (1982), pp. 191–210
21. [21] F.A. Williams; The mathematics of combustion; XIAM, Philadelphia PA (1985) 97–131
22. [22] N. Peters; Turbulent combustion; Cambridge University Press (2000)
23. [23] T.K. Bose; High temperature gas dynamics; (1^st Edition)Springer-Verlag Berlin Heidelberg, New York (2004)
24. [24] D.A. Kontogeorgos, E.P. Keramida, M.A. Founti; Assessment of simplified thermal radiation models for engineering calculations in natural gas-fired furnace; Int J Heat Mass Transfer, 50 (2007), pp. 5260–5268
25. [25] C.A. Estrada; Simulación de una cámara de combustión para una microturbina de gas utilizando el programa de dinámica de fluidos Fluent; Scientia et Technica, 34 (2007), pp. 255–260
26. [26] H. Copete, A. Amell, F. Cadavid; Simulación numérica de una cámara de combustion de alta velocidad con dos configuraciones de inyección de combustible; Dyna, 156 (2008), pp. 109–120
27. [27] D. Liu, F. Ding, H. Zhang, W. Zheng; Numerical simulation of high temperature air combustion in aluminium hydroxide gas suspension calcinations; Trans Nonferrous Met Soc China, 19 (2009), pp. 259–266
28. [28] C.F. Taylor; The internal combustion engine in theory and practice; (2^nd Edition)MIT Press (1985)

Notes

1. Nótese que una consecuencia de que la función level set sea considerada una distancia es que el módulo de G resulta unitario. Sin embargo, en las Ecs. (7) y (8) se ha utilizado la definición estricta de curvatura. En el resto del artículo se ha trabajado con el módulo de G calculado puesto que la level set solo se ha reinicializado una franja alrededor del frente de llama en lugar de la totalidad del dominio.