Documento 6-06-09

ÍNDICE GENERAL DE LA TESIS DOCTORAL CAPÍTULO I. PRESENTACIÓN Y OBJETIVOS I.1 Antecedentes e introducción I.2 Objetivos y desarrollo I.3 Perspectivas CAPÍTULO II. FUNDAMENTOS TEÓRICOS II.1 Introducción II.2 Transmisión de calor por conducción II.2.1 La ecuación de Fourier o ecuación de conducción II.2.2 Medios multicapas II.2.3 Conducción de calor en aletas II.2.3.1 Introducción II.2.3.2 Ecuación diferencial de la aleta 1D. Hipótesis simplificadoras II.2.4 Condiciones iniciales y de frontera II.2.5 Soluciones analíticas II.3 El Método de simulación por redes II.3.1 Idea del método. Tipos de monopuertas II.3.2 El MESIR como método numérico II.4 Modelos en red II.4.1 La analogía termoeléctrica clásica (antecedentes) II.4.2 Modelos para medios homogéneos y para condiciones de contorno II.5 El software PSpice-Orcad II.5.1 Introducción y aplicaciones II.5.2 Simulación. Presentación de resultados II.6 El programa C# CAPÍTULO III. EL PROGRAMA PROCCA-09 III.1 Introducción III.2 Estructura del programa PROCCA-09 III.3 Creación de archivos de modelos III.3.1 Presentación del programa III.3.1.1 Pantallas del módulo CONCBA III.3.1.2 Pantallas del módulo CONCAL III.4 Criterios para la numeración de celdas, nodos y elementos del modelo III.5 Estructura de los archivos de texto de modelos III.6 Pantallas de presentación de resultados III.7 Ejemplos de archivos de modelo III.7.1 Placa rectangular con una oquedad III.7.2 Cilindro hueco de dos capas III.7.3 Aleta rectangular 1-D CAPÍTULO IV. APLICACIONES DOCENTES Y DE INVESTIGACIÓN IV.1 Introducción IV.2 Aplicaciones docentes IV.2.1 Transitorios en medios 1-D y 2-D de geometría rectangular Placa 1-D bajo enfriamiento convectivo Placa 2-D con condiciones de contorno armónicas IV.2.2 Placas 1-D bajo condiciones de convección y radiación IV.2.3 Estudio de aletas simples IV.3 Aplicaciones de investigación IV.3.1 Aislamiento de tanques esféricos multicapas IV.3.2 Introducción al diseño de aletas y conjuntos aleta-pared rectangulares CAPÍTULO V. CONTRIBUCIONES CAPÍTULO VI. REFERENCIAS

1

I. Presentación y objetivos

Capítulo I

PRESENTACIÓN Y OBJETIVOS

I.1

Antecedentes e introducción

2

I.2

Objetivos y desarrollo

5

I.3

Perspectivas

8

2

I. Presentación y objetivos

I.1

ANTECEDENTES E INTRODUCCIÓN

Durante los años xxx-xxxx de mi estancia como profesor TEU en la Escuela de Ingeniería Industrial de la actual Universidad Politécnica en Cartagena, como profesor del Área de Máquinas y Motores Térmicos, la elaboración de una tesis doctoral en el campo de la termodinámica o en el de la transmisión de calor era, casi, una tarea impensable pues no existían doctores ni laboratorios adecuados en esa escuela para abordar este objetivo. La solución era investigar en otras Universidades, Murcia o Valencia, lo que para muchos suponía un esfuerzo extraordinario. Poco antes de mi traslado a la Universidad de Murcia se formó en la UPCT el grupo de investigación “Simulación por Redes”, dirigido por el profesor Carlos González, con el objetivo inicial de aplicar las técnicas del Método de redes a problemas de transmisión del calor. Yo empecé a interesarme en esta técnica y fruto de ello fueron mis primeros trabajos en este campo que se tradujeron en publicaciones en congresos nacionales e internacionales y en revistas científicas, Alcaraz y col. [2001], Alhama y col. [2001], Alarcón y col [2001] y Zueco y col. [2001], Alarcón y col. [2002a, 2003a] Zueco y col. [2002b] y Alcaraz y col [2003b]. Tras varios años en los que mi tarea fundamental en la UMU derivó hacia las Ciencias de la Educación, sin abandonar mi trabajo docente en el campo termodinámico ni abandonar mi contacto con el grupo de Cartagena y mi intención de hacer mi tesis doctoral, surge un punto de encuentro: combinar las Nuevas Tecnologías con la docencia propia de mi área. La elaboración de esta memoria es el fruto de ese encuentro. [Paco: Incluir un párrafo sobre los aspectos educativos de la tesis (necesidad de una enseñanza técnica moderna en la que el ordenador juegue un papel esencial); incluir tres referencias bibliográficas como máximo), 10 líneas en total como máximo para este párrafo]. Con todo, en el año 2004, el responsable del grupo de investigación me invitó a trabajar en la elaboración de dos programas básicos, uno para la conducción de calor en medios multicapas rectangulares y otro para el cálculo y diseño de aletas simples, con el objetivo de que ambos fueran el germen de una posible memoria de tesis doctoral tras extenderlos a otras geometrías y ampliar los objetivos en el campo general de conducción de calor, materia del área de conocimiento a la que pertenezco. En 2005 salieron las primeras versiones educativas de estos programas, el segundo de los cuales ya se ha registrado en su versión definitiva, “PRODASIM 1” [,2007] y el primero, PESCAR está en fase de registro. Como fruto de este trabajo se han presentado varias comunicaciones a congresos, nacionales e internacionales y en revistas científicas, tanto de carácter técnico como educativo, Del Cerro Velázquez y col. [2004a, 2004b, 2004c, 2004d, 2005a, 2005b, 2008]. Durante este periodo también he desarrollado trabajos con el método de simulación por redes dentro del grupo sobre modelos específicos aplicados a educación y otros de carácter científico, Del Cerro Velázquez y col. [2044e, 2004f], Del Cerro Velázquez y Alhama [2004g, 2004h].

3

I. Presentación y objetivos

Tras la finalización de estos programas, cuyos derechos de explotación fueron cedidos a la UPCT, el director de esta memoria me invitó a trabajar en un programa más amplio y refinado que constituyera el núcleo de una tesis doctoral. Tras marcar los objetivos de este programa, de nombre PROCCA-09, ahora en fase de registro para su explotación, después de algunos años de trabajo, hemos concluido con la realización de la presente memoria. La conducción de calor es, sin duda, uno de los campos más relevantes en ingeniería pues son muchos los procesos industriales en los que se da este fenómeno físico. Conocer los aspectos esenciales del mismo es, pues, un objetivo fundamental de la docencia en las escuelas de ingeniería, en particular en la carrera de ingeniero químico. Para cubrir este objetivo, además de una adecuada preparación del profesorado y una planificación adecuada de la enseñanza, es preciso el concurso del material didáctico necesario que incluye los textos de estudio de teoría y problemas, las prácticas de laboratorio y, cuando es posible, el uso de programas de ordenador específicos. Las ventajas de esta última herramienta, en contraposición con las prácticas de laboratorio, son muchas: es generalmente barata y no tiene riesgos, portátil, permite el autoaprendizaje, se adapta al ritmo personal de cada alumno, se auto corrige, etc. En esta línea, la presente tesis doctoral se orienta a la elaboración de un programa de conducción de calor que cubra los objetivos de docencia en este campo, dentro del programa exigido por el ministerio en sus orientaciones troncales para las diferentes titulaciones de grado de todas las ingenierías. El programa de conducción de calor, “PROCCA-09” [2009], que se constituye entonces en el objetivo fundamental de esta tesis, se estructura en dos módulos: el primero, “CONCBA” (módulo básico de conducción de calor), de carácter general, se orienta a la resolución numérica de la ecuación lineal transitoria del calor en diferentes geometrías y medios homogéneos (monocapa) y heterogéneos (multicapas), 1-D y 2-D, con condiciones de contorno arbitrarias (lineales o no) y con extensiones particulares que permiten la localización de dominios anómalos en la región estudiada; el segundo, “CONCAL (módulo de aplicación al cálculo y optimización de aletas)”, se orienta específicamente al diseño y optimización de algunos tipos de aletas y conjuntos aleta-pared rectangulares, un tema esencial en el diseño de equipos térmicos. Aunque algunos de los problemas lineales estudiados con estos programas tienen solución analítica, si bien frecuentemente mediante series complicadas de convergencia lenta (particularmente en geometrías cilíndricas y esféricas), la solución presentada por los programas se basa en cálculos numéricos bien contrastados. La herramienta usada a este fin es el método de simulación por redes (MESIR) que hace uso de la similitud entre el transporte de calor y el de carga eléctrica, si bien no se trata de la conocida analogía termoeléctrica contenida en muchos libros de texto (Holman [1981], Incropera y de Witt [1981], Chapman [1984], Grigull y Sandner [1984], Kakaç y Yener [1985], Özisik [1993], Lienhard [1987], White [1988], Brodkey 4

I. Presentación y objetivos

y Hershey [1988], Bayazitoglu y Özisik [1988], Thomas [1992] y Mills [1995]) cuyo objetivo, además de orientarse exclusivamente a problemas lineales, es frecuentemente académico o formal, resolviendo problemas en los que la analogía se establece entre partes finitas del medio en el que tiene lugar la transmisión del calor que, a su vez, están conectadas en serie o paralelo con otras partes del mismo u otro medio. En los años 60 la analogía termo-eléctrica llegaba a implementar físicamente los circuitos y tuvo un cierto auge pues los cálculos numéricos no se habían desarrollado suficientemente. Sin embargo, tal simulación analógica tuvo no pocas dificultades que propiciaron su abandono como herramienta de cálculo; entre ellas, los relativamente elevados márgenes de tolerancia en la fabricación de componentes y su elevado número en problemas complejos. El enorme desarrollo en capacidad y velocidad que experimentaron los equipos informáticos hizo que estos desbancaran sin esfuerzo las técnicas analógicas de simulación. El Método de simulación por redes aplicado a los problemas de conducción térmica, objeto de los programas presentados en esta memoria, pretende aunar todo el potencial existente en la analogía termo-eléctrica con la potencia de los modernos ordenadores. En este sentido, el software empleado para la resolución de los modelos en red tiene las ventajas de trabajar con dispositivos ideales, disponer de amplias librerías de componentes, aportar soluciones con errores tan pequeños como se soliciten y requerir tiempos de ejecución relativamente pequeños. Además el número de reglas requeridas para elaboración de los archivos (modelos) contenidos en los programas objeto de esta tesis es muy pequeño ya que, por un lado, son muy pocos los componentes eléctricos diferentes que necesitan y por otro, el tipo de operaciones y resultados requeridos son siempre los mismos (determinación de temperaturas y flujos de calor en el medio). En cuanto a los criterios para el desarrollo del programa presentado en esta memoria, mencionaremos que se ha desarrollado en un ambiente (o interfaz), bajo Windows, agradable para el usuario al que sólo se le exige unas nociones fundamentales de conducción de calor. La introducción sucesiva y ordenadas de datos, la manipulación y tratamiento de los archivos de programa generados, el acceso al módulo de cálculo y la obtención y tratamiento de los resultados gráficos y numéricos son aspectos que se han elaborado cuidadosamente y con la información de ayuda necesaria. El programa dispone, así mismo de herramientas autocorrectoras que informan al usuario de los límites en los valores de los rangos numéricos y de posibles errores en la introducción de datos. Por último diremos que el presente trabajo se enmarca dentro de las líneas de investigación del grupo “Simulación por Redes” de la Universidad politécnica de Cartagena y pretende abrir un nuevo frente en sus objetivos: la elaboración de programas docentes y de investigación basados en los modelos diseñados por este grupo y aplicados a diferentes campos de la ciencia e 5

I. Presentación y objetivos

ingeniería: transmisión de calor, flujo de fluidos con transporte de solutos, mecánica del rozamiento, elasticidad, etc. I.2

OBJETIVOS Y DESARROLLO

El objetivo de esta memoria es la elaboración de un amplio programa, “PROCCA-09”, de ordenador para la resolución de problemas transitorios y estacionarios de conducción de calor en medios lineales, incluido el cálculo y optimización de aletas, que pueda ser aplicado tanto al campo docente como de investigación. El software de programación usado es C# (Archer [2004]), un código que aprovecha las ventajas de la programación en C++. La presentación bajo entorno Windows del programa, muy intuitiva y con ayudas, permite al usuario seguir los pasos lógicos, introducción de datos para la elaboración del modelo asociado al problema, comprobación y/o modificaciones del modelo, simulación en el correspondiente software numérico (PSpice [1994]) y tratamiento y presentación de los resultados en forma tabulada o gráfica, usando para ello el programa MATLAB [1997] consecución del problema y presentación de las soluciones en forma tabulada o gráfica. La finalidad del programa es doble: i) docencia, problemas ya organizados y otros propuestos ii) investigación, problemas nuevos y diseño de equipos. Como se ha mencionado en la introducción, el programa se ha estructurado en dos módulos denominados, CONCBA y CONCAL. El primer módulo, CONCBA, es de carácter fundamental y permite la elaboración y simulación de modelos para los siguientes problemas: ⇒ Problemas transitorios y estacionarios, 1-D, en geometría rectangular, cilíndrica (variable el radio del cilindro) y esférica (variable el radio de la esfera), multicapas (con capas de materiales homogéneos) y condiciones de contorno arbitrarias. Pueden adoptarse una geometría hueca tanto para las coordenadas cilíndricas como para las esféricas. ⇒ Problemas transitorios y estacionarios, 2-D, en geometría rectangular y cilíndrica (con el radio y la altura del cilindro como variables), multicapas (con capas de materiales homogéneos) y condiciones de contorno arbitrarias. Puede adoptarse también una geometría hueca para el cilindro. Las condiciones de contorno son arbitrarias y pueden extenderse a cada una de las cuatro regiones que definen el contorno o a dominios más pequeños (y hasta puntuales) dentro de estas regiones. ⇒ En relación con las condiciones de contorno pueden adoptarse las siguientes: -

Isoterma

-

Temperatura dependiente del tiempo (rectangular, sinusoidal, polinomial, a tramos, pulso…) 6

I. Presentación y objetivos

-

De flujo constante

-

De flujo dependiente del tiempo (rectangular, sinusoidal, polinomial, a tramos, pulso…).

-

Adiabática

-

De convección

-

De radiación

-

De convección y radiación simultáneas

⇒ Con las condiciones dadas, el programa es capaz de realizar los cálculos requeridos para el diseño de paredes multicapas para la optimización de recintos de aislamiento térmico tales como tanques industriales, paredes de edificación, recintos refrigerados, etc. El módulo específico de elaboración y simulación de modelos para el cálculo y optimización de aletas, CONCAL, es de carácter aplicado y permite abordar los siguientes problemas: ⇒ Problemas transitorios y estacionarios, 1-D, de aletas simples de cualquier sección y condiciones de contorno arbitrarias, ⇒ Problemas transitorios y estacionarios, 1-D, de aletas anulares (geometría cilíndrica hueca), ⇒ Problemas transitorios y estacionarios, 2-D, de conjuntos aleta-pared recta rectangular, con condiciones de contorno arbitrarias. Se incluye el caso de la aleta rectangular con pared desnuda, es decir, el caso de aletas rectas rectangulares 2-D, ⇒ En relación con las condiciones de contorno, se pueden aplicar las mismas que se han mencionado para el otro módulo (CONCBA). ⇒ Con las posibilidades mencionadas es posible abordar el diseño de todos los conjuntos mencionados, es decir resolver el problema de optimización: dado un volumen de material encontrar el tipo de aleta que permite la máxima transferencia de calor (o, dado el dato de calor disipado, encontrar la menor cantidad de material de aleta que permite tal disipación, El funcionamiento de los diferentes subprogramas que interactúan entre si, PROCCA-09 (para la elaboración de modelo), PSpice (para simulación numérica y presentación de resultados) y MATLAB (para presentaciones gráficas alternativas) esta coordinado por el programa principal, PROCCA-09, el cual, mediante los comandos oportunos es capaz de arrancar la simulación, en sus diferentes opciones, cuando sea preciso, ejecutar el entorno gráfico de salida de datos de

7

I. Presentación y objetivos

PSpice y solicitar y tratar sus salidas numéricas de datos para hacer representaciones complementarias a conveniencia mediante MATLAB. Las rutinas creadas a tal fin se incorporan en el programa principal, PROCCA-09. En este sentido se constituyen en objetivos parciales del programa principal los siguientes: ⇒ Creación de las subrutinas necesarias para corregir los modelos, ⇒ Creación de un archivo de texto que se integrará como un programa general de ayuda integrado en el programa principal, ⇒ Creación de las subrutinas de arranque de programas auxiliares, tratamiento de datos, representaciones gráficas y otras opciones, ⇒ Creación de las opciones de simulación que permitan ejecutar simulaciones múltiples para optimización y tratar los datos en estos casos, Se resume, a continuación el desarrollo de los Capítulos de esta memoria. En el Capítulo II se presentan los aspectos teóricos y metodológicos que hemos considerado necesarios para enmarcar adecuadamente los contenidos de la misma. En primer lugar los aspectos relativos a los fundamentos de conducción de calor: las ecuaciones de transporte lineales en diferentes geometrías, 1-D y 2-D, bajo la hipótesis de Fourier. En segundo lugar los aspectos asociados a la aplicación de la técnica numérica usada para la simulación, el método de simulación por redes (MESIR), exponiendo brevemente los contenidos de esta técnica. En el mismo capítulo se incluyen tanto los fundamentos del programa PSpice usado para la ejecución (simulación) de los modelos como una introducción a la programación en C# y al programa de cálculo matemático MATLAB, usado para ciertas representaciones gráficas. El Capítulo III se dedica a la exposición del programa objeto de esta memoria, PROCCA-09. Se muestra en detalle el funcionamiento del programa en todos los aspectos: definición del problema, entrada de datos, creación del modelo, modificaciones, opciones de simulación y búsqueda de resultados. Se exponen todas las posibilidades de PROCCA-09 aunque creemos sinceramente que es posible arrancar otras opciones al programa no incluidas en este capítulo, particularmente fruto de las posibilidades de las nuevas versiones de PSpice en las que se incorporan nuevas opciones de simulación múltiple. El Capítulo IV se dedica a las aplicaciones docentes y de investigación de PROCCA-09. Dado que el enfoque docente como programa para realizar prácticas de conducción de calor es uno de los objetivos del programa, se exponen dos aplicaciones que pueden desarrollarse como prácticas guiadas de laboratorio. También se incluyen dos aplicaciones, una de diseño de paredes multicapas y otra de optimización de aletas, que por su contenido podrían considerarse objetivos de investigación. Las contribuciones de esta memoria se recogen al final, en un apartado específico. 8

I. Presentación y objetivos

I.3

PERSPECTIVAS

El camino abierto con la elaboración del programa PROCCA-09 es, de hecho, una nueva línea de trabajo del grupo de investigación “Simulación por redes” de la UPCT a la que dedicaremos un gran esfuerzo en los próximos años. Creemos que la elaboración de programas puede convertirse en una actividad interesante que, sobre la base de que el diseño previo es la actividad fundamental de dicho grupo, permite cerrar el trabajo del grupo haciéndolo completo. La elaboración de otros programas de transmisión de calor que incluyan sistemas de ecuaciones acopladas (movimiento de fluidos con transporte de calor simultáneo), problemas de radiación entre paredes, protocolos para la resolución de problemas inversos, optimización multivariable, etc., junto con programas de otro tipo que han sido o son objetivos de investigación del grupo, tales como flujo de aguas subterráneas con transporte de soluto, modelos de rozamiento, problemas de elasticidad, caos, etc., proporciona un campo de trabajo casi ilimitado y muy fructífero. Es aquí donde tengo el propósito de proyectar mi futuro trabajo de investigación.

9

II. Fundamentos teóricos

Capítulo II

FUNDAMENTOS TEÓRICOS

II.1

Introducción

II.2

Transmisión de calor por conducción

II.2.1

La ecuación de Fourier o ecuación de conducción

II.2.2

Medios multicapa

II.2.3

Conducción de calor en aletas

II.2.3.1

Introducción

II.2.3.2

Ecuación diferencial de la aleta 1D. Hipótesis simplificadoras

II.2.4

Condiciones iniciales y de frontera

II.2.5

Soluciones analíticas

II.3

El Método de simulación por redes

II.3.1

Idea del método. Tipos de monopuertas

II.3.2

El MESIR como método numérico

II.4

Modelos en red

II.4.1

La analogía termoeléctrica clásica (antecedentes)

II.4.2

Modelos para medios homogéneos y para condiciones de contorno

II.5

El software PSpice-OrCAD

II.5.1

Introducción y aplicaciones

II.5.2

Simulación. Presentación de resultados

II.6

El programa C#

10

II. Fundamentos teóricos

II.1 INTRODUCCIÓN Los objetivos cubiertos en este capítulo son varios. Se exponen, en primer lugar, los fundamentos de la ciencia de conducción del calor, particularizando al caso de los problemas en medios lineales objeto del programa de esta memoria. Se presenta la ecuación transitoria de conducción, 1-D y 2-D, en medios homogéneos y heterogéneos multicapas, en las geometrías rectangular, cilíndrica y esférica, así como las ecuaciones de conducción en aletas simples, aletas anulares y en conjuntos aleta-pared rectangulares. Se completan los modelos matemáticos con las condiciones de contorno habituales, que incluyen la condición no lineal de convección, radiación y convección más radiación. Las soluciones analíticas, a estos problemas de conducción de calor, cuando existen, están habitualmente constituidas por desarrollos en serie de convergencia lenta, particularmente en geometrías no rectangulares. A continuación, se introducen brevemente los fundamentos del Método de Simulación por Redes (MESIR), herramienta de cálculo para la simulación numérica de los problemas que el programa propuesto es capaz de abordar. El MESIR es un método versátil y potente, muy extendido en la literatura científica, capaz de modelar, en principio, cualquier problema matemático definido mediante un conjunto de ecuaciones de gobierno y de condiciones de contorno. Los modelos usados en esta memoria que se introducen en este capítulo, han sido ya publicados en la literatura científica por investigadores del grupo de simulación por redes y han demostrado, en efecto, ser suficientemente precisos ya que los errores quedan reducidos a valores del 0,5 ó 1 % para problemas transitorios lineales, 1-D y 2-D, con mallados relativamente pequeños (Alhama [1999], Alarcón [2001] y Alarcón y col. [2000a y 2000b]). La aplicación del MESIR precisa de un programa de resolución de circuitos eléctricos. De entre los existentes en el mercado se ha adoptado OrCAD (derivado del antiguo software PSpice [1994]) en sus diferentes versiones. Un objetivo de este capítulo se dedica a las posibilidades de análisis y simulación de este programa. Por último, se hace una introducción al software comercial C#, código de programación usado para la elaboración del programa.

II.2

TRANSMISIÓN DE CALOR POR CONDUCCIÓN

II.2.1 La ecuación de Fourier o ecuación de conducción Para un medio homogéneo e isótropo, es decir, de conductividad independiente de la posición y de la dirección espacial, la relación entre la densidad de flujo calorífico, j (W/m2), energía por segundo que cruza la unidad de superficie isoterma, y el gradiente térmico ∇ T (K/m), vector normal a la superficie o línea isoterma, viene dada por la ley de Fourier (introducida con anterioridad por Biot [Biot, 1816]),

11

II. Fundamentos teóricos

j(r, t) = - k∇[T(r, t) ]

(2.1)

donde r es la posición (m), t el tiempo (s) y k la conductividad térmica (W K-1 m-1), magnitud escalar positiva. Se trata de una ley fenomenológica, es decir, no deducida de principios fundamentales sino derivada a partir de observaciones directas. Como es sabido, la conducción es uno de los tres modos de transmisión del calor en el que el intercambio energético tiene lugar en sólidos, o en fluidos en reposo (ausencia de movimientos convectivos resultantes del desplazamiento de porciones macroscópicas del medio), debido a la presencia de una diferencia de temperaturas. El balance energético local en un medio en reposo entre la energía almacenada (energía térmica interna), la energía en tránsito (calor) y la generada o consumida por el propio medio (fuentes o sumideros), permite ser formulado matemáticamente en términos de una ecuación en derivadas parciales denominada ecuación de conducción del calor, Figura 2.1. Para un medio homogéneo e isótropo dicha ecuación tiene la forma

( ∇ ⋅ j) + g ( r, t ) = ρ ce  ∂ T 

(2.2)

 ∂t 

Figura 2.1 Balance energético local

En esta ecuación g(r,t), W/m3, es la densidad de potencia (potencia generada por unidad de volumen en un punto del medio de posición r), ρ es la densidad (kg m-3), y ce el calor específico por unidad de masa a presión constante (J kg-1 K-1). Sustituyendo la densidad de flujo de calor de la ecuación (ecuación (2.1)) en la ecuación anterior, la ecuación de conducción puede expresarse en términos de una única variable dependiente, la temperatura, y dos independientes, la posición y el tiempo:  ∂T  ∇⋅ [ k∇T( r, t ) ] + g ( r, t ) = ρ ce    ∂t 

(2.3)

Si se trata de un medio lineal (conductividad y calor específico constantes, es decir, independientes de la temperatura y la posición) la ecuación de conducción se reduce a

12

II. Fundamentos teóricos

∇2T ( r, t ) + g(r, t) = α −1

∂T ( r, t ) ∂t

(2.4)

donde α = k/(ρ ce) es la llamada difusividad térmica (m2/s), asociada con la rapidez global de propagación del calor en el medio. En ausencia de generación interna ni sumideros de calor la ecuación se simplifica a

 ∂T ( r, t )  ∇2T ( r, t ) = α −1    ∂t 

(2.5)

En coordenadas rectangulares, Figura 2.2a, la ecuación de conducción tiene la forma ∂  ∂T  ∂  ∂T  ∂  ∂T  ∂T k  k +  + ∂z k ∂z  = ρ c e ∂t ∂x  ∂x  ∂y  ∂ y    

(2.6a)

o bien ∂2T ∂2T ∂2T 1 ∂T + + 2 = 2 2 α ∂t ∂x ∂y ∂z

(2.6b)

En coordenadas cilíndricas, Figura 2.2b, la ecuación de conducción tiene la forma 1 ∂ ∂T  1 ∂  ∂T  ∂  ∂T  ∂T k  k r +  + ∂z k ∂z  = ρ c e ∂t r ∂r  ∂r  r 2 ∂φ  ∂ φ    

(2.7a)

o bien 1 ∂  ∂T  1 ∂2T ∂2T 1 ∂T + = r + 2 2 2 r ∂r  ∂r  r ∂φ α ∂t ∂z

(2.7b)

mientras que el coordenadas esféricas, Figura 2.2c, la forma de la ecuación es 1 ∂  2 ∂T  1 ∂ ∂T  1 ∂  ∂T  ∂T k  k r + 2 k sen θ + 2 2 2  = ρ ce ∂t ∂r  r sen θ ∂θ  ∂θ  r sen θ ∂φ  ∂ φ r ∂r   

(2.8a) o bien

13

II. Fundamentos teóricos

1 ∂  2 ∂T  1 ∂ ∂T  1 ∂2T 1 ∂T r + sen θ + =     2 2 2 2 2 ∂θ  r sen θ ∂φ α ∂t r ∂r  ∂r  r sen θ ∂θ 

(2.8b)

14

II. Fundamentos teóricos

Figura 2.2a Geometría rectangular

Figura 2.2b Geometría cilíndrica

Figura 2.2c Geometría esférica

Los problemas abordados estudiados en esta memoria abordan sólo tratan con geometrías 2-D rectangulares, 2-D cilíndricas, con r y z como variables independientes, y 1-D esféricas, con r como variable espacial. Así, las ecuaciones anteriores se simplifican para estas coordenadas, respectivamente, a las siguientes ∂  ∂T  ∂  ∂T  ∂  ∂T  ∂T k  k +  + ∂z k ∂z  = ρ ce ∂t ∂x  ∂x  ∂y  ∂ y    

(2.9a)

1 ∂  ∂T  ∂  ∂T  ∂T k r + k  = ρ ce r ∂r  ∂r  ∂z  ∂z  ∂t

(2.9b) 1 ∂  2 ∂T  ∂T k r  = ρ ce 2 ∂r  ∂t r ∂r 

(2.10a)

o bien, en medios con características térmicas constantes, ∂2T ∂2T 1 ∂T + = 2 2 α ∂t ∂x ∂y

(2.10b)

15

II. Fundamentos teóricos

1 ∂  ∂T  ∂2T 1 ∂T r + 2 = r ∂r  ∂r  ∂z α ∂t

(2.11a)

16

II. Fundamentos teóricos

1 ∂  2 ∂T   1  ∂T  r  =    r 2 ∂r  ∂r   α  ∂t 

(2.11b)

La deducción de estas ecuaciones puede encontrarse en numerosos libros de texto entre los que cabe mencionar, por su carácter pedagógico, a Chapman [1984], Özisik (1993), Bejan [1995], Mills [1995] y Incropera y DeWitt [1996].

II.2.2 Medios multicapa La geometría de estos medios, mostrada en la Figura 2.3 para el caso de geometría cilíndrica hueca 2-D, es frecuente en multitud de procesos de transmisión de calor tales como: i) equipos propios de la industria de producción de frío y calor, ii) recintos aislantes de almacenamiento de productos alimenticios y gases y líquidos industriales, iii) aislamiento térmico de edificaciones, etc. Los medios suelen estar formados por varias capas una de las cuales suele ser la capa térmicamente no conductora (o aislante), situada en la zona central, mientras que las otras suelen tener otras propiedades mecánicas o químicas requeridas por el diseño. En cada capa se satisface la ecuación de calor, (2.9b), (2.10b) ó (2.11b), según la geometría correspondiente, rectangular, cilíndrica y esférica, respectivamente, de forma que el conjunto global de ecuaciones, es de la forma

