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 [[#bib0005|[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 [[#bib0010|[2]]]  and [[#bib0015|[3]]] , grietas en materiales [[#bib0020|[4]]] , movimiento de olas debido al viento [[#bib0025|[5]]] , solidificación [[#bib0030|[6]]]  and [[#bib0035|[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 [[#bib0040|[8]]]  and [[#bib0045|[9]]] , y en cuanto a combustión en motores, se pueden encontrar en la literatura los trabajos de [[#bib0050|[10]]] , [[#bib0055|[11]]] , [[#bib0060|[12]]] , [[#bib0065|[13]]] , [[#bib0070|[14]]]  and [[#bib0075|[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 [[#bib0080|[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 [[#tbl0005|tabla 1]] .

<span id='tbl0005'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

Tabla 1.

Especificaciones técnicas

|-

! Parámetro
! Valor
|-

| Tipo de motor
| Dos tiempos, Otto
|-

| Potencia (kW)
| 7,5
|-

| Velocidad (rpm)
| 6.500
|-

| Cilindrada (cm<sup>3</sup> )                                                    
| 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[[#tblfn0005|<sup>*</sup>]]
|}

*. 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), [[#fig0005|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, [[#fig0005|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, [[#fig0005|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, [[#fig0005|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.

<span id='fig0005'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_354709251-1-s2.0-S0213131512000338-gr1.jpg|center|393px|Funcionamiento del motor.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 1.

Funcionamiento del motor.

</span>
|}

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, [[#fig0005|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 [[#fig0010|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.

<span id='fig0010'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_354709251-1-s2.0-S0213131512000338-gr2.jpg|center|289px|Inicio de la combustión en un motor de encendido por chispa.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 2.

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

</span>
|}

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 [[#fig0015|figura 3]] .

<span id='fig0015'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_354709251-1-s2.0-S0213131512000338-gr3.jpg|center|335px|Presión en el interior del cilindro medida experimentalmente.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 3.

Presión en el interior del cilindro medida experimentalmente.

</span>
|}

==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  [[#fig0020|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.

<span id='fig0020'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_354709251-1-s2.0-S0213131512000338-gr4.jpg|center|249px|Frente de llama e isocontornos de level set alrededor del mismo.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 4.

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

</span>
|}

Supóngase que ''x''<sub>''f''</sub>''(t)''  es el trayecto de un punto en el frente de llama, es decir, donde la distancia es cero:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>G(x_f(t),t)=0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 1)
|}

Derivando respecto al tiempo resulta la siguiente expresión:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial G}{\partial t}+\nabla G\cdot \frac{d{\overset{\rightarrow}{x}}_f}{dt}=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 2)
|}

Esta última ecuación fue la introducida por Osher y Sethian [[#bib0005|[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 [[#fig0025|figura 5]]  (a) y (b) para los instantes de tiempo 1 s y 5 s respectivamente.

<span id='fig0025'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_354709251-1-s2.0-S0213131512000338-gr5.jpg|center|301px|Crecimiento de una esfera de radio inicial 1 m y velocidad de crecimiento radial ...]]


|-
| <span style="text-align: center; font-size: 75%;">

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.

</span>
|}

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.  [[#bib0085|[17]]]  introdujeron un procedimiento iterativo para reinicializar ''G'' , el cual se basa en la siguiente ecuación:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial G}{\partial \tau }=sign(G_0)\left(1-\right. </math><math>\left. \left|\nabla G\right|\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 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''<sub>''0''</sub>  el campo de la función level set al principio de cada iteración en este tiempo artificial introducido y sign(''G''<sub>0</sub> ) el signo de ''G''<sub>''0''</sub> , dado por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>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\}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 4)
|}

Como el signo vale cero cuando la función level set es cero, ''G''  tiene el mismo cero que ''G''<sub>''0''</sub> . 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:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial G}{\partial \tau }+\overset{\rightarrow}{w}\cdot \nabla G=</math><math>sign(G_0)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 5)
|}

siendo:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\overset{\rightarrow}{w}=sign(G_0)\frac{\nabla G}{\left|\nabla G\right|}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 6)
|}

