Desarrollo mediante UNITY 3D del modelo para un gemelo digital de un captador solar tipo Fresnel
Abstract
Este trabajo desarrolla un modelo virtual de un Colector Solar Fresnel (FSC) comercial instalado en una planta de refrigeración por absorción. El objetivo es sentar las bases para la creación de un Digital Twin (DT) de última generación, utilizando los modelos validados previamente creados bajo el proyecto DENiM y adaptándolos al lenguaje de programación C# y al entorno de simulación Unity3D. La idea es que el gemelo digital ayude en la optimización de arranque y apagado, así como en el diseño de control y fiabilidad. El modelo seleccionado utiliza un enfoque fenomenológico basado en Ecuaciones Diferenciales Parciales (EDP) y la identificación de parámetros. Los conjuntos de datos de medición disponibles abarcan 25 días de operación con un tiempo de muestreo de 20 segundos, que, tras la eliminación de valores atípicos, filtrado y tratamiento, resultaron en 108416 muestras. La validación se realizó considerando seis días de operación por separado.
Full text
Proyecto Fin de Carrera Ingeniería de Telecomunicación Formato de Publicación de la Escuela Técnica Superior de Ingeniería Autor: F. Javier Payán Somet Tutor: Juan José Murillo Fuentes Dep. Teoría de la Señal y Comunicaciones Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2013 Trabajo Fin de Grado Grado en Ingeniería Electrónica, Robótica y Mecatrónica Desarrollo mediante UNITY 3D del modelo para un gemelo digital de un captador solar tipo Fresnel Autor: Pablo Carmona Benito Tutor: Juan Manuel Escaño González Cotutor: William David Chicaiza Salazar Dpto. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2023
Trabajo Fin de Grado Grado en Ingeniería Electrónica, Robótica y Mecatrónica Desarrollo mediante UNITY 3D del modelo para un gemelo digital de un captador solar tipo Fresnel Autor: Pablo Carmona Benito Tutor: Juan Manuel Escaño González Profesor Titular Cotutor: William David Chicaiza Salazar Dpto. Ingeniería de Sistemas y Automática Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2023
Trabajo Fin de Grado: Desarrollo mediante UNITY 3D del modelo para un gemelo digital de un captador solar tipo Fresnel Autor: Pablo Carmona Benito Tutor: Juan Manuel Escaño González Cotutor: William David Chicaiza Salazar El tribunal nombrado para juzgar el trabajo arriba indicado, compuesto por los siguientes profesores: Presidente: Vocal/es: Secretario: acuerdan otorgarle la calificación de: El Secretario del Tribunal Fecha:
Resumen Es te trabajo desarrolla un modelo virtual de un Colector Solar Fresnel (FSC) comercial instalado en una planta de refrigeración por absorción. El objetivo es sentar las bases para la creación de un Digital Twin (DT) de última generación, utilizando los modelos validados previamente creados bajo el proyecto DENiM y adaptándolos al lenguaje de programación C# y al entorno de simulación Unity3D. La idea es que el gemelo digital ayude en la optimización de arranque y apagado, así como en el diseño de control y fiabilidad. El modelo seleccionado utiliza un enfoque fenomenológico basado en Ecuaciones Diferenciales Parciales (EDP) y la identificación de parámetros. Los conjuntos de datos de medición disponibles abarcan 25 días de operación con un tiempo de muestreo de 20 segundos, que, tras la eliminación de valores atípicos, filtrado y tratamiento, resultaron en 108416 muestras. La validación se realizó considerando seis días de operación por separado. I
Abstract Th is work presents the development of a virtual model for a commercial Fresnel Solar Collector (FSC) that is installed in an absorption cooling plant. The main objective is to establish the foundation for creating a state-of-the-art Digital Twin (DT). The DT will be based on validated models that were previously created under the DENiM project. These models will be adapted for the C# programming language and the Unity3D simulation environment. The purpose of the digital twin is to assist in optimizing start-up and shut-down processes and to improve control system reliability. The selected model employs a phenomenological approach, using Partial Differential Equations (PDE) and parameters identification. To validate the model, a dataset of measurement data spanning 25 days of operation was collected, with a sampling time of 20 seconds. After removing outliers and applying filtering and data treatment techniques, the dataset was reduced to 108,416 samples. The validation process involved six separate operating days. III
2Capítulo 1. Introducción • Toma de decisiones y optimización: proporcionan información valiosa para la toma de decisiones y la optimización de procesos. Basados en análisis y simulaciones, los DT ayudan a identificar oportunidades de mejora, predecir fallos o problemas, y optimizar el rendimiento del sistema físico. En resumen, los DT son réplicas digitales en tiempo real de sistemas físicos, que permiten el análisis, la simulación y la toma de decisiones basadas en datos y escenarios del mundo real. Estos modelos digitales se actualizan continuamente utilizando datos obtenidos en tiempo real para reflejar el estado y el comportamiento de su contraparte física. 1.1.2 Los gemelos digitales en la actualidad Hoy en día, la ingeniería y la fabricación utilizan los DT predominantemente para proporcionar representaciones virtuales precisas de objetos y simulaciones de procesos productivos. Las aplicaciones de los DT en la gestión de operaciones y cadenas de suministro, especialmente en términos de trazabilidad de operaciones, mantenimiento de transporte, asistencia remota, visualización de activos y personalización de diseño, se revisan en publicaciones relacionadas [ 14 ]. La tecnología está lista para cumplir con sus numerosas promesas en otras industrias, como automotriz, aeroespacial, construcción, agricultura, minería, servicios públicos, comercio minorista, atención médica, militar, recursos naturales y seguridad pública. Daremos algunos ejemplos como botón de muestra del desarrollo y la implantación de los DT en la actualidad. •SoFi Stadium de Los Angeles Rams tiene un DT diseñado por Willow Inc. • McKinsey diseñó y produjo el DT para la nave de vela del equipo Neozelandés en la Copa de América. • fábricas, centros de distribución y logística, llevan años empleando DT para la optimización de sus procesos. • el Foro Económico Mundial, junto con la Academia China de las Tecnologías de la Información y la Comunicación, llevan dos años embarcados en el proyecto de facilitar el desarrollo de DT para ciudades completas. • el concepto de Gemelo Digital nació en los años 60 en la NASA, en el contexto de la misión Apollo; hasta el día de hoy, usan la tecnología de los DT para sus aeronaves. SpaceX también la pone en práctica para su "Dragon", de la mano de SIEMENS. • un grupo de investigadores están desarrollando un DT del sistema inmune, con el objetivo de mejorar el tratamiento de enfermedades y detectar posibles complicaciones antes de que se den. 1.2 Captadores Fresnel: estado del arte 1.2.1 La tecnología Fresnel Las plantas solares Fresnel son un tipo de sistema de energía solar concentrada (CSP, por sus siglas en inglés) que utiliza espejos para concentrar la luz solar en un receptor y generar calor. Un receptor, compuesto por tubos y rodeado por un reflector secundario, contiene el fluido que transporta el calor. Los concentradores solares, formados por espejos, dirigen la radiación solar hacia el receptor. Un sistema de control, como un autómata programable, orienta los concentradores para seguir la trayectoria del sol. El sistema también incluye motores para mover los concentradores y una estructura que sostiene todos los componentes. Así, se logra generar calor para su uso en una máquina de absorción o, alternativamente, vapor de agua para producir electricidad a través de una turbina.
1.2 Captadores Fresnel: estado del arte 3 La idea del sistema Fresnel es lograr un efecto similar al conseguido mediante un captador cilindro-parabólico (Figura 1.1), pero ganando en algunos aspectos. Figura 1.1 Principio de funcionamiento de un captador cilindro-parabólico. Este tipo de concentradores pueden considerarse como un particionamiento (aproximación de primer orden, lineal) de un perfil parabólico, y los segmentos son colocados en una superficie plana. Figura 1.2 Principio de funcionamiento de un captador Fresnel. Esta tecnología ofrece varias ventajas, incluyendo un coste más bajo en comparación con otras tecnologías CSP, un diseño e instalación más sencillos y una mayor eficiencia de conversión de energía solar a térmica. El uso de una distribución de lentes planas, que funcionan en su conjunto como una lente Fresnel, permite un ángulo de aceptación más amplio, lo que facilita el seguimiento solar y aumenta la recolección de energía. Estas plantas han sido estudiadas e implementadas en diversas regiones del mundo, demostrando su potencial como una solución de energía renovable viable y eficiente. Por ejemplo, la investigación realizada por He et al. (2020) resaltó la mejora en el rendimiento logrado por las plantas solares Fresnel en términos de eficiencia térmica y generación de electricidad. Además, estudios como los realizados por Zoulias y Lymberopoulos (2009) y Wang et al. (2018) han explorado los aspectos técnicos y la viabilidad económica de las plantas solares Fresnel, demostrando su papel prometedor en el futuro de la generación de energía solar. Para ver con más detalle el funcionamiento de esta tecnología, acudir a Zhu et al. (2014).
4Capítulo 1. Introducción 1.3 Objetivos del proyecto El objetivo perseguido con este trabajo es el desarrollo del modelo físico para un gemelo digital de un captador solar mediante Unity 3D. Con ello se pretende, en primer lugar, contribuir al depósito de conocimientos de la Escuela, y, en particular, del Departamento, sobre el modelado de DT. Aunque no está oficialmente vinculado a él, de alguna manera este trabajo ha sido motivado por el proyecto DENiM. Como segundo objetivo, se busca dar las herramientas para, en el momento oportuno, volver a poner en funcionamiento la planta solar Fresnel de la ETSI. Finalmente, y de modo más concreto, se procura proporcionar un modelo eficiente programado en C# e integrado en Unity3D, dos tecnologías de uso comercial más extendido que Matlab, programa en que se basaba el modelo anterior. Precisamos que ha habido dos actividades principales en el transcurso del proyecto: una, el modelado en 3D del sistema Fresnel. Aunque ya existía un modelo en 3D en Unity3D, era parcial. Se ha pretendido dar un modelo que abarca todo el captador. Por cuestiones de tiempo y falta de fichas técnicas y planos más precisos, no se ha modelado hasta el detalle, dejando ese punto a posteriores trabajos de investigación. La segunda actividad ha sido la "traducción" y adaptación del modelo desde Matlab a C#, y hacerlo compatible con la plataforma Unity3D. 1.4 Descripción de la planta solar de tipo Fresnel de la ETSI La Escuela Técnica Superior de Ingeniería de Sevilla (ETSI) cuenta con una instalación de refrigeración solar (PRS), que utiliza concentradores lineales de Fresnel para aprovechar la radiación solar y generar frío para sistemas de climatización. Esta tecnología fue desarrollada a partir de un proyecto presentado por la empresa "Gas Natural" en colaboración con la ETSI y la Corporación Tecnológica de Andalucía en 2006. La instalación piloto, puesta en marcha en 2008, demostró la viabilidad técnica de las aplicaciones de absorción de doble efecto con concentradores de Fresnel. Figura 1.3 Vista general de la planta. El objetivo principal era establecer parámetros técnicos, económicos y medioambientales para el diseño de instalaciones de este tipo en el sector terciario. La PRS de la ETSI consta de 352 m2 de concentradores lineales de Fresnel que calientan un caudal de agua para generar energía térmica utilizada en calefacción y refrigeración. La instalación incluye un sistema de almacenamiento de energía para manejar los excesos de potencia generados. Los estudios posteriores confirmaron la fiabilidad y seguridad de los captadores solares de Fresnel, y se observó la necesidad de implementar dichos sistemas de almacenamiento para maximizar su eficiencia y rendimiento.
1.4 Descripción de la planta solar de tipo Fresnel de la ETSI 5 Tabla 1.1 Características de la planta Fresnel de la ETSI. Parámetro (símbolo) Valor Unidad Área total 512 m2 Área total del espejo (At) 352 m2 Longitud del absorbedor 64 m Altura del absorbedor 4 m Longitud unitaria del espejo 4 m Ancho unitario del espejo 0.5 m Abertura total del espejo (G) 5.5 m Filas de espejos 11 - Número de espejos 176 - Masa específica del tubo absorbedor (Pm) 8027 kg/m3 Diámetro exterior del tubo absorbedor 0.07 m Diámetro interior del tubo absorbedor 0.066 m Calor específico del tubo absorbedor (cm) 500 J/(kg °C) Abertura de la cavidad receptora 0.3 m Fluido portador del calor Agua - Temperatura nominal 180 °C Presión nominal 13 bar Potencia térmica nominal 120 kW En nuestra planta, cada espejo concentrador es movido por un motor, a su vez controlado por el PLC conectado a la planta, que tiene el objetivo de concentrar la radiación solar en el tubo receptor, para lograr capturar la mayor cantidad de energía posible. El espejo reflector secundario asiste a los reflectores primarios concentrando los haces de luz que se desvían en el tubo, para no desperdiciar energía. Figura 1.4 Dibujo de concepto de la planta con algunas medidas. Las coordenadas de la instalación son 6°de longitud hacia el oeste y 37,41°de latitud. Su orientación es 12°3’1" de la dirección sur. Para una descripción más detallada del sistema Fresnel instalado en la ETSI, ver González Parrilla (2019).
6Capítulo 1. Introducción Figura 1.5 Vista general de la planta.
2 Plataforma Unity 3D Co mo hemos mencionado, una parte fundamental de un DT es la representación virtual del sistema modelado. La tendencia es usar representaciones 3D, ya que la tecnología actual permite construir modelos precisos y mostrarlos e interactuar con ellos sin demasiado esfuerzo. Hemos diseñado un modelo virtual en 3D del captador Fresnel que queremos modelar. Para ello empleamos el software Unity 3D. La elección del mismo ha sido motivada por la decisión del grupo de investigación de emplear esa herramienta en proyectos anteriores, debido a su adecuación al propósito de la investigación. Para un estudio detallado de los softwares disponibles y la justificación del uso de Unity 3D, se puede consultar el trabajo de González (2021) [4]. 2.1 Utilidad Unity 3D es una plataforma de desarrollo de videojuegos en 2D y en 3D. Es compatible con los grandes sistemas operativos, y permite exportar los proyectos a una gran variedad de soportes: PC, Smartphone, Televisión, aparatos de realidad virtual, consolas de videojuegos... La motivación principal para usar Unity 3D es su motor físico. Permite la simulación realista del comportamiento de la materia. Sobre todo, trata de modos muy deseables la interacción gravitatoria, la interacción por colisión entre dos o más cuerpos, y la interacción de las superficies con la luz. Permite también tener en cuenta momentos e inercias, y facilita la programación del comportamiento físico de un objeto. 2.2 Breve introducción a Unity 3D Para obtener el software se puede acceder a la página https://unity.com/ y completar el registro, para posteriormente descargarlo. Unity funciona por proyectos, así que, una vez instalado en nuestro ordenador, podemos abrir el programa y seleccionar, en la interfaz de UnityHub, la opción de crear un nuevo proyecto. Para nuestro caso, se trata de un proyecto 3D. Una vez hecho esto, se nos presentará la interfaz de Unity 3D. Veremos los elementos principales del software para orientarnos en el proceso de creación del modelo: • Scene: pretende ser un nivel dentro del juego. En un proyecto pueden darse varias escenas, aunque nosotros, en esta ocasión, trabajaremos con una sola escena. Encontramos diversos gameobjects. Por defecto, aparecen un foco de luz y una cámara (o punto de vista). En el juego se muestra la parte de la escena que esté en el rango visual de la cámara. • Game object: elemento básico en el que se pueden incluir componentes que modifiquen sus propiedades. 7
8Capítulo 2. Plataforma Unity 3D Figura 2.1 Pantalla de selección del tipo de proyecto en Unity Hub. • Camera: su posición y orientación, modificables, determinan la vista de juego o game view: su campo de visión es lo que se renderiza. Viene dado por defecto. • Light: es la fuente de luz por defecto de la escena. Sus parámetros son accesibles ymodificables. Se pueden añadir puntos de luz de tipos variados. • Component: definen propiedades, modos de comportarse o acciones de los game objects. Hay un buen número de tipos de componentes en Unity, además de los que se pueden crear por cuenta propia. • Asset: elementos del proyecto, que pueden ser texturas, objetos importados, prototipos, scripts... • Prefab: se trata de un game object guardado en algún directorio del proyecto, habitualmente un modelo en 2D o 3D, que se puede instanciar, generando copias dependientes del prefab (los cambios se propagan del prefab a sus instancias, no al revés). • Transform: contiene la posición, rotación y escalado del game object al que pertenece. Puede ser absoluto (respecto del "mundo") o relativo, si el objeto en cuestión es "hijo" de otro. • Script: archivo de código que permite acceder a los elementos del proyecto y modificar sus propiedades, definir su comportamiento... también en tiempo de juego (durante la simulación, en nuestro caso). Nosotros programamos en C# (C-Sharp). • Tag: quizá no es de los elementos fundamentales, pero interesa para nuestro proyecto. Se trata de lo que su nombre indica: una etiqueta que marca un game object, y sirve para buscar en la escena, interactuar con el objeto mediante código... 2.2.1 Estructura de un script C# para Unity 3D La estructura básica de un script común para Unity 3D escrito en C# es la siguiente: •Inclusión de librerías o paquetes. •Declaración de la clase: nombre y tipo de la clase. • Start / Awake: dentro de la clase, definimos el comportamiento inicial del elemento al que va ligado el script, en su caso, o de algún determinado fenómeno o condición en el juego. El programa incluido dentro de la función Start / Awake, se ejecuta al iniciarse el juego o la simulación.
2.2 Breve introducción a Unity 3D 9 • Update / FixedUpdate: dentro de la clase se define también su comportamiento en los sucesivos frames del juego. Para generar la imagen y las acciones, se ejecuta el código contenido en esta sección cada tiempo determinado 1 2using UnityEngine; 3using System; 4 5 6class ejemplo 7{ 8 9// Use this for initialization 10 void Start () { 11 12 // Acciones que se ejecutan una vez al principio del juego / simulación,→ 13 14 } 15 16 // Update is called once per frame 17 void Update () { 18 19 // Acciones que se ejecutan en cada "frame" del juego; cada vez que se refresca,→ 20 21 } 22 23 } 24 Para el manejo del software hemos recurrido a los cursos que ofrecen en la plataforma UnityLearn.
3 Modelo de la planta El objetivo del proyecto que nos ocupa ha evolucionado notablemente en el transcurso del tiempo, hasta acabar concretándose en el presente desarrollo. Al principio, cuando empezamos a pensar sobre ello, la idea era hacer un SCADA de la planta de absorción solar del tejado de la Escuela de Ingenieros de la Universidad de Sevilla, sirviéndonos para ello de herramientas más rudimentarias de modelado 3D, y centrándonos en la manipulación de parámetros y entradas de la planta empleando como interfaz una pantalla táctil que instalaríamos en la azotea. Para ello dedicamos unas semanas a la formación en SCADA con Arduino, Python y HTML, con idea de generar un SCADA online que permitiese la modificación remota de parámetros. Desistimos pronto de aquel camino, porque la pantalla estaba deteriorada y soportaba un software anticuado. El siguiente enfoque fue contribuir a la realización de un DT de la planta, empleando Unity 3D, pero sin incorporar el modelo matemático (limitándonos a un modelado 3D). Aquel concepto fue enriqueciéndose y decidimos concretarlo de la siguiente forma: un modelo de fluidos de la planta, que modelase, mediante plugins de Unity, el comportamiento del agua caloportadora que recorre el tubo del captador en la planta. Este planteamiento requirió un período de investigación para encontrar el mejor modo de simular fluidos en Unity 3D. Dimos con un paquete que parecía adecuado: Zibra Liquids de Zibra AI. Se trata de un plugin con versión gratuita de prueba: tras investigar la dinámica del fluido y el aspecto de la simulación, supimos que el grupo de investigación había comprado ya, por otros motivos, dos licencias. De manera que accedimos a una y comenzamos a hacer pruebas para cerciorarnos de la completa adecuación al propósito. Sin embargo, topamos con una dificultad: el paquete no ha conseguido modelar satisfactoriamente la dinámica del agua por un conducto, teniendo en cuenta presiones, cambios de caudal... Sí que permite la simulación y renderización de fluido libre o el efecto de la gravedad. Pero al interactuar el fluido con conductos cerrados, en los que entra en juego el transporte de masas (revisar si es correcto, y, si no, corregir la idea), el plugin no interactua adecuadamente con la superficie del conducto, no detecta correctamente las normales. Tras acudir reiteradamente a varias plataformas de resolución de problemas que posee la compañía, desistimos de usar el plugin. Así, acudimos al departamento de Ingeniería de Fluidos de la Escuela, para informarnos sobre posibles modos de llevar a cabo la simulación. La conversación no fue muy animante, pero salimos con la idea de intentar generar un modelo de partículas que permitiese simular y renderizar el fluido en Unity 3D. Con esto en mente, encontramos dos publicaciones interesantes por Jos Stam [ 11 ] [ 12 ] sobre el modelado de fluidos para videojuegos. Sin embargo, viendo la complejidad que suponía, decidimos optar otra vía menos escabrosa, dejando este proyecto a futuros desarrollos. De este modo acabamos encauzando el trabajo hacia donde estamos ahora. 11
18 Capítulo 3. Modelo de la planta caracteristicasdelfluido.m Simplemente toma la temperatura actual del fluido en cada punto del circuito y el caudal, y con eso calcula la densidad del fluido, su calor específico, su número de Prandtl y el coeficiente de transmisión de calor por convección metal-fluido. perdidasmetal.m Halla el coeficiente de pérdidas metal-ambiente, que es distinto según el tramo de tubo/tubería que se recorra. interpola.m Esta es una función auxiliar que sirve para interpolar los datos como hemos explicado anteriormente. Emplea interpolación lineal. Error y tiempo de simulación El modelo en que se inspira nuestro trabajo da unos resultados bastante aceptables, teniendo en cuenta el rango de operación de la temperatura de salida de la planta (salida del tubo receptor): un error máximo de -2.82ºC en el caso presentado en la figura 3.8 para un rango de operación entre 64ºC y 160ºC. Figura 3.7 Resultados del modelo PDE en Matlab para las muestras del día 29 de junio de 2009. Así pues, creemos acertado basarnos en este modelo para construir la parte física del nuestro. Por otro lado, el modelo en Matlab tarda en simular un set de datos una media de 13 segundos con algunas centésimas. Se trata de un tiempo de simulación estupendo para comprender varias horas de funcionamiento de la planta y dada la finalidad del experimento.
3.2 Modelo de parámetros distribuidos 19 Figura 3.8 Error del modelo PDE en Matlab para las muestras del día 29 de junio de 2009. 3.2.2 Estructura del modelo en C# Trasladar el modelo de ecuaciones en diferencias parciales a C# requiere una adaptación considerable. El ventajoso trato de matrices de Matlab no existe en C#, por lo que recorrer los vectores de datos y operar con ellos no es tan sencillo. Por otro lado, la discretización del modelo, que heredamos del propuesto en Matlab, ha de sufrir varios pasos de adaptación para ser ejecutable y práctico en C# y compatible con Unity 3D. Figura 3.9 Estructura del modelo en C#. En Unity 3D, lo que llamamos simulación es el producto de la renderización de una escena 3D y su interacción con los scripts que se hayan programado para la escena, todo ello realizado en cada
20 Capítulo 3. Modelo de la planta frame de la simulación. Basándonos en lo descrito en la sección 2.2.1 podemos inicializar los valores necesarios en la función de Start (en nuestro caso usamos Awake ) de Simulacion.cs, y ejecutar el bucle de simulación del modelo en la función Update (en nuestro caso, usamos FixedUpdate), que se ejecutará para cada frame de la simulación. 3.2.3 Conceptos clave del lenguaje C# Este lenguaje de programación tiene algunas particularidades en las que nos apoyamos para construir el modelo. Las nociones más relevantes en este contexto son: • Clases: las clases son los bloques de construcción básicos de la programación orientada a objetos en C#. Se utilizan para crear objetos que tienen datos y comportamientos. • Objetos: los objetos son instancias de clases. Tienen los datos y comportamientos definidos en su clase. Nos aprovechamos, además, del concepto de clase estática. Una clase estática en C# es una clase que no puede ser instanciada. Esto significa que no se puede crear una variable de tipo "clase estática". Las clases estáticas se utilizan para almacenar datos y funciones que son globales y no están asociados con ninguna instancia particular de una clase. Para declararlas se usa la palabra clave static. En la figura 3.9 hemos reflejado la estructura del modelo. En turquesa hemos señalado el script principal, Simulacion.cs; en un color más claro, FresnelSupport.cs y Fresnel.cs, que contienen clases estáticas con funciones que se emplean en los cálculos. En amarillo, las funciones auxiliares que hemos creado, junto con el script CSVReader.cs, que contiene otra función auxiliar para leer los datos de que disponemos. En rosa señalamos el fichero que almacena los datos experimentales, de los cuales nos servimos para probar el modelo. Salvo el script principal y el fichero de datos (que es simplemente un archivo .CSV), el resto de elementos del modelo son clases estáticas o funciones estáticas contenidas en ellas. 3.2.4 Simulacion.cs El programa principal de nuestro modelo se limita a declarar algunas variables y a llamar a las funciones contenidas en otros scripts. En la función Awake, primeramente modifica el tiempo de refresco, o tiempo entre frames, a 25ms, para que sea múltiplo del tiempo de muestreo del modelo 2 . void Awake() { Debug.Log("Simulacion ha empezado"); // Modifica el tiempo fijado de repetición del bucle. Para que tuviera sentido del todo,,→ // habría que ejecutar los cálculos cada 10 iteraciones del bucle (así, el tiempo de,→ // muestreo es realmente 0.25 s, como es en el modelo en matlab). De todas formas, los,→ // resultados son aceptables sin llevar a cabo esta corrección sugerida. Time.fixedDeltaTime = 0.025f; 2 El tiempo entre frames por defecto en Unity 3D es de 20ms, en tanto que el tiempo de muestreo de nuestro modelo es de 250ms
3.2 Modelo de parámetros distribuidos 21 ... Acto seguido invoca el Setup del script Fresnel.cs, y toma los valores de temperatura del fluido y del metal que se han generado en dicha función (luego veremos que se inicializan los vectores de valores a una misma cantidad). También se toman de Fresnel.cs los valores generados para el minuto y la hora de comienzo de la simulación (o sea, el momento en que comienzan las muestras sobre las que trabajaremos), el valor de eficiencia media calculado para el sistema según sus características físicas, y el número total de muestras de que disponemos, teniendo en cuenta las muestras originales y la interpolación hecha para aumentar la densidad de datos. // Se pone en marcha el sistema: Fresnel.Setup(); TemperaturaFluidoSimulacion = Fresnel.Tf; TemperaturaMetalSimulacion = Fresnel.Tm; // A hora y minuto de la simulación se le dan los valores del primer elemento de las muestras,,→ // porque a partir del primer paso, el mismo bucle de simulación las calculará.,→ minutoSimulacion = Fresnel.minuto1[0]; // Atencion al espíritu o la intención de la declaración en el modelo original: no es lo que estamos,→ // haciendo: (original:) x = find(abs(hora-horacomienzo)==min(abs(hora-horacomienzo))); hora = hora(x); ,→ ,→ horaSimulacion = Fresnel.hora1[0]; EficienciaMediaSimulacion = Fresnel.Eficienciamedia; j=0; Ntotal = (int)Fresnel.Ntotal; La función FixedUpdate actualiza las variables de trabajo correspondientes a la entrada a los valores de la iteración presente. FixedUpdate se repetirá hasta que ya no queden más muestras en el archivo de datos (entonces acabará la simulación); y cada vez que se ejecuta, corresponde a una muestra. Por ello, en cada iteración, asignamos: • La temperatura de entrada del fluido en la tubería correspondiente a la muestra actual, a la variable instantánea de temperatura del fluido 3. •La irradiancia incidente según la muestra actual, a la variable instantánea de irradiancia. • La temperatura ambiente recogida en la muestra actual, a la variable instantánea de temperatura ambiente. 3 asignamos ese valor a la primera posición del array que contiene la temperatura de cada parte discretizada del fluido en ese instante de la simulación
22 Capítulo 3. Modelo de la planta • El caudal correspondiente a la muestra actual, a la variable instantánea de caudal por la tubería. TemperaturaFluidoSimulacion[0] = Fresnel.Tentrada1[j]; // Inicializamos el primer valor (la temperatura de entrada del fluido),→ IrradianciaSimulacion = Fresnel.I1[j]; // Es constante para todo el sistema,→ TemperaturaAmbienteSimulacion = Fresnel.Tambiente1[j];// Es constante para todo el sistema,→ qSimulacion = Fresnel.caudal1[j]; // Es constante para toda la tubería,→ Luego calcula la variable de porcentaje del tubo en que incide radiación solar. Hemos simplificado ese valor a 1 o 0, según si los espejos están activos (en ese caso, consideramos que todos están activos y apuntan al tubo) o inactivos (sucede cuando la temperatura del fluido alcanzada en la salida supera los 180ºC, umbral a partir del cual, por seguridad, la planta debe cesar su captación de energía. También corrige el valor de caudal si, por cálculo, saliera negativo (en ese caso lo iguala a 0). PorcentajeTuboSimulacion = 1; // Si la temperatura de salida es mayor que 180ºC if(TemperaturaFluidoSimulacion[Fresnel.numeroPartesDiscretasTuboCaptador] > (double)180) PorcentajeTuboSimulacion = 0;,→ if(qSimulacion < 0) qSimulacion = 0; Cada vez que transcurre el tiempo marcado para recalcular el término de eficiencia (lo veremos más adelante al revisar las funciones dentro de Fresnel.cs en la sección 3.2.7), se actualiza la cuenta de tiempo y el factor de eficiencia según la función eficienciaGeoYSombras, que veremos en la sección 3.2.5. if(j % Fresnel.tactualizacion == 0) { minutoSimulacion = minutoSimulacion + (double)10 / (double)60; if(minutoSimulacion >= (double)60) { horaSimulacion++; if(horaSimulacion == 24) horaSimulacion = 0; minutoSimulacion -= (double)60; } TiempoSolarSimulacion = FresnelSupport.CalculoHoraSolar(horaSimulacion, minutoSimulacion, Fresnel.mes, Fresnel.angulodiario1);,→ FactorSombraSimulacion = FresnelSupport.eficienciaGeoYSombras(Fresnel.angulodiario1, TiempoSolarSimulacion); ,→ ,→ }
3.2 Modelo de parámetros distribuidos 23 Por fin, las operaciones más centrales del modelo son las relativas al cálculo de la temperatura del metal de la tubería y el fluido portador. Estos se llevan a cabo en la función CalculoTemperatura (ver sección 3.2.6). if(j % Fresnel.tactualizacion == 0) { minutoSimulacion = minutoSimulacion + (double)10 / (double)60; if(minutoSimulacion >= (double)60) { horaSimulacion++; if(horaSimulacion == 24) horaSimulacion = 0; minutoSimulacion -= (double)60; } TiempoSolarSimulacion = FresnelSupport.CalculoHoraSolar(horaSimulacion, minutoSimulacion, Fresnel.mes, Fresnel.angulodiario1);,→ FactorSombraSimulacion = FresnelSupport.eficienciaGeoYSombras(Fresnel.angulodiario1, TiempoSolarSimulacion); ,→ ,→ } Y con esto obtenemos ya los valores de temperatura de cada una de las partes del fluido y el metal de la tubería en el instante actual. Estos valores deben ser guardados para posterior estudio o correcciones del modelo: los almacena en un archivo .txt. Además, comprueba que queden datos en el archivo para finalizar o no la ejecución. // Store data for later plotting // In each iteration, append the value to the file File.AppendAllText(filePath, TemperaturaFluidoSimulacion[64].ToString() + Environment.NewLine);,→ j++; if(j > Ntotal) UnityEditor.EditorApplication.isPlaying = false; 3.2.5 Funciones del modelo óptico Obtener energía térmica del sol implica que la óptica va a jugar un papel importante en nuestro trabajo. Conseguir que el tubo reciba siempre el máximo de radiación solar posible (respetando unos límites de seguridad) tiene su complejidad: hay que conocer la posición del sol en todo momento, para reflejar los haces de luz sobre el tubo; hay que conocer el nivel de irradiancia en función de sombras de la propia estructura y de las nubes; basándose en eso y en la situación y orientación geográfica de la planta, hay que calcular el movimiento de los espejos. Todas las funciones del modelo óptico están alojadas en el script FresnelSupport.cs, dentro de su clase estática pública ("global"). Previamente al cálculo de los parámetros del modelo óptico que nos permitirán conocer la transferencia de energía al tubo y al fluido caloportador, hemos de calcular la hora solar, que nos brinda una información más veraz para el cálculo de la posición del sol que la hora oficial.
24 Capítulo 3. Modelo de la planta Figura 3.10 Esquema de las funciones del modelo óptico. CalculoDiaJuliano La función recibe como entrada la fecha de la simulación (año, mes y día) y devuelve, en double, el correspondiente día del año (contando entre 1 y 365). public static double CalculoDiaJuliano(double anyo, double mes, double dia) { int[] tabla = new int[] {31,28,31,30,31,30,31,31,30,31,30, 31};,→ int[] tablabisiesto = new int[] {31,29,31,30,31,30,31,31,30,31, 30,31};,→ int suma, i; double diajuliano; if( anyo % 4== 0&& ( anyo % 100 >0|| anyo % 400 == 0)) { tabla[1] = tablabisiesto[1]; } suma = 0; for(i=0;i<mes-1; i++) { suma = suma + tabla[i];
3.2 Modelo de parámetros distribuidos 25 } return diajuliano = suma + dia; } La única cuestión a destacar de este código, puesto que su funcionamiento es bastante sencillo, es el modo de determinar cuando un año es bisiesto 4 . Un posible método es comprobar la divisibilidad de la cifra del año por 4 ypor 100 ó400. Si cumple esa condición, se trata de un año bisiesto. CalculoAnguloDiario Ni siquiera merece la pena comentar una función tan sencilla. Está hallando la parte proporcional de circunferencia que se ha recorrido en los días del año pasados, si cada día se recorriese una parte alícuota. public static double CalculoAnguloDiario(double Diajuliano) { double anguloDiario; return anguloDiario = (2*Math.PI/365)*(Diajuliano - 1); } CalculoHoraSolar Se corrige la hora oficial tomando en consideración tres factores: la diferencia originada por la posición de la planta dentro de la zona horaria en que se halla (no está justo sobre un meridiano); la diferencia debida al horario de verano o de invierno, sumando la hora de desfase que tiene España (GMT+1 de base); y, por último, la ecuación de tiempo de Spencer. Con todo ello, la función recibe la hora y el minuto de las muestras de la simulación, el mes en que se realizaron las medidas, y el ángulo diario calculado; y devuelve la hora corregida en la variable Ts. public static double CalculoHoraSolar(double hora_HoraSolar, double minuto_HoraSolar, double mes_HoraSolar, double angulodiario_HoraSolar) { double Longitudmeridianocentralhuso, diferenciaporposicion, diferenciapormes, diferenciatotal;,→ double horaactualoficial, Et, Ts; horaactualoficial = hora_HoraSolar; // Hora oficial en decimal horaactualoficial = horaactualoficial + minuto_HoraSolar/60; Longitudmeridianocentralhuso = 0;// Meridiano de Greenwich // Diferencia por posición en minutos diferenciaporposicion = Math.Abs((Fresnel.Longitud - Longitudmeridianocentralhuso)*4);,→ 4Para más información acudir a https://learn.microsoft.com/en-us/office/troubleshoot/excel/determine-a-leap-year
26 Capítulo 3. Modelo de la planta // diferencia por horario de verano if(mes_HoraSolar < 4|| mes_HoraSolar > 10) diferenciapormes = 60;// 1 hora si no estamos entre abril y octubre,→ else diferenciapormes = 120;// 2 horas entre abril y octubre // Ecuación del tiempo Et = (0.000075+0.001868*Math.Cos(angulodiario_HoraSolar) -0.032077*Math.Sin(angulodiario_HoraSolar) -0.014615*Math.Cos(2*angulodiario_HoraSolar) -0.04089*Math.Sin(2*angulodiario_HoraSolar))*229.18;// En minutos // Diferencia total en minutos diferenciatotal = diferenciapormes + diferenciaporposicion - Et; return Ts = horaactualoficial - (diferenciatotal/60); } eficienciaGeoYSombras Esta función llama a otras funciones para obtener un factor unificado que sirva para valorar el efecto de las sombras, la geometría y otros aspectos del modelo óptico, sobre el resultado de la radiación incidente en el tubo. La función recibe el ángulo diario calculado y el tiempo solar corregido. public static double eficienciaGeoYSombras(double AnguloDiario_EficienciaGeoYSombras,,→ double Ts_EficienciaGeoYSombras) { El array XE contiene las coordenadas en dirección horizontal de cada fila de espejos. Definimos también la anchura de cada espejo, que es 25cm, y otros parámetros físicos de la planta. double[] XE = {-3.5, -2.8, -2.1, -1.4, -0.7,0,0.7,1.4,2.1,2.8,3.5}; double Anchuraespejos = 0.25; int Ytubo = 4, Areatotal = 352, longitudPlanta = 64; La primera función a la que llama es CalculoAzimut: su entrada es la hora solar y devuelve el ángulo azimutal del sol en ese instante y la altura solar. Con esa información, más las coordenadas horizontales de los espejos y la altura del tubo captador, podemos calcular la posición que deben tener los espejos para reflejar correctamente la luz sobre el tubo. La función a la que llama a continuación se ocupa de hacer ese cálculo en dos dimensiones. Obtenemos el vector que une el centro del espejo con el centro del tubo captador. resultadoAzimut = CalculoAzimut(Ts_EficienciaGeoYSombras); alturasolar = resultadoAzimut[0];
3.2 Modelo de parámetros distribuidos 27 azimut = resultadoAzimut[1]; resultadoVectorPosicionEspejos = CalculoVectorPosicion2DEspejos(alturasolar, azimut, XE, Ytubo);,→ Teniendo esto, calculamos los ángulos que han de formar los espejos: su inclinación respecto del sistema de coordenadas (función CalculoInclinacionEspejos). Esto nos permite también prever qué sombras se arrojan unos sobre otros, que puede ser resumido en un factor (función CalculoFactorSombra). Y dado eso, sacamos un coeficiente de eficiencia geométrica. resultadoCalculoInclinacionEspejos = CalculoInclinacionEspejos(XE, i2d, r2d, modi2d, gamma);,→ alpha = resultadoCalculoInclinacionEspejos[0]; theta = resultadoCalculoInclinacionEspejos[1]; beta = resultadoCalculoInclinacionEspejos[2]; kappa = resultadoCalculoInclinacionEspejos[3]; inclinacion = resultadoCalculoInclinacionEspejos[4]; factorsombra = resultadoCalculoInclinacionEspejos[5]; factorsombra = CalculoFactorSombra(XE, inclinacion, Anchuraespejos, gamma, factorsombra);,→ resultadoCalculoEficienciaGeometrica = CalculoEficienciaGeometrica(factorsombra, Areatotal, inclinacion, kappa, gamma); ,→ ,→ Aglutinando todos estos factores, sacamos un coeficiente que representa la eficiencia del captador contando con su geometría y las sombras que proyectan unos espejos sobre otros (función CalculoParteSombreada). Factorsombra = CalculoParteSombreada(inclinacion, alturasolar, i3d, longitudPlanta, Ytubo, eficienciageometrica, Pps, azimut, XE);,→ return Factorsombra; } CalculoAzimut El azimut o ángulo azimutal es la magnitud angular en grados, generalmente medida en sentido horario desde el norte, que describe la dirección de un objeto o punto en relación con el norte en un plano horizontal. Es una medida comúnmente utilizada en astronomía, geología, navegación y otras disciplinas para indicar la dirección o posición de un objeto o punto en el horizonte. El azimut solar se ilustra en la figura 3.11. Para el cálculo del azimut y la altura solar, sólo necesitamos saber la hora solar, que hemos calculado previamente. private static double[] CalculoAzimut(double Ts)
34 Capítulo 3. Modelo de la planta resultado[0, l] = pf[l]; resultado[1, l] = Cf[l]; resultado[2, l] = Ht[l]; } CalculoPerdidasMetal Esta función se corresponde a una redistribución de los parámetros de la ecuación 3.1, ˙ Qa=Dmπha(Tm(k,n)−Ta(k))3+b(Tm(k,n)−Ta(k))i(3.1) Dm es el diámetro equivalente de las paredes del tubo, a y b son los coeficientes de pérdidas de calor de la función polinómica del absorbedor, que quedan reunidos en el resultado de esta función, Hl. for(m=0; m < tamanyoT_PerdidasMetal; m++) { Hl[m] = (0.45659247191149893/Sup)*(T[m] - Tambiente) - 1.062045206593238/Sup;,→ } for(m=64;m<=164; m++) { Hl[m] = 0.0140251565842975*(T[m] - Tambiente)*0.09 + 7.328725247698581*0.09;,→ } for(m=165; m<tamanyoT_PerdidasMetal; m++) { Hl[m] = 0*(T[m] - Tambiente) + 0.290223*5; } return Hl; Donde el valor de Sup es el área de la superficie total de espejos, en metros. Los metros 0-63 de la tubería se corresponden al tubo receptor; los metros 64 en adelante, al resto de la tubería y el intercambiador de calor. De ahí los índices de los bucles for, que aplican fórmulas distintas según la zona del sistema para la que se esté haciendo el cálculo. CalculoTemperatura Esta función responde a las ecuaciones 3.2 y 3.3: Tm(k+1,n) = Tm(k,n)+∆t ˙ Qsol(k) ρmcmAm − ˙ Qa(k,n) ρmcmAm − ˙ Qf(k,n) ρmcmAm!(3.2) Tf(k+1,n) = Tf(k,n)+∆t q(k) Af Tf(k,n)−Tf(k,n−1) ∆x+ ˙ Qf(k,n) ρfcfAf!(3.3) En primer lugar invocamos las funciones previamente explicadas, de cálculo de pérdidas en el metal y características del fluido portador.
3.2 Modelo de parámetros distribuidos 35 El bucle for se rige por el número de columnas; y, como los datos calculados los recibo de la función de cálculo de las características por filas, tengo que recorrer cada fila entera recogiendo todos los datos del mismo tipo. Hl_local = CalculoPerdidasMetal(TactualPerdidas, Tambiente_CalcTemp); aux2 = CalculoCaracteristicasDelFluido(TactualFluido_CalcTemp, Caudal_CalcTemp);,→ pf_local = new double[aux2.GetLength(1)]; Cf_local = new double[aux2.GetLength(1)]; Ht_local = new double[aux2.GetLength(1)]; for(n=0; n < aux2.GetLength(1); n++) { pf_local[n] = aux2[0, n]; Cf_local[n] = aux2[1, n]; Ht_local[n] = aux2[2, n]; } Ahora calculamos el vector de radiación, usando el factor de sombra que ha resultado de todas las operaciones del modelo óptico. Incorporamos también al cálculo la eficiencia del sistema y la variable Porcentajetubo que indicaba si los espejos apuntan al tubo receptor o no (por seguridad en casos de sobrecalentamiento). Salta a la vista que sólo lo calculamos para las 135 primeras partes del modelo: se debe a que los tubos del interior del intercambiador no reciben radiación solar; y los trozos correspondientes a la tubería después del receptor están aislados y consideramos despreciable la radiación incidente. vectorradiacion = new double[135]; for(n=0;n<64; n++) { vectorradiacion[n] = I_CalcTemp*G*Eficienciamedia_CalcTemp*Factorsombra_CalcTemp*Porcentajetubo;,→ } for(n=64;n<135; n++) { vectorradiacion[n] = 0; } Calculamos la temperatura del metal del tubo receptor y actualizamos la variable de salida. for(n=0;n<64; n++) { Tmetal_CalcTemp[n] = TactualMetal_CalcTemp[n] + (tint_CalcTemp/(pm*Cm*Am)),→ *(vectorradiacion[n] - Hl_local[n]*G*(TactualMetal_CalcTemp[n] - Tambiente_CalcTemp),→ - L*Ht_local[n]*(TactualMetal_CalcTemp[n] - TactualFluido_CalcTemp[n]));,→ TemperaturasCalcTemp[0, n] = Tmetal_CalcTemp[n]; }
36 Capítulo 3. Modelo de la planta Calculamos la temperatura del fluido cuando el caudal es nulo, y actualizamos el primer elemento del array de salida correspondiente. for(n=0;n<64; n++) { TfluidoSinCorregir[n] = TactualFluido_CalcTemp[n] + ((L*Ht_local[n]*tint_CalcTemp)/(Af*pf_local[n]*Cf_local[n])),→ *(TactualMetal_CalcTemp[n] - TactualFluido_CalcTemp[n]);,→ } TemperaturasCalcTemp[1,0] = TfluidoSinCorregir[0]; Calculamos, para el resto de elementos del tubo captador en el array de salida, la temperatura del fluido, teniendo en cuenta la energía por transporte de caudal [13] [8]. for(n=1;n<64; n++) { Tfluido_CalcTemp[n] = TfluidoSinCorregir[n] - (q*tint_CalcTemp)/(Af*incx[n-1]),→ *(TfluidoSinCorregir[n] - TfluidoSinCorregir[n-1]);,→ // Debug.Log("Tfluido captador: "+Tfluido_CalcTemp[n]); TemperaturasCalcTemp[1, n] = Tfluido_CalcTemp[n]; } Calculamos ahora la temperatura del metal para el resto de la tubería. Nótese que las ecuaciones incorporan los términos correspondientes. for(n=64;n<135; n++) { Tmetal_CalcTemp[n] = TactualMetal_CalcTemp[n] + (tint_CalcTemp/(pm*Cm*Am2)),→ *(vectorradiacion[n] - Hl_local[n]*Math.PI*0.16*(TactualMetal_CalcTemp[n] - Tambiente_CalcTemp)*0.1 ,→ ,→ - L2*0.01*Ht_local[n]*(TactualMetal_CalcTemp[n] - TactualFluido_CalcTemp[n]));,→ TemperaturasCalcTemp[0, n] = Tmetal_CalcTemp[n]; } for(n=135;n<254; n++) { Tmetal_CalcTemp[n] = TactualMetal_CalcTemp[n]; TemperaturasCalcTemp[0, n] = Tmetal_CalcTemp[n]; } A continuación, se obtiene la temperatura del fluido como antes, pero para la parte del sistema distinta del receptor.
3.2 Modelo de parámetros distribuidos 37 for(n=64;n<254; n++) { TfluidoSinCorregir[n] = TactualFluido_CalcTemp[n] - L2*Hl_local[n]*(TactualFluido_CalcTemp[n] - Tambiente_CalcTemp),→ *tint_CalcTemp/(Af2*pf_local[n]*Cf_local[n]); } Y de modo análogo, se corrige el cálculo en estado estacionario, teniendo en cuenta la energía por transporte de caudal. for(n=64;n<254; n++) { Tfluido_CalcTemp[n] = TfluidoSinCorregir[n] - (q*tint_CalcTemp)/(Af2*incx[n-1]),→ *(TfluidoSinCorregir[n] - TfluidoSinCorregir[n-1]);,→ TemperaturasCalcTemp[1, n] = Tfluido_CalcTemp[n]; } 3.2.7 Funciones auxiliares A parte de los cálculos propios del modelo óptico y térmico de la planta, necesitamos de otras funciones que solucionen algunas necesidades. CrossProduct3D Por ejemplo, para el cómputo del producto vectorial del vector sol-espejo con el vector inclinación del espejo, necesitamos una función que lo calcule, puesto que C# implementa un producto vectorial entre vectores (arrays de tipo vector), pero nosotros no trabajamos con vectores (por comodidad). private static double[] CrossProduct3D(double[] a, double[] b) { double[] result = new double[3]; result[0] = a[1] * b[2] - a[2] * b[1]; result[1] = a[2] * b[0] - a[0] * b[2]; result[2] = a[0] * b[1] - a[1] * b[0]; return result; } Fresnel.Setup La función de Setup dentro del script Fresnel.cs deja todo listo para poder comenzar la simulación. Invoca Inicializacion yObtenerDatos, y llamando a Interpola, aumenta la densidad de muestras. También inicializa las variables del día juliano, el ángulo diario, el tiempo solar, la temperatura del fluido y la del metal (no mostramos aquí el código). Inicializacion(); ObtenerDatos();
38 Capítulo 3. Modelo de la planta numeroperiodosint = (int)Math.Ceiling(incrementodetiempo / tint); hora1 = Interpola(numeroperiodosint, hora1); caudal1 = Interpola(numeroperiodosint, caudal1); Tentrada1 = Interpola(numeroperiodosint, Tentrada1); Tsalida1 = Interpola(numeroperiodosint, Tsalida1); I1 = Interpola(numeroperiodosint, I1); Tambiente1 = Interpola(numeroperiodosint, Tambiente1); ... Hallamos otros tres parámetros: • tiempomax: el número máximo de iteraciones que hará el bucle principal. Se corresponde con la longitud del archivo de datos (podríamos haber puesto una variable cualquiera de las medidas en las muestras). • Ntotal: el número de iteraciones del bucle principal de simulación. Lo calculamos en base a los tiempos que se indican en la función Inicializacion. En caso de que sea mayor que tiempomax, se sobreescribe con su valor (de otro modo, llegaría un punto en que no habría datos para continuar con la simulación). • tactualizacion: indica la frecuencia con que se actualiza FactorSombra, que aglutina, como se ha explicado antes, la eficiencia geométrica y los efectos de las sombras en el captador. tiempomax = Tentrada1.Length; Ntotal = Math.Floor(Fresnel.tiemposimulacion / Fresnel.tint); if(Ntotal > tiempomax) Ntotal = tiempomax; tactualizacion = Math.Ceiling(Fresnel.tacteficiencias / Fresnel.tint);
3.2 Modelo de parámetros distribuidos 39 Fresnel.Inicializacion Define multitud de variables necesarias para la operación. Primeramente, la fecha correspondiente al set de datos utilizado. Después, el tiempo de integración (tint), el tiempo de muestreo (incrementodetiempo), el tiempo que transcurre desde el primer hasta el último dato de la muestra (tiemposimulacion), y el tiempo que esperaremos entre actualizaciones de las eficiencias (tacteficiencias). Después define algunos parámetros físicos de la planta. Tabla 3.5 Parámetros definidos en la función de inicialización. Parámetro (símbolo) Valor Unidad Tiempo de integración 0.25 s Tiempo de actualización de las eficiencias 10 s Diámetro interior de la tubería que continúa el circuito 0.07699 m Diámetro interior de la tubería que continúa el circuito 0.16 m Absortancia 0.94 - Transmisividad 0.96 - Reflectividad 0.8 - Reflectividad secundaria 0.77 - Factor de ensuciamiento 1a0.8 - Factor de ensuciamiento 2b0.8 - aAsumiendo que el espejo está sucio bAsumiendo que el espejo está sucio Por último, realiza una serie de cálculos: Halla las eficiencia media del espejo en su reflexión, y nuevos parámetros físicos que usaremos en los cálculos. Eficienciamedia = Absortancia*Reflectividad*Transmisividad*Reflectividadsecundaria,→ *factorensuciamiento1*factorensuciamiento2; // Para Tubería 1 // Diametro propio de la tuberia [m] L = Math.PI*diametrointeriortubo; // Áreas transversales del fluido[m2] Af = Math.PI*Math.Pow(diametrointeriortubo, 2)/4; Am = Math.PI*(Math.Pow(diametroexteriortubo, 2) - Math.Pow(diametrointeriortubo, 2))/4;,→ // Para Tubería 2 // Diametro propio de la tuberia 2[m] L2 = Math.PI*diametrointeriortubo2; // Áreas transversales del fluido [m2] Af2 = Math.PI*Math.Pow(diametrointeriortubo2, 2)/4; Am2 = Math.PI*(Math.Pow(diametroexteriortubo2, 2) - Math.Pow(diametrointeriortubo2, 2))/4;,→
40 Capítulo 3. Modelo de la planta CSVReader.cs Lo primero que debemos hacer es hallar el número de líneas que contiene el archivo de datos (es decir, el número de muestras tomadas). Queremos averiguar también el número de columnas, pues así es como están separados nuestros datos (cada variable forma una columna): para ello, identificamos una línea cualquiera (hemos tomado la primera, por comodidad) y contamos el número de valores separados por comas. public static class CSVReader { public static string[,] LeerCSV(string filePath) { int i=0, j; string[] values; string[,] valor; int numberOfLines = File.ReadLines(filePath).Count(); var firstLine = File.ReadLines(filePath).FirstOrDefault(); var numberOfColumns = firstLine.Split(',').Length; valor = new string[numberOfLines, numberOfColumns]; ... Comprobamos que exista el archivo y usamos la clase StreamReader 5 de la librería System.IO de C# para recorrer el archivo. En el bucle while leemos la línea en curso y separamos las cantidades almacenadas siguiendo las comas que delimitan el archivo. Obtenemos tantos string como comas haya en el archivo de datos, de modo que podemos tomar cada string y guardarlo en el array creado a tal efecto. ... if (File.Exists(filePath)) { using (StreamReader reader = new StreamReader(filePath)) { while (!reader.EndOfStream) { string line = reader.ReadLine(); values = line.Split(','); j=0; foreach (string value in values) { valor[i, j] = value; j++; } i++; } } 5Consultar [7]
3.2 Modelo de parámetros distribuidos 41 } else { Debug.LogError("File not found: " + filePath); } return valor; } } Fresnel.ObtenerDatos Cargamos la dirección del archivo que almacena los datos; declaramos las variables; en el bucle for, convertimos los datos leídos de caracteres (formato en el que salen de la función de lectura) a double y se almacena en la variable pertinente. Mostramos el caso de una variable, pero se repite el procedimiento para todas (ver Anexo). private static void ObtenerDatos() { string[,] DatosCSV = CSVReader.LeerCSV(@"C:\Users\Pablo\OneDrive - UNIVERSIDAD DE SEVILLA\TFG\0_Integracion\Assets\DataFiles\datosMedidos.csv"); ,→ ,→ bool success; int longitudBucle = DatosCSV.GetLength(0); hora1 = new double[longitudBucle]; ... for(int contador = 0; contador < longitudBucle; contador++) { // Convertimos los datos leídos de caracteres a [double] el formato que tenga la variable de destino.,→ success = Double.TryParse(DatosCSV[contador, 0], out hora1[contador]); // h,→ if(!success) Debug.Log("Error al convertir hora a formato double, en el indice: " + contador);,→ ... } } Fresnel.Interpola Se implementa una interpolación lineal con el cálculo previo del primer valor. private static double[] Interpola(int numeroMuestras, double[] datosOrigen) { int tam = datosOrigen.Length; int tamvect = numeroMuestras; int n = tamvect;
42 Capítulo 3. Modelo de la planta ... for(int i = 1; i <= tamvect; i++) { vectorsalida[i - 1] = datosOrigen[0] + ((datosOrigen[1] - datosOrigen[0]) / numeroMuestras) * i;,→ } for(int i = 1; i < tam; i++) { for(int j = 1; j <= tamvect; j++) { if(i >= tam - 1) { vectorsalida[n] = datosOrigen[i - 1] + ((datosOrigen[i] - datosOrigen[i - 1]) / numeroMuestras) * j;,→ } else vectorsalida[n] = datosOrigen[i] + ((datosOrigen[i + 1] - datosOrigen[i]) / numeroMuestras) * j;,→ n++; } } ... 3.2.8 Movimiento del sol Para comprobar el correcto funcionamiento del modelo 3D que explicamos en la sección 3.3, y a falta de la conexión entre las variables del modelo físico y los elementos del modelo 3D, hemos creado unos script adicionales para simular el movimiento del sol y demostrar que los espejos cumplen su función adecuadamente. Se trata de un cálculo simplificado respecto del modelo óptico de la planta Fresnel que se ha explicado anteriormente. La primera simplificación es que tan sólo considera movimiento del sol en un plano perpendicular al eje longitudinal de la planta. Esto facilita bastante los cálculos, convirtiéndolos en bidimensionales. En el Awake del script principal (SolarMovement.cs) se identifica la orientación de la planta: void Awake() { fresnelSystem = SolarMovementSupport.FindFresnelObjects(); rotacionYplanta = (float)GameObject.Find("Planta").transform.localRotation.eulerAngles.y;,→ rotacionZplanta = (float)GameObject.Find("Planta").transform.localRotation.eulerAngles.z;,→ } En el FixedUpdate: •Se rota el sol un ángulo fijo cada frame para simular su movimiento. •Se detecta su ángulo de rotación.
3.2 Modelo de parámetros distribuidos 43 • Se calcula el ángulo de rotación correspondiente de cada espejo (en esta versión del código hacemos el cálculo para cada espejo, aunque se podría reducir a un cálculo por filas en este modelo simplificado). •Se comprueba que el sol no haya completado media circunferencia. void FixedUpdate() { fresnelSystem.Sun.Rotate(SolarMovementSupport.variacionMovimientoSolar, 0,0);,→ rotacionActualSol = (float)GameObject.Find("Sun").transform.localRotation.eulerAngles.x;,→ for(int i = 0; i < fresnelSystem.Mirrors.GetLength(0); i++) { anguloRotacionEspejo = (float)SolarMovementSupport.MirrorAngle(i + 1, anguloRotacionAnterior);,→ rotacionActualEspejo = (float)fresnelSystem.Mirrors[i, 0].transform.localRotation.eulerAngles.x;,→ // Recorremos todas las filas y todas las columnas para girar los espejos,→ for(int j = 0; j < fresnelSystem.Mirrors.GetLength(1); j++) fresnelSystem.Mirrors[i, j].transform.rotation = Quaternion.Euler(anguloRotacionEspejo, rotacionYplanta, rotacionZplanta); ,→ ,→ anguloRotacionAnterior = anguloRotacionEspejo; } if(GameObject.Find("Sun").transform.localRotation.eulerAngles.x > 180) UnityEditor.EditorApplication.isPlaying = false;,→ } En el script de soporte hemos incluido la función que calcula la rotación que deben tener los espejos. Para ello, primero encontramos cada uno de los espejos: public static FresnelSystem FindFresnelObjects() { GameObject[] aux; int i = 0,j=0; FresnelSystem fresnelSystem; fresnelSystem.Sun = GameObject.Find("Sun").transform; aux = GameObject.FindGameObjectsWithTag("Mirror"); fresnelSystem.Mirrors = new Transform[aux.Length/numberOfMirrorsInLine, numberOfMirrorsInLine];,→ foreach(GameObject mirror in aux)
50 Capítulo 3. Modelo de la planta Figura 3.22 Barras laterales. Figura 3.23 Cubierta del tubo receptor. 3.4 Conexión de los modelos El verdadero potencial de este modelo se despliega cuando se combinan el modelo físico con la visualización 3D. En este trabajo no se ha llegado a explotar completamente, puesto que las funciones de cálculo no toman los datos de Unity, sino de los archivos de datos recogidos. Sin embargo, ofrecemos una muestra con las funciones que modelan el comportamiento del sol y los espejos. 3.4.1 Importar el modelo 3D en Unity Para importar el modelo a Unity 3D del modo más eficiente posible, lo hemos diseñado entero y "ensamblado" en Blender, de manera que en Unity queda como un solo prefab. Esto permite, por otro lado, mejor compatibilidad para las herramientas de control de versiones (hemos estado trabajando en un repositorio de GitHub) y que los otros usuarios puedan visualizar el modelo en su
3.4 Conexión de los modelos 51 ordenador. El problema de esto es que, de este modo, no es posible acceder a los componentes del modelo por separado. Para solucionarlo hay que desempaquetar el prefab. Otra de las cuestiones importantes es la relación padre-hijo: el objeto Planta es el padre de todo el resto de objetos que conforman la planta, por ejemplo. Figura 3.24 Esquema de relaciones de los objetos del modelo. Por último, nos hemos aprovechado de la instanciación de prefabs para la propagación de cambios entre objetos. Por ejemplo, al convertir la superficie de los espejos a plano, bastó con cambiarla en uno de los espejos de cada fila, puesto que eran instancias del mismo prefab. 3.4.2 Incluir los scripts Como decíamos en la sección sobre el modelado de los espejos, es esencial preservar el orden y la coherencia en los nombres y la estructuración de los objetos del modelo. Esto se acentúa más cuando hay scripts que actúan sobre objetos de la escena. Para esa interacción, el script que hemos usado se fija en las etiquetas que asignamos a los objetos del mismo tipo. Figura 3.25 Cada espejo del está etiquetado para ser identificado por el script.
52 Capítulo 3. Modelo de la planta El único script del modelo físico que está vinculado a un GameObject es Simulacion.cs, ya que desde ahí se llama a todas las demás funciones. Esto facilita la subsanación de errores. Figura 3.26 Vista del inspector que muestra los scripts de Simulacion.cs y SolarMovement.cs asociados al objeto Planta.
4 Conclusiones y trabajos futuros Hemos sido capaces de implementar un modelo de parámetros distribuidos de una planta solar de tipo Fresnel usando C# y Unity 3D, basándonos en otro modelo diseñado en Matlab; y hemos desarrollado también un modelo en 3D del captador solar, funcional y listo para su integración con el modelo computacional. Para ello, previamente, hemos llevado a cabo un trabajo de investigación de varias etapas en torno al concepto de Gemelo Digital. Primero, procurando aprender a diseñar un SCADA de una planta física. Después, probando varias aproximaciones a la simulación en 3D de fluidos empleando Unity 3D. Por último, concentrándonos en el modelo de parámetros distribuidos que hemos desarrollado. A lo largo de este trabajo, hemos ido detectando aspectos que nos gustaría mejorar o profundizar a posteriori, pero, por limitación de tiempo o desconocimiento, no hemos sido capaces de acometer ahora. Entre ellas, están, además de las mencionadas en secciones anteriores: •La posibilidad de simular nubes, jugando con el nivel de irradiancia. • En vez de ver valores en Inspector, generar visualizadores para build y game mode, de manera que un usuario pueda ver algunos valores, y manipular otros. • Considerar la irradiancia individualmente en cada espejo para dar lugar a un modelo más preciso. • En el código, parametrizar el número de variables medidas, de manera que identifique el número y asigne cada columna a una variable poniéndole el nombre de la primera fila de la columna (habría que tratar distinto la primera fila y avisar de que los archivos .csv deberían llevar los nombres de cada variable en la primera fila de cada columna). • El archivo de datos debería poder direccionarse, o insertarse, al ejecutar la simulación (y no incluirlo en el código). • Aprovechando que tenemos un modelo en Unity, hacer un build del modelo, de manera que se pueda correr la simulación desde distintos sitios, disponiendo de los ficheros de datos. • Que el caudal que es una de las entradas al modelo, se tome de Zibra Liquids (usando el paquete). 53
Apéndice A Código comentado A.1 Fresnel.cs 1using UnityEngine; 2using System; 3using System.Linq; 4using System.IO; 5using System.Collections; 6using System.Collections.Generic; 7 8// Código de inicialización de parámetros y variables y asignación de valores.,→ 9 10 public static class Fresnel 11 { 12 // Aquí se incluye todo lo previo al bucle de simulacion. Es el Setup del proyecto.,→ 13 14 ///////////////////////////////////// 15 // Variables para el Setup 16 ///////////////////////////////////// 17 public static int numeroperiodosint; 18 public static double diajuliano1, angulodiario1, Ts1; 19 public static double tiempomax, Ntotal, tactualizacion; 20 public static double[] caudal1, Tentrada1, Tsalida1, I1, Tambiente1, minuto1, hora1;,→ 21 public static double[] Tf, Tm; 22 23 ///////////////////////////////////// 24 // Variables para Inicializacion 25 ///////////////////////////////////// 26 public static double Latitud, oriplanta, Longitud, Eficienciamedia;,→ 27 public static int incrementodetiempo, tacteficiencias; 28 public static double tint, diametrointeriortubo, diametroexteriortubo,,→ 55
56 Capítulo A. Código comentado 29 diametrointeriortubo2, diametroexteriortubo2; 30 public static float Absortancia, Transmisividad, Reflectividad, Reflectividadsecundaria,,→ 31 factorensuciamiento1, factorensuciamiento2; 32 public static long tiemposimulacion; 33 public static double Af, Am, G, L, pm, Cm, Am2, Af2, L2; 34 public static int anyo, mes, dia, i; 35 public static int numeroPartesDiscretasSistema = 254; 36 public static int numeroPartesDiscretasTuboCaptador = 64; 37 38 public static void Setup() 39 { 40 // Primero inicializamos los datos del sistema del captador. 41 Inicializacion(); 42 43 // Ahora preparamos el sistema para operación. 44 45 46 //--1--// Cargamos los datos experimentales 47 48 ObtenerDatos(); 49 // Aprovechamos que son variables públicas, no tenemos ni que recogerlas.,→ 50 51 52 //--2--// Interpolamos los datos de las medidas, para tener una base de muestras mayor,→ 53 54 numeroperiodosint = (int)Math.Ceiling(incrementodetiempo / tint);,→ 55 hora1 = Interpola(numeroperiodosint, hora1); 56 caudal1 = Interpola(numeroperiodosint, caudal1); 57 Tentrada1 = Interpola(numeroperiodosint, Tentrada1); 58 Tsalida1 = Interpola(numeroperiodosint, Tsalida1); 59 I1 = Interpola(numeroperiodosint, I1); 60 Tambiente1 = Interpola(numeroperiodosint, Tambiente1); 61 62 63 //--3--// Inicialización de las variables. 64 65 diajuliano1 = FresnelSupport.CalculoDiaJuliano(anyo, mes, dia);,→ 66 67 angulodiario1 = FresnelSupport.CalculoAnguloDiario(diajuliano1);,→ 68 69 minuto1 = new double[hora1.Length]; 70 for(i=0; i < hora1.Length; i++) 71 {
A.1 Fresnel.cs 57 72 minuto1[i] = Math.Floor(Math.Abs(hora1[i] - Math.Truncate(hora1[i]))*60);,→ 73 } 74 75 // Aplicamos, para inicializar el parámetro Ts1, la función, al primer elemento de los datos experimentales,→ 76 Ts1 = FresnelSupport.CalculoHoraSolar(hora1[0], minuto1[0], mes, angulodiario1);,→ 77 78 Tf=new double[numeroPartesDiscretasSistema]; // Tiene la longitud de las partes de fluido que consideramos en el sistema.,→ 79 for(i=0; i < numeroPartesDiscretasSistema; i++) 80 { 81 Tf[i] = Tentrada1[0]; // Inicializar la temperatura del fluido,→ 82 } 83 84 Tm=new double[Tf.Length]; 85 for(i=0; i < numeroPartesDiscretasSistema; i++) 86 { 87 Tm[i] = Tf[i] + 7;// Inicializar la temperatura del metal,→ 88 } 89 90 91 //--4--// Cálculos necesarios para el funcionamiento del programa.,→ 92 93 // Cálculo del número de iteraciones, correspondientes a los tiempos de,→ 94 // simulación e integración que se han definido en el fichero de,→ 95 // configuración. También se calcula el tiempo máximo de simulación.,→ 96 //-> tiempomax es el número máximo de iteraciones que hará el bucle principal.,→ 97 // Es así porque se corresponde con la longitud del archivo de datos (podríamos haber,→ 98 // puesto cualquier otra de las variables que contiene, aquí usamos "Tentrada1").,→ 99 tiempomax = Tentrada1.Length; 100 101 //-> Ntotal es el número de iteraciones del bucle principal de simulación,,→ 102 // calculado en base a los tiempos que hemos indicado en Inicializacion.,→ 103 // Si es mayor que tiempomax, se sobreescribe, pues no tendríamos datos para hacer,→ 104 // la simulación. 105 Ntotal = Math.Floor(Fresnel.tiemposimulacion / Fresnel.tint);
58 Capítulo A. Código comentado 106 if(Ntotal > tiempomax) Ntotal = tiempomax; 107 108 // Cada cuántas muestras se actualizan las eficiencias 109 //-> tactualizacion marca la frecuencia de actualización de FactorSombra, que aglutina bastantes cálculos,→ 110 // sobre la eficiencia geométrica y las sombras posibles en el captador.,→ 111 tactualizacion = Math.Ceiling(Fresnel.tacteficiencias / Fresnel.tint);,→ 112 113 Debug.Log("Fresnel se ha cargado"); 114 } 115 116 private static void Inicializacion() 117 { 118 Debug.Log("Los datos del sistema se han cargado"); 119 120 // Introducir fecha correspondiente a los datos que se usarán para la simulación,→ 121 // En realidad estos datos tenemos que cargarlos de un fichero, estas declaraciones no me sirven,→ 122 123 // anyo = 2009; 124 // mes = 06; 125 // dia = 16; 126 127 anyo = 2009; 128 mes = 06; 129 dia = 29; 130 131 // anyo = 2009; 132 // mes = 07; 133 // dia = 06; 134 135 // anyo = 2009; 136 // mes = 07; 137 // dia = 10; 138 139 // anyo = 2009; 140 // mes = 08; 141 // dia = 17; 142 143 // anyo = 2009; 144 // mes = 08; 145 // dia = 22; 146 147 148 149 // ------------------->>PARTE CONFIGURABLE DEL FICHERO>>--------------------,→
A.1 Fresnel.cs 59 150 //__________________________________________________________________________,→ 151 // Tiempo de integración, Tiempo de muestreo de sensores, Simulación,,→ 152 // ---------------------Representación en segundos--------------------------,→ 153 // Tiempo de integración 154 tint = 0.25; 155 156 // Tiempo de muestreo de sensores. 157 // Si se dan cada 15 minutos son 900 segundos. 158 incrementodetiempo = 60; 159 160 // Tiempo de simulación en segundos 161 tiemposimulacion = 7*3600; 162 163 // El tiempo para actualizar las eficiencias en segundos. 164 tacteficiencias = 10; 165 166 167 // Asignacion de valor de variables dimensionales 168 169 // Tuberia 1 [m] 170 diametrointeriortubo = 0.068;//[m] 171 diametroexteriortubo = 0.071;//[m] 172 173 // Tuberia 2 [m] 174 diametrointeriortubo2 = 0.076986659489532;//[m] 175 diametroexteriortubo2 = 0.16;//[m] 176 177 // Apertura del colector [m] 178 G=0.5*11; 179 180 // Características del metal DIN 1.4404 181 pm=8027;//Densidad del metal (DIN 1.4404) [kg/m3] 182 Cm=500;//Calor específico del metal [J/kg°C] 183 184 // MODELO ÓPTICO DEL CAPTADOR SOLAR FRESNEL 185 //__________________________________________________________________________,→ 186 // Variables de orientación y posición:Latitud y longitud 187 // y orientacion e inclinación del campo 188 189 Longitud = -6; 190 Latitud = 37.41; 191 oriplanta = (12 +3/600 +11/3600)*Math.PI/180; 192 // ¿Qué sentido tiene esta transformación? 193 // Entiendo que se debe a que la planta tiene desviacion de,→
66 Capítulo A. Código comentado 135 resultadoCalculoEficienciaGeometrica = CalculoEficienciaGeometrica(factorsombra, Areatotal, inclinacion, kappa, gamma); ,→ ,→ 136 Pps = resultadoCalculoEficienciaGeometrica[0][0];,→ 137 eficienciageometrica = resultadoCalculoEficienciaGeometrica[1];,→ 138 139 Factorsombra = CalculoParteSombreada(inclinacion, alturasolar, i3d, longitudPlanta, Ytubo, eficienciageometrica, Pps, azimut, XE); ,→ ,→ 140 141 return Factorsombra; 142 } 143 144 private static double[] CalculoAzimut(double Ts) 145 { 146 // Compmutes the azimuth and "solar height" for given data "solar time",→ 147 double[] resultado; 148 double angulodiario = Fresnel.angulodiario1; // This is to be updated on every iteration (in fact, only when angulodiario changes) ,→ ,→ 149 double alturasolar, azimut, sinazimut, declinacion, anghorario, latitud;,→ 150 151 latitud = Fresnel.Latitud; 152 153 // declinación calculada por Spencer 154 declinacion = 0.006918 -0.399912*Math.Cos(angulodiario) 155 +0.070257*Math.Sin(angulodiario) 156 -0.006758*Math.Cos(2*angulodiario) 157 +0.000907*Math.Sin(2*angulodiario) 158 -0.002697*Math.Cos(3*angulodiario) 159 +0.00148*Math.Sin(3*angulodiario); 160 161 // ángulo horario 162 anghorario = ((Ts - 12)*15)*(Math.PI/180); // en radianes 163 164 // cálculo de la altura y azimut solar 165 alturasolar = Math.Asin(Math.Sin(declinacion)*Math.Sin((latitud*Math.PI)/180),→ 166 + Math.Cos(declinacion)*Math.Cos((latitud*Math.PI)/180)*Math.Cos(anghorario));,→ 167 168 azimut = Math.Acos((Math.Sin(alturasolar)*Math.Sin((latitud*Math.PI)/180),→ 169 - Math.Sin(declinacion))/(Math.Cos(alturasolar)*Math.Cos((latitud*Math.PI)/180)));,→ 170
A.2 FresnelSupport.cs 67 171 // Cálculo del signo del azimut solar 172 sinazimut = (Math.Cos(declinacion)*Math.Sin(anghorario))/(Math.Cos(alturasolar));,→ 173 174 // si la cantidad es positiva el ángulo es negativo 175 if(sinazimut > 0) azimut = - azimut; 176 177 return resultado = new double[] {alturasolar, azimut}; 178 } 179 180 private static double[][] CalculoVectorPosicion2DEspejos(double alturasolar, double azimut, double[] XE, int Ytubo),→ 181 { 182 double[][] resultado; 183 double[] i2d, i3d, r2d0, r2d1; 184 double oriplanta = Fresnel.oriplanta; 185 double gamma, modi2d; 186 int numeroFilasEspejos, j; 187 188 numeroFilasEspejos = XE.Length; // This is not changing (for this Solar plant),→ 189 resultado = new double[6][]; 190 r2d0 = new double[numeroFilasEspejos]; 191 r2d1 = new double[numeroFilasEspejos]; 192 i2d = new double[2]; 193 i3d = new double[3]; 194 195 // Modelo 2D, cálculo del vector i2d solar 196 // No usamos vectores, sino arrays 197 if(azimut > 0) 198 { 199 i2d = new double[] {Math.Cos(alturasolar)*Math.Cos(azimut + oriplanta), Math.Sin(alturasolar)};,→ 200 i3d = new double[] {i2d[0], i2d[1], Math.Cos(alturasolar)*Math.Sin(azimut + oriplanta)};,→ 201 } 202 else if(azimut<0&& Math.Abs(azimut)<=oriplanta) 203 { 204 i2d = new double[] {Math.Cos(alturasolar)*Math.Cos(oriplanta - Math.Abs(azimut)), Math.Sin(alturasolar)}; ,→ ,→ 205 i3d = new double[] {i2d[0], i2d[1], Math.Cos(alturasolar)*Math.Sin(oriplanta - Math.Abs(azimut))};,→ 206 } 207 else if(azimut<0&& Math.Abs(azimut)>oriplanta) 208 { 209 i2d = new double[] {Math.Cos(alturasolar)*Math.Cos(Math.Abs(azimut) - oriplanta), Math.Sin(alturasolar)}; ,→ ,→
68 Capítulo A. Código comentado 210 i3d = new double[] {i2d[0], i2d[1], Math.Cos(alturasolar)*Math.Sin(Math.Abs(azimut) - oriplanta)};,→ 211 } 212 213 modi2d = Math.Sqrt(Math.Pow(i2d[0], 2) + Math.Pow(i2d[1], 2)); 214 215 for(j=0; j < i3d.Length; j++) 216 { 217 i3d[j] = i3d[j]*(1/Math.Sqrt(Math.Pow(i3d[0], 2) + Math.Pow(i3d[1], 2) + Math.Pow(i3d[2], 2)));,→ 218 } 219 220 gamma = Math.Atan(i2d[1]/i2d[0]); 221 222 if(alturasolar <= 0) gamma = 0;//gamma=alturasolar; 223 224 // Cálculo del vector r2d de posiciones de espejos. 225 for(j=0; j < numeroFilasEspejos; j++) 226 { 227 r2d0[j] = - XE[j]/Math.Sqrt(Math.Pow(XE[j], 2) + Math.Pow(Ytubo, 2));,→ 228 r2d1[j] = Ytubo/Math.Sqrt(Math.Pow(XE[j], 2) + Math.Pow(Ytubo, 2));,→ 229 } 230 231 resultado[0] = i2d; 232 resultado[1] = i3d; 233 resultado[2] = r2d0; // r2d.CopyTo(resultado, 2); 234 resultado[3] = r2d1; 235 resultado[4] = new double[] {modi2d}; 236 resultado[5] = new double[] {gamma}; 237 238 return resultado; 239 } 240 241 private static double[][] CalculoInclinacionEspejos(double[] XE, double[] i2d, double[,] r2d, double modi2d, double gamma),→ 242 { 243 // Cálculo del ángulo que forma cada espejo con la horizontal beta.,→ 244 245 double[][] resultado; 246 double[] alpha, theta, beta, kappa, inclinacion, factorsombra;,→ 247 int numeroFilasEspejos = XE.Length, j; 248 249 resultado = new double[6][]; 250 alpha = new double[numeroFilasEspejos]; theta = new double[numeroFilasEspejos];,→
A.2 FresnelSupport.cs 69 251 beta = new double[numeroFilasEspejos]; kappa = new double[numeroFilasEspejos];,→ 252 inclinacion = new double[numeroFilasEspejos]; factorsombra = new double[numeroFilasEspejos];,→ 253 254 for(j=0; j < numeroFilasEspejos; j++) 255 { 256 beta[j] = Math.Atan(r2d[1, j]/Math.Abs(r2d[0, j])); 257 } 258 259 // Cálculo del ángulo alpha y theta 260 for(j=0; j < numeroFilasEspejos; j++) 261 { 262 alpha[j] = Math.Acos((r2d[0, j]*i2d[0] + r2d[1, j]*i2d[1])/modi2d);,→ 263 theta[j] = alpha[j]/2; 264 } 265 266 // Distintos casos, para cada espejo, primero para las 267 // primeras 5 filas. Cálculo de inclinación 268 for(j=0; j < beta.Length; j++) 269 { 270 inclinacion[j] = 1; 271 factorsombra[j] = 0;// Inicializacion de factorsombra. Lo usamos mas adelante.,→ 272 kappa[j] = 1;// Inicializacion de kappa. Lo usamos mas adelante.,→ 273 } 274 275 for(j=0; j<5; j++) 276 { 277 if(gamma >= beta[j]) 278 { 279 inclinacion[j] = Math.PI/2- beta[j] - theta[j]; 280 } 281 if(gamma < beta[j]) 282 { 283 inclinacion[j] = Math.PI/2- beta[j] + theta[j]; 284 } 285 } 286 inclinacion[5] = theta[5]; 287 for(j=6; j < numeroFilasEspejos; j++) 288 { 289 inclinacion[j] = beta[j] + theta[j] - Math.PI/2; 290 } 291 292 resultado[0] = alpha; 293 resultado[1] = theta; 294 resultado[2] = beta; 295 resultado[3] = kappa;
70 Capítulo A. Código comentado 296 resultado[4] = inclinacion; 297 resultado[5] = factorsombra; 298 299 return resultado; 300 } 301 302 private static double[] CalculoFactorSombra(double[] XE, 303 double[] inclinacion, 304 double Anchuraespejos, 305 double gamma, 306 double[] factorsombra) 307 { 308 // Una vez calculados los ángulos necesarios, se calculan 309 // los factores de sombra y la eficiencia geométrica 310 311 double[] D, E; 312 double Pii, pj, dx, dy, modulo, w, a, b; 313 int j, k; 314 315 for(k=0;k<10; k++) // ¿Por que se hacen 10 iteraciones? 316 { 317 j=k+1; 318 Pii = inclinacion[k]; 319 pj = inclinacion[j]; 320 321 // analizamos los casos: 322 // 1ºSi ambas inclinaciones son positivas 323 // 2ºSi la primera positiva y la segunda negativa 324 // 3ºSi las dos son negativas 325 if(Pii > 0&&pj>0) 326 { 327 D=new double[] {XE[j] - Anchuraespejos*Math.Cos(pj), Anchuraespejos*Math.Sin(pj)};,→ 328 E=new double[] {XE[k] + Anchuraespejos*Math.Cos(Pii), -Anchuraespejos*Math.Sin(Pii)};,→ 329 330 // ángulo w 331 dx = Math.Abs(D[0] - E[0]); 332 dy = Math.Abs(D[1] - E[1]); 333 modulo = Math.Sqrt(Math.Pow(dx, 2) + Math.Pow(dy, 2)); 334 w = Math.Atan(dy/dx); 335 a = w - gamma; 336 b = Pii + gamma; 337 338 if(gamma < w) 339 { 340 factorsombra[k] = modulo*Math.Sin(a)/Math.Sin(b); 341 } 342 } 343 if(Pii > 0&&pj<0)
A.2 FresnelSupport.cs 71 344 { 345 D=new double[] {XE[j] + Anchuraespejos*Math.Cos(pj), Anchuraespejos*Math.Sin(pj)};,→ 346 E=new double[] {XE[k] + Anchuraespejos*Math.Cos(Pii), -Anchuraespejos*Math.Sin(Pii)};,→ 347 348 // ángulo w 349 dx = Math.Abs(D[0] - E[0]); 350 dy = Math.Abs(D[1] - E[1]); 351 modulo = Math.Sqrt(Math.Pow(dx, 2) + Math.Pow(dy, 2)); 352 w = Math.Atan(dy/dx); 353 a = w - gamma; 354 b = Pii + gamma; 355 356 if(gamma < w) 357 { 358 factorsombra[k] = modulo*Math.Sin(a)/Math.Sin(b); 359 } 360 } 361 if(Pii < 0&&pj<0) 362 { 363 factorsombra[k] = 0; 364 } 365 } 366 367 return factorsombra; 368 } 369 370 private static double[][] CalculoEficienciaGeometrica( double[] factorsombra,,→ 371 int Areatotal,,→ 372 double[] inclinacion,,→ 373 double[] kappa,,→ 374 double gamma),→ 375 { 376 double[][] resultado; 377 double[] eficienciageometrica; 378 double Pps; 379 int j; 380 381 eficienciageometrica = new double[kappa.Length]; 382 resultado = new double[2][]; 383 384 for (j = 0; j < factorsombra.Length; j++) 385 { 386 factorsombra[j] = Math.Abs(factorsombra[j])*32;
72 Capítulo A. Código comentado 387 } 388 Pps = 1- factorsombra.Sum()/Areatotal; 389 390 // Se halla ahora la parte del tubo que está 391 // en sombra. Angulo kappa 392 for(j=0;j<11; j++) 393 { 394 if(inclinacion[j] > 0) 395 { 396 kappa[j] = Math.PI/2- gamma - inclinacion[j]; 397 } 398 if(inclinacion[j] < 0) 399 { 400 kappa[j] = Math.PI/2+ inclinacion[j] - gamma; 401 } 402 } 403 for (j = 0; j < kappa.Length; j++) 404 { 405 eficienciageometrica[j] = Math.Abs(Math.Cos(kappa[j])); 406 } 407 408 resultado[0] = new double[] {Pps}; 409 resultado[1] = eficienciageometrica; 410 411 return resultado; 412 } 413 414 private static double CalculoParteSombreada(double[] inclinacion, 415 double alturasolar, 416 double[] i3d, 417 int longitudPlanta, 418 double Ytubo, 419 double[] eficienciageometrica,,→ 420 double Pps, double azimut, double[] XE),→ 421 { 422 double[,] ene, C; 423 double[] aux, ff, partesombra; 424 double Factorsombra, tubosombra, fulit, Eficienciatubosombra;,→ 425 int ZR, j; 426 int numeroFilasEspejos = XE.Length; 427 428 ene = new double[3, numeroFilasEspejos]; 429 C = new double[3, numeroFilasEspejos]; 430 ff = new double[numeroFilasEspejos]; 431 partesombra = new double[numeroFilasEspejos]; 432 433 // Modelo 3D, parte del tubo que está en sombra
A.2 FresnelSupport.cs 73 434 for(j=0; j < numeroFilasEspejos; j++) // Justificar por que 11 iteraciones aqui,→ 435 { 436 ene[0, j] = Math.Sin(inclinacion[j]); 437 ene[1, j] = Math.Cos(inclinacion[j]); 438 ene[2, j] = 0; 439 } 440 441 if(azimut > 0) ZR = longitudPlanta; 442 else ZR = 0; 443 444 for(j=0; j < numeroFilasEspejos; j++) 445 { 446 aux = new double[] {ene[0, j], ene[1, j], ene[2, j]}; 447 C[0, j] = (double)CrossProduct3D(i3d, aux).GetValue(0); 448 C[1, j] = (double)CrossProduct3D(i3d, aux).GetValue(1); 449 C[2, j] = (double)CrossProduct3D(i3d, aux).GetValue(2); 450 ff[j] = Math.Pow((C[0, j]/Math.Cos(inclinacion[j])), 2); // Explicar por que C[0],→ 451 partesombra[j] = ZR - Math.Sqrt((Math.Pow(XE[j], 2) + Math.Pow(Ytubo, 2))*ff[j]/(1- ff[j]));,→ 452 } 453 fulit = partesombra.Average(); 454 455 // Comprobamos si hay sombra efectivamente o no. 456 if(ZR == 0) 457 { 458 tubosombra = Math.Abs(fulit/longitudPlanta); 459 } 460 if(ZR == longitudPlanta) 461 { 462 tubosombra = (longitudPlanta - fulit)/longitudPlanta; 463 } 464 // Hay que poner una alternativa, pues no es posible que tubosombra no tome ningun valor.,→ 465 // Vamos a considerar que, si no se da ninguno de los dos casos previos, tubosombra = 0;,→ 466 else tubosombra = 0; 467 468 Eficienciatubosombra = 1- tubosombra; 469 470 if(Eficienciatubosombra > 1) 471 { 472 Eficienciatubosombra = 1; 473 } 474 if(Eficienciatubosombra < 0) 475 { 476 Eficienciatubosombra = 0; 477 } 478
74 Capítulo A. Código comentado 479 Factorsombra = Eficienciatubosombra*eficienciageometrica.Max()*Pps;,→ 480 if(Factorsombra < 0|| alturasolar < 0) 481 { 482 Factorsombra = 0; 483 } 484 if(Factorsombra > 1) 485 { 486 Factorsombra = 1; 487 } 488 489 return Factorsombra; 490 } 491 492 private static double[] CrossProduct3D(double[] a, double[] b) 493 {// Performs cross-product for 3D arrays 494 double[] result = new double[3]; 495 result[0] = a[1] * b[2] - a[2] * b[1]; 496 result[1] = a[2] * b[0] - a[0] * b[2]; 497 result[2] = a[0] * b[1] - a[1] * b[0]; 498 return result; 499 } 500 501 public static double[,] CalculoTemperatura( double[] TactualMetal_CalcTemp,,→ 502 double[] TactualFluido_CalcTemp,,→ 503 double Tambiente_CalcTemp,,→ 504 double I_CalcTemp, 505 double Factorsombra_CalcTemp,,→ 506 double Caudal_CalcTemp,,→ 507 double Eficienciamedia_CalcTemp,,→ 508 double Porcentajetubo) 509 // CalculoTemperatura recibe los parametros del sistema para cada instante de simulacion,→ 510 // Los dos primeros tienen un valor para cada discretizacion del tubo. Los demas, son escalares:,→ 511 // los suponemos constantes en toda la longitud del sistema. 512 { 513 // I en realidad es un array con valor para cada muestra. 514 double Af, Af2, Am, Am2, G, L, L2, pm, Cm, tint_CalcTemp, q; 515 double[] incx, Hl_local, pf_local, Cf_local, Ht_local, TfluidoSinCorregir;,→ 516 int n; 517 double[] Tfluido_CalcTemp, Tmetal_CalcTemp, TactualPerdidas, vectorradiacion;,→
A.2 FresnelSupport.cs 75 518 double[,] aux2; 519 double[,] TemperaturasCalcTemp; 520 521 // Debug.Log("Temperatura fluido 33: "+TactualFluido_CalcTemp[33]);,→ 522 523 Af = Fresnel.Af; Af2 = Fresnel.Af2; Am = Fresnel.Am; Am2 = Fresnel.Am2;,→ 524 G = Fresnel.G; L = Fresnel.L; L2 = Fresnel.L2; 525 pm = Fresnel.pm; Cm = Fresnel.Cm; 526 tint_CalcTemp = Fresnel.tint; 527 528 incx = new double[TactualFluido_CalcTemp.Length]; 529 for(n=0; n < TactualFluido_CalcTemp.Length; n++) 530 { 531 if(n<64) incx[n] = 1; 532 else incx[n] = 0.9; 533 } 534 535 // Características del fluido y coeficientes de pérdidas y transmisión de,→ 536 // calor correspondientes a cada trozo de metal y fluido 537 538 TactualPerdidas = new double[TactualFluido_CalcTemp.Length]; 539 for(n=0;n<64; n++) 540 { 541 TactualPerdidas[n] = TactualMetal_CalcTemp[n]; 542 } 543 for(n=64; n < TactualFluido_CalcTemp.Length; n++) 544 { 545 TactualPerdidas[n] = TactualFluido_CalcTemp[n]; 546 } 547 548 Hl_local = CalculoPerdidasMetal(TactualPerdidas, Tambiente_CalcTemp);,→ 549 550 aux2 = CalculoCaracteristicasDelFluido(TactualFluido_CalcTemp, Caudal_CalcTemp);,→ 551 pf_local = new double[aux2.GetLength(1)]; 552 Cf_local = new double[aux2.GetLength(1)]; 553 Ht_local = new double[aux2.GetLength(1)]; 554 for(n=0; n < aux2.GetLength(1); n++) // Quiero como limite del bucle el numero de columnas,→ 555 {// Ojo que aqui los datos se organizan por filas, y no por columnas (como es el caso del fichero Datos.csv),→ 556 pf_local[n] = aux2[0, n]; 557 Cf_local[n] = aux2[1, n]; 558 Ht_local[n] = aux2[2, n]; 559 } 560
82 Capítulo A. Código comentado 15 double[,] TemperaturasSimulacion; 16 double[] TemperaturaFluidoSimulacion, TemperaturaMetalSimulacion; 17 // double[] TemperaturaSalidaFluido, TemperaturaTuberia, TemperaturaEntradaFluido,,→ 18 // CaudalSalida, IrradianciaSalida; 19 double IrradianciaSimulacion, PorcentajeTuboSimulacion, TiempoSolarSimulacion, FactorSombraSimulacion,,→ 20 TemperaturaAmbienteSimulacion, qSimulacion, EficienciaMediaSimulacion,,→ 21 minutoSimulacion, horaSimulacion, RepeticionesSubcicloSimulacion;,→ 22 int j, i, n, Ntotal, contador10; 23 24 // Archivo donde almacenaremos los resultados de la temperatura de salida,→ 25 string filePath = @"C:\Users\Pablo\OneDrive - UNIVERSIDAD DE SEVILLA\TFG\FSCmodel\FSCmodelUnity\Assets\Scripts\Debugging_data_060923.txt";,→ 26 27 void Awake() 28 { 29 Debug.Log("Simulacion ha empezado"); 30 31 // Modifica el tiempo fijado de repetición del bucle. Para que tuviera sentido del todo,,→ 32 // habría que ejecutar los cálculos cada 10 iteraciones del bucle (así, el tiempo de,→ 33 // muestreo es realmente 0.25 s, como es en el modelo en matlab). De todas formas, los,→ 34 // resultados son aceptables sin llevar a cabo esta corrección sugerida.,→ 35 Time.fixedDeltaTime = 0.025f; 36 37 38 // Se pone en marcha el sistema: 39 Fresnel.Setup(); 40 41 TemperaturaFluidoSimulacion = Fresnel.Tf; 42 TemperaturaMetalSimulacion = Fresnel.Tm; 43 44 // A hora y minuto de la simulación se le dan los valores del primer elemento de las muestras,,→ 45 // porque a partir del primer paso, el mismo bucle de simulación las calculará.,→ 46 minutoSimulacion = Fresnel.minuto1[0]; 47 // Atencion al espíritu o la intención de la declaración en el modelo original: no es lo que estamos,→ 48 // haciendo: (original:) x = find(abs(hora-horacomienzo)==min(abs(hora-horacomienzo))); hora = hora(x); ,→ ,→ 49 horaSimulacion = Fresnel.hora1[0];
A.3 Simulacion.cs 83 50 51 EficienciaMediaSimulacion = Fresnel.Eficienciamedia; 52 j=0; 53 54 Ntotal = (int)Fresnel.Ntotal; 55 56 57 // TemperaturaSalidaFluido = new double[Ntotal]; 58 // TemperaturaTuberia = new double[Ntotal]; 59 // TemperaturaEntradaFluido = new double[Ntotal]; 60 // CaudalSalida = new double[Ntotal]; 61 // IrradianciaSalida = new double[Ntotal]; 62 63 contador10 = 1; 64 65 } 66 67 void FixedUpdate() 68 { 69 TemperaturaFluidoSimulacion[0] = Fresnel.Tentrada1[j]; // Inicializamos el primer valor (la temperatura de entrada del fluido) ,→ ,→ 70 IrradianciaSimulacion = Fresnel.I1[j]; // Es constante para todo el sistema,→ 71 TemperaturaAmbienteSimulacion = Fresnel.Tambiente1[j];// Es constante para todo el sistema,→ 72 qSimulacion = Fresnel.caudal1[j]; // Es constante para toda la tubería,→ 73 74 PorcentajeTuboSimulacion = 1; 75 // Si la temperatura de salida es mayor que 180ºC 76 if(TemperaturaFluidoSimulacion[Fresnel.numeroPartesDiscretasTuboCaptador] > (double)180) PorcentajeTuboSimulacion = 0; ,→ ,→ 77 78 if(qSimulacion < 0) qSimulacion = 0; 79 80 if(j % Fresnel.tactualizacion == 0) 81 { 82 minutoSimulacion = minutoSimulacion + (double)10 / (double)60;,→ 83 if(minutoSimulacion >= (double)60) 84 { 85 horaSimulacion++; 86 if(horaSimulacion == 24) horaSimulacion = 0; 87 minutoSimulacion -= (double)60; 88 } 89
84 Capítulo A. Código comentado 90 TiempoSolarSimulacion = FresnelSupport.CalculoHoraSolar(horaSimulacion, minutoSimulacion, Fresnel.mes, Fresnel.angulodiario1); ,→ ,→ 91 92 FactorSombraSimulacion = FresnelSupport.eficienciaGeoYSombras(Fresnel.angulodiario1, TiempoSolarSimulacion); ,→ ,→ 93 } 94 95 TemperaturasSimulacion = FresnelSupport.CalculoTemperatura( TemperaturaMetalSimulacion,,→ 96 TemperaturaFluidoSimulacion,,→ 97 TemperaturaAmbienteSimulacion,,→ 98 IrradianciaSimulacion,,→ 99 FactorSombraSimulacion,,→ 100 qSimulacion,,→ 101 EficienciaMediaSimulacion,,→ 102 PorcentajeTuboSimulacion);,→ 103 104 for(int indice = 0; indice < TemperaturasSimulacion.Length/2; indice++),→ 105 { 106 // El vector conjunto tiene longitud 2xlongitud de uno de los vectores de temperaturas,→ 107 TemperaturaMetalSimulacion[indice] = TemperaturasSimulacion[0, indice];,→ 108 TemperaturaFluidoSimulacion[indice] = TemperaturasSimulacion[1, indice];,→ 109 } 110 111 // Store data for later plotting 112 // In each iteration, append the value to the file 113 File.AppendAllText(filePath, TemperaturaFluidoSimulacion[64].ToString() + Environment.NewLine);,→ 114 // Guardar también los valores de temperatura del metal 115 116 j++; 117 118 if(j > Ntotal) UnityEditor.EditorApplication.isPlaying = false;,→ 119 } 120 }
A.4 SolarMovement.cs 85 A.4 SolarMovement.cs 1using System; 2using System.Collections; 3using System.Collections.Generic; 4using UnityEngine; 5using UnityEditor; 6 7public class SolarMovement : MonoBehaviour 8{ 9// Struct to save the position of the sun and the mirrors 10 SolarMovementSupport.FresnelSystem fresnelSystem; 11 float rotacionActualSol, rotacionActualEspejo, anguloRotacionEspejo, anguloRotacionAnterior;,→ 12 float rotacionYplanta, rotacionZplanta; 13 14 // Awake to find every mirror and the position of the sun 15 void Awake() 16 { 17 fresnelSystem = SolarMovementSupport.FindFresnelObjects(); 18 rotacionYplanta = (float)GameObject.Find("Planta").transform.localRotation.eulerAngles.y;,→ 19 rotacionZplanta = (float)GameObject.Find("Planta").transform.localRotation.eulerAngles.z;,→ 20 } 21 22 // Fixed Update to sims every instant the movement of the sun and the automatic control of every mirror,→ 23 void FixedUpdate() 24 { 25 fresnelSystem.Sun.Rotate(SolarMovementSupport.variacionMovimientoSolar, 0,0); ,→ ,→ 26 27 rotacionActualSol = (float)GameObject.Find("Sun").transform.localRotation.eulerAngles.x;,→ 28 29 // Lo hacemos por espejos, aunque para lo que estamos haciendo daría lo mismo hacerlo por filas,,→ 30 // con idea de implementarlo con más precisión en un futuro. 31 for(int i = 0; i < fresnelSystem.Mirrors.GetLength(0); i++) 32 { 33 anguloRotacionEspejo = (float)SolarMovementSupport.MirrorAngle(i + 1, anguloRotacionAnterior); ,→ ,→ 34 rotacionActualEspejo = (float)fresnelSystem.Mirrors[i, 0].transform.localRotation.eulerAngles.x;,→ 35
86 Capítulo A. Código comentado 36 // Recorremos todas las filas y todas las columnas para girar los espejos,→ 37 for(int j = 0; j < fresnelSystem.Mirrors.GetLength(1); j++),→ 38 fresnelSystem.Mirrors[i, j].transform.rotation = Quaternion.Euler(anguloRotacionEspejo, rotacionYplanta, rotacionZplanta); ,→ ,→ 39 40 // Debug.Log("RotacionEspejo: "+anguloRotacionEspejo); 41 42 anguloRotacionAnterior = anguloRotacionEspejo; 43 } 44 45 if(GameObject.Find("Sun").transform.localRotation.eulerAngles.x > 180) UnityEditor.EditorApplication.isPlaying = false; ,→ ,→ 46 } 47 } A.5 SolarMovementSupport.cs 1using System; 2using System.Collections; 3using System.Collections.Generic; 4using UnityEngine; 5using UnityEditor; 6 7public class SolarMovementSupport : MonoBehaviour 8{ 9// Parameters 10 const int numberOfMirrorsInLine = 16; 11 const int fresnelCollectorHeight = 4;//metros 12 public const float variacionMovimientoSolar = 0.5f;//grados 13 //public static double anguloEspejo; 14 static int flagSol = 0; 15 16 // Struct to save the position of the sun and the mirrors 17 public struct FresnelSystem 18 { 19 public Transform Sun; 20 public Transform[,] Mirrors; 21 22 } 23 24 // Function that find every mirror and the position of the sun 25 public static FresnelSystem FindFresnelObjects() 26 { 27 GameObject[] aux;
A.5 SolarMovementSupport.cs 87 28 int i = 0,j=0; 29 30 FresnelSystem fresnelSystem; 31 fresnelSystem.Sun = GameObject.Find("Sun").transform; 32 aux = GameObject.FindGameObjectsWithTag("Mirror"); 33 fresnelSystem.Mirrors = new Transform[aux.Length/numberOfMirrorsInLine, numberOfMirrorsInLine]; ,→ ,→ 34 35 foreach(GameObject mirror in aux) 36 { 37 fresnelSystem.Mirrors[i, j] = mirror.transform; 38 j++; 39 if(j == numberOfMirrorsInLine) 40 { 41 i++; 42 j=0; 43 } 44 } 45 46 return fresnelSystem; 47 } 48 49 public static double MirrorAngle(int rowNumber, float anguloRotacionAnterior),→ 50 { 51 // La funcion recibe el numero de fila del espejo, sabiendo que empezamos a contar por el espejo que en Unity,→ 52 // se corresponde con Fila1 53 double anguloEspejo, h, angulo1; 54 55 int coeficiente = Math.Abs(6- rowNumber); 56 57 // Calculo de la longitud del vector que une el centro del espejo con el tubo colector,→ 58 h = Math.Sqrt(Math.Pow(fresnelCollectorHeight, 2) + Math.Pow(coeficiente*0.7,2));,→ 59 60 if(Math.Asin(coeficiente*0.7 /h)<=90) 61 angulo1 = 180*Math.Asin(coeficiente*0.7 / h)/Math.PI; 62 else 63 angulo1 = 180 -180*Math.Asin(coeficiente*0.7 / h)/Math.PI;,→ 64 65 double angulo2 = 180 -90 - angulo1; 66 double angulo3 = 90 + angulo1; 67 68 // Antes de hallar el angulo que debe girar el espejo, necesito saber qué angulo forma el sol,→
88 Capítulo A. Código comentado 69 double anguloSol = GameObject.Find("Sun").transform.localRotation.eulerAngles.x;,→ 70 71 if(flagSol == 1)if(anguloSol < 90) anguloSol = 180 - anguloSol;,→ 72 73 if(rowNumber > 5) anguloEspejo = (angulo2 + anguloSol)/2; 74 else anguloEspejo = (angulo3 + anguloSol)/2; 75 76 if(anguloSol == 90) flagSol = 1; 77 78 return -anguloEspejo; 79 } 80 81 } A.6 CSVReader.cs 1using UnityEngine; 2using System; 3using System.Linq; 4using System.IO; 5using System.Collections; 6using System.Collections.Generic; 7 8public static class CSVReader 9{ 10 public static string[,] LeerCSV(string filePath) 11 { 12 int i = 0, j; // Contadores para recorrer el array donde vamos a almacenar los valores.,→ 13 string[] values; // Valores de las cadenas de cada línea del archivo .csv.,→ 14 string[,] valor; // Array 2D de strings que contiene los valores del .csv organizados en matriz.,→ 15 16 // Primero calculamos el número de líneas del archivo. En el caso,→ 17 // de nuestro proyecto, eso equivale al número de muestras reales.,→ 18 int numberOfLines = File.ReadLines(filePath).Count(); 19 // Ahora identificamos la primera línea. 20 var firstLine = File.ReadLines(filePath).FirstOrDefault(); 21 // Con el objetivo de calcular el número de columnas que tiene el,→ 22 // archivo. Es decir, el número de variables medidas, en 23 // nuestro caso. 24 var numberOfColumns = firstLine.Split(',').Length;
A.6 CSVReader.cs 89 25 // Inicializamos las dimensiones de nuestra variable "valor". 26 valor = new string[numberOfLines, numberOfColumns]; 27 28 if (File.Exists(filePath)) 29 { 30 using (StreamReader reader = new StreamReader(filePath)) 31 { 32 // Mientras queden datos por leer: 33 while (!reader.EndOfStream) 34 { 35 // Leemos la línea en curso, 36 string line = reader.ReadLine(); 37 // separamos los valores que contiene según las comas (importante que estén puestas),,→ 38 values = line.Split(','); 39 // reseteamos j a 0 antes de entrar en el bucle de asignación.,→ 40 j=0; 41 // Como hemos separado los valores por las comas, ahora tenemos tantos strings,→ 42 // como columnas tenga el archivo. 43 // Así que, cogemos cada string y lo asignamos a una posición del array 2D.,→ 44 foreach (string value in values) 45 { 46 valor[i, j] = value; 47 j++; 48 } 49 i++; 50 } 51 } 52 } 53 else 54 { 55 // Si no existe el archivo cuyo path nos han pasado, sacamos error por pantalla.,→ 56 Debug.LogError("File not found: " + filePath); 57 } 58 59 return valor; 60 } 61 }
Índice de Figuras 1.1 Principio de funcionamiento de un captador cilindro-parabólico 3 1.2 Principio de funcionamiento de un captador Fresnel 3 1.3 Vista general de la planta 4 1.4 Dibujo de concepto de la planta con algunas medidas 5 1.5 Vista general de la planta 6 2.1 Pantalla de selección del tipo de proyecto en Unity Hub 8 3.1 Fotograma renderizado de una escena promocional de Zibra Liquids Pro 12 3.2 Ejemplo de aproximación Lagrangiana (basada en partículas) 13 3.3 Concepto de simulación Euleriana de fluidos (basada en malla o matriz) 13 3.4 13 3.5 Las partículas no respetan las superficies "complejas" 13 3.6 Estructura del modelo en Matlab 15 3.7 Resultados del modelo PDE en Matlab para las muestras del día 29 de junio de 2009 18 3.8 Error del modelo PDE en Matlab para las muestras del día 29 de junio de 2009 19 3.9 Estructura del modelo en C# 19 3.10 Esquema de las funciones del modelo óptico 24 3.11 Representación gráfica del azimut y la altura solares 28 3.12 Comparación de los resultados de los modelos en Matlab y C# con las medidas 45 3.13 La componente fija del error es prácticamente nula 45 3.14 A medida que avanza la simulación, la temperatura calculada por el modelo se aleja más de la salida real 46 3.15 Al final de la simulación el error deja de ser despreciable 46 3.16 Modelo de los espejos de la planta Fresnel del tejado de la ETSI 47 3.17 Espejo con superficie reflectante renderizado 47 3.18 Sistema de coordenadas de un espejo en el modelo 3D 48 3.19 Un espejo insertado en la estructura del modelo 48 3.20 Estructura de barras del modelo 49 3.21 Patas de la planta 49 3.22 Barras laterales 50 3.23 Cubierta del tubo receptor 50 3.24 Esquema de relaciones de los objetos del modelo 51 3.25 Cada espejo del está etiquetado para ser identificado por el script 51 3.26 Vista del inspector que muestra los scripts de Simulacion.cs y SolarMovement.cs asociados al objeto Planta 52 91