Simulación numérica de la combustión de carbón pulverizado

XX Congreso de Ecuaciones Diferenciales y Aplicaciones X Congreso de Matem´ atica Aplicada Sevilla, 24-28 septiembre 2007 (pp. 1–8) Simulaci´on num´erica de la combusti´on de carb´on pulverizado A. Berm´ udez 1, J.L. Ferr´ ın1, A. Li˜ n´ an2, L. Saavedra1 1Dpto.Matem´atica Aplicada, Universidad de Santiago de Compostela, Aptdo. 15782, Santiago de Compostela. E-mails: [email protected], [email protected], [email protected]. 2E.T.S. Ingenieros Aeron´auticos, Universidad Polit´ecnica de Madrid, Aptdo. 28040, Madrid. E-mail: [email protected]. Palabras clave: carb´on pulverizado, elementos finitos Resumen El objetivo de esta comunicaci´on es presentar un modelo matem´atico de la combusti´on de carb´on pulverizado y el correspondiente algoritmo para su resoluci´on. El modelo que se ha desarrollado consta de dos fases fuertemente acopladas: la fase gaseosa, para la que se utilizar´a una descripci´on Euleriana, que proporciona las distribuciones de temperatura y de las concentraciones de las distintas especies en el gas y la fase s´olida, para la que se utiliza una descripci´on Lagrangiana, que trata los procesos simult´aneos de evaporaci´on de la humedad y la devolatilizaci´on, junto con las reacciones heterog´eneas de gasificaci´on del “char”. El acoplamiento de estas fases viene determinado por el hecho de que las part´ıculas de carb´on son fuentes y sumideros de masa y energ´ıa, mientras que la fase gaseosa determina el movimiento y la atm´osfera en la que se produce la combusti´on de las part´ıculas. Ese modelo de combusti´on se va a incorporar al programa SC3D (Simulaci´on de Calderas en 3 Dimensiones). Se trata de un programa de CFD (Mec´anica de Fluidos Computacional) adaptado para la simulaci´on de la zona del hogar de una caldera de carb´on pulverizado utilizada en una Central T´ermica. Para la resoluci´on de las ecuaciones en derivadas parciales del modelo de fase gaseosa se utiliza un m´etodo de elementos finitos combinado con el m´etodo de las caracter´ısticas. Los sistemas de ecuaciones diferenciales ordinarias y algebraicas del modelo de fase s´olida se resuelven utilizando los m´etodos de Euler impl´ıcito o expl´ıcito y el m´etodo de Newton discretizado. 1. Modelo de combusti´on Los procesos f´ısico-qu´ımicos a los que est´a sometida una part´ıcula de carb´on durante su combusti´on son de una gran complejidad. Ello obliga a desarrollar modelos simplificados 1 A. Berm´udez, J.L. Ferr´ın, A. Li˜n´an, L. Saavedra para la generaci´on de vol´atiles, la gasificaci´on del “char” y la posterior combusti´on de las especies gaseosas liberadas. El modelo cin´etico simplificado que consideramos consiste en los siguientes procesos f´ısico-qu´ımicos en el interior de las part´ıculas porosas, expresados por las reacciones heterog´eneas: 1CO2+C(s)→2CO + (q1) 21 2O2+C(s)→CO + (q2) 3H2O+C(s)→CO +H2+ (q3) 4V(s)→V(g)+ (q4) 5H2O(s)→H2O(g)+ (q5) y por las siguientes reacciones de oxidaci´on en fase gaseosa: 6CO +1 2O2→CO2+ (q6) 7V(g)+ν1O2→ν2CO2+ν3H2O+ν4SO2+ (q7) 8H2+1 2O2→H2O+ (q8) donde el ´ındice sdenota la fase s´olida y gla fase gaseosa, mientras que qies el calor liberado en la reacci´on ipor unidad de masa gasificada. Por simplicidad, se considerar´an los vol´atiles como mol´eculas ´unicas, V(g)=Cκ1Hκ2Oκ3Sκ4,de masa molecular Mvol,donde los coeficientes κ1, κ2, κ3yκ4se obtienen del an´alisis elemental del carb´on. Las tasas de gasificaci´on, wi, i = 1, ..., 5 para las reacciones heterog´eneas se modelan usando leyes de Arrhenius globales y la tasa de gasificaci´on del char, wC, por unidad de volumen, wC=w1+w2+w3, queda determinada por las velocidades globales de las reacciones 1, 2 y 3. Para la derivaci´on de este modelo se ha generalizado el an´alisis de Burke-Schumann para tener en cuenta la competici´on por el ox´ıgeno entre las especies combustibles. As´ı, se supone que las reacciones 6-8 o no ocurren o lo hacen infinitamente r´apidas en una llama delgada (llama de difusi´on) que tiene lugar en el interior de la part´ıcula, en el gas en la vecindad de la part´ıcula, o en el gas lejos de la part´ıcula. El tipo de combusti´on que ocurra depende de la temperatura y de las concentraciones locales de O2,CO, vol´atiles y H2en el entorno gaseoso y ser´a el propio modelo el que decidir´a en que situaci´on nos encontramos. Entonces el modelo da la posibilidad de que existan zonas en las que se acabe el ox´ıgeno, que denotaremos por ΩF, a diferencia de las zonas donde s´ı hay ox´ıgeno que denotamos por ΩO. Cuando la part´ıcula llega a ΩO, los vol´atiles, CO yH2que se han producido por las reacciones de gasificaci´on se queman completamente, en una llama de difusi´on, dentro de la part´ıcula o fuera en el entorno de la misma. En este caso las part´ıculas no representan fuentes de CO, vol´atiles o H2para la fase gaseosa. Cuando la part´ıcula est´a en ΩF, al no haber ox´ıgeno, las reacciones homog´eneas 6, 7 y 8 no tienen lugar. Entonces los vol´atiles, CO yH2se unen a la fase gaseosa sin quemarse. Como suponemos que las reacciones homog´eneas son infinitamente r´apidas, estas especies se quemar´an cuando encuentren ox´ıgeno en una llama de difusi´on, ondulada por la turbulencia; es la superficie ΓFen la figura 1, que separa la regi´on ΩF, sin ox´ıgeno, de la regi´on ΩOdonde los vol´atiles, H2yCO se encuentran s´olo en el entorno de las part´ıculas de carb´on dentro de la llama delgada que las rodea. 2 Simulaci´on num´erica de la combusti´on de carb´on pulverizado gas/carb´on ΩF ΩO aire secundario ΓF ΩO Figura 1: Llamas de difusi´on. 1.1. Modelo para la fase gaseosa El modelo de combusti´on de part´ıculas de carb´on que presentamos est´a acoplado a una fase gaseosa de dos formas. Por un lado, el modelo de la fase gaseosa establece la temperatura, velocidad, presi´on y concentraciones que determinan la atm´osfera en la que se mueven y queman las part´ıculas. Por otro lado, las part´ıculas debido a su combusti´on aparecen como fuentes o sumideros de masa y energ´ıa. Sea Lgel operador diferencial definido por Lg(u) = ∂(ρgu) ∂t +∇ ·(ρguvg)−∇·(ρgD∇u), donde Des el coeficiente de difusi´on de la fase gaseosa que, por simplicidad, se considerar´a igual para todas las especies e igual a la difusitividad t´ermica. Entonces, las ecuaciones de conservaci´on de la masa, de la masa de cada una de las especies gaseosas y de la energ´ıa, conocidas las fuentes debidas a las part´ıculas, cuyas expresiones se pueden ver en Berm´udez et al [2], incluyen los t´erminos w6,w7yw8debidos a las reacciones qu´ımicas 6, 7 y 8, que tienen lugar en la fase gaseosa. Al considerar la hip´otesis de Burke-Schuman, lo cual implica la no coexistencia de CO, vol´atiles y H2con O2,se pueden considerar los siguientes escalares conservados o combinaciones lineales de Shvab-Zeldovich: βg 1=Yg O2−4 7Yg CO −32ν1 Mvol Yg V−8Yg H2,(1) βg 2=Yg CO2+11 7Yg CO +44ν2 Mvol Yg V,(2) βg 3=Yg H2O+18ν3 Mvol Yg V+ 9Yg H2,(3) βg 4=Yg SO2+64ν4 Mvol Yg V,(4) Hg=hg T+q6Yg CO +q7Yg V+q8Yg H2,(5) y de esta forma se pueden eliminar esos t´erminos de las reacciones gaseosas, obteni´endose Lg(βg 1) = fm O2−4 7fm CO −32ν1 Mvol fm V−8fm H2,(6) Lg(βg 2) = fm CO2+11 7fm CO +44ν2 Mvol fm V,(7) Lg(βg 3) = fm H2O+18ν3 Mvol fm V+ 9fm H2,(8) Lg(βg 4) = fm SO2+64ν4 Mvol fm V,(9) Lg(Hg) = fe+q6fm CO +q7fm V+q8fm H2−∇·qrg.(10) 3 A. Berm´udez, J.L. Ferr´ın, A. Li˜n´an, L. Saavedra As´ı pues, gracias a la hip´otesis de Burke-Schuman, adem´as de resolver menos ecuaciones y m´as sencillas, el propio modelo dice en qu´e regi´on nos encontramos. As´ı, si Xg 1>0 habr´a ox´ıgeno y nos encontraremos en el dominio ΩOy si Xg 1<0 no habr´a ox´ıgeno y estaremos en la regi´on ΩF.Las reacciones 6, 7 y 8 tienen lugar en una llama infinitamente delgada situada donde Xg 1= 0,con w6,w7yw8actuando como deltas de Dirac. 1.2. Modelo para la combusti´on de una part´ıcula El modelo de gasificaci´on del char que proponemos es v´alido para part´ıculas con un alto contenido en cenizas. Adem´as, suponemos que la fracci´on de cenizas que tienen las part´ıculas de carb´on no se pierde durante la devolatilizaci´on ni la oxidaci´on del char, es decir, no se produce fragmentaci´on durante su combusti´on. De esta forma, consideraremos que cada part´ıcula mantiene constante su radio y su densidad de cenizas, ρash, aunque la densidad de H2O, vol´atiles y char disminuye con el tiempo. Para llevar a cabo un tratamiento m´as simple de los efectos de las reacciones de gasificaci´on del char 1, 2 y 3, tendremos en cuenta que las energ´ıas de activaci´on de estas reacciones son grandes. Los n´umeros de Damk¨ohler, que pueden ser definidos para estas reacciones como Dai= (a2/De)Bie−Ei/RTp, i = 1,2,3,determinar´an la etapa de la combusti´on en la que nos encontramos: Primera etapa. Los n´umeros de Damk¨ohler son menores que 1 y se pueden considerar las reacciones 1, 2 y 3 congeladas, por lo tanto debemos resolver el sistema: dρV dt =−3ρgD a2λ4, λ4=a2 3ρgDB4e−E4/RTpρV, dρH2O dt =−3ρgD a2λ5, λ5=a2 3ρgDB5e−E5/RTpρH2O. (11) Segunda etapa. Los n´umeros de Damk¨ohler son mayores que 1 y se consideran las reacciones 1, 2 y 3 infinitamente r´apidas. En esta etapa se pueden dar varias situaciones que determinan las ecuaciones para λ1yλ3. 1. Cuando la part´ıcula est´a en la regi´on ΩOy adem´as Ys O2es mayor que cero (es decir, rc< rf≤a) debemos resolver eλD De(a/rf−1)+λ=ϕ+ 1, 11 3 λ1 λ−·22 3 λ1 λ+11 3 λ3 λ+44ν2 Mvol λ4 λ¸eλD De(a/rf−a/rc) =½Yg CO2−11 3µλ1 λ+λ3 λ¶−44ν2 Mvol λ4 λ¾eλD De(1−a/rc)−λ, 3 2 λ3 λ−λ5 λ−µ3 2 λ3 λ+18ν3 Mvol λ4 λ¶eλD De(a/rf−a/rc) =½Yg H2O−λ5 λ−18ν3 Mvol λ4 λ¾eλD De(1−a/rc)−λ, ρ0 C ρgaDr2 c drc dt =−(λ1+λ3), (12) 4 Simulaci´on num´erica de la combusti´on de carb´on pulverizado con ϕ=eλa/rf−1. La ecuaci´on para la temperatura en este caso 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 +14 3q6λ1+7 3q6λ3+q7λ4+1 6q8λ3).(13) 2. Si la part´ıcula est´a en ΩOy adem´as Ys O2= 0 (es decir, rf> a) debemos resolver 11 3 λ1 λ=½Yg CO2+11 3 λ1 λ+µ22 3 λ1 λ+11 3 λ3 λ+44ν2 Mvol λ4 λ¶ϕ¾eλD De(1−a/rc)−λ, 3 2 λ3 λ−λ5 λ=½Yg H2O+3 2 λ3 λ−λ5 λ+µ3 2 λ3 λ+18ν3 Mvol λ4 λ¶ϕ¾eλD De(1−a/rc)−λ, ρ0 C ρgaDr2 c drc dt =−(λ1+λ3). (14) Es este caso para obtener la temperatura debemos resolver 4 3πa3ρpcs dTp dt = 4πa2(q00 p+q00 r)+4πρgaD(q1λ1+q3λ3+q4λ4+q5λ5 +14 3q6λ1+7 3q6λ3+q7λ4+1 6q8λ3).(15) 3. Si la part´ıcula est´a en ΩOy el ox´ıgeno llega a la superficie del “core”de la part´ıcula, debido a que las reacciones en fase gaseosa est´an congeladas, 11 3 λ1 λ=½Yg CO2+11 3 λ1 λ¾eλe(1−/rc)−λ, 4 3 λ2 λ=½Yg O2+4 3 λ2 λ¾eλe(1−/rc)−λ, 3 2 λ3 λ−λ5 λ=½Yg H2O+3 2 λ3 λ−λ5 λ¾eλe(1−/rc)−λ, ρ0 C ρg r2 c drc dt =−(λ1+λ2+λ3) (16) Para la temperatura tenemos 4 3πa3ρpcs dTp dt = 4πa2(q00 p+q00 r)+4πρgaD(q1λ1+q2λ2+q3λ3+q4λ4+q5λ5 +14 3q6λ1+7 3q6λ3+q7λ4+1 6q8λ3).(17) 4. Finalmente, cuando la part´ıcula est´a en la regi´on ΩF, debemos resolver 11 3 λ1 λ=½Yg CO2+11 3 λ1 λ¾eλD De(1−a/rc)−λ, 3 2 λ3 λ−λ5 λ=½Yg H2O+3 2 λ3 λ−λ5 λ¾eλD De(1−a/rc)−λ, ρ0 C ρgaDr2 c drc dt =−(λ1+λ3) (18) 5 A. Berm´udez, J.L. Ferr´ın, A. Li˜n´an, L. Saavedra junto con la ecuaci´on 4 3πa3ρpcs dTp dt = 4πa2(q00 p+q00 r)+4πρgaD(q1λ1+q3λ3+q4λ4+q5λ5).(19) Por ´ultimo, las fuentes homogeneizadas en la fase gaseosa por unidad de volumen y tiempo, en un punto xde la caldera, se obtienen a partir de las fuentes individuales de cada part´ıcula mediante la expresi´on fα(x) = Ne X j=1 Np X i=1 ˜qj pij 100 Ztij f 0 Fα ij(t)δ(x−xij p(t))dt (20) donde Fα ij(t) es la contribuci´on de una part´ıcula de tipo iintroducida por la entrada j, en el instante t(cuyas expresiones pueden verse en Berm´udez et al [2]), xij s(t) denota la posici´on ocupada por esa part´ıcula en el instante t,δ(x) es la medida de Dirac en el punto 0, tij fes el tiempo que necesita la part´ıcula para quemarse completamente o abandonar la caldera, ˜qjes el caudal m´asico de carb´on que entra por j,pij es el porcentaje de part´ıculas ique entran por j, y NeyNpson el n´umero de entradas a la caldera y de part´ıculas, respectivamente. 2. Resoluci´on num´erica 2.1. Algoritmo y m´etodos num´ericos empleados En esta secci´on se describe el algoritmo empleado para la resoluci´on num´erica de las ecuaciones del modelo descrito anteriormente. Mientras que para las ecuaciones de la fase gaseosa se utiliza una descripci´on Euleriana, para la fase s´olida se sigue una descripci´on Lagrangiana. El algoritmo desarrollado sigue los siguientes pasos: 1. Se resuelven las ecuaciones (6)-(10) del modelo de fase gaseosa utilizando un m´etodo de elementos finitos P2 combinado con el m´etodo de las caracter´ısticas, de orden 2, para el tratamiento del t´ermino convectivo. Dependiendo de si se est´a en una zona con O2o sin ´el, el c´alculo de la temperatura y de las fracciones m´asicas de las especies que componen la mezcla gaseosa se realiza de la siguiente forma: a) En la regi´on ΩF,w6=w7=w8= 0 y, por tanto: Yg CO2=Xg 2−11 7Yg CO −44ν2 Mvol Yg V,(21) Yg H2O=Xg 3−18ν3 Mvol Yg V−9Yg H2,(22) Yg SO2=Xg 4−64ν4 Mvol Yg V.(23) Adem´as, para calcular Yg CO, Y g VyYg H2, la ecuaci´on (6) para Xg 1tiene que ser complementada con dos ecuaciones de conservaci´on para Yg CO, Y g VoYg H2.Por ejemplo, con las ecuaciones Lg(Yg V) = fm VyLg(Yg H2) = fm H2en ΩF,que deben ser integradas usando las condiciones de contorno Yg V= 0 y Yg H2= 0 en la superficie ΓFdada por Xg 1= 0, que separa las regiones ΩOy ΩF. 6 Simulaci´on num´erica de la combusti´on de carb´on pulverizado b) En la regi´on ΩOse tiene que Yg CO =Yg V=Yg H2= 0 y, adem´as, Yg O2=Xg 1, Yg CO2=Xg 2,Yg H2O=Xg 3,Yg SO2=Xg 4yhg T=Hg. 2. Se calculan los datos de la fase gaseosa que determinan la concentraci´on de O2,H2O yCO2en el entorno de la part´ıcula y, utilizando la ecuaci´on del movimiento de la part´ıcula, se calcula su posici´on. 3. Se resuelven las ecuaciones (11) utilizando un m´etodo de diferencias finitas de un paso (Euler impl´ıcito), obteni´endose los valores de λn+1 4yλn+1 5. 4. Una vez calculados los valores de ρn+1 V,ρn+1 H2O,λn+1 4yλn+1 5, si rn c>0 (lo que equivale a decir que todav´ıa hay carbono fijo en la part´ıcula) hay dos posibilidades: a) Si los n´umeros de Damk¨ohler Dai,i= 1,2,3, son menores que 1 se considera que las reacciones 1, 2 y 3 no tienen lugar y, por lo tanto λn+1 1=λn+1 2=λn+1 3= 0. b) Si los n´umeros de Damk¨ohler Dai,i= 1,2,3 son mayores que 1, las reacciones 1, 2 y 3 se consideran infinitamente r´apidas y, pueden ocurrir dos casos dependiendo de la atm´osfera en la que se encuentra la part´ıcula: Si Yg O2= 0 la part´ıcula est´a en la regi´on ΩFy se resuelve el sistema (12). Si Yg O26= 0 la part´ıcula est´a en la regi´on ΩO. Se resuelve en primer lugar el sistema (14) correspondiente a la situaci´on en la que la llama est´a fuera de la part´ıcula (rf> a). Si al resolver este sistema resulta que rf≤aes que nos encontramos en la situaci´on en la que la llama est´a dentro de la part´ıcula y se resolver´ıa el sistema (18), comprobando que, efectivamente rf≤a(y por tanto la consistencia en el modelo). Para resolver cualquiera de estos tres sistemas se utiliza el mismo m´etodo. En primer lugar discretizamos la ecuaci´on que determina la evoluci´on de rc utilizando el m´etodo de Euler expl´ıcito con la condici´on inicial r0 c=a. En segundo lugar se utiliza el m´etodo de Newton discretizado para resolver el sistema de ecuaciones correspondiente, que nos proporciona λn+1 1yλn+1 3. Si rn c= 0 se ha consumido todo el carbono fijo y, por lo tanto, no tendr´ıan lugar las reacciones de gasificaci´on del char, lo que implica que λn+1 1=λn+1 2=λn+1 3= rn+1 c=rn+1 f= 0. 5. Se resuelve la ecuaci´on de la energ´ıa que corresponda. Para su discretizaci´on se utiliza el m´etodo de Euler impl´ıcito. Debido al t´ermino de calor por radiaci´on, representado por la ley de Stefan-Boltzmann, queda una ecuaci´on no lineal que se resuelve mediante el m´etodo de Newton. 6. Se calculan las fuentes para la fase gaseosa dadas por la expresi´on (20). 2.2. Algunos resultados obtenidos Hemos aplicado la resoluci´on del modelo de combusti´on a la simulaci´on de una caldera de carb´on pulverizado de una Central T´ermica (ver Saavedra [5]). En primer lugar mostraremos algunos resultados obtenidos para una part´ıcula de 500 micras que entra en la caldera (Figura 2). A continuaci´on se muestra en un plano de la 7 A. Berm´udez, J.L. Ferr´ın, A. Li˜n´an, L. Saavedra caldera la fuente de energ´ıa homogeneizada y se compara con la obtenida mediante Fluent (Figuras 3 y 4). Figura 2: Velocidades adimensionales de gasificaci´on del “char” Contours of sensible-enthalpy-source FLUENT 6.2 (3d, segregated, pdf19, ske) Mar 22, 2007 3.000e+06 2.850e+06 2.700e+06 2.550e+06 2.399e+06 2.249e+06 2.099e+06 1.949e+06 1.799e+06 1.649e+06 1.499e+06 1.348e+06 1.198e+06 1.048e+06 8.979e+05 7.478e+05 5.976e+05 4.475e+05 2.973e+05 1.472e+05 -3.000e+03 ZY X Figura 3: SC3D Figura 4: Fluent Agradecimientos Parte de este trabajo ha sido financiado por el MEC a trav´es del proyecto ENE200509190-C04-01/CON. Referencias [1] A. Berm´udez de Castro. Continuum thermomechanics. Birkh¨auser, Verlag, Berlin, 2005. [2] A. Berm´udez de Castro, J. L. Ferr´ın y A. Li˜n´an. The modelling of the generation of volatiles, H2 and CO, and their simultaneous diffusion controlled oxidation, in pulverised furnaces. Aceptado en Combust. Theory Model. [3] S.P. Burke and T.E.W. Schumann. Diffusion flames. Ind. and Eng. Chemistry, 20:998-1004 (1928). [4] J. L. Ferr´ın. Algunas contribuciones a la modelizaci´on matem´atica de procesos de combusti´on de carb´on. Tesis. Universidade de Santiago de Compostela, 1999. [5] L. Saavedra. Simulaci´on num´erica de la combusti´on de part´ıculas de carb´on y simulaci´on num´erica en Mec´anica de Fluidos. Trabajo de Investigaci´on Tutelado. USC, 2006. 8