El vector <math display="inline">\overset{\rightarrow}{w}</math>  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  [[#fig0030|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  [[#fig0030|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 [[#fig0030|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 [[#fig0030|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 [[#fig0030|figura 6]] ) a cada lado del frente de llama.

<span id='fig0030'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_354709251-1-s2.0-S0213131512000338-gr6.jpg|center|299px|Procedimiento de reinicialización de la level set para una esfera.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 6.

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

</span>
|}

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[[#fn0005|<sup>1</sup>]] .

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\overset{\mbox{ˆ}}{n}=\frac{\nabla G}{\left|\nabla G\right|}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 7)
|}

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\kappa =-\nabla .\overset{\mbox{ˆ}}{n}=-\nabla .\left(\frac{\nabla G}{\left|\nabla G\right|}\right)</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 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 [[#bib0040|[8]]] , [[#bib0045|[9]]] , [[#bib0050|[10]]] , [[#bib0055|[11]]] , [[#bib0060|[12]]] , [[#bib0065|[13]]] , [[#bib0070|[14]]] , [[#bib0075|[15]]]  and [[#bib0090|[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:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{d\overset{\rightarrow}{x}f}{dt}=\overset{\rightarrow}{u}f+</math><math>\overset{\rightarrow}{n}u_L</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 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:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial G}{\partial t}+\overset{\rightarrow}{u}f.\nabla G=</math><math>uL\left|\nabla G\right|</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 10)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>u_L=u_L^0-D\kappa </math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 11)
|}

siendo <math display="inline">{\mu }_L^0</math>  la velocidad que tendría una llama plana (sin corrugar), obtenida mediante el modelo de Metghalchi y Keck [[#bib0100|[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:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial G}{\partial t}+\overset{\rightarrow}{u}f.\nabla G=</math><math>u_L^0\left|\nabla G\right|-D\kappa \left|\nabla G\right|</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 12)
|}

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial G}{\partial t}+\overset{\rightarrow}{u}f.\nabla G\frac{{\rho }_n}{\rho }u_T\left|\nabla G\right|-</math><math>D_T\kappa \left|\nabla G\right|</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 13)
|}

siendo ''ρ''<sub>''n''</sub>  la densidad de los gases no quemados, ''ρ  ''  la densidad del gas en la llama definida por <math display="inline">G(\overset{\rightarrow}{x},t)=</math><math>0</math>  y ''u''<sub>''T''</sub>  la velocidad de llama turbulenta. La relación entre la velocidad de llama laminar y turbulenta que propuso Peters [[#bib0110|[22]]]  viene dada por la siguiente expresión:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>u_T=u_L+{\upsilon }^'\left\{-\frac{a_4b_3^2}{2b_1}Da+\right. </math><math>\left. {\left[{\left(\frac{a_4b_3^2}{2b_1}Da\right)}^2+\right. \right. </math><math>\left. \left. a_4b_3^2Da\right]}^{1/2}\right\}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 14)
|}

siendo ''a''<sub>''4''</sub> , ''b''<sub>''1''</sub> , y ''b''<sub>''3''</sub>  constantes de valor 0,78, 2 y 1 respectivamente, ''υ''<sup>'</sup>  la intensidad turbulenta y ''Da''  el número de DamKohler, dado por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>Da=\frac{u_Ll}{{\upsilon }^'l_F}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 15)
|}

siendo ''l''  la longitud de escala turbulenta y ''l''<sub>''F''</sub>  el espesor de llama, dado por:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>l_F=\frac{\lambda /c_p}{\rho u_L}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 16)
|}

siendo ''λ''  la conductividad térmica y ''c''<sub>''p''</sub>  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):

<span id='eq0085'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial \rho }{\partial t}+\frac{\partial }{\partial x_i}(\rho u_i)=</math><math>0</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 17)
|}

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

<span id='eq0090'></span>
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial }{\partial t}(\rho u_i)+\frac{\partial }{\partial x_j}(\rho u_iu_j)=</math><math>-\frac{\partial p}{\partial x_i}+\frac{\partial {\tau }_{ij}}{\partial x_j}\frac{\partial }{\partial x_j}(-</math><math>\rho \overline{{u^{{'}}}_i{u^{{'}}}_j})</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 18)
|}

