Software de simulación para la generación fotovoltaica

Equation Chapter 1 Section 1 Trabajo Fin de Máster Máster en Ingeniería Electrónica, Robótica y Automática Software de simulación para la generación fotovoltaica Dpto. de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Autor: Alejandro Trujillo Quevedo Tutor: Alejandro del Real Torres Sevilla, 2021 iii Trabajo Fin de Máster Máster en Ingeniería Electrónica, Robótica y Automática Software de simulación para la generación fotovoltaica Autor: Alejandro Trujillo Quevedo Tutor: Alejandro del Real Torres Profesor asociado Dpto. de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2021 v Trabajo Fin de Máster: Software de simulación para la generación fotovoltaica Autor: Alejandro Trujillo Quevedo Tutor: Alejandro del Real Torres El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2021 El Secretario del Tribunal vii A mi familia A mis maestros ix Agradecimientos Me gustaría, antes de todo, agradecer a mi tutor, Alejandro del Real Torres, por la oportunidad brindada, así como la ayuda y el apoyo prestados durante la realización de este proyecto. A mi familia por su inquebrantable confianza y su inestimable apoyo durante todos y cada uno de estos años, en los que me alentaban a no rendirme y seguir adelante. A mis amigos, los de toda la vida y las grandes amistades formadas durante estos años, por el apoyo incondicional y la ayuda prestada cuando más lo necesitaba. Modelo simplificado 35 Maximum Power Point Tracking 36 7 Matriz de módulos fotovoltaicos 37 Salida en CC 37 Salida en CA 38 8 Implementación en Julia 41 Diagrama de flujo 41 Documentación 43 8.1.1 Cuantitativo 43 8.1.2 Módulo cualitativo (usando el nombre) 45 8.1.3 Módulo cualitativo (usando eficiencia) 47 8.1.4 Módulo (usando el nombre) e inversor cualitativos 49 8.1.5 Módulo (usando la eficiencia) e inversor cualitativos 50 8.1.6 Módulo (usando el nombre), inversor y pérdidas cualitativas 51 8.1.7 Módulo (usando la eficiencia), inversor y pérdidas cualitativas 52 8.1.8 Inversor cualitativo 53 8.1.9 Inversor y pérdidas cualitativas 55 8.1.10 Pérdidas cualitativas 56 Descripción del paquete 58 9 Resultados 61 Modelo cualitativo 61 Modelo cuantitativo 67 Modelo cualitativo con GHI 72 Panel solar fotovoltaico 76 10 Conclusiones 79 11 Líneas de trabajo futuras 81 Referencias 83 Glosario 91 Apéndice I: Resumen del modelo 93 Esquema general 93 Entradas 93 Salidas 93 Módulo: Angles for horizontal coordinate system 94 Entradas 94 Salidas 94 Parámetros 94 Módulo: Horizontal Irradiance 95 Entradas 95 Salidas 95 Parámetros 95 Módulo: Plane of array irradiance 96 Entradas 96 Salidas 96 Parámetros 96 Módulo: Soiling derate 98 Entradas 98 Salidas 98 Parámetros 98 Módulo: Effective irradiance 99 Entradas 99 xvii Salidas 99 Parámetros 99 Módulo: Module power 100 Entradas 100 Salidas 100 Parámetros 100 Módulo: Net DC Array Power 102 Entradas 102 Salidas 102 Parámetros 102 Módulo: Net AC Array Power 105 Entradas 105 Salidas 105 Parámetros 105 Apéndice II: Modelos REST2 y Engerer2 107 REST2 107 Irradiancia directa 107 Irradiancia difusa y global 110 Engerer2 112 Apéndice III: Códigos en Julia 115 Main.jl 115 Show_Doc.jl 122 Load_Data.jl 123 Sun_Angles.jl 149 Weather_File_Data.jl 156 REST2.jl 157 Engerer2.jl 161 Tilted_Irradiance.jl 162 Shading_Soiling.jl 164 Module_Irradiance.jl 165 Module_Power.jl 169 Model.jl 171 Gross_DC_Power.jl 173 Net_DC_Power.jl 174 Gross_AC_Power.jl 175 Net_AC_Power.jl 176 Save_Data.jl 177 xix ÍNDICE DE TABLAS Tabla 1: Parámetros típicos del modelo REST2 25 Tabla 2: Matriz de pruebas del estándar IEC-61853 34 Tabla 3: Parámetros del modo cuantitativo 44 Tabla 4: Parámetros del módulo cualitativo (usando el nombre) 46 Tabla 5: Parámetros del módulo cualitativo (usando la eficiencia) 48 Tabla 6: Parámetros del módulo (usando el nombre) e inversor cualitativos 49 Tabla 7: Parámetros del módulo (usando la eficiencia) e inversor cualitativos 50 Tabla 8: Parámetros del módulo (usando el nombre), inversor y pérdidas cualitativas 51 Tabla 9: Parámetros del módulo (usando la eficiencia), inversor y pérdidas cualitativas 52 Tabla 10: Parámetros del inversor cualitativo 54 Tabla 11: Parámetros del inversor y pérdidas cualitativas 55 Tabla 12: Parámetros de las pérdidas cualitativas 57 Tabla 13: Lista de archivos en la carpeta Function_codes 59 Tabla 14: Entradas de la prueba cualitativa 61 Tabla 15: Parámetros de la prueba cualitativa 63 Tabla 16: Análisis percentil del APE de la prueba cualitativa 65 Tabla 17: Análisis percentil para la prueba cualitativa (tras eliminar amanecer y anochecer) 67 Tabla 18: Entradas de la prueba cuantitativa 68 Tabla 19: Parámetros de la prueba cuantitativa 70 Tabla 20: Análisis percentil para la prueba cuantitativa (tras eliminar amanecer y anochecer) 72 Tabla 21: Entradas de la prueba cualitativa con GHI 72 Tabla 22: Parámetros de la prueba cualitativa con GHI 74 Tabla 23: Análisis percentil para la prueba cualitativa con GHI (tras eliminar amanecer y anochecer) 76 Tabla 24: Métricas de validación del modelo de los paneles solares 78 Tabla 25: Entradas generales del modelo 93 Tabla 26: Salidas generales del modelo 93 Tabla 27: Entradas del módulo Angles for horitontal coordinate system 94 Tabla 28: Salidas del módulo Angles for horitontal coordinate system 94 Tabla 29: Parámetros del módulo Angles for horitontal coordinate system 94 Tabla 30: Entradas del módulo Horizontal Irradiance 95 Tabla 31: Salidas del módulo Horizontal Irradiance 95 Tabla 32: Parámetros del módulo Horizontal Irradiance 95 Tabla 33: Entradas del módulo Plane Of Array Irradiance 96 Tabla 34: Salidas del módulo Plane Of Array Irradiance 96 Tabla 35: Parámetros del módulo Plane Of Array Irradiance 97 Tabla 36: Entradas del módulo Soiling derate 98 Tabla 37: Salidas del módulo Soiling derate 98 Tabla 38: Parámetros del módulo Soiling derate 98 Tabla 39: Entradas del módulo Effective Irradiance 99 Tabla 40: Salidas del módulo Effective Irradiance 99 Tabla 41: Parámetros del módulo Effective Irradiance 99 Tabla 42: Entradas del módulo Module Power 100 Tabla 43: Salidas del módulo Module Power 100 Tabla 44: Parámetros del módulo Module Power 101 Tabla 45: Entradas del módulo Net DC Array Power 102 Tabla 46: Salidas del módulo Net DC Array Power 102 Tabla 47: Parámetros del módulo Net DC Array Power 103 Tabla 48: Entradas del módulo Net AC Array Power 105 Tabla 49: Salidas del módulo Net AC Array Power 105 Tabla 50: Parámetros del módulo Net AC Array Power 106 Tabla 51: Parámetros del modelo Engerer2 113 xxi ÍNDICE DE FIGURAS Figura 1: Evolución de la eficiencia de las células solares para cada tecnología 2 Figura 2: Potencia solar fotovoltaica en España 3 Figura 3: Partes de una célula fotovoltaica 7 Figura 4: Efecto de un fotón en una unión P-N 7 Figura 5: Tabla con los días del año para años bisiestos y no bisiestos 10 Figura 6: Husos horarios 10 Figura 7: Sistema de coordenadas eclíptico 11 Figura 8: Rotación (R), Precesión (P) y Nutación (N) 12 Figura 9: Sistema de coordenadas ecuatorial 13 Figura 10: Sistema de coordenadas ecuatorial (utilizando el ángulo horario) 14 Figura 11: Sistema de coordenadas horizontal 15 Figura 12: Ecuación del Tiempo 16 Figura 13: Irradiancia espectral 19 Figura 14: Principio de funcionamiento de las células multiunión 20 Figura 15: Eficiencia en función del número de uniones 20 Figura 16: Masa de aire 21 Figura 17: Variación de la constante solar a lo largo del año 22 Figura 18: Descomposición de la irradiancia extraterrestre 23 Figura 19: Efecto de la atmósfera en la irradiancia espectral 24 Figura 20: Ángulo de incidencia sobre un plano inclinado 25 Figura 21: Componentes de la irradiancia difusa 26 Figura 22: Sombreado de los paneles 27 Figura 23: Circuito equivalente para una celda solar 32 Figura 24: Estructura de un módulo 35 Figura 25: Curva V-I (roja) y curva V-P (azul) 36 Figura 26: Disposición de una matriz de paneles fotovoltaicos 37 Figura 27: Diagrama de flujo principal 41 Figura 28: Diagrama de flujo del cálculo de la irradiancia 42 Figura 29: Diagrama de flujo del cálculo de la energía 42 Figura 30: Ejemplo de archivo meteorológico 60 Figura 31: Histograma del AE para la prueba cualitativa 63 Figura 32: Histograma del APE para la prueba cualitativa 64 Figura 33: Histograma del SE para la prueba cualitativa 64 Figura 34: Histograma del AE para la prueba cualitativa (tras eliminar amanecer y anochecer) 66 Figura 35: Histograma del APE para la prueba cualitativa (tras eliminar amanecer y anochecer) 66 Figura 36: Histograma del SE para la prueba cualitativa (tras eliminar amanecer y anochecer) 67 Figura 37: Histograma del AE para la prueba cuantitativa (tras eliminar amanecer y anochecer) 70 Figura 38: Histograma del APE para la prueba cuantitativa (tras eliminar amanecer y anochecer) 71 Figura 39: Histograma del SE para la prueba cuantitativa (tras eliminar amanecer y anochecer) 71 Figura 40: Histograma del AE para la prueba cualitativa con GHI (tras eliminar amanecer y anochecer) 74 Figura 41: Histograma del APE para la prueba cualitativa con GHI (tras eliminar amanecer y anochecer) 75 Figura 42: Histograma del SE para la prueba cualitativa con GHI (tras eliminar amanecer y anochecer) 76 Figura 43: Validación del modelo del panel solar en Cocoa 77 Figura 44: Validación del modelo del panel solar en Eugene 77 Figura 45: Validación del modelo del panel solar en Golden 78 Figura 46: Esquema general del modelo 93 Figura 47: Esquema del módulo Angles for horitontal coordinate system 94 Figura 48: Esquema del módulo Horizontal Irradiance 95 Figura 49: Esquema del módulo Plane Of Array Irradiance 96 Figura 50: Esquema del módulo Soiling derate 98 Figura 51: Esquema del módulo Effective Irradiance 99 Figura 52: Esquema del módulo Module Power (cuantitativo) 100 Figura 53: Esquema del módulo Module Power (cualitativo 1) 100 Figura 54: Esquema del módulo Module Power (cualitativo 2) 100 Figura 55: Esquema del módulo Net DC Array Power (cuantitativo) 102 Figura 56: Esquema del módulo Net DC Array Power (cualitativo) 102 Figura 57: Esquema del módulo Net AC Array Power (cuantitativo) 105 Figura 58: Esquema del módulo Net AC Array Power (cualitativo) 105 xxiii Notación 𝑡𝑢𝑡𝑐 Tiempo en UTC ℎ𝑟 Hora local 𝑚𝑖𝑛 Minuto local 𝑡𝑧 Huso horario 𝑗𝑑𝑜𝑦 Fecha ordinal 𝑗𝑢𝑙𝑖𝑎𝑛 Fecha juliana 𝑚𝑛𝑙𝑜𝑛𝑔 Longitud media 𝑚𝑛𝑎𝑛𝑜𝑚 Anomalía media 𝑒𝑐𝑙𝑜𝑛𝑔 Longitud de la eclíptica 𝑜𝑏𝑙𝑒𝑞 Oblicuidad de la eclíptica 𝑟𝑎 Ascensión recta 𝛿 Declinación 𝑔𝑚𝑠𝑡 Tiempo sideral medio en Greenwich 𝑙𝑚𝑠𝑡 Hora sideral media local 𝐻𝐴 Ángulo horario 𝛼 Ángulo de elevación 𝛾 Ángulo azimutal 𝜃𝑧 Ángulo cenital 𝐸𝑂𝑇 Ecuación del tiempo 𝑡𝑠𝑢𝑛𝑟𝑖𝑠𝑒 Tiempo de amanecer 𝑡𝑠𝑢𝑛𝑠𝑒𝑡 Tiempo de atardecer 𝐴𝑀 Masa de aire 𝐺𝑠𝑐 Constante solar 𝐹𝑛 Variación estacional 𝐺𝑜𝑛 Irradiancia extraterrestre normal 𝐷𝑁𝐼 Irradiancia directa normal 𝐷𝐻𝐼 Irradiancia difusa horizontal 𝐺𝐻𝐼 Irradiancia global horizontal 𝛾𝑠 Ángulo acimutal del panel solar 𝛽𝑠 Ángulo de inclinación del panel solar 𝜃 Ángulo de incidencia 𝐺𝑏𝑖 Irradiancia directa en el plano del panel 𝜌 Albedo 𝐺𝑔𝑖 Irradiancia reflejada sobre el panel 𝐺𝑑𝑖 Irradiancia difusa sobre el panel 𝐺𝑖 Irradiancia total sobre el panel 𝑂 Pérdidas por suciedad 𝐺𝑏 Irradiancia directa en el plano del panel tras pérdidas por suciedad 𝐺𝑑 Irradiancia difusa en el plano del panel tras pérdidas por suciedad 𝐺𝑟 Irradiancia reflejada en el plano del panel tras pérdidas por suciedad 𝐺 Irradiancia global en el plano del panel tras pérdidas por suciedad 𝐾𝑏 Modificador del ángulo de incidencia de la irradiancia directa 𝐾𝑑 Modificador del ángulo de incidencia de la irradiancia difusa 𝐾𝑟 Modificador del ángulo de incidencia de la irradiancia reflejada 𝐺𝑎𝑏𝑠 Irradiancia total tras pérdidas por modificador del ángulo de incidencia 𝜏𝛼 Producto transmitancia-absorbancia 𝑀 Modificador de la masa de aire 𝐺𝑒𝑓𝑓 Irradiancia efectiva finalmente absorbida por la célula 𝑣𝑤 Velocidad del viento 𝑇𝑁𝑂𝐶𝑇 Temperatura NOCT 𝑇𝑎 Temperatura ambiente 𝜂𝑟𝑒𝑓 Eficiencia de la célula en condiciones STC 𝑇𝑐 Temperatura de la célula 𝐼 Corriente de la célula 𝑉 Tensión de la célula 𝐼𝐿 Corriente fotogenerada 𝑎 Factor de idealidad 𝐼𝑜 Corriente de saturación del diodo 𝑅𝑠 Resistencia en serie 𝑅𝑠ℎ Resistencia en paralelo 𝑘 Constante de Boltzmann 𝑞 Carga del electrón 𝐼𝑠𝑐 Corriente de cortocircuito 𝑉𝑜𝑐 Tensión de circuito abierto 𝐸𝑔 Energía de banda prohibida 𝑃 Energía de la célula 𝑁𝑠 Número de células en serie 𝑁𝑝 Número de células en paralelo 𝜇𝐼𝑆𝐶 Coeficiente de variación de la corriente de cortocircuito con la temperatura 𝜇𝑉𝑂𝐶 Coeficiente de variación de la tensión de circuito abierto con la temperatura 𝜇𝑃 Coeficiente de variación de la energía con la temperatura 𝑉𝑑𝑐,𝑛 Tensión de la submatriz 𝑁𝑚,𝑠𝑡𝑟𝑖𝑛𝑔 Número de paneles por cadena 𝐼𝑑𝑐,𝑛 Corriente de la submatriz 𝑁𝑠𝑡𝑟𝑖𝑛𝑔𝑠 Número de cadenas 𝑃𝑑𝑐,𝑔𝑟𝑜𝑠𝑠,𝑛 Energía bruta generada en CC en la submatriz 𝑁𝑚𝑜𝑑𝑢𝑙𝑒𝑠 Número total de paneles 𝐹𝑑𝑐 Factor de pérdidas en CC 𝐿𝑚𝑖𝑠𝑚𝑎𝑡𝑐ℎ Pérdidas por desajuste (mismatch) 𝐿𝑛𝑎𝑚𝑒𝑝𝑙𝑎𝑡𝑒 Pérdidas nominales (nameplate) 𝐿𝑑𝑐𝑤𝑖𝑟𝑖𝑛𝑔 Pérdidas en el cableado en CC 𝑃𝑑𝑐,𝑛𝑒𝑡,𝑛 Energía neta generada en CC en la submatriz xxv 𝑃𝑑𝑐,𝑛𝑒𝑡 Energía neta generada en CC 𝑃𝑚𝑎𝑥𝑖𝑛𝑣 Energía máxima de entrada al inversor 𝑁𝑖𝑛𝑣 Número de inversores 𝑃𝑖𝑛𝑣 Energía de entrada al inversor 𝜂𝑖𝑛𝑣 Eficiencia del inversor 𝑃𝑎𝑐 Energía a la salida del inversor 𝑃𝑠,0 Energía consumida por el inversor debido a su funcionamiento 𝑃𝑎𝑐,0 Energía máxima a la salida del inversor 𝑃𝑛𝑡 Energía nocturna consumida por el inversor 𝑃𝑎𝑐𝑔𝑟𝑜𝑠𝑠 Energía bruta generada en CA 𝐿𝐴𝐶𝑤𝑖𝑟𝑖𝑛𝑔 Pérdidas en el cableado en CA 𝐿𝑡𝑟,𝑛𝑜 𝑙𝑜𝑎𝑑 Pérdidas en vacío del transformador 𝐿𝑡𝑟,𝑙𝑜𝑎𝑑 Pérdidas en carga del transformador 𝑃𝑎𝑐𝑛𝑒𝑡 Energía neta generada en CA Hardware y Software 6 6 definiciones de funciones como propiedades de una estructura. Esto hace posible la herencia en Julia y supone un gran beneficio para extender los paquetes. • Hecha para ML: Julia, a diferencia de Python, fue concebida para ser utilizada en estadística y ML. Además, el uso de álgebra lineal y el uso de las matrices es mucho más sencillo e intuitivo. Python puede ser optimizado para estas tareas, pero Julia ya viene optimizada de fábrica. • Gestor de paquetes: El gestor de paquetes de Julia (Pkg) es muy superior al de Python (Pip). Este viene con su propio REPL (Read Eval Print Loop) y permite construir, añadir, eliminar e instanciar los distintos paquetes. 7 3 CONSTRUCCIÓN DE UNA CÉLULA SOLAR La construcción de células solares se basa en el uso de un semiconductor que normalmente es silicio cristalino. La primera capa es un semiconductor tipo-N y se la conoce como capa emisora. La segunda capa, por el contrario, es un semiconductor tipo-P y se llama capa base. Como ambas capas están unidas se forma una unión P-N entre ellas [38] [39] [40]. Figura 3: Partes de una célula fotovoltaica Los átomos de silicio están compuestos por 14 electrones repartidos en 3 capas diferentes. La capa más exterior o de valencia tiene 4 electrones, por lo que un átomo intentará recombinarse con otros átomos cercanos para poder completar su capa de valencia. Si se aplica la suficiente energía al silicio puro, es posible arrancar uno de los electrones de la capa de valencia para que abandone el átomo. Estos electrones suelen llamase portadores libres y se moverán aleatoriamente por el silicio buscando un hueco con el que recombinarse. La capa emisora se dopa con fósforo, el cual tiene 5 electrones en la capa de valencia. De esta manera, se requiere menos energía para poder arrancar ese electrón adicional ya que no está unido a ningún átomo. Como resultado, se obtienen muchos más portadores libres que si se usara silicio puro. Un razonamiento análogo se utiliza para la capa base. En este caso se dopa con boro, el cual tiene solamente tres electrones de valencia. Nótese que los átomos de la capa N presentan ahora un hueco. Al unir estas dos capas, algunos de los electrones de la capa emisora se moverán a la capa base recombinándose con los huecos disponibles. De esta forma, se crea una zona de deplexión donde no hay ni electrones ni huecos libres. Debido al desplazamiento de los electrones, la capa emisora está cargada positivamente y la capa base negativamente, creándose, por tanto, un campo eléctrico entre ambas. Figura 4: Efecto de un fotón en una unión P-N Construcción de una célula solar 8 8 Los fotones provenientes del Sol están en constante interacción con la unión P-N. Si alguno de estos fotones presenta una energía superior a la de la banda prohibida, promocionará uno de los electrones de la banda de valencia a la banda de conducción, resultando a su vez en un hueco en la banda de valencia. Debido al campo eléctrico que existe entre ambas capas, los electrones se verán atraídos hacia la capa emisora mientras que los huecos se moverán a la capa base. Ahora que la concentración de electrones y huecos se ha incrementado, existe un potencial eléctrico que hará que los electrones fluyan si se conecta una carga. Tras completar su camino, los electrones se recombinarán con los huecos de la capa base. De esta manera, se consigue la generación de energía eléctrica mediante células solares. Normalmente la capa emisora es fina y fuertemente dopada mientras que la capa base es más ancha y débilmente dopada. Esto se hace así para aumentar la eficiencia de la célula. La zona de deplexión formada bajo esta configuración es mayor, haciendo posible una mayor generación de pares electrón-hueco y, por tanto, una mayor corriente eléctrica. 9 4 ÁNGULOS SOLARES Un requisito indispensable para el estudio de los cielos es determinar dónde se encuentran los objetos en él tanto temporalmente como espacialmente. Para tener en el primero se suele tomar un marco temporal de referencia para los objetos celestes. Con relación al segundo, los astrónomos han desarrollado varios sistemas de coordenadas para especificar la posición de los objetos en la bóveda celeste. La mayoría de estos sistemas se expresan en coordenadas esféricas utilizando una cuadrícula proyectada sobre la Esfera Celeste, es decir, una esfera imaginaria de gran radio centrada en la Tierra. Todos los objetos que pueden verse en la bóveda celeste pueden considerarse como situados en la superficie de esta esfera. La diferencia principal entre los sistemas de coordenadas es la elección del plano fundamental, el cual divide la Esfera Celeste en dos hemisferios iguales a lo largo de un gran círculo. En el caso particular de este proyecto, el cuerpo que hay que localizar en la Esfera Celeste es el Sol. Marco temporal El primer paso que se necesita es especificar el formato temporal que se va a utilizarla. Las horas estarán en formato de reloj de 24 horas comenzando a las 00:00 y con su fin a las 23:00. Además, la frecuencia temporal, aunque suele ser horaria, puede ser variable, tomándose como referencia el punto medio entre dos horas, excepto para la salida y la puesta del sol. Los archivos meteorológicos que se van a utilizar en las pruebas de validación son archivos de años típicos, lo que significa que las marcas temporales no tienen por qué utilizar siempre el mismo año. Por ejemplo, todas las marcas de tiempo de enero pueden ser de 1991, mientras que todas las marcas de tiempo de febrero pueden ser de 1995. Todos los archivos meteorológicos proporcionan un año típico no bisiesto (365 días), aunque los años de los meses de marzo a diciembre pueden corresponderse a años bisiestos. La fecha juliana del año jdoy (estrictamente llamada la fecha ordinal [41] [42]) es el número de días desde el 1 de enero y cambia para los años bisiestos. La Figura 5 muestra una tabla que relaciona el día del año con el día de cada mes. Sin embargo, el tiempo es relativo a la ubicación del observador. Debido a los husos horarios, la hora de un español será diferente a la de un australiano. Por lo tanto, es necesario utilizar una referencia horaria para todas las zonas horarias. Esta referencia es el UTC (Tiempo Universal Coordinado) [43], el cual es el principal estándar de tiempo por el que el mundo regula los relojes y la hora. Se puede calcular como la hora local del observador menos la zona horaria del observador. 𝑡𝑢𝑡𝑐=ℎ𝑟+𝑚𝑖𝑛 60−𝑡𝑧 (4-1) Esta ecuación puede dar lugar a valores fuera del rango [0,24), en cuyo caso se aplican las siguientes correcciones: {𝑡𝑢𝑡𝑐=𝑡𝑢𝑡𝑐+24 𝑗𝑑𝑜𝑦=𝑗𝑑𝑜𝑦−1 𝑖𝑓 𝑡𝑢𝑡𝑐<0 {𝑡𝑢𝑡𝑐=𝑡𝑢𝑡𝑐−24 𝑗𝑑𝑜𝑦=𝑗𝑑𝑜𝑦+1 𝑖𝑓 𝑡𝑢𝑡𝑐>24 (4-2) Hay que destacar que, en caso de que el nuevo jdoy sea cero, al año se le debe restar una unidad y, dependiendo de si el nuevo año es bisiesto o no, jdoy valdrá 366 o 365 respectivamente. Por el contrario, si el nuevo jdoy es mayor que 366 o 365 (dependiendo del año), debe ponerse a uno y el año debe ser incrementado en una unidad. Ángulos solares 10 10 Figura 5: Tabla con los días del año para años bisiestos y no bisiestos La Figura 6 muestra los husos horarios de cada país. Figura 6: Husos horarios La fecha actual no se va a expresar mediante el calendario gregoriano, sino con el calendario juliano, es decir, el número de días transcurridos desde el 1 de enero de 4713 a.C. [44].Este calendario está basado en el periodo 11 11 Software de simulación para la generación fotovoltaica juliano (7980 años), que es donde coinciden el ciclo solar de 28 años, el ciclo lunar de 19 años y la indicción (un impuesto) romana cada 15 años. La razón para utilizar este formato de fechas es que los calendarios civiles son poco manejables para medir largos intervalos de tiempo (los días del mes varían, bisiestos, etc.) [45] [46]. Dentro de este calendario, se ha marcado un punto de referencia o época llamando J2000, que es el número de días transcurridos desde el 1 de enero del 2000. Teniendo en cuenta que dicha referencia expresada en el calendario juliano es 2451545, la fecha actual se puede calcular como: 𝑗𝑢𝑙𝑖𝑎𝑛=2432916.5+365(𝑦𝑟−1949)+⌊𝑦𝑟−1949 4⌋+𝑗𝑑𝑜𝑦+𝑡𝑢𝑡𝑐 24−2451545 (4-3) Donde se ha tomado 1949 (2432916.5 en calendario juliano) arbitrariamente para poder calcular el número de días desde J2000. Marco especial Los ángulos solares más relevantes para la producción fotovoltaica terrestre son la elevación, el azimut y la declinación, cuyo cálculo se detallará en las subsecciones siguientes. Las ecuaciones presentadas en este apartado para la posición solar están basadas en el método de Michalsky [47], el cual a su vez está basado en el algoritmo del Almanaque Astronómico. Este algoritmo permite obtener una precisión de 0.01° desde 1950 hasta el 2050. Cabe destacar que se han incluido las modificaciones de Iqbal [48] para corregir los ángulos de azimut en localizaciones del hemisferio sur. 4.1.1 Sistema de coordenadas eclíptico El sistema de coordenadas eclíptico utiliza como plano fundamental la Eclíptica, es decir, la trayectoria que el Sol parece seguir a través del cielo a lo largo de un año [49]. También puede ser la proyección de la órbita de la Tierra sobre la Esfera Celeste. El Ecuador Celeste se obtiene siguiendo el mismo proceso, o sea, proyectando el ecuador de la Tierra sobre la Esfera Celeste. Los puntos de intersección entre la Eclíptica y el Ecuador Celeste se denominan Equinoccio Vernal o Equinoccio de Primavera y Equinoccio de Otoño. Las coordenadas de este sistema son la longitud y latitud eclípticas o celestes. La longitud eclíptica mide la distancia angular de un objeto a lo largo de la eclíptica desde la dirección primaria, es decir, la línea que une el centro de la Tierra y el equinoccio de primavera. La longitud eclíptica se mide en sentido positivo hacia el este en el plano fundamental. La latitud eclíptica mide la distancia angular de un objeto desde la eclíptica hacia el polo norte (positivo) o sur (negativo) de la eclíptica. Figura 7: Sistema de coordenadas eclíptico La latitud eclíptica puede asumirse como nula ya que tanto la Tierra como el Sol están en el plano eclíptico. Por tanto, sólo hay que calcular una coordenada, que es la longitud eclíptica. Para obtener esta longitud, es necesario calcular la anomalía media y la longitud media. Ángulos solares 12 12 Longitud media en grados (0≤𝑚𝑛𝑙𝑜𝑛𝑔<360): 𝑚𝑛𝑙𝑜𝑛𝑔=280.46+0.9856474 𝑗𝑢𝑙𝑖𝑎𝑛 (4-4) La anomalía media en radianes (0≤𝑚𝑛𝑎𝑛𝑜𝑚<2𝜋): 𝑚𝑛𝑎𝑛𝑜𝑚= 𝜋 180(357.528+0.9856003 𝑗𝑢𝑙𝑖𝑎𝑛) (4-5) La longitud eclíptica en radianes (0≤𝑒𝑐𝑙𝑜𝑛𝑔<2𝜋) es, por tanto: 𝑒𝑐𝑙𝑜𝑛𝑔= 𝜋 180(𝑚𝑛𝑙𝑜𝑛𝑔+1.915sin(𝑚𝑛𝑎𝑛𝑜𝑚)+0.02sin(2 𝑚𝑛𝑎𝑛𝑜𝑚)) (4-6) La oblicuidad de la eclíptica, el ángulo que forman los planos de la eclíptica y el ecuatorial celeste, será necesaria para posteriores transformaciones y viene dada por: 𝑜𝑏𝑙𝑒𝑞= 𝜋 180(23.439−0.0000004 𝑗𝑢𝑙𝑖𝑎𝑛) (4-7) 4.1.2 Sistema de coordenadas ecuatorial El sistema de coordenadas ecuatoriales suele ser el preferido por los astrónomos para hacer un seguimiento de la posición de los cuerpos en el cielo. Este sistema puede expresar la posición de una estrella o galaxia con independencia de dónde esté el observador en la Tierra [50]. Sin embargo, debido a que este sistema se ve afectado por los movimientos de precesión y nutación, se ha comenzado el estudio a partir de las coordenadas eclípticas. La precesión es el movimiento del eje de rotación de un cuerpo astronómico por el que el eje traza lentamente un cono. En el caso de la Tierra, este tipo de precesión se conoce también como "precesión de los equinoccios". La Tierra atraviesa un ciclo de precesión de este tipo en un periodo de aproximadamente 26.000 años, durante el cual las posiciones de los astros cambian lentamente tanto en las coordenadas ecuatoriales como en la longitud de la eclíptica. La nutación es un pequeño cambio en el ángulo de inclinación de la Tierra con respecto al Sol que cambia la ubicación de los principales círculos de latitud. Este movimiento tiene un periodo de 18,6 años. Figura 8: Rotación (R), Precesión (P) y Nutación (N) Ambos movimientos se representan en la Figura 8. En este proyecto se utilizan las coordenadas eclípticas como punto de partida para tener en cuenta estos dos movimientos adicionales de la Tierra. Sin embargo, dado que estos movimientos cambian la posición del equinoccio de primavera en el cielo a lo largo de los siglos, se 13 13 Software de simulación para la generación fotovoltaica necesita un punto temporal de referencia o época [51]. La época estándar que se utiliza actualmente es la que se comentó anteriormente en el apartado 0, es decir, la "época juliana J2000.0", la cual se define como el número de días transcurridos desde el mediodía del 1 de enero de 2000. El plano fundamental es el Ecuador Celeste y las coordenadas de este sistema son la ascensión recta y la declinación. Además, la dirección primaria de este sistema es el ecuador. La declinación es el ángulo latitudinal y mide la distancia angular de un objeto perpendicular al Ecuador Celeste, positiva hacia el norte y negativa hacia el sur. La ascensión recta es la coordenada longitudinal y mide la distancia angular de un objeto hacia el este a lo largo del Ecuador Celeste desde el Equinoccio Vernal. Figura 9: Sistema de coordenadas ecuatorial Los parámetros calculados para el sistema eclíptico permiten calcular las coordenadas ecuatoriales de manera sencilla. La ascensión recta en radianes es: 𝑟𝑎= { atan(𝑐𝑜𝑠(𝑜𝑏𝑙𝑒𝑞)𝑠𝑖𝑛(𝑒𝑐𝑙𝑜𝑛𝑔) 𝑐𝑜𝑠(𝑒𝑐𝑙𝑜𝑛𝑔))+𝜋 𝑖𝑓cos(𝑒𝑐𝑙𝑜𝑛𝑔)<0 atan(𝑐𝑜𝑠(𝑜𝑏𝑙𝑒𝑞)𝑠𝑖𝑛(𝑒𝑐𝑙𝑜𝑛𝑔) 𝑐𝑜𝑠(𝑒𝑐𝑙𝑜𝑛𝑔))+2𝜋 𝑖𝑓cos(𝑜𝑏𝑙𝑒𝑞)sin(𝑒𝑐𝑙𝑜𝑛𝑔)<0 (4-8) El ángulo de la declinación solar en radianes viene dado por: 𝛿=asin(𝑠𝑖𝑛(𝑜𝑏𝑙𝑒𝑞)𝑠𝑖𝑛(𝑒𝑐𝑙𝑜𝑛𝑔)) (4-9) Sin embargo, en lugar de la ascensión recta, se suele utilizar otro ángulo: el ángulo horario. Para definir este nuevo ángulo, hay que explicar en primer lugar el tiempo sideral. El tiempo sideral significa "tiempo de las estrellas" y es diferente del tiempo utilizado en nuestra vida cotidiana, que es el tiempo solar. La unidad fundamental del Tiempo Solar es un Día, que se define como el tiempo que tarda el Sol en recorrer 360 grados alrededor del cielo debido a la rotación de la Tierra. Sin embargo, hay un problema con el Tiempo Solar: la Tierra no gira realmente 360 grados en un Día Solar. La Tierra está en órbita alrededor del Sol y, en el transcurso de un día, se mueve aproximadamente un grado a lo largo de su órbita. Por lo tanto, la Tierra tiene que girar 361 grados para que parezca que el Sol ha viajado 360 grados alrededor del Cielo. Este periodo de rotación se denomina Día Sideral y es 4 minutos más corto que un Día Solar por término medio debido al grado extra [52]. La ascensión recta de un objeto indica el tiempo sideral en el que transitará por el meridiano local, un círculo imaginario en la esfera celeste que pasa por los polos celestes y un punto directamente por encima de la cabeza de un observador. Este punto directamente por encima de la cabeza se llama cénit y el punto a 180º de Ángulos solares 14 14 distancia se llama nadir. Figura 10: Sistema de coordenadas ecuatorial (utilizando el ángulo horario) Por otro lado, el Ángulo Horario del objeto indica el tiempo que ha pasado desde que el objeto estaba en el meridiano local. También es la distancia angular entre el objeto y el meridiano, medida en horas. Los ángulos horarios negativos indican el tiempo que falta para el próximo tránsito por el meridiano local y los positivos el tiempo transcurrido desde que el objeto estuvo en el meridiano local. Por supuesto, un Ángulo Horario de cero significa que el objeto está actualmente en el meridiano local. La Figura 10 permite ilustrar este sistema de coordenadas. Para calcular este ángulo, se necesita la ascensión recta y el tiempo sideral medio local. El tiempo sideral medio de Greenwich en horas se puede obtener como 𝑔𝑚𝑠𝑡=6.697375+0.0657098242 𝑗𝑢𝑙𝑖𝑎𝑛+𝑡𝑢𝑡𝑐 (4-10) La hora sideral media local no es más que un desplazamiento de la hora sideral media de Greenwich hasta el meridiano local: 𝑙𝑚𝑠𝑡=𝑔𝑚𝑠𝑡+𝑙𝑜𝑛𝑔 15 (4-11) La longitud anterior debe expresarse en grados y la hora sideral media local obtenida se mide en horas. El ángulo horario en radianes (−𝜋<𝐻𝐴<𝜋) viene dado por: 𝑏=15 𝜋 180𝑙𝑚𝑠𝑡−𝑟𝑎 (4-12) 𝐻𝐴={𝑏 𝑖𝑓 −𝜋<𝑏<𝜋 𝑏+2𝜋 𝑖𝑓 𝑏≤−𝜋 𝑏−2𝜋 𝑖𝑓 𝑏≥𝜋 (4-13) 4.1.3 Sistema de coordenadas horizontal El sistema de coordenadas ecuatoriales utiliza el horizonte local del observador como plano fundamental [53]. Esto divide convenientemente el cielo en el hemisferio superior que se puede ver y el hemisferio inferior que no se puede ver. Como ya se ha mencionado, el polo del hemisferio superior se llama cenit y el polo del hemisferio inferior se llama nadir. El ángulo de un objeto por encima o por debajo del horizonte se llama elevación y el ángulo alrededor del horizonte (medido desde el punto Norte, hacia el Este) se llama acimut. Como este sistema de coordenadas 15 15 Software de simulación para la generación fotovoltaica está fijado a la Tierra, la elevación y el acimut de un objeto cambian con el tiempo. Además, como el sistema horizontal está definido por el horizonte local de cada uno, el mismo objeto visto desde diferentes lugares de la Tierra en el mismo momento tendrá diferentes valores de Elevación y Acimut. Figura 11: Sistema de coordenadas horizontal El Sistema de coordenadas horizontal puede ser calculado a partir de las coordenadas del sistema ecuatorial y haciendo uso de algunas relaciones angulares. El ángulo de elevación solar en radianes es: 𝑎=sin(𝛿)sin(𝜋 180𝑙𝑎𝑡𝑖𝑡𝑢𝑑𝑒)+cos(𝛿)cos(𝜋 180𝑙𝑎𝑡𝑖𝑡𝑢𝑑𝑒)cos(𝐻𝐴) (4-14) 𝛼0= { asin(𝑎) 𝑖𝑓−1≤𝑎≤1 𝜋2 𝑖𝑓 𝑎>1 −𝜋2 𝑖𝑓 𝑎<−1 (4-15) El ángulo azimutal 𝛾 en radianes: 𝑎=sin(𝛼0)sin(𝜋 180𝑙𝑎𝑡𝑖𝑡𝑢𝑑𝑒)−sin(𝛿) cos(𝛼0)cos(𝜋 180𝑙𝑎𝑡𝑖𝑡𝑢𝑑𝑒) (4-16) 𝑏={acos(𝑎) 𝑖𝑓−1≤𝑎<1 𝜋 𝑖𝑓cos(𝛼0)=0 𝑜𝑟 𝑎<−1 0 𝑖𝑓 𝑎>1 (4-17) 𝛾={𝑏 𝑖𝑓 𝐻𝐴<−𝜋 𝜋−𝑏 𝑖𝑓−𝜋≤𝐻𝐴≤0 𝑜𝑟 𝐻𝐴≥0 𝜋+𝑏 𝑖𝑓 0<𝐻𝐴<𝜋 (4-18) Sin embargo, cuando se observa cerca del horizonte y durante el amanecer y el atardecer, es necesario tener en cuenta la extensión angular del Sol, pues se utiliza el limbo superior del disco solar como punto de referencia. En estos casos, la refracción no puede ser obviada y es necesario corregir el ángulo de elevación en consecuencia. 𝛼0𝑑=180 𝜋𝛼0 (4-19) 𝑟={3.51561(0.1594+0.0196𝛼0𝑑+0.00002𝛼0𝑑 2 1+0.505𝛼0𝑑+0.0845𝛼0𝑑 2) 𝑖𝑓 𝛼0𝑑>−0.56 0.56 𝑖𝑓 𝛼0𝑑≤−0.56 (4-20) Análisis de la irradiancia 22 22 Las variaciones de la TSI son de aproximadamente un 0.1% en cada ciclo solar (11 años) debido a las manchas solares. El nuevo valor aceptado es de 1360.8±0.5 𝑊/𝑚2, el cual se midió durante el mínimo solar ocurrido en 2008 [63]. El valor de TSI que se estableció en 1990 y que estuvo en vigor hasta la actualización del mismo en 2008 era 1365.4±1.3 𝑊/𝑚2. Es importante conocer el valor antiguo, ya que se sigue encontrando en la literatura. El valor de la constante solar se puede obtener de diferentes formas. Una primera aproximación es calcular la integral de la irradiancia espectral sobre todo el espectro. Otra forma es aplicando la ley de Stephan-Boltzmann para calcular la potencia emitida por el Sol y multiplicarla por la relación de áreas entre la superficie del Sol y la de una esfera de radio una unidad astronómica. Sin embargo, los métodos anteriores son solo meras aproximaciones. Actualmente, la constante solar se obtiene gracias a las mediciones de los satélites. Por ejemplo, la NASA lleva a cabo mediciones denominadas estado del arte en el satélite de su misión Solar Radiation and Climate Experiment (SORCE) [64]. La irradiancia extraterrestre es la cantidad de irradiancia total que llega al punto más exterior de la atmósfera. Aunque pueda parecer que dicha irradiancia es el valor de la constante solar, esto no es del todo cierto. Debido a que la órbita terrestre es elíptica, la distancia que hay entre el Sol y la Tierra varía a lo largo del año. Esto resulta en una variación de la irradiancia extraterrestre de aproximadamente un ±3.3% que se puede visualizar en la Figura 17 [65]. Figura 17: Variación de la constante solar a lo largo del año Matemáticamente, la irradiancia extraterrestre en un plano normal a los rayos solares contando con la variación estacional viene dada por [65]: 𝐺𝑜𝑛=𝐺𝑠𝑐𝐹𝑛=𝐺𝑠𝑐(1+0.033cos360𝑛 365) (5-2) Donde 𝐺𝑠𝑐 es la constante solar y 𝑛 es el día del año según la Figura 5. Si se necesita una mayor precisión en la variación estacional: 𝐹𝑛=1.000110+0.034221cos(Γ)+0.001280sin(𝛤)+0.000719cos(2Γ) +0.000077sin(2𝛤) (5-3) 23 23 Software de simulación para la generación fotovoltaica Con Γ=2𝜋(𝑛−1) 365 Para obtener la irradiancia extraterrestre en un plano horizontal, solo es necesario tener en cuenta el ángulo cenital: 𝐺𝑜=𝐺𝑜𝑛cos𝜃𝑧 (5-4) Descomposición de la irradiancia extraterrestre La irradiancia extraterrestre debe atravesar toda la atmósfera para poder llegar a la superficie terrestre. Durante dicho proceso, la irradiancia se atenuará debido a la interacción de los diferentes componentes de la atmósfera con determinadas bandas del espectro. Una parte de la irradiancia total que llega a la superficie terrestre es llamada irradiancia difusa y proviene de todas las direcciones de la bóveda celeste. La otra parte se llama irradiancia directa y alcanza al observador en una línea recta desde el Sol. Existe un tercer tipo de irradiancia llamada reflejada o albedo que es la que llega al observador tras reflejarse en el suelo o en el agua. La suma de las tres es conocida como irradiancia global. En la Figura 18 se puede observar este proceso de descomposición. Figura 18: Descomposición de la irradiancia extraterrestre La reducción del espectro de la irradiancia se debe a los diferentes componentes de la atmósfera debido a absorción y dispersión [56]. El efecto que produce cada uno de dichos componentes es: • Ozono: absorbe la irradiancia en el espectro de 200-300 nm y ligeramente en el espectro visible • Vapor de agua, dióxido de carbono y oxígeno: absorbe la irradiancia cercana a la banda infrarroja • Dispersión Rayleigh y aerosoles: reduce la irradiancia en el espectro visible por medio de dispersión. La Figura 19 muestra cada uno de los efectos anteriores sobre el espectro de la irradiancia extraterrestre. Esta última es la más importante en aplicaciones fotovoltaicas. La distinción entre estas componentes es de vital importancia para calcular la potencia de salida y se pueden calcular mediante diversos sensores. Como se puede observar, la cantidad de parámetros atmosféricos que intervienen en la descomposición de la irradiancia es realmente alto. Es por ello por lo que normalmente las distintas componentes de la irradiancia se obtienen mediante sensores. Análisis de la irradiancia 24 24 Figura 19: Efecto de la atmósfera en la irradiancia espectral Dependiendo de la componente que se desea medir, se utiliza un sensor u otro [66]. Por ejemplo, para medir la irradiancia global, ya sea horizontal (GHI) o sobre cualquier plano, se utiliza un piranómetro. Si se busca medir la irradiancia directa normal (DNI) se utiliza un dispositivo llamado pirheliómetro montado sobre un seguidor solar. La irradiancia difusa horizontal (DHI) necesita de ambos sensores. La relación entre las distintas componentes es la siguiente: 𝐺𝐻𝐼=𝐷𝑁𝐼cos𝜃𝑧+𝐷𝐻𝐼 (5-5) Aunque lo ideal es contar con sensores que permitan conocer todas las componentes de la irradiancia, normalmente el único dispositivo disponible es el piranómetro. En estos casos se pueden usar modelos matemáticos como el Engerer2 [67] combinado con el REST2 [68] [69] para poder descomponer la GHI en las componentes directa y difusa. Ambos modelos están detallados en el Apéndice II: Modelos REST2 y Engerer2. Los modelos anteriores se han implementado para poder dar la opción al usuario de usar solamente la componente GHI. Sin embargo, ambos métodos requieren de una gran cantidad de parámetros que deben ser obtenidos de distintas fuentes, por ejemplo, el AOD (Aerosol Optical Depth) se puede obtener mediante imágenes por satélites o de AERONET (Aerosol Robotic Network) [70]. Con el objetivo de mantener el modelo lo más simple posible, se van a utilizar valores típicos por defecto, los cuales se pueden consultar en la Tabla 1. Parámetro Rango Valor típico Presión 𝑝∈(300,100) 𝑚𝑏𝑎𝑟 𝑝=1013.25 𝑚𝑏𝑎𝑟 Ozono 𝑢𝑜∈(0,0.6) 𝑎𝑡𝑚−𝑐𝑚 𝑢𝑜=0.3 𝑎𝑡𝑚−𝑐𝑚 Dióxido de nitrógeno 𝑢𝑛∈(0,0.03) 𝑎𝑡𝑚−𝑐𝑚 𝑢𝑛=2∙10−4 𝑎𝑡𝑚−𝑐𝑚 Precipitaciones 𝑤∈(0,10) 𝑐𝑚 𝑤=1.5 𝑐𝑚 Coeficientes de turbiedad de Ångström 𝛼1,𝛼2∈(0,2.5) 𝛽∈(0,1.1) 𝛼1=𝛼2=1.3 𝛽=0.1 25 25 Software de simulación para la generación fotovoltaica Tabla 1: Parámetros típicos del modelo REST2 Irradiancia sobre un plano inclinado Las irradiancias medidas por los sensores del apartado anterior hacen referencia a las distintas componentes medidas sobre la superficie terrestre sobre un plano horizontal (GHI y DHI) o perpendicular (DNI) a los rayos solares. Sin embargo, las placas solares nunca están dispuestos en estos ángulos, por lo que es preciso calcular las componentes anteriores para un ángulo genérico. 5.1.1 Ángulo de incidencia La descripción de un plano con una inclinación y orientación arbitraria (refiriéndose al panel) se puede hacer usando dos ángulos: azimut 𝛾𝑠∈[−𝜋,𝜋] e inclinación 𝛽𝑠∈[0,𝜋/2]. La Figura 20 describe perfectamente las relaciones trigonométricas para calcular dicho ángulo, el cual viene dado por: 𝑎=cos(𝜃)=sin(𝜃𝑧)cos(𝛾−𝛾𝑠)sin(𝛽𝑠)+cos(𝜃𝑧)cos(𝛽𝑠) (5-6) 𝜃={𝜋 𝑖𝑓 𝑎<−1 0 𝑖𝑓 𝑎>1 𝑜 𝑠𝑒𝑔𝑢𝑖𝑚𝑖𝑒𝑛𝑡𝑜 𝑑𝑒 𝑑𝑜𝑠 𝑒𝑗𝑒𝑠 arccos𝑎 𝑖𝑓−1≤𝑎≤1 (5-7) Donde 𝜃 es el ángulo de incidencia, 𝜃𝑧 el ángulo cenital, 𝛾 el ángulo acimutal, 𝛾𝑠 el ángulo acimutal del plano inclinado y 𝛽𝑠 la inclinación del propio plano Figura 20: Ángulo de incidencia sobre un plano inclinado Es necesario destacar que algunos colectores solares poseen la capacidad de realizar el seguimiento del Sol, ya sea mediante uno o dos ejes. Si el seguimiento se realiza acimutalmente, 𝛾𝑠=𝛾. Si, por el contrario, el seguimiento se realiza en dos ejes, 𝛾𝑠=𝛾 y 𝛽𝑠=𝜃𝑧 5.1.2 Irradiancia directa La relación entre la irradiancia directa en un plano inclinado y en una superficie horizontal es básicamente un factor geométrico puro que depende de los ángulos de incidencia: 𝑅𝑏=max(0,𝑐𝑜𝑠𝜃 𝑐𝑜𝑠𝜃𝑧) (5-8) Análisis de la irradiancia 26 26 Por tanto, la componente directa se puede calcular como: 𝐺𝑏𝑖=𝑅𝑏𝐷𝑁𝐼 (5-9) 5.1.3 Irradiancia debida a la reflexión La irradiancia reflejada es aquella que alcanza la superficie inclinada después de haberse reflejado en el suelo. La irradiancia reflejada es función de la DNI, el ángulo cenital, la DHI y la reflectividad del suelo, también conocida como albedo [71]: 𝐺𝑔𝑖=𝜌(𝐷𝑁𝐼cos𝜃𝑧+𝐷𝐻𝐼)1−cos𝛽 2 (5-10) Donde 𝜌 es el albedo. 5.1.4 Irradiancia difusa La irradiancia difusa es aquella que ha sido dispersada por la atmósfera terrestre antes de alcanzar la superficie. Esta irradiancia, a su vez, está compuesta por otras tres componentes: isotrópica, circunsolar y brillante horizontal. La irradiancia circunsolar es una región alrededor del disco solar en la que esta presenta pequeñas variaciones comparada con la directa. La irradiancia brillante horizontal es la que hace que el cielo parezca más claro en el horizonte y toma especial importancia en áreas nevadas. La irradiancia isotrópica es la parte dominante de la difusa. Un esquema de las tres se puede observar en la Figura 21. Figura 21: Componentes de la irradiancia difusa Una sencilla aproximación es suponer que la parte circunsolar y la brillante horizontal son despreciables con respecto a la isotrópica [71]. La irradiancia difusa que llega al colector es, por tanto, solo una fracción de la irradiancia difusa que alcanza a una superficie horizontal: 𝐺𝑑𝑖=𝐺𝑑𝑖 𝑖𝑠𝑜=(1+cos𝛽 2)𝐷𝐻𝐼 (5-11) Por supuesto, existen modelos matemáticos más complejos que tienen en cuenta todas las partes de la irradiancia difusa. Un ejemplo de ellos son el HDKR [65] o el modelo de Pérez [72] [73]. 27 27 Software de simulación para la generación fotovoltaica Sombras y suciedad La irradiancia efectiva en una superficie inclinada se puede definir como la irradiancia incidente menos las pérdidas debidas a las sombras de objetos cercanos, el auto-sombreado y a la suciedad. Esta es, por tanto, la irradiancia que finalmente llegaría al módulo fotovoltaico. Las sombras causadas por objetos cercanos se traducen en una reducción de la irradiancia que incide sobre una superficie inclinada. Ejemplos de ello pueden ser edificios, árboles, colinas, etc. Esto reduce tanto la irradiancia directa como la difusa. Sin embargo, cuando se habla de parques solares, estas sombras no deben considerarse en absoluto asumiendo que, si el parque está bien diseñado, no existirán objetos cercanos que puedan proyectar sombras sobre los paneles. El auto-sombreado es el sombreado de los módulos fotovoltaicos de una fila sobre los módulos de las filas vecinas. Esto, como en el caso anterior, provocaría una reducción de la irradiancia que llega al módulo. Afortunadamente, no es necesario tener en cuenta estas pérdidas suponiendo de nuevo que la central está bien diseñada, ya que las filas estarán lo suficientemente alejadas como para evitar este efecto. Esta distancia mínima viene ilustrada en la Figura 22 y puede calcularse fácilmente según la siguiente ecuación [74]: 𝑑=𝑑1+𝑑2= 𝑧 tan (ℎ0)+𝑧 tan (𝛽𝑠) (5-12) Figura 22: Sombreado de los paneles La variable z se puede calcular como 𝑊𝑠𝑖𝑛(𝛽𝑠), donde 𝑊 es el ancho del módulo. Por último, ℎ0 es el mínimo valor del ángulo de elevación solar a lo largo del año, el cual ocurre el 21 de diciembre. Dicho valor viene dado por: ℎ0=(90−𝜙)+𝛿 (5-13) En el peor de los casos (21 de diciembre) el valor de 𝛿 es −23.45°. Por último, pero no menos importante, las pérdidas por suciedad causan una reducción de la irradiancia incidente debido a la presencia polvo y tierra en la superficie de los paneles. Estas pérdidas se expresarán en términos de porcentaje: 𝑂=1− 𝐿 100 (5-14) Por tanto, la irradiancia global que alcanza la superficie del módulo tras tener en cuenta las pérdidas por sombras y suciedad se puede calcular como: 𝐺=𝐺𝑖𝑂 (5-15) Análisis de la irradiancia 28 28 Modificador del ángulo de incidencia y de la masa de aire Una vez se han obtenido las componentes de la irradiancia que inciden en un plano inclinado, o lo que es lo mismo, en un colector solar, es necesario calcular cuánto de dicha irradiancia es realmente absorbida. La transmitancia es función de la irradiancia incidente y del grosor, índice de refracción y coeficiente de extinción del material [65] [75]. 𝑛=1.526 (índice de refracción del cristal) 𝐿=0.002 (grosor de la cubierta de cristal) 𝐾=4 (constante de proporcionalidad) (5-16) Utilizando la ley de Snell, el ángulo de refracción viene dado por: 𝜃𝑟=arcsin(1𝑛𝑠𝑖𝑛(𝜃𝑖)) (5-17) Donde 𝜃𝑖 es el ángulo de incidencia y 𝑛 el índice de refracción. Los ángulos de incidencia de las diferentes componentes de la irradiancia vienen dados por 𝜃𝑏=𝜃 (Á𝑛𝑔𝑢𝑙𝑜 𝑑𝑒 𝑖𝑛𝑐𝑖𝑑𝑒𝑛𝑐𝑖𝑎) 𝜃𝑑=59.7−0.1388𝛽𝑠+0.001497𝛽𝑠2 𝜃𝑔=90−0.5788𝛽𝑠+0.002693𝛽𝑠2 (5-18) La transmitancia como función del ángulo de incidencia viene dado por: 𝜏(𝜃)=𝑒−𝐾𝐿 cos(𝜃𝑟)(1−12(sin2(𝜃𝑟−𝜃) sin2(𝜃𝑟+𝜃)+tan2(𝜃𝑟−𝜃) tan2(𝜃𝑟+𝜃))) (5-19) Donde 𝐾 es el coeficiente de extinción y 𝐿 el grosor del material. En caso de que el ángulo de incidencia sea 0°, la ecuación anterior se reduce a: 𝜏(0)=lim 𝜃→0𝜏(𝜃)=𝑒−𝐾𝐿(1−(1−𝑛 1+𝑛)2) (5-20) El modificador del ángulo de incidencia se puede calcular como el ratio entre las transmitancias para los ángulos 𝜃 y 0°: 𝐾𝑏=𝜏(𝜃𝑏) 𝜏(0) 𝐾𝑑=𝜏(𝜃𝑑) 𝜏(0) 𝐾𝑔=𝜏(𝜃𝑔) 𝜏(0) (5-21) La irradiancia absorbida por la célula fotovoltaica es, por tanto: 𝐺𝑎𝑏𝑠=𝐺𝑏𝐾𝑏+𝐺𝑑𝐾𝑑+𝐺𝑔𝐾𝑔 (5-22) Además, el producto transmitancia-absorbancia necesario para calcular la temperatura de la célula es: 29 29 Software de simulación para la generación fotovoltaica 𝜏𝛼=0.9𝐺𝑎𝑏𝑠 𝐺 (5-23) El modificador de la masa de aire recoge los efectos de la masa de aire en la distribución espectral y se puede calcular como: 𝑀=𝑎0+𝑎1𝐴𝑀+𝑎2𝐴𝑀2+𝑎3𝐴𝑀3+𝑎4𝐴𝑀4 (5-24) Donde 𝐴𝑀 ha sido corregida con un factor exponencial para tener en cuenta la altura sobre el nivel del mar del observador: 𝐴𝑀=(cos(𝜃𝑧)+0.5057(96.08−𝜃𝑧)−1.634)−1𝑒−0.0001184𝐴 (5-25) Los coeficientes dependen del tipo de material [76]. Por ejemplo, para el silicio policristalino que es bastante común son: 𝑀=0.918093+0.086257𝑚−0.024459𝑚2+0.002816𝑚3−0.000126𝑚4 (5-26) Finalmente, la irradiancia efectiva que absorbe la célula solar es: 𝐺𝑒𝑓𝑓=𝑀𝐺𝑎𝑏𝑠 (5-27) Análisis de la irradiancia 30 30 31 6 PRINCIPIO DE FUNCIONAMIENTO Temperatura de la célula fotovoltaica La temperatura de operación de una célula fotovoltaica se determina mediante un balance energético [65]. Parte de la energía solar que es absorbida por una célula se convierte en energía térmica y la otra parte en energía eléctrica que se extrae de la célula gracias a un circuito externo. La temperatura de funcionamiento nominal de la célula (NOCT) se define como aquella que se alcanza con la célula montada de manera normal bajo condiciones de trabajo de una irradiancia de 800 𝑊/𝑚2, la velocidad del viento es 1 𝑚/𝑠 y la temperatura ambiente es de 20 °𝐶. La velocidad del viento ha de ajustarse dependiendo a la altura que se coloquen los paneles solares del suelo: 𝑣𝑤,𝑎𝑑𝑗={0.51𝑣𝑤 3.3𝑚 𝑜 𝑚𝑒𝑛𝑜𝑠 0.61𝑣𝑤 6.6𝑚 𝑜 𝑚á𝑠 (6-1) Del mismo modo, la temperatura NOCT también ha de ajustarse dependiendo del soporte de montaje de los paneles: 𝑇𝑁𝑂𝐶𝑇,𝑎𝑑𝑗= { 𝑇𝑁𝑂𝐶𝑇 𝑖𝑛𝑡𝑒𝑔𝑟𝑎𝑑𝑜,𝑚𝑎𝑦𝑜𝑟 𝑑𝑒 3.5′𝑜 𝑚𝑜𝑛𝑡𝑎𝑑𝑜 𝑒𝑛 𝑠𝑢𝑒𝑙𝑜 𝑇𝑁𝑂𝐶𝑇+2 𝑑𝑒 2.5′ 𝑎 3.5′ 𝑇𝑁𝑂𝐶𝑇+6 𝑑𝑒 1.5′𝑎 2.5′ 𝑇𝑁𝑂𝐶𝑇+11 𝑑𝑒 0.5′ 𝑎 1.5′ 𝑇𝑁𝑂𝐶𝑇+18 𝑚𝑒𝑛𝑜𝑠 𝑑𝑒 0.5′ (6-2) La temperatura de la célula en grados Celsius viene dada por: 𝑇𝑐=𝑇𝑎+𝐺𝑒𝑓𝑓 800(𝑇𝑁𝑂𝐶𝑇,𝑎𝑑𝑗−20)(1−𝜂𝑟𝑒𝑓 𝜏𝛼)( 9.5 5.7+3.8𝑣𝑤,𝑎𝑑𝑗) (6-3) Donde 𝜏𝛼 es el producto transmitancia-absorbancia y 𝜂𝑟𝑒𝑓 es la eficiencia de la célula trabajando en el punto de energía máxima. Este último parámetro puede calcularse como: 𝜂𝑟𝑒𝑓=𝐼𝑚𝑉𝑚 𝐺𝑟𝑒𝑓𝐴𝑚 (6-4) En la ecuación anterior 𝐼𝑚 y 𝑉𝑚 son la corriente y la tensión en el punto de máxima energía bajo condiciones de referencia (STC). Dichas condiciones son concretamente una irradiancia de 1000 𝑊/𝑚2 y una temperatura de 25°𝐶. 𝐴𝑚 es área de la célula. Modelo fotovoltaico El comportamiento de células fotovoltaicas puede modelarse con un circuito como el de la Figura 23 [65] [75]. Dicho circuito está compuesto por una resistencia en serie y un diodo en paralelo con una resistencia. Dicho circuito puede usarse para modelar una sola celda solar o un panel fotovoltaico compuesto por varias de estas. Matriz de módulos fotovoltaicos 38 38 La energía calculada anteriormente no tiene en cuenta todas las pérdidas que puede presentar un sistema real. Algunas de estas pérdidas son [85] [86] [87] [88]: • Pérdidas por temperatura: estas pérdidas son modeladas por el coeficiente térmico de energía y están en torno al 0,4-0,5% por cada grado que aumenta la temperatura, teniendo como referencia la temperatura a STC (25ºC). Las pérdidas por este efecto son mayores en verano. Afortunadamente, esta pérdida ya ha sido considerada en el modelo • Pérdidas nominales (nameplate): los módulos fotovoltaicos se fabrican gracias a un proceso industrial y, por tanto, va a haber pequeñas diferencias entre ellos. Esto significa que la energía va a ser ligeramente diferente si se comparan dos modelos del mismo módulo. Este dato lo da el fabricante como tolerancia. • Pérdidas por desajuste (mismatch): estas pérdidas se producen por la interconexión de módulos solares que no tienen propiedades idénticas (como se ha explicado antes) o que experimentan condiciones diferentes entre sí. El módulo en serie con menor energía restringirá la corriente y lo mismo ocurre con la tensión en los módulos en paralelo. Estas pérdidas se sitúan en torno al 1-4% y sólo pueden determinarse experimentalmente. • Pérdidas por sombreado: las sombras en los módulos producen pérdidas de energía debido a la menor irradiación incidente. Estas pérdidas ya han sido consideradas. • Pérdidas por suciedad: estas pérdidas dependen de la ubicación del parque solar y de la frecuencia de la lluvia, entre otros factores. Estas pérdidas ya han sido consideradas. • Pérdidas angulares: la energía nominal de un módulo se da en condiciones STC siempre que los rayos solares lleguen a la superficie en una dirección normal al mismo. Si los rayos no llegan perpendicularmente al plano del módulo, la irradiancia absorbida por el módulo será menor y, por tanto, también lo será la energía. Estas pérdidas ya se han tenido en cuenta en el modelo. • Pérdidas en el cableado: en las partes de corriente continua de la instalación, habrá pérdidas debidas a la resistividad del conductor. El factor de pérdidas es, por tanto: 𝐹𝑑𝑐=(1−𝐿𝑚𝑖𝑠𝑚𝑎𝑡𝑐ℎ)(1−𝐿𝑛𝑎𝑚𝑒𝑝𝑙𝑎𝑡𝑒)(1−𝐿𝑑𝑐𝑤𝑖𝑟𝑖𝑛𝑔) (7-2) La energía neta en CC de cada submatriz se calcula como la energía bruta a la que se le sustraen las pérdidas anteriores: 𝑃𝑑𝑐,𝑛𝑒𝑡,𝑛=𝑃𝑑𝑐,𝑔𝑟𝑜𝑠𝑠,𝑛𝐹𝑑𝑐 (7-3) La energía neta de salida en CC de la matriz no es más que la suma de las energías netas de cada una de las submatrices: 𝑃𝑑𝑐,𝑛𝑒𝑡= ∑ 𝑃𝑑𝑐,𝑛𝑒𝑡,𝑛 𝑁𝑠𝑢𝑏𝑎𝑟𝑟𝑎𝑦𝑠 𝑛=1 (7-4) Salida en CA La energía neta en CC debe convertirse en CA mediante el uso de inversores para poder ser transmitida por la red de transporte o distribución hasta los consumidores finales. Por tanto, lo primero que es necesario hacer es calcular el número de inversores necesarios para manejar la energía obtenida de la matriz fotovoltaica. Las ecuaciones que definen al inversor vienen dadas por [89], particularizando para el caso en el que los datos se toman de la hoja de especificaciones del fabricante, es decir, haciendo los coeficientes empíricos 𝐶𝑖 nulos. 39 39 Software de simulación para la generación fotovoltaica 𝑁𝑖𝑛𝑣=⌈𝑃𝑑𝑐𝑛𝑒𝑡 𝑃𝑚𝑎𝑥𝑖𝑛𝑣⌉ (7-5) El número de inversores viene dado por la relación entre la energía total existente y la energía máxima que el inversor puede manejar (redondeado al entero siguiente). Con el número de inversores, se puede distribuir la energía en CC de forma equitativa entre todos ellos de la siguiente manera: 𝑃𝑖𝑛𝑣=𝑃𝑑𝑐𝑛𝑒𝑡 𝑁𝑖𝑛𝑣 (7-6) La energía en CA a la salida del inversor depende de la eficiencia de este, dato el cual viene dado por el fabricante: 𝑃𝑎𝑐=𝑃𝑖𝑛𝑣∙𝜂𝑖𝑛𝑣 (7-7) Sin embargo, el inversor también tiene un consumo propio mientras está funcionando. Al tener en cuenta el propio consumo del inversor, la ecuación anterior queda reducida a: 𝑃𝑎𝑐=𝑃𝑎𝑐,0 𝑃𝑎𝑐,0 𝜂𝑖𝑛𝑣−𝑃𝑠,0(𝑃𝑑𝑐𝑛𝑒𝑡−𝑃𝑠,0∙𝑁𝑖𝑛𝑣) (7-8) No solo eso, sino que, durante la noche, aunque los paneles no estén produciendo energía, los inversores sí están consumiendo un mínimo de energía que es necesario tener en cuenta: 𝑃𝑛𝑡=𝑃𝑛𝑡,0𝑁𝑖𝑛𝑣 (7-10) Una vez tenidas todas las consideraciones anteriores en cuenta, la energía bruta en CA se calcula como: 𝑃𝑎𝑐𝑔𝑟𝑜𝑠𝑠=𝑃𝑎𝑐𝑁𝑖𝑛𝑣 (7-11) De forma análoga a la energía en CC, es necesario sustraer las pérdidas que se daría en un sistema real, entre ellas [85] [86] [87] [88]: • Pérdidas en el transformador: las pérdidas debidas al transformador son dos. Por un lado, las pérdidas en carga, debidas a fenómenos como la histéresis, las corrientes de Foucault en el núcleo o las pérdidas dieléctricas. Por otro lado, las pérdidas en vacío, que se producen por la resistencia eléctrica del bobinado o por las corrientes de Foucault dentro de los conductores del mismo. • Pérdidas en el cableado: en las partes de CA de la instalación, habrá pérdidas debidas a la resistividad del conductor. Las pérdidas anteriores vienen dadas en porcentajes y la forma en la que influyen en la energía bruta en CA es la siguiente: 𝑃𝑎𝑐𝑛𝑒𝑡=𝑃𝑎𝑐𝑔𝑟𝑜𝑠𝑠∙(1−𝐿𝐴𝐶𝑤𝑖𝑟𝑖𝑛𝑔)−𝐿𝑡𝑟,𝑛𝑜 𝑙𝑜𝑎𝑑∙𝑃𝑎𝑐,0∙𝑁𝑖𝑛𝑣−𝐿𝑡𝑟,𝑙𝑜𝑎𝑑𝑃𝑎𝑐𝑔𝑟𝑜𝑠𝑠∙(1−𝐿𝐴𝐶𝑤𝑖𝑟𝑖𝑛𝑔) 𝑃𝑎𝑐,0∙𝑁𝑖𝑛𝑣 (7-12) Esta energía es la que finalmente podría distribuirse a la red de transporte o distribución en caso de provenir de una central solar fotovoltaica. Matriz de módulos fotovoltaicos 40 40 41 8 IMPLEMENTACIÓN EN JULIA Diagrama de flujo Con el fin de explicar la implementación en Julia de este modelo matemático, se proporciona un diagrama de flujo en la Figura 27 que describe su flujo de trabajo interno. Figura 27: Diagrama de flujo principal Implementación en Julia 42 42 En primer lugar, el usuario debe introducir un cierto número de parámetros según el modo en el que se desea trabajar. Después de esto, el programa inicializa automáticamente varias estructuras y variables necesarias durante la operación, por lo que esto es transparente para el usuario. Si no se proporciona ningún parámetro, el programa mostrará la documentación completa con todos los modos de uso. El primer parámetro debe ser siempre el modo de uso (Ver documentación en el apartado 0) para que el programa verifique que se ha introducido un modo válido. El programa mostrará la documentación completa si el primer parámetro no se corresponde con uno válido. Además, si el primer parámetro es correcto pero el número de parámetros del modo no coincide con los requeridos, el programa mostrará la documentación de ese modo específico. Tras comprobar que los parámetros introducidos por el usuario son correctos, el programa los utilizará para obtener los datos necesarios para el cálculo de la energía. En caso de que algún parámetro esté en un formato incorrecto (por ejemplo, un porcentaje fuera de los límites), se notificará al usuario con el error correspondiente para que pueda solucionarlo e introducir ese parámetro correctamente. Una vez realizadas todas las comprobaciones, el programa realiza los cálculos necesarios para obtener la energía del parque solar. Este cálculo se puede dividir en dos bloques: el primero calcula la irradiancia finalmente absorbida por un solo módulo y el segundo obtiene la energía utilizando esa irradiancia, entre otros datos. Por último, pero no menos importante, el programa mostrará la energía anual en CA y almacenará la energía obtenida en un archivo .CSV situado en la misma ruta que el código. El bloque de cálculo de la irradiancia utiliza únicamente las ecuaciones detalladas en los apartados anteriores y lo mismo ocurre con el bloque de cálculo de la energía neta de CA. Los cálculos internos de estos bloques se representan en la siguiente figura: Figura 28: Diagrama de flujo del cálculo de la irradiancia Figura 29: Diagrama de flujo del cálculo de la energía 43 43 Software de simulación para la generación fotovoltaica Documentación Como se ha mencionado anteriormente, hay ciertos casos en los que los parámetros no son correctos, ante lo que el programa responde mostrando la documentación. Los diferentes modos de uso se encuentran listados a continuación, explicando los parámetros que cada uno de los modos requiere. 8.1.1 Cuantitativo En este modo, el modelo matemático se calcula totalmente de forma cuantitativa. Esto significa que los parámetros del módulo, del inversor y del transformador se obtienen de la hoja de datos del fabricante y se siguen todas las fórmulas que se han explicado para cada uno de los bloques. La energía del módulo se calcula según las ecuaciones del apartado 0. La lista completa de parámetros es: Nombre Descripción Mode Nombre del modo ("Quantitative" en este caso) Source Ruta y nombre del archivo con los datos meteorológicos a usar Albedo Albedo entre 0 y 1 (si se desconoce, marcar a 0.2) ModuleVmpp Voltaje máximo del módulo en voltios ModuleImpp Corriente máxima del módulo en amperios ModuleKisc Coeficiente de corriente de cortocircuito en porcentaje ModuleKvoc Coeficiente de voltaje de circuito abierto en porcentaje ModuleKp Coeficiente de energía en porcentaje Am Área del módulo en metros cuadrados Tnoct Temperatura NOCT del módulo en grados Track_opt Seguimiento del parque: 0 para fijo, 1 para azimut y 2 para ambos ejes beta_0 Ángulo de inclinación de los módulos del parque. Entre 0 y 90 grados azimuth_0 Ángulo de orientación de los módulos del parque. Entre 0 y 180 grados Nsub Número de submatrices del parque: Máximo de 4 Nstr Número de cadenas por cada submatriz Nmod Número de módulos por cadena Inv_Pdc Energía de entrada máxima en CC del inversor en vatios Inv_Pac Energía de salida máxima en CA del inversor en vatios. Debe ser menor que Inv_Pac Inv_Vmax Voltaje máximo en CC en voltios. Debe ser mayor que el voltaje de circuito abierto de una cadena Inv_MinMPPT Voltaje MPPT mínimo en voltios. Debe ser menor que el voltaje máximo de una cadena Implementación en Julia 44 44 Inv_MaxMPPT Voltaje MPPT máximo en voltios. Debe ser mayor que el voltaje máximo de una cadena Inv_Ps0 Consumo del inversor durante la operación en vatios Inv_Pnt0 Consumo nocturno del inversor en vatios O Pérdidas por suciedad en porcentaje Lnameplate Pérdidas nominales (nameplate) en porcentaje Lmismatch Pérdidas por desajuste (mismatch) en porcentaje L_dcwiring Pérdidas en CC en el cableado en porcentaje L_acwiring Pérdidas en CA en el cableado en porcentaje Ltrans_load Pérdidas del transformador en carga en porcentaje Ltrans_noload Pérdidas del transformador en vacío en porcentaje Material Material del módulo. “Mono-c-Si” para monocristalino, cualquier otra cosa será considerada policristalino IrrComp Disponibilidad de las componentes de la irradiancia. 1 para disponible, 0 para desconocidas (se aplicará Engerer2) Tabla 3: Parámetros del modo cuantitativo 45 45 Software de simulación para la generación fotovoltaica 8.1.2 Módulo cualitativo (usando el nombre) En este modo, el único parámetro cualitativo es el nombre del módulo solar. Una vez proporcionado el nombre del módulo, los parámetros del este se obtendrán de una base de datos y las ecuaciones para calcular la energía de salida pertenecen a la sección 0. El resto de los parámetros sigue siendo el mismo y la lista completa es: Nombre Descripción Mode Nombre del modo ("Module_Qualitative1" en este caso) Source Ruta y nombre del archivo con los datos meteorológicos a usar Albedo Albedo entre 0 y 1 (si se desconoce, marcar a 0.2) ModuleName Nombre completo del módulo solar Track_opt Seguimiento del parque: 0 para fijo, 1 para azimut y 2 para ambos ejes beta_0 Ángulo de inclinación de los módulos del parque. Entre 0 y 90 grados azimuth_0 Ángulo de orientación de los módulos del parque. Entre 0 y 180 grados Nsub Número de submatrices del parque: Máximo de 4 Nstr Número de cadenas por cada submatriz Nmod Número de módulos por cadena Inv_Pdc Energía de entrada máxima en CC del inversor en vatios Inv_Pac Energía de salida máxima en CA del inversor en vatios. Debe ser menor que Inv_Pac Inv_Vmax Voltaje máximo en CC en voltios. Debe ser mayor que el voltaje de circuito abierto de una cadena Inv_MinMPPT Voltaje MPPT mínimo en voltios. Debe ser menor que el voltaje máximo de una cadena Inv_MaxMPPT Voltaje MPPT máximo en voltios. Debe ser mayor que el voltaje máximo de una cadena Inv_Ps0 Consumo del inversor durante la operación en vatios Inv_Pnt0 Consumo nocturno del inversor en vatios O Pérdidas por suciedad en porcentaje Lnameplate Pérdidas nominales (nameplate) en porcentaje Lmismatch Pérdidas por desajuste (mismatch) en porcentaje L_dcwiring Pérdidas en CC en el cableado en porcentaje L_acwiring Pérdidas en CA en el cableado en porcentaje Ltrans_load Pérdidas del transformador en carga en porcentaje Implementación en Julia 46 46 Ltrans_noload Pérdidas del transformador en vacío en porcentaje IrrComp Disponibilidad de las componentes de la irradiancia. 1 para disponible, 0 para desconocidas (se aplicará Engerer2) Tabla 4: Parámetros del módulo cualitativo (usando el nombre) 47 47 Software de simulación para la generación fotovoltaica 8.1.3 Módulo cualitativo (usando eficiencia) En este modo, los únicos parámetros cualitativos son la energía nominal y la superficie del módulo solar, que se utilizarán para calcular la eficiencia del panel. Con esta eficiencia, el programa buscará en la base de datos el módulo con la eficiencia más cercana y utilizará los parámetros correspondientes para calcular la energía. La lista de parámetros de este modo es: Nombre Descripción Mode Nombre del modo ("Module_Qualitative2" en este caso) Source Ruta y nombre del archivo con los datos meteorológicos a usar Albedo Albedo entre 0 y 1 (si se desconoce, marcar a 0.2) ModulePower Energía del módulo solar en vatios ModuleArea Área del módulo en metros cuadrados Track_opt Seguimiento del parque: 0 para fijo, 1 para azimut y 2 para ambos ejes beta_0 Ángulo de inclinación de los módulos del parque. Entre 0 y 90 grados azimuth_0 Ángulo de orientación de los módulos del parque. Entre 0 y 180 grados Nsub Número de submatrices del parque: Máximo de 4 Nstr Número de cadenas por cada submatriz Nmod Número de módulos por cadena Inv_Pdc Energía de entrada máxima en CC del inversor en vatios Inv_Pac Energía de salida máxima en CA del inversor en vatios. Debe ser menor que Inv_Pac Inv_Vmax Voltaje máximo en CC en voltios. Debe ser mayor que el voltaje de circuito abierto de una cadena Inv_MinMPPT Voltaje MPPT mínimo en voltios. Debe ser menor que el voltaje máximo de una cadena Inv_MaxMPPT Voltaje MPPT máximo en voltios. Debe ser mayor que el voltaje máximo de una cadena Inv_Ps0 Consumo del inversor durante la operación en vatios Inv_Pnt0 Consumo nocturno del inversor en vatios O Pérdidas por suciedad en porcentaje Lnameplate Pérdidas nominales (nameplate) en porcentaje Lmismatch Pérdidas por desajuste (mismatch) en porcentaje L_dcwiring Pérdidas en CC en el cableado en porcentaje L_acwiring Pérdidas en CA en el cableado en porcentaje Implementación en Julia 54 54 IrrComp Disponibilidad de las componentes de la irradiancia. 1 para disponible, 0 para desconocidas (se aplicará Engerer2) Tabla 10: Parámetros del inversor cualitativo 55 55 Software de simulación para la generación fotovoltaica 8.1.9 Inversor y pérdidas cualitativas En este modo, el usuario debe introducir tanto el grado de eficiencia del inversor como el grado de pérdidas en el modelo. La lista de parámetros es: Nombre Descripción Mode Nombre del modo ("Inverter_Losses_Qualitative" en este caso) Source Ruta y nombre del archivo con los datos meteorológicos a usar Albedo Albedo entre 0 y 1 (si se desconoce, marcar a 0.2) ModuleVmpp Voltaje máximo del módulo en voltios ModuleImpp Corriente máxima del módulo en amperios ModuleKisc Coeficiente de corriente de cortocircuito en porcentaje ModuleKvoc Coeficiente de voltaje de circuito abierto en porcentaje ModuleKp Coeficiente de energía en porcentaje Am Área del módulo en metros cuadrados Tnoct Temperatura NOCT del módulo en grados Track_opt Seguimiento del parque: 0 para fijo, 1 para azimut y 2 para ambos ejes beta_0 Ángulo de inclinación de los módulos del parque. Entre 0 y 90 grados azimuth_0 Ángulo de orientación de los módulos del parque. Entre 0 y 180 grados Nsub Número de submatrices del parque: Máximo de 4 Nstr Número de cadenas por cada submatriz Nmod Número de módulos por cadena Inv_eff Grado de eficiencia del inversor. De 1 a 5 en orden de menor a mayor O Pérdidas por suciedad en porcentaje Leff Grado de pérdidas. De 1 a 5 en orden de menor a mayor Material Material del módulo. “Mono-c-Si” para monocristalino, cualquier otra cosa será considerada policristalino IrrComp Disponibilidad de las componentes de la irradiancia. 1 para disponible, 0 para desconocidas (se aplicará Engerer2) Tabla 11: Parámetros del inversor y pérdidas cualitativas Implementación en Julia 56 56 8.1.10 Pérdidas cualitativas El único parámetro cualitativo de este modo son las pérdidas, que es, una vez más, un grado de pérdidas proporcionado por el usuario. La lista de parámetros es: Nombre Descripción Mode Nombre del modo ("Losses_Qualitative" en este caso) Source Ruta y nombre del archivo con los datos meteorológicos a usar Albedo Albedo entre 0 y 1 (si se desconoce, marcar a 0.2) ModuleVmpp Voltaje máximo del módulo en voltios ModuleImpp Corriente máxima del módulo en amperios ModuleKisc Coeficiente de corriente de cortocircuito en porcentaje ModuleKvoc Coeficiente de voltaje de circuito abierto en porcentaje ModuleKp Coeficiente de energía en porcentaje Am Área del módulo en metros cuadrados Tnoct Temperatura NOCT del módulo en grados Track_opt Seguimiento del parque: 0 para fijo, 1 para azimut y 2 para ambos ejes beta_0 Ángulo de inclinación de los módulos del parque. Entre 0 y 90 grados azimuth_0 Ángulo de orientación de los módulos del parque. Entre 0 y 180 grados Nsub Número de submatrices del parque: Máximo de 4 Nstr Número de cadenas por cada submatriz Nmod Número de módulos por cadena Inv_Pdc Energía de entrada máxima en CC del inversor en vatios Inv_Pac Energía de salida máxima en CA del inversor en vatios. Debe ser menor que Inv_Pac Inv_Vmax Voltaje máximo en CC en voltios. Debe ser mayor que el voltaje de circuito abierto de una cadena Inv_MinMPPT Voltaje MPPT mínimo en voltios. Debe ser menor que el voltaje máximo de una cadena Inv_MaxMPPT Voltaje MPPT máximo en voltios. Debe ser mayor que el voltaje máximo de una cadena Inv_Ps0 Consumo del inversor durante la operación en vatios 57 57 Software de simulación para la generación fotovoltaica Inv_Pnt0 Consumo nocturno del inversor en vatios O Pérdidas por suciedad en porcentaje Leff Grado de pérdidas. De 1 a 5 en orden de menor a mayor Material Material del módulo. “Mono-c-Si” para monocristalino, cualquier otra cosa será considerada policristalino IrrComp Disponibilidad de las componentes de la irradiancia. 1 para disponible, 0 para desconocidas (se aplicará Engerer2) Tabla 12: Parámetros de las pérdidas cualitativas Implementación en Julia 58 58 Descripción del paquete El lenguaje de programación utilizado para desarrollar este modelo es Julia en su versión 0.6.1. El programa se encuentra en la carpeta Function_codes, que contiene los siguientes archivos: Nombre Descripción Documentation.txt Este archivo contiene toda la documentación en caso de que el programa necesite mostrarla por pantalla. Engerer2.jl Este archivo contiene la función Engerer2, la cual calcula las componentes de la irradiancia. Gross_AC_Power.jl Este archivo contiene la función Gross_AC_Power, la cual calcula la energía bruta en CA de la matriz. Gross_DC_Power.jl Este archivo contiene la función Gross_DC_Power, la cual calcula la energía bruta en CC de la matriz. Load_Data.jl Este archivo contiene la función Load_Data, la cual obtiene todos los datos necesarios una vez se han introducido los parámetros. Main.jl Este es el archivo principal del programa y es el que llama al resto de funciones. Model.jl Este archivo contiene la función Model, la cual calcula la energía del módulo cuando un modo cualitativo es seleccionado. Module_Irradaince.jl Este archivo contiene la función Module_Irradiance, la cual calcula la irradiancia finalmente absorbida por el panel solar. Module_Power.jl Este archivo contiene la función Module_Power, la cual calcula la energía del panel solar para todos los modos. Net_AC_Power.jl Este archivo contiene la función Net_AC_Power, la cual calcula la energía neta en CA de la matriz. Net_DC_Power.jl Este archivo contiene la función Net_DC_Power, la cual calcula la energía neta en CC de la matriz. REST2.jl Este archivo contiene la función REST2, la cual calcula la irradiancia en un día despejado con un modelo de dos bandas. Results.csv Este archivo contiene la energía de salida de la matriz. Save_Data.jl Este archivo contiene la función Save_Data, la cual guarda la energía de salida en el archivo Results.csv. Shading_Soiling.jl Este archivo contiene la función Shading_Soiling, la cual calcula la irradiancia tras sombras y suciedad. Show_Doc.jl Este archivo contiene la función Show_Doc, que muestra por pantalla la parte correspondiente del archivo Documentation.txt. 59 59 Software de simulación para la generación fotovoltaica Sun_Angles.jl Este archivo contiene la función Sun_Angles, la cual calcula las componentes del sistema de coordenadas horizontal. Tilted_Irradiance.jl Este archivo contiene la función Tilted_Irradiance, la cual calcula la irradiancia en el plano de la matriz. Weather_File_Data.jl Este archivo contiene la función Weather_File_Data, la cual lee y guarda información del archive meteorológico introducido como parámetro. Tabla 13: Lista de archivos en la carpeta Function_codes La carpeta restante es Modules y en ella se encuentra un archivo con una enorme base de datos de numerosos módulos solares. Esta carpeta es necesaria en caso de que se seleccione cualquier modo que utilice un módulo cualitativo. El archivo contiene no sólo todos los datos proporcionados en la ficha técnica de cada módulo, sino también algunos parámetros que se han calculado por separado y que son necesarios para el modelo del apartado 0. En cuanto al archivo meteorológico dado como parámetro de entrada (que debe estar en formato .CSV), debe contener al menos lo siguiente • Una celda con la etiqueta Latitude y en la celda inferior el valor de la latitud local. • Una celda con la etiqueta Longitude y en la celda inferior el valor de la longitud local. • Una celda con la etiqueta Time Zone y en la celda inferior del valor de la zona horaria local. • Una celda con la etiqueta Elevation y en la celda inferior el valor de la altitud local. • Una celda con la etiqueta Year y en la celda inferior el valor del año. • Una celda con la etiqueta Month y en las celdas inferior el valor del mes. • Una celda con la etiqueta Day y en las celdas inferior del valor del día. • Una celda con la etiqueta Hour y en las celdas inferior del valor de la hora. • Una celda con la etiqueta DNI o Beam y en las celdas inferiores el valor de la irradiancia directa. • Una celda con la etiqueta DHI o Diffuse y en las celdas inferiores el valor de la irradiancia difusa. • Una celda con la etiqueta Tdry o Temperature y en las celdas inferiores el valor de la temperatura local. • Una celda con la etiqueta Wspd o Wind Speed y en las celdas inferiores el valor de la velocidad del viento. En la Figura 30 se muestra un ejemplo de archivo con las especificaciones requeridas. Si las componentes DNI y DHI se desconocen y deben ser estimadas, en vez de las respectivas celdas, el archivo debe contener: • Una celda con la etiqueta GHI o Global y en las celdas inferiores el valor de la irradiancia global. Implementación en Julia 60 60 Figura 30: Ejemplo de archivo meteorológico 61 9 RESULTADOS Tras haber explicado en profundidad los principios matemáticos de la generación fotovoltaica, así como los archivos involucrados para la ejecución del programa, se procederá a presentar los resultados obtenidos a fin de conocer la fiabilidad del modelo desarrollado. Modelo cualitativo Para poder comprobar la validez del modelo implementado en Julia, los resultados obtenidos por este se van a comparar con los del programa desarrollado por NREL: System Advisor Model (SAM) [90]. SAM es un software tecnológico y económico gratuito para realizar simulaciones de diversas tecnologías renovables tales como fotovoltaica o eólica. Al descargarse uno el programa anterior, se provee al usuario de forma automática con archivos .CSV que contienen datos meteorológicos de un año típico para un gran número de ciudades. Entre dichos archivos se encuentra la ciudad de Sevilla, que es la que se utilizará para esta prueba de validación. Del mismo modo, también se provee al usuario de una base de datos de paneles fotovoltaicos con los parámetros necesarios para el utilizar el modelo de circuito equivalente con un solo diodo (sección 0). Las entradas y parámetros que se han configurado tanto en SAM como en el modelo desarrollado son las siguientes (las entradas, salidas y parámetros de cada módulo están detallados en el Apéndice I): ENTRADAS Nombre Descripción Valor [unidades] 𝑻𝒔𝒕𝒂𝒎𝒑 Marca temporal Horaria [ℎ] 𝑫𝑵𝑰 Irradiancia directa normal Año típico [𝑊/𝑚2] 𝑫𝑯𝑰 Irradiancia difusa horizontal Año típico [𝑊/𝑚2] 𝑻𝒂𝒎𝒃 Temperatura ambiente Año típico [°𝐶] 𝑾𝒔𝒑𝒆𝒆𝒅 Velocidad del viento Año típico [𝑚/𝑠2] Tabla 14: Entradas de la prueba cualitativa PARÁMETROS Módulo 1: Angles for horizontal coordinate system Nombre Descripción Valor [unidades] Lat Latitud local 37.42 [°] Long Longitud local -5.9 [°] A Altitud local 31 [𝑚] 𝑡𝑧 Zona horaria local 1 [ℎ] Resultados 62 62 Módulo 3: Plane Of Array Irradiance Name Description Value [units] 𝜌 Albedo 0.2 [Adimensional] 𝛾𝑠 Ángulo azimuthal del panel 180 [°] 𝛽𝑠 Ángulo de inclinación del panel 30 [°] 𝑇𝑟𝑎𝑐𝑘𝑜𝑝𝑡 Opción de seguimiento 0 [Adimensional] Módulo 4: Soiling derate Name Description Value [units] 𝑂 Pérdidas por suciedad 5 [%] Módulo 5: Effective Irradiance Name Description Value [units] 𝐴 Altitud local 31 [𝑚] 𝛽𝑠 Ángulo de inclinación del panel 30 [°] 𝑇𝑟𝑎𝑐𝑘𝑜𝑝𝑡 Opción de seguimiento 0 [Adimensional] Módulo 6: Module Power Name Description Value [units] 𝑛𝑎𝑚𝑒 Modelo de panel solar SunPower SPR-X21-335- BLK [Adimensional] Módulo 7: Net DC Array Power Name Description Value [units] 𝑁𝑠𝑢𝑏 Número de submatrices 4 [Adimensional] 𝑁𝑠𝑡𝑟 Número de cadenas por submatriz 250 [Adimensional] 𝑁𝑚 Número de módulos por cadena 6 [Adimensional] 𝐿𝑛𝑎𝑚𝑒𝑝𝑙𝑎𝑡𝑒 Pérdidas nominales (nameplate) 0 [%] 𝐿𝑚𝑖𝑠𝑚𝑎𝑡𝑐ℎ Pérdidas por desajuste (mismatch) 2 [%] 𝐿𝐷𝐶𝑤𝑖𝑟𝑖𝑛𝑔 Pérdidas DC por cableado 2 [%] Módulo 8: NET AC Array Power Name Description Value [units] 𝑁𝑖𝑛𝑣 Número de inversores 4 [Adimensional] 𝜑 Factor de potencia 1 [Adimensional] 𝜂𝑖𝑛𝑣 Eficiencia del inversor 96 [%] 𝑃𝑎𝑐,0 Máxima energía AC del inversor 500000 [𝑊] 𝑃𝑠,0 Consumo del inversor durante la 4000 [𝑊] 63 63 Software de simulación para la generación fotovoltaica operación 𝑃𝑛𝑡,0 Consumo nocturno del inversor 125 [𝑊] 𝐿𝑡𝑟,𝑙𝑜𝑎𝑑 Pérdidas en carga del transformador 2.5 [%] 𝐿𝑡𝑟,𝑛𝑜 𝑙𝑜𝑎𝑑 Pérdidas en vacío del transformador 0.86 [%] 𝐿𝐴𝐶𝑤𝑖𝑟𝑖𝑛𝑔 Pérdidas AC por cableado 1.67 [%] Tabla 15: Parámetros de la prueba cualitativa Las métricas utilizadas para analizar el rendimiento de los paneles son el error medio absoluto (MAE), el error medio absoluto porcentual (MAPE) y el error cuadrático medio (MSE): 𝑀𝐴𝐸=1𝑛∑|𝑃𝑚𝑜𝑑𝑒𝑙−𝑃𝑆𝐴𝑀| 𝑛 𝑖=1 𝑀𝐴𝑃𝐸=1𝑛∑|𝑃𝑚𝑜𝑑𝑒𝑙−𝑃𝑆𝐴𝑀 𝑃𝑆𝐴𝑀 | 𝑛 𝑖=1 𝑀𝑆𝐸=1𝑛∑(𝑃𝑚𝑜𝑑𝑒𝑙−𝑃𝑆𝐴𝑀)2 𝑛 𝑖=1 (9-1) Las gráficas obtenidas son las siguientes: Figura 31: Histograma del AE para la prueba cualitativa Resultados 70 70 𝑃𝑠,0 Consumo del inversor durante la operación 4000 [𝑊] 𝑃𝑛𝑡,0 Consumo nocturno del inversor 125 [𝑊] 𝐿𝑡𝑟,𝑙𝑜𝑎𝑑 Pérdidas en carga del transformador 2.5 [%] 𝐿𝑡𝑟,𝑛𝑜 𝑙𝑜𝑎𝑑 Pérdidas en vacío del transformador 0.86 [%] 𝐿𝐴𝐶𝑤𝑖𝑟𝑖𝑛𝑔 Pérdidas AC por cableado 1.67 [%] Tabla 19: Parámetros de la prueba cuantitativa De la misma forma que en la prueba anterior, se van a utilizar las mismas métricas. Además, las gráficas solo tienen en cuenta aquellas horas en las que la irradiancia sea nula (noche) o mayor de 100 W/m2 (día sin amanecer ni atardecer) para eliminar los errores introducidos al amanecer y atardecer que son despreciables a efectos prácticos. Como se puede apreciar en los distintos histogramas, las métricas ahora alcanzan valores más altos. Esto era de esperar debido que se están comparando dos modelos completamente diferentes en las que uno de ellos (SAM) tiene un número de parámetros significativamente más alto. Figura 37: Histograma del AE para la prueba cuantitativa (tras eliminar amanecer y anochecer) 71 71 Software de simulación para la generación fotovoltaica Figura 38: Histograma del APE para la prueba cuantitativa (tras eliminar amanecer y anochecer) Los errores medios derivados de las gráficas anteriores son 3.20 kW, 0.83 % y 26.97 kW2 para el MAE, MAPE y MSE, respectivamente. Esto supone un incremento en las métricas del 110,5 %, 336.8 % y 244.4 % con respecto al modelo de circuito equivalente. Sin embargo, hay que tener en cuenta que el número de parámetros que modelan el panel fotovoltaico se ha reducido de cinco a uno. A pesar de todo lo comentado anteriormente, el análisis estadístico con percentiles recogido en la Tabla 20 muestra que el 99 % de las muestras presenta un error inferior a 16.10 kW, 1.70 % y 103.73 kW2 para el MAE, MAPE y MSE, respectivamente. Si bien es cierto que las métricas no son especialmente altas para el sistema en cuestión teniendo en cuenta que solo hay un parámetro que modele el panel, se aconseja dejar este modo para aquellos consumidores domésticos o comerciales que no necesiten una gran precisión en los resultados. Para aquellos usuarios cuya operación financiera este intrínsecamente ligada con la producción fotovoltaica, tales como un parque solar, se aconseja utilizar el modelo del circuito equivalente. Figura 39: Histograma del SE para la prueba cuantitativa (tras eliminar amanecer y anochecer) Resultados 72 72 Concepto AE APE SE Muestras 7186 7186 7186 Media 3.20 0.83 26.97 Desviación estándar 4.09 1.68 42.45 Percentil del 25% 0.00 0.00 0.00 Percentil del 50% 0.00 0.00 0.00 Percentil del 75% 6.44 0.68 41.49 IQR 6.44 0.68 41.49 -1.5*IQR -9.66 -1.02 -62.24 +1.5*IQR 16.10 1.70 103.73 Tabla 20: Análisis percentil para la prueba cuantitativa (tras eliminar amanecer y anochecer) Modelo cualitativo con GHI Como se ha comentado en el presente documento, existen casos en los que solamente se cuenta con un piranómetro para medir la componente de la irradiancia GHI y el resto no están disponible. Es por este motivo que se han incorporado los modelos de REST2 y Engerer2 para poder obtener la DNI y DHI a partir de la GHI. Si bien es cierto que se han utilizado los valores típicos de los parámetros involucrados en las transmitancias, esto sirve para tener una idea general de la producción energética que se va a obtener. Se va a proceder a repetir la primera prueba, pero, esta vez, usando los modelos de REST2 y Engerer2 para ver como de diferente resulta la salida. Una vez más, los resultados obtenidos mediante el software SAM se tomarán como valores reales con los que comparar el modelo desarrollado y poder estudiar las métricas necesarias. Las entradas y parámetros que se han configurado en el modelo desarrollado (las de SAM no han cambiado) son las siguientes (las entradas, salidas y parámetros de cada módulo están detallados en el Apéndice I): ENTRADAS Nombre Descripción Valor [unidades] 𝑻𝒔𝒕𝒂𝒎𝒑 Marca temporal Horaria [ℎ] 𝑮𝑯𝑰 Irradiancia global horizontal Año típico [𝑊/𝑚2] 𝑻𝒂𝒎𝒃 Temperatura ambiente Año típico [°𝐶] 𝑾𝒔𝒑𝒆𝒆𝒅 Velocidad del viento Año típico [𝑚/𝑠2] Tabla 21: Entradas de la prueba cualitativa con GHI 73 73 Software de simulación para la generación fotovoltaica PARÁMETROS Módulo 1: Angles for horizontal coordinate system Nombre Descripción Valor [unidades] Lat Latitud local 37.42 [°] Long Longitud local -5.9 [°] A Altitud local 31 [𝑚] 𝑡𝑧 Zona horaria local 1 [ℎ] Módulo 3: Plane Of Array Irradiance Name Description Value [units] 𝜌 Albedo 0.2 [Adimensional] 𝛾𝑠 Ángulo azimuthal del panel 180 [°] 𝛽𝑠 Ángulo de inclinación del panel 30 [°] 𝑇𝑟𝑎𝑐𝑘𝑜𝑝𝑡 Opción de seguimiento 0 [Adimensional] Módulo 4: Soiling derate Name Description Value [units] 𝑂 Pérdidas por suciedad 5 [%] Módulo 5: Effective Irradiance Name Description Value [units] 𝐴 Altitud local 31 [𝑚] 𝛽𝑠 Ángulo de inclinación del panel 30 [°] 𝑇𝑟𝑎𝑐𝑘𝑜𝑝𝑡 Opción de seguimiento 0 [Adimensional] Módulo 6: Module Power Name Description Value [units] 𝑛𝑎𝑚𝑒 Modelo de panel solar SunPower SPR-X21-335- BLK [Adimensional] Módulo 7: Net DC Array Power Name Description Value [units] 𝑁𝑠𝑢𝑏 Número de submatrices 4 [Adimensional] 𝑁𝑠𝑡𝑟 Número de cadenas por submatriz 250 [Adimensional] 𝑁𝑚 Número de módulos por cadena 6 [Adimensional] 𝐿𝑛𝑎𝑚𝑒𝑝𝑙𝑎𝑡𝑒 Pérdidas nominales (nameplate) 0 [%] 𝐿𝑚𝑖𝑠𝑚𝑎𝑡𝑐ℎ Pérdidas por desajuste (mismatch) 2 [%] 𝐿𝐷𝐶𝑤𝑖𝑟𝑖𝑛𝑔 Pérdidas DC por cableado 2 [%] Resultados 74 74 Módulo 8: NET AC Array Power Name Description Value [units] 𝑁𝑖𝑛𝑣 Número de inversores 4 [Adimensional] 𝜑 Factor de potencia 1 [Adimensional] 𝜂𝑖𝑛𝑣 Eficiencia del inversor 96 [%] 𝑃𝑎𝑐,0 Máxima energía AC del inversor 500000 [𝑊] 𝑃𝑠,0 Consumo del inversor durante la operación 4000 [𝑊] 𝑃𝑛𝑡,0 Consumo nocturno del inversor 125 [𝑊] 𝐿𝑡𝑟,𝑙𝑜𝑎𝑑 Pérdidas en carga del transformador 2.5 [%] 𝐿𝑡𝑟,𝑛𝑜 𝑙𝑜𝑎𝑑 Pérdidas en vacío del transformador 0.86 [%] 𝐿𝐴𝐶𝑤𝑖𝑟𝑖𝑛𝑔 Pérdidas AC por cableado 1.67 [%] Tabla 22: Parámetros de la prueba cualitativa con GHI De la misma forma que en la primera prueba, se van a utilizar las mismas métricas. Además, las gráficas solo tienen en cuenta aquellas horas en las que la irradiancia sea nula (noche) o mayor de 100 W/m2 (día sin amanecer ni atardecer) para eliminar los errores introducidos al amanecer y atardecer que son despreciables a efectos prácticos. Como se puede apreciar en los distintos histogramas, las métricas ahora alcanzan valores más altos. Esto era de esperar debido que se están estimando las componentes de la irradiancia a partir de la irradiancia global con un modelo que depende de un gran número de parámetros. Figura 40: Histograma del AE para la prueba cualitativa con GHI (tras eliminar amanecer y anochecer) 75 75 Software de simulación para la generación fotovoltaica Figura 41: Histograma del APE para la prueba cualitativa con GHI (tras eliminar amanecer y anochecer) Los errores medios derivados de las gráficas anteriores son 5.99 kW, 1.11 % y 168.64 kW2 para el MAE, MAPE y MSE, respectivamente. Esto supone un incremento en las métricas del 296.7 %, 484.2 % y 2053.8 % con respecto a usar todas las componentes de la irradiancia. Sin embargo, hay que tener en cuenta que se han modelado dichas componentes a partir de valores típicos de todos los parámetros involucrados en las transmitancias en vez de usar los valores reales. A pesar de todo lo comentado anteriormente, el análisis estadístico con percentiles recogido en la Tabla 23 muestra que el 99 % de las muestras presenta un error inferior a 17.53 kW, 2.20 % y 122.98 kW2 para el MAE, MAPE y MSE, respectivamente. Las métricas son bastante más altas en esta prueba, especialmente para el MSE, por lo que se aconseja dejar este modo para aquellos consumidores domésticos o comerciales que no necesiten una gran precisión en los resultados. Para aquellos usuarios cuya operación financiera este intrínsecamente ligada con la producción fotovoltaica, tales como un parque solar, se aconseja utilizar todas las componentes de la irradiancia. En caso de que esto no sea posible, la opción más viable es colocar el piranómetro en la misma dirección que el ángulo de incidencia de los paneles solares y utilizar la irradiancia medida como efectiva para calcular la salida en potencia. Resultados 76 76 Figura 42: Histograma del SE para la prueba cualitativa con GHI (tras eliminar amanecer y anochecer) Concepto AE APE SE Muestras 7186 7186 7186 Media 5.99 1.11 168.64 Desviación estándar 11.52 2.76 614.86 Percentil del 25% 0.00 0.00 0.00 Percentil del 50% 0.00 0.00 0.00 Percentil del 75% 7.01 0.88 49.19 IQR 7.01 0.88 49.19 -1.5*IQR -10.52 -1.32 -73.79 +1.5*IQR 17.53 2,20 122.98 Tabla 23: Análisis percentil para la prueba cualitativa con GHI (tras eliminar amanecer y anochecer) Panel solar fotovoltaico El panel solar es seguramente la parte más crítica del modelo fotovoltaico planteado y, es por ello, que su funcionamiento necesita ser validado de forma correcta. NREL ha puesto a disposición pública datos recabados de instalaciones solares reales con el objetivo concreto de poder validar los distintos modelos que se fueran desarrollando a lo largo del tiempo [92] [93] [94]. Estas instalaciones pertenecen a tres localizaciones diferentes: Cocoa, Florida; Eugene, Oregón; y Golden, colorado. Estas localizaciones atienen a tres tipos de climas diferentes según la clasificación de Köppen [95]: subtropical húmedo (Cfa), mediterráneo (Csb) y continental templado (Dfb), respectivamente. Además, la base 77 77 Software de simulación para la generación fotovoltaica de datos contiene módulos de diversas tecnologías, es decir, silicio monocristalino, silicio policristalino o teluro de cadmio, entre otros. Se han utilizado los datos pertenecientes a los paneles de silicio monocristalino y policristalino en las tres localizaciones provistas. La razón de esto es que es la tecnología más ampliamente utilizada tanto por consumidores domésticos como en grandes centrales fotovoltaicas. Del mismo modo, como en los casos mencionados anteriormente se suele contar únicamente con los parámetros proporcionados en la hoja de datos del fabricante, se ha utilizado el modo cuantitativo en la validación. Las métricas utilizadas para analizar el rendimiento de los paneles son el error medio absoluto porcentual (MAPE), el sesgo porcentual (MBPE), la raíz del error cuadrático medio porcentual (RMSPE) y el coeficiente de correlación de Pearson (r2): 𝑛𝑀𝐴𝑃𝐸= 𝑀𝐴𝑃𝐸 𝑃𝑛𝑜𝑚𝑖𝑛𝑎𝑙=1 𝑃𝑛𝑜𝑚𝑖𝑛𝑎𝑙1𝑛∑|𝑃𝑚𝑜𝑑𝑒𝑙−𝑃𝑣𝑎𝑙𝑖𝑑𝑎𝑡𝑖𝑜𝑛 𝑃𝑣𝑎𝑙𝑖𝑑𝑎𝑡𝑖𝑜𝑛 | 𝑛 𝑖=1 𝑛𝑀𝐵𝑃𝐸= 𝑀𝐵𝑃𝐸 𝑃𝑛𝑜𝑚𝑖𝑛𝑎𝑙=1 𝑃𝑛𝑜𝑚𝑖𝑛𝑎𝑙1𝑛∑(𝑃𝑚𝑜𝑑𝑒𝑙−𝑃𝑣𝑎𝑙𝑖𝑑𝑎𝑡𝑖𝑜𝑛 𝑃𝑣𝑎𝑙𝑖𝑑𝑎𝑡𝑖𝑜𝑛 ) 𝑛 𝑖=1 𝑛𝑅𝑀𝑆𝑃𝐸=𝑅𝑀𝑆𝑃𝐸 𝑃𝑛𝑜𝑚𝑖𝑛𝑎𝑙=1 𝑃𝑛𝑜𝑚𝑖𝑛𝑎𝑙√1𝑛∑(𝑃𝑚𝑜𝑑𝑒𝑙−𝑃𝑣𝑎𝑙𝑖𝑑𝑎𝑡𝑖𝑜𝑛 𝑃𝑣𝑎𝑙𝑖𝑑𝑎𝑡𝑖𝑜𝑛 )2 𝑛 𝑖=1 𝑟2=1−∑(𝑃𝑚𝑜𝑑𝑒𝑙−𝑃𝑣𝑎𝑙𝑖𝑑𝑎𝑡𝑖𝑜𝑛)2 𝑛𝑖=1 ∑(𝑃𝑚𝑜𝑑𝑒𝑙−𝑃𝑚𝑜𝑑𝑒𝑙          )2 𝑛𝑖=1 (9-2) A continuación, se muestran los resultados obtenidos para la localización de Cocoa (Figura 43), Eugene (Figura 44) y Golden (Figura 45), así como una tabla (Tabla 24) que recoge los resultados. Figura 43: Validación del modelo del panel solar en Cocoa Figura 44: Validación del modelo del panel solar en Eugene Resultados 78 78 Figura 45: Validación del modelo del panel solar en Golden Como se puede observar, los resultados generales son muy buenos. De las tres localizaciones, aquella con menores errores es Eugene, y la que cuenta con los errores más altos, Golden. Del mismo modo, los resultados son ligeramente mejores para la tecnología de silicio policristalino, independientemente de la localización. El coeficiente de correlación de Pearson es del 0.98 y cercano al 0.99 en algunos casos y los errores globales son generalmente bajos: 4.24% (nMAPE), -1.68% (nMBPE) y 5.95% (nRMSPE). Localización Clima Tecnología nMAPE (%) nMBPE (%) nRMSPE (%) r2 Cocoa Cfa Monocristalino 4,10 -2,74 5,73 0,988 Cocoa Cfa Policristalino 3,74 -1,07 5,14 0,989 Eugene Csb Monocristalino 3,60 -1,35 5,59 0,986 Eugene Csb Policristalino 3,70 0,06 5,41 0,986 Golden Dfb Monocristalino 5,47 -3,59 7,34 0,979 Golden Dfb Policristalino 4,80 -1,38 6,46 0,980 Media Global 4,24 -1,68 5,95 0,98 Media - Monocristalino 4,39 -2,56 6,22 0,98 Media - Policristalino 4,08 -0,80 5,67 0,99 Tabla 24: Métricas de validación del modelo de los paneles solares 79 10 CONCLUSIONES Tal y como se ha podido comprobar a lo largo del presente proyecto, la energía fotovoltaica se encuentra en uno de sus puntos más álgidos a lo largo de la historia y será una de las mayores conductoras de la transición energética que tendrá lugar durante los próximos años. Además, cada vez son más los usuarios domésticos y comerciales que están optando por este tipo de tecnología para poder abastecerse completa o parcialmente y reducir la dependencia de la red eléctrica. Con el objetivo de que los usuarios puedan realizar estimaciones de la cantidad de energía que se puede generar con un sistema fotovoltaico, se ha implementado un modelo matemático completo en el lenguaje de programación Julia. Dicho modelo permite configurar las distintas partes del sistema, por ejemplo, panel, inversor o ángulos, de acuerdo con las necesidades concretas de los usuarios. Los resultados del modelo implementado se han comparado con el software SAM desarrollado por NREL, uno de los mayores exponentes en desarrollo fotovoltaico. Los errores medios obtenidos son del 1.52 kW, 0.19 % y 7.83 kW2 para el MAE, MAPE y MSE, respectivamente. Por lo tanto, a nivel de resultados, los conseguidos con el modelo desarrollado son comparables, haciendo viable el modelo planteado en el presente proyecto. No solo eso, sino que se han incorporado diversas funcionalidades adicionales para dotar de aún más versatilidad al modelo desarrollado. Las más relevantes son la inclusión de un modelo simplificado para el panel fotovoltaico, ya que la mayoría de los usuarios no van a poder contar con todos los parámetros necesarios para el modelo de circuito equivalente. Además, también es posible estimar las componentes de la irradiancia en caso de que el usuario solo disponga de la irradiancia global. Hay que destacar que los resultados obtenidos en las pruebas para estos casos son de menor calidad que aquella en la que se utiliza el modelo conociendo todas las variables y parámetros. Sin embargo, esto era de esperar dado que, o bien se utilizan menos parámetros, o bien se utilizan valores típicos de estos. Además, estos modelos no están destinados a aquellos usuarios cuya retribución económica dependa de la generación fotovoltaica, sino para aquellos usuarios domésticos y comerciales que busquen conocer una estimación de su producción para gestionar sus recursos y que generalmente no dispondrán de toda la información para utilizar el modelo completo. Con este proyecto se pretendía brindar una introducción a la generación fotovoltaica, su funcionamiento y las matemáticas detrás de ellas. Gracias a este proyecto, el lector podrá tener un conocimiento bastante profundo de esta tecnología y podrá ayudarlo a decidir si subirse al barco del autoconsumo o convertirse en productor de energía eléctrica para ofrecer dicha energía a la red. Referencias 86 86 Available: https://www.electrical4u.com/working-principle-of-photovoltaic-cell-or-solar-cell/. [Último acceso: June 2021]. [40] «Solar Cell Construction & Working Principle,» Electrical Engineering 123, [En línea]. Available: https://electricalengineering123.com/solar-cell-construction-working-principle/. [Último acceso: June 2021]. [41] «ISO 8601,» Wikipedia, [En línea]. Available: https://en.wikipedia.org/wiki/ISO_8601. [Último acceso: September 2021]. [42] «ISO 8601:2004,» ISO, [En línea]. Available: https://www.iso.org/standard/40874.html. [Último acceso: September 2021]. [43] «Time Conventions,» National Aeronautics and Space Administration (NASA), [En línea]. Available: https://solarsystem.nasa.gov/basics/chapter2-3/. [Último acceso: September 2021]. [44] «Julian Date,» European Space Agency (ESA), 23 February 2012. [En línea]. Available: https://gssc.esa.int/navipedia/index.php/Julian_Date. [Último acceso: September 2021]. [45] «A Brief Note on Time Systems,» The Ohio State University, 10 January 2011. [En línea]. Available: http://www.astronomy.ohio-state.edu/~pogge/Ast350/timesys.html. [Último acceso: September 2021]. [46] «Julian Date,» Wolfram Research, [En línea]. Available: https://scienceworld.wolfram.com/astronomy/JulianDate.html. [Último acceso: September 2021]. [47] J. J. Michalsky, «The Astronomical Almanac's algorithm for approximate solar position (1950-2050),» Solar Energy, vol. 40, nº 3, pp. 227-235, 1988. [48] M. Iqbal, An Introduction to Solar Radiation, New York: NY: Academic Press, 1983. [49] «Ecliptic coordinates,» Museo Galileo, [En línea]. Available: https://catalogue.museogalileo.it/indepth/EclipticCoordinates.html. [Último acceso: September 2021]. [50] «Cosmic Coordinates,» Las Cumbres Observatory, [En línea]. Available: https://lco.global/spacebook/sky/equatorial-coordinate-system/. [Último acceso: September 2021]. [51] B. King, «Right Ascension & Declination: Celestial Coordinates for Beginners,» Sky & Telescope, 26 February 2019. [En línea]. Available: https://skyandtelescope.org/astronomy-resources/right-ascension- declination-celestial-coordinates/. [Último acceso: September 2021]. [52] «Solar Time vs Sidereal Time,» Las Cumbres Observatory, [En línea]. Available: https://lco.global/spacebook/sky/sidereal-time/. [Último acceso: September 2021]. [53] k. Bikos, «Altitude & Azimuth: The Horizontal Coordinate System,» Time and Date, [En línea]. Available: https://www.timeanddate.com/astronomy/horizontal-coordinate-system.html. [Último acceso: September 2021]. [54] D. Ford, «Equation Of Time,» In-The-Sky, 19 June 2020. [En línea]. Available: https://in-the- sky.org/article.php?term=equation_of_time. [Último acceso: September 2021]. 87 87 Software de simulación para la generación fotovoltaica [55] «Solar Irradaince,» National Aeronautics and Space Administration (NASA), January 2008. [En línea]. Available: https://www.nasa.gov/mission_pages/sdo/science/Solar%20Irradiance.html. [Último acceso: June 2021]. [56] «Introduction to Solar Radiation,» NewPort, [En línea]. Available: https://www.newport.com/t/introduction-to-solar-radiation. [57] «Multijunction III-V Photovoltaics research,» Office of Energy Efficiency & Renewable Energy (U.S. Department of Energy), [En línea]. Available: https://www.energy.gov/eere/solar/multijunction-iii-v- photovoltaics-research. [Último acceso: June 2021]. [58] «Air Mass,» PV Performance Modeling Collaborative, [En línea]. Available: https://pvpmc.sandia.gov/modeling-steps/1-weather-design-inputs/irradiance-and-insolation-2/air-mass/. [Último acceso: June 2021]. [59] C. Riordan y R. Hulston, «What is an air mass 1.5 spectrum? (solar cell performance calculations),» IEEE Conference on Photovoltaics Specialists, vol. 2, pp. 1085-1088, 1990. [60] «Reference Air Mass 1.5 Spectra,» National Renewable Energy Laboratory (NREL), [En línea]. Available: https://www.nrel.gov/grid/solar-resource/spectra-am1.5.html. [Último acceso: June 2021]. [61] F. Kasten y A. T. Young, «Revised optical air mass tables and approximation formula,» Applied optics, vol. 28, pp. 4735-4738, 1989. [62] «Solar Radiation,» National Aeronautics and Space Administration (NASA), [En línea]. Available: https://earth.gsfc.nasa.gov/climate/research/solar-radiation. [Último acceso: June 2021]. [63] G. Kopp y J. L. Lean, «A new, lower value of total solar irradiance: Evidence and climate significance,» Geophysical Research Letters, vol. 38, 2011. [64] «Welcome to the Home Page of the SOlar Raditaion and Climate Experiment,» National Aeronautics and Space Administration (NASA), [En línea]. Available: https://lasp.colorado.edu/home/sorce/. [Último acceso: June 2021]. [65] J. A. Duffie y W. A. Beckman, Solar Engineering of Thermal Processes, Wiley, 2013. [66] «How to measure solar radiation,» Hukseflux - Thermal Sensors, [En línea]. Available: https://www.hukseflux.com/applications/solar-energy-pv-system-performance-monitoring/how-to- measure-solar-radiation. [Último acceso: September 2021]. [67] N. A. Engerer, «Minute resolution estimates of the diffuse fraction of global irradiance for southeastern Australia,» Solar Energy, vol. 116, pp. 215-237, 2015. [68] C. A. Gueymard, «REST2: High-performance solar radiation model for cloudless-sky irradiance, illuminance, and photosynthetically active radiation - Validation with a benchmark dataset,» Solar Energy, vol. 82, nº 3, pp. 272-285, 2008. [69] C. A. Gueymard, «Direct solar transmittance and irradiance predictions with broadband models. Part I: detailed theoretical performance assessment,» Solar Energy, vol. 74, nº 5, pp. 355-378, 2003. [70] «Home,» Aerosol Robotic NEtwork (AERONET), [En línea]. Available: https://aeronet.gsfc.nasa.gov. [Último acceso: September 2021]. Referencias 88 88 [71] B. Liu y R. Jordan, «A Rational Procedure for Predicting The Long-term Average Performance of Flatplate Solar-energy Collectors,» Solar Energy, vol. 7, nº 2, pp. 53-74, 1963. [72] R. Perez, P. Ineichen, R. Seals, J. Michaslky y R. Stewart, «Modeling Daylight Availability and Irradiance Components from Direct and Global Irradiance,» Solar Energy, vol. 44, nº 5, pp. 271-289, 1990. [73] «Perez Sky Diffuse Model,» PV Performance Modeling Collaborative, [En línea]. Available: https://pvpmc.sandia.gov/modeling-steps/1-weather-design-inputs/plane-of-array-poa- irradiance/calculating-poa-irradiance/poa-sky-diffuse/perez-sky-diffuse-model/. [Último acceso: June 2021]. [74] «How to calculate the minimum distance between PV panels,» EasySolar, 20 April 2021. [En línea]. Available: https://easysolar-app.com/en/ufaqs/how-to-calculate-the-minimum-distance-between-pv- panels/. [Último acceso: September 2021]. [75] W. de Soto, «Improvement and Validation of a Model for Photovoltaic Array Performance,» University of Wisconsin-Madison, Madison, 2004. [76] W. de Soto, S. A. Klein y W. A. Beckman, «Improvement and validation of a model for photovoltaic array performance,» Solar Energy, vol. 80, pp. 78-88, 2006. [77] A. P. Dobos y S. M. MacAlpine, «Procedure for applying IEC-61853 test data to a single diode model,» 2014 IEEE 40th Photovoltaic Specialist Conference (PVSC), pp. 2846-2849, 2014. [78] «Module Circuit Design,» PVEducation, [En línea]. Available: https://www.pveducation.org/pvcdrom/modules-and-arrays/module-circuit-design. [Último acceso: September 2021]. [79] S. Murmson, «The Output Watts of Solar Panels,» Sciencing, 10 March 2018. [En línea]. Available: https://sciencing.com/average-photovoltaic-system-efficiency-7092.html. [Último acceso: September 2021]. [80] D. T. Cotfas, P. A. Cotfas y O. M. Machidon, «Study of Temperature Coefficients for Parameters of Photovoltaic Cells,» International Journal of Photoenergy, vol. 2018, nº 5945602, 2018. [81] E. Skoplaki y J. A. Palyvos, «On the temperature dependence of photovoltaic module electrical performance: A review of efficiency/power correlations,» Solar Energy, vol. 83, nº 5, pp. 614-624, 2009. [82] M. Seyedmahmoudian, B. Horan, T. Kok Soon , R. Rahmani, A. M. Muang Than Oo y A. Stojcevski, «State of the art artificial intelligence-based MPPT techniques for mitigating partial shading effects on PV systems - A review,» Renewable and Sustainable Energy Reviews, pp. 435-455, 2016. [83] «GRID-CONNECTED PC: OPERATIONS AND MAINTENANCE GUIDELINE,» Pacific Power Association, 2019. [En línea]. Available: https://www.ppa.org.fj/wp-content/uploads/2019/08/OM-PV- Systems-Guideline-V2-D2.1.pdf. [Último acceso: September 2021]. [84] T. Ma, H. Yang y L. Lu, «Solar photovoltaic system modeling and performance prediction,» Renewable and Sustainable Energy Reviews, vol. 36, pp. 304-315, 2014. 89 89 Software de simulación para la generación fotovoltaica [85] V. Ionin, «System Loss Diagram,» Aurora, September 2021. [En línea]. Available: https://help.aurorasolar.com/hc/en-us/articles/235994088-System-Loss-Diagram. [Último acceso: September 2021]. [86] C. Hopper, «System Losses,» Aurora, September 2021. [En línea]. Available: https://help.aurorasolar.com/hc/en-us/articles/220450107-System-Losses. [Último acceso: September 2021]. [87] P. Chauhan, «Solar PV Power-Simulation and Designing INTERSHIP REPORT (in partial fulfilment on VCE internship program),» 2020. [88] K. Omkar, M. V. Srikanth, K. P. Swaroop y R. Rama, «Performance evaluation of 50KWp rooftop solar PV plant,» de 2015 International Conference on Industrial Instrumentation and Control (ICIC), Pune, India, 2015. [89] D. L. King, S. González, G. M. Galbraith y W. E. Boyson, «Performance Model for Grid-Connected Photovoltaic Inverters,» Sandia National Laboratories, Albuquerque, 2007. [90] «System Advisor Model (SAM),» National Renewable Energy Laboratory (NREL), [En línea]. Available: https://sam.nrel.gov. [Último acceso: September 2021]. [91] SunPower, «SunPower X-Series: X21-350-BLK,» [En línea]. Available: https://cdn.enfsolar.com/Product/pdf/Crystalline/5ab85dd1ba85b.pdf. [Último acceso: September 2021]. [92] B. Marion, A. Anderberg, C. Deline, M. Muller, G. Perrin, J. Rodríguez, S. Rummel, T. Silverman, F. Vignola y S. Barkaszi, «Data for Validating Models for PV Module Performance,» U.S. Department of Energy - Office of Scientific and Technical Information, United States, 2021. [93] B. Marion, M. G. Deceglie y T. J. Silverman, «Analysis of measured photovoltaic module performance for Florida, Oregon, and Colorado locations,» Solar Energy, vol. 110, pp. 736-744, 2014. [94] W. Marion, A. Anderberg, C. Deline, S. Glick, M. Muller, G. Perrin, J. Rodríguez, S. Rummel, K. Terwilliger y T. J. Silvermand, «User’s Manual for Data for Validating Models for PV Module Performance,» National Renewable Energy Laboratory (NREL), Golden, Colorado, United States, 2014. [95] J. A. Arnfield, «Köppen climate classification,» Encyclopedia Britannica, 11 November 2020. [En línea]. Available: https://www.britannica.com/science/Koppen-climate-classification. [Último acceso: September 2021]. Referencias 90 90 91 GLOSARIO AC Altern Current AE Absolute Error AERONET Aerosol Robotic Network AM0 Air Mass 0 AM1 Air Mass 1 AOD Aerosol Optical Depth APE Absolute Percentage Error BOE Boletín Oficial del Estado CA Corriente Alterna CC Corriente Continua CIGS Cobre, Indio, Galio y Selenio CSV Comma Separated Value DC Direct Current DHI Diffuse Horizontal Irradiance DNI Direct Normal Irradiance ESA European Space Agency EU European Union FER Fuentes de Energía Renovables GHI Global Horizontal Irradiance H2020 Horizon 2020 HDKR Hay and Davies, Klucher and Reindl IEC International Electrotechnical Commission IMM Inverted Metamorphic IQR Interquartile Range JIT Just In Time MAE Mean Absolute Error MAPE Mean Absolute Percentage Error MBPE Mean Bias Percentage Error ML Machine Learning MPPT Maximum Power Point Tracking MSE Mean Squared Error NASA National Aeronautics and Space Administration NOCT Nominal Operative Cell Temperature NREL National Renewable Energy Laboratory PV Photovoltaics RD Real Decreto REE Red Eléctrica de España REPL Read Eval Print Loop RES Renewable Energy Sources REST2 Reference Evaluation of Solar Transmittance, 2 bands RMSPE Root Mean Squared Percentage Error SAM System Advisor Model SE Squared Error SI International System SORCE Solar Radiation and Climate Experiment Glosario 92 92 STC Standard Test Conditions TIOBE The Importance Of Being Earnest TSI Total Solar Irradiance UE Unión Europea UNEF Unión Española Fotovoltaica UTC Universal Coordinated Time 93 APÉNDICE I: RESUMEN DEL MODELO Esquema general Figura 46: Esquema general del modelo Entradas Nombre Descripción Unidades 𝑻𝒔𝒕𝒂𝒎𝒑 Marca temporal Decisión del usuario 𝑮 Irradiancia Vatios por metro cuadrado (W/m2) 𝑻𝒂𝒎𝒃 Temperatura ambiente Grados (°) 𝑾𝒔𝒑𝒆𝒆𝒅 Velocidad del viento Metros por segundo(m/s) Tabla 25: Entradas generales del modelo Salidas Nombre Descripción Unidades 𝑷𝒂𝒄𝒏𝒆𝒕 Energía neta en CA Vatios (W) Tabla 26: Salidas generales del modelo Apéndice I: Resumen del modelo 94 94 Módulo: Angles for horizontal coordinate system Figura 47: Esquema del módulo Angles for horitontal coordinate system Entradas Nombre Descripción Unidades 𝑻𝒔𝒕𝒂𝒎𝒑 Marca temporal Decisión del usuario Tabla 27: Entradas del módulo Angles for horitontal coordinate system Salidas Nombre Descripción Unidades 𝒁 Ángulo cenital Grados (°) 𝜶 Ángulo de elevación Grados (°) 𝜸 Ángulo acimutal Grados (°) Tabla 28: Salidas del módulo Angles for horitontal coordinate system Parámetros Nombre Descripción Unidades 𝑳𝒂𝒕 Latitud local Grados (°) 𝑳𝒐𝒏𝒈 Longitud local Grados (°) 𝑨 Altitud local Metros (m) 𝒕𝒛 Huso horario Horas (h) Tabla 29: Parámetros del módulo Angles for horitontal coordinate system 95 95 Software de simulación para la generación fotovoltaica Módulo: Horizontal Irradiance Figura 48: Esquema del módulo Horizontal Irradiance Entradas Nombre Descripción Unidades 𝒁 Ángulo cenital Grados (°) 𝑨 Altitud local Metros (m) 𝑮𝟎𝒉 Irradiancia extraterrestre horizontal Vatios por metro cuadrado (𝑊/𝑚2) Tabla 30: Entradas del módulo Horizontal Irradiance Salidas Nombre Descripción Unidades 𝑮𝒃𝒉 Irradiancia directa horizontal Vatios por metro cuadrado (𝑊/𝑚2) 𝑮𝒅𝒉 Irradiancia difusa horizontal Vatios por metro cuadrado (𝑊/𝑚2) Tabla 31: Salidas del módulo Horizontal Irradiance Parámetros Nombre Descripción Unidades 𝑰𝒓𝒓𝑪𝒐𝒎𝒑 Disponibilidad de las componentes de la irradiancia Vatios por metro cuadrado (𝑊/𝑚2) Tabla 32: Parámetros del módulo Horizontal Irradiance Apéndice I: Resumen del modelo 102 102 Módulo: Net DC Array Power Figura 55: Esquema del módulo Net DC Array Power (cuantitativo) Figura 56: Esquema del módulo Net DC Array Power (cualitativo) Entradas Nombre Descripción Unidades 𝑷𝒅𝒄𝒎 Energía producida por el panel Vatios (W) 𝑽𝒅𝒄𝒎 Tensión del panel Voltios (V) 𝑰𝒅𝒄𝒎 Corriente del panel Amperios (A) Tabla 45: Entradas del módulo Net DC Array Power Salidas Nombre Descripción Unidades 𝑷𝒅𝒄𝒏𝒆𝒕 Energía neta del sistema en CC producida Vatios (W) 𝑽𝒅𝒄𝒏𝒆𝒕 Tensión del sistema en CC Voltios (V) 𝑰𝒅𝒄𝒏𝒆𝒕 Corriente del sistema en CC Amperios (A) Tabla 46: Salidas del módulo Net DC Array Power Parámetros Nombre Descripción Unidades 𝑵𝒔𝒖𝒃 Número de submatrices Adimensional 𝑵𝒔𝒕𝒓 Número de cadenas por submatriz Adimensional 𝑵𝒎 Número de módulos por cadena Adimensional 𝑳𝒏𝒂𝒎𝒆𝒑𝒍𝒂𝒕𝒆 Pérdidas nominales (nameplate) Porcentaje (%) 103 103 Software de simulación para la generación fotovoltaica 𝑳𝒅𝒄𝒘𝒊𝒓𝒊𝒏𝒈 Pérdidas por cableado en CC Porcentaje (%) 𝑳𝒎𝒊𝒔𝒎𝒂𝒕𝒄𝒉 Pérdidas por desajuste (mismatch) Porcentaje (%) 𝒅𝒆𝒈𝒓𝒆𝒆_𝑫𝑪 Grado de pérdidas Adimensional Tabla 47: Parámetros del módulo Net DC Array Power Apéndice I: Resumen del modelo 104 104 105 Módulo: Net AC Array Power Figura 57: Esquema del módulo Net AC Array Power (cuantitativo) Figura 58: Esquema del módulo Net AC Array Power (cualitativo) Entradas Nombre Descripción Unidades 𝑷𝒅𝒄𝒏𝒆𝒕 Energía neta del sistema en CC producida Vatios (W) Tabla 48: Entradas del módulo Net AC Array Power Salidas Nombre Descripción Unidades 𝑷𝒂𝒄𝒏𝒆𝒕 Energía neta del sistema en CA producida Vatios (W) 𝑰𝒂𝒄𝒏𝒆𝒕 Corriente del sistema en CA Amperios (A) Tabla 49: Salidas del módulo Net AC Array Power Parámetros Nombre Descripción Unidades 𝑵𝒊𝒏𝒗 Número de inversores Adimensional 𝝋 Factor de potencia Adimensional 𝜼𝒊𝒏𝒗 Eficiencia del inversor Porcentaje (%) 𝑷𝒂𝒄,𝟎 Energía nominal del inversor en CA Vatios (W) 𝑷𝒔,𝟎 Consumo del inversor durante la operación Vatios (W) 𝑷𝒏𝒕,𝟎 Consumo nocturno del inversor Vatios (W) 𝑳𝒕𝒓,𝒍𝒐𝒂𝒅 Pérdidas del transformador en carga Porcentaje (%) 𝑳𝒕𝒓,𝒏𝒐 𝒍𝒐𝒂𝒅 Pérdidas del transformador en vacío Porcentaje (%) 𝑳𝒂𝒄𝒘𝒊𝒓𝒊𝒏𝒈 Pérdidas por cableado en CA Porcentaje (%) Error! Use the Home tab to apply Título 1 to the text that you want to appear here. 106 106 𝒅𝒆𝒈𝒓𝒆𝒆_𝑨𝑪 Grado de pérdidas Adimensional Tabla 50: Parámetros del módulo Net AC Array Power 107 APÉNDICE II: MODELOS REST2 Y ENGERER2 REST2 REST2 es un modelo que calcula la irradiancia para un día despejado utilizando un modelo de dos bandas, ya que este tipo funciona mejor que los modelos de banda ancha más simples. La separación de bandas se sitúa en 0,7μm. La primera banda cubre el ultravioleta y el visible (0,29 a 0,7μm) y se caracteriza por una fuerte absorción por parte del ozono en el ultravioleta y una fuerte dispersión por parte de las moléculas y los aerosoles en toda la banda. La segunda banda cubre el infrarrojo cercano, de 0,7 a 4μm, y se caracteriza por una fuerte absorción por parte del vapor de agua, el dióxido de carbono y otros gases, junto con una dispersión limitada. La parte correspondiente a cada una de las bandas de la constante solar son 𝐺0𝑛1=635,4 𝑊/𝑚2 y 𝐺0𝑛2=709,7 𝑊/𝑚2. Irradiancia directa Para cada una de las dos bandas, la irradiancia directa normal se puede obtener como producto de las transmitancias individuales: 𝐺𝑏𝑛𝑖=𝑇𝑅𝑖𝑇𝑔𝑖𝑇𝑜𝑖𝑇𝑛𝑖𝑇𝑤𝑖𝑇𝑎𝑖𝐺0𝑛𝑖 (II-1) Donde 𝑇𝑅𝑖, 𝑇𝑔𝑖, 𝑇𝑜𝑖, 𝑇𝑛𝑖, 𝑇𝑤𝑖 y 𝑇𝑎𝑖 son las transmitancias de la banda para la dispersión de Rayleigh, la absorción de los gases uniformemente mezclados, la absorción del óxido de nitrógeno, la absorción del vapor de agua y la extinción de los aerosoles, respectivamente. En lugar de usar una sola masa de aire para caracterizar el camino recorrido por los rayos solares a lo largo de la atmósfera, se utilizarán masas ópticas individuales. De forma análoga a las transmitancias, 𝑚𝑅, 𝑚𝑜, 𝑚𝑤 y 𝑚𝑎 se usan para la dispersión de Rayleigh y la absorción de los gases uniformemente mezclados, la absorción del óxido de nitrógeno, la absorción del vapor de agua y la extinción de los aerosoles, respectivamente. Estas masas ópticas dependen del ángulo cenital y vienen dadas por: { 𝑚𝑅=[cos(𝜃𝑧)+0.48353∙𝜃𝑧0.095846∙(96.741−𝜃𝑧)−1.754]−1 𝑚𝑅′=𝑚𝑅(𝑝/𝑝0) 𝑚𝑜=[cos(𝜃𝑧)+1.0651∙𝜃𝑧0.6379∙(101.8−𝜃𝑧)−2.2694]−1 𝑚𝑤=[cos(𝜃𝑧)+0.10648∙𝜃𝑧0.11423∙(93.781−𝜃𝑧)−19203]−1 𝑚𝑎=[cos(𝜃𝑧)+0.16851∙𝜃𝑧0.18198∙(95.318−𝜃𝑧)−1.9542]−1 (II-2) Donde 𝑚𝑅′ es la masa de aire óptica para la dispersión Rayleigh corregida para la presión local con el valor típico 𝑝0=1013.25 𝑚𝑏𝑎𝑟. El resto de las entradas atmosféricas para calcular las transmitancias son la cantidad de ozono, 𝑢𝑜 (atm-cm); la cantidad total de óxido de nitrógeno, 𝑢𝑛 (atm-cm); precipitaciones, w (cm), y los coeficientes de turbiedad de Ångström, 𝛼1, 𝛼2 and 𝛽. La transmitancia para la dispersión Rayleigh para ambas bandas depende solo de la masa óptica de Rayleigh y puede obtenerse como: 𝑇𝑅1=1+1.8169𝑚𝑅′−0.033454𝑚𝑅′2 1+2.063𝑚𝑅′+0.31978𝑚𝑅′2 𝑇𝑅2=1−0.010394𝑚𝑅′ 1−0.00011042𝑚𝑅′2 (II-3) Apéndice II: Modelos REST2 y Engerer2 108 108 Del mismo modo, la dispersión debida a los gases uniformemente mezclados: 𝑇𝑔1=1+0.95885𝑚𝑅′+0.012871𝑚𝑅′2 1+0.96321𝑚𝑅′+0.015455𝑚𝑅′2 𝑇𝑔2=1+0.27284𝑚𝑅′−0.00063699𝑚𝑅′2 1+0.30306𝑚𝑅′ (II-4) La transmitancia para la absorción de ozono en la banda 1 depende de 𝑢𝑜, mientras que es constante en la banda 2: { 𝑓1=𝑢𝑜10.979−8.5421𝑢𝑜 1+2.0115𝑢𝑜+40.189𝑢𝑜2 𝑓2=𝑢𝑜−0.027589−0.005138𝑢𝑜 1−2.4857𝑢𝑜+13.942𝑢𝑜2 𝑓3=𝑢𝑜10.995−5.5001𝑢𝑜 1+1.6784𝑢𝑜+42.406𝑢𝑜2 𝑇𝑜1=1+𝑓1𝑚𝑜+𝑓2𝑚𝑜2 1+𝑓3𝑚𝑜 𝑇𝑜2=1 (II-5) La transmitancia para la absorción del óxido de nitrógeno es similar, pero usando 𝑢𝑛 en su lugar: { 𝑔1=0.17499+41.654𝑢𝑛−2146.4𝑢𝑛2 1+22295𝑢𝑛2 𝑔2=𝑢𝑛−1.2134+59.324𝑢𝑛 1+8847.8𝑢𝑛2 𝑔3=0.17499+61.658𝑢𝑛+9196.4𝑢𝑛2 1+74109𝑢𝑛 𝑇𝑛1=min(1,1+𝑔1𝑚𝑤+𝑔2𝑚𝑤 2 1+𝑔3𝑚𝑤) 𝑇𝑛2=1 (II-6) La transmitancia para la absorción de precipitaciones es diferente para cada banda: 109 109 Software de simulación para la generación fotovoltaica { ℎ1=𝑤0.065445+0.00029901𝑤 1+1.2728𝑤 ℎ2=𝑤0.065687+0.0013218𝑤 1+1.2728𝑤 𝑇𝑤1=1+ℎ1𝑚𝑤 1+ℎ2𝑚𝑤 { 𝑐1=𝑤19.566−1.6506𝑤+1.0672𝑤2 1+5.4248𝑤+1.6005𝑤2 𝑐2=𝑤0.50158−0.14732𝑤+0.047584𝑤2 1+1.1811𝑤+1.0699𝑤2 𝑐3=𝑤21.286−0.39232𝑤+1.2692𝑤2 1+4.8318𝑤+1.412𝑤2 𝑐4=𝑤0.70992−0.23155𝑤+0.096514𝑤2 1+0.44907𝑤+0.75425𝑤2 𝑇𝑤2=1+𝑐1𝑚𝑤+𝑐2𝑚𝑤 2 1+𝑐3𝑚𝑤+𝑐4𝑚𝑤 2 (II-7) Por último, la única transmitancia que queda por calcular es la asociada a los aerosoles. Los aerosoles son pequeñas partículas sólidas y líquidas suspendidas en la atmósfera, como el polvo arrastrado por el viento, las sales marinas, las cenizas volcánicas o la contaminación de las fábricas. Para calcular esta transmitancia se necesitan los coeficientes de turbidez de Ångström, que se relacionan con la linealización de la ley de Ångström: ln(𝜏𝑎𝜆)=ln(𝛽𝑖)−𝛼𝑖ln(𝜆) (II-8) Donde 𝜆 es la longitud de onda y 𝜏𝑎𝜆 es la Longitud Óptica del Aerosol (AOD) para esa longitud de onda en específico. Los coeficientes ln(𝛽𝑖) y −𝛼𝑖 representan la intersección y la pendiente (negativa) de la recta ajustada a los 𝑛 pared de datos (𝜏𝑎𝜆,𝜆), respectivamente. Dicho ajuste puede realizarse de forma independiente para longitudes de onda inferiores y superiores a 700 nm para obtener (𝛼1,𝛽1) y (𝛼2,𝛽2) con las siguientes restricciones: {𝛽1=𝛽0.7𝛼1−𝛼2 𝛽2=𝛽 (II-9) Las parametrizaciones anteriores son válidas al menso para los rangos especificados en la Tabla 1. Con los coeficientes de turbiedad de Ångström, la AOD media en cada banda puede calcularse como: 𝜏𝑎𝑖=𝛽𝑖𝜆𝑒𝑖 −𝛼𝑖 (II-10) Siendo 𝜆𝑒𝑖 la longitud de onda efectiva para toda la banda, la cual se puede calcular como: Apéndice II: Modelos REST2 y Engerer2 110 110 { 𝑑0=0.57664−0.024743𝛼1 𝑑1=0.093942−0.2269𝛼1+0.12848𝛼12 1+0.6418𝛼1 𝑑2=−0.093819+0.36668𝛼1−0.12775𝛼12 1−0.11651𝛼1 𝑑3=𝛼10.15232−0.087214𝛼1+0.012664𝛼12 1−0.90454𝛼1+0.26167𝛼12 𝜆𝑒1=𝑑0+𝑑1𝑢𝑎+𝑑2𝑢𝑎2 1+𝑑3𝑢𝑎2 { 𝑒0=1.183−0.022989𝛼2+0.020829𝛼22 1+0.11133𝛼2 𝑒1=−0.50003−0.18329𝛼2+0.23835𝛼22 1+1.6756𝛼2 𝑒2=−0.50001+1.1414𝛼2+0.0083589𝛼22 1+11.168𝛼2 𝑒3=−0.70003−0.73587𝛼2+0.51509𝛼22 1+4.7665𝛼2 𝜆𝑒2=𝑒0+𝑒1𝑢𝑎+𝑒2𝑢𝑎2 1+𝑒3𝑢𝑎 (II-11) De hecho, 𝜆𝑒𝑖 es básicamente una función de 𝑢𝑎=ln(1+𝑚𝑎𝛽𝑖) La transmitancia de aerosoles en cada banda es, por tanto: 𝑇𝑎𝑖=𝑒−𝑚𝑎𝜏𝑎𝑖 (II-12) Con todas estas transmitancias, la irradiancia directa normal para cada banda se puede obtener con la siguiente ecuación: 𝐺𝑏𝑛𝑖=𝑇𝑅𝑖𝑇𝑔𝑖𝑇𝑜𝑖𝑇𝑛𝑖𝑇𝑤𝑖𝑇𝑎𝑖𝐺0𝑛𝑖 (II-13) La irradiancia directa normal es, por tanto: 𝐺𝑏𝑛=𝐺𝑏𝑛1+𝐺𝑏𝑛2 (II-14) La extinción de los aerosoles es causada principalmente por dispersión y por absorción de la parte restante. La transmitancia de los aerosoles asociada a la dispersión queda expresada como: 𝑇𝑎𝑠𝑖=𝑒−𝑚𝑎𝜔𝑖𝜏𝑎𝑖 (II-15) Donde 𝜔𝑖 es la dispersión individual debida al albedo. Su valor depende del tipo dominante de los aerosoles, aunque para un uso general se pueden utilizar los valores típicos 𝜔1=0.92 y 𝜔2=0.84. Esta transmitancia tendrá un impacto al calcula la irradiancia difusa horizontal, pero no afecta a la irradiancia directa normal. Irradiancia difusa y global El formalismo aquí se basa también en un esquema de dos capas. Se supone que la capa superior es la fuente de toda la dispersión de Rayleigh, así como de toda la absorción de ozono y gases mixtos. Del mismo modo, se asume que la capa inferior es la fuente de toda la dispersión de aerosoles, así como de la absorción de 111 111 Software de simulación para la generación fotovoltaica aerosoles, vapor de agua y dióxido de nitrógeno. Una vez que se ha producido la dispersión en la capa superior, se asume que la irradiancia difusa descendente se comporta como irradiancia directa con una masa de aire efectiva 𝑚′=1.66. Este es el valor de la masa de aire que se utiliza en el cálculo de las transmitancias (indicadas con una prima en lo que sigue) que tratan las absorciones en la capa inferior. La irradiancia difusa incidente sobre un suelo perfectamente absorbente (es decir, con albedo cero) se define como: 𝐺𝑑𝑝𝑖=𝑇𝑜𝑖𝑇𝑔𝑖𝑇𝑛𝑖 ′𝑇𝑤𝑖 ′[𝐵𝑅𝑖(1−𝑇𝑅𝑖)𝑇𝑎𝑖 0.25+𝐵𝑎 𝐹𝑖𝑇𝑅𝑖(1−𝑇𝑎𝑠𝑖 0.25)](𝐺0𝑛𝑖cos(𝜃𝑧)) (II-16) La función 𝐹𝑖 es un factor de corrección introducido para compensar los efectos de la dispersión múltiple y otras deficiencias en el enfoque simple de la transmitancia utilizado aquí. Se puede calcular como: { 𝑔0=3.715+0.368𝑚𝑎+0.036294𝑚𝑎2 1+0.0009391𝑚𝑎2 𝑔1=−0.164−0.72567𝑚𝑎+0.20701𝑚𝑎2 1+0.0019012𝑚𝑎2 𝑔2=−0.052288+0.31902𝑚𝑎+0.17871𝑚𝑎2 1+0.0069592𝑚𝑎2 𝐹1=𝑔0+𝑔1𝜏𝑎1 1+𝑔2𝜏𝑎1 { ℎ0=3.4352+0.65267𝑚𝑎+0.00034328𝑚𝑎2 1+0.034388𝑚𝑎1.5 ℎ1=1.231−1.63853𝑚𝑎+0.20667𝑚𝑎2 1+0.1451𝑚𝑎1.5 ℎ2=0.8889−0.55063𝑚𝑎+0.50152𝑚𝑎2 1+0.14865𝑚𝑎1.5 𝐹2=ℎ0+ℎ1𝜏𝑎2 1+ℎ2𝜏𝑎2 (II-17) 𝐵𝑅1 y 𝐵𝑅2 son las fracciones de dispersión en dirección frontal para la extinción Rayleigh. En ausencia de dispersión múltiple, serían exactamente 0.5 porque las moléculas se dispersan por igual hacia delante y hacia atrás. La dispersión múltiple es despreciable en la banda 2, pero no en la banda 1. Por lo tanto: 𝐵𝑅1=0.5(0.89013−0.0049558𝑚𝑅+0.000045721𝑚𝑅2) 𝐵𝑅2=0.5 (II-18) El factor de dispersión frontal de los aerosoles, 𝐵𝑎, viene dado por: 𝐵𝑎=1−𝑒−0.6931−1.8326cos(𝜃𝑧) (II-19) La irradiancia difusa de banda ancha que incide sobre una superficie horizontal e idealmente negra viene dada por: 𝐺𝑑𝑝=𝐺𝑑𝑝1+𝐺𝑑𝑝2 (II-20) Sin embargo, hay que añadir una contribución de retrodispersión debido a la interacción entre la superficie terrestre reflectante y las capas de dispersión de la atmósfera. Esta contribución suele ser pequeña (por ejemplo, <10% de 𝐺𝑔𝑝), pero puede llegar a ser mucho más significativa en regiones nevadas. El albedo del suelo considerado aquí, 𝜌𝑔, se refiere a una media sobre una amplia zona de 5-50 km de radio alrededor del lugar analizado. Para cada banda, el albedo del cielo, 𝜌𝑠𝑖, se obtiene en función de 𝛼𝑖 y 𝛽𝑖: Apéndice III: Códigos en Julia 118 118 end_line=192; doc=1; end elseif (ARGS[1]=="Module_Inverter_Losses_Qualitative2") if (length(ARGS)==15) doc=0; else start_line=194; end_line=214; doc=1; end elseif (ARGS[1]=="Inverter_Qualitative") if (length(ARGS)==26) doc=0; else start_line=216; end_line=247; doc=1; end elseif (ARGS[1]=="Inverter_Losses_Qualitative") if (length(ARGS)==21) doc=0; else start_line=249; end_line=275; doc=1; end elseif (ARGS[1]=="Losses_Qualitative") if (length(ARGS)==27) doc=0; else start_line=277; end_line=309; doc=1; end else start_line=1; end_line=309; doc=1 end else start_line=1; end_line=309; doc=1 end # If the input parameters are correct, compute the power output if (doc==0) 119 119 Software de simulación para la generación fotovoltaica Data=Load_Data(Location_data, Module_data, Park_data, Inverter_data, Losses_data); doc=Data[1]; msg=Data[2]; Location=Data[3]; Module=Data[4]; Park=Data[5]; Inverter=Data[6]; Losses=Data[7]; IrrComp=Data[8]; # If there was a problem while loading the data, show message explaining it if (doc==1) println("",msg); else println("",msg); # Computation of Sun Angles Data=Sun_Angles(Location); alpha=Data[1]; azimuth=Data[2]; omega=Data[3]; Z=Data[4]; M=Data[5]; A=Data[6]; ts=Data[7]; jdoy=Data[8]; G0h=Data[9]; EOT=Data[10]; ## Extraction of Irradiance Components after entering atmosphere Data=Weather_File_Data(Location, Z, M, IrrComp); Tamb=Data[1]; Wspeed=Data[2]; if IrrComp BHI=Data[3]; DHI=Data[4]; GHI=BHI+DHI; Gbh=BHI; Gdh=DHI; else GHI=Data[3]; Data=REST2(Z); Ghc=Data; Data=Engerer2(Location, jdoy, GHI, Ghc, G0h, EOT); Gbh=Data[1]; Gdh=Data[2]; end ## Extraction of Irradiance Components on Tilted Surface Apéndice III: Códigos en Julia 120 120 Data=Tilted_Irradiance(M, Park, azimuth, Z, Gbh, Gdh, Location); Gbi=Data[1]; Gdi=Data[2]; Gri=Data[3]; AOI=Data[4]; beta_s=Data[5]; ## Irradiance Components on Tilted Surface after Shading and Soiling Data=Shading_Soiling(Gbi, Gdi, Gri); Gb=Data[1]; Gd=Data[2]; Gr=Data[3]; ## Irradiance absorbed by the PV Module Data=Module_Irradiance(M, Park, AOI, Gb, Gd, Gr, Z, A, beta_s); Geff=Data[1]; tau_alpha=Data[2]; ## Power Output of a single PV Module Data=Module_Power(M, Module, Wspeed, tau_alpha, Tamb, Geff); Tc=Data[1]; P=Data[2]; V=Data[3]; I=Data[4]; ## Gross DC Power Output of the array Data=Gross_DC_Power(M, P, V, I, Park); Pdc_gross=Data[1]; Vdc=Data[2]; Idc=Data[3]; Pdc_gross_n=Data[4]; ## Net DC Power Output of the array Data=Net_DC_Power(M, Pdc_gross_n, Losses, Park); Pdc_net=Data; ## Gross AC Power Output of the array Data=Gross_AC_Power(M, Inverter, Pdc_net); Pac_gross=Data; ## Net AC Power Output of the array Data=Net_AC_Power(Pac_gross, Losses, Inverter); Pac_net=Data; # Print the annual DC and AC energy EnergyDC=round(sum(Pdc_net)/1e3); @printf "The DC energy produced throughout a year is %f kWh/year\n" EnergyDC; 121 121 Software de simulación para la generación fotovoltaica EnergyAC=round(sum(Pac_net)/1e3); @printf "The AC energy produced throughout a year is %f kWh/year\n" EnergyAC; # Save power ouput Save_Data(ts, Z, alpha, azimuth, Gbh, Gdh, AOI, Geff, Tc, Pdc_gross_n, Pac_net); println("\nData saved successfully"); end # If the parameters were not intruduced correctly else Show_Doc(start_line, end_line); end Apéndice III: Códigos en Julia 122 122 Show_Doc.jl Script Show_Doc.jl Descripción Función que muestra la documentación si los parámetros no se introdujeron de forma correcta function Show_Doc(start_line, end_line) if (start_line!=0) f = open("Documentation.txt"); lines = readlines(f); lines = lines[start_line:end_line]; if (start_line>1) println("\nThe number of arguments does not match. Please, check the list below"); end for l in lines println("$l") end close(f); else println("\nNo matching method. Please, enter help as argument to show documentation"); end end 123 123 Software de simulación para la generación fotovoltaica Load_Data.jl Script Load_Data.jl Descripción Función que inicializa las estructuras de acuerdo con los parámetros introducidos function Load_Data(Location_data, Module_data, Park_data, Inverter_data, Losses_data) # Flag in case the data was introduced in wrong format doc=0; msg=""; # Structures' initialization Location=Location_data(); Module=Module_data(); Park=Park_data(); Inverter=Inverter_data(); Losses=Losses_data(); # Location data directly from arguments Location.filename=ARGS[2]; Location.albedo=parse(Float64,ARGS[3]); if (Location.albedo<0 || Location.albedo>1) doc=1; msg="\nAlbedo must be between 0 and 100\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end # List with all the modules of the database module_list=readdlm("$(pwd())/../Modules/CEC Modules.csv",','); # Flag indicating if the irradiance components are available IrrComp=parse(Int,ARGS[end]); if (IrrComp==1) IrrComp=true; else IrrComp=false; end # Switch depending on the slected mode if (ARGS[1]=="Quantitative") # Module parameters Module.material=ARGS[31]; Module.Vmpp=parse(Float64,ARGS[4]); Module.Impp=parse(Float64,ARGS[5]); Module.Kisc=parse(Float64,ARGS[6])/100; if (Module.Kisc<0 || Module.Kisc>1) doc=1; msg="\nThe short circuit coefficient must be between 0 and 100\n"; Apéndice III: Códigos en Julia 124 124 return doc, msg, Location, Module, Park, Inverter, Losses; end Module.Kvoc=parse(Float64,ARGS[7])/100; if (Module.Kvoc<-1 || Module.Kvoc>1) doc=1; msg="\nThe open circuit coefficient must be between -100 and 100\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end Module.Kp=parse(Float64,ARGS[8])/100; if (Module.Kp<-1 || Module.Kp>1) doc=1; msg="\nThe power coefficient must be between -100 and 100\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end Module.Am=parse(Float64,ARGS[9]); Module.Tnoct=parse(Float64,ARGS[10])+273.15; # Park parameters Park.Track_opt=parse(Int,ARGS[11]); if (Park.Track_opt<0 || Park.Track_opt>2) doc=1; msg="\nThe tracking options are 0: Fixed, 1: Azimuth tracking and 2: Two- Axis tracking\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end Park.beta_0=parse(Float64,ARGS[12])*pi/180; if (Park.beta_0<0 || Park.beta_0>pi/2) doc=1; msg="\nThe tilt angle must be between 0 and 90\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end Park.azimuth_0=parse(Float64,ARGS[13])*pi/180; if (Park.azimuth_0<0 || Park.azimuth_0>pi) doc=1; msg="\nThe orientation angle must be between 0 and 180\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end Park.nsub=parse(Int,ARGS[14]); Park.nstr=parse(Int,ARGS[15]); Park.nmod=parse(Int,ARGS[16]); # Inverter parameters 125 125 Software de simulación para la generación fotovoltaica Inverter.Pdc=parse(Float64, ARGS[17]); Inverter.Pac=parse(Float64, ARGS[18]); if (Inverter.Pac > Inverter.Pdc) doc=1; msg="\nThe AC power of the inverter cannot be greater than the DC power of the invert er\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end Inverter.efficiency=Inverter.Pac/Inverter.Pdc; Inverter.Ninv=ceil((Module.Vmpp*Module.Impp*Park.nsub*Park.nstr*Park.nmod)/Inverter.Pdc); Inverter.Vdc_max=parse(Float64, ARGS[19]); if (isdefined(Module, 6) && Inverter.Vdc_max < (Module.Voc*Park.nmod)) doc=1; msg="\nThe maximum DC voltage of the inverter is not greater than Open Circuit Voltag e of a single string\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end Inverter.Vmppt_min=parse(Float64, ARGS[20]); if (Inverter.Vmppt_min > (Module.Vmpp*Park.nmod)) doc=1; msg="\nThe minimum MPPT voltage is not lower than maximum voltage of a single string\ n"; return doc, msg, Location, Module, Park, Inverter, Losses; end Inverter.Vmppt_max=parse(Float64, ARGS[21]); if (Inverter.Vmppt_max < (Module.Vmpp*Park.nmod)) doc=1; msg="\nThe maximum MPPT voltage is not greater than maximum voltage of a single strin g\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end Inverter.Ps0=parse(Float64, ARGS[22]); Inverter.Pnt0=parse(Float64, ARGS[23]); # Losses parameters Losses.soiling=parse(Float64, ARGS[24])/100; if (Losses.soiling < 0 || Losses.soiling > 1) doc=1; msg="\nThe soiling loss must be between 0 and 100\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end Losses.nameplate=parse(Float64, ARGS[25])/100; if (Losses.nameplate < 0 || Losses.nameplate > 1) doc=1; Apéndice III: Códigos en Julia 126 126 msg="\nThe nameplate loss must be between 0 and 100\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end Losses.mismatch=parse(Float64, ARGS[26])/100; if (Losses.mismatch < 0 || Losses.mismatch > 1) doc=1; msg="\nThe mismatch loss must be between 0 and 100\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end Losses.DC_wiring=parse(Float64, ARGS[27])/100; if (Losses.DC_wiring < 0 || Losses.DC_wiring > 1) doc=1; msg="\nThe DC wiring loss must be between 0 and 100\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end Losses.AC_wiring=parse(Float64, ARGS[28])/100; if (Losses.AC_wiring < 0 || Losses.AC_wiring > 1) doc=1; msg="\nThe AC wiring loss must be between 0 and 100\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end Losses.transformer_load=parse(Float64, ARGS[29])/100; if (Losses.transformer_load < 0 || Losses.transformer_load > 1) doc=1; msg="\nThe transformer load loss must be between 0 and 100\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end Losses.transformer_noload=parse(Float64, ARGS[30])/100; if (Losses.transformer_noload < 0 || Losses.transformer_noload > 1) doc=1; msg="\nThe transformer no load loss must be between 0 and 100\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end elseif (ARGS[1]=="Module_Qualitative1") # Module parameters Module.name=ARGS[4]; index=find(x -> x==Module.name, module_list); Module.material=module_list[index,22][1]; Module.Tnoct=module_list[index,4][1]+273.15; Module.Am=module_list[index,5][1]; Module.Ns=module_list[index,6][1]; 127 127 Software de simulación para la generación fotovoltaica Module.Isc=module_list[index,7][1]; Module.Voc=module_list[index,8][1]; Module.Impp=module_list[index,9][1]; Module.Vmpp=module_list[index,10][1]; Module.Kisc=(module_list[index,11][1]/module_list[index,7][1]); Module.Kvoc=(module_list[index,12][1]/module_list[index,8][1]); Module.aref=module_list[index,13][1]; Module.Iph_ref=module_list[index,14][1]; Module.Io_ref=module_list[index,15][1]; Module.Rs=module_list[index,16][1]; Module.Rsh_ref=module_list[index,17][1]; Module.Kp=module_list[index,19][1]/100; # Park parameters Park.Track_opt=parse(Int,ARGS[5]); if (Park.Track_opt<0 || Park.Track_opt>2) doc=1; msg="\nThe tracking options are 0: Fixed, 1: Azimuth tracking and 2: Two- Axis tracking\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end Park.beta_0=parse(Float64,ARGS[6])*pi/180; if (Park.beta_0<0 || Park.beta_0>pi/2) doc=1; msg="\nThe tilt angle must be between 0 and 90\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end Park.azimuth_0=parse(Float64,ARGS[7])*pi/180; if (Park.azimuth_0<0 || Park.azimuth_0>pi) doc=1; msg="\nThe orientation angle must be between 0 and 180\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end Park.nsub=parse(Int,ARGS[8]); Park.nstr=parse(Int,ARGS[9]); Park.nmod=parse(Int,ARGS[10]); # Inverter parameters Inverter.Pdc=parse(Float64, ARGS[11]); Inverter.Pac=parse(Float64, ARGS[12]); if (Inverter.Pac > Inverter.Pdc) doc=1; msg="\nThe AC power of the inverter cannot be greater than the DC power of the invert er\n"; return doc, msg, Location, Module, Park, Inverter, Losses; end