∂ 2Ti

+

∂ 2Ti

 1   ∂T  =    i   αi   ∂t 

i = 1, 2, …n

(2.12)

1 ∂  ∂Ti  ∂ 2Ti 1 ∂T  ri  + = ⋅ i 2 ri ∂ri  ∂ri  ∂zi αi ∂t

i = 1, 2, …n

(2.13)

1 ∂  2 ∂Ti  1 ∂Ti  ri = ⋅ 2 ri ∂ri  ∂ri  αi ∂t

i = 1, 2, …n

(2.14)

∂x i

2

∂yi

2

donde ‘n’ es el número de capas. Además de estas ecuaciones deben cumplirse las condiciones de frontera entre capas; en el caso de contacto térmico perfecto estas ecuaciones son la de conservación del flujo calorífico entre medios y la de continuidad del valor de la temperatura al pasar de un medio a otro, es decir

− k1∂T − k 2 ∂T pared = pared ∂n ∂n

(2.15)

T ( x i = Li ) = T( x i +1 = 0)

(2.16)

17

II. Fundamentos teóricos

donde Li es el tamaño de capa i (en esta ecuación el origen de coordenadas se ha tomado al principio de cada capa).

18

II. Fundamentos teóricos

Figura 2.3 Medio multicapas. Geometría cilíndrica hueca 2-D

II.2.3 Conducción de calor en aletas II.2.3.1 Introducción Dentro del campo de la ingeniería, las aletas son elementos adicionales que se adosan a la superficie de un cuerpo cuando se desea eliminar calor de éste. Pueden ser del mismo o distinto material que la pared del cuerpo al que están adosadas. Las aletas forman parte esencial de dispositivos tan variados y comunes como intercambiadores de calor, compresores, motores térmicos o eléctricos refrigerados por aire, disipadores de calor en dispositivos eléctricos y electrónicos, etc. La adición de la aleta trae como consecuencia un aumento del área por la que se intercambia calor entre el cuerpo y el medio. Sin embargo, dado que el flujo de calor ha de atravesar mayor cantidad de material (el material de la propia aleta) se produce un aumento de la resistencia térmica. Así, aunque en muchas aplicaciones las aletas se emplean para disipar calor, es decir, para aumentar la transmisión de calor hacia el entorno, más frío, también pueden realizar la función inversa, es decir, aumentar la ganancia térmica de un objeto. En general, la relativamente gran longitud relativa y pequeño espesor de la aleta proporcionan una gran superficie de contacto con el fluido que baña la superficie exterior del cuerpo, superficie a través de la cual se disipa el calor que entra en la aleta por su base. El mecanismo más frecuente de intercambio térmico a lo largo de toda esta superficie exterior es la convección, en las distintas formas que ésta adopta. En este caso se habla de aletas convectivas.

19

II. Fundamentos teóricos

En otros casos, cuando el salto térmico es importante o no existe fluido exterior, la radiación puede ser el mecanismo de disipación. También puede ser importante el flujo disipado por la pared desnuda por lo que es necesario considerar ésta en muchas aplicaciones. La terminología de las aletas rectas de sección constante se muestra en la Figura 2.4. Las aletas longitudinales o rectas se caracterizan porque la base de la aleta es plana y se extiende a lo largo de la superficie primaria también plana o casi-plana. Se denomina base de la aleta a la superficie en contacto con la pared sobre la que descansa la aleta o superficie primaria; el extremo o punta de la aleta es la superficie más alejada de la base; los flancos o superficies laterales son las principales superficies de disipación de calor y los bordes son las superficies que cierran el volumen de la aleta, cuyo efecto sobre la transmisión de calor es despreciable en ciertos tipos de aletas longitudinales. Se incluyen también en los programas propuestos aletas de sección cuadrada, circular y sección arbitraria, Figura 2.5. Las dimensiones principales de las aletas rectas rectangulares son la longitud o altura, La, el espesor, 2e, y la anchura, b. Se denomina sección transversal a la sección en el plano YZ, plano perpendicular al flujo térmico, mientras que la sección longitudinal es la del plano XY que define el perfil de la aleta. El plano XZ es, en este caso, un plano de simetría en aletas rectas y contiene al eje de la aleta de dirección axial, eje X, dirección principal en la que fluye el calor. En relación con las aletas tipo aguja, Figura 2.5, el programa incluye todos los tipos de aleta, cuya denominación atiende a la forma de la sección transversal: rectangulares, cuadradas, circulares, etc.

bordes

flancos

y

x

base

z

2e

b

La

punta

Figura 2.4 Nomenclatura de la aleta recta longitudinal

20

II. Fundamentos teóricos

b)

a)

c) Figura 2.5 Aletas tipo aguja: Fig. 2.5a Aleta rectangular; Fig. 2.5b Aleta circular; Fig. 2.5c Aleta de sección elíptica.

El modelo físico más frecuente que se encuentra en la literatura técnica es el de aleta aislada (Kern y Kraus [1972]; Chapman [1984], etc.) Figura 2.6, es decir una aleta que se considera independiente de la pared soporte, asumiendo en su base una temperatura conocida, dependiente o no del tiempo, u otra condición de contorno tipo convectivo, de radiación. y x

T∞

2e

b

z

Tb La

Figura 2.6 Aleta longitudinal o recta aislada de perfil rectangular o de sección constante

21

II. Fundamentos teóricos

II.2.3.2 Ecuación diferencial de la aleta 1D. Hipótesis simplificadoras Los primeros estudios formales provienen de principios del siglo pasado, Harper y Brown, [1922], y abarcan un significativo número de geometrías. Estos estudios clásicos se basan en las llamadas hipótesis de Murray-Gardner (Gardner [1945]), o hipótesis simplificadoras (limiting asumptions) (Kern y Kraus [1972]). Aún hoy se asumen la mayor parte de estas hipótesis, tanto en los libros de texto especializados como en gran parte de la literatura científica, para el diseño práctico de estos conjuntos. Estas hipótesis son: i)

el calor a través de la aleta y su temperatura permanecen constantes a lo largo del tiempo (proceso estacionario)

ii) iii)

el medio es uniforme e isótropo, con conductividad térmica constante el coeficiente de convección es constante y uniforme en toda la superficie exterior de la aleta

iv) la temperatura del medio exterior que rodea la aleta es uniforme v)

el espesor es tan pequeño comparado con la longitud o altura de la aleta que es despreciable el gradiente térmico en la dirección del espesor de la aleta OY (unidimensional en la dirección axial de la aleta)

vi)

la temperatura en la base de la aleta es conocida y uniforme, lo que implica la consideración de la aleta aislada y despreciando su posible interacción con la pared a la que está unida, o sin formar un conjunto con otras aletas

vii)

no existe resistencia térmica de contacto entre la aleta y la pared que la sustenta o superficie primaria

viii)

no existen fuentes o sumideros internos de calor

ix)

aleta suficientemente ancha (dirección OZ) como para despreciar los efectos de borde

x)

el calor intercambiado por la aleta es proporcional a la sobretemperatura entre la aleta y el medio exterior (no existe radiación),

Otra condición que a menudo se impone en muchos estudios es: xi) extremo adiabático Las aletas objeto de estudio en esta tesis pueden asumir, pues, una condición complementaria que transciende de estas hipótesis clásicas, la condición de radiación tanto en la superficie de la aleta como en el extremo de la misma. Esta condición, además, puede simultanearse o no con la de convección y su inclusión impide la obtención de soluciones analíticas o semianalíticas del problema. Además, los modelos aleta recta rectangular y conjunto aleta pared recta rectangular son 2-D por lo que no necesitan satisfacer la hipótesis aludida en relación con el valor del número de Biot. Algunas de estas hipótesis (conducción 1-D, efecto de la pared y sólo disipación por convección) no son necesariamente asumidas por el programa que es capaz de dar soluciones a problemas menos restrictivos. 22

II. Fundamentos teóricos

El modelo unidimensional proporciona resultados suficientemente precisos para el cálculo y diseño de aquellas aplicaciones en las que el número de Biot, Bi =

ha2e , está por debajo de ka

0,1 y se cumplen sustancialmente las hipótesis clásicas, Irey [1968]. La ecuación de conducción, 1-D, para aletas rectas (Figura 2.7), viene dada por

d 2θ − m 2θ = 0 2 dx donde m =

(2.17)

hC es el parámetro de la aleta (inverso de la longitud característica). kA

L

x

x

C =π d A=

π d2 4

L

C=2 A=ωx1

Figura 2.7 Aleta recta de espesor uniforme o aguja de sección transversal constante

Para aletas anulares, Figura 2.8, esta ecuación es d 2θ 1 + 2 r dr

dθ dr

−

2h θ =0 kw

23

(2.18)

II. Fundamentos teóricos

re rb

tf

W

h

r

W

k

t0

δ

Figura 2.8 Aleta anular de espesor constante

A estas ecuaciones hay que añadir las condiciones de contorno. Por último, para aletas rectas 2D las ecuaciones de gobierno son ∂2θ ∂x 2

+

 ∂θ −k   ∂y 

∂θ = 0, ∂y θ =1,

∂θ = 0, ∂y

∂2θ ∂y 2

= 0;

  =hθ 

y = 0,

0 < x < l;

y = e;

0 < y < e;

(2.19)

0 < z < l;

(2.20a)

0 < z < l;

(2.20b)

x = 0;

(2.20c)

0 < y < e,

0 < y < e,

x = l;

(2.20d)

o  ∂θ  − k  ∂y   = h θ,  

0 < y < e,

x = l;

(2.20e)

En estas ecuaciones se ha utilizado la temperatura adimensional θ =

( T − Tref )

( Tb − Tref )

donde Tb es la temperatura de la base y T ref la del fluido externo. La ecuación (2.20d) representa una condición adiabática, caso particular de la ecuación (2.20e) que es una condición convectiva. Para los conjuntos aleta-pared a estas ecuaciones hay que añadir las de conducción 2-D en la pared y sus condiciones de contorno, ecuaciones similares a las expuestas anteriormente.

II.2.4 Condiciones iniciales y de frontera La ecuación de conducción requiere, para su solución, la especificación de unas condiciones de contorno asociadas a los límites del medio o fronteras y de unas condiciones iniciales que proporcionan la información sobre el campo inicial de temperaturas en el medio y el origen de

24

II. Fundamentos teóricos

la coordenada tiempo. Estas condiciones se expresan mediante ecuaciones matemáticas simples para los problemas estudiados en esta memoria.

25

II. Fundamentos teóricos

Si llamamos a la superficie de la frontera A, se define la densidad de flujo calorífico de conducción (W/m2) que cruza la superficie A como j n ( x frontera ) =

−k ∂T | pared . Esta densidad ∂n

de flujo de calor, de acuerdo con el tipo de interacción entre el medio y el exterior, puede ser debida a diferentes causas. Entre éstas cabe distinguir: jfuente, asociada a una fuente incidente externa. jconv, representa fuentes o sumideros de tipo convectivo debida al contacto de la superficie

(

)

A con un fluido exterior, j conv = h Tpared − T∞ , donde Tpared es la temperatura en la pared, T∞ la del fluido exterior lejos de la superficie y h el llamado coeficiente de transferencia de calor o coeficiente de película; h depende, en general, del tipo de flujo (laminar, turbulento, etc.), de la geometría y aspecto superficial de la pared, de la diferencia de temperaturas “Tpared - T∞”, de las propiedades físicas del fluido, etc. El rango de valores de este coeficiente es muy amplio y en muchos textos de ingeniería pueden encontrarse tabulaciones de h para cada tipo de aplicación. jrad, representa las pérdidas (o ganancias si jrad > 0) de calor por radiación hacia el medio exterior,

donde ε

es la emisividad de la superficie, σ la

[ ( ) − (T ) ] constante de Stefan-Boltzmann y T la temperatura de referencia del medio exterior para j rad = ε σ Tpared

4

r

4

r

la radiación. El balance energético a través de la superficie A es de la forma jn(xfrontera) = jfuente + jconv + jrad que usando las expresiones anteriores, puede escribirse en función de una única variable, la temperatura en la pared, mediante la ecuación no lineal:

(

)

[(

−k ∂T(r) | T pared = − j fuente + h Tpared − T∞ + ε σ Tpared ∂n

) 4 − ( Tr ) 4 ]

(2.21)

Con todo, las condiciones de contorno se clasifican en los textos habituales de acuerdo con la terminología siguiente: Condición de contorno de primera clase: Se especifica la temperatura en la frontera como una función generalmente dependiente de la posición y el tiempo, T ( rfrontera, t ) . El caso

T ( rfrontera, t ) = 0 se llama condición de contorno homogénea de primera clase. Condición de contorno de segunda clase: Se conoce el flujo externo en la superficie de la pared, que debe coincidir con el flujo calorífico de conducción,

26

−k ∂T(r) | T pared . El caso ∂n

II. Fundamentos teóricos

especial

∂T(r) | T pared = 0 , condición homogénea de contorno de segunda clase, ∂n

corresponde a una pared adiabática.

27

II. Fundamentos teóricos

Condición de contorno de tercera clase: Se trata de una condición convectiva; el flujo de conducción en la superficie de la pared coincide con el de convección del fluido exterior,

(

)

− k ∂T |pared = h Tpared − T∞ . El caso T∞ = 0 se llama condición homogénea de contorno ∂n

de tercera clase. Las condiciones de primera y segunda clase pueden derivarse de esta última haciendo k = 0 ó h = 0, respectivamente. Condición de radiación: Algunos autores (Mills [1955]), se refieren a ésta como condición de cuarta clase cuando la emisividad, ε , y la temperatura de referencia, Tr, son constantes. La condición de radiación

ecuación lineal

− k ∂T( x ) |pared = + ε ∂ Tpared ∂n

[(

(

) − ( T ) ] puede aproximarse a la 4

4

r

)

− k ∂T |T pared = hr Tpared − Tr , con hr = 4ε σTr 3 , cuando entre las ∂n

temperaturas Tpared y Tr se cumple

(T

pared

−Tr

Tr

)

< <1.

Condición de contacto entre medios: Se da cuando dos medios de igual o diferente conductividad se someten a un contacto térmico imperfecto. El balance energético en la frontera

requiere

que

se

cumpla

una

ecuación

lineal

del

tipo

− k1∂T − k 2∂T |pared,1 = h1,2 Tpared,1 − Tpared,2 = |pared,2 , donde h1,2 es la llamada ∂n ∂n

(

)

conductancia térmica de contacto entre medios, que depende de la rugosidad de las superficies, tipo de material, presión de contacto, material de separación de los medios, etc. Para un contacto perfecto, conductancia infinita, la ecuación anterior se convierte en

− k1∂T − k 2 ∂T |pared = |pared . ∂n ∂n

II.2.5 Soluciones analíticas Muchos de los problemas que pueden simularse con el programa de conducción propuesto en esta tesis tienen solución analítica, particularmente los problemas con condiciones de contorno lineales, aunque estas soluciones sean de difícil tratamiento por tratarse de series matemáticas de lenta convergencia que han de ser resueltas con programas de cálculo matemático. Remitimos a los textos Carslaw y Jeager [1959], Chapman [1984], Incropera y DeWit [1981] y Özisik [1993] para estas soluciones. Para otros problemas no lineales pueden encontrarse, en ocasiones, soluciones semianalíticas y en cualquier caso todos pueden tratarte numéricamente. La herramienta numérica con la que se 28

II. Fundamentos teóricos

resuelven los problemas del programa, el MESIR, es muy precisa ya que ha sido contrastada en numerosas aplicaciones tanto en el campo de transmisión de calor como en otros campos, por lo que sus soluciones pueden considerarse fiables. II.3

EL MÉTODO DE SIMULACIÓN POR REDES

II.3.1 Idea del método. Tipos de monopuertas El Método de simulación por redes (MESIR o NSM, Network Simulation Method) (GonzálezFernández y col. [2002]) es una técnica para el estudio de diferentes procesos físicos que puedan definirse mediante un conjunto de ecuaciones, o modelo matemático. Partiendo de éstas el procedimiento consiste, en primer lugar, en elaborar un “modelo en red” o circuito eléctrico equivalente al proceso, y en segundo lugar en simular dicho proceso, obteniendo la solución del modelo mediante un programa adecuado de resolución de circuitos eléctricos. Un modelo en red se considera equivalente a un determinado proceso cuando, en su descripción, las ecuaciones del modelo matemático discretizadas, y las ecuaciones del modelo en red para un elemento del volumen o celda elemental, correspondientes a variables análogas, coinciden. La técnica de elaboración del modelo consiste en reticular el espacio en elementos de volumen o celdas elementales; al aplicar a estas reticulaciones las ecuaciones diferenciales, se obtiene un conjunto de ecuaciones en diferencias finitas que se constituyen en el punto de partida para la obtención del modelo en red correspondiente a cada celda elemental; una seleccionada correspondencia entre variables dependientes del problema y variables eléctricas, tensiones e intensidades, permite interpretar los resultados de la simulación en términos del proceso que se modela. La asociación de celdas, de acuerdo con la geometría del problema, configura el modelo en red correspondiente a todo el medio finito, que es tanto más preciso cuanto mayor sea el número de estas celdas. Las condiciones de contorno e iniciales se incorporan al modelo de manera simple mediante dispositivos eléctricos adecuados. En el caso de los procesos de transmisión de calor, la posibilidad de elaborar modelos en red representativos de los mismos, es decir, el hecho de que admitan redes eléctricas equivalentes, supone no sólo la equivalencia matemática sino, también, la equivalencia física entre las variables características de unos y otros procesos (térmicos y eléctricos). Además la equivalencia física permite, en casos muy concretos, determinar cualitativa y cuantitativamente ciertas magnitudes asociadas a la red que pueden jugar un papel en la descripción del fenómeno de transporte, similar al correspondiente en el transporte de carga eléctrica, como es el caso de la impedancia, por ejemplo. A partir del modelo matemático y siguiendo los planteamientos de lo que se conoce como “Teoría de Redes” de Peusner [1986,1987], se obtiene un grafo equivalente al proceso cuya simulación (solución) se lleva a cabo mediante un adecuado programa de ordenador. Con el método de simulación por redes, el grafo es del tipo eléctrico (un circuito eléctrico), y la 29

II. Fundamentos teóricos

simulación se realiza mediante el software de simulación de circuitos. En este trabajo se ha utilizado PSpice (PSpice [1994], Nagel [1975 y 1977] y Vladimirescu [1994]).

30

II. Fundamentos teóricos

El método de simulación por redes presenta diferencias notables respecto de los métodos numéricos clásicos. Desde el punto de vista conceptual supone la sustitución de un complicado sistema de ecuaciones diferenciales, que ya no es necesario manipular, por un circuito eléctrico equivalente, de cuya solución se encarga PSpice, y donde resulta fácil visualizar la interconexión entre flujos y fuerzas, y relacionar los procesos físicos locales con la evolución de las variables en los componentes eléctricos que simulan el medio (modelo en red). En cuanto a la reticulación, sólo requiere una división de la variable espacio, como en los llamados métodos de líneas (Berezin y Zhidkov [1965] y Rukos [1978]). Por otro lado, en tanto que la continuidad de la corriente eléctrica (1ª ley de Kirchhoff) y la propiedad asociada al potencial eléctrico, derivado de un campo conservativo, (2ª ley de Kirchhoff) son leyes exigidas a los circuitos, no es necesario hacer las aproximaciones o tanteos que muchos métodos numéricos exigen para cumplir estos requerimientos; PSpice advierte cuando alguna de estas reglas no ha sido respetada en el diseño del modelo en red. Cabe señalar la ventaja que supone un buen conocimiento de la teoría de circuitos a la hora de implementar la red; no obstante es preciso poco esfuerzo para la familiarizarse con este aspecto del método, ya que son bastante reducidas las agrupaciones de términos de las expresiones matemáticas que se convierten en elementos o partes del circuito. En cuanto a la manipulación y elaboración del programa podemos afirmar que las dificultades son mínimas. Su presentación bajo “windows” (tanto para PC como para estación de trabajo) supone una guía constante para el manipulador, y la nomenclatura utilizada, organizada de nudos y componentes eléctricos, permite el acceso inmediato a las tensiones y corrientes que constituyen las variables fundamentales en el circuito. El MESIR, que utiliza la teoría de redes para modelar el proceso físico objeto de estudio, es un método de simulación en tanto que incluye la resolución numérica del modelo en red obtenido mediante la reticulación. Así, las variables flujo y fuerza características del mismo deben satisfacer las leyes de Kirchhoff y sus relaciones determinarán los elementos de circuito correspondientes. Ahora bien, en cada proceso concreto y una vez elegidas las variables conjugadas, la información de qué elementos de circuito intervienen en el modelo en red y cómo se conectan entre sí, se obtiene del modelo matemático y no de consideraciones de tipo físico acerca del papel que juegan estas variables. En síntesis, en la teoría de redes, la viabilidad de un modelo en red supone: i) La existencia de una red independiente del tiempo, ii) La existencia de una magnitud jN-N´ llamada flujo, asociada a cada rama que conecta los

nudos N-N´ y que va de N a N´. jN-N´ obedece las leyes de Kirchhoff para corrientes (LCK), iii) La existencia de una magnitud, φ , asociada a cada nudo, tal que la diferencia XN-N´ = φ

φ

, llamada fuerza, obedece la ley de los voltajes de Kirchhoff (LVK).

N´

31

N

-

II. Fundamentos teóricos

El hecho de que las leyes LCK y LVK se refieran, respectivamente, a corrientes eléctricas y voltajes no es significativo. En el caso de procesos de transporte es normal establecer una correspondencia entre variables flujo por un lado (densidad de corriente eléctrica con flujo de calor, flujo de masa, ...) y variables tipo potencial por otro (potencial eléctrico con temperatura, concentración, ...), pero es posible establecer otras analogías aún en procesos físicos que describan el transporte de una determinada magnitud. Las relaciones entre flujo y fuerza asociados a una rama y sus (dos) nudos límite, que pueden incluir o no variaciones temporales de estas variables que se dicen conjugadas, definen los elementos concretos del circuito equivalente a esa rama. La relación causa-efecto entre las variables conjugadas es completamente arbitraria con tal que sea consistente con ii) y iii). Volviendo a considerar un proceso de transporte caracterizado por las variables flujo y fuerza que obedecen las leyes de Kirchhoff y que, por consiguiente, son variables LCK y LVK respectivamente, tales leyes dan cuenta de la topología de la red relativa al proceso. A la red se le asocia un conjunto de flujos que obedecen a una ley de balance local y un conjunto de fuerzas que satisfacen la condición de unicidad. Las propiedades topológicas dependen únicamente de la asignación de conexiones entre los diferentes puntos o de las posibles combinaciones de trayectorias que unen un nudo dado con otros nudos. Son independientes de las medidas y, desde un punto de vista topológico, dos grafos son iguales o isomorfos si las asignaciones de vértices y ramas son las mismas. Las leyes de Kirchhoff establecen relaciones entre flujos y fuerzas por separado, pero no expresan ningún tipo de relación entre flujos y fuerzas entre sí. Las relaciones entre el par conjugado flujo-fuerza se conocen como ecuaciones constitutivas o fenomenológicas y definen los elementos de circuito que expresan características específicas de cada proceso. Se dice que dos grafos son geométricamente iguales si los potenciales y flujos de cada par de puntos y su rama correspondiente son iguales para cualquier conjunto de valores que puedan ser elegidos para los flujos o las fuerzas. Las propiedades geométricas de la red, es decir, sus características métricas, se siguen de las relaciones constitutivas. Las relaciones constitutivas se pueden establecer entre las variables de un par flujo-fuerza, en cuyo caso se habla de monopuerta. Las monopuertas pueden ser pasivas, que disipan o almacenan energía, y activas o fuentes, que generan potencia de acuerdo con una ley preestablecida. Las monopuertas pasivas son las resistencias, bobinas y condensadores mientras que las activas son las fuentes de tensión y/o corriente, controladas o no. En función de la relación expresa existente entre las variables LCK y LVK las monopuertas pasivas tienen nombre específicos. Las usadas en esta memoria son:

32

II. Fundamentos teóricos

Monopuerta resistiva. Es un elemento de circuito asociado a una relación entre las derivadas temporales de las variables flujo y fuerza de una misma rama, mediante una función independiente del tiempo que llamaremos resistencia, R, que puede depender o no del flujo o de la fuerza:

dX(t) dX(t) dJ(t) =R , o bien R = . dJ(t) dt dt

Una monopuerta resistiva es lineal cuando la relación entre las variables X(t) y J(t) lo es, es decir X(t) = R J(t); naturalmente R es una constante en este caso. Su acción es instantánea, no importa cual sea su estado anterior, en este sentido carecen de memoria. En su analogía física representan efectos disipativos, fricciones, efectos viscosos, energías de reacción, etc., y desde el punto de vista termodinámico son elementos generadores de entropía. Las monopuertas resistivas no lineales se definen a través de las funciones que las caracterizan y constituyen, en definitiva, fuentes controladas de corriente o tensión, respectivamente. Su representación simbólica de una monopuerta resistiva se muestra en la Figura 2.9. La traducción al modelo en red es una resistencia eléctrica de valor R ohmios para el caso lineal o una fuente controlada de corriente o tensión para el caso no lineal.

j(t)

jo

A B a) Lineal

b ) No lineal j(t)=F

(VAB )

R

c) No lineal V(t) =

v(t) + _

FR-1(jo)

Figura 2.9 Representación simbólica de monopuertas resistivas

Monopuerta capacitiva. Es un elemento de circuito asociado a una relación entre la variable flujo y la derivada temporal de la variable fuerza, de una misma rama, mediante una función no dependiente del tiempo que designaremos como capacidad, C, J(t) = C

dX(t) . En estas monopuertas se produce algún tipo de almacenamiento, sin dt

