Control Óptimo Aplicado a Campos de Colectores Solares Distribuidos

Revista Iberoamericana de Automática e Informática Industrial 15 (2018) 327-338 www.revista-riai.org Control ´ Optimo Aplicado a Campos de Colectores Solares DistribuidosI Francisco R. Rubioa,∗, Sergio J. Navasa, Pedro Ollerob, Jo˜ ao M. Lemosc, Manuel G. Ortegaa aDepartamento de Ingenier´ıa de Sistemas y Autom´atica, Universidad de Sevilla, Spain bDepartamento de Ingenier´ıa Qu´ımica y Ambiental, Universidad de Sevilla, Spain cINESC-ID, Instituto Superior T´ecnico, University of Lisbon, Portugal Resumen En este trabajo se describe el control de una planta solar de colectores cil´ ındrico parab´ olicos cuando la radiaci´ on solar est´ a sometida a variaciones debidas al paso de nubes. El objetivo de las estrategias de control desarrolladas es optimizar la potencia generada, a diferencia de otras estrategias que persiguen el mantenimiento de la temperatura de salida del campo solar. El modelo desarrollado de la planta solar incluye tanto el campo solar con todos sus lazos, como el sistema de generaci´ on de potencia y el sistema de almacenamiento. As´ ı mismo se modela el paso de las nubes con diferentes tama˜ nos de cobertura del campo solar. Se desarrollan diferentes estrategias de control para maximizar la potencia generada y al mismo tiempo intentar producir dicha potencia el m´ aximo de tiempo posible y con las menores variaciones. Palabras Clave: Control de procesos, Control predictivo, Optimizaci´ on, Energ´ ıa solar. Optimal Control Applied to Distributed Solar Collector Fields Abstract This work describes the optimal control of a parabolic trough solar plant when the solar radiation is subject to variations due to the passage of clouds. The objective of the control strategies developed is to optimize the generated power, unlike other strategies that pursue the maintenance of the outlet oil temperature of the solar field. The solar plant model developed includes the solar field with all its loops, as well as the power generation system and the storage system. It also models the passage of clouds with different sizes of coverage of the solar field. Different control strategies are developed to maximize the power generated and at the same time to try to produce that power as long as possible and with the smallest variations. Keywords: Process Control, Predictive control, Optimization, Solar energy. 1. Introducci´ on Desde hace alg´ un tiempo se ha puesto de manifiesto un creciente inter´ es por el aprovechamiento de fuentes energ´ eticas no exploradas hasta ahora, o que lo han sido insuficientemente. Tal es el caso de la energ´ ıa solar, y concretamente de las plantas el´ ectricas termosolares. Las Centrales El´ ectricas Termosolares (CET), son sistemas empleados para obtener la energ´ ıa el´ ectrica a partir de la energ´ ıa solar mediante la transformaci´ on previa de ´ esta en energ´ ıa t´ ermica. Una de estas plantas son las de energ´ ıa solar concentrada (CSP). Los colectores parab´ olicos, las torres solares, el colector de Fresnel y los platos solares son las tecnolog´ ıas m´ as utilizadas para concentrar la energ´ ıa solar. Este trabajo se centra en los campos solares parab´ olicos, que consisten en un campo de colectores (Figura 1), un ciclo de potencia y elementos auxiliares como bombas, tuber´ ıas, tanques de almacenamiento y v´ alvulas. El campo de colectores solares recoge la radiaci´ on solar y la enfoca en un tubo por el que circula un fluido caloportador, usualmente aceite sint´ etico. El aceite se calienta y luego se utiliza en el ciclo de potencia para producir vapor de alta presi´ on en una caldera y se produce la electricidad mediante la expansi´ on en un turbo-generador. En los ´ ultimos a˜ nos tambi´ en IEste trabajo ha sido soportado por los proyectos DPI2013-44135-R y DPI2015-70973-R del Ministerio Espa˜ nol de Ciencia e Innovaci´ on. ∗Autor para correspondencia: [email protected] To cite this article: Francisco R. Rubio, Sergio J. Navas, Pedro Ollero, João M. Lemos, Manuel G. Ortega. 2018. Optimal Control Applied to Distributed Solar Collector Fields. Revista Iberoamericana de Automática e Informática Industrial 15, 327-338. https://doi.org/10.4995/riai.2018.8944 Attribution-NonCommercial-NoDerivatives 4,0 International (CC BY-NC-ND 4,0) se han a˜ nadido a las plantas tanques de almacenamiento para prolongar el tiempo de conexi´ on a la red el´ ectrica y/o para fines de regulaci´ on. Actualmente la posibilidad de almacenar energ´ ıa es una de las mayores ventajas de las plantas CET. Figura 1: Campo de colectores distribuidos ACUREX La diferencia fundamental entre una planta convencional y una planta solar es que la energ´ ıa primaria no puede ser manipulada por el hombre, es intermitente y cuando se encuentra disponible resulta cara de transformar. Por todo ello, en una planta solar se requiere un sofisticado sistema de control capaz de mantener unas especificaciones m´ as rigurosas que en el caso de un control cl´ asico. Durante los ´ ultimos treinta a˜ nos, muchos investigadores han dedicado un esfuerzo considerable para mejorar la eficiencia de las plantas de energ´ ıa solar t´ ermica con colectores distribuidos desde los puntos de vista de control y optimizaci´ on (Camacho, 1997), (Camacho, 2012), (Lemos, 2014). El objetivo principal de un campo solar parab´ olico es recoger la energ´ ıa solar con el fin de producir la mayor cantidad posible de energ´ ıa el´ ectrica. Normalmente, la mayor´ ıa de las centrales t´ ermicas solares tratan de lograr este objetivo manteniendo la temperatura del aceite de salida del campo alrededor del valor m´ aximo permisible, que en este caso es de 400 ◦C, impuesto para evitar la degradaci´ on del aceite. Sin embargo, algunos estudios como (Lippke, 1995) muestran que esta forma de operar el campo no produce los mejores resultados de la energ´ ıa el´ ectrica generada. En este trabajo se sugiere que la estrategia ´ optima se basa en adaptar la temperatura de salida del aceite a la radiaci´ on solar incidente, manteniendo constante la temperatura de sobrecalentamiento del vapor, mientras que en (Camacho, 2013) se propuso cambiar la temperatura de salida seg´ un el valor de la radiaci´ on solar. En (Lippke, 1995) la variable controlada es la temperatura de sobrecalentamiento, mientras que en (Camacho, 2013) es la temperatura de salida del aceite. La novedad principal de este art´ ıculo es que el efecto de la radiaci´ on solar sobre la temperatura ´ optima de salida se estudia con un modelo de ciclo de potencia completo. Adem´ as se estudia el caso cuando se tiene el campo solar parcialmente cubierto por las nubes, lo cual puede hacer que tengamos unos lazos con radiaci´ on solar y otros totalmente sin radiaci´ on. Por otro lado la planta solar se complementa con un sistema de almacenamien- to de energ´ ıa mediante dos tanques de sales fundidas, lo cual permitir´ a tener otra variable de control que es caudal de aceite (y por ende la energ´ ıa) que se deriva para almacenamiento. En este trabajo se propone el uso de una funci´ on polinomial para calcular la energ´ ıa el´ ectrica generada en funci´ on del caudal m´ asico y de la temperatura de salida del aceite del campo, en lugar de depender s´ olo de esa temperatura. Tambi´ en se utiliza un modelo del campo solar completo para tener en cuenta no s´ olo la radiaci´ on incidente total, sino tambi´ en su distribuci´ on entre cada uno de los lazos que constituyen el campo solar. Con este modelo es posible simular cada lazo del campo, en lugar de simular s´ olo uno de ellos suponiendo que el comportamiento de todo el campo es el mismo. En a˜ nos anteriores se ha realizado el modelado de campos solares cil´ ındrico parab´ olicos por varios autores (Camacho, 1988), (Meaburn, 1993), (Rubio, 2006), (Lemos, 2006), (Camacho, 2013), (Romera, 2017). El uso de un modelo del campo solar que tiene en cuenta individualmente todos sus lazos fue propuesto en (Abutayeh, 2014), pero se us´ o para probar una estrategia de control basada en maximizar la temperatura de salida del aceite del campo, que no es la forma ´ optima de producir la m´ axima potencia el´ ectrica. En este trabajo se utiliza este tipo de modelo del campo para comparar las estrategias de control cuyo principal objetivo es maximizar la potencia el´ ectrica, especialmente durante los d´ ıas con cobertura parcial. Como se ver´ a el objetivo de potencia m´ axima no coincide con el de mantener la temperatura de salida a un valor determinado. Las dos primeras estrategias consisten en un controlador predictivo multivariable (MPC), que utiliza predicciones de las futuras nubes junto con el modelo del campo de colectores y del ciclo de potencia, aunque difieren en su n´ umero de variables manipuladas. Mientras que en una de las estrategias propuestas se manipula el flujo total de aceite, que despu´ es se distribuye igualmente entre los lazos, la otra manipula individualmente el flujo de aceite que circula a trav´ es de cada lazo. Con ambas estrategias, se consigue una mejora de la potencia el´ ectrica generada, en comparaci´ on con la estrategia de mantener constante la temperatura de salida del aceite; sin embargo, se ver´ a que la estrategia que manipula individualmente el flujo de cada lazo no produce una mejora notable en comparaci´ on con la que manipula el flujo total. Figura 2: Campo solar cubierto parcialmente Utilizando el software Engineering Equation Solver c , los autores han simulado un campo solar parab´ olico conectado a un ciclo de potencia. Estas simulaciones se han realizado para diferentes valores de radiaci´ on solar, calculando en el r´ egimen permanente, para cada una el valor ´ optimo de la temperatura de salida que produce la mayor cantidad de energ´ ıa el´ ectrica. Este valor ´ optimo tambi´ en corresponde al valor ´ optimo del flujo de aceite, debido a que ambas variables son dependientes para un valor dado de Irradiancia (I). El ciclo de potencia utilizado en este trabajo es un ciclo Rankine con una temperatura m´ axima 328 Francisco R. Rubio et al. / Revista Iberoamericana de Automática e Informática Industrial 15 (2018) 327-338 de 374◦C, una presi´ on m´ axima de 70,5 bar y un rango de potencia de 800-2330 kW. La selecci´ on de un ciclo de Rankine se debe al hecho de que es el mismo utilizado en el campo ACUREX; sin embargo, el rango de potencia se ha incrementado de acuerdo con el aumento de tama˜ no en el campo utilizado en este trabajo y la temperatura de trabajo ha sido aumentada para coincidir con los utilizadas por las plantas comerciales (Navas, 2017). El art´ ıculo est´ a organizado de la siguiente manera: En la secci´ on 2 se describen los modelos del campo solar, del paso de nubes, del ciclo de potencia y el del sistema de almacenamiento utilizados con fines de simulaci´ on. En la secci´ on 3 se describen las estrategias de control probadas: la estrategia global, que consiste en un controlador MPC con una sola variable manipulada (flujo total de aceite), la estrategia distribuida, que consiste en un controlador MPC con 24 variables manipuladas (flujo de aceite a trav´ es de cada lazo) y la estrategia de control ´ optimo del campo con almacenamiento de sales. La secci´ on 5 muestra los resultados obtenidos en simulaci´ on. Por ´ ultimo, el documento llega a su fin con algunas observaciones finales (secci´ on 6). 2. Modelo del Sistema En este apartado se presenta el modelo de cada una de las partes que se han utilizado para simular la operaci´ on de un campo solar durante d´ ıas con cobertura parcial. Estas partes son: el campo de colectores solares, el paso de las nubes, el sistema de almacenamiento y el ciclo de generaci´ on de energ´ ıa. 2.1. Modelo del Campo de Colectores Solares El modelo del campo del colector solar es el mismo utilizado en (Navas, 2016), siendo al mismo tiempo una ligera modificaci´ on del modelo propuesto por (Camacho, 1997)-(Camacho, 2012)-(Carmona, 1985) para el campo ACUREX (Figura 1). Este modelo puede usarse para simular campos solares parab´ olicos mediante la selecci´ on de par´ ametros como el n´ umero de zonas activas (las partes donde la radiaci´ on solar llega al tubo) y pasivas (articulaciones y otras partes no alcanzadas por la radiaci´ on solar concentrada), la longitud de cada zona, o la abertura del colector. El campo solar simulado en este trabajo est´ a modelado utilizando datos de radiaci´ on solar que corresponden a la localizaci´ on de la Escuela T´ ecnica Superior de Ingenier´ ıa de Sevilla. Se compone de 24 lazos y tiene las dimensiones de 144x240 m. Cada lazo ( j) est´ a modelado por el siguiente sistema de ecuaciones diferenciales parciales que describen el balance energ´ etico: Zonas Activas: Para un elemento de metal: ρmCmAm ∂Tm j ∂t=Ijµ0G−GHl j(Tm j −Tamb)−d Ht j(Tm j −Ta j) (1) Para un elemento de fluido (aceite): ρa jCa jAa ∂Ta j ∂t+ρa jCa j ˙qa j ∂Ta j ∂x=d Ht j(Tm j −Ta j) (2) donde se han obtenido los coeficientes, en funci´ on de la temperatura, a partir de los datos suministrados por el productor del aceite, ρa j =903−0,672 Ta j (kg/m3); Ca j =1820+3,478 Ta j(J/kgoC) (3) Zonas Pasivas: Para un elemento de metal: ρmCmAm ∂Tm j ∂t=−GHp j(Tm j −Tamb)−d Ht j(Tm j −Ta j) (4) En las ecuaciones anteriores el sub´ ındice mse refiere al metal, el aal fluido de aceite y jal lazo. Los par´ ametros del modelo y sus unidades pueden verse en la tabla 1. Tabla 1. Par´ ametros del modelo y variables S´ ımbolo Descripci´ on y unidad t: tiempo (s) x: espacio (m) ρ: densidad (kg/m3) C: calor espec´ ıfico (J/kgK) A:´ area transversal (m2) T: temperatura (oC) ˙q: flujo volum´ etrico impulsado por la bomba (m3/s) I: irradiaci´ on solar (W/m2) µo: eficiencia ´ optica de los espejos G: apertura ´ optica de los espejos (m) Tamb : temperatura ambiente (K, oC) Hl: coeficiente global de p´ erdidas t´ ermicas, zonas activas(W/moC, W/mK) Hp: coeficiente global de p´ erdidas t´ ermicas, zonas pasivas(W/moC, W/mK) Ht: coeficiente de transmisi´ on metal-fluido (W/m2oC, W/m2K) d: di´ ametro interior de los tubos (m) El detalle de los diferentes par´ ametros puede encontrase en (Camacho, 1997)-(Camacho, 2012)-(Carmona, 1985). El colector tiene una longitud de 480 m(432 mde zonas activas y 48 m de zonas pasivas) y un ´ area de secci´ on transversal de 5,3∗10−4 m2. La abertura del colector es de 1.82 my la eficiencia ´ optica es 0,675. El coeficiente de p´ erdidas t´ ermicas para las zonas activas, Hl, y para las zonas pasivas, Hp, tiene un valor de 0,49 y 0,24 respectivamente. Se supone que estos valores son constantes, aunque var´ ıan con las temperaturas y como la temperatura ambiente tiene muy poca influencia en la eficiencia del campo solar (Manzolini, 2012), esta suposici´ on es aceptable para los prop´ ositos del presente trabajo. La soluci´ on de este sistema de ecuaciones en derivadas parciales se ha programado como una ecuaci´ on de diferencias finitas en dos etapas, considerando cada segmento de 1 mpara las zonas pasivas y de 3 mpara las zonas activas (Camacho, 1997) Francisco R. Rubio et al. / Revista Iberoamericana de Automática e Informática Industrial 15 (2018) 327-338 329 Francisco R. Rubio et al. / Revista Iberoamericana de Automática e Informática Industrial 15 (2018) 327-338 329 Figura 3: Diagrama de la planta solar con el ciclo de potencia y el sistema de almacenamiento resolviendo (1)-(2)-(4). Por ello, el modelo tiene 432/3=144 trozos en la zona activa y 48 en la zona pasiva. Una vez resuelto el sistema de ecuaciones diferenciales, se puede calcular el caudal volum´ etrico de aceite total (˙qat), la temperatura de salida del campo (Ta) y el flujo m´ asico de aceite total ( ˙mat), ˙ma j =˙qa j ρa j; ˙mat = 24 X j=1 ρa j ˙qa j;Ta= 24 X j=1 ˙ma j Ta j ˙mat (5) Este modelo de campo solar est´ a conectado con un modelo del ciclo de potencia (Figura 3) como se comenta en la secci´ on 1 y se describe m´ as adelante. As´ ı mismo se ha modelado un sistema de almacenamiento de energ´ ıa en dos tanques de sales fundidas. Se dispone de una v´ alvula de tres v´ ıas que puede ser controlada para derivar el caudal total entre el ciclo y el almacenamiento ( ˙mat =˙mac +˙maa). Los modelos tienen que ser simulados simult´ aneamente dado que algunas de sus variables son compartidas, siendo la simulaci´ on de los modelos un proceso iterativo. 2.2. Modelado del paso de las nubes En este trabajo se estudia c´ omo afecta el paso de nubes sobre el campo solar y c´ omo controlarlo en ese caso. La idea es modelar el paso de nubes que no cubren todo el campo y con ello conocer c´ omo se distribuye la radiaci´ on solar en todo el campo. La propuesta es realizar un modelo por trozos, para ello se propone utilizar una matriz que represente la extensi´ on del campo completo. A cada elemento de la matriz se le asigna el valor de la radiaci´ on solar incidente en esa secci´ on en cada momento. El campo solar tiene dimensiones de 144x240 m, por lo que si el campo est´ a dividido en elementos de 3x3 m, se necesita una matriz de 48x80. La matriz representa el campo de tal manera que cada elemento contiene una fracci´ on de un elemen- to activo o pasivo. La figura 4 muestra la fracci´ on de toda la matriz que cubre un lazo de un campo gen´ erico. Figura 4: Ejemplo de la matriz para un lazo del campo El valor de la radiaci´ on en cada elemento de la matriz depende de los siguientes par´ ametros: La radiaci´ on solar incidente global en el campo. La direcci´ on de la ruta seguida por la nube (el ´ angulo formado por esta direcci´ on y la direcci´ on seguida por el fluido t´ ermico, siendo 0 grados la direcci´ on en que coinciden y 90 grados la perpendicular). La velocidad de la nube (determinada por el n´ umero de elementos de la matriz supuestamente recorrida por la nube cada 39 segundos, que es el tiempo de muestreo del campo). 330 Francisco R. Rubio et al. / Revista Iberoamericana de Automática e Informática Industrial 15 (2018) 327-338 El tama˜ no de la nube (se supone una forma rectangular, definida por el n´ umero de filas y columnas que cubre). El factor de atenuaci´ on, que es el valor que multiplica la radiaci´ on en las zonas donde est´ a presente la nube. Este valor var´ ıa entre 0 y 1 dependiendo de la radiaci´ on que la nube permita llegar al campo, siendo 0 el caso en el que ninguna radiaci´ on alcance el campo, y 1 el caso en el que toda la radiaci´ on lo alcanza. Una vez establecidos estos valores, el programa calcula para cada colector el valor medio de radiaci´ on de cada elemento de la matriz relacionado con ´ el. Este valor medio se asigna a la variable Ide (1). 2.3. Modelo del ciclo de potencia El ciclo de potencia de Rankine, simulado en este trabajo, consiste en un economizador, una caldera y un sobrecalentador (cada uno indicado por el sub´ ındice E,BySrespectivamente) seguido por una turbina de vapor. Siguiendo el recorrido realizado por el flujo de aceite (Figura 3), primero entra en el sobrecalentador donde el flujo de vapor se sobrecalienta a una cierta temperatura. Este sobrecalentador fue modelado por las ecuaciones ˙macCPaTae −QS=˙macCPaTas, ˙mvCPvTve +QS=˙mvCPvTvs, QS=USASDT LMS,(6) que forman un sistema de tres ecuaciones con tres inc´ ognitas. Las inc´ ognitas son las temperaturas del aceite de salida y del vapor de salida, Tas yTvs, respectivamente, y el calor intercambiado en el sobrecalentador QS. Los intercambiadores de calor (E,ByS) se han dise˜ nado siguiendo la regla t´ ıpica de que la diferencia de temperatura en cada salida de cada intercambiador de calor para el punto de funcionamiento nominal de energ´ ıa (900W/m2, radiaci´ on solar y 393 ◦C de temperatura de salida del aceite), sea mayor de 10◦C. Esto permite establecer los valores del producto de la transmitancia t´ ermica Uy del ´ area de intercambio Ade cada intercambiador (tabla 2). Se supone que estos valores son constantes, asumiendo los errores asociados debido a la ligera variaci´ on de U. El calor espec´ ıfico del aceite CPa y el del vapor CPv dependen de la temperatura del fluido, siendo el primero calculado por la funci´ on polin´ omica dada en (Camacho, 1997) y el segundo utilizando las tablas de propiedades de vapor. El caudal m´ asico ˙mac y la temperatura de entrada Tae del aceite se obtienen del modelo del campo solar y el caudal m´ asico ˙mvy la temperatura de entrada del vapor Tve del modelo de la caldera que se mostrar´ a a continuaci´ on. Finalmente, la diferencia de temperatura media logar´ ıtmica DT LM necesaria en una de las ecuaciones se calcula usando la ecuaci´ on DT LM =(Tae −Tvs)−(Tas −Tve) ln (Tae−Tvs) (Tas−Tve).(7) Tabla 2: Valores de UA Par´ ametro Valor (kJ/s◦C) USAS20,15 UBAB96,95 UEAE25 Despu´ es del sobrecalentador, el caudal de aceite se utiliza para calentar el caudal de agua para producir vapor saturado en la caldera. Este proceso se lleva a cabo a una presi´ on flotante para maximizar la cantidad de calor intercambiado. Las ecuaciones utilizadas para modelar la caldera son ˙macCPaTae −QB=˙macCPaTas, ˙mwCPwTwe +QB=˙mvHv+˙mwCPwTws, QB=UBABDT LMB.(8) De nuevo, estas ecuaciones forman un sistema de tres ecuaciones con tres inc´ ognitas, que son las temperaturas de salida del aceite y la del vapor, Tas yTvs respectivamente, y el calor intercambiado en la caldera, QB. Debido al hecho de que se considera un modelo de estado estacionario, el flujo de agua ˙mwes igual al flujo de vapor ˙mv. El calor de vaporizaci´ on Hvse puede encontrar en las tablas de vapor o calculado por funciones polin´ omicas dependientes de la presi´ on o temperatura. Las variables restantes se pueden obtener como en el p´ arrafo anterior. Finalmente, el flujo de aceite se introduce en un economizador para precalentar el flujo de agua, y luego se recircula al campo. El economizador es modelado por las ecuaciones ˙macCPaTae −QE=˙macCPaTas, ˙mwCPwTwe +QE=˙mwCPwTws, QE=UEAEDT LME.(9) En este caso, como en los anteriores, tenemos un sistema de tres ecuaciones y tres inc´ ognitas, que son las temperaturas de salida del aceite y la del agua de salida, Tas yTws, respectivamente, y el calor intercambiado en el economizador QE. Las otras variables se calculan como se explica en el sobrecalentador. El vapor generado en el sobrecalentador fluye hacia una turbina de alta presi´ on. Para modelar esta turbina se utiliza el M´etodo de L´ınea de Willan descrito en (Smith, 2005). En primer lugar, la presi´ on de salida de la turbina se establece en un valor de 5,63 kPa, que es un valor t´ ıpico para las turbinas de condensaci´ on, ya que permite el uso de agua a temperatura ambiente para condensar el vapor. Despu´ es se utiliza la siguiente ecuaci´ on para calcular la presi´ on de entrada, Pe=k˙mv(10) Esta ecuaci´ on es una modificaci´ on de la ecuaci´ on de Stodola (Stodola, 1945) cuando la presi´ on de entrada es mucho m´ as alta que la de salida, como es el caso. Con estos valores de presi´ on es posible calcular las temperaturas de saturaci´ on del caudal de vapor a la entrada y a la salida, de manera que la diferencia de temperatura de saturaci´ on ∆Tspodr´ ıa calcularse y posteriormente utilizada en a=0,662∆Ts b=1,191 +0,000759∆Ts(11) para calcular los par´ ametros aybsi la potencia generada por la turbina es menor o igual a 2000 kW, y en a=−463 +3,53∆Ts b=1,220 +0,000148∆Ts(12) si es superior a 2000 kW. Francisco R. Rubio et al. / Revista Iberoamericana de Automática e Informática Industrial 15 (2018) 327-338 331 Francisco R. Rubio et al. / Revista Iberoamericana de Automática e Informática Industrial 15 (2018) 327-338 331 A continuaci´ on, es necesario calcular la diferencia isentr´ opica de la entalp´ ıa ∆Hiy establecer el par´ ametro Lque es la tasa de interceptaci´ on que depende de las caracter´ ısticas de la turbina, siendo su valor t´ ıpico entre 0,05 y 2. Con estos dos par´ ametros y con aybpodemos calcular la pendiente y la intercepci´ on de linea de Willans (nyWint) usando n=(L+1) b(∆Hi−a/˙mv), Wint =L b(∆Hi˙mv−a),(13) y con ellos, utilizando W=n˙mv−Wint,(14) podemos calcular la potencia el´ ectrica generada por la turbina. El modelo se ha simulado utilizando el software Engineering Equation Solver c y luego con estos resultados de simulaci´ on se obtuvieron dos funciones polinomiales de la temperatura de salida del aceite y el caudal de masa para calcular la potencia el´ ectrica generada (W) y la temperatura del aceite devuelto al campo (Tin) respectivamente. W=8230 −49,96 ˙mac −2,7 ˙m2 ac −47,15Ta +0,068T2 a+0,54 ˙macTa,(15) Tin =340 +1,78 ˙mac −0,155 ˙m2 ac −Ta +0,0011T2 a+0,022 ˙macTa,(16) Esta simplificaci´ on fue necesaria para reducir el tiempo de simulaci´ on, manteniendo el error en torno al 1 %. Adem´ as, se asumi´ o que la din´ amica del ciclo era como un modelo de primer orden con una constante de tiempo de 100 segundos, es decir, la din´ amica de la parte m´ as lenta del ciclo, la caldera. El ciclo de potencia se model´ o utilizando el software Aspen Hysys c y la constante de tiempo se obtuvo mediante simulaci´ on. Se asume que la constante de tiempo del ciclo es la misma que la de la caldera ya que es la din´ amica m´ as lenta, y ´ esta se obtiene mediante simulaci´ on con el software de simulaci´ on de procesos mencionado. El uso de ecuaciones est´ aticas m´ as la adici´ on de din´ amica es de uso com´ un en este tipo de sistemas. La potencia el´ ectrica neta producida por el campo es el resultado de sustraer la potencia consumida por la bomba a la potencia generada por la turbina. Por lo tanto es necesario calcular el consumo de la bomba, la cual se ha obtenido utilizando la ecuaci´ on de Darcy para estimar la p´ erdida de presi´ on Plj=9806,65 8f l˙q2 a j gπ2d5.(17) El n´ umero de Reynolds Re y el coeficiente de fricci´ on de Barr fson calculados por Rej=ρf˙qa jd Afµfj=0,25 1 (log10( r 3,7d+5,74 Re0,9 j ))2(18) El consumo de energ´ ıa Wbomba depende de la eficiencia de la bomba ηbomba, del flujo de aceite ˙qy de la perdida de presi´ on hpl Wbomba =P24 j=1˙qa j Plj ηbomba .(19) Por ´ ultimo, la potencia el´ ectrica neta es Wnet =W−Wbomba 1000 (kW).(20) Un resumen de los par´ ametros del modelo y sus unidades se muestra en la Tabla 3. Los valores num´ ericos pueden encontrarse en (Navas, 2017). Tabla 3: Par´ ametros del modelo de ciclo de alimentaci´ on y variables S´ ımbolo Descripci´ on Unidades QCalor intercambiado kW ˙mFlujo de masa kg/s CCapacidad calor´ ıfica espec´ ıfica kJ/(K kg) AArea de intercambio m2 TTemperatura ◦C UTransmitancia t´ ermica kW/(m2◦C) DT LM Diferencia de temperatura media logar´ ıtmica ◦C HEntalp´ ıa espec´ ıfica kJ/kg PPresi´ on kPa aPar´ ametro de correlaci´ on para turbinas de vapor kW bPar´ ametro de correlaci´ on para turbinas de vapor adimensional LTasa de interceptaci´ on adimensional nPendiente de la l´ ınea de Willans kW/kg Wint Intercepci´ on de l´ ınea de Willans kW ρfDensidad del fluido kg/m3 ˙qFlujo volum´ etrico de aceite m3/s AfArea de la secci´ on transversal m2 µViscosidad din´ amica kg/(m s) dDi´ ametro del tubo m rRugosidad relativa m gGravedad m/s2 lLongitud del lazo m Pl Ca´ ıda de presi´ on Pa La validaci´ on del modelo propuesto en (Camacho, 1997) y que se ha utilizado aqu´ ı como base, se realiz´ o con el campo ACUREX de la planta de Tabernas (Almer´ ıa). El modelo propuesto aqu´ ı no ha sido validado directamente al no corresponder con una planta real. Se considera que al estar basado en el mencionado anteriormente es v´ alido a los efectos de prueba de los controladores presentados. 2.4. Modelo del sistema de almacenamiento El sistema de almacenamiento consiste en un intercambiador de calor y dos tanques de sales fundidas, cuya composici´ on es 40 % de NaNO3y 60 % de KNO3. Uno de los tanques contiene sal a 291◦C, que es la sal que se calentar´ a en el intercambiador de calor utilizando el aceite de salida del campo solar ˙maa, cuya energ´ ıa t´ ermica est´ a destinada a almacenarse. El flujo de sal, despu´ es de haber sido calentado en el intercambiador de calor, se env´ ıa al otro tanque, donde se almacena hasta que se necesite calentar el aceite de la planta solar para mantener la producci´ on de energ´ ıa el´ ectrica cuando no hay suficiente radiaci´ on solar. El esquema del sistema de almacenamiento se puede ver en la figura 3. 332 Francisco R. Rubio et al. / Revista Iberoamericana de Automática e Informática Industrial 15 (2018) 327-338 Cada tanque est´ a modelado por el balance de materia, ρs dV dt =˙mse −˙mss,(21) y un balance de energ´ ıa ρsCsVdTs dt =Cs˙mse(Tse −Ts)−Cs˙mss(Tss −Ts)−Hs(Ts−Tamb), (22) donde Ves el volumen del tanque, ˙mes el caudal m´ asico, Tses la temperatura de sal fundida en el tanque, Tamb es la temperatura ambiente, Hses el coeficiente de p´ erdidas t´ ermicas, ρses la densidad y Cses la capacidad calor´ ıfica espec´ ıfica y pueden ser calculados por, ρs=2090 −0,636Ts,Cs=1443 −0,172Ts.(23) El coeficiente de p´ erdidas t´ ermicas Hstiene un valor de 0.015 W/◦C. Este valor se ha calculado suponiendo que el tanque a la capacidad m´ axima de volumen pierde 1◦C de temperatura por d´ ıa. El intercambiador de calor ha sido modelado por las ecuaciones ˙maaCPaTae −QS t =˙maaCPaTas, ˙msCsTse +QS t =˙msCsTss, QS t =US tAS tDT LMS t,(24) que forman un sistema de tres ecuaciones y tres inc´ ognitas. Los valores desconocidos son las temperaturas del aceite de salida y de la sal, Tas yTss respectivamente, y el intercambio de calor QS t. El valor del producto de la transmitancia t´ ermica US t y el ´ area de intercambio AS t como se ha descrito anteriormente en el modelo del ciclo. Este valor es 161.32 kJ/s◦C. 3. Estrategias de control sin tanques de almacenamiento En primer lugar se aborda el problema sin los tanques de almacenamiento con dos estrategias consistentes en un controlador MPC con los mismos horizontes de predicci´ on y control, modelo de predicci´ on y restricciones, pero con diferentes variables manipuladas. El primero, denominado por los autores como la estrategia global (GS), manipula el flujo total de aceite que circula por el campo, suponiendo que este flujo est´ a igualmente dividido entre los 24 lazos del campo. El otro, denominado estrategia distribuida (DS), manipula por separado el flujo de aceite que circula a trav´ es de cada lazo, lo que significa que tiene 24 variables manipuladas. En ambas estrategias se dispone de una bomba principal m´ as v´ alvulas en cada lazo; manual en el caso GS y autom´ atico en el caso DS. Los lazos se tratan como no acoplados. 3.1. Estrategia global La estrategia de control global (GS) se basa en el uso de un controlador MPC, que calcula el valor ´ optimo del caudal de salida de aceite que maximiza la potencia el´ ectrica generada, teniendo en cuenta las limitaciones de temperatura y de caudal de aceite. Para reducir los tiempos de c´ alculo, el modelo de campo solar utilizado por el MPC es una simplificaci´ on del utilizado para simular el campo real (1)-(2). Se supone que el tubo solo se divide en 6 partes de 80 metros en lugar de las divisiones de 1 metro para zonas pasivas y 3 metros para las zonas activas. La predicci´ on de las nubes se ha realizado suponiendo que se conoce la posici´ on exacta de la nube para cada intervalo de muestreo del horizonte de predicci´ on, aunque se supone que el valor de la radiaci´ on solar incidente es constante para cada intervalo de muestreo, lo que agrega cierto grado de error entre el valor predicho y el real. La estructura de control se puede ver en la figura 5. Figura 5: Esquema de la estrategia de control global La energ´ ıa el´ ectrica generada se considera a lo largo del per´ ıodo de tiempo (k+1)hhasta (k+H)h, donde kdenota el tiempo discreto actual, Hes el horizonte de predicci´ on, y hes el intervalo de muestreo, que tiene un valor de 39 segundos. Por lo tanto, el objetivo es seleccionar el caudal que maximice la potencia neta producida (ecuaci´ on 20), mediante la funci´ on de costos siguiente, J(˙qk+H−1 atk )= H X i=1 Wnet(˙qat((k+i−1)h)).(25) sujeto a las siguientes condiciones: ˙qa j(min)≤˙qa j(k+i)≤˙qa j(max)∀i=0, . . . , H Ta j ≤Ta j(max)∀i=0, . . . , H La maximizaci´ on de Jse realiza con una estrategia de horizonte deslizante, utilizando un horizonte de predicci´ on de H =12 y un horizonte de control de 3. El valor del horizonte de predicci´ on se eligi´ o teniendo en cuenta el tiempo de establecimiento de la planta y despu´ es de hacer algunas simulaciones para diferentes valores del mismo y seleccionar el que permiti´ o producir la mayor cantidad de energ´ ıa el´ ectrica durante un d´ ıa de operaci´ on. Los resultados de estas simulaciones se pueden ver en la figura 6. Se seleccion´ o un horizonte de control de 3 debido al hecho de que con uno superior, la mejora lograda en la cantidad de energ´ ıa el´ ectrica era inferior al 0,1 % y tambi´ en aument´ o el tiempo de c´ alculo del optimizador. Francisco R. Rubio et al. / Revista Iberoamericana de Automática e Informática Industrial 15 (2018) 327-338 333 Francisco R. Rubio et al. / Revista Iberoamericana de Automática e Informática Industrial 15 (2018) 327-338 333 Horizonte de Predicción 1 2 3 4 5 6 7 8 9 10 11 12 Potencia Eléctrica (MW) 837 838 839 840 841 842 843 844 845 846 Figura 6: Potencia El´ ectrica generada para diferentes valores del horizonte de predicci´ on El procedimiento de optimizaci´ on din´ amica se puede expresar como el siguiente algoritmo: Se miden todas las variables necesarias del proceso. Se utiliza el modelo del campo solar simplificado para calcular la temperatura del aceite de salida. Es decir se utiliza el mismo modelo dado por las ecuaciones (1 a 4), pero en vez de simular el tubo de aceite en 144 trozos, solo se simulan 6 trozos, como se ha indicado anteriormente. El caudal de aceite (˙qat) es la variable independiente utilizada por el optimizador. Con los valores de la temperatura de salida de aceite y el caudal m´ asico de aceite, se usan las ecuaciones (15) y (16) para calcular la potencia el´ ectrica generada por la turbina (W) y la temperatura del aceite de entrada (Tin). Finalmente, con las ecuaciones (18), (17), (19) y (20), se calcula la potencia el´ ectrica neta (Wnet). Por lo tanto, el MPC tiene que maximizar el valor de la energ´ ıa el´ ectrica neta sujeta a restricciones en la temperatura del aceite (Ta≤400◦C) y el flujo de aceite (0,133 l/s≤˙qa j≤1,58 l/s) para cada ciclo. La optimizaci´ on se lleva a cabo para cada paso de integraci´ on mediante la funci´ on f mincon en MATLAB. El valor ´ optimo del flujo de aceite se env´ ıa a continuaci´ on como el nuevo punto de ajuste del controlador de caudal. Para el c´ alculo, se asume que la din´ amica del lazo de control de caudal se comporta como un sistema de segundo orden, con din´ amica mucho m´ as r´ apida que los 39 segundos del periodo de muestro que se utiliza. 3.2. Estrategia distribuida La estrategia de control distribuido (DS) es similar a la anterior, pero con la diferencia de que en lugar de calcular el flujo total de aceite del campo, el MPC calcula el valor del flujo de aceite que debe circular a trav´ es de cada lazo para maximizar la potencia el´ ectrica generada, teniendo en cuenta las restricciones de la temperatura y caudal de aceite. El esquema de control se puede ver en la figura 7. La optimizaci´ on tambi´ en se lleva a cabo de la misma manera que en la estrategia global, pero la salida del MPC ahora consiste en 24 variables, cada una correspondiente al valor del flujo de aceite de cada lazo. El aumento en el n´ umero de variables lleva a un aumento en el tiempo de c´ alculo en cada paso de integraci´ on del optimizador. El c´ alculo de cada paso de integraci´ on de todo el programa (incluyendo el MPC, campo solar, nubes pasantes y ciclo de energ´ ıa) requiere alrededor de 20 segundos para simular un per´ ıodo de 39 segundos y, por lo tanto, puede implementarse en una planta real. Figura 7: Estrategia de control distribuido 4. Estrategia de control para planta solar con almacenamiento En esta secci´ on se aborda el problema de optimizaci´ on del suministro de energ´ ıa de una planta con almacenamiento mediante tanques de sales fundidas, intentando maximizar las horas de almacenamiento para su posterior suministro a la red, sujeto a dar la potencia deseada en todo momento. El objetivo en este caso es seleccionar el flujo total de aceite ˙qat y el flujo m´ asico de aceite que se deriva ˙maa para el almacenamiento de forma que se maximice la funci´ on de costos siguiente (ver figura 8), J= H X i=1 [−λ1(Wre f −Wnet i)2−λ2(∆˙qat)2+(Wmax −Wnet i)].(26) donde Wmax es la m´ axima potencia que se obtendr´ ıa sin almacenamiento, o sea, si todo el caudal de aceite fuese al ciclo de potencia y nada a almacenamiento. El ´ ındice propuesto pondera el error entre la potencia que se produce en cada instante respecto a la referencia de potencia (Wre f ) que se pretende suministrar, penalizando el esfuerzo de control para evitar variaciones bruscas del caudal de aceite y por otro lado se incluye un ´ ultimo t´ ermino para maximizar la energ´ ıa almacenada en el tanque de sales (Wmax −Wnet i). Figura 8: Estrategia de control para planta con tanques de almacenamiento 334 Francisco R. Rubio et al. / Revista Iberoamericana de Automática e Informática Industrial 15 (2018) 327-338 La temperatura Taque se obtiene en cada paso es la correspondiente a la m´ axima energ´ ıa producida. Los valores de λ1 yλ2se eligen de forma que los tres t´ erminos del ´ ındice est´ en en el mismo rango num´ erico y posteriormente se ajustan para conseguir el comportamiento deseado. La maximizaci´ on de Jse realiza con una estrategia de horizonte deslizante, utilizando un horizonte de predicci´ on de H = 12 y un horizonte de control de 3. 5. Resultados de simulaci´ on En esta secci´ on se muestran los resultados obtenidos por simulaci´ on con los modelos descritos en las secciones anteriores. En primer lugar, el objetivo es comparar el rendimiento de las estrategias global y distribuida durante las situaciones de cobertura parcial del campo solar al pasar nubes. Las nubes se han simulado con el modelo descrito en la secci´ on 2.2 con los siguientes par´ ametros: una matriz de 16x16, una velocidad de 2 elementos de matriz por paso de integraci´ on, una direcci´ on de 45◦y una factor de atenuaci´ on de 0. Para cada simulaci´ on el tiempo de paso ha sido modificado. Este tiempo de paso se ha definido como el n´ umero de pasos de integraci´ on donde no hay nube sobre el campo despu´ es de que el ´ ultimo lo haya abandonado; eso significa que durante toda la simulaci´ on hay nubes que entran y salen del campo, separadas por este tiempo. El primer objetivo de estas simulaciones es determinar a partir de qu´ e valor de la radiaci´ on solar es mejor operar el campo a una temperatura m´ as baja y con un caudal superior de aceite en lugar de trabajar a la m´ axima temperatura pero con un caudal inferior. Por lo tanto, estas simulaciones se llevaron a cabo calculando la temperatura de salida ´ optima para cada valor de radiaci´ on solar. El valor m´ aximo permitido de la temperatura de salida es 400 ◦C para evitar la degradaci´ on del aceite, y el flujo m´ asico de aceite est´ a limitado entre 3,7 y 37 kg/s de acuerdo con la velocidad m´ ınima y m´ axima permitida a trav´ es de las tuber´ ıas (0,5-5 m/s) (ver (Navas, 2017)). Radicacion Solar (W/m 2) 500 550 600 650 700 750 800 850 900 950 1000 Tiempo(h) 10 11 12 13 14 15 16 17 Temperatura (ºC) 310 320 330 340 350 360 370 380 390 400 410 Figura 9: Evoluci´ on de la temperatura ´ optima durante un d´ ıa con nubes La figura 9 presenta el caso de un d´ ıa con nubes. En este caso, est´ a claro que, dependiendo de la cantidad de nubes y la cantidad del campo solar que est´ e cubierto, tendremos que cambiar el punto de operaci´ on (la referencia de temperatura de salida del aceite), muchas veces durante el d´ ıa para obtener la cantidad m´ axima de energ´ ıa el´ ectrica. Figura 10: radiaci´ on solar Las siguientes simulaciones est´ an encaminadas a probar las estrategias control global (GS) y control distribuido (DS). Estas simulaciones se han llevado a cabo con un tiempo de paso de nubes de 40 tiempos de muestreo, lo que da como resultado la curva de radiaci´ on solar que se muestra en la figura 10. Con esta curva de radiaci´ on, se evaluaron las estrategias de control junto con la utilizada para mantener constante la temperatura de salida del campo, para ver si la estrategia distribuida logr´ o una mejora en la potencia el´ ectrica generada. Pero, como se puede ver en las figuras 11, 12, y 13 ambas estrategias tienen un comportamiento similar, la mejora lograda por la distribuida es insignificante. Figura 11: Temperatura de salida del aceite con las estrategia de control global y distribuido Figura 12: Potencia el´ ectrica con las estrategias de control global y distribuido Francisco R. Rubio et al. / Revista Iberoamericana de Automática e Informática Industrial 15 (2018) 327-338 335 Francisco R. Rubio et al. / Revista Iberoamericana de Automática e Informática Industrial 15 (2018) 327-338 335