siendo ''τ''<sub>ij</sub>  el tensor de tensiones y el término <math display="inline">\left(-\right. </math><math>\left. \rho \overline{{u^{{'}}}_i{u^{{'}}}_j}\right)</math>  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:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial }{\partial t}(\rho H)+\frac{\partial }{\partial x_i}(\rho u_iH)=</math><math>\frac{\partial }{\partial x_i}\left[\left(\frac{\mu }{\sigma }+\right. \right. </math><math>\left. \left. \frac{{\mu }_t}{{\sigma }_t}\right)\frac{\partial H}{\partial x_i}\right]+</math><math>S_{rad}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 19)
|}

siendo ''H''  la entalpía, ''μ''  la viscosidad, ''μ''<sub>''t''</sub>  la viscosidad turbulenta, ''σ''  el número de Prandtl, ''σ''<sub>''t''</sub>  el número de Prandtl turbulento y ''S''<sub>''rad''</sub>  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<sup>-1</sup> , Bose [[#bib0115|[23]]] . Otros trabajos que han utilizado el modelo P1 aplicado a combustión son los de Kontogeorgos et al. [[#bib0120|[24]]] , Estrada [[#bib0145|[25]]] , Copete et al. [[#bib0130|[26]]]  y Liu et al. [[#bib0135|[27]]] .      

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

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>H=H_{ref}+{\int }_{T_{ref}}^Tc_pdT+{\alpha }_pH_c</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 20)
|}

siendo ''T''  la temperatura, ''T''<sub>''ref''</sub>  la temperatura de referencia, ''H''<sub>''ref''</sub>  la entalpía de referencia, ''c''<sub>''p''</sub>  el calor específico, ''α''<sub>''p''</sub>  la fracción volumétrica de los productos y ''H''<sub>''c''</sub>  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:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>q=h(T_{gas}-T_{agua})</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 21)
|}

siendo ''q''  la tasa de transferencia de calor, ''T''<sub>''gas''</sub>  la temperatura en el interior del cilindro, ''T''<sub>''agua''</sub>  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  [[#bib0140|[28]]] :

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>h=10,4kb^{-1/4}{\left(u_{piston}/\nu \right)}^{3/4}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 22)
|}

siendo ''b''  la carrera, ''k''  la conductividad térmica, ''u''<sub>''piston''</sub>  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<sup>2</sup> K.      