pérdidas (no hay disipación energética), y su estado, que no cambia instantáneamente, tiene en cuenta todas las operaciones llevadas a cabo en el pasado (tiene memoria). En su analogía, representa procesos físicos en los que se produce algún tipo de almacenamiento como condensadores, tanques, etc. La relación constitutiva anterior puede expresarse en términos de la capacidad

C=

dq dF (X) = c , que es constante cuando la dependencia q = Fc (X) es lineal, dX dX

33

II. Fundamentos teóricos

C=

q . Las dependencias q = Fc (X ) no lineales deben estudiarse particularmente en X

cada caso. Su representación simbólica de la monopuerta capacitiva lineal se muestra en la Figura 2.10. La traducción al modelo en red es un condensador eléctrico de valor C faradios.

34

II. Fundamentos teóricos

C

Figura 2.10 Representación simbólica de una monopuerta capacitiva lineal

Los procesos de almacenamiento y disipación de energía, bajo la hipótesis de continuidad del medio, se originan en todo los puntos del sistema. Los elementos R y C se identifican sin embargo con regiones pequeñas pero finitas del medio y sus conexiones con las otras puertas se realizan con enlaces ideales de energía, es decir, con conductores de resistencia nula. El que cada elemento pueda ser caracterizado por un par de variables conjugadas con una única ecuación constitutiva entre ellas es una hipótesis básica en el MESIR que deriva de la teoría de redes. Físicamente equivale a decir que es posible elegir un elemento de volumen lo suficientemente pequeño como para que su tiempo de relajación interna sea mucho menor que el del sistema global, pero suficientemente grande como para que las fluctuaciones de las variables que describe el sistema en él sean despreciables. En las monopuertas activas se produce una aportación o extracción de energía al sistema. Las empleadas dentro de los programas objeto de esta memoria son: Fuentes constantes. Son monopuertas definidas de acuerdo con las expresiones

FJ ( J ) = 0 y Fx ( X ) = 0 , según se trate de fuentes de flujo o de fuerza, respectivamente. Tienen asignado un sentido (o signo) que indica la dirección en que fluye la energía. La representación simbólica es la de la Figura 2.11a; eléctricamente, se corresponden a pilas o generadores de corriente constante. Fuentes dependientes del tiempo. La relación constitutiva entre las variables tiene la misma forma de las fuentes constantes; además, X =X(t ) y J =J(t) según se trate de fuentes de fuerza o de flujo. Ejemplos de representación simbólica se muestran en la Figura 2.11b. Fuentes controladas. Se trata de monopuertas especiales asociadas a relaciones constitutivas entre variables, conjugadas o no, expresadas mediante cualquier función que no contiene explícitamente el tiempo. Se trata de elementos de entradas múltiples con una única salida que corresponde a un flujo o una fuerza que depende funcionalmente de otros flujos o fuerzas de distintas ramas y nudos del mismo o diferente circuito. Estas fuentes van a permitir especificar acoplamientos energéticos de distinto tipo. Existen cuatro tipos de fuentes controladas por una sola variable: a) fuentes de tensión controladas por tensión, b) fuentes de tensión controladas por corriente, c) fuentes de corriente controladas por tensión y d) fuentes de corriente controladas por corriente (VERIFICAR, NO COINCIDE EN LA FIGURA DE LA PÁGINA SIGUIENTE).

35

II. Fundamentos teóricos

La acción de control puede ser ejercida por más de una variable y las funciones de control pueden ser complejas. Aunque la monopuerta puede especificarse arbitrariamente, su implementación como elemento de circuito puede no ser posible en tanto que no esté contenida en las librerías del software elegido. La teoría de circuitos permite, mediante circuitos auxiliares, resolver prácticamente todos los casos de diseño de la red eléctrica que se necesiten para cualquier tipo complejo de fuente controlada. La representación simbólica de estas fuentes se muestra en la Figura 2.11c. a) Fuente de tensión controlada por tensión

X

J

b) Fuente de tensión controlada por corriente

X( t )

J (t)

X( t )

X ( t ) = x 0 sen ( t )

J  t  pulso

J ( t ) = j0 e x p( t ) d) Fuente de corriente controlada por dos tensiones

c) Fuente de corriente controlada por tensión

X = FX ( X C )

X = FJ ( J C )

J = FX ( X C1 , X C2 )

{XC}

{J C }

{ X C1 , X C 2 }

J = FX ( X C )

J = FJ ( J C )

{ JC }

{XC}

Figura 2.11 Representación simbólica de monopuertas activas: 2.11a Fuentes constantes; 2.11b Fuentes dependientes del tiempo; Fig. 2.11c Fuentes controladas por una variable; 2.11d Ejemplo de fuente controlada por varias variables.

El potencial de estas monopuertas activas para establecer los modelos en red de sistemas fuertemente no lineales es inmenso ya que su uso permite imponer a la monopuerta el valor de una variable (en función de variables de otras monopuertas) sin influir en la otra variable, cuyo valor, se ajusta a la topología y geometría del modelo en red. En términos de componentes eléctricos el software elegido en esta memoria para la simulación, PSpice [1994], es capaz de reconocer un largo catálogo de componentes eléctricos: resistencias (R), condensadores (C), fuentes constantes de tensión y corriente, fuentes de tensión y corriente dependientes del tiempo, fuentes controladas de tensión y corriente y fuentes controladas no lineales de tensión y corriente. 36

II. Fundamentos teóricos

II.3.2 El MESIR como método numérico Desde el punto de vista de cálculo, el MESIR es, efectivamente, un método numérico en tanto que los procedimientos usados para obtener una solución aproximada de los circuitos eléctricos (que, a su vez, son los modelos en red de los problemas estudiados) son estrictamente numéricos. Se trata de procedimientos muy sofisticados y perfeccionados fruto de una continuada investigación en el campo eléctrico. Digamos que el MESIR es, en esencia, un usuario de estas técnicas ya que su aportación se limita a diseñar unos circuitos cuyas ecuaciones diferenciales en diferencias finitas son equivalentes a las correspondientes, discretizadas en el espacio, del modelo matemático del problema. Esta tarea, sin embargo, no es trivial y depende de la complejidad del modelo. Son numerosos los diferentes tipos de problemas que el MESIR ha abordado con éxito y que han sido publicados en la literatura científica. Entre estos cabe destacar: Problemas de difusión a través de membranas (Horno y col. [1990]), ii) ídem. de procesos electroquímicos (González-Fernández y col. [1995]), iii) ídem. de determinación de propiedades superficiales de coloides (Horno y col. [1993], LópezGarcía y col. [1996]), iv) ídem. de transferencia de calor (Alhama [1999], Alarcón [2001], Alarcón y col. [2002a, 2002b y Zueco y Alhama [2007]), v) ídem. de flujo de fluidos (Soto y col. [2007a, 2007b]), vi) ídem. mecánicos (cadenas cinemáticas, problemas tribológícos, modos de vibración en vigas…, Moreno y col. [2004,2007]), vii) ídem. de ecología (López Sánchez y col. [2005], viii) problemas inversos (Zueco y Alhama [2005,2006] y Zueco y col. [2006a,2006b]), ix) procesos caóticos (Alcover y col. [2006]), etc. A diferencia de otros métodos numéricos clásicos, como los de formulaciones en diferencias finitas explícita, implícita e híbrida, métodos de elementos finitos, métodos variacionales, métodos iterativos específicos y de autovalores (Crank y Nicolson [1947], Crandall [1956], Wood [1977], Orivuori [1979], Seeniraj y Arunachalam [1979], Lin [1979], Vujanovic, [1994], Fried y Metzler [1978] y Patankar [1978], Varoglu y Liam Finn [1980], Chung [1981] y Palmieri y Rathjem [1978]), los cuales requieren una preparación matemática considerable y modificaciones sustanciales en el software al cambiar cualquier dependencia paramétrica o condición de contorno, el MESIR es una técnica que aprovecha los potentes algoritmos integrados en los modernos software de resolución de circuitos los cuales, quizás, son los más sofisticados y capaces ya que tienen que afrontar señales fuertemente no lineales (pulsos) y frecuencias muy altas. Además, conserva en cierto modo la visualización del proceso de transporte que siempre queda oscurecida entre todo el aparato matemático en los otros métodos. El punto de partida es siempre el modelo matemático del proceso, esto es, un conjunto de ecuaciones en derivadas parciales (EDP) espacio-temporales; la discretización de la variable espacial permite establece el modelo en red o red eléctrica equivalente. Esta es la única manipulación directa que se hace de las ecuaciones. 37

II. Fundamentos teóricos

El modelo en red es, pues, el formato que se da al modelo matemático para que pueda ser utilizado como entrada (fichero) en un programa de resolución de circuitos eléctricos tal como PSpice (Nagel [1975], PSpice [1994], Vladimirescu [1994]). Este software es el que resuelve las ecuaciones de la red y proporciona la solución numérica del modelo matemático. A continuación exponemos las diferencias de estrategias más notables al compararlo con otros métodos numéricos. Cuando en una ecuación en derivadas parciales se hace una doble reticulación, espacial y temporal, se reemplazan de hecho las derivadas parciales por aproximaciones algebraicas, lo que conduce a un conjunto de ecuaciones algebraicas que aproximan las EDP. Para la solución numérica de éstas se utiliza un software adecuado de matemáticas. Este procedimiento es la base de los bien conocidos métodos numéricos de diferencias finitas, elementos finitos y volúmenes finitos para la solución de las EDP. Como ya se ha comentado, la elaboración del modelo en red pasa por la reticulación espacial, pero no temporal. Se parte, pues, de un sistema de ecuaciones en derivadas parciales cuya reticulación espacial las convierte en ecuaciones diferenciales ordinarias en el tiempo, que son las del circuito correspondiente a una celda elemental. La diferencia esencial es, pues, que en los métodos numéricos convencionales se realiza una reticulación simultánea de las dos variables independientes, espacio y tiempo, mientras que en el MESIR la reticulación es sucesiva; 1ª etapa, una reticulación espacial de la que se obtiene el modelo en red y 2ª etapa, una reticulación temporal, realizada por el propio software en el proceso de simulación. Alhama [1999] ha demostrado que la precisión de los resultados de la simulación o error respecto de la solución exacta, en problemas lineales, depende del tamaño de la reticulación, pero son suficientes reticulaciones del orden de 40 a 60 elementos de volumen para reducir estos errores a valores por debajo del 0.5-0.2%. Cuando se trata de problemas fuertemente no lineales, por ejemplo problemas de cambio de fase con frontera móvil, basta duplicar el tamaño de la retícula para obtener soluciones con errores del mismo orden.

II.4

MODELOS EN RED

II.4.1 La analogía termoeléctrica clásica Una amplia y detallada discusión sobre este punto, con numerosas referencias bibliográficas, se puede ver en Alhama [1999]. Es frecuente el uso de analogías eléctricas de procesos simples de transmisión de calor por un interés meramente académico; sencillamente porque las ecuaciones algebraicas del proceso de transporte (no diferenciales) aplicadas a medios finitos son exactamente iguales a las que relacionan la intensidad y la tensión (ley de Ohm) en los componentes pasivos de los circuitos eléctricos.

38

II. Fundamentos teóricos

Textos recientes como Grigull y Sandner [1984], Chapman [1984], Kakaç y Yener [1985], Lienhard [1987], Siegel y Howell [1992] y Mills [1995], no olvidan este aspecto tanto en procesos de conducción como de radiación; Holman [1981], presenta una discusión muy didáctica de la aplicación de la analogía eléctrica a problemas más complicados de transmisión de calor, en particular en el campo de radiación, con numerosos ejemplos. Otros textos clásicos y modernos no hacen referencia alguna a la analogía termo-eléctrica (Carslaw y Jaeger [1959]; Özisik [1993] y Gebhart [1993]) o la mencionan muy de pasada (Bennett y Myers [1982]; White [1988], Bejan [1993], Taine y Petit [1993]). Textos modernos, específicos de tratamiento numérico, (Patankar [1980], Shih [1984], Ozisik [1989] y Gebhart [1993]) no hacen, obviamente, referencia alguna a analogías eléctricas ni mencionan el “network method” de Oppenheim [1956]. En todo caso las referencias en estos textos siempre aparecen en capítulos que tratan problemas lineales y con condiciones de contorno simples, primera y segundan clase. Ningún texto ni artículo, con excepción de los publicados por miembros del grupo de simulación por redes, menciona la analogía termo-eléctrica como método de cálculo numérico. Un ejemplo importante es el uso de la analogía eléctrica empleado por Chapman [1984]. Dedica un importante apartado a este tema pero siempre dentro de los procesos de transmisión de calor lineales, cuya solución sólo requiere resistencias térmicas, condensadores para el almacenamiento y fuentes constantes para las condiciones de contorno. Hemos de mencionar un caso especial de analogías a las que se hace referencia como “lumpedcapacity systems”, “lumped-thermal-capacity problems” o “lumped system formulation” (Lienhard [1987], Özisik [1993], Mills [1995], Thomas, 1992...). Se trata de problemas de interés exclusivamente académico. Las aportaciones reales de interés científico de la analogía termo-eléctrica se reducen a las derivadas de la técnica “resistance-network model”, que mediante el uso exclusivo de resistencias se llega a la construcción de complicados circuitos en los que manipulando el intervalo de tiempo permiten simular problemas no lineales, que incluyen procesos de cambio de fase; Liebmann [1954a, 1954b, 1956a y 1956b] estudia la elevación térmica en cuerpos de turbina durante el transitorio, y Bonilla y Strupezewsky [1965] que estudian el comportamiento térmico de carcasas de reactores ante pruebas de fuego. A ellas hay que añadir las derivadas del MESIR ya mencionadas. El paréntesis de casi tres décadas que separa estas publicaciones es debido a las dificultades inmensas en la realización práctica de los circuitos para obtener pruebas fiables; en la década actual los programas de resolución de circuitos disponibles en el mercado han salvado este obstáculo. Como dice Özisik [1993] “El método analógico de redes eléctricas, usado a menudo en los primeros años (Kreith y Romie [1955], Otis [1956], Liebmann [1956a y 1956b] y Baxter

39

II. Fundamentos teóricos

[1961]) ha sido desplazado por la aplicación de métodos numéricos puros a causa de la gran potencialidad de los ordenadores digitales”, o Lienhard [1987], “... los computadores digitales actuales hacen posible la solución de complicados problemas por aplicación directa de métodos numéricos”. De todo ello se puede deducir que el “método de simulación por redes” es algo sustancialmente diferente a la analogía termo-eléctrica clásica esencialmente por su capacidad de abordar cualquier tipo de problemas lineales o no, acoplados o no, y con condiciones de contorno arbitrarias. II.4.2 Modelos para medios homogéneos y para las condiciones de contorno El modelo en red de una celda elemental correspondiente a un medio 1-D, homogéneo, de espesor ∆ x y superficie unidad (normal a la propagación del calor) es el indicado en la Figura 2.12, González y Alhama [2002]. La celda se ha organizado simétricamente disponiendo condensador de almacenamiento en el punto central de la misma. La pared completa estará constituida por N celdas iguales conectadas en serie. Los valores de los componentes resistencia y condensador son

R izq = R der =

Δx1 Δx = 2kS 2N1kS

C = ρ ce Δx 1 =

(2.22)

ρ ce l0 N1

donde lo es la longitud total del dominio. El modelo asociado a medios multicapa formados por capas homogéneas de diferente naturaleza es una conexión en serie del modelo de capa simple. Cuando el contacto entre capas no es perfecto hay que añadir una resistencia entre capas. El número total de celdas es N 1 + N2 + … Nn. Para resistencias de contacto constantes, un medio de n capas de espesores l1, l2, … ln, superficie unidad y resistencia de contacto entre capas despreciable, el modelo en red se muestra en la Figura 2.13.

Ti - ∆ ji - ∆

Ti Rizq

T1+∆

j i, c

C

Rder

Figura 2.12 Modelo en red de la celda elemental

40

ji +∆

II. Fundamentos teóricos

Tf

Celda 1

Celda 2

x = ∆x

x =0

Celda 3

x = 2∆x

Celda N

x = N∆x

Figura 2.13 Modelo en red de un medio multicapas

Para geometrías cilíndrica y esférica 1-D, con variable espacial r, los modelos tienen la misma configuración pero el valor de las resistencias y condensador son, Figura 2.13, Geometría cilíndrica:

R izq =

R der =

Δri   Δr   2k  2π ri + i   4    (2.23)

Δri   Δr   2k  2π ri + 3 i   4   

C = ρ ⋅ ce ⋅ 2π ri Δri

Geometría esférica:

R izq =

R der =

Δri Δr   2k ⋅ 4π ri + i  4  

2

Δri Δr   2k ⋅ 4π ri + 3 i  4  

(2.24)

2

C = ρ ce ⋅ 4π ri 2Δri

41

II. Fundamentos teóricos

42

Figura 2.14 Nomenclatura de la celda en las geometrías cilíndrica y esférica

La celda rectangular 2-D se muestra en la Figura 2.15. Las resistencias horizontales y verticales pueden ser diferentes si se adoptan valores ∆ x≠ ∆ y. Los valores de las resistencias y condensador son

R izq = R der =

Δx i 2k ⋅ Δyi

R inf = R sup =

Δyi 2k ⋅ Δx i

(2.25)

C = ρ ⋅ ce Δx ⋅ Δy

Para la geometría cilíndrica 2-D el modelo es el mismo. Los valores de las resistencias y condensador son

R izq =

R der =

Δri  Δr    2k  2ππ  ri + i   4    Δri  Δr    2k  2ππ  ri + 3 i   4   

(2.26)

C = ρ ⋅ ce ⋅ 2 π ri Δri H

Figura 2.15 Modelo en red de la celda elemental 2-D. Geometría rectangular

En relación con los modelos de aletas simples, 1-D, los valores de las resistencias y condensador son R izq = R der =

Δz 2ks t

(2.27)

C = ρ ⋅ ce ⋅ Δz ⋅ s t

donde ∆ z es el espesor de la celda elemental y S t la sección transversal de la aleta. La condición de convección modifica el modelo. Esta condición se implementa introduciendo un generador controlado por voltaje en el centro de la celda, Figura 2.16. La corriente de este generador, especificada mediante programación es el flujo de convección dado por la ley de newton, jcon = h ( Ti − Tref ) , donde Ti es la temperatura en el centro de la celda, Tref la del fluido exterior y h el coeficiente de convección. La polaridad de estos generadores debe ser la adecuada para que el flujo de corriente circule en el sentido de mayor a menor temperatura. Si se trata de la condición de radiación, el modelo no cambia pero la corriente viene definida por la ley de Stefan de la radiación,

(

)

4 4 jrad = ε σSi T1 − Tref , donde ε es la emisividad de la superficie, σ la constante de Stefan

y Si la superficie exterior del elemento de volumen.

Figura 2.16 Modelo en red de la celda elemental de una aleta simple

La existencia de condición simultánea de convección y radiación puede implementarse separadamente (o conjuntamente) por medio de dos generadores unidos al centro de la celda. En aletas anulares 1-D, los valores de las resistencias y condensador son R izq = R der =

Δz 2kπR 2

(2.28)

donde ∆ z es el espesor de la celda elemental y R el radio de la sección transversal de la aleta. El modelo en red es el mismo que el mostrado en la figura anterior. En relación con las aletas 2-D, conjunto aleta-pared rectangular, puede usarse el modelo 2-D de la Figura 2.15 para las celdas interiores en combinación con el de la Figura 2.17 para las celdas del contorno, con 1 ó 2 generadores para una o las dos condiciones de contorno simultáneas (convección y radiación).

Figura 2.17 Modelo en red de una celda elemental de contorno en aletas 2-D

En relación con las condiciones de contorno (aparte de lo indicado anteriormente para la convección y radiación, cuyos modelos en el caso de aletas se han integrado en las celdas de contorno), su implementación se muestra en la Figura 2.18. La condición isoterma se simula mediante una fuente de voltaje de valor constante; otras fuentes dependientes del tiempo tales como voltajes sinusoidales, en rampa, continua a trazos, rectangular, etc., pueden implementarse fácilmente pues estos componentes eléctricos están integrados en las librerías de los programas de simulación), Figura 2.18a y 2.18b. Para el flujo constante o dependiente del tiempo se usan fuentes de corriente, Figura 2.18c y 2.18d. La condición adiabática se implementa con una resistencia de valor elevado, teóricamente de valor infinito, Figura 2.18e. Las condiciones de convección radiación se implementan mediante fuentes controladas de corriente, Figura 2.18f y 2.18g.

Figura 2.18 Componentes para implementar las condiciones de contorno. Fig. 2.18a Isoterma; Fig. 2.18b De temperatura dependiente del tiempo; Fig. 2.18c De flujo de calor constante; Fig. 2.18d De flujo de calor dependiente del tiempo; Fig. 2.18e adiabática; Fig. 2.18f De convección o de radiación.

II.5

EL SOFTWARE PSpice-OrCAD

II.5.1 Introducción y aplicaciones Una vez obtenido el modelo en red se procede a su análisis mediante su simulación. Para ello hemos buscado un software adecuado para la solución de circuitos eléctricos tal como PSpice u OrCAD [1994]. El segundo es, en realidad, la versión actual del primero. PSpice ha sido utilizado por otros autores para resolver problemas de otras disciplinas. Baker y Shortt [1990] estudia el comportamiento de componentes integrados en diferentes rangos de temperatura, Bello [1991] lo aplica a la resolución de problemas mecánicos, Herbert [1992] lo aplica a la resolución de ecuaciones diferenciales ordinarias, Hamill [1993], a problemas estadísticos y relacionados con el caos, etc. En el proceso de simulación el circuito se presenta al ordenador como un conjunto de ecuaciones matemáticas y éste, mediante procedimientos de análisis numérico, proporciona toda la información solicitada por el investigador para cada tipo de análisis. De esta forma se obtienen los datos correspondientes a medidas típicas de laboratorio con un margen de error despreciable y sin afectar al circuito; más aún, pueden alterarse las condiciones iniciales, de contorno, y las características térmicas del medio con sencillos cambios en el programa, y el análisis puede aportar datos sobre el comportamiento del circuito más allá de los límites que virtualmente se pueden obtener con medidas reales. La simulación está estructurada en cinco subprogramas principales, que interaccionan entre ellos a través de una estructura de datos que es almacenada en un área común del programa. Estos subprogramas son: entrada, organización, análisis, salida y utilidades, Figura 2.19. El subprograma de entrada lee el archivo de entrada, construye una estructura de datos y chequea el circuito. El de organización, una vez que el programa se ha ejecutado con éxito, construye las estructuras adicionales de datos que serán requeridas en el programa de análisis, parte esencial de la simulación. El subprograma de salida genera y organiza, en la memoria central o en discos, los resultados solicitados por el usuario en forma tabular o gráfica. Las utilidades son aspectos secundarios no relacionados directamente con la simulación; éstas permiten, por ejemplo, almacenar componentes o partes de modelos para ser compartidos por otros usuarios. El subprograma análisis es la parte importante del programa de simulación. Ejecuta los análisis del circuito requeridos, de acuerdo con las indicaciones del archivo de entrada; la información resultante se almacena en la memoria central o en discos para su posterior procesamiento en los archivos de salida. Mientras que la facilidad de uso del programa reside en los subprogramas de entrada y salida, el programa de análisis, que contiene algoritmos más complejos y consume la fracción mayor del tiempo de computación, determina la eficiencia de la simulación.

CONTROL

ENTRADA

ORGANIZACI ÓN

ANÁLISIS

SALIDA

UTILIDADES

ESTACIONARIO

TRANSITORIO

DE PEQUEÑA SEÑAL

Figura 2.19 Diagrama bloques del programa de simulación de circuitos PSpice

En el proceso de simulación, se obtiene la solución numérica de la representación matemática del modelo en red. Esta contiene i) las ecuaciones matemáticas de los diferentes tipos de monopuertas, ii) las ecuaciones correspondientes a las restricciones impuestas por las leyes de Kirchhoff, propias de la teoría de circuitos, que han de satisfacerse entre las ramas y nudos del circuito, y iii) la información particular sobre la interconexión de los diferentes componentes eléctricos de cada modelo. Toda esta información compone un extenso sistema de ecuaciones algebraico-diferenciales. El conjunto de tareas que componen el proceso de simulación puede ser agrupado en los siguientes tópicos (o algoritmos de computación): i) formulación de las ecuaciones, ii) solución de ecuaciones lineales, iii) solución de ecuaciones no lineales y iv) integración numérica. PSpice es miembro de la familia de programas de simulación de circuitos SPICE2 (Nagel, [1975]); mucho más potente y rápido que sus predecesores, fue desarrollado en la Universidad de California en los años setenta y utiliza algoritmos numéricos más refinados con formatos de entrada-salida idénticos. En el análisis de continua PSpice determina el punto de trabajo, es decir, los valores de polarización de sus componentes en ausencia de excitaciones alternas. Para este cálculo se elimina la acción de los condensadores y bobinas, los primeros quedan como circuitos abiertos y las bobinas se cortocircuitan. Para el análisis transitorio PSpice parte del intervalo de tiempo (0,t) solicitado, que puede ser menor o mayor que la duración del transitorio, y facilita los datos en forma de listado o mediante gráficas. Si los resultados se quieren en forma tabular el usuario debe indicar el instante inicial, el final, el paso temporal y el número de variables listadas; si se solicitan en forma gráfica, una sentencia de programa permite organizarlos y almacenarlos para ser utilizados con ese propósito en cada momento.

La formulación de las ecuaciones del circuito se realiza mediante el método conocido como Análisis Nodal Modificado La solución transitoria se determina computacionalmente extrayendo del intervalo temporal un conjunto discreto de instantes (0, t 1, t2, ..., T). En cada uno de ellos, empezando por 0, el tiempo se incrementa una pequeña porción o paso, δ t, y, mediante métodos de integración (algoritmo implícito de Backward-Euler) y procesos de iteración hasta conseguir la convergencia, se resuelven las ecuaciones algebraicas equivalentes de las monopuertas que contienen derivadas temporales; cada iteración requiere de la linealización de las ecuaciones del modelo y de su solución; el método de linealización es el de Newton-Raphson que utiliza una serie de Taylor truncada después del término de primer orden. Para la solución del sistema matricial de ecuaciones lineales se utiliza el método de factorización LU que elimina directamente las incógnitas (este método descompone la matriz de coeficientes en producto de matrices triangulares, “lower and upper, LU”, cuyas inversas no precisan ser calculadas, lo que redunda en un menor esfuerzo computacional). Para minimizar el esfuerzo de cálculo, las ecuaciones se reordenan usando el método de Markowitz. Los métodos de integración implantados en PSpice incorporan un proceso de variación dinámica del paso del tiempo de integración para mantener una razonable exactitud en la solución y un tiempo mínimo de computación. PSpice utiliza unos métodos de integración polinomiales basados en el análisis de error de truncamiento local y en la estabilidad (propiedades contrapuestas). Debido a que ciertos circuitos (que contienen constantes de tiempo de valores muy diferentes) pueden dar lugar a un sistema de ecuaciones “stiff”, es conveniente que el algoritmo de integración sea “stiff-estable”. Para conseguir este objetivo se utilizan métodos de integración trapezoidal y Gear de orden variable de dos a seis. Tras conseguir la convergencia, la solución se almacena y se reinicia el proceso para el instante siguiente. El paso δ t es, en consecuencia, variable de unos tramos del intervalo a otros; el programa los ajusta en función de la precisión exigida a los resultados de manera que el tiempo de computación sea el mínimo. Los datos de simulación correspondientes a tiempos fuera del conjunto discreto de instantes 0, t1, t2, ..., T se obtienen por interpolación. La Figura 2.20 representa un diagrama de flujo que ilustra los cuatro algoritmos de computación que tienen lugar en la simulación de un proceso transitorio (para simplificar se ha supuesto un δ t constante).

t=0

t=t+h

Integración numérica

