Resolución eficiente de modelos de combustión

Trabajo Fin de Grado Resolución eficiente de modelos de combustión Héctor Manuel de los Santos Porto 2018/2019 UNIVERSIDAD DE SANTIAGO DE COMPOSTELA GRADO DE MATEMÁTICAS Trabajo Fin de Grado Resolución eficiente de modelos de combustión Héctor Manuel de los Santos Porto Julio, 2019 UNIVERSIDAD DE SANTIAGO DE COMPOSTELA Trabajo propuesto Área de Conocimiento: Matemática Aplicada Título: Resolución eficiente de modelos de combustión Breve descripción del contenido: En este trabajo se propone la resolución mediante métodos numéricos de las ecuaciones que componen el modelo de combustión BFL compuesto por ecuaciones diferenciales ordinarias y ecuaciones algebraicas. Además, se comparará la eficiencia del mismo comparando los resultados obtenidos con los utilizados para la validación del modelo matemático publicado en Combustion and Flame. Recomendaciones: Buenas capacidades de programación. Otras observaciones: iii Índice general Índice de figuras vii Índice de cuadros ix Resumen xi Introducción xiii 1. Modelo de combustión 1 1.1. Modelo de combustión del carbón . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2. Modelo de la fase sólida . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 1.2.1. Modelo de gasificación de la partícula . . . . . . . . . . . . . . . . . 4 1.2.2. Fuentes homogéneas de la fase gaseosa . . . . . . . . . . . . . . . . . 8 1.2.3. Modelo del movimiento de la partícula . . . . . . . . . . . . . . . . . 10 2. Métodos numéricos 11 2.1. Problema de valor inicial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2. Métodosde1-paso ................................ 12 2.2.1. MétodosdeEuler............................. 13 2.2.2. Métodos Runge-Kutta . . . . . . . . . . . . . . . . . . . . . . . . . . 14 2.2.3. Métodos Runge-Kutta adaptativos . . . . . . . . . . . . . . . . . . . 17 2.3. MétododeNewton................................ 20 2.3.1. Método de Newton para funciones de una variable . . . . . . . . . . 20 2.3.2. Método de Newton para funciones de varias variables . . . . . . . . . 22 2.3.3. Método de Newton discretizado para funciones de varias variables . . 23 3. Resultados numéricos 25 3.1. Problematest................................... 25 3.2. Resultados..................................... 27 v vi ÍNDICE GENERAL 3.2.1. Esquemas de paso fijo . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.2.1.1. Esquema 1: Sistema desacoplado, RK4 . . . . . . . . . . . . 28 3.2.1.2. Esquema 2: Sistema acoplado, Euler explícito . . . . . . . . 33 3.2.1.3. Esquema 3: Sistema acoplado, Euler implícito . . . . . . . . 36 3.2.1.4. Esquema 4: Sistema acoplado, RK4 . . . . . . . . . . . . . 38 3.2.2. Esquemas de paso adaptativo . . . . . . . . . . . . . . . . . . . . . . 40 3.2.2.1. Esquema 5: Sistema desacoplado, RKF45 . . . . . . . . . . 41 3.2.2.2. Esquema 6: Sistema acoplado, RK4 . . . . . . . . . . . . . 44 4. Conclusiones 49 Nomenclatura 51 A. Datos del problema test 53 Bibliografía 55 Índice de figuras 1.1. Llamadedifusión. ................................ 2 1.2. Esquema de la combustión. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4 3.1. Diagrama de flujo del esquema 1. . . . . . . . . . . . . . . . . . . . . . . . . 29 3.2. Aproximaciones de ρVyλ4............................ 29 3.3. Aproximaciones de ρH2Oyλ5........................... 30 3.4. Aproximaciones de λ1,λ2yλ3del esquema 1. . . . . . . . . . . . . . . . . . 31 3.5. Aproximación de rcdelesquema1........................ 31 3.6. Aproximación de Tpdelesquema1........................ 32 3.7. Diagrama de flujo para los esquemas 2, 3 y 4. . . . . . . . . . . . . . . . . . 33 3.8. Aproximaciones de λ1,λ2yλ3del esquema 2. . . . . . . . . . . . . . . . . . 34 3.9. Aproximación de rcdelesquema2........................ 34 3.10. Aproximación de Tpdelesquema2........................ 35 3.11. Aproximaciones de λ1,λ2yλ3del esquema 3. . . . . . . . . . . . . . . . . . 36 3.12. Aproximación de rcdelesquema3........................ 36 3.13. Aproximación de Tpdelesquema3........................ 37 3.14. Aproximaciones de λ1,λ2yλ3del esquema 4. . . . . . . . . . . . . . . . . . 38 3.15. Aproximación de rcdelesquema4........................ 38 3.16. Aproximación de Tpdelesquema4........................ 39 3.17. Diagrama de flujo de los esquemas 5 y 6. . . . . . . . . . . . . . . . . . . . . 40 3.18. Aproximaciones de ρVyλ4delesquema5.................... 41 3.19. Aproximaciones de ρH2Oyλ5del esquema 5. . . . . . . . . . . . . . . . . . 42 3.20. Aproximaciones de λ1,λ2yλ3del esquema 5. . . . . . . . . . . . . . . . . . 42 3.21. Aproximación de rcdelesquema5........................ 43 3.22. Aproximación de Tpdelesquema5........................ 43 3.23. Aproximaciones de ρVyλ4delesquema6.................... 44 3.24. Aproximaciones de ρH2Oyλ5del esquema 6. . . . . . . . . . . . . . . . . . 45 vii xiv INTRODUCCIÓN teniendo especial cuidado en algunas hipótesis cruciales del modelo. A continuación se describen las ecuaciones que modelan matemáticamente los diferentes procesos que ocurren en la caldera para la fase sólida. En el segundo capítulo se expondrán los métodos numéricos que emplearemos en la resolución del modelo BFLs1 descrito en el capítulo anterior. La resolución de las ecuaciones diferenciales ordinarias (ODE’s) que forman el modelo se lleva a cabo mediante el método Runge-Kutta clásico (paso fijo) o mediante el método Runge-Kutta-Fehlberg (paso adaptativo). Para el sistema de ecuaciones algebraico-diferenciales (DAE’s) realizaremos una semidiscretización en tiempo, utilizando un método de Euler explícito, implícito o el Runge-Kutta de orden 4, para posteriormente resolver el sistema de ecuaciones no lineales utilizando el método de Newton discretizado. En el tercer capítulo se muestran los resultados numéricos obtenidos al resolver el problema test utilizando la metodología de soluciones manufacturadas, es decir, añadiendo al problema real términos fuente para que la solución exacta sea la que se propone. A partir de los datos especificados en el apéndice, hemos resuelto el problema con los distintos métodos numéricos y, apoyándonos en la gráficas incluidas, discutiremos su comportamiento para el test. Los diferentes métodos numéricos han sido implementados en el software comercial MatLab R2017a. Por último, el cuarto capítulo está dedicado a la redacción de las conclusiones obtenidas en este trabajo de Fin de Grado. Capítulo 1 Modelo de combustión En este primer capítulo introduciremos el modelo matemático, referido como BFLs1 en [2], que describe los fenómenos físico-químicos que ocurren en la vecindad de una partícula de carbón pulverizado durante su combustión en una caldera de una Central Térmica. La validez de este modelo se limita al caso que trata partículas de un tamaño pequeño (menor a100 µm de diámetro) y un alto contenido en cenizas. Nos centramos en este modelo por ser el que mejor se adapta al experimento que se describe y simula en [2]. 1.1. Modelo de combustión del carbón Comenzaremos explicando el proceso que ocurre en la caldera, donde se introducen las partículas de carbón pulverizado a través de distintos conductos. Las partículas son arrastradas por una corriente de gases de recirculación (aire primario) y se quemarán con el oxígeno aportado por el aire (aire secundario) que se introduce a través de otros conductos, formando una llama de difusión como se puede ver en el esquema de la Figura 1.11. Cuando éstas se encuentran dentro de la atmósfera de la caldera, comienzan a intercambiar calor con ésta y se producen las distintas reacciones que llevan a su combustión, las cuales expondremos en la siguiente sección. Dicha combustión da lugar a una llama de difusión en el dominio de la caldera (es una zona infinitamente delgada y que se ha denotado por ΓF) dividiéndolo en dos regiones: ΩF, que no contiene oxígeno y ΩOdonde sí hay oxígeno. 1Las figuras que aparecen en este capítulo las hemos obtenido del artículo [2] y del Trabajo de Fin de Grado [5] que aparecen en la bibliografía. 1 2CAPÍTULO 1. MODELO DE COMBUSTIÓN Figura 1.1: Llama de difusión. El modelo de combustión que estamos tratando se basa en una descripción lagrangiana para estudiar la evolución de cada partícula de carbón individualmente. Además, hay que tener en cuenta que la atmósfera de la caldera cambiará debido a la reacciones de volatilización y gasificación de las mismas. Los volátiles emitidos por las partículas, el monóxido de carbono y el hidrógeno resultantes de la gasificación reaccionarán con el oxígeno presente en la fase gaseosa para producir vapor de agua y dióxido de carbono. Una de las hipótesis cruciales del modelo BFLs1 es que las reacciones de oxidación de los volátiles, hidrógeno y monóxido de carbono que ocurren en la fase gaseosa son irreversibles e infinitamente rápidas en comparación con el resto de procesos que suceden en el horno con lo cual estas reacciones ocurrirán simultáneamente en la llama de difusión, ΓF. De esta manera se obtiene una buena explicación de la temperatura y de las concentraciones de las distintas especies. Por otro lado, el modelo se basa en un cálculo de la temperatura y densidad de cada partícula de carbón pulverizado a lo largo de su trayectoria. Entonces supondremos que las partículas no interactúan entre ellas, basándonos en el hecho de que se verifican las siguientes desigualdades: L >> lc>> lp>> a, donde Les la longitud del horno, lces la longitud de la celda computacional, lpes la distancia entre las partículas y ael radio de la partícula de carbón. Además, se supone que las partículas son esféricas y que las densidades de carbono fijo y cenizas de las mismas se mantienen constantes durante el proceso de combustión. Las reacciones de oxidación también pueden tener lugar en pequeñas llamas de difusión alrededor de las partículas o en el interior de las mismas en la región ΩO. Para que esto ocurra el tamaño de la partícula debe ser grande comparado con el espesor de la llama. 1.2. MODELO DE LA FASE SÓLIDA 3 Sin embargo, en el BFLs1 esto no es posible ya que las partículas son demasiado pequeñas como para mantenerla por si mismas, aunque si podrían mantenerla en un entorno en el que no afecte a las reacciones de gasificación. El carbono fijo se gasifica mediante las reacciones con el vapor de agua, el dióxido de carbono y el oxígeno que deben penetrar en las partículas desde la atmósfera gaseosa. Cada una de estas reacciones posee una temperatura de activación distinta, aunque en nuestro caso supondremos que difieren poco entre ellas y la denotaremos por Tc. Por debajo de ésta, podremos obviar las reacciones de gasificación, mientras que por encima suceden de manera instantánea por difusión. En las partículas con un alto contenido en cenizas, las reacciones de gasificación suceden en el núcleo de radio rc, menor que el radio original a. Las moléculas de dióxido de carbono, vapor de agua y oxígeno penetran a través de la capa de ceniza hasta alcanzar la capa de carbono fijo con una velocidad que determinará el decrecimiento del radio del núcleo. Por último, la capa de ceniza restante después de la gasificación mantiene el tamaño de la partícula de radio aconstante, puesto que es una especie inerte. 1.2. Modelo de la fase sólida Nuestro punto de partida será la exposición del modelo introducido en el artículo [2]. Basándonos en éste consideraremos un modelo cinético simplificado consistente en las siguientes reacciones físico-químicas que ocurren dentro de la partícula porosa: CO2(g) + C(s)→2CO(g)+(q1)(1.1) 1 2O2(g) + C(s)→CO(g)+(q2)(1.2) H2O(g) + C(s)→CO(g) + H2(g)+(q3)(1.3) V(s)→V(g)+(q4)(1.4) H2O(s)→H2O(g)+(q5)(1.5) y las siguientes reacciones de oxidación que ocurren en la fase gaseosa: CO(g) + 1 2O2(g)→CO2(g)+(q6)(1.6) V(g) + ν1O2(g)→ν2CO2(g) + ν3H2O(g) + ν4S2O(g)+(q7)(1.7) H2(g) + 1 2O2(g)→H2O(g)+(q8)(1.8) Observación 1.1.Nótese que estas tres últimas reacciones no se tratarán en este trabajo debido a que nos centraremos en el modelo de la fase sólida. 4CAPÍTULO 1. MODELO DE COMBUSTIÓN Notación 1.2.En las reacciones anteriores, hemos indicado que los reactivos y productos se encuentran en fase sólida o gaseosa mediante (s)y(g), respectivamente. Mientras que las reacciones (1.1),(1.3) y (1.5) son endotérmicas, la reacción (1.2) es exotérmica. Por otro lado, tomaremos q4= 0 debido a que desconocemos el comportamiento termoquímico de la reacción (1.4) por tratarse de una molécula idealizada, como describimos a continuación. Los volátiles considerados en las reacciones anteriores se toman como una única molécula de la siguiente manera: V=CK1HK2OK3SK4(1.9) cuya masa molecular se denota por Mvol y los coeficientes Kicon i= 1,...,4vienen determinados por un análisis del carbón. Los coeficientes estequiométricos de la reacción (1.7) se calculan a partir de los términos Kicomo se puede ver en [2]. A continuación se detallará el modelo que rige durante la combustión de una partícula de carbón, a través de la conservación de su masa y energía. 1.2.1. Modelo de gasificación de la partícula Dividiremos el modelo de la gasificación de la partícula de carbón en función de la temperatura de la misma. Por este motivo podemos describir este proceso en distintas etapas que se han representado en la Figura 1.2. V(g) H2O(g) H 2 CO H2O(g) O 2 CO 2 Tp>T1 Tp>Tc Y v =Y H2O =Y c =0 Figura 1.2: Esquema de la combustión. En primer lugar, la partícula entra en la caldera y comienza a intercambiar calor con el entorno, hasta alcanzar una temperatura T1que marca el inicio de la segunda etapa. Posteriormente, se dan las reacciones de evaporación y volatilización, al mismo tiempo que 1.2. MODELO DE LA FASE SÓLIDA 5 sigue incrementándose la temperatura de la partícula. La tercera etapa comienza cuando Tpes mayor que Tcpor lo que se darán las reacciones de gasificación del carbono fijo. La última etapa comienza cuando se han consumido todos los materiales combustibles y durante ella las cenizas restantes seguirán intercambiando calor hasta salir del horno. Por supuesto, todo el proceso se interrumpe si la partícula sale del horno. A continuación se describirán de forma exhaustiva cada una de estas etapas: 1. Primera etapa: en ella la temperatura de la partícula es menor que T1por lo que no se produce ninguna de las reacciones anteriormente expuestas. Por tanto, podemos tomar la velocidad de la i-ésima reacción λi= 0 con i= 1,...,5y el tamaño del radio invariante. Además, si suponemos que la temperatura de la partícula es la misma en todos sus puntos tenemos que su variación a lo largo del tiempo viene dada por: 4 3πa3ρpcs dTp dt = 4πa2(q00 p+q00 r),(1.10) donde los dos miembros de la parte derecha de la ecuación representan los flujos de calor que llegan a la superficie de la partícula por conducción y radiación, respectivamente. Los términos q00 pyq00 rse expresan de la siguiente forma: q00 p=kdT dr |r=a+=k a(Tg−cs cp Tp),(1.11) donde en esta última igualdad se ha utilizado la Ley de Fourier, y q00 r=εp1 4ZS2 I(x, ω)dω −σT4 p,(1.12) donde I(x, ω)es la intensidad de radiación en la dirección de ωen la posición de la partícula denotada por x,εpes la emisividad de la partícula y S2es la esfera unidad. 2. Segunda etapa: esta etapa comienza cuando la temperatura de la partícula es mayor o igual que T1. Diferenciaremos dos casos en función de la cantidad de volátiles y humedad presentes en la partícula. Si ambas cantidades denotadas por YVeYH2Oson nulas entonces λi= 0 con i= 4,5y resolveríamos de nuevo la ecuación para la temperatura (1.10). Por el contrario, si alguna de las dos cantidades es mayor que cero y asumiendo las desigualdades entre las distintas escalas anteriormente mencionadas, podemos modelar la evolución en el tiempo de las densidades de ambas componentes mediante las siguientes ecuaciones: dρV dt =−B4exp −E4 RTpρV,(1.13) 6CAPÍTULO 1. MODELO DE COMBUSTIÓN dρH2O dt =−B5exp −E5 RTpρH2O.(1.14) Estas relaciones de tipo Arrhenius se basan en las suposiciones de que las densidades de las diferentes especies y la temperatura en la partícula son uniformes y el volumen de la partícula se mantiene constante. Por otra parte, se introducen las velocidades adimensionales de reacción de los volátiles y la humedad, como la relación entre la masa de las componentes y el flujo de difusión característico, dadas por: λ4=a2 3ρgDB4exp −E4 RTpρV,(1.15) λ5=a2 3ρgDB5exp −E5 RTpρH2O.(1.16) Por último, podemos considerar que las reacciones de gasificación del carbono fijo están congeladas; en consecuencia no se produce monóxido de carbono ni hidrógeno en la partícula. De esta forma la variación de la temperatura en nuestro modelo BFLs1 se puede escribir de la siguiente forma: 4 3πa3ρpcs dTp dt = 4πa2(q00 p+q00 r)+4πρgaD(q4λ4+q5λ5),(1.17) donde q00 pse puede escribir como: q00 p=k acp (hg T−hs T)λ eλ−1,(1.18) yλ=λ4+λ5. Observación 1.3.Nótese que en esta última expresión se ha despreciado el efecto del movimiento relativo de la partícula y la atmósfera gaseosa en la vecindad de la partícula. 3. Tercera etapa: se caracteriza por el hecho de que la temperatura de la partícula es mayor que la temperatura crítica, Tc, por lo que ocurren, controladas por difusión, las reacciones de gasificación en la superficie de la partícula de carbón. Por tanto, YCO2=YO2=YH2O= 0 si r≤rc, siendo rcel radio del núcleo de la partícula de carbón. En otro caso, es decir, si r > rcentonces ρC= 0. En esta etapa, queremos determinar como varían las cantidades de monóxido de carbono, oxígeno y vapor de agua debidas a las reacciones (1.1),(1.2) y (1.3). Para ello estudiaremos sus velocidades asociadas, dadas por λicon i= 1,2,3, aunque su expresión no es tan simple como para las vistas anteriormente debido a que dichas 1.2. MODELO DE LA FASE SÓLIDA 7 especies no son uniformes en la partícula. Para describir sus efectos en la gasificación del carbono fijo contenido en la partícula se ha considerado una simplificación, que la energía de activación de las tres reacciones es la misma y además grande. Adicionalmente, se supuso que la densidad de carbono fijo se mantiene constante cuando ocurren las reacciones de gasificación mientras que, el volumen varía. Notación 1.4.Denotaremos por λ= 5 P i=1 λial caudal másico de gas generado en la partícula dividido por 4πaρgD. De manera similar a como hicimos en la segunda etapa diferenciaremos varios casos: Si la cantidad de carbono fijo en la partícula es cero entonces no se pueden dar dichas reacciones, y por ello resolveremos como en la segunda etapa. En cambio, si la cantidad de carbono fijo denotada por YCes no nula tendremos que tener en cuenta otra característica de la atmósfera gaseosa: la cantidad de oxígeno. Esto nos permite diferenciar qué sistema de ecuaciones no lineales y ecuaciones diferenciales ordinarias tenemos que resolver, dado que en cada caso estaremos en un dominio distinto de la caldera. •Si la partícula se encuentra en ΩFentonces hay ausencia de oxígeno en un entorno de la partícula y únicamente las reacciones (1.1) y (1.3) intervienen en la gasificación. Por tanto, tenemos que resolver el siguiente sistema: 11 3 λ1 λ=Yg CO2+11 3 λ1 λexp λD De1−a rc−λ 3 2 λ3 λ−λ5 λ=Yg H2O+3 2 λ3 λ−λ5 λexp λD De1−a rc−λ ρ0 C ρgaDr2 c drc dt =−(λ1+λ3)                (1.19) Mientras, la temperatura de la partícula viene dada por: 4 3πa3ρpcs dTp dt = 4πa2(q00 p+q00 r)+4πρgaD(q1λ1+q3λ3+q4λ4+q5λ5),(1.20) donde q00 pyq00 rvienen dados por (1.18) y (1.12), respectivamente. •Si la partícula está en ΩOentonces las reacciones (1.1), (1.2) y (1.3) contribuyen a la gasificación del carbono fijo mediante las especies CO2,O2y 8CAPÍTULO 1. MODELO DE COMBUSTIÓN H2O. El modelo que rige la gasificación de la partícula resulta ser: 11 3 λ1 λ=Yg CO2+11 3 λ1 λexp λD De1−a rc−λ 4 3 λ2 λ=Yg O2+4 3 λ2 λexp λD De1−a rc−λ 3 2 λ3 λ−λ5 λ=Yg H2O+3 2 λ3 λ−λ5 λexp λD De1−a rc−λ ρ0 C ρgaDr2 c drc dt =−(λ1+λ2+λ3)                        (1.21) Por otro lado, la temperatura de la partícula evoluciona de la siguiente manera: 4 3πa3ρpcs dTp dt = 4πa2(q00 p+q00 r) + 4πρgaD(q1λ1+q2λ2+q3λ3+q4λ4+q5λ5), (1.22) donde q00 pyq00 rvienen dados por (1.18) y (1.12), respectivamente. Observación 1.5.En esta tercera etapa, donde se cumple que Tp> Tc> T1, el lector debe tener en cuenta que los procesos de la segunda etapa pueden seguir ocurriendo mientras las cantidades de volátiles y humedad sean distintos de cero. Por último, si la partícula sigue estando dentro de la zona de combustión pero todas las especies se han consumido debido a las reacciones de gasificación, tenemos una nueva situación debida al alto contenido en cenizas, que es una especie inerte. En ese caso, la partícula no cambia de tamaño, manteniendo su radio igual al inicial a; con lo cual podríamos considerar que es una nueva etapa hasta que la partícula abandone la caldera donde ésta sigue intercambiando calor con la atmósfera gaseosa. Su comportamiento se puede describir con las mismas ecuaciones de la primera etapa. 1.2.2. Fuentes homogéneas de la fase gaseosa El modelo de combustión del carbón pulverizado consiste en dos modelos acoplados: uno para la fase gaseosa (el cual no trataremos en nuestro trabajo y determina la atmósfera en la que se encuentra nuestra partícula) y otro para la fase sólida. La resolución del modelo de fase sólida va a permiter obtener las fuentes de masa y energía que las partículas aportan a la fase gaseosa. Las fuentes homogeneizadas en la fase gaseosa por unidad de volumen y tiempo, en un punto x, se calculan mediante las aportaciones de cada partícula que está en la posición x en un instante t. Para ello, utilizaremos la expresión: fα(x) = Ne X j=1 ˜qjZtj f 0 Fα j(t)δ(x−xj p(t))dt, (1.23) 1.2. MODELO DE LA FASE SÓLIDA 9 donde Fα j(t)es la fuente de masa o energía de una partícula introducida por la entrada jen el instante t,xj p(t)es la posición que ocupa la partícula en el instante t,δ(x)es la medida de Dirac en el punto 0, tj fes el tiempo de la partícula para consumirse por completo o para salir del dominio de la caldera, ˜qjes el flujo de masa por la entrada jyNees el número de entradas, respectivamente. En particular, para el modelo simplificado que estamos considerando, las fuentes de masa para cada especie debidas a una partícula vienen dadas por: 1. Si estamos en ΩF, es decir, donde hay ausencia de oxígeno tenemos: Fm O2= 0,(1.24) Fm SO2= 0,(1.25) Fm CO2=4πak cp−11 3λ1,(1.26) Fm H2O=4πak cpλ5−3 2λ3,(1.27) Fm CO =4πak cp14 3λ1+7 3λ3,(1.28) Fm V=4πak cp λ4,(1.29) Fm H2=4πak cp 1 6λ3.(1.30) 2. En cambio, si nos encontrásemos en ΩOtendríamos: Fm O2=4πak cp−4 3λ2,(1.31) Fm CO2=4πak cp−11 3λ1,(1.32) Fm H2O=4πak cpλ5−3 2λ3,(1.33) Fm SO2= 0,(1.34) Fm CO =4πak cp14 3λ1+7 3λ2+7 3λ3,(1.35) Fm V=4πak cp λ4,(1.36) Fm H2=4πak cp 1 6λ3.(1.37) Finalmente, podemos expresar las fuentes de masa total y de energía de la siguiente manera: Fm=4πak cp λ, (1.38) 16 CAPÍTULO 2. MÉTODOS NUMÉRICOS Para estudiar la estabilidad de los métodos Runge-Kutta mencionaremos un resultado de los métodos multipaso que se puede ver en [9]. Definición 2.12. Un método lineal de kpasos se corresponde con la siguiente fórmula general:      k P j=0 αjyi+j=h k P j=0 βjfi+jcon i= 0, . . . , n −k y0, . . . , yk−1dados donde fi=f(ti, yi). Manteniendo esta notación, podemos definir el primer polinomio característico como ρ(r) = k P j=0 αjrj. Teorema 2.13. (Dahlquist) Una condición necesaria y suficiente para que un método multipaso sea estable es que todas las raíces de ρ(r)=0sean de módulo menor o igual que 1 y que las de módulo 1 sean simples. Proposición 2.14. Un método Runge-Kutta es estable ya que aplicando el teorema (2.13) la única raíz de ρ(r) = r+ 1 es −1que es simple y de módulo unidad. Proposición 2.15. Un método Runge-Kutta es consistente si y sólo si s P l=1 bl= 1. Proposición 2.16. Todo método de Runge-Kutta consistente es convergente 1. En este trabajo no abordaremos como establecer el orden de los métodos Runge-Kutta puesto que no es uno de nuestros objetivos, aunque puede consultarse en [6, pág.143-154]. Aún así, enunciaremos una propiedad relevante relacionada con el orden de los métodos Runge-Kutta. Propiedad 2.17. El orden de un método Runge-Kutta explícito de setapas no puede ser mayor que s. Además, no existen métodos Runge-Kutta de setapas explícitos con orden s≥5. A continuación, mostraremos los máximos órdenes alcanzables por los Runge-Kutta explícitos en la siguiente tabla: etapas 1 2 3 4 6 7 9 11 orden 1 2 3 4 5 6 7 8 node condiciones 1 2 4 8 17 37 85 200 1El teorema de Lax-Richtmyer establece la equivalencia entre estabilidad y convergencia para los métodos consistentes de problemas (PVI) bien planteados. 2.2. MÉTODOS DE 1-PASO 17 2.2.3. Métodos Runge-Kutta adaptativos En esta subsección expondremos los métodos Runge-Kutta adaptativos o encajados que consisten en adaptar el número y la posición de los nodos de la discretización de la variable independiente. Su principal objetivo es asegurar que el error local no supere una cierta cota prefijada de antemano. Para llevar a cabo dicho objetivo consideraremos simultáneamente dos métodos Runge- Kutta de setapas con orden pyq2, respectivamente, los cuales emplearán los mismos coeficientes para la condición de la suma por filas. De esta manera, podemos formularlos mediante las siguientes expresiones: ki,l =f xi,l, yi+h s X j=1 aljki,j , l = 1, . . . , s, (2.13) yi+1 =yi+h s X l=1 blki,l,(2.14) ˆyi+1 =yi+h s X l=1 ˆ blki,l.(2.15) De forma equivalente, se puede hacer a través de la siguiente tabla de Butcher: c2a21 c3a31 a32 . . .. . .. . .... csas1as2· · · ass−1 b1b2· · · bs−1bs ˆ b1ˆ b2· · · ˆ bs−1ˆ bs e1e2· · · es−1es donde el método de orden pviene dado por las matrices A,Cy el vector bmientras que el método de orden qse identifica con A,Cy ˆ b. Además, denotaremos por eal vector columna cuyas componentes son el=ˆ bl−blcon l= 1, . . . , s. A continuación, tomando las diferencias entre las soluciones aproximadas obtenidas mediante ambos métodos para un nodo de la discretización, podemos proporcionar una estimación del error local en dicho punto para el esquema de orden p. Esta estimación se 2Tendremos en cuenta que q > p y concretamente en este trabajo consideraremos dos esquemas de métodos Runge-Kutta cumpliendo que q=p+ 1. 18 CAPÍTULO 2. MÉTODOS NUMÉRICOS puede expresar mediante la siguiente fórmula: ξi(h) = ˆyi−yi=h s X l=1 elki,l con i= 1, . . . , n, (2.16) aprovechándonos de que los coeficientes ki,l coinciden en ambas formulaciones. Por otra parte, hemos supuesto que las aproximaciones obtenidas a partir de los métodos Runge-Kutta son de orden pyp+ 1, con lo cual podemos escribir las siguientes expresiones: yi=y(ti) + O(hp+1), ˆyi=y(ti) + O(hp+2), con i= 1, . . . , n, lo que conlleva que la diferencia de las dos aproximaciones cumpla la siguiente condición: ˆyi−yi=y(ti)−yi+O(hp+2)con i= 1, . . . , n. A partir de esta expresión se puede concluir que cuanto menor sea el paso hmayor será la precisión de la estimación local del error en el esquema de orden p. A continuación concretaremos los coeficientes de las tablas de Butcher para los métodos Runge-Kutta que utilizaremos en la resolución numérica del modelo matemático que estamos estudiando. Los esquemas que se muestran en los Cuadros 2.1 y 2.2 se deben a Fehlberg y utilizan simultáneamente un método de orden 4 emparejado con otro de orden 5. 0 2 9 2 9 1 3 1 12 1 4 3 4 69 128 −243 128 135 64 1−17 12 27 4−27 5 16 15 5 6 65 432 −5 16 13 16 4 27 5 144 bt1 909 20 16 45 1 12 0 ˆ bt47 450 012 25 32 225 1 30 6 25 et−1 150 03 100 −16 75 −1 20 6 25 Cuadro 2.1: Esquema RKF4(5). 2.2. MÉTODOS DE 1-PASO 19 0 1 4 1 4 3 8 3 32 9 32 12 13 1932 2197 −7200 2197 7296 2197 1439 216 −83680 513 −845 4104 1 2−8 27 2−3544 2565 1859 4104 −11 40 bt25 216 01408 2565 2197 4104 −1 50 ˆ bt16 135 06656 12825 28561 56430 −9 50 2 55 et1 360 0−128 4275 −2197 75240 1 50 2 55 Cuadro 2.2: Esquema RKF45. Observación 2.18.El esquema RKF45 tiende a subestimar el error local en el método de orden p. Como tal, su uso no es completamente fiable cuando el tamaño del paso hes grande como se puede ver en [9]. Por último, haremos referencia a una característica determinante de los Runge-Kutta adaptativos que es el control del tamaño del paso h. Para cualquier paso inicial hel esquema calculará dos soluciones aproximadas yieˆyiy después determinará una estimación del error local en el nodo ti. Esta estimación queremos que verifique que: |ˆyi,l −yi,l| ≤ m´ax(|yi−1,l|,|yi,l|)τ=sclcon l= 1, . . . , s ei= 1, . . . , n, donde τes la tolerancia prefijada por el usuario para el error relativo. Para medir el error emplearemos la siguiente expresión: eh,i =v u u t 1 s s X l=1 ˆyl,i −yl,i scl2 , que estudiaremos si se aproxima a 1 para encontrar el paso óptimo. De esta manera, el paso óptimo sigue la siguiente fórmula: hopt =hant 1 eh,i  1 q+ 1 .(2.17) Estas expresiones son correctas teóricamente, pero a la hora de implementar el código en el ordenador debemos de ser más precisos. Por ello, no permitiremos un crecimiento o decrecimiento drástico del paso hmediante la utilización de la siguiente expresión: h=hant m´ın   ˆγ, m´ax   ˜γ, 0.91 eh,i  1 q+ 1     ,(2.18) 20 CAPÍTULO 2. MÉTODOS NUMÉRICOS donde ˆγy˜γson dos factores que permiten asegurar la condición anterior. 2.3. Método de Newton En esta sección explicaremos el método de Newton, el cual utilizaremos para resolver los sistemas de ecuaciones no lineales que aparecen en (1.19) y (1.21). El método de Newton es uno de los métodos iterativos más conocidos para la búsqueda de las raíces de una función, debido a su orden de convergencia y su bajo coste computacional. 2.3.1. Método de Newton para funciones de una variable Dada f: [a, b]⊂R→R, trataremos de encontrar un punto fijo de la función: g(x) = x−f(x) f0(x).(2.19) Para ello plantearemos un algoritmo iterativo, que construirá una sucesión de puntos de la siguiente forma:    x0dado, xi+1 =xi−f(xi) f0(xi)con i≥0.(2.20) Observación 2.19.Nótese que la función fdebe ser derivable y su derivada no puede anularse en ninguno de los elementos de la sucesión. Para su implementación en el ordenador debemos construir un test de parada, ya que es necesario un criterio que pueda asegurar que se ha alcanzado una buena aproximación de la solución. Dada una tolerancia, que denotaremos por ¯τ, aplicaremos el algoritmo para construir la sucesión de elementos hasta que se cumpla alguna de las siguientes condiciones: Dos iterantes consecutivos están muy próximos si: |xi+1 −xi|<¯τcon i= 0, . . . , n −1. La distancia relativa entre dos iterantes es cercana a cero, es decir: |xi+1 −xi| |xi+1|<¯τcon i= 0, . . . , n −1. El valor de la función en el iterante es próximo a cero, lo que expresaremos mediante: |f(xi)|<¯τcon i= 0, . . . , n. 2.3. MÉTODO DE NEWTON 21 En particular, para este trabajo usaremos la tercera opción siguiendo el criterio escogido en [5]. A continuación enunciaremos una serie de resultados sobre la convergencia global y local del método de Newton, los cuales pueden verse con más detenimiento en [10]. Teorema 2.20. Sea f: [a, b]⊂R→R,f∈C2([a, b]) cumpliendo: a) f(a)f(b)<0, b) f0(x)6= 0,∀x∈[a, b], c) f00(x)≤0of00(x)≥0,∀x∈[a, b], d) Si c∈ {a, b}denota el extremo de [a, b]en el que |f0|es más pequeño, se tiene que  f(c) f0(c) ≤b−a. Entonces, la ecuación f(x)=0 tiene una única raíz α∈(a, b)y para cualquiera que sea el iterante inicial x0∈[a, b]el método de Newton-Rahpson converge a α. Teorema 2.21. Sea f: [a, b]⊂R→Rcon una raíz α∈(a, b). Si fes derivable en un entorno de αde la forma (α−δ1, α +δ1)⊂[a, b],f0es continua en αyf0(α)6= 0. Entonces: a) La aplicación g(x) = x−f(x) f0(x)está definida en un entorno (α−δ2, α +δ2)tal que δ2< δ1. b) ges derivable en αyg0(α)=0. c) Existe un entorno (α−δ3, α +δ3)tal que para todo x0∈(α−δ3, α +δ3)la sucesión obtenida aplicando el algoritmo converge a αy además xi∈(α−δ3, α +δ3)para todo i≥0. d) Para todo x0∈(α−δ3, α +δ3)la convergencia es superlineal: l´ım i→∞ |xi+1 −α| |xi−α|= 0. Teorema 2.22. (Estimación asintótica del error) Sea αuna raíz simple de f: [a, b]⊂R→Rtal que f0(α)6= 0 y supongamos que existe un entorno (α−δ, α+δ)donde fes dos veces continuamente derivable. Entonces si x06=α y se verifica que xi6=αpara i≥0definiendo ei=xi−αse verifica que l´ım i→∞ ei+1 e2 i =f00(α) 2f0(α). de modo que el método es al menos orden 2. 22 CAPÍTULO 2. MÉTODOS NUMÉRICOS Observación 2.23.Nótese que el algoritmo depende de hacer la evaluación de la derivada de la función fy que si el iterante inicial no es suficientemente cercano a la solución exacta el método podría no converger. 2.3.2. Método de Newton para funciones de varias variables En esta subsección trataremos de generalizar el planteamiento del método de Newton para funciones de varias variables, donde hemos seguido [3], puesto que en este trabajo los sistemas no lineales (1.19) y (1.21) son de mecuaciones con mincógnitas. Sea F:R ⊂ Rm→Rm, con m≥1, una función vectorial real de varias variables donde Res un abierto de Rm. Queremos resolver el sistema de ecuaciones planteado por: F(x) = 0⇔             F1(x1, . . . , xm)=0 F2(x1, . . . , xm)=0 . . . Fm(x1, . . . , xm) = 0 Para resolverlo construiremos una función G:R ⊂ Rm→Rmdada por: G(x) = x−J(x)−1F(x),(2.21) donde J(x)denota la matriz jacobiana asociada a la función vectorial F, cuyas componentes las podemos escribir de la siguiente manera: J(x)ij =∂Fi(x) ∂xj =∂jFi(x). Siguiendo el proceso del caso unidimensional construiremos un algoritmo iterativo que podemos expresar como: (x0∈ R dado, xi+1 =xi−J(xi)−1F(xi) = G(xi)con i≥0,(2.22) junto con un criterio de parada que nos asegure que alcanzamos una buena aproximación de la solución exacta. Para finalizar, enunciaremos un resultado sobre la convergencia del método de Newton para funciones vectoriales reales de varias variables. Teorema 2.24. Supongamos que αes solución de la ecuación x=F(x). Si existe un número δque cumpla que: a) ∂Fi ∂xj es continua en un entorno de αdado por Nδ={x∈ R/||x−α|| < δ}, para todo i= 1, . . . , m yj= 1, . . . , m. 2.3. MÉTODO DE NEWTON 23 b) ∂2Fi(x) ∂xj∂xk es continua, y  ∂2Fi(x) ∂xj∂xk ≤Mpara alguna constante M y para cualquier x∈Nδ con i, j, k = 1, . . . , m. c) ∂Fi(α) ∂xj = 0 para todo i= 1, . . . , m yj= 1, . . . , m. Entonces existe un número ˜ δ < δ tal que la sucesión generada por el algoritmo del método de Newton converge de forma cuadrática a αpara cualquier elección del iterante inicial cumpliendo que ||x0−α|| < δ. Además, ||xi−α||∞≤m2M 2≤ ||xi−1−α||2 ∞,para cada i≥1. 2.3.3. Método de Newton discretizado para funciones de varias variables Por último, tenemos que resaltar que el método de Newton presenta dificultades para su implementación, debido a que debemos evaluar la matriz jacobiana de la función F. Por ello en la práctica utilizaremos el método de Newton discretizado. El método de Newton discretizado consiste en aproximar las derivadas de la función F mediante cocientes incrementales como la siguiente aproximación: Jjk(xi, Hi) = Fj(xi+hi jkek)−Fj(xi) hi jk ≈∂kFj(xi),(2.23) donde hi jk son parámetros de discretización conocidos y ekes el k-ésimo vector de la base canónica de Rm. Además, {Hi}i≥0es una sucesión de matrices que se aproximan a la matriz nula. Por tanto, podemos reescribir el sistema de la siguiente forma: (x0∈ R dado, xi+1 =xi−Jjk(xi, Hi)−1F(xi) = G(xi)con i≥0,(2.24) junto con el test de parada correspondiente que nos asegure que alcanzamos una buena aproximación de la solución exacta. Finalmente, siguiendo la estructura de las secciones anteriores, enunciaremos un resultado sobre la convergencia del método de Newton discretizado siguiendo la referencia [9]. Teorema 2.25. Sea F:R ⊂ Rm→Rmuna función C1yRun abierto convexo tal que α∈ R. Si existen dos constantes positivas yhcumpliendo que: a) x0∈B(α, ) = {x∈Rm/||x−α||1< }, b) 0<|hi jk|< h para todo j, k = 0, . . . , m 24 CAPÍTULO 2. MÉTODOS NUMÉRICOS Entonces la sucesión dada por (2.24) esta bien definida y converge de manera lineal a α. Además, si existe una constante positiva Ctal que m´ax j,k=1,...,n |hi jk| ≤ C||xi−α||1, entonces la sucesión converge de forma cuadrática. Capítulo 3 Resultados numéricos En este capítulo expondremos los resultados obtenidos utilizando los diferentes códigos implementados a lo largo del trabajo sobre un ejemplo test. Primero presentaremos el problema test que queremos resolver y, posteriormente, los resultados a los cuales hemos llegado mediante la resolución del mismo, usando los diferentes métodos descritos en el capítulo anterior. 3.1. Problema test El problema test se basa en considerar el problema dado por las ecuaciones (1.13), (1.14), (1.21) y (1.22), y añadir unos términos fuente, de forma que la solución exacta del problema sea conocida. De esta manera podremos calcular los errores que cometen los códigos que hemos implementado y verificar si su comportamiento es correcto. En este test hemos considerado que el dominio temporal es [0,2], y tomaremos el siguiente problema describiendo el comportamiento de una partícula de carbón pulverizado sufriendo los procesos descritos en el capítulo 1 para la tercera etapa. El test seleccionado cumplirá que, en el intervalo de tiempo mencionado anteriormente, las cantidades de humedad, volátiles y carbono fijo serán no nulas, la fracción másica de oxígeno en la atmósfera será distinta de cero y la temperatura de la partícula será mayor que T1, para que ocurran las reacciones de volatilización (1.13) y (1.14), y también mayor que Tcpara que ocurran las reacciones de gasificación del carbono fijo dadas por el sistema (1.21). Con ello, el modelo matemático que hay que resolver resulta: dρV dt =−B4exp −E4 RTpρV+gρV(t),(3.1) 25 32 CAPÍTULO 3. RESULTADOS NUMÉRICOS 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 300 310 320 330 340 350 360 370 380 390 400 Temperatura Sol. exacta Sol. aproximada Figura 3.6: Aproximación de Tpdel esquema 1. En el Cuadro 3.3 debemos destacar que todas las variables presentan reducción del error a medida que aumentamos el número de nodos. Además, en las tres primeras columnas correspondientes a la velocidades de gasificación, se observa que el método de Newton presenta un orden ligeramente superior a 1, lo cual estropea el orden 4 del Runge-Kutta aplicado al radio del núcleo de carbono fijo. En cambio, esto no sucede al resolver la temperatura. Nodos Errorλ1Errorλ2Errorλ3ErrorrcErrorTp 25 0.010753 0.025541 0.014378 8.1464·10−60.0137 50 0.0052437 0.012558 0.0068502 3.6586·10−60.00096819 100 0.0025627 0.0061524 0.0033447 1.7463·10−61.5733·10−5 250 0.0010128 0.0024314 0.0013197 6.8011·10−71.5733·10−6 500 0.00050434 0.0012108 0.00065675 3.3719·10−79.7661·10−8 1000 0.00025163 0.00060414 0.00032761 1.6788·10−76.0782·10−9 Cuadro 3.3: Errores del esquema 1. 3.2. RESULTADOS 33 3.2.1.2. Esquema 2: Sistema acoplado, Euler explícito Figura 3.7: Diagrama de flujo para los esquemas 2, 3 y 4. Los tres próximos esquemas se basan en semidiscretizar la derivada temporal asociada a la variable rcdel sistema (3.3) mediante los métodos vistos en el capítulo 2 y resolver el sistema no lineal obtenido mediante el método de Newton. Este proceso se ha sintetizado de manera gráfica en la Figura 3.7. En el esquema 2 semidiscretizamos la variación temporal de la variable rcmediante el método de Euler explícito. En las Figuras 3.8, 3.9 y 3.10, podemos observar las aproximaciones de las distintas variables, en las cuales podemos destacar que cualitativamente la resolución de las velocidades parece empeorar. Por otra parte, en la solución aproximada del radio de carbono fijo observamos la acumulación del error, un hecho distintivo de los métodos de Euler. 34 CAPÍTULO 3. RESULTADOS NUMÉRICOS 012 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1Velocidad lambda1 012 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3Velocidad lambda2 012 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 3.2Velocidad lambda3 Sol. exacta Sol. aproximada Figura 3.8: Aproximaciones de λ1,λ2yλ3del esquema 2. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 6 7 8 9 10 11 12 10-5 Radio carbono fijo Sol. exacta Sol. aproximada Figura 3.9: Aproximación de rcdel esquema 2. 3.2. RESULTADOS 35 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 300 310 320 330 340 350 360 370 380 390 400 Temperatura Sol. exacta Sol. aproximada Figura 3.10: Aproximación de Tpdel esquema 2. En el Cuadro 3.4 mostraremos los errores de la evolución a lo largo del tiempo de las distintas soluciones aproximadas donde debemos reseñar que son peores cuantitativamente para las velocidades, y por otra parte, mejoran para el radio del núcleo de carbono fijo y la temperatura. En todos los casos se sigue reduciendo el error al aumentar el número de nodos, mientras observamos orden algo mejor que 1 en las variables del sistema y cercano a orden 4 en la resolución de la temperatura. Nodos Errorλ1Errorλ2Errorλ3ErrorrcErrorTp 25 0.1277 0.16564 0.46694 2.8414·10−60.012064 50 0.065821 0.081378 0.23196 1.387·10−60.00091421 100 0.033371 0.040334 0.11496 6.8557·10−76.0185·10−5 250 0.013454 0.016049 0.045748 2.7233·10−71.555·10−6 500 0.0067445 0.0080106 0.022832 1.3585·10−79.7115·10−8 1000 0.0033766 0.0040018 0.011405 6.7847·10−86.0609·10−9 Cuadro 3.4: Errores del esquema 2. Por otra parte, aunque el esquema 1 presente ciertas oscilaciones en el radio de carbono fijo, de manera global parece cometer menos error que el esquema 2 donde semidiscretizamos mediante Euler explícito la ecuación diferencial asociada a la variable rc. 36 CAPÍTULO 3. RESULTADOS NUMÉRICOS 3.2.1.3. Esquema 3: Sistema acoplado, Euler implícito En el esquema 3 semidiscretizamos la derivada temporal de la variable rcmediante el método de Euler implícito. En las Figuras 3.11, 3.12 y 3.13 mostraremos la evolución a lo largo del tiempo de la solución exacta y aproximada de las variables λ1,λ2,λ3,rcyTp. 012 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1Velocidad lambda1 012 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3Velocidad lambda2 012 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 3.2Velocidad lambda3 Sol. exacta Sol. aproximada Figura 3.11: Aproximaciones de λ1,λ2yλ3del esquema 3. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 6 7 8 9 10 11 12 10-5 Radio carbono fijo Sol. exacta Sol. aproximada Figura 3.12: Aproximación de rcdel esquema 3. 3.2. RESULTADOS 37 En ellas podemos observar como se ajustan mejor que en el esquema 1, aunque ya no percibimos el comportamiento de acumular el error en el radio de carbono fijo sino una oscilación más suave como en el esquema 1. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 300 310 320 330 340 350 360 370 380 390 400 Temperatura Sol. exacta Sol. aproximada Figura 3.13: Aproximación de Tpdel esquema 3. En el Cuadro 3.5 se muestran los errores asociados a cada variable por columnas; donde debemos reseñar la reducción del error en las velocidades y el radio por el nuevo método para discretizar la derivada. En él seguimos observando el descenso del error al aumentar el número de nodos y un orden mayor que 1 en las variables asociadas al sistema (3.3). Nodos Errorλ1Errorλ2Errorλ3ErrorrcErrorTp 25 0.0020461 0.006094 0.0042057 1.5841·10−60.013651 50 0.0010288 0.0030689 0.0021252 7.9948·10−70.00096705 100 0.00051543 0.0015397 0.0010657 4.0112·10−76.1954·10−5 250 0.00020644 0.00061666 0.00042714 1.6083·10−71.5731·10−6 500 0.00010326 0.00030843 0.00021371 8.0471·10−89.7653·10−8 1000 5.1614·10−50.00015424 0.00010688 4.0249·10−86.078·10−9 Cuadro 3.5: Errores del esquema 3. 38 CAPÍTULO 3. RESULTADOS NUMÉRICOS 3.2.1.4. Esquema 4: Sistema acoplado, RK4 En el esquema 4 semidiscretizamos la variación respecto al tiempo de la variable rc mediante el método de Runge-Kutta clásico cuya tabla de Butcher se puede ver en la Figura 3.1. En las Figuras 3.14, 3.15 y 3.16 mostramos la solución exacta y aproximada para las velocidades de gasificación, el radio de carbono fijo y la temperatura de la partícula. 012 0.5 0.6 0.7 0.8 0.9 1 1.1Velocidad lambda1 012 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 3.2Velocidad lambda2 012 1.2 1.4 1.6 1.8 2 2.2 2.4 2.6 2.8 3 3.2Velocidad lambda3 Sol. exacta Sol. aproximada Figura 3.14: Aproximaciones de λ1,λ2yλ3del esquema 4. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 6 7 8 9 10 11 12 10-5 Radio carbono fijo Sol. exacta Sol. aproximada Figura 3.15: Aproximación de rcdel esquema 4. 3.2. RESULTADOS 39 De manera cualitativa podemos asegurar que mediante el esquema 4 obtenemos las mejores aproximaciones entre los esquemas de paso constante, por la manera que se ajustan las curvas. 0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 300 310 320 330 340 350 360 370 380 390 400 Temperatura Sol. exacta Sol. aproximada Figura 3.16: Aproximación de Tpdel esquema 4. En el Cuadro 3.6 mostramos los errores cometidos en cada aproximación de las respectivas variables dependiendo del número de nodos. A partir de él podemos asegurar que el esquema 4 es el que menor error comete, viendo como se reduce a medida que aumentamos el número de nodos y observando orden cercano a 2 en las variables resueltas mediante Newton discretizado y orden 4 para la resolución de la temperatura. Nodos Errorλ1Errorλ2Errorλ3ErrorrcErrorTp 25 0.00070463 0.0018073 0.0010677 6.2088·10−70.013668 50 0.0001683 0.00043367 0.00025573 1.503·10−70.00096743 100 4.1358·10−50.00010637 6.2676·10−53.6875·10−86.1966·10−5 250 6.5552·10−61.6862·10−59.9322·10−65.8467·10−91.5732·10−6 500 1.6419·10−64.2225·10−62.4858·10−61.463·10−99.7655·10−8 1000 4.1454·10−71.0657·10−66.2679·10−73.6879·10−10 6.078·10−9 Cuadro 3.6: Errores del esquema 4. 40 CAPÍTULO 3. RESULTADOS NUMÉRICOS 3.2.2. Esquemas de paso adaptativo Aunque los métodos de paso fijo vistos anteriormente aproximan satisfactoriamente nuestro ejemplo test, presentan una gran desventaja al intentar resolver un problema real que modele la combustión de carbón pulverizado. Esta desventaja radica en cómo son las soluciones exactas para nuestras variables de interés, ya que experimentan gradientes muy fuertes, por lo que con la utilización de métodos de paso fijo, o bien perderíamos información, o bien tendríamos que utilizar pasos de tiempo demasiado pequeños en intervalos de tiempo donde no sería necesario. Para atajar este problema consideramos métodos de paso variable en los cuales intentamos ajustar el paso de tiempo, y de esta manera explicar mejor las zonas donde se dan cambios bruscos en el gradiente. Por ello en esta subsección explicaremos dos esquemas que difieren en la metodología utilizada para resolver el sistema (3.3). En ellos hemos implementado un método de paso adaptativo como es el Runge-Kutta-Fehlberg (RKF45) expuesto en el capítulo 2, cuya tabla de Butcher se puede ver en el Cuadro 2.2. Como en la sección anterior, dicho método se utilizará para resolver las ODE’s que intervienen en el problema test expuesto. Figura 3.17: Diagrama de flujo de los esquemas 5 y 6. 3.2. RESULTADOS 41 Para la implementación de la selección de paso hemos utilizado la referencia [3], puesto que es más sencilla que la expuesta en el capítulo 2 y requiere menos cálculos. Su expresión es la siguiente: h=hant m´ın     ˆγ, m´ax     ˜γ, 0.84    τ |1 hant ξi+1(hant)|   1/4       ,(3.5) con i= 0, . . . , n −1. Nótese que en (3.5) hemos mantenido la notación para las variables τ,ˆγ,˜γyξi+1(hant). Además, en la Figura 3.17 podemos ver el diagrama de flujo, dando una idea de la estructura del código implementado. 3.2.2.1. Esquema 5: Sistema desacoplado, RKF45 El esquema 5 se basa en resolver de manera “desacoplada” el sistema (3.3). El procedimiento utilizado es análogo al propuesto en el esquema 1, y ambos esquemas únicamente se diferenciarán en la resolución de las ODE’s. Mientras que en el esquema 1 se resolvía el comportamiento de las variables ρV,ρH2O,rcyTpmediante el método RK4, en este nuevo esquema usaremos el método RKF45. A continuación, en las Figuras 3.18, 3.19, 3.20, 3.21 y 3.22 se muestra la evolución en el tiempo de las variables que intervienen en nuestro modelo y podremos valorar de forma cualitativa las soluciones exactas y aproximadas para nuestro problema test. 0 0.5 1 1.5 2 0 10 20 30 40 50 60 70 80 90 100 Densidad de volatiles Sol. exacta Sol. aproximada 0 0.5 1 1.5 2 0 1 2 3 4 5 610-4Velocidad lambda4 Figura 3.18: Aproximaciones de ρVyλ4del esquema 5. Capítulo 4 Conclusiones En este trabajo de Fin de Grado hemos realizado un estudio sobre la resolución de modelos de combustión. En primer lugar, hemos descrito el modelo de combustión BFLs1 presentado en [2], centrándonos en el desarrollo de su fase sólida. Posteriormente, hemos hecho una revisión de los métodos para la resolución de ecuaciones diferenciales y sistemas de ecuaciones no lineales, puesto que a partir de ellos podemos aproximar el comportamiento de los distintos procesos que suceden en la combustión de las partículas de carbón pulverizado. A continuación, utilizando estos métodos hemos implementado seis códigos diferentes para modelar la combustión del carbón, y mediante la resolución de un test, hemos podido validar su correcto funcionamiento. Además, por los resultados vistos en el capítulo 3 los esquemas que resuelven el sistema de DAE’s de forma conjunta aproximan mejor el modelo, debido a que ajustan de manera más adecuada el radio de carbono fijo. De esta manera podemos asegurar que, tanto para los métodos de paso fijo como para los de paso adaptativo, la semidiscretización mediante el método RK4 de la derivada temporal del radio es el que mejor aproxima las soluciones exactas para las distintas variables. Además, hemos comprobado que los métodos que adaptan el paso por el comportamiento de la solución resuelven de manera más eficiente las ecuaciones que intervienen en el modelo de combustión. Por último, como trabajo futuro podríamos desarrollar los códigos que hemos implementado para tener en cuenta la trayectoria de cada partícula de carbón pulverizado dentro de la caldera debido a que es un modelo tridimensional como vimos en el capítulo 1. Además, si implementásemos códigos de localización y supiésemos como es la atmósfera gaseosa podríamos ser capaces de afrontar un problema real. 49 Nomenclatura A continuación se especificará la notación usada a lo largo del trabajo de Fin de Grado. 1. Símbolos relativos a los gases. Notación Descripción Unidades (SI) cpcalor específico a presión constante J/kgK Dcoeficiente de difusión de la mezcla de gases m2/s kconductividad térmica W/mK Rconstante universal de los gases J/molK σconstante de Stefan-Boltzmann de los gases W/m2K4 ρgdensidad de la atmósfera kg/m3 hTentalpía térmica J/kg Yifracción másica de la especie i-ésima fefuente de energía homogeneizada en la fase gaseosa J/m3s fmfuente de masa homogeneizada procedente de las partículas kg/m3s ggravedad m/s2 Mmasa molecular kg Tgtemperatura de la atmósfera K vvelocidad m/s µviscosidad molecular kg/ms 2. Símbolos relativos a las reacciones químicas. Notación Descripción Unidades (SI) qicalor de la reacción i-ésima J/kg Eienergía de activación de la reacción i-ésima J/mol Bifactor de frecuencia de la reacción i-ésima 1/s λivelocidad adimensional de la reacción i-ésima 51 52 NOMENCLATURA 3. Símbolos relativos a la partícula. Notación Descripción Unidades (SI) cscalor específico J/kgK Decoeficiente de difusión efectiva a través de los poros m2/s ρpdensidad kg/m3 ρidensidad de la especie i kg/m3 εpcoeficiente de emisividad q00 pflujo de calor por conducción J/m2s q00 pflujo de calor por radiación J/m2s Fefuente de energía procedente de la partícula J/m3s Fm ifuente de masa de la especie i-ésima procedente de las partículas kg/m3s mpmasa kg xpposición aradio inicial m rcradio del núcleo de carbono fijo m T1temperatura de volatilización K Tctemperatura de gasificación K Tptemperatura K vpvelocidad m/s Apéndice A Datos del problema test A continuación se especificarán los datos utilizados a lo largo del trabajo de Fin de Grado. 1. Datos relativos a los gases. Notación Valor Unidades (SI) cp1000 J/kgK D2.88·10−15 m2/s k0.0454 W/mK R8.314472·103J/molK σ5.67·10−8W/m2K4 ρg0.25 kg/m3 YCO20.4 YO20.4 YH2O0.2 YSO20 rad 0 W/m2 Tg1750 K 53 54 APÉNDICE A. DATOS DEL PROBLEMA TEST 2. Datos relativos a las reacciones químicas. Notación Valor Unidades (SI) q1-1.43568·107J/kg q29.20247·106J/kg q3-1.09305·107J/kg q40 J/kg q50 J/kg E11.24·106J/mol E22·104J/mol E32·104J/mol E42021 J/mol E58.32·104J/mol B12.46·1081/s B21.367·1081/s B31.37·1081/s B43.11·1081/s B53.228·1081/s 3. Datos relativos a la partícula. Notación Valor Unidades (SI) cs1000 J/kgK De2.88·10−5m2/s ρp900 kg/m3 ρc579 kg/m3 ρa167 kg/m3 εp0.9 a60·10−6m T1300 K Tc300 K Bibliografía [1] Ascher, U.M.; Petzold, L.R., Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations, SIAM cop., Filadelfia, 1998. [2] Bermúdez, A.; Ferrín, J.L.; Liñán, A. y Saavedra, L., Numerical simulation of group combustion of pulverized coal, Combust. and Flame 158(9)(2011), 1852-1865. [3] Burden, R.L. y Douglas, J., Numerical Analysis, 3aed., PWS Publishers, Boston, 1985. [4] Ferrín, J.L., Algunas contribuciones a la modelización matemática de procesos de combustión de carbón(Tesis doctoral), USC, 1999. [5] Ferrín, J.L. y Martínez, S., Resolución numérica de DAE’s en modelos de combustión(Trabajo de Fin de Grado), USC, 2018. [6] Hairier, E.; Nørset, S.P. y Wanner, G., Solving Ordinary Differential Equations I: Nonstiff Problems, 2aed. revisada, Springer series in Computational Mathematics, 8, Springer-Verlag, Berlin, 1993. [7] Hairier, E. y Wanner, G., Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2aed. revisada, Springer series in Computational Mathematics, 14, Springer-Verlag, Berlin, 1996. [8] Lambert, J.D., Numerical methods for ordinary differential systems: the initial value problem, John Wiley and Sons, Chichester, 1991. [9] Quarteroni, A.; Sacco, R. y Saleri, F., Numerical Mathematics, Text in Applied Mathematics 37, Springer-Verlag, Nueva York, 2000. [10] Viaño, J.M., Lecciones de métodos numéricos 2: Resolución de ecuaciones numéricas, Tórculo edicións, Santiago de Compostela, 1997. 55