Respecto a las condiciones iniciales, la presión en el instante correspondiente a 30° antes del PMS es de 3,85 bar, [[#fig0015|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 [[#bib0080|[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 [[#fig0035|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.

<span id='fig0035'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_354709251-1-s2.0-S0213131512000338-gr7.jpg|center|355px|Movimiento y adaptación de la malla para el comienzo y final de la simulación y ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 7.

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

</span>
|}

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 [[#fig0040|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.

<span id='fig0040'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_354709251-1-s2.0-S0213131512000338-gr8.jpg|center|347px|Presión en el interior del cilindro obtenida numérica y experimentalmente.]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 8.

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

</span>
|}

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 [[#fig0040|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:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\rho \frac{\partial c}{\partial t}+\rho {\overset{\rightarrow}{u}}_f\cdot \nabla c=</math><math>\nabla \cdot \left(\frac{{\mu }_t}{{\sigma }_t}\nabla c\right)+</math><math>{\rho }_uu_T\left|\nabla c\right|</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 23)
|}

siendo ''μ''<sub>''t''</sub>  la viscosidad turbulenta y ''σ''<sub>''t''</sub>  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 [[#fig0040|figura 8]] . Como en la escala de la gráfica no se puede apreciar el error obtenido con demasiada claridad, en la [[#tbl0010|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).

<span id='tbl0010'></span>

{| class="wikitable" style="min-width: 60%;margin-left: auto; margin-right: auto;"
|+

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 [[#fig0045|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.

<span id='fig0045'></span>

{| style="text-align: center; border: 1px solid #BBB; margin: 1em auto; max-width: 100%;" 
|-
|


[[Image:draft_Content_354709251-1-s2.0-S0213131512000338-gr9.jpg|center|301px|Campo de temperaturas (K) y frente de llama para varios valores de ángulo de ...]]


|-
| <span style="text-align: center; font-size: 75%;">

Figura 9.

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

</span>
|}

Esta caracterización del campo de temperaturas es muy importante para cuantificar las emisiones de NO<sub>x</sub> , 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<sub>x</sub> . 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<sub>x</sub>  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 [[#eq0085|(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 <math display="inline">\Gamma =</math><math>\frac{\mu }{\sigma }+\frac{{\mu }_t}{{\sigma }_h}</math> :
* ''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 [[#eq0090|(18)]] , se ha reordenado como sigue:

{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;" 
|-
| <math>\frac{\partial }{\partial t}(\rho u_i)+\frac{\partial }{\partial x_j}(\rho u_iu_j)+</math><math>\frac{\partial }{\partial x_j}(\rho \overline{{u^{{'}}}_i{u^{{'}}}_j}-</math><math>{\tau }_{ij})=-\frac{\partial p}{\partial x_i}</math>
|}
| style="width: 5px;text-align: right;white-space: nowrap;" | ( 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==

<ol style='list-style-type: none;margin-left: 0px;'><li><span id='bib0005'></span>
[[#bib0005|[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</li>
<li><span id='bib0010'></span>
[[#bib0010|[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</li>
<li><span id='bib0015'></span>
[[#bib0015|[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</li>
<li><span id='bib0020'></span>
[[#bib0020|[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</li>
<li><span id='bib0025'></span>
[[#bib0025|[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</li>
<li><span id='bib0030'></span>
[[#bib0030|[6]]] L. Tan, N. Zabaras; A level set simulation of dentritic solidification of multi-component alloys; J Comput Phys, 221 (2007), pp. 9–40</li>
<li><span id='bib0035'></span>
[[#bib0035|[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</li>
<li><span id='bib0040'></span>
[[#bib0040|[8]]] V. Smiljanovski, V. Moser, R. Klein; A capturing-tracking hybrid scheme for deflagration discontinuities; J Combust Theory Model, 2 (1997), pp. 183–215</li>
<li><span id='bib0045'></span>
[[#bib0045|[9]]] R. Fedkiw, T. Islam, S. Xu; The ghost fluid method for deflagration and detonation discontinuities; J Comput Physics, 154 (1999), pp. 393–427</li>
<li><span id='bib0050'></span>
[[#bib0050|[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</li>
<li><span id='bib0055'></span>
[[#bib0055|[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</li>
<li><span id='bib0060'></span>
[[#bib0060|[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</li>
<li><span id='bib0065'></span>
[[#bib0065|[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</li>
<li><span id='bib0070'></span>
[[#bib0070|[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</li>
<li><span id='bib0075'></span>
[[#bib0075|[15]]] J. Ewald; A level set based flamelet model for the prediction of combustion in homogenous charge and direct injection spark ignition engines;  (1<sup>st</sup>  Edition)Curvillier Verlag, Göttingen (2006)                                        </li>
<li><span id='bib0080'></span>
[[#bib0080|[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</li>
<li><span id='bib0085'></span>
[[#bib0085|[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</li>
<li><span id='bib0090'></span>
[[#bib0090|[18]]] H.K. Versteeg, W. Malalasekera; An introduction to computational fluid dynamics the finite volume method;  (2<sup>nd</sup>  Edition)Pearson Prentice Hall (2007)                                        </li>
<li><span id='bib0095'></span>
[[#bib0095|[19]]] G. Markstein; Nonsteady Flame Propagation; Pergamon Press, Oxford (1964)</li>
<li><span id='bib0100'></span>
[[#bib0100|[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</li>
<li><span id='bib0105'></span>
[[#bib0105|[21]]] F.A. Williams; The mathematics of combustion; XIAM, Philadelphia PA (1985) 97–131</li>
<li><span id='bib0110'></span>
[[#bib0110|[22]]] N. Peters; Turbulent combustion; Cambridge University Press (2000)</li>
<li><span id='bib0115'></span>
[[#bib0115|[23]]] T.K. Bose; High temperature gas dynamics;  (1<sup>st</sup>  Edition)Springer-Verlag Berlin Heidelberg, New York (2004)                                        </li>
<li><span id='bib0120'></span>
[[#bib0120|[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</li>
<li><span id='bib0145'></span>
[[#bib0145|[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</li>
<li><span id='bib0130'></span>
[[#bib0130|[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</li>
<li><span id='bib0135'></span>
[[#bib0135|[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</li>
<li><span id='bib0140'></span>
[[#bib0140|[28]]] C.F. Taylor; The internal combustion engine in theory and practice;  (2<sup>nd</sup>  Edition)MIT Press (1985)                                        </li>
</ol>

== Notes ==

<span id='fn0005'></span>

[[#bfn0005|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.