Linealización

Formulación de las ecuaciones

Soluc ión de las ecuaciones lineales No ¿Converge? Sí Almacenar la solución No

¿Fin de l transitorio? Sí Archivo de salida

Figura 2.20 Operaciones en el análisis de circuitos

Los algoritmos utilizados en PSpice, que se documentan en la tesis de Nagel, son el resultado de implementaciones, modificaciones y comparaciones cuidadosas de los métodos numéricos existentes en el contexto especial de la simulación de circuitos. El objeto de la tesis es seleccionar los métodos de simulación de circuitos más exactos y eficaces, con una mínima interacción por parte del usuario.

II.5.2 Simulación. Presentación de resultados El software PSpice se programa en su forma clásica por sentencias, elaborando archivos de texto, en un lenguaje relativamente simple (alternativamente es posible elaborar archivos por medio de la opción gráfica ‘schematics’ seleccionando directamente los elementos de circuito y conectándolos eléctricamente entre si en forma de esquela eléctrico). La sintaxis de entrada no requiere especiales disposiciones ordenadas de datos, su estilo puede catalogarse más bien como libre y dispone de una razonable fuente de datos que se adjudican por omisión a los componentes cuando éstos no se especifican en detalle. También realiza un buen número de chequeos para asegurar que el circuito ha sido introducido correctamente y el resto de las sentencias de programa están bien escritas, advirtiendo al programador de posibles errores mediante mensajes previos a la ejecución. En definitiva, un usuario principiante necesita especificar un número mínimo de parámetros y controles de simulación para extraer unos resultados de simulación aceptables. El programa, por fin, se estructura como un listado que contiene todos los componentes eléctricos del circuito (existe la posibilidad de organizar el programa mediante subcircuitos), resistencias, condensadores, fuentes, interruptores, etc., que se introducen uno por uno indicando el nombre, valor, nudos de conexión y otros parámetros característicos. El programa PSpice (como, en general, cualquier otro software de resolución de circuitos eléctricos) ofrece muchas posibilidades para el estudioso de los sistemas térmicos. A continuación se destacan algunas utilizadas en esta memoria: •

Permite conocer directamente el estacionario del sistema térmico (BIAS POINT), mediante el análisis en continua del circuito, antes referido. La opción TRANS proporciona el transitorio del proceso,

•

Con el análisis en alterna se obtiene de forma inmediata el análisis de respuesta en frecuencia del sistema térmico,

•

La aplicación Probe muestra de forma gráfica los resultados de la simulación con la máxima precisión que da el programa. Esta aplicación permite la representación de funciones resultado de operaciones entre variables de la simulación. Por ejemplo, las curvas de la admitancia o la impedancia del sistema (cociente entre corriente y tensión o viceversa) pueden ser directamente obtenidas de Probe,

•

El software admite la parametrización del modelo en red (sentencia PARAM), lo que constituye un modo ventajoso de utilizar la técnica de cambiar de valores los componentes del circuito para obtener soluciones de problemas similares,

•

Las sentencias PARAM y STEP combinadas obtienen la variación secuencial de la respuesta del sistema ante la variación de un parámetro, lo que es una herramienta muy útil para problemas sencillos de optimización (una o dos variables),

•

La aplicación Stimulus permite la confección de fuentes de tensión o corriente de prácticamente cualquier forma, que pueden representar cualquier estímulo térmico del sistema,

•

El programa admite la ejecución sucesiva de programas, técnica que permite arrancar indefinidamente el PSpice por otro programa y resolver el circuito para una amplia gama de valores de los componentes. En este caso el programa actúa como un solver, cuya misión es resolver las ecuaciones diferenciales del sistema, mientras que al otro programa se le confía la resolución de un problema más amplio.

En relación con la presentación de resultados, el programa permite acceder a los resultados de la simulación (temperaturas y flujos de calor en todo el medio) de dos formas: i) directamente usando el entorno gráfico de PSpice, muy intuitivo y potente, o accediendo a los archivos de salida de datos los cuales muestran los resultados en forma tabulada; en general estos resultados vienen dados usando como variable independiente en tiempo por lo que son muy útiles en problemas transitorios pero no tanto en problemas estacionarios, y ii) mediante representaciones gráficas alternativas, bien en el entorno directo de PROCCA-09 bien por medio de MATLAB [1997], usando opciones del programa. Estas representaciones alternativas, muy útiles en problemas estacionarios, se han elaborado para dar mayor información pues consisten en tipos de representación no proporcionados por el entorno gráfico PSpice (por ejemplo mapas de temperatura en función de la posición). Además pueden obtenerse representaciones animadas de isotermas en problemas transitorios tanto en el entorno gráfico del propio programa como en MATLAB. Todo ello con acceso directo desde las propias ventanas del PROCCA-09 mediante comandos de botón que dan acceso al usuario a todas las opciones de representación incorporadas al programa. Finalmente, es posible extraer los datos en forma tabulada para ser manipulados mediante hoja de cálculo.

II.6

El PROGRAMA C#

C# (pronunciado en inglés “C Sharp” y en español “C Almohadilla”) es el nuevo lenguaje sencillo e intuitivo diseñado por Microsoft específicamente para su plataforma .NET (Archer [2004], Gunnerson [2002], Albahari y col. [2001] y Robinson y col. [2001]). Se trata de un lenguaje orientado a objetos sencillo, moderno, amigable, intuitivo y fácilmente legible que ha sido diseñado por Microsoft con el ambicioso objetivo de recoger las mejores características de muchos otros lenguajes, fundamentalmente Visual Basic, Java y C++, y combinarlas en uno solo en el que se unen la alta productividad y facilidad de aprendizaje de Visual Basic con la potencia de C++.

C# incorpora muchos elementos de los que carecen otros lenguajes, tales como sistema de tipos homogéneo, propiedades, indexadores, tablas multidimensionales, operadores redefinibles, etc.) y tiene una velocidad de ejecución muy competitiva. Las principales características de este lenguaje son: -

Dispone de todas las características propias de cualquier lenguaje orientado a objetos: encapsulación, herencia y polimorfismo,

-

Permite definir estructuras, que son clases un tanto especiales: sus objetos se almacenan en pila de acceso rápido,

-

Es un lenguaje fuertemente tipado, lo que significa se controla que todas las conversiones entre tipos se realicen de forma compatible, lo que asegura que nunca se acceda fuera del espacio de memoria ocupado por un objeto. Así se evitan frecuentes errores de programación y se consigue que los programas no puedan poner en peligro la integridad de otras aplicaciones,

-

Tiene a su disposición un recolector de basura que libera al programador de la tarea de tener que eliminar las referencias a objetos que dejen de ser útiles, encargándose de ello éste y evitándose así que se agote la memoria porque al programador olvide liberar objetos inútiles o que se produzcan errores porque el programador libere áreas de memoria ya liberadas y reasignadas,

-

Incluye soporte nativo para eventos y delegados. Los delegados son similares a los punteros de funciones de otros lenguajes como C++ aunque más cercanos a la orientación a objetos, y los eventos son mecanismos mediante los cuales los objetos pueden notificar de la ocurrencia de sucesos. Los eventos suelen usarse en combinación con los delegados para el diseño de interfaces gráficas de usuario, con lo que se proporciona al programador un mecanismo cómodo para escribir códigos de respuesta a los diferentes eventos que puedan surgir a lo largo de la ejecución de la aplicación. (pulsación de un botón, modificación de un texto, etc.),

-

Incorpora propiedades, que son un mecanismo que permite el acceso controlado a miembros de una clase tal y como si de campos públicos se tratasen,

-

Permite la definición del significado de los operadores básicos del lenguaje (+, -, *, &, =, etc.) para nuestros propios tipos de datos, lo que facilita enormemente tanto la legibilidad de las aplicaciones como el esfuerzo necesario para escribirlas. Es más, se puede incluso definir el significado del operador [] en cualquier clase, lo que permite acceder a sus objetos tal y como si fuesen tablas. A la definición de este último operador se le denomina indizador, y es especialmente útil a la hora de escribir o trabajar con colecciones de objetos.

-

Admite unos elementos llamados atributos que no son miembros de las clases sino información sobre éstas que podemos incluir en su declaración. Por ejemplo, indican si un miembro de una clase ha de aparecer en la ventana de propiedades de Visual Studio.NET, cuáles son los valores admitidos para cada miembro en ésta, etc.

-

Visual Studio.Net permite diseñar la interfaz de la aplicación de manera visual, sin más que arrastrar con el ratón los elementos que necesitemos (botones, lista de selección, etc.) sobre las posiciones adecuadas en la ventana de nuestra aplicación. También incluye otras facilidades para el desarrollo, como una ventana de propiedades desde la que se puede modificar los valores de las propiedades de cada objeto sin tener que escribir código, un depurador de código gráfico, un editor de códigos inteligente que puede detectar nuestros errores de sintaxis instantáneamente, etc.

III. Programa PROCCA-09

Capítulo III

El programa PROCCA-09

III.1 Introducción III.2 Estructura del programa PROCCA-09 III.3 Creación de archivos de modelos III.3.1 Presentación del programa III.3.1.1 Descripción del módulo CONCBA III.3.1.2 Descripción del módulo CONCAL III.4 Criterios para la numeración de celdas, nodos y elementos del modelo III.5 Estructura de los archivos de texto de modelos III.6 Pantallas de presentación de resultados III.7 Ejemplos de archivos de modelo III.7.1 Placa rectangular con oquedad III.7.2 Cilindro hueco de dos capas III.7.3 Aleta rectangular 1-D III.1 INTRODUCCIÓN 57

III. Programa PROCCA-09

Este capítulo contiene una descripción detallada del manejo del programa PROCCA-09, desde la generación de modelos de archivo hasta su ejecución y presentación de resultados en sus diferentes opciones. Se describen, sucesivamente, las pantallas de introducción de datos con explicación de los detalles que contienen, la estructura de los archivos de texto de los modelos, su manipulación, la ejecución de estos en PSpice y la consiguiente manipulación de resultados de acuerdo con los diferentes entornos de salida. Para aclarar en lo posible el uso del programa se describen al final de este capítulo tres ejemplos que cubren una gran parte del espectro posible de posibilidades de PROCCA-09. III.2 ESTRUCTURA DEL PROGRAMA PROCCA-09 PROCCA-09 es, en principio, uno más entre los muchos programas comerciales educativos aplicados a ciencias e ingeniería. Sin embargo, su diseño tiene ciertas peculiaridades que creemos le diferencia y distingue de otros programas básicamente orientados al cálculo, con un diseño numérico propio, y que funcionan a modo de cajas negras de contenido, en general, inaccesible a usuario. El objetivo de PROCCA-09 es tanto científico como didáctico, en consecuencia, no solamente realiza el cálculo numérico necesario para la simulación de los problemas de conducción térmica para los que ha sido desarrollado sino que permite aprender, simultáneamente, los contenidos básicos de esta ciencia en lo referente tanto a problemas de conducción de calor en medios multicapa de diferentes geometrías (módulo CONCBA) como a problemas de análisis y diseño de aletas (módulo CONCAL). El programa objeto de esta tesis, que hace uso de la analogía o equivalencia entre el transporte eléctrico y la conducción térmica, se presenta al usuario mediante un entorno de comunicación ameno, tipo “window”, que dirige paso a paso las acciones y opciones posibles del usuario, tales como selección del tipo de problema, entrada de datos, creación y manipulación de archivos de modelos, simulación, presentación de resultados… Los archivos de modelos en red se ejecutan en PSpice [1994] y los resultados de simulación se ofrecen directamente en el entorno, gráfico o tabulado, de salida de PSpice o bien, mediante manipulaciones adecuadas incorporadas en PROCCA-09, en forma gráfica bien directamente en el entorno gráfico del propio programa como (en mayor detalle) en el entorno gráfico del software MATLAB [1997] por medio de rutinas auxiliares incorporadas el programa. PROCCA-09 permite también presentar soluciones animadas de las isotermas en problemas transitorios. La Figura 3.1 presenta un esquema del funcionamiento básico del programa. Su puesta en marcha permite seleccionar el módulo de trabajo (CONCBA o CONCAL),..0. en donde se incardina el problema objeto de simulación o solución, y se procede directamente a la entrada de datos: geometría de la reticulación, características térmicas, condiciones de contorno, etc. Una vez completada la especificación se puede crear un archivo básico del modelo que permite su 58

Figura 3.1 Diagrama de funcionamiento de PROCCA-09

manipulación directa y su modificación. La introducción de datos complementarios relativos al tipo de simulación, tales como la precisión en los cálculos, número de dígitos, tiempo de simulación, opciones de presentación y otros, se introduce paralelamente o al final de la creación del archivo del modelo. La ejecución o simulación y consiguiente solución del modelo da acceso al entorno de salida gráfico de PSpice el cual permite representar simultáneamente las variables de salida más comunes, a saber, los flujos de calor y las temperaturas (corrientes y tensiones) en los distintos elementos y nudos del medio elegidos por el usuario. La asignación de nombres a los elementos del modelo en red así como la asignación de nodos, tanto en problemas 1-D como 2-D sigue una regla sencilla e intuitiva independiente del tipo de geometría, permitiendo al usuario localizar inmediatamente el elemento, sección o punto del modelo del que se desea obtener información acerca del valor de las variables dependientes, flujos de calor y temperaturas. El programa PROCCA-09 incorpora un gráfico directo en el que se muestra la disposición de las celdas o volúmenes de control mostrando la leyenda de los nudos centrales de cada celda, que sirve también para identificar todos los componentes de la celda (resistencias y condensadores). En este gráfico también se muestran las celdas en las que se incorporan condiciones de contorno particulares. Por último, un archivo de ayuda accesible desde cualquier paso del programa da información al usuario de cómo resolver e interpretar las dificultades que surgen con el uso del programa. III.3 CREACIÓN DE ARCHIVOS DE MODELOS III.3.1 Presentación del programa La ejecución del programa PROCCA-09 cuyo icono aparece en la Figura 3.2, da entrada al mismo por medio de la pantalla de la Figura 3.3, “Diseñador de Modelos”, con la que se inicia el diseño del modelo buscado o se elije un modelo anteriormente diseñado y archivado en el programa. Los botones de acceso al usuario son tres: “Archivo” Vista” y “Ayuda”. Con el primero, Figura 3.4, se opta bien por iniciar el diseño de un nuevo modelo a través de los módulos CONCBA y CONCAL (opciones “Nuevo modelo CONCBA, conducción de calor básica” y “Nuevo modelo CONCAL, conducción de calor en aletas”), bien por importar modelos ya existentes en otras carpetas (opción “Cargar”), para lo cual se hace clic directamente en el directorio de archivo de los mismos que se presenta en la pantalla, Figura 3.5. Las opciones “Nuevo modelo CONCBA” y “Nuevo modelo CONCAN” presentan, respectivamente, las pantallas mostradas en la Figura 3.6 a) y b).

Figura 3.2 Icono del programa PROCCA-09

Figura 3.3 Pantalla de inicio del programa PROCCA-09

Figura 3.4 Opciones de inicio de la barra de herramientas de PROCCA

Figura 3.5 Pantalla de búsqueda de modelos ya diseñados

Figura 3.6 Opciones de selección del problema. a) Módulo CONCBA

Figura 3.6 Opciones de selección del problema. b) Módulo CONCAL

III.3.1.1 Descripción del módulo CONCBA En la pantalla correspondiente al módulo CONCBA se selecciona el “Tipo de modelo” con arreglo a la geometría del mismo, “Rectangular”, “Cilíndrica” y “Esférica”, Figura 3.7. La selección de cada una de estas opciones da paso a las pantallas de la Figura 3.7 a), b) y c), respectivamente. En estas pantallas se introduce el número de celdas horizontales y verticales adoptado para la simulación. Además, para geometrías cilíndrica y esférica se introduce también el radio interior (0 para medios no huecos).

Figura 3.7 Opciones del “tipo de modelo” (Módulo CONCBA). a) geometría rectangular

Figura 3.7 Opciones del “tipo de modelo” (Módulo CONCBA). b) geometría cilíndrica

Figura 3.7 Opciones del “tipo de modelo” (Módulo CONCBA). c) geometría esférica

Una vez introducidos los datos de esta pantalla se pulsa el botón “Siguiente” para dar paso a la pantalla de la Figura 3.8, que corresponde a un modelo rectangular, que permite introducir los datos de información de cada capa. Simultáneamente en la pantalla principal “Diseñador de modelos” aparece una retícula con el número de celdas seleccionado, 15(horizontal)× 10(vertical) en la Figura 3.9.

Figura 3.8 Pantalla de introducción de datos por capas en geometría rectangular

Figura 3.9 Retícula de celdas del medio

La introducción de datos se realiza capa a capa definiendo por medio de la pantalla de la Figura 3.10 a), 3.10 b)…, el nombre de la capa, las celdas que ocupa, el tamaño de la capa (anchura y altura de la misma), las propiedades térmicas del material (densidad, conductividad y calor específico de la capa) y la temperatura inicial. Todas las unidades se introducen en el SI. La información ancho y alto de la celda de la capa actual la calcula el programa y la presenta en la pantalla para información al usuario. La elección de este método de introducción de datos tiene la ventaja de poder asociar a cada región o capa tanto un tamaño diferente de las celdas que la definen como una condición inicial también diferente, lo que permite ampliar el campo de aplicación del programa al poder reticular más finamente aquellas regiones en las que los procesos locales transitorios son más cambiantes, amén de otras ventajas de aspecto educativo (influencia del tamaño de la retícula en la solución, análisis de resultados detallados en regiones seleccionadas, etc. El programa puede contener los datos de una serie de materiales metálicos y no metálicos con propiedades conocidas por lo que es posible seleccionar el tipo directamente y el programa presenta las propiedades del material seleccionado ya almacenadas en su memoria. También es posible definir nuevos materiales con la opción “Definir nuevo material”; los materiales así definidos se almacenarán como nuevos materiales en la memoria del programa. Una vez introducidos los datos de una capa es preciso pulsar el botón “Añadir capa” para finalizar correctamente la introducción de los datos de la capa. Cada vez que se introduce una nueva capa la pantalla de la retícula se colorea con un color de define la capa introducida. Así para las dos capas verticales definidas en las pantallas de la Figura 3.10 a) y b) la retícula se presenta de la forma indicada en la pantalla de la Figura 3.11.

Figura 3.10 Pantalla de introducción de datos por capas en geometría rectangular. a) capa 1

Figura 3.10 Pantalla de introducción de datos por capas en geometría rectangular. b) capa 2

Figura 3.11 Retícula asociada a las capas definidas en la Figura 3.10 a) y 3.10 b)

Finalizada la definición de capas, el botón “Siguiente” da acceso a la pantalla de “Condiciones de contorno”, Figura 3.12. La introducción de las condiciones de contorno se hace definiendo la celda o conjunto de celdas a las que se aplica la misma condición, seleccionando tal condición e introduciendo los datos complementarios asociados a dicha condición y especificando la posición de la celda a la que se aplica (extremo inferior, superior, izquierdo o derecho). Las condiciones posibles son adiabática (no requiere datos complementarios), isoterma con temperatura fija (requiere el valor de la temperatura), flujo de calor constante (valor de dicho flujo), de convección (coeficiente de convección más temperatura de referencia), radiación (emisividad más temperatura de referencia), convección más radiación (requiere el coeficiente de convección y la temperatura de referencia para la convección, la emisividad y la temperatura de referencia para la radiación), isoterma función del tiempo y flujo de calor función del tiempo (que requieren información sobre la dependencia temporal de la temperatura y del flujo de calor, respectivamente). La condición de temperatura y flujo dependiente del tiempo permite la implementación directa de varios tipos de dependencias: preestablecidas, lineal a tramos, sinusoidal, rectangular y pulso. Las pantallas de la Figura 3.13 muestran los cuadros de introducción de estos datos. Es necesario pulsar el botón “Establecer condiciones de contorno” cada vez, para asegurar que la condición definida se ha implementado en el modelo. Si en una celda se han definido varias condiciones de contorno en la misma posición, el programa fija la última condición establecida.

Figura 3.12 Pantalla para la introducción de las condiciones de contorno. a) selección del tipo de condición

Figura 3.12 Pantalla para la introducción de las condiciones de contorno. b) selección de la posición, dentro de las celdas, a la que se aplica la condición

Figura 3.13 Pantallas típicas mostrando los diferentes tipos de condiciones de contorno a) condición de convección

Figura 3.13 Pantallas típicas mostrando los diferentes tipos de condiciones de contorno b) condición de radiación

Figura 3.13 Pantallas típicas mostrando los diferentes tipos de condiciones de contorno c) condición de convección más radiación

Figura 3.13 Pantallas típicas mostrando los diferentes tipos de condiciones de contorno d) condición de temperatura función del tiempo

Al finalizar la introducción de las condiciones de contorno, la pantalla de reticulación de celdas colorea las condiciones separadamente e indica en qué posición se han impuesto tales condiciones para que el usuario pueda hacer comprobaciones inmediatas. Ver, por ejemplo la Figura 3.14.

Figura 3.14 Pantalla de reticulación definiendo las capas más las condiciones de contorno

La siguiente pantalla, Figura 3.15a permite introducir la ventana de simulación junto con otros tiempos característicos asociados a la presentación de resultados tabulados y al paso de tiempo interno para los cálculos numéricos de PSpice. La simulación simultánea de un parámetro variable en un rango dado se especifica también en esta pantalla proporcionando el parámetro o variable elegida, Figura 3.15b, su rango de variación y el incremento en su valor. Finalizado este paso, la siguiente pantalla (Figura 3.16) permite introducir, en las primeras líneas del archivo de texto del modelo, una descripción del problema que se simula tan detallada como se precise. Esta descripción permite al usuario reconocer el modelo de manera detallada al abrir el archivo. Con esto se completa la especificación del archivo de modelo y el programa presenta la pantalla de la Figura 3.17, que contiene las opciones de visualización gráfica de PROCCA-09. Estas opciones están básicamente preparadas para modelos 3-D pendientes de integrar en el programa y se refieren al modo de presentar la retícula. Contiene opciones para integrar o no en la retícula la numeración de celdas y el nombre de la capa. El modelo ahora puede guardarse en una carpeta de directorio, Figura 3.18, u opcionalmente presentar el texto del mismo para ser modificado, y posteriormente ejecutado, en una nueva pantalla, Figura 3.19.

Figura 3.15a Pantalla de ventanas de simulación y simulación de un parámetro variable

Figura 3.15b Pantalla de opciones de selección de parámetros

Figura 3.16 Pantalla de descripción del problema modelo

Figura 3.17 Pantalla de opciones de visualización

Figura 3.18 Pantalla de guardar modelo

Figura 3.19 Pantalla de archivo “.cir” del modelo

Una ruta similar de pantallas, que consideramos no es necesario presentar, se da para geometrías cilíndricas y esféricas. Para terminar diremos que el programa integra una ayuda completa que contiene para cada tipo de problema una descripción de sus fundamentos teóricos junto con un listado de posibles errores y fallos del mismo.

III.3.1.2 Descripción del módulo CONCAL Este módulo está dedicado al cálculo y diseño de aletas pudiendo trabajar con aletas simples, anulares, rectangulares y conjuntos rectangulares aleta-pared. La aleta simple 1-D de sección arbitraria, sólo requiere introducir los datos de sección transversal y perímetro y el número de celdas horizontales; las anulares 1-D requieren también el radio interior y las rectangulares, la longitud, el espesor, y los números de celdas horizontales y verticales (permitiendo trabajar en 2-D aunque en la mayor parte de los casos prácticos el comportamiento es 1-D). Por último, los conjuntos más complejos aleta-pared requieren información sobre la geometría y reticulación de la pared y de la aleta, separadamente. Las pantallas de presentación para introducir datos, de cada una de estas opciones, se muestran en la Figura 3.20a-d. Así mismo, un ejemplo de retícula para un conjunto aleta-pared se muestra en la Figura 3.21.

Figura 3.20 a) Pantalla aleta simple

Figura 3.20 b) Pantalla aleta anular

Figura 3.20 c) Pantalla aleta rectangular

Figura 3.20 d) Pantalla aleta pared

Figura 3.21 Retícula de un conjunto aleta-pared

III.4 CRITERIOS PARA LA NUMERACIÓN DE CELDAS, NODOS Y ELEMENTOS DEL MODELO PROCCA-09 genera automáticamente la numeración de las celdas siguiendo un criterio lógico consistente en atribuir a cada una un conjunto de 4 dígitos, de los cuales los dos primeros indican la posición horizontal de la celda y los dos siguientes la posición vertical (Figura 3.22). Esta misma numeración se asigna al subcircuito correspondiente a la celda. El origen para la numeración se sitúa en la posición izquierda-inferior de la geometría. El nudo correspondiente al centro de la celda se define igual que la propia celda, los nudos de los bordes izquierdo y derecho llevan una x al final de su denominación mientras que los nudos inferior y superior llevan una y. Los nudos izquierdo e inferior de la celda tienen la misma numeración que el central mientras que el derecho y superior tienen una unidad más en la coordenada correspondiente (x e y, respectivamente). De esta forma es inmediato establecer una correspondencia entre nudos y posiciones locales del medio. Un detalle de esta numeración se muestra en la misma Figura 3.22. De esta forma es fácil identificar la posición relativa de cada punto del mallado a partir del número de celdas que contiene (tanto en geometría 1-D como 2D) y solicitar los datos de temperatura en los puntos requeridos una vez realizada la simulación. Esta definición de nudos es muy útil cuando se trata de buscar los errores o fallos de programación del archivo usando directamente los resultados de la simulación mostrados por medio de la opción “Schematics” de PSpice. En relación con la denominación de elementos del modelo, estos se definen con una letra inicial que los identifica (R, resistencia; C condensador; E, generador de tensión controlado por

tensión, etc.) seguida de los números correspondientes a la celda a la que pertenecen. Debido al diseño simétrico de la celda, se añade “izq o “der” a las resistencias para identificar su posición relativa en la celda. Con esta identificación intuitiva de elementos el usuario puede encontrar fácilmente el correspondiente a la posición requerida para solicitar el flujo de calor.

Figura 3.22 Numeración de celdas y nodos

Finalmente, en relación con los elementos de contorno se sigue igualmente una regla lógica para identificarlos. Cuando se trata de condiciones adiabáticas las resistencias que las implementan van conectadas entre el nudo correspondiente del contorno y masa pero cuando se trata de condiciones isotermas, convectivas o radiativas los generadores que las implementan van conectados a un nudo común que a su vez se une con masa por medio de un generador (pila) de tensión nula que actúa como integrador de las corrientes de todo el contorno o de la parte del entorno sometida a igual condición. De esta forma disponemos de este valor directamente leyendo la corriente de la pila. Estos nudos de referencia comunes de los generadores tiene la denominación NrefSuperiorIso NrefIzquierdaConv NrefIzquierdaRad

para la condición isoterma para la condición convectiva para la condición radiativa

En relación con la denominación de los generadores se ha seguido la regla VnulaIzquierdaiso VnulaDerechaconv VnulaInferiorrad

generador común a la pared isoterma superior generador común a la pared convectiva derecha generador común a la pared inferior radiativa

III.5 ESTRUCTURA DE LOS ARCHIVOS DE TEXTO DE MODELOS La estructura de estos archivos está dividida en los bloques siguientes, por este orden: -

Nombre del archivo (con la opción de incluir, a continuación del nombre, una descripción general del problema),

-

Parámetros físicos, geométricos y de reticulación,

-

Descripción de los subcircuitos de las celdas correspondientes a cada capa donde se especifican los componentes de los mismos y la denominación de sus nudos internos,

-

Listado de interconexión de subcircuitos, especificando el tipo de subcircuito y la numeración de nudos externos,

-

Listado de elementos que implementan las condiciones de contorno, indicando el tipo de elemento, su valore y los nudos de conexión,

-

Listado de variables a imprimir

-

Sentencias de opciones de simulación

El archivo, pues, está encabezado por el nombre que lo identifica y, opcionalmente, una descripción no limitada del problema a que se refiere (descripción introducida por el usuario en el diseño del modelo). Todas las líneas de esta sección contienen doble asterisco “**” al principio de las mismas lo que constituye una clave para ser obviadas por el programa de simulación PSpice. Ejemplo de encabezamiento: ** Ejemplo 1 ** Conducción 2-D en una placa homogénea, rectangular, que contiene un agujero rectangular ** de paredes adiabáticas. ** Los lados derecho e inferior son adiabáticos, el superior isotermo y el izquierdo convectivo. **

La sección siguiente del programa está formada por un listado de las variables que usa. Éstas se refieren a los parámetros físicos (conductividades, calores específicos, densidades de cada capa), parámetros geométricos del problema (longitudes del medio). parámetros asociados a las condiciones de contorno (coeficientes de convección, emisividades, temperaturas de referencia…) y tamaño de las celdas de cada capa (ancho y alto). La denominación de estas variables, que toman el valor dado en la especificación del problema o lo deducen de los datos de entrada, es una abreviatura de su nombre completo. Ejemplo de listado de variables: ** .PARAM .PARAM .PARAM .PARAM .PARAM .PARAM .PARAM .PARAM .PARAM

PI = 3.1415 TempRefIsoSuperior = 500 ConvCoefIzquierda = 50 ConvTempRefIzquierda = 300 capa1Ancho = 0.0025 capa1Alto = 0.002 capa1Cond = 200 capa1Dens = 6000 capa1CalorEsp = 500

**

El tercer bloque define los subcircuitos. Cada uno de ellos se corresponde a una capa, cuando todos los subcircuitos de la misma son iguales, o bien o a una celda cuando los subcircuitos de la capa son diferentes como ocurre en coordenadas cilíndricas y esféricas. La primera línea define el nombre del subcircuito seguido de la numeración de nodos internos del mismo. En problemas 1-D cada subcircuito o celda contiene 3 nudos (izquierdo, derecho y central, por este orden) más el nudo de potencial nulo (masa), mientras que en problemas 2-D los nudos son 6 (izquierdo, derecho, inferior, superior, central y masa, por este orden). Los componentes que contiene para el caso general son dos resistencias en problemas 1-D (cuatro en problemas 2-D) dispuestas simétricamente y un condensador conectado entre el centro de la celda y el nudo de referencia (masa). La denominación de estos componentes ya se indicó en el Capítulo II. Para especificar su valor se escribe la fórmula o función de las variables de las que depende entre llaves. Ejemplo de descripción de subcircuitos ** .SUBCKT capa1 1 2 3 4 5 6 Rizq 1 5 {capa1Ancho / (2 * capa1Cond * capa1Alto)} Rder 5 2 {capa1Ancho / (2 * capa1Cond * capa1Alto)} Rinf 3 5 {capa1Alto / (2 * capa1Cond * capa1Ancho)} Rsup 5 4 {capa1Alto / (2 * capa1Cond * capa1Ancho)} Ccell 5 6 {capa1Ancho * capa1Alto * capa1Dens * capa1CalorEsp} IC=300 .ENDS capa1 **

Un cuarto bloque se refiera al listado de interconexión entre subcircuitos (o celdas) en el que aparece el número total de subcircuitos existentes (perteneciente a alguno de los ya definidos en el bloque anterior) y la interconexión entre ellos con arreglo a la numeración de nodos ya explicada en la sección III.3.1. Cada subcircuito contiene seis nodos que se escriben en el orden nudo izquierda, nudo derecha (que terminan en la letra y), nudo inferior, nudo superior (que terminan en la letra x) nudo central y masa. A continuación se escribe el nombre del subcircuito asignado en el bloque anterior. El listado se organiza por bloques de subcircuitos correspondientes a la misma columna vertical, siguiendo un orden desde la primera columna hasta la última. Ejemplo de listado de interconexión entre subcircuitos para un total de 40 (horizontales)× 20(verticales) celdas o subcircuitos ** *Listado de interconexión entre celdas ** X0101 0101y 0201y 0101x 0102x 0101 X0102 0102y 0202y 0102x 0103x 0102 X0103 0103y 0203y 0103x 0104x 0103 … X0120 0120y 0220y 0120x 0121x 0120 X0201 X0202 … X0220 … … X4001 X4002 … X4020 **

0 capa1 0 capa1 0 capa1 0 capa1

0201y 0301y 0201x 0202x 0201 0 capa1 0202y 0302y 0202x 0203x 0202 0 capa1 0220y 0320y 0220x 0221x 0220 0 capa1 4001y 4101y 4001x 4002x 4001 0 capa1 4002y 4102y 4002x 4003x 4002 0 capa1 4020y 4120y 4020x 4021x 4020 0 capa1

El siguiente bloque de la estructura del archivo de modelo es el correspondiente a las condiciones de contorno de las celdas sometidas a esta condición. En modelos 2-D los contornos son los bordes izquierdo, derecho, superior e inferior del modelo más los bordes correspondientes a los huecos existentes dentro del medio. Los componentes que implementan estas condiciones, siguiendo una numeración ordenada, se especifican uno a uno de acuerdo con las reglas mencionadas en el Capítulo II de esta memoria. Dado que las líneas de programa tienen un número limitado de dígitos, es frecuente que el listado se separe en bloques de variables ocupando un gran número de líneas. El siguiente ejemplo corresponde a condiciones isotermas en el borde superior, adiabáticas en los bordes derecho e inferior y de convección en el borde izquierdo. ** *Listado de condiciones de contorno *Izquierda GIzquierdaconv01 0101y NrefIzquierdaConv VALUE={ConvCoefIzquierda*{capa1alto}*(V(0101y,0)-{ConvTempRefIzquierda})} GIzquierdaconv02 0102y NrefIzquierdaConv VALUE={ConvCoefIzquierda*{capa1alto}*(V(0102y,0)-{ConvTempRefIzquierda})} … GIzquierdaconv20 0120y NrefIzquierdaConv VALUE={ConvCoefIzquierda*{capa1alto}*(V(0120y,0)-{ConvTempRefIzquierda})} *Derecha RDerechaadi01 2113y 0 1E15 RDerechaadi02 2114y 0 1E15 … RDerechaadi25 4120y 0 1E15 *Superior VSuperioriso01 0121x NrefSuperiorIso {TempRefIsoSuperior} VSuperioriso02 0221x NrefSuperiorIso {TempRefIsoSuperior} … VSuperioriso20 2021x NrefSuperiorIso {TempRefIsoSuperior} *Inferior RInferioradi01 0101x 0 1E15 RInferioradi02 0201x 0 1E15 ... RInferioradi56 4001x 0 1E15

**

Con objeto de poder determinar, en todos los casos, no sólo el calor (la corriente) de cada uno de los elementos que definan las condiciones de contorno (es decir las corrientes que entran o salen por los contornos de las celdas sometidas a estas condiciones) sino el calor total que entra o sale porcada una de los bordes del medio, se implementan circuitos auxiliares que permiten determinar estos valores directamente. Así todos los componentes de contorno de cada borde se conectan a un generador de tensión nula, en serie con la conexión de masa, que integra las corrientes individuales. Este generador actúa como amperímetro en el circuito. El listado que implementa estos integradores siempre aparece en el modelo con independencia de que se den o no las condiciones isotermas, de convección, de radiación o mixtas, o de otro tipo. Para las condiciones dadas en el bloque anterior los circuitos integradores son los siguientes:

*Elementos de integración para las superficies isotermas Vnulaizquierdaiso NrefIzquierdaIso 0 0 Vnuladerechaiso NrefDerechaIso 0 0 Vnulasuperioriso NrefSuperiorIso 0 0 Vnulainferioriso NrefInferiorIso 0 0 * Elementos de integración para las superficies de convección Vnulaizquierdaconv NrefIzquierdaConv 0 0 Vnuladerechaconv NrefDerechaConv 0 0 Vnulasuperiorconv NrefSuperiorConv 0 0 Vnulainferiorconv NrefInferiorConv 0 0 * Elementos de integración para las superficies de radiación Vnulaizquierdarad NrefIzquierdaRad 0 0 Vnuladerecharad NrefDerechaRad 0 0 Vnulasuperiorrad NrefSuperiorRad 0 0 Vnulainferiorrad NrefInferiorRad 0 0 **

El bloque en el que se listan las variables cuyos resultados de simulación se desean obtener en forma tabulada en el archivo “.out” constituye la siguiente sección en el archivo de modelo. Por defecto, siempre se solicita la impresión de la tensión en todos los centros de las celdas durante el transitorio, de acuerdo con el intervalo de tiempo solicitado para la impresión. Cualquier otro valor que desee ser tabulado debe solicitarse añadiendo al archivo las sentencias adecuadas. Para un modelo de 40× 20 celdas el listado que aparece por defecto es el siguiente: *Listado de variables a imprimir en el archivo .out” .PRINT .PRINT .PRINT … .PRINT .PRINT .PRINT … .PRINT .PRINT .PRINT … .PRINT

TRAN V(0101,0) TRAN V(0102,0) TRAN V(0103,0) TRAN V(0120,0) TRAN V(0201,0) TRAN V(0202,0) TRAN V(0220,0) TRAN V(0301,0) TRAN V(0302,0) TRAN V(4020,0)

El último bloque contiene un grupo de sentencias fijas. La que define la ventana del transitorio de tiempos de la simulación, ya introducidos al diseñar el modelo, sentencia “.TRAN”; la asociada a la precisión requerida en los cálculos, sentencia “.OPTIONS RELTOL”; la que define el número de dígitos con que se presentan los resultados tabulados, sentencia “.OPTIONS NUMDIG”; sentencia que activa el entorno gráfico de PSpice, sentencia “.PROBE”

y, finalmente, la sentencia de cierre del archivo modelo, “.END”. Un ejemplo de

estas sentencias es ** .TRAN 10ms 2s 0 UIC .OPTIONS RELTOL 0.0001 .OPTIONS NUMDGT 4 .PROBE .END

El primero de los valores de la sentencia .TRAN define el intervalo de tiempo de impresión de los datos tabulados de salida, el segundo valor se refiere al tiempo total solicitado en la simulación y el tercero al valor a partir del cual se imprimen los resultados.

III.6 PANTALLAS DE PRESENTACIÓN DE RESULTADOS Una vez creado el archivo de modelo su simulación es inmediata pulsando el botón “Simular” que aparece en la pantalla “Archivo de texto del modelo”, que a su vez se ha originado al pulsar “Mostrar CIR” en la pantalla que contiene la retícula final del modelo. Esta acción arranca PSpice y ejecuta el modelo. Durante el tiempo que dura la simulación PSpice muestra la pantalla de la Figura 3.23 en la que aparece, por defecto, en la parte inferior derecha el intervalo de tiempo pedido a la simulación y el paso de tiempo que PSpice adopta en cada instante para la evaluación de los cálculos numéricos. PSpice varía continuamente este paso de tiempo de cálculo de acuerdo con la tendencia uniforme o cambiante de los resultados actuales, para reducir al máximo (sin merma de la precisión en los resultados) los tiempos de computación totales.

Figura 3.23 Pantalla de simulación en curso (entorno PSpice)

Si se ha incluido la opción gráfica de presentación de resultados (sentencia “.PROBE”), una vez finalizada la simulación se muestra directamente el entorno gráfico PSpice, Figura 3.24. Este entorno gráfico consiste en una cuadrícula en la que irán apareciendo la evolución temporal de las variables solicitadas (o sus valores cuando se trata de soluciones estacionarias) a través del comando “ADD TRACE” de dicha pantalla, Figura 3.25. Al pulsar este comando aparece la pantalla mostrada en la Figura 3.26 en la que se listan todas las variables posibles que pueden solicitarse (tensiones –temperaturas– en todos los nudos del modelo y corrientes –flujos de calor– en todos los elementos del circuito). Basta pulsar la variable solicitada para que ésta pase a formar parte de la lista de variables que quiere representarse. Es posible hacer operaciones entre variables, por ejemplo obtener directamente flujos totales de calor, diferencias de temperatura, etc.

Figura 3.24 Entorno gráfico PSpice

Figura 3.25 Comando de selección de variables

Figura 3.26 Pantalla de selección de variables para el entorno gráfico PSpice

El uso de opciones avanzadas permite representar en el entorno PSpice los resultados simultáneos de una variable para todos los valores del parámetro elegido (siempre que se haya usado esta opción en el diseño del modelo). También es posible, mediante el botón Append de la regleta del entorno gráfico de PSpice, Figura 3.25, representar simultáneamente resultados de una misma variable perteneciente a modelos diferentes con objeto de comparar soluciones de modelos del mismo tipo o problema. Estas representaciones son muy útiles para el diseño y optimización de conjuntos térmicos. No obstante lo anterior la forma en que PSpice presenta los resultados en su entorno gráfico está limitada ya que presenta sólo los resultados de cada variable o grupo de variables en función del tiempo, sin posibilidad de hacer representaciones instantáneas de un grupo de variables u otro tipo de representaciones tipo 2-D. Los resultados tabulados pueden obtenerse al llamar, desde PSpice o bien desde Word, el archivo “.out”. Este archivo, que también puede solicitarse para buscar los fallos de diseño del archivo de modelo ya que vienen listados al principio del mismo, contiene ordenadamente un listado de todas las variables solicitadas para la solución transitoria o estacionaria del problema. El intervalo de tiempo usado es el especificado en el diseño del modelo (sentencia .PRINT TRAN). Estos resultados pueden ser transportados con relativa facilidad a una hoja de cálculo para su manipulación y nuevas representaciones gráficas de perfiles, curvas de temperatura constante, etc. Un ejemplo de este listado para las variables temperatura en los nodos (0101) y (0204) se muestra en la Figura 3.27 donde el paso de tiempo es de 0.01 s. ****

TRANSIENT ANALYSIS

TEMPERATURE = 27.000 DEG C

**************************************************************************** TIME

V(0101,0)

0.000E+00 4.928E-04 1.000E-02 9.216E-01 2.000E-02 1.773E+00 … TIME

V(0204,0)

0.000E+00 -1.549E-04 1.000E-02 4.452E-02 2.000E-02 1.098E-01 …

Figura 3.27 Representación resultados tabulados en el archivo de texto “.out” El entorno gráfico de PROCCA-09, aunque limitado, permite acceder a isolíneas de temperatura obtenidas por interpolación de los resultados tabulados anteriores. Estas líneas están representadas por puntos y muestran la lectura directa de la temperatura en el instante solicitado del transitorio. El acceso a estas representaciones es directo una vez simulado el modelo (aparece la pantalla gráfica de PROCCA-09) o bien, si se ha cargado un archivo ya simulado a través de la ruta “Modelo” → “Procesar archivo .out de este modelo” de la pantalla “Diseñador de modelos”.

Para acceder a representaciones gráficas de mayor precisión, realizadas con MATLAB, hay que usar la pantalla mostrada en la Figura 3.29 que arranca directamente MATLAB mostrando la rutina de programación de la Figura 3.30 (en principio no manipulable) en la que se especifica los valores de temperatura que se quieren representar. Un ejemplo de representación de isotermas para un recinto rectangular 1-D, se muestra en la Figura 3.31. Finalmente PROCCA-09 incorpora una rutina de programación que permite representar animaciones de evolución de isotermas, tanto en su propio entorno como en el entorno gráfico de MATLAB. Para generarlas se usa la pantalla de la Figura 3.32, que da acceso a la pantalla de introducción de datos de la animación (Figura 3.33), valores de temperatura de las isotermas y tiempo de paso entre fotogramas de la animación. La rutina de programación se muestra en la Figura 3.34.

Figura 3.28 Curvas de isotemperaturas del entorno gráfico de PROCCA-09

Figura 3.29 Pantalla de generación de gráficos MATLAB

Figura 3.30 Rutina de programación de gráficos MATLAB

Figura 3.31 Curvas de isotemperaturas del entorno gráfico de MATLAB

Figura 3.32 Pantalla de generación de animaciones MATLAB

Figura 3.33 Pantalla de introducción de datos de animaciones MATLAB

Figura 3.34 Rutina de programación de animaciones MATLAB

III.7 EJEMPLOS DE ARCHIVOS DE MODELO Para mejor entender la estructura de los archivos, su simulación y la presentación de resultados se presentan tres ejemplos.

III.7.1 Placa rectangular con oquedad Las dimensiones de la placa son 0.1× 0.04 m, Figura 3.35. Las características físicas y geométricas del medio se muestran en la Tabla 3.1. La ecuación que rige el problema es ∂ 2T/∂x2 + ∂ 2T/∂y2 = α (∂T/∂t) y las condiciones de contorno vienen dadas por:

T(x,y=0.04 = Tsup

(borde superior isotermo)

j(0,y) = h(Ts-Tref)

(borde izquierdo convectivo)

j(x,0) = 0, j(0.1,y) = 0

(bordes inferior y derecho de la placa adiabáticos)

j(0.05,0.024
(bordes de la oquedad adiabáticos)

α = k/ρ ce es la difusividad térmica, h es el coeficiente de convección y Tref la temperatura de referencia para la convección.

Tsup

h

d f

adiab.

Jconv.

y x

a adiab. Figura 3.35 Geometría del Ejemplo 1

a: 0,1

∆x (ancho de la celda):

0,0025

b: 0,04

∆y (alto de la celda):

0,002

c: 0,01

k (conductividad térmica):

200

d: 0,006

Ce (calor específico):

500

e: 0,04

ρ (densidad):

6000

f: 0,01

Tsup (temperatura superior):

500

h (coeficiente de convección):

50

Nº celdas horizontales:

40

Tconv (temp. ref. convección):

300

Nº celdas verticales:

20

Tini:

300

Tabla 3.1 Parámetros físicos y geométricos del Ejemplo 1 (unidades en el SI)

0 3 0 1 y

c

b

e

0 2 0 1 y

adiabático

Se ha adoptado una reticulación de celdas rectangulares del mismo tamaño. La anchura de las celdas es ∆ x = 2.5 mm y la altura ∆ y = 2 mm. Las Figuras 3.36 a) y b) muestran las retículas con los nombres y numeración de celdas más las condiciones de contorno. Una sección del archivo de texto que contiene parte de las condiciones de contorno se muestra en la Figura 3.37.

Figura 3.36a Retículas del Ejemplo 1 (mostrando la denominación y numeración de celdas)

Figura 3.36b Retículas del Ejemplo 1 (mostrando las condiciones de contorno)

Figura 3.37 Archivo de texto del modelo del Ejemplo 1

Una vez simulado el modelo, que ha requerido un tiempo de paso de computación de 1.464 s en el tramo final, según se muestra en la parte inferior de la pantalla PSpice, Figura 3.38, los resultados de algunas de las variables se muestran en las Figuras 3.39 a 3.41. Una sección del archivo tabulado de salida “ejemplo1.out” se muestra en la Figura 3.42.

Figura 3.38 Sección de la pantalla de PSpice mostrando la ventana de simulación y el tiempo de computación

Figura 3.39 Transitorio de temperaturas en puntos típicos cercanos a la diagonal. Centros de las celdas 0716, 1412, 2407 y 3602

Figura 3.40 Temperaturas en puntos cercanos al centro de los lados del hueco. Centros de las celdas 2912, 2918, 2015 y 3715

Figura 3.41 Flujos de calor en celdas típicas del contorno isotermo, celdas 1020, 3420, y convectivo, celdas 0105 y 0115

El entorno gráfico de PROCCA-09 permite visualizar los puntos de las curvas isotermas. Las Figuras 3.43a y b muestran estas isolíneas para t = 20 y 50 s, respectivamente.

Figura 3.42 Archivo de salida “.out” del Ejemplo 1 (detalle)

Figura 3.43 Isolíneas. a): t = 20 s, b) : 50 s

Figura 3.43 b) Isolíneas. t=20 s

III.7.2 Cilindro hueco de dos capas La geometría del problema, 2-D, se muestra en la Figura 3.44. Las características físicas y geométricas del medio se muestran en la Tabla 3.2. La coordenada angular está degenerada. La ecuación de conducción es 1 ∂  ∂T  1 ∂2 T ∂2 T 1 ∂T + = r + r ∂r  ∂r  r 2 ∂φ 2 α ∂t ∂z 2

Las condiciones de contorno vienen dadas por T(r=0.02,0
(superficie interior isoterma)

j(r=0.09,z) = j(r,0) = j(r,0.5) = j(r=0.02,0.25
(resto de superfícies adiabáticas) adiab. adiab.

adiab. adiab.

Capa 1

Capa 2

03 01 y

01 y

03Rint 01 y 02

02 01 y R1 03 0102 y 01 y

R2

T0

Figura 3.44 Geometría del Ejemplo 2

Jconv

H (altura): Rint: R1: R2: Nº celdas Capa 1 (N1): Nº celdas Capa 2 (N2): Nº celdas verticales (Nz): T0: Tref.conv.

0,05 0,02 0,007 0,009 10 20 10 1000 300

k (conductividad térmica Capa 1): k (conductividad térmica Capa 2): Ce (calor específico Capa 1): Ce (calor específico Capa 2):

10 100 300 200

ρ

(densidad):

1000

ρ Capa 2 (densidad): ΔR Celda Capa 1: ΔR Celda Capa 2: ΔH:

5000 0,005 0,001 0,005

Capa 1

Tabla 3.2 Parámetros físicos y geométricos del Ejemplo 2 (unidades en el SI)

Se ha adoptado una reticulación de 30 (radiales)× 10(verticales) celdas. La capa interior está formada por 10 del mismo espesor radial y la misma altura, ∆ r=0.005 y ∆ h=0.05; para las 20 celdas de la capa exterior ∆ r = 0.001 y ∆ h = 0.05. La Figura 3.45 a) y b) muestra dos vistas de la retícula mientras que la Figura 3.46 muestra una sección del archivo de texto del modelo en el que puede apreciarse que se han definido un total de 30 subcircuitos ya que por ser geometría cilíndrica los elementos (resistencias y condensadores) de cada columna son diferentes (de acuerdo con lo indicado en el Capítulo II).

Figura 3.45 Retícula del Ejemplo 2. a) Nomenclatura y numeración de celdas

Figura 3.45 Retícula del Ejemplo 2. b) condiciones de contorno

Figura 3.46 Sección del archivo de texto del Ejemplo 2

Las Figuras 3.47 a 3.49 muestran los resultados de la simulación. En la primera se muestra los transitorios de temperatura en posiciones típicas del cilindro; la segunda recoge los flujos de calor en la sección de pared isoterma (ambas en el entorno gráfico de PSpice). La tercera muestra isolíneas típicas en el entorno de PROCCA-09 para un tiempo de 316 s. Finalmente, la Figura 3.50a y 3.50b muestra isolíneas obtenidas con MATLAB en dos tiempos característicos t = 50 s y t = 316 s, respectivamente.

Figura 3.47 Temperaturas en posiciones típicas de las bases del cilindro. Base inferior: centro de las celdas 0501, 0503, 0505, 0507, y 0509; base superior: ídem 1501, 1503, 1505, 1507 y 1509

Figura 3.48 Flujos en las celdas con condición de contorno isoterma. Centros de las celdas 0101, 0102, 0103, 0104 y 0105

Figura 3.49 Isolíneas en el instante t = 316 s en el entorno gráfico de PROCCA-09

Figura 3.50 a) Isolíneas con MATLAB t = 50 s.

Figura 3.50 b) Isolíneas con MATLAB t = 316 s.

III.7.3 Aleta rectangular 1-D La geometría se muestra en la Figura 3.51. Las características físicas y geométricas del medio se muestran en la Tabla 3.3. La ecuación de conducción viene dada por el conjunto de ecuaciones (2.19) a (2.20e) del Capítulo II.

jconv.

Tb

Tref.conv.

ρ , k, ce, Tini

j=0

e

y

x l

j=0

Figura 3.51 Geometría de la aleta

Las condiciones de contorno (por simetría se ha adoptado la mitad superior de la aleta) vienen dadas por

T(x=0,t) = To j(superficie de simetría,t) = 0 j(superficie lateral) = h (Ts - Tref) j(extremo) = h (Ts - Tref)

(superficie interior isoterma) (superficie de simetría adiabática) (superficie lateral convectiva) (extremo convectivo)

l = 0.020

h = 2000

e =0.002

Tini = 300

k = 100

Trefconv = 300

ρ = 2000

Tb = 500

ce = 300

Nx=40, Ny = 4

Tabla 3.3 Parámetros físicos y geométricos del Ejemplo 3

Las celdas son cuadrados de tamaño ∆ x = 0.0005 y ∆ y = 0.005. La retícula del modelo se muestra en la Figura 3.52 y una sección del archivo de testo del modelo en la Figura 3.53

Figura 3.52 Retícula y condiciones de contorno del Ejemplo 3

Figura 3.53 Sección del archivo de texto del Ejemplo 3

Figura 3.54 Transitorio de temperaturas en posiciones típicas. Centros de las celdas 0110, 0120, 0130 y 0140

Los resultados de la simulación se muestran en las Figuras 3.54 y 3.55, transitorio de temperaturas y flujos de calor en posiciones típicas de la aleta. Una representación de las isotermas para t = 0.2 s, en el entorno gráfico de PROCCA-09, se muestra en la Figura 3.56 (puede apreciarse en las mismas el comportamiento 1-D de la aleta). Por último, la Figura 3.57 muestra los isotermas para t = 100 s usando MATLAB.

Figura 3.55 Flujos de calor de convección en las celdas xx, xx, xx y xx

Figura 3.56 Isotermas de la aleta para t = 0.2 s

Figura 3.57 Isotermas en el entorno MATLAB para t = 100 s.

IV. Aplicaciones docentes y de investigación

Capítulo IV

Aplicaciones docentes y de investigación

IV.1 Introducción IV.2 Aplicaciones docentes IV.2.1 Transitorios en medios 1-D y 2-D de geometría rectangular Placa 1-D bajo enfriamiento convectivo Placa 2-D con condiciones de contorno armónicas IV.2.2 Placas 1-D bajo condiciones de convección y radiación IV.2.3 Estudio de aletas simples IV.3 Aplicaciones de investigación IV.3.1 Aislamiento de tanques esféricos multicapa IV.3.2 Introducción al diseño de aletas y conjuntos aleta-pared rectangulares

104

Condición de contorno 1

y

Condición de contorno 2

ρ , k, ce, Tini b

IV. Aplicaciones docentes y de investigación

x IV.1

INTRODUCCIÓN

Condición de contorno 3

En este capítulo se presentan aplicaciones del programa PROCCA-09 en su doble vertiente de docencia e investigación. PROCCA-09, al igual que su antecesor PRODASIM, enfocado básicamente a cálculos de aletas simples, nació básicamente con vocación docente. El objetivo era sustituir algunas de las prácticas de laboratorio de transmisión de calor impartidas en las carreras de ingeniería, muy costosas y de carácter fundamental, por prácticas con ordenador, de coste mínimo y mayor complejidad, aunque con las desventajas inherentes asociadas a la falta de contacto directo del alumno con equipos e instrumentos de medida. PROCCA-09, mucho más potente que su antecesor PRODASIM, permite el diseño de un amplio conjunto de prácticas de conducción de calor. En el presente capítulo mostramos tres aplicaciones docentes que, aunque no organizadas como prácticas de laboratorio, muestran las posibilidades del programa en relación con este objetivo: conducción en medios multicapas rectangulares, placas bajo condición simultánea de convección y radiación, y estudio de aletas simples. La exposición de las mismas, no obstante, sigue un enfoque típicamente docente: objetivos, esquema de la aplicación, descripción, realización, obtención y tratamiento de resultados, elaboración de tablas y gráficos, cuestiones y trabajos a desarrollar por el alumno, y conclusiones finales. Nuestra intención es elaborar para el curso próximo 2009-2010 un cuaderno de prácticas, con un total de 6 a 8 prácticas de dos horas de duración, que cubra parcialmente los objetivos de la enseñanza práctica de esta materia en las carreras de Ingeniería Industrial (básicamente Mecánicos y Químicos). En relación con las aplicaciones de investigación consideramos que PROCCA-09 es una herramienta muy útil en este campo aunque su alcance se ciña al campo de la conducción térmica con condiciones de contorno (no lineales) prácticamente arbitrarias. Las posibilidades del módulo CONCBA de poder usar medios 2-D, multicapas, de geometrías tanto rectangulares como cilíndricas y esféricas, lo convierte en una herramienta de optimización muy interesante en múltiples aplicaciones en los campos industrial y de la construcción, particularmente para el diseño de recintos aislados, acondicionados térmicamente; del mismo modo el módulo CONCAL es potencialmente aplicable al diseño y optimización de superficies extendidas incluyendo el estudio de la influencia de la pared desnuda, un tema cuyo tratamiento en la literatura científica es matemáticamente muy complejo y sólo accesible a investigadores especializados. Se presentan dos aplicaciones: diseño de aislamiento de un tanque cilíndrico multicapa bajo diferentes condiciones en sus paredes y estudio de optimización de aletas y de la influencia de la pared desnuda en el diseño de conjuntos aleta-pared.

105

a

IV. Aplicaciones docentes y de investigación

IV.2 APLICACIONES DOCENTES IV.2.1 Transitorios en medios 1-D y 2-D de geometría rectangular Los objetivos de esta aplicación son los siguientes: -

Estudiar, separadamente y en conjunto, la influencia de las características térmicas conductividad, k, y calor específico, ce, en el transitorio y estacionario final de un proceso de conducción de calor 2-D

-

Ídem para el coeficiente de convección. Efecto de la dependencia del coeficiente de convección con la temperatura superficial

-

Estudio del coeficiente de difusividad térmica

-

Ídem el efecto de condiciones de contorno dependientes del tiempo estudiando la influencia de los valores máximos y del periodo de la condición

-

Familiarizarse con las distribuciones de potenciales y flujos en el medio

-

Representar perfiles de temperaturas en el medio

-

Representar las isotermas durante el transitorio con los entornos gráficos de PROCCA y MATLAB

-

Hacer representaciones animadas de las isotermas con las herramientas de MATLAB

El esquema general de esta aplicación se presenta en la Figura 4.1: una placa rectangular homogénea (de densidad ρ o, conductividad térmica k y calor específico ce), con temperatura inicial constante, Tini, se somete a diferentes tipos de condiciones de contorno dependientes o no del tiempo (adiabática, isoterma y de convección) en parte o la totalidad de sus bordes. Los valores o rangos de los parámetros se muestran en la Tabla 4.1.

Condición de contorno 1 Condición de contorno 2

b

y

ρ , k, ce, Tini

x a Figura 4.1a Esquema de la geometría

106

Condición de contorno 3

IV. Aplicaciones docentes y de investigación

Figura 4.1b Modelo en red de la celda h = 100, 300 y 1000

p = 1, 1.2 y 1.3,

k = 50, 100 y 500

Trefconv = 300

ρ ce,1 = 5E5 y 1E6

Tini = 500 Tabla 4.1 Valores y rangos de los parámetros (SI)

Placa 1-D bajo enfriamiento convectivo. Por simetría, el esquema de este problema y modelo en red se muestran en la Figura 4.2 a) y b), respectivamente. Las condiciones de contorno son: ∂2 T 1 ∂T = 2 α ∂t ∂x 1 ∂  ∂T  ∂  ∂T  ∂T k r + k  = ρ ce r ∂r  ∂r  ∂z  ∂z  ∂t

T-Δx

j(0,t) = 0

Rizq

(condición adiabática en el extremo izquierdo)

jx(l1,t) = h(Ts- Trefconcv)p

(condición de convección en el extremo derecho)

Tini (x,0) = 500

(condición inicial)

La longitud de la placa es lo = 0.5 m y el espesor 0.001 m; el número de celdas es N = 50. Con todo, ∆ x = 0.01 m. La simulación se realizará para el intervalo temporal [0,20000] y los datos tabulados se imprimen cada 10 s.

Jent.x adiabático

y

Tref.conv.

ρ , k, ce j=0 (adiabático)

Tini

Jconv = h (Ts - Tref.conv.) Ts

x

adiabático l0

Figura 4.2a Esquema del enfriamiento convectivo

107

C

IV. Aplicaciones docentes y de investigación

1 mm

T-Δx

Rizq

Rder

T

Jent.x

T+Δx

Jsal.x C

Figura 4.2b Modelo en red de la celda

Los valores de las resistencias y condensador del modelo son (tomando la unidad de longitud para la profundidad): Rizq = Rder = ∆ x / (2× 0.001× k) = 10lo / k, C = ρ ce ∆ x = ρ ce lo / 50. Con estos datos el alumno deberá diseñar un modelo general con el módulo CONCBA, en el que irá modificando los valores de k, ρ , ce, h y p para obtener los modelos particulares que ejecutará en PSpice. Los resultados de la simulación se muestran en las siguientes figuras del entorno gráfico de PSpice que el alumno comentará de acuerdo con las guías establecidas en el cuaderno de prácticas. Las influencias separadas del k, h y ρ ce en los resultados del transitorio de temperaturas se muestran en las Figuras 4.3 a 4.5, respectivamente. Se han representado las temperaturas en los centros de las celdas primera y última que corresponden muy aproximadamente a los extremos izquierdo y derecho, respectivamente, de la placa. Puede observarse que un aumento de k o h acorta la duración del transitorio de enfriamiento mientras que un aumento de ρ ce lo alarga. Las curvas de flujo de calor en la superficie convectiva, para algunos casos, se muestran en la Figura 4.6. Es interesante la comparación entre estas curvas que se cruzan para determinados valores de las características térmicas, así como la comparación de la duración del transitorio, que es función de estos valores (el alumno debe sacar conclusiones generales y particulares acerca de la influencia de los parámetros a la vista de los resultados). La construcción de estas figuras se logra con el uso del comando APPEND del entorno gráfico de PSpice, que permite integrar en un mismo gráfico de tiempos, variables con el mismo nombre pertenecientes a diferentes modelos.

108

IV. Aplicaciones docentes y de investigación

Figura 4.3 Influencia de k. Temperaturas en posiciones típicas x = 0.01 (centro de la primera celda) y

0.49 (centro de la última). k = 50, 100 y 500. ρ ce = 5E5, h = 100, p = 1

Figura 4.4 Influencia de h. Temperaturas en posiciones típicas x = 0.01 (centro de la primera celda) y 0.49 (centro de la última). k = 50, ρ ce = 5E5, h = 100, 300 y 1000, p = 1

109

IV. Aplicaciones docentes y de investigación

Figura 4.5 Influencia de ρ ce. Temperaturas en posiciones típicas x = 0.01 (centro de la primera celda) y 0.49 (centro de la última). k = 50, ρ ce = 2E5, 5E5 y 1E6, h= 100, p = 1

110

Figura 4.6 Flujos de convección, p= 1. a) k = 50, h = 100, ρ ce = 5E5; b) k = 500, h = 100, ρ ce = 5E5; c) k = 50, h = 1000, ρ ce = 5E5; d) k = 50, h = 100, ρ ce = 1E6;

En relación con la dependencia del coeficiente de convección con la diferencia de temperaturas Ts – Trefconv, adoptando la expresión típica h = ho(Ts-Trefconc)p, con p una constante cercana a la unidad, es necesario manipular el generador de convección para incluir esta potencia, lo que se hace en el archivo del modelo directamente. La sentencia que implementa esta función es GDerechaconv01 5101y NrefDerechaConv VALUE={ConvCoefDerecha*0.1*((V(5101y,0){ConvTempRefDerecha})^p)}

La Figura 4.7 muestra los resultados del transitorio de temperaturas en las celdas primera y última para tres valores de p, parámetro muy determinante en la duración del transitorio y en la distribución de temperaturas en el medio. Los flujos de convección para los mismos parámetros se muestran en la Figura 4.8.

Figura 4.7 Influencia de p. Temperaturas en posiciones típicas x = 0.01 (celda 0101) y 0.49 (celda 0140). k = 50, ρ ce = 5E5, h = 100, p = 1, 1.2 y 1.4

Figura 4.8 Flujos de convección. k = 50, ρ ce = 5E5, h = 100, p = 1, 1.2 y 1.4

Una primera visión de la evolución temporal de las isotermas puede obtenerse mediante el entorno gráfico de PROCCA-09. Para dos instantes dados, estas isotermas se muestran en la Figura 4.9.

Figura 4.9 Isotermas en el entorno PROCCA. k = 50, ρ ce = 5E5, h = 100, p = 1. a) t = 100

Figura 4.9 Isotermas en el entorno PROCCA. k = 50. ρ ce = 5E5, h = 100, p = 1. b) t = 200

El alumno deberá comprobar aspectos complementarios de la aplicación tales como: -

estudiar modelos con el mismo valor de difusividad (y valores diferentes de k y ρ ce) proporcionan los mismos resultados de temperatura y flujos de calor en el transitorio

-

determinar el tiempo característico del problema, t* = lo/α y comprobará que es del orden de magnitud de la duración del transitorio

-

construir perfiles de temperatura para los diferentes modelos y representarlos

-

estudiar la influencia del número de celdas en la precisión de los resultados

-

utilizar el entorno gráfico de PROCCA-09 para observar la evolución de las isotermas

-

ídem de MATLAB para una visión más detallada de la evolución anterior y para obtener representaciones animadas

-

comentar muchos de los resultados anteriores contestando a una serie de preguntas guiadas y establecerá las conclusiones finales

Placa 2-D con condiciones de contorno armónicas. El esquema de este problema se muestra en la Figura 4.10, el modelo en red es el de la Figura 4.1b. Los valores de los parámetros se muestran en la Tabla 4.2. y T = T (t) adiabático

T0

l2

ρ , k, ce, Tini

adiabático

isotermo x adiabático l1 Figura 4.10 Esquema de la placa 2-D bajo condiciones de contorno armónicas

k = 100 ρ = 1000 ce = 300 Tini = To = 300

l1 = 0.4 l2 = 0.2 m Nx = 50 Ny = 50 Tabla 4.2 Parámetros del problema (SI)

Con todo, los valores de las resistencias y condensador de cada celda son Rizq = Rder = 0.01, Rinf = Rsup = 0.0025 y C = 9.6. Las condiciones de contorno son: jy(x,0) = -k (∂T/∂y) = 0

(condición adiabática en el borde inferior)

jx(l1,y) = -k (∂T/∂x) = 0

(condición adiabática en el borde derecho)

jy(0<x
(condición adiabática en la mitad del borde superior)

T(0,y) = To

(isoterma en el borde izquierdo)

T(l1/2<x
(temperatura función del tiempo en la mitad del borde superior)

Se asumen dos dependencias temporales armónicas de la función T(t), sinusoidal y rectangular. Los parámetros de estas dependencias se muestran en la Figura 4.11 a y b; el valor medio de temperatura es Tmed = Tmed = 300 K y la amplitud Tamp,sen = Tamp,sen = 100; los valores del periodo son to,sen = 500 s y to,rec = 1000 s. La simulación se realizará para el intervalo temporal [0,1000] s los valores de impresión del archivo tabulado de salida se obtendrán cada 1 s. T

Tmedia

Tmáx

t0

t Figura 4.11a Dependencias temporales (sinusoidal)

T t1

Tmáx

Tmín

t0

t Figura 4.11b Dependencias temporales (rectangular)

Con estos datos el alumno deberá diseñar dos modelos correspondientes a las dos dependencias temporales de la condición isoterma anterior; para ello usará la ventana que permite definir directamente este tipo de dependencias, Figura 4.12.

Figura 4.12 Ventana (CONCBA) de introducción de condiciones de contorno dependientes del tiempo

Algunos de los resultados de la simulación, temperaturas y flujos de calor en posiciones típicas del medio, se muestran en las Figuras 4.13, 4.19 (para la condición sinusoidal) y 4.20 y 4.26 (para la onda rectangular). Puede apreciarse cómo el retraso de la onda depende de la distancia del punto a a la condición de contorno armónica. Es importante observar también la forma de onda de los flujos de calor en las celdas en donde se aplican las dos condiciones de contorno diferentes. El alumno comentará éstas y otras gráficas en diferentes puntos del medio para determinar la profundidad de penetración en la placa del efecto armónico de la condición de contorno debido tanto a la amplitud máxima como al periodo de dicha condición, así como la onda de calor producida. También, mediante el entorno gráfico de PROCCA comentará la evolución (armónica) global de las temperaturas en la placa; este efecto se aprecia con mayor precisión en las representaciones gráficas de MATLAB a partir de los datos del archivo de salida .out. Las Figuras 4.27a, 4.27b y 4.27c muestran tres isotermas en la placa, para la onda sinusoidal, correspondientes a los instantes t = 200, 400 y 600 s.

Figura 4.13 Condición de contorno sinusoidal. Temperaturas en posiciones típicas cercanas a la diagonal. Centros de las celdas 2501, 2510, 2520, 2530, 2540 y 30 50

Figura 4.14 Condición de contorno sinusoidal. Temperaturas en puntos de la vertical x = 0.076. Centros de las celdas 1001, 1010, 1020, 1030, 1040 y 1050

Figura 4.15 Condición de contorno rectangular. Temperaturas en puntos de la vertical x = 0.336. Centros de las celdas 4001, 4010, 4020, 4030, 4040 y 4050

Figura 4.16 Condición de contorno sinusoidal. Temperaturas en puntos de la horizontal y = 0.038. Centros de las celdas 0110, 1010, 2010, 3010, 4010 y 5000

Figura 4.17 Condición de contorno sinusoidal. Temperaturas en puntos de la horizontal y = 0.158. Centros de las celdas 0140, 1040, 2040, 3040, 4040 y 5040

Figura 4.18 Condición de contorno sinusoidal. Flujos de calor en las fuentes dependientes del tiempo. Celdas 2650, 3450, 4250 y 5050

Figura 4.19 Condición de contorno sinusoidal. Flujos de calor en las fuentes isotermas. Celdas 0101, 0110, 0120, 0130, 0140 y 0150

Figura 4.20 Condición de contorno rectangular. Temperaturas en posiciones típicas cercanas al eje de simetría. Centros de las celdas 2501, 2510, 2520, 2530, 2540 y 3050

Figura 4.21 Condición de contorno rectangular. Temperaturas en posiciones típicas de la línea vertical x = 0.076. Centros de las celdas 1001, 1010, 1020, 1030, 1040 y 1050

Figura 4.22 Condición de contorno rectangular. Temperaturas en posiciones típicas de la línea vertical x = 0.336. Centros de las celdas 4001, 4010, 4020, 4030, 4040 y 4050

Figura 4.23 Condición de contorno rectangular. Temperaturas en puntos de la horizontal y = 0.038. Centros de las celdas 0110, 1010, 2010, 3010, 4010 y 5000

Figura 4.24 Condición de contorno rectangular. Temperaturas en puntos de la horizontal y = 0.158. Centros de las celdas 0140, 1040, 2040, 3040, 4040 y 5040

Figura 4.25 Condición de contorno rectangular. Flujos de calor en las fuentes dependientes del tiempo. Celdas 2650, 3450, 4250 y 5050

Figura 4.26 Condición de contorno rectangular. Flujos de calor en las fuentes isotermas. Celdas 0101, 0110, 0120, 0130, 0140 y 0150

Figura 4.27 Condición de contorno sinusoidal. a) Isotermas para t = 200

Figura 4.27 Condición de contorno sinusoidal. b) Isotermas para t = 400

Figura 4.27 Condición de contorno sinusoidal. c) Isotermas para t = 600

Otros tareas a realizar por el alumno son: -

trabajar con ondas de valor medio cero para ver sus efectos

-

cambiar las condiciones de contorno en los bordes de condiciones fijas

-

construir perfiles de temperatura

-

utilizar el entorno gráfico de PROCCA-09 para observar la evolución de las isotermas para condiciones de contorno diversas

-

ídem de MATLAB para una visión más precisa y para representaciones animadas

-

comentar muchos de los resultados anteriores contestando a una serie de preguntas guiadas y establecerá las conclusiones finales

IV.2.2 Placas 1-D bajo condiciones de convección y radiación Los objetivos de esta aplicación son los siguientes: -

Comparar entre sí los efectos de las condiciones de contorno de convección y radiación para diferentes rangos de los parámetros que las definen: coeficiente de convección, emisividad y temperaturas de referencia

-

Estudiar los rangos para los que dichos efectos (flujos de convección y radiación en la superficie) son de un orden de magnitud similar

-

Para una aplicación simultánea de ambas condiciones estudiar los rangos de valores de los parámetros para los que una condición se hace despreciable frente a la otra

-

Representar las isotermas durante el transitorio con los entornos gráficos de PROCCA y MATLAB

-

Hacer representaciones animadas de las isotermas con las herramientas de MATLAB

El esquema general de esta aplicación se presenta en la Figura 4.28. Una placa 1-D homogénea (de densidad ρ o, conductividad térmica k y calor específico ce), con temperatura inicial constante, Tini, y temperatura en el borde izquierdo constante, T o, se somete a diferentes condiciones de contorno, convección, radiación y convección más radiación en el otro extremo. El modelo en red es el de la Figura 4.2b. Los valores o rangos de los parámetros se muestran en la Tabla 4.3.

y l0

ρ , k, Jconv = h (Ts - Tref.conv.)

T0 (isotermo)

Jred = ε T (Ts - Tref.red.) Tref.conv., Tref.rad.

x Figura 4.28 Esquema de la geometría de la placa

k = 100 ρ = 3000 ce = 300 h1 = 10, 50, 100 y 500 ε 1 = 0.3, 0.5, 0.7 y 0.9

lo = 0.1 e = 0.01 Tini = 300, To = 1000 Trefconv = 300 Trefrad = 300 N = 50 Tabla 4.3 Valores y rangos de los parámetros (SI)

El modelo matemático y condiciones de contorno son: T(0,t) = To

(isoterma en el borde izquierdo)

j(lo,t) = h(Ts – Trefconv)

(condición convectiva en el borde derecho)

j(lo,t) = σ ε (Ts4 – Trefrad4)

(condición de radiación en el borde derecho)

j(lo,t) = h(Ts – Trefconv)+ σ ε (Ts4 – Trefrad4)

(condición de convección más radiación en el borde derecho)

Los valores de las resistencias y condensador de cada celda son Rizq = Rder = 0.001, C = 1.8. Enfriamiento por convección y radiación separadamente. El intervalo temporal de simulación es [0,25] s y el intervalo de tiempo de impresión 0.1 s. Con todo lo anterior el alumno diseñará un conjunto de modelos para cada conjunto de valores de los parámetros o un único modelo en el que manipulará los valores de los parámetros. La simulación en PSpice proporciona los siguientes resultados. La Figura 4.29 muestra el transitorio de temperaturas en posiciones típicas del medio (a) y flujos de calor en la superficie (b) para los diferentes valores del coeficiente de convección, h, mientras que la Figura 4.30 muestra los mismos resultados para la condición de radiación. Los flujos de calor en la superficie isoterma se muestran en la Figura 4.31. Puede apreciarse como estos flujos, en el estacionario, coinciden con los de las superficies de convección y radiación para cada una de las condiciones de contorno. Finalmente, la Figura 4.32 muestra varias representaciones de isotermas en el entorno gráfico de PROCCA-09 para el caso de convección solamente; en ellas puede observarse que el estacionario ya se alcanza a los diez segundos.

Figura 4.29 Condición de convección. h = 10, 50, 100 y 500. a) Temperaturas en las celdas 2501 (centro de la placa) y 5001 (extremo)

Figura 4.29 Condición de convección. h = 10, 50, 100 y 500. b) flujos de calor en la superficie de convección

Figura 4.30 Condición de radiación. ε = 0.3, 0.5, 0.7 y 0.9. a) Temperaturas en las celdas 2501 y 5001

Figura 4.30 Condición de radiación. ε = 0.3, 0.5, 0.7 y 0.9. b) flujos de calor en la superficie de radiación

Figura 4.31 Flujos en la pared isoterma. a) h = 10, 50, 100 y 500 (sólo convección); b) ε = 0.3, 0.5, 0.7 y 0.9 (sólo radiación)

Figura 4.32 Isotermas (sólo convección). a) t = 1 s

Figura 4.32 Isotermas (sólo convección). b) t = 10 s

Figura 4.32 Isotermas (sólo convección). c) t = 20 s

Enfriamiento por convección y radiación simultáneamente. El intervalo temporal de simulación es [0,100] s y el intervalo de tiempo de impresión 0.1 s. Los resultados de la simulación se muestran en las Figuras 4.33 a 4.36. La Figura 4.33 muestra el transitorio de temperaturas en posiciones típicas del medio para diferentes valores del coeficiente de convección y una misma emisividad mientras que los flujos de convección radiación para esos mismos valores se muestra en la Figura 4.34. Estos mismos resultados para un mismo h y diferentes valores de ε se muestran en las Figuras 4.35 (temperaturas en el medio) y 4.36 (flujos de convección y radiación). Puede observarse que la influencia de la emisividad en el campo térmico es muy pequeña (excepto para los flujos de radiación) en comparación con la influencia del coeficiente de convección. Las isotermas de MATLAB para diferentes tiempos se muestran en la Figura 4.37.

Figura 4.33 Condición de convección más radiac ión. h = 100, 500 y 1000, ε = 0.5. Temperaturas en las celdas 2501 (extremo isotermo), 2525 (centro de la placa) y 5001 (extremo de convección y radiación)

Figura 4.34 Condición de convección más radiación. h = 100, 500 y 1000, ε = 0.5. Flujos de convección y radiación

Figura 4.35 Condición de convección más radiación. h = 10, ε = 0.5, 0.7 y 0.9. Temperaturas en las celdas 2501 (extremo isotermo), 2525 (centro de la placa) y 5001 (extremo de convección y radiación)

Figura 4.36 Condición de convección más radiación. h = 10, ε = 0.5, 0.7 y 0.9. Flujos de convección y radiación

Figura 4.37a Isotermas. Sólo convección, h = 100. t = xx y xx s

Figura 4.37b Isotermas. Sólo radiación, ε = 0.7, t = xx y xx s

Figura 4.37c Isotermas. Convección más radiación, h = 100, ε = 0.5 y 0.7, t = xx y xx s

El alumno completará el aprendizaje de esta aplicación estudiando: -

la diferencia cualitativa y cuantitativa de los resultados de las dos condiciones para interpretar la influencia de los valores de los parámetros (y otros nuevos)

-

la influencia de la conductividad y del producto ρ ce en la radiación

-

la diferencia cuantitativa entre los resultados de la condición simultánea convección más radiación y los resultados de ambas por separado discutidos en el caso anterior

-

una ampliación de los rangos de valores para ver cuándo la influencia de una condición es despreciable frente a otra en el caso de su aplicación simultánea

-

los perfiles de temperatura en el medio (que representará)

-

los tiempos característicos

-

resultados complementarios, tales como flujos de calor en el interior del medio, cualitativamente diferentes según el valor de los parámetros del problema

-

la influencia del número de celdas en la precisión de los resultados

-

la evolución de las isotermas mediante el entorno gráfico de PROCCA-09, dicha evolución proporciona una información muy interesante de la que derivan los flujos de calor en el medio aludidos anteriormente

-

representaciones detalladas de las isotermas mediante MATLAB y representaciones animadas para estudiar sobre ellas muchas de las conclusiones obtenidas

-

el alumno, finalmente, contestará un conjunto de preguntas sobre la aplicación y establecerá las conclusiones finales

IV.2.3 Estudio de aletas simples Los objetivos de esta aplicación son los siguientes: -

estudiar el proceso de disipación térmica de una aleta simple de sección arbitraria sometida a convección, y la influencia del valor de los diferentes parámetros que lo regulan: conductividad, coeficiente de convección, longitud, sección transversal y perímetro de la aleta

-

comprender el concepto de longitud característica como parámetro de referencia para adimensionalizar la longitud real de la aleta

-

obtener los coeficientes de eficiencia y efectividad

-

estudiar la influencia de las condiciones del extremo

-

representar perfiles de temperaturas y flujos de calor en la aleta

El esquema general de esta aplicación se presenta en la Figura 4.38a. Una aleta de sección arbitraria So, perímetro Po y longitud lo (espín), que se supone trabaja en condiciones 1-D, se somete a convección en su superficie lateral, condición isoterma en la superficie de contacto con la pared, To, y condición convectiva o adiabática en el extremo; Tini es la temperatura inicial. El modelo en red para un elemento de volumen cualquiera se muestra en la Figura 4.38b y los valores de los parámetros en la Tabla 4.4.

l0 Figura 4.38a Esquema de la geometría

T0 x

S0 Figura 4.38b Modelo en red de la celda

k = 100

Trefconv = Tini = 300

ρ = 4000

lo = 0.02

ce = 300

So = 3E-6

h = 2000

Po = 6E-3

To = 500

N = 50 Tabla 4.4 Valores y rangos de los parámetros (SI)

El modelo matemático y condiciones de contorno son:

d 2∂ T − Tref.conv − m 2θ = 0, θ = ,m = 2 T0 − Tref.conv. dx

h P0 k S0

T(0,t) = To

(isoterma en el extremo derecho)

j(lo,t) = h(Te – Trefconv) o j(lo,t) = 0

(condición convectiva o adiabática en el extremo derecho)

j(x,t) = h(Tl(x) – Trefconv)

(condición de convección en la superficie lateral)

Te y Tl son las temperaturas en el extremo y en la superficie lateral, respectivamente. Con estos

Rizq

datos, ∆ x = 4E-4, y las resistencias y condensador de cada celda, Rizq = Rder = ∆ x/(2kSo) = lo(N2k So) = 4, C = So ∆ x ρ ce = 1.44E-3. El intervalo temporal de simulación es [0,1] s y el intervalo de tiempo de impresión 0.001 s.

T

jc

Con todo lo anterior el alumno diseñará dos modelos, uno con extremo adiabático y otro con extremo convectivo. La simulación en PSpice proporciona los siguientes resultados mostrados en las siguientes figuras. Las Figuras 4.39 a 4.42 muestran el transitorio de temperaturas en posiciones típicas para diferentes rangos de los parámetros h, k, So y lo, respectivamente, mientras que las Figuras 4.43 a 4.46 muestran los flujos de calor convectivos. La influencia de cada uno de los parámetros es patente y se manifiesta en condiciones estacionarias de funcionamiento. Los parámetros influyentes de la longitud de la aleta, Figura 4.42, y la sección transversal, Figura 4.45.

Figura 4.39 Influencia de la conductividad térmica, k=50, 100 y 150. Temperatura en las celdas 0101, 1001 y 3001

Figura 4.40 Influencia del coeficiente de convección, h = 1000, 1500 y 2000. Temperaturas en las celdas 1001, 2001, 3001 y 5001

Figura 4.41 Influencia de la sección transversal, So = 3E-6, 4E-6 y 5E-6. Temperaturas en las celdas 0101, 2501 y 5001

Figura 4.42 Influencia de la longitud de la aleta, lo = 0.01, 0.015 y 0.02. Temperaturas en las celdas 0101, 2501 y 5001

Figura 4.43 Influencia de la conductividad térmica, k=50, 100 y 150. Flujos de convección en las celdas 0101, 1001 y 3001

Figura 4.44 Influencia del coeficiente de convección, h = 1000, 1500 y 2000. Flujos de convección en las celdas 1001, 2001, 3001 y 5001

Figura 4.45 Influencia de la sección transversal, So = 3E-6, 4E-6 y 5E-6. Flujos de convección en las celdas 0101, 2501 y 5001

¿Otra Figura 4.45?

Figura 4.46 Influencia de la longitud de la aleta, lo = 0.01, 0.015 y 0.02. Flujos de convección en las celdas 0101, 2501 y 5001

Los perfiles de temperatura y flujos convectivos, para algunos valores concretos de los parámetros se muestran en las Figuras 4.47 y 4.48. La longitud característica se define como la parte de una aleta muy larga para la que la temperatura ha caído a un valor de To/e , con T0 la temperatura de la base y e la base de los logaritmos neperianos. A partir de las gráficas anteriores es inmediato determinar la longitud característica de las aletas correspondientes, longitud que se utiliza como parámetro de referencia para adimensionalizar la longitud real de una aleta. La influencia de convección o no en el extremo suele ser despreciable (a menos que se trate de aletas cortas); la Figura 4.49 muestra las temperaturas en celdas cercanas al extremo y flujos de convección en la última celda (y en el extremo) para un caso general de aletas en donde puede apreciarse el pequeño efecto de la convección. 6,000E+02

Temperatura (K)

5,000E+02

4,000E+02

3,000E+02

2,000E+02

1,000E+02

Curva a Curva b Curva c Curva d

0,000E+00 0,0001 0,0009 0,0017 0,0025 0,0033 0,0041 0,0049 0,0057 0,0065 0,0073 0,0081 0,0089 0,0097

Figura 4.47 Perfiles de temperatura en la aleta. a) k = 150, h = 2000, So = 3E-6 y lo = 0,02; b) k = 100, h = 2000, So = 3E-6 y lo = 0,02; c) k = 100, h = 2000, So = 5E-6 y lo = 0,02; d) k = 100, h = 2000, So = 3E-6 y lo = 0.01

Figura 4.48 Perfiles de flujos de convección. a) k = 100, h = 2000, So = 3E-6 y lo = 0,02; b) k = 100, h = 2000, So = 3E-6 y lo = 0,02; c) k = 100, h = 2000, So = 3E-6 y lo = 0,02; d) k = 100, h = 2000, So = 3E-6 y lo = 0,02

Figura 4.49 Influencia de la convección en el extremo. k = 150, h = 2000, So = 3E-6 y lo = 0.02 Curvas inferiores: temperaturas en las celdas 4501 y 5001; Superiores: flujo en la última celda y en el extremo

El cálculo de los parámetros de eficiencia puede realizarse a partir de los resultados anteriores, así como la influencia de la convección en el extremo en relación con la condición adiabática generalmente adoptada. El alumno completará el aprendizaje de esta aplicación estudiando: -

diferentes tipos de sección comparando los resultados entre sí

-

una ampliación de los rangos de valores para determinar la expresión de la longitud característica en función de los parámetros

-

la representación de diferentes perfiles de temperatura para mejor comprender el papel de la longitud característica

-

la influencia del número de celdas en la precisión de los resultados

-

la eficiencia de un grupo amplio de aletas para construir gráficos de eficiencia

-

ídem la efectividad

-

contestará un conjunto de preguntas sobre la aplicación y establecerá las conclusiones finales

IV.3 APLICACIONES DE INVESTIGACIÓN Hemos mencionado ya que PROCCA-09 puede ser aplicado para diseño e investigación. Aunque en su versión actual sólo trata de materiales homogéneos y por tanto lineales, si bien las condiciones de contorno que admite pueden ser fuertemente no lineales, modificaciones relativamente pequeñas en el programa permitirían que éste simulara medios no lineales; por ejemplo, medios con conductividad y/o calor específico dependientes de la temperatura, problemas cuyo diseño del modelo en red ya ha sido objeto de estudio en otras tesis doctorales [Alhama, 1999; Alarcón, 2001; Zueco, 2004 y Moreno, 2004] y numerosos artículos de investigación [Horno y col., 1993; González y col., 1995; Alhama y col., 1997; Alarcón y col., 2002a y 2002b; Zueco y Alhama, 2005, 2006, 2007; Zueco y col. 2006a, 2006b; Moreno Soto y col., 2007; Soto y col., 2007]. Es nuestra intención, de hecho, incorporar esta posibilidad y otras tales como problemas 3-D, de reticulación variable, etc. en futuras versiones de PROCCA-09. Consideramos, no obstante que el programa, en su versión actual, es capaz de abordar ciertos problemas de diseño e investigación en transmisión de calor, de los cuales presentamos a continuación dos aplicaciones. IV.3.1 Aislamiento de tanques esféricos multicapa Este es un problema de gran interés en la industria química y energética en general. Nos ceñimos aquí al caso de un tanque constituido por tres capas de características térmicas diferentes de las que una de ellas, la capa central, actúa básicamente como aislante mientras que las otras dos, de espesor constante, están definidas por requerimientos mecánicos, químicos o de otra índole. La Figura 4.50 representa un esquema físico del problema cuyos datos se muestran en la Tabla 4.5. El modelo en red es el de la Figura 4.1b. La pared interna está sometida a valores de temperatura conocidos mientras que el exterior se enfría por convección al ambiente. El objetivo es encontrar el espesor de aislamiento adecuado para unas condiciones dadas de temperatura en la pared exterior y otras asociadas con criterios de aislamiento externo.

Tref.conv.

Rint

e1

e2

e3

T(t)

Figura 4.50 Esquema físico del problema

jconv.

Rint = 0,1

ρ 3 = 3000

e1 = 0,005

ce,3 = 200

e2 = 0,05

k3 = 200

e3 = 0,005

N1 = 5

ρ 1 = 7500

N2 = 40

ce,1 = 200

N3 = 5

k1 = 100

h = 10,

ρ 2 = 100

Tini = 300

ce,2 = 50

Tref.conv. = 300

k2 = 10 Tabla 4.5 Valores de los parámetros (SI)

El modelo matemático está formado por las ecuaciones: Ec. (2.14)

Ecuación de conducción, i = 1,2 y 3

Tint = T(t)

(pared interior)

T1 = T2, r = Rint + e1

(contacto entre capas 1 y 2) (contacto entre capas 2 y 3)

T2 = T2, r = Rint + e1 +e2 jext = h (Ts – Trefconv)

(pared exterior)

T(r,t=0) = Tini

(condición inicial)

Pared interna isoterma y temperatura limitada en la pared exterior por el coeficiente de convección La temperatura de la pared interior del tanque es T int = 1000 K mientras que la temperatura exterior no debe exceder de un cierto valor adaptando el coeficiente de convección, h; Tini = 300. El objetivo, pues, es determinar la temperatura exterior del tanque en función de h. El modelo en red se diseña para un conjunto de valores de h utilizando la sentencia .STEP PARAM. En esta sentencia se especifica el conjunto de valores de h dando el valor inicial (200), el final (5000) y el salto entre valores (200). La estructura de la sección del archivo de texto de esta sentencia es la siguiente: .PARAM ConvCoefDerecha = 10 .STEP PARAM ConvCoefDerecha 200 5000 200

La Figura 4.51 muestra las temperaturas en el centro de la última celda, 5001, para cada valor del coeficiente de convección. Los flujos de convección se muestran en la Figura 4.52. Puede observarse la no linealidad entre h y los valores de temperatura y flujo. La gráfica de la Figura 4.53 muestra mejor la relación h = h (Tcelda 5001)

Figura 4.51 Temperaturas en el centro de la última celda (5001). h = 200, 400, 60… 5000

Figura 4.52 Flujos de convección. h = 200, 400, 600…, 5000

Figura 4.53. h = h(Tcelda 5001).

Pared interna isoterma y temperatura limitada en la pared exterior por el espesor de la capa aislante Se simula el programa para un conjunto de espesores de la capa aislante y los valores constantes del resto de los parámetros dados en la Tabla 4.5. La estructura de la sentencia de parámetro variable es .PARAM capa2Ancho = 0.00125 .STEP PARAM capa2Ancho 0.00125 0.01 0.00125

La Figura 4.54 representa las temperaturas transitorias y estacionarias de la pared externa del recipiente para un conjunto de espesores en torno al valor idóneo. Los flujos de calor por convección están representados en la Figura 4.55 para el mismo rango de espesores. La gráfica de valores estacionarios de temperatura T ext = f(e2) que facilita mejor la elección de e2 se muestra en la Figura 4.56.

Figura 4.54 Temperaturas en la celda de convección (5501), e2 = 0.00125, 0.0025, 0.00375…, 0.01

Figura 4.55 Flujos de calor de convección. e2 = 0.00125, 0.0025, 0.00375…, 0.01

e2 T celda

Figura 4.56 Temperaturas en la celda 45501 en función del espesor de la capa 2

Pared interna isoterma y temperatura limitada en la pared exterior por la conductividad térmica del aislante Se simula el modelo para h = 3000 y los valores 1, 2, 3…, 5 de la conductividad de la capa aislante, k2. El resto de los valores se mantiene constante, Tabla 4.5. La estructura de la sentencia del parámetro variable es .PARAM capa2Ancho = 0.00125 .STEP PARAM capa2Cond 1 5 1

La Figura 4.57 muestra la temperatura de la pared externa del recipiente para cada k. Los flujos de calor por convección están representados en la Figura 4.58 para el mismo rango de valores. La gráfica de valores estacionarios de temperatura Text = f(k2), que facilita mejor la elección de k2, se representa en la Figura 4.59.

Figura 4.57 Temperaturas en la celda de convección (5501). h = 3000, k2 = 1, 2, 3, 4 y 5

Figura 4.58 Flujos de calor de convección

PARAM CAPA2COND = {1, 2, 3, 4, 5}

Figura 4.59 Temperaturas en la celda 5001 en función de k2. k2 = 1, 2, 3, 4 y 5

Temperatura limitada en la pared exterior para un cambio armónico en la pared interna PULSE m Se asume el cambio de temperatura armónico representado en la Figura 4.60 en la pared interna del recipiente. L especificación de esta temperatura puede hacerse con la opción “Temperatura” dependiente del tiempo.

T Duración del pulso

Tmin

Tmax

retardo

t

periodo

tiempo de subida

tiempo de bajada

VSuperior26 0101y 0 PULSE(300 1000 0 10 10 90 200) Figura 4.60 Variación de la temperatura en la pared interna del recipiente

La influencia del espesor de la capa aislante en la temperatura de la superficie convectiva, con h = 200 y los valores de la Tabla 4.5 para el resto de los parámetros, se muestra en la Figura 4.61. La temperatura en la frontera entre las capas 1 y 2 (celda 0501) y en el centro del cilindro (celda 2501) se representa en la Figura 4.62; en ella puede apreciarse como el retraso de la onda aumenta con el espesor mientras que su amplitud disminuye con éste (efecto más apreciable al alejarse de la superficie interior).

Figura 4.61 Temperatura en el extremo para e2 = 0.00125, 0.00250, 0.00375…, 0.018, h = 200

Figura 4.62 Temperatura en las celdas 0501 y 2501, e2 = 0.00125, 0.00250, 0.00375…, 0.018, h = 200

La frecuencia de la onda pulsante influye apreciablemente en el resultado. Su efecto se muestra en las figuras para 4.63 a 4.66 de acuerdo con lo indicado en la Tabla 4.6. El aumento del tiempo de duración del pulso y del periodo de la onda causa mayores oscilaciones en las temperaturas del medio, si bien este efecto se desvanece al alejarnos de la pared interna y al aumentar el espesor de la capa aislante, pudiéndose ser despreciable en la superficie de convección.

90

Period o 200

4.61 y 4.62

50

950

2000

4.63 y 4.64

50

1950

4000

4.65 y 4.66

50

450

2000

-------------

Onda pulsante

Tiempo de subida

Tiempo de bajada

Duración del pulso

1

10

10

2

50

3

50

4

59

Figuras

Tabla 4.6 Parámetros de la onda armónica

Figura 4.63 Temperatura en el extremo para e2 = 0.00125, 0.00250, 0.00375…, 0.018. Onda pulsante 2, h = 200

Figura 4.64 Temperatura en las celdas 0501 y 2501, e2 = 0.00125, 0.00250, 0.00375…, 0.018. Onda pulsante 2, h = 200

Figura 4.65 Temperatura en el extremo para e2 = 0.00125, 0.00250, 0.00375…, 0.018. Onda pulsante 3, h = 200

Figura 4.66 Temperatura en las celdas 0501 y 2501, e2 = 0.00125, 0.00250, 0.00375…, 0.018. Onda pulsante 3, h = 200

Para el espesor más pequeño de la capa aislante (e 2 = 0.00125) y la onda pulsante 2, la influencia del coeficiente de convección se muestra en las Figuras 4.67 y 4.68 (el resto de los parámetros adquieren los valores dados en la Tabla 4.5). En la primera se muestra la temperatura de las celdas 0501 y 2501 mientras que en la segunda se muestra la temperatura cerca de la superficie de convección (celda 5001). Puede apreciarse que, debido al pequeño espesor del aislante, las variaciones de temperatura no son pronunciadas ni siquiera para valores muy grandes del coeficiente de convección. Por último, la influencia de la conductividad de la capa aislante en las temperaturas en las celdas 0501 y 2501, y 5001, se muestra en las Figuras 4.69 y 4.70, respectivamente. Para estas gráficas se ha tomado el espesor mínimo, e2 = 0.00125, h = 1000 y la onda pulsante 4 (Tabla 4.5). Puede apreciarse la enorme influencia de k2 en la amplitud de la temperatura en el extremo convectivo, amplitud relativamente pequeña para el menor valor de la conductividad.

Figura 4.67 Temperatura en las celdas 0501 y 2501, e2 = 0.00125 y onda pulsante 2. Tabla 4.5: h = 200, 400, 600…, 5000

Figura 4.68 Temperatura en la celda 5001. e2 = 0.00125 y onda pulsante 2. h = 200, 400, 600…, 5000

Figura 4.69 Temperatura en las celdas 0501 y 2501. e2 = 0.00125, h= 1000 y onda pulsante 4. k2 = 1, 2, 3, 4 y 5

Figura 4.70 Temperatura en el extremo convectivo (celda 5001). e2 = 0.00125, h = 1000 y onda pulsante 4. k2 = 1, 2, 3, 4 y 5

Todo lo anterior representa un estudio completo de este problema que puede extenderse a otros requisitos de diseño de carácter no técnico, por ejemplo, requisitos de carácter económico tales como coste mínimo de aislante o una combinación de éste con el mínimo volumen de recipiente utilizando diferentes aislantes. Todo ello para unos valores de temperatura límite en la superficie de convección o un flujo de convección máximo establecido. El diseño de modelos por medio de PROCCA-09 junto con las posibilidades de programación avanzada de PSpice permite simular el problema con tiempos de computación mínimos y seleccionar la opción de diseño más conveniente.

IV.3.2 Introducción al diseño de aletas y de conjuntos aleta-pared rectangulares Los elementos de refrigeración son una parte esencial del diseño de sistemas térmicos en los que el ahorro energético es una de los objetivos clave; tal ocurre por ejemplo en los dispositivos de refrigeración de componentes electrónicos de potencia o en los motores eléctricos o de gasolina de pequeña potencia. Las aletas son los dispositivos más comunes de refrigeración pero en su diseño es frecuente no considerar el efecto de la disipación de la pared desnuda. Los parámetros de eficiencia generalmente adoptados son la eficiencia y la efectividad pero ambos tienen limitaciones y se aplican sólo a la aleta (sin la pared); el primero no es coherente en cuanto que una mejora de la eficiencia no supone, en ocasiones, un aumento del calor disipado. En cuanto al segundo, no tiene un valor de referencia óptimo ya que puede ser incrementado en cualquier circunstancia. En relación con el diseño de conjuntos aleta-pared, la asunción de los efectos de la pared desnuda complica considerablemente los diseños de estos conjuntos y, con frecuencia, las referencias de la literatura científica que han dedicado un enorme esfuerzo a este tema [Razelos, 1979 y 2003; Aziz, 1992; Chung y Iyer, 1993; Wood y col., 1996; Rong-Hua, 1997; Kalman, 2000; Krauss y col., 2001] no proporcionan una guía clara y fácil de aplicar. El módulo CONCAL de PROCCA-09 puede ser, en este escenario, una herramienta muy útil pues permite diseñar modelos numéricos muy fiables que se computan con tiempo despreciable. Con un protocolo de programación en el que se determine la influencia de alguna o algunas de las variables en el calor disipado, es posible conseguir un diseño óptimo, tanto en aletas simples como en conjuntos aleta-pared. Aleta recta rectangular El problema de diseño se plantea en los términos de encontrar la geometría rectangular que disipa el máximo calor por unidad de gradiente térmico, para un volumen dado de aleta, problema equivalente al de encontrar el volumen mínimo de aleta que permite la disipación de un calor previamente especificado. Así, partiremos de un volumen dado Vo de aleta y adoptaremos una profundidad de aleta de valor unidad, buscando cuál es la relación de aspecto l/e (longitud/semi-espesor), que permite la máxima disipación de calor por unidad de gradiente térmico. La geometría del problema, que aprovecha su simetría, se muestra en la Figura 4.71. jconv Tref.conv.

jconv

z0

Tb

e

j=0

y x

l Figura 4.71 Geometría de la aleta rectangular

El modelo matemático, 2-D, viene dado por el conjunto de ecuaciones (2.19) a (2.20d) más la ecuación  ∂θ  −k   ∂y   = 0, o < y < e, x = l que corresponde a extremo adiabático.  

Los valores numéricos del problema se muestran en la Tabla 4.7. Por simplicidad, el gradiente térmico entre Tb y Trefconv es la unidad, ya que la solución es independiente de las unidades de temperatura. También se ha trabajado en 2-D a pesar de que la mayor parte de aletas que trabajan de forma eficiente pueden aproximarse aceptablemente con modelos 1-D. Como se sabe, en relación con los flujos estacionarios en los que estamos interesados, los valores de ρ y ce no influyen en sus resultados. Vo = 1E-4 z=1 k = 50 ρ = 7000 ce = 300

h = 100 Tb = 1 Tref.conv. = 0 Nx = 50 Ny = 10 Tabla 4.7 Parámetros del problema (SI)

Una vez diseñado el modelo mediante el módulo CONCAL, con valores arbitrarios de e y l (tales que el producto 2e× l vale Vo, por ejemplo e = 0.0025, l = 0.02), se corrigen las sentencias que especifican los parámetros e y l para definir estos en función del volumen, el cual se especifica mediante una sentencia añadida ya que este parámetro no está incluido en la entrada de datos del módulo. La Figura 4.72 presenta, en primer lugar los flujos de calor transitorios totales de la pared isoterma y los flujos disipados, para diferentes longitudes de aleta y un mismo espesor. Puede observarse que ambos convergen en el estacionario, según el detalle de la Figura 4.73. Las temperaturas cerca del extremo (celda 5005), para las mismas longitudes de aleta y espesor se muestran en la Figura 4.74.

Figura 4.72 Flujos transitorios de la base isoterma y flujos disipados para distintas longitudes de aleta, l = 50× 0.0001, 50× 0.0012 a intervalos de 50 × 0.0001,. e = 0.0025

Figura 4.73 Flujos estacionarios de la pared isoterma y disipados. l = 50× 0.0001, 50× 0.0012 a intervalos de 50 × 0.0001, e = 0.0025

Figura 4.74 Temperaturas en la celda del extremo (celda 5005). l = 50× 0.0001, 50× 0.0012 a intervalos de 50 × 0.0001, e = 0.0025

Para tener una idea de la longitud característica de la aleta definida como la coordenada de una aleta muy larga para la que el valor de la temperatura desciende a Tb/e, la Figura 4.75 muestra el perfil temperaturas de esta aleta mientras que las isotermas representadas por medio de MATLAB se muestran en la Figura 4.76. El valor de la longitud característica es aproximadamente 0,036 m. Por otra parte, el comportamiento de la aleta, cuyo Biot transversal es Bit = 0,005, puede aproximarse a 1-D como muestra la Figura 4.77 (temperaturas en puntos de una misma sección de la mitad de la aleta).

Figura 4.75 Perfil de temperaturas de una aleta larga, l = 0.2 y e = 0.0025 (resto de los parámetros Tabla 4.6)

Figura 4.76 Isotermas de la aleta larga, longitud 0.2 (aleta larga) y e = 0.0025. (Dejar espacio para 2)

Figura 4.77 Temperaturas en puntos de la misma coordenada, celdas 2501, 2502, 2503…, 2510, longitud 0.2 (aleta larga) y e = 0.0025

Para encontrar el perfil óptimo se trabaja con aletas de un mismo volumen (proceso de optimización). La estructura de las sentencias del archivo de modelo, que afectan a los parámetros e y l, pasan a ser definidas por medio de las dimensiones de la celda mediante funciones dependientes del volumen. Este conjunto de sentencias tiene la forma .PARAM volumen = 0.00005 .PARAM capa1Ancho = 0.0001 .PARAM capa1Alto = {volumen*0.001/capa1ancho} .STEP PARAM capa1ancho 0.0001 0.0018 0.0001

La Figura 4.78 muestra el calor disipado, estacionario y transitorio, para diferentes relaciones de aspecto; el eje de abscisas es la longitud de la aleta. Puede apreciarse que al aumentar la longitud, al principio aumenta el flujo disipado, pero a partir de una cierta longitud de aleta el flujo disminuye. Esta representación, que permite elegir la aleta óptima (lopt ≅ 0.03, eopt ≅ 0.00166) puede obtenerse mejor de la curva de flujos de la Figura 4.79 (punto de máxima disipación).

Figura 4.78 Flujos de convección totales aleta de volumen constante (Tabla 4.6)

I(VnulaArrConv)

Figura 4.79 Perfil de flujos disipados en función de la longitud de la aleta de volumen constante (Tabla 4.6)

Para la longitud óptima los perfiles de temperatura y flujos de convección se muestran en las Figuras 4.80 y 4.81. Puede observarse que la temperatura en el extremo no es despreciable en relación con salto térmico Tb Trefconv. En estas mismas figuras se representan las temperaturas y flujos de convección de la misma aleta bajo condición de convección en el extremo, con objeto de apreciar la influencia relativa de cada condición en los resultados. Como se ve, esta influencia es mínima.

Figura 4.80 Perfil de temperaturas en la aleta óptima, l = 0.03. a) extremo adiabático (curva superior). b) extremo convectivo (curva inferior)

a)

Figura 4.81 Perfil de flujos de convección en la aleta óptima, l = 0.03. extremo adiabático (curva superior). b) extremo convectivo (curva inferior)

Siguiendo la línea expuesta, es inmediato estudiar la influencia de los parámetros k y h (conjuntamente o por separado) en la determinación de la relación de aspecto óptima. Las Figuras 4.82 y 4.83 muestran los flujos disipados para diferentes relaciones de aspecto con k = 100 (Figura 4.82) y h = 200 (Figura 4.83). En ambas Figuras las longitudes varían entre 0.005 y 0.09 en pasos de 0.005 m (sentencia .STEP PARAM capa1ancho 0.0001 0.0018 0.0001). Puede apreciarse que la relación óptima para estos casos es muy similar a la encontrada anteriormente, lo que permite deducir que los valores individuales de h y k no cambian apenas el valor de la longitud óptima.

Figura 4.82 Flujos de convección totales aleta de volumen constante (Tabla 4.6), k = 100

Figura 4.83 Flujos de convección totales aleta de volumen constante (Tabla 4.6), h = 200

Conjuntos aleta-pared rectangulares La consideración de la pared complica el problema por el número de variables en juego. Entre las nuevas variables que pueden considerarse están: la conductividad de la pared, su espesor, la longitud de pared desnuda y el coeficiente de convección de la pared desnuda; todas ellas pueden incorporarse al módulo CONCAL y seleccionar el número de celdas adecuado en cada problema. Esta selección, en relación con el número de celdas de la aleta y la pared desnuda, no es independiente pues la anchura y altura de todas las celdas debe ser la misma a fin de sus bordes se acoplen en las líneas de intersección entre la pared y la aleta y las secciones de la pared separadas por la horizontal de la superficie de la aleta. La optimización en este caso no siempre es posible pues depende de los valores de los parámetros del problema. Además es preciso comprobar siempre que la aleta óptima es efectivamente eficiente (el conjunto aleta-pared debe disipar más con la aleta que sin ella). Nos ceñimos aquí a la posibilidad de estudiar estos conjuntos mediante el módulo CONCAL, presentado algunos resultados. La geometría del problema se muestra en la Figura 4.84, las ecuaciones de gobierno son las (2.19) a (2.20d) del Capítulo II más la condición adiabática del extremo, y la lista de valores de las variables en la Tabla 4.8. La Figura 4.85 muestra la geometría de las celdas del conjunto (a) y las condiciones de contorno (b).

b

ep y

J=0 j=0

Figura 4.84 Geometría del conjunto aleta-pared

ce,p = 300

ancho de celdaa = ancho de celdap = alto de celdaa = alto de celdap = 0.0005,

ce,a = 300

a

z=1

ep = ea = 0.002

hp = ha = 100 Tb = 1

lp = 0.005

Trefconv = 0

la = 0.018

Nx,p = 4

kp = ka = 50

ρ p = ρ a = 7000

jconv.

Nx,a = 36

Tb

Ny,p = 10 Ny,a = 4

Tabla 4.8 Valores de las variables (los subíndices ‘p’ y ‘a’ se refieren a pared y aleta, respectivamente)

j=0

Figura 4.85. a) Estructura de celdas del modelo

lt

Figura 4.85. b) Condiciones de contorno

La simulación proporciona los siguientes resultados. La Figura 4.86a muestra el transitorio de temperaturas en la superficie de convección de la aleta (centros de las celdas superficiales). Puede observarse que las temperaturas cerca del extremo tienen valores apreciables, en relación con el gradiente térmico, para esta longitud de aleta. El transitorio en nudos de la superficie de convección de la pared desnuda se muestra en la Figura 4.87b. Los valores estacionarios son altos en comparación con los anteriores por la cercanía de estos nudos con la pared interna isoterma. Los flujos de convección en las celdas de ambas paredes se muestran en la Figura 4.87 mientras que los flujos totales de convección se muestran en la Figura 4.88 donde puede observarse que la aleta disipa más calor que la pared desnuda. Por último, la Figura 4.89 muestra las temperaturas en una sección de la aleta para apreciar el pequeño efecto 2-D del conjunto.

Figura 4.86 a) Temperaturas en la superficie de convección de la aleta, celdas 0504 a 4004

Figura 4.86 b) Temperaturas en la superficie de convección de la pared desnuda, nudos 0405 a 0410

Figura 4.87 a) Flujos de convección en las celdas (0504 a 4004) de la aleta

Figura 4.87 b) Flujos de convección en las celdas de la pared desnuda, celdas 0405 a 0410

Figura 4.88 Flujos totales de convección en la aleta y en la pared desnuda

Figura 4.89 temperaturas en la columna 18 de la aleta, celdas 1801, 1802, 1803 y 1804

La Figura 4.90 representa los flujos totales de convección (el de la aleta, que se estabiliza más tarde, siempre es mayor que el de la pared desnuda) para diferentes relaciones entre el alto de la celda en la aleta (incluyendo la zona no desnuda de la pared) y el alto de la celda de la pared desnuda. Ambos flujos de calor crecen al amentar esta relación. El efecto 2-D para estas mismas relaciones puede apreciarse en las temperaturas de las celdas del extremo de la aleta, Figura 4.91. Este efecto es pequeño aunque aumenta al crecer el espesor de la aleta (o disminuir el alto de la pared desnuda). Las temperaturas de la pared desnuda crecen al alejarnos de la aleta aunque son temperaturas altas, cercanas a las de la cara interior de la pared; la influencia de la relación “alto de la celda en la aleta/alto de la celda de la pared desnuda” se muestra en la Figura 4.92.

Figura 4.90 Flujos totales de convección para diferentes cocientes “altura de celda de aleta/altura de celda pared desnuda”: 0.00020.0008; 0.00025-0.00075; 0.0003-0.0007; 0.00035-0.00065; 0.0004-0.0006; 0.00045-0.00055; 0.0005-0.0005; 0.00055-0.00045; 0.0006-0.0004; 0.00065-0.00035; 0.0007-0.0003

Figura 4.91 Temperaturas en el extremo de la aleta (celdas 4001, 4002, 4003 y 4004) para los cocientes “altura de celda de aleta/altura de celda pared desnuda” anteriores

Figura 4.92 Temperaturas en dos celdas típicas de la superficie convectiva de la pared, 0405 (más cercana a la aleta) y 0410 (más alejada)

La influencia del ancho de la celda (que afecta más, lógicamente, a la longitud total de la aleta que a la de la pared) en los flujos totales de convección de la aleta (superior) y de la pared (inferior) se muestra en la Figura 4.93. La disipación de la aleta (curvas superiores) aumenta más que la disipación de la pared desnuda (curvas inferiores) al incrementarse este parámetro. La influencia del ancho en la temperatura de la superficie de convección, para dos puntos típicos de la aleta se muestra en la Figura 4.94, mientras que la temperatura en las celdas del extremo, para los mismos valores del parámetro, se muestra en la Figura 4.95.

Figura 4.93 Flujos de convección de aleta (superior) y pared (inferior) para diferentes anchos de celda: 0.0005, 0.0006, 0.0007, 0.0008, 0.0009, 0.001

Figura 4.94 Temperaturas en celdas típicas (1004 y 4004) de la superficie convectiva de la aleta para los mismos anchos de celda

Figura 4.95 Temperaturas en las celdas del extremo, 4001, 4002, 4003 y 4004, para los mismos valores de anchos de la figura anterior. Las aletas más largas (donde el efecto 2-D es menor) son las curvas inferiores.

La influencia de la conductividad térmica en los flujos de convección en la aleta (curvas superiores) y en la pared desnuda (curvas inferiores) se muestra en la Figura 4.96. Un aumento de k hace crecer la disipación. La influencia de k en las temperaturas de las celdas del extremo, en celdas típicas de la superficie convectiva de la aleta y en celdas convectivas de la pared desnuda se muestra en las Figuras 4.97, 4.98 y 4.99, respectivamente. Con toda esta información y de acuerdo con los criterios y parámetros establecidos es posible el diseño de estos conjuntos siguiendo protocolos de optimización adecuados.

Figura 4.96. Flujos de convección en la aleta (curvas superiores) y en la pared desnuda (curvas inferiores). Influencia de la conductividad k = 50, 100, 150 y 200

Figura 4.97 Temperaturas en el extremo derecho (celdas 4001, 4002, 4003 y 4004) para los mismos valores de k

Figura 4.98 Temperaturas en nudos típicos (1004 y 4004) de la superficie convectiva para los mismos valores de k

Figura 4.99. Temperaturas en nudos típicos de la superficie convectiva de la pared (0505y y 0510y) para los k anteriores

V. Contribuciones

Capítulo V

CONTRIBUCIONES

164

V. Contribuciones

V. CONTRIBUCIONES

1.- La contribución fundamental de esta memoria es la realización de un programa de simulación, PROCCA-09, que conjuga la técnica del método de redes, capaz de formular ecuaciones diferenciales en diferencias finitas (discretizadas en el espacio) en términos de ecuaciones formalmente equivalentes de circuitos eléctricos ideales (e implementar estos circuitos) , con la potencia de cálculo numérico de programas comerciales de resolución de estos circuitos. PROCCA-09 permite: - diseñar modelos en red, - modificar modelos mediante el acceso directo a los mismos, lo que permite mediante ligeras modificaciones implementar características del problema que no se recogen explícitamente en PROCCA-09 mediante nuevas líneas de programa, - simularlos los modelos en PSpice, - representar los resultados de la simulación en forma gráfica o tabulada. 2.- El manejo de PROCCA-09 se ha desarrollado bajo entorno “windows” usando el programa C# y permitiendo así al usuario acceder al diseño de modelos de una manera intuitiva y ordenada, disponiendo de una ayuda suficiente. El programa contiene dos módulos diferenciados, uno de carácter general y otro específico de aletas, CONCBA y CONCAL, respectivamente. 3.- La simulación de modelos es directa desde PROCCA-09 ya que éste arranca automáticamente el simulador OrCAD o PSpice, permitiendo e acceder a su entorno gráfico y tabulado una vez realizada la simulación. 4.- PROCCA-09 dispone asimismo, además de su propio entorno gráfico capaz de representar las isotermas en cualquier instante del transitorio así como animaciones de estas representaciones, rutinas auxiliares que arrancan MATLAB para acceder a soluciones gráficas más precisas, incluyendo también animaciones de isotermas. Con todo el usuario puede observar y aprender los fenómenos asociados a la conducción térmica con gran detalle 5.- La extensión del programa PROCCA-09 incluye: i) en relación con el tipo de problema -

la solución transitoria y estacionaria, permitiendo elegir las ventanas de simulación,

-

el uso de geometrías rectangulares 1-D y 2-D, cilíndricas 1-D y 2-D y esféricas 1-D,

165

V. Contribuciones

-

la implementación de medios multicapa extendidos a toda o parte de la geometría del problema, así como la implementación de regiones de particulares de geometría casi arbitraria de característica térmicas diferentes (incluyendo oquedades en el medio),

-

la posibilidad de reticular con un número de celdas, global y para cada capa, elegido por el usuario,

ii) en relación con las condiciones de contorno, el programa permite implementar las siguientes condiciones: -

Isoterma

-

De flujo de calor constante

-

Adiabática

-

De convección

-

De radiación

-

De convección más radiación

-

Temperatura dependiente del tiempo (diferentes tipos de dependencia: sinusoidal, rectangular, lineal a tramos y pulsos)

-

Flujo de calor dependiente del tiempo (diferentes tipos de dependencia: sinusoidal, rectangular, lineal a tramos y pulsos)

-

Otras condiciones implementadas por programación en el archivo del modelo

iii) en relación con las posibilidades de programación PROCCA-09 aprovecha las ventajas de programación avanzada de OrCAD o PSpice, que incluye -

posibilidad de presentar soluciones conjuntas de diferentes problemas,

-

ídem de ejecutar simulaciones correspondientes a diferentes valores de uno más parámetros,

-

la posibilidad de representación gráfica directa de funciones arbitrarias cuyos argumentos son las variables (flujos y temperaturas) en cualquier región del medio,

-

todas las posibilidades de representación gráfica contenidas en OrCAD o PSpice (elección de rangos de representación, marcadores numéricos, colores de gráficos, etc.

iv) en relación con la presentación de resultados de la simulación -

permite representar las temperaturas y flujos de calor en topo el medio (centros y bordes de los elementos de volumen),

-

también los valores de estas variables para todos los componentes del medio y de las condiciones de contorno,

-

flujos integrados en cualquier región del contorno definida por las mismas condiciones

166

V. Contribuciones

6.- Dado que entendemos las aplicaciones del programa pueden enfocarse bien a docencia o a investigación, como contribución última de esta tesis se han presentado aplicaciones en ambas vertientes que justifican ampliamente estos objetivos. En estas aplicaciones se ha cubierto un gran número de posibilidades de aplicación del programa. Entre ellas, en su vertiente docente se pueden citar: -

influencia de las características térmicas conductividad, calor específico y densidad en las soluciones,

-

influencia de las condiciones de contorno, coeficiente de convección, emisividad, etc,

-

influencia separada de las condiciones de convección y radiación, para discernir la influencia relativa de cada condición en la solución del problema,

-

influencia de los parámetros asociados a condiciones de contorno definidas por temperaturas dependientes del tiempo, amplitud, frecuencia y forma de onda,

-

estudio de los retrasos causados en problemas de conducción armónica,

-

estudio de los parámetros que influyen en la disipación térmica de una aleta simple: perímetro, sección transversal y coeficiente de convección. Influencia de la condición del extremo, etc., mientras que en su vertiente de investigación citaremos:

-

diseño de aislamiento de tanques cilíndricos multicapa bajo diferentes supuestos de temperatura en la cavidad interior del tanque, incluyendo cambios de tipo armónico. Limitación de la temperatura de la pared externa del tanque controlando diferentes parámetros (espesor y conductividad térmica de la capa aislante y coeficiente de convección de la pared externa),

-

optimización de aletas (geometría óptima de un volumen de aleta para disipar la máxima cantidad de calor),

-

estudio de la influencia de los diferentes parámetros del conjunto aleta-pared en la disipación térmica del conjunto.

167

- REFERENCIAS: Alarcón, M. Transporte de calor en sistemas con aletas. Coeficientes de rendimiento y red de transferencia. Tesis doctoral, Universidad Politécnica de Cartagena (2001) Alarcón, M., Alhama, F. y González-Fernández, C.F. Transient heat conduction in a fin-wall assembly under harmonic excitation. Network thermal admittance. Heat Transfer Engineering, 23: 31-43 (2002a) Alarcón, M., Alhama, F. y González-Fernández, C.F. Time-dependent heat transfer in a fin-wall assembly. New performance coefficient: Thermal reverse admittance. Int. J. Therm. Sci., 41, 386-395, (2002b) Alarcón, M., del Cerro Velázquez, F. y Alhama, F. Modelización de procesos térmicos en máquinas frigoríficas mediante el método de simulación por redes. I CYTEF, (Cartagena) Libro de resúmenes, 34. (2002) Alarcón, M., del Cerro Velázquez, F. y Alhama, F. Modelización de procesos térmicos en máquinas frigoríficas mediante el método de simulación por redes. I CYTEF, Cartagena. Libro de resúmenes, 34 (2002a). Publicado en: Avances en ciencias y técnicas del frío Alarcón, A., del Cerro Velázquez, F. y Alhama F. Modelización de procesos térmicos en máquinas frigoríficas mediante el método de redes. Avances en ciencias y técnicas del frío, 109-115. UPCT y Sociedad Española de Ciencias y Técnicas del Frío (2003a) ´ Albahari, B., Drayton, P y Merril, B. C# Essentials, Editorial O’Reilly (2001) Alcaraz D., Zueco J., Alhama F., y del Cerro Velázquez, F. Un modelo simple para la solución de problemas de transmisión de calor en medios ortotrópicos bajo condiciones de contorno arbitrarias. Información Tecnológica 14(2), 27-34, May. (2003b) Alcaraz, D., Alhama, F., Zueco, y del Cerro Velázquez, F. Solución de problemas de transmisión de calor en medios compuestos bajo condiciones de contorno arbitrarias. V CIDIM, IV CONIM, (Venezuela), Proceedings, 1, 1333-1338 (2001) Alcover, V., Mareca, P. and F. Alhama. A network implementation of the Chua oscillator. Characterization of the Chua phase coherence. Proceedings. Dynamics Days. Greek-Italia (2006) Alhama, F. Estudio de respuestas térmicas transitorias en procesos no lineales de transmisión de calor mediante el método de simulación por redes. Tesis doctoral, Universidad de Murcia (1999) Alhama, F., López-Sánchez, J. F., y González-Fernández, C. F. “Heat conduction through a multilayered wall with variable boundary conditions”. Energy, 22, 797-803. (1997)

168

Alhama, F., Zueco, J., Alcaraz, D. y del Cerro Velázquez, F. Un modelo simple para la solución de problemas de transmisión de calor en medios ortotrópicos bajo condiciones de contorno arbitrarias V CIDIM, IV CONIM, (Venezuela), Proceedings, 1, 867-872 (2001) Archer, T. Inside C#, Editorial Microsoft (2004) Aziz, A. Optimum dimensions of extended surfaces operating in a convective environment. Appl. Mech. Rev. 45 (5), pp. 155-180, (1992). Baker, W. E. y Shortt, D. J., “Integrated electrical/thermal component modeling”. Naval Res. Lab., Washington (1990) Baxter, D. C., “The fision times of slabs and cylinders”. ASME Paper No. 61-WA-179, (1961) Bayazitoglu, Y. Y Özisik, M.N., “Elements of heat transfer”, Caps. 4 y 11, McGraw-Hill, New York (1988) Bejan, A., “Heat transfer”. John Wiley & Sons, Inc. New York (1993) Bejan, A., Convection Heat Transfer, 2ª ed, John Wiley and Sons, New York (1995) Bello, V. G., “Electrical models of mechanical units widen simulator´s scope”. Electronics Design News, March (1991) Bennett, C. O. Y Myers, J. E., “Momentum, heat and mass transfer”. 3ª ed. McGraw-Hill International Editions, Singapore (1982) Berezin, I.S. y Zhidkov, N.P. “Computing methods II”. Pergamon, Oxford (1995) Biot, J. B., “Traité de Physique”. 4 518, París (1816) Bonilla, C. F. Y Strupczewsky, A. L., “An electric analog computer for nuclear fuel shipping cask fire tests”. Nuclear Structural Eng., 2, 40-47 Amsterdam (1965) Brodkey, S, Robert y Hershey, C. Harry, “Transport phenomena”, Cap. 11, McGraw-Hill, New Grigull, U. y Sandner, H., “Heat conduction”.Springer-Verlag, Berlin (1984) Carslaw, H. S. y Jaeger, J. C., “Conduction of heat in solids”. 2ª ed., Oxford Science Publications (1959) Chapman, A. J., “Heat transfer”. 4ª ed. Macmillan Publishing Company, New York (1984) Crandall, S.H., “Engineering analysis”. pp.376-380, McGraw Hill, New York (1956) Crank, J. y Nicolson, P., “A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type”. Proc. Camb. Philos. Soc., pp. 189-191, Wiley-Interscience, New York (1967) Chung, B.T.F. and Iyer, R., Optimum design of longitudinal rectangular fins and cylindrical spines with variable heat transfer coefficient, Heat transfer engineering, vol. 14, no. 1, pp. 31-42, 1993. 169

Chung, K. S., “The fourth-dimension concept in the finite element analysis of transient heat transfer problems”. Int. J. Numer. Methods Eng., 17, 315-325 (1981) Del Cerro Velázquez, F., Gómez-Lopera, S. y Alhama, F. A powerful and versatile educational software to simulate transient heat transfer processes in simple fins. CAEE, Computer Application in Engineering Education 16, 72-82 (2008) Del Cerro Velázquez, Gómez-Lopera, S. y Alhama F. Solución de problemas lineales de transmisión de calor en medios multicapa mediante la analogía RC y el código Pspice. XVI CNIM Proceedings. León (2004a) Del Cerro Velázquez, Gómez-Lopera, S. y Alhama F. Software para el diseño de aletas simples basado en la analogía RC. XVI CNIM Proceedings, León (2004b) Del Cerro Velázquez, Gómez-Lopera, S. y Alhama F. Programa para la resolución numérica de problemas lineales de transmisión de calor en medios multicapa mediante la analogía RC y el código Pspice. Anales de Ingeniería Mecánica, Vol. 1, 749-756, (2004c) Del Cerro Velázquez, Gómez-Lopera, S. y Alhama F. Simulador para el diseño de aletas simples basado en la analogía RC. Anales de Ingeniería Mecánica, Vol. 1, 743-748. (2004d) Del Cerro Velázquez, F., A. Campo and F. Alhama. The Teaching of Unsteady Heat Conduction Using the Thermoelectric and the Code PSPICE. Proc. of International Network on Engineering Education and Research. INEER, Praga (2004e) Del Cerro Velázquez, F., Alhama, F. and Campo, A. The Teaching of Unsteady Heat Conduction Using the Thermoelectric and the Code PSPICE. Engineering Education and Research “Progress Through Partnership”. VSB-TUO, Ostrava. 79-86 (2004f) Del Cerro Velázquez, F. and Alhama, F. Teaching Coupled Differential Equations by Network Method. Proc. of International Network of Engineering Education and Research. INEER, Praga (2004g) Del Cerro Velázquez, F. and Alhama, F. Teaching Coupled Differential Equations by Network Method. Engineering Education and Research “Progress Through Partnership”. VSB-TUO, Ostrava. 815-821 (2004h) Fried, I. Y Metzler, J., “Conjugate gradients in a finite element discretization”. Int. J. Numer. Methods Eng., 12 1329-1342 (1978) García Hernández, M. T., “Simulación digital de procesos de electrodo. Método de redes”. Tesis Doctoral, Dpto. Física Aplicada, Universidad de Granada (1994) Gardner, K.A., “Heat exchanger tube sheet of temperatures”. Refiner. Nat. Gas. Manuf., 21, 7177 (1942) Gebhart, B., “Heat conduction and mass diffusion”. McGraw-Hill Book Co., Singapore (1993) 170

Gómez-Lopera S.A., Del Cerro Velázquez, F. y Alhama, F. Analogías de redes eléctricas como herramienta educativa para el estudio de diversos fenómenos físicos. VII Congreso Int. Sobre Investigación en Didáctica de las Ciencias. Granada (2005a) Gómez-Lopera, S.A., Del cerro Velázquez, F. y Alhama, F. Analogías de redes eléctricas como herramienta educativa para el estudio de diversos fenómenos físicos. Rev. Enseñanza de las Ciencias (Revista de investigación y experiencias didácticas), ICE (UAB), Número Extra. P. 1-4 (2005b) González-Fernández, C. F., García Hernández, M. T. y Horno, J., “Computer simulation of a square scheme with reversible and irreversible charge transfer by the network method”. J. Electroanal. Chem., 395, 39-44 (1995) González-Fernández C. F. and Alhama, F. Heat Transfer and the Network Simulation Method. Horno J. Ed. Research Signpost, Kerala (2002) Gunnerson, E. A programmer’s introduction to C#, Ed. Apress (2002) Grigull, U. Y Sandner, H., “Heat conduction· , Cap. 6, Springer-Verlag, Berlin (1984) Hamill, D. C., “Learning about chaotic circuits with pspice”. IEEE Transactions on education, 36, 28-35 (1993) Handbook of Chemistry and Physics. 73 ed., Lide, D. R. Editor. CRC Press (1992) Harper, D.R. and Brown, W.B., “Mathematical equations for heat conduction in fins of air cooled engines”. NACA Report 158, (1922) Herbert, D. B., “Simulations differential equations with pspice2”. IEEE Circuits and devices, 8, 11-14 (1992) Holman, J. P., “Heat transfer”. 5º ed. McGraw-Hill, New York (1981) Horno, J., González Caballero, A., Hayas, A. y González-Fernández, C. F., “The effect of previous convective flux on the nonstationary diffusion through membranes”. J. Membr. Sci, 48, 67-77 (1990) Horno, J., García Hernández, M. T., y González-Fernández, C. F., “Digital simulation of electrochemical processes by network approach”. J. Electroanal. Chem., 352, 83-97 (1993) Horno, J., García Hernández, M. T., y González-Fernández, C. F., “Digital simulation of electrochemical processes by network approach”. J. Electroanal. Chem., 352, 83-97 (1993) Incropera, F.P. y DeWitt, D.P. Fundamentals of heat and mass transfer. John Wiley & Sons, New York (1996) Irey, R.K., “Errors in the One-Dimensional Fin Solution”. J. Heat Trans-T. ASME, 90, 175-176 (1968) 171

Kakaç, S. y Yener, Y., “Heat conduction”, 2ª ed. Hemisphere Publishing Corporation, New York (1985) Kalman, H., in Recent advances in analysis of heat transfer for fin type surfaces, ed. B. Sundén and P.J. Heggs, pp. 97-143, Wit Press, Southampton, 2000. Kern, D. Q. y Kraus A. Extended surface heat transfer. McGraw-Hill Book Company. New York (1972) Kraus, A. D., Aziz, A. and Welty, J. Extended Surface Heat Transfer, John Wiley and Sons, New York (2001). Kreith, F. Y Romie, F. E., “A study of the thermal diffusion equation with boundary condition corresponding to solidification or melting of materials initially at the fusion temperature”. Proc. Phys. Soc. Lond. 68(B), 383, (1955) Liebmann, G., “Resistance-network analogues with unequal meshes or subdivides meshes”. Brirish J. Appl. Phys., 5, 362-266 (1954a) Liebmann, G., “Note on the resistance network analogue solution of field problems of spherical simmetry”.Brirish J. Appl. Phys, 5, 412 (1954b) Liebmann, G., “A new electrical analog method for the solution of transient heat conduction problems”. Trans. ASME, 78, 655-665 (1956a) Liebmann, G., “Solution of transient heat-transfer problems by the resistance-network analog method”. Trans. ASME, 78, 1267-1272 (1956b) Lienhard J. H., “A heat transfer textbook”. 2º ed. Prentice-Hall Inc. New Jersey (1987) Lin, S.H., “Transient heat conduction in a composite slab with variable thermal conductivity”. Int. J. Numer. Methods Eng., 14, 1726-1731 (1979) López-García, J. J., Moya, A. A., Horno, J., Delgado, A. Y González-Caballero, F., “A network model of the electrical double layer around a colloid particle”. J Colloid Interface Sci., 183, 124-130 (1996) López Sánchez, J. F., Alhama, F. and González-Fernández, C.F. Introduction and performance of species in a diffusive Lotka-Volterra system with time dependent coefficients. Ecological Modelling 183, 1-9 (2005) MATLAB 6. MathWorks, Natick, MA (1997) Mathematica, Wolfram Research (UK) Ltd., Abington, Oxon OX13 6TG. Reino Unido Mills, A. F., “Heat and mass transfer”. Richard D. Irwin, Inc., Chicago (1995) Moreno Nicolás, J.A., Hernández Grau, J. y Alhama, F. Modelo de cadena cinemática de un motor diesel de media velocidad. Anales de Ingeniería Mecánica, Vol. 3, 1685-1691 (2004)

172

Moreno Nicolás, J. A., Gómez de León Hyjes y F. Alhama. Solution of temperature fields in hydrodynamics bearing by the numerical network model. Tribology International 40, 139145 (2007) Nagel, L.W., “SPICE2: A computer program to simulate semiconductor circuits”. Memo. Nº. UCB/ERL M520. Electronic Research Laboratory, Univ. de California, Berkeley, CA 94720 (1975). Nagel, L.W., “SPICE (Simulation program with integrated circuit emphasis)”, Berkeley, CA; University of California, Electronics Res. Lab., ERL-M382, (1977) Oppenhein, A. K., “Radiation analysis by the network Method”. Trans. ASME, 78, 725 (1956) Orivuori, S., “Efficient method for solution of nonlinear heat conduction problems”. Int. J. Numer. Methods Eng., 14, 1461-1476 (1979) Otis, D. R. “Solving the melting problem using the eletric analog to heat conduction”.Heat Transfer and Fluid Mechanics Institute, Stanford Univ., (1956) Özisik, M.N., “Heat conduction”, 2ª ed. John Wiley & Sons, Inc., New York (1993) Palmieri, J. V. y Rathjem, K. A., “CAVE 3-A General transient heat transfer computer code utilizing eigenvectors and eigenvalues”. NASA Contract. Rept. 145290 (1978) Patankar, S. V., “A numerical method for conduction in composite materials, flow in irregular geometries and conjugate heat transfer”.Proc. 6th Int. Heat Transfer Conf. Toronto, 3, 297 (1978) Patankar, S. V., “Numerical heat transfer and fluid flow”. Hemisphere Publishing Corporation, New York (1980) Peusner, L., “Studies in network thermodynamics”, Elseveier, Amsterdam (1986) Peusner, L., “The principles of network thermodynamics: Theory and biophysical applications”, Entropy, Lincoln, Massachusetts (1987) PROCCA-09. Programa de conducción de calor, 2009. Del Cerro Velázquez y col. En fase de Registro (2009) Prodasim 1. Programa de diseño de aletas simples. Gómez Lopera, S., Alhama, F. y del Cerro Velázquez, F. NAR: 08 / 2005 / 544 (12/09/2005). © UPCT (2007) PSPICE, versión 6.0: Microsim Corporation, 20 Fairbanks, Irvine, California 92718 (1994) Razelos, P. and Krikkis, R.N., On the optimum thermal design of individual longitudinal fins with rectangular profile, International Communications Heat and Mass Transfer, vol. 30, no. 3, pp. 349-358, 2003. 173

Razelos, P., A note on the heat transfer in convective fins, Wärme und Stoffübertragung, vol 12, pp. 133-119, 1979. Robinson, S, y col. Professional C#, Wrox Press (2001) Rong-Hua Yeh, An analytical study of the optimum dimensions of rectangular fins and cylindrical pin fins, International journal heat mass transfer, vol. 40, no. 15, pp. 3607-3615, 1997. Rukos, E.A., “Continuous elements in the finite element method”. Int. J. Numer. Methods Eng., 12, 11-33 (1978) Seeniraj, V. y Arunachalam, M., “Heat conduction in a semi-infinite solid with variable thermophysical properties”. Int. J. Heat Mass Transf., 22, 1455-1456 (1979) Shih, T. M., “Numerical heat transfer”. Hemisphere Publishing Corporation, New York (1984) Siegel, R. y Howell, J. R., “Thermal radiation heat transfer”. Hemisphere Publishing Corporation, 3º ed., New York (1992) Soto Meca, A., Alhama, F. y González-Fernández, C. F.bAn efficient model for solving density driven groundwater flow problems based on the network simulation method. J. Hidrology 339, 39-53 (2007) Soto, A., Alhama, F. y González-Fernández, C.F. Density-driven flow and solute transport problems. A 2-D numerical model based on the network simulation method. Computer physics communications 177, 720-728 (2007) Taine, J. y Petit, J.P., “Heat transfer”. Prentice Hall International, New York (1993) Thomas, L.C., “Heat transfer”, Caps 3 y 5, Prentice Hall, New Jersey (1992) Varoglu, E. y Liam Finn, W. D., “Finite elements incorporating characteristic for onedimensional diffusion-convection equation”. J. Comput. Phys., 34, 371-389 (1980) Vladimirescu, A., “The spice book”. Joun Wiley & Sons, Inc., New York (1994) Vujanovic, B. D. y Baclic, B. S., “advances in heat transfer”. Academic Press, Inc., San Diego (1994) White, F. M., “Heat and mass transfer”. Addison-Wesley Publishing Company, New York (1988) Wood, W.L., “Control of Crank-Nicolson noise in the numerical solution of the heat conduction equation”. Int. J. Numer. Methods Eng., 11, 1059-1065 (1977) Wood, A. S., Tupholme, G. E., Bhatti, M. I. H., and Heggs, P. J. Performance indicators for steady-state heat transfer through fin assemblies. ASME J. Heat Transfer, 118, pp. 310-316, (1996). 174

Zueco, J., Alhama, F. y del Cerro Velázquez, F. El problema inverso de la determinación del número de Nusselt en convección natural. Solución mediante el método de simulación por redes. V CIDIM, IV CONIM, (Venezuela), Proceedings, 1, 1315-1320 (2001) Zueco, J., Alhama, F. y del Cerro Velázquez, F. El problema inverso de la determinación del número de Nusselt en convección natural. Solución mediante el Método de simulación por redes. Ciencia e Ingeniería 23(3) 37-41, Nov. (2002) Zueco, J. and F. Alhama. Estimación de propiedades térmicas en aleaciones como problema inverso.Revista Española de Metalurgia, 41, 227-232 (2005) Zueco, J. y F. Alhama. Inverse estimation of temperature dependent emissivity of solid metals. Journal of Quantitative Spectroscopy and Radiative Transfer 101, 73-86 (2006) Zueco, J., Alhama, F. and González-Fernández, C.F. Inverse determination of heat generation sources on two dimensional homogeneous solids. Application to orthotropic medium. Int. Comms. Heat Mass Transfer 32, 49-55 (2006a) Zueco, J., F. Alhama y C.F. González Fernández. Inverse determination of temperature dependent thermal conductivity using the network simulation method. Journal of Material Processing and Technology 174, 137-144 (2006b) Zueco, J. y Alhama, F. Simultaneous inverse determination of temperature-dependent thermophysical properties in fluids using the network simulation method. Int. J. Heat Mass Transfer 50, 3234-3243 (2007)

175