Introducción a la resolución de problemas de Transmisión de Calor con MATLAB (2).pdf

Introducción a la resolución de problemas de Transmisión de Calor mediante MATLAB Jaime López Gutiérrez, 2012

Tabla de Contenidos Tabla de Contenidos Introducción Cálculo del calor intercambiado por radiación entre un gas y una pared con qrad.m Cálculo de las propiedades convectivas para agua en ebullición en película con filmboil.m Cálculo del coeficiente de película de convección-radiación para ebullición en película con htotal.m Cálculo de absortividades de gases participativos con absort.m Cálculo de funciones de radiación del cuerpo negro con planck.m Cálculo de emisividades/absortividades medias con emismedia.m Resolución paramétrica de un problema completo con ayuda de MATLAB

Introducción a la resolución de problemas de Transmisión de Calor mediante mediante MATLAB Página 2 de 22

2 3 5 7 9 12 14 17 18

Introducción El presente documento pretende servir de ayuda a los alumnos que cursen la asignatura de “Ampliación de Transmisión de Calor” del Plan de Estudios de Ingeniero Industrial (2000), o sus equivalentes en los Grados y Másteres adaptados al Espacio Europeo de Educación superior en fase de implantación, a la hora de enfrentarse a la resolución académica y/o profesional de problemas ingenieriles en el ámbito de la Transmisión de Calor.

Dado que se trata de una disciplina eminentemente aplicada, y en la cual es frecuente encontrarse con problemas parcialmente indeterminados, los procesos de resolución iterativa cobran especial importancia a la hora de abordar numerosos problemas tanto en el ámbito académico como en el profesional. Los cálculos de eficiencia de cambiadores de calor mediante el método de determinación de la tempereatura logarítmica media o los balances convectivo-radiantes en sistemas con fuentes o sumideros de calor son buenos ejemplos de ello. Además, muchos de los cálculos que se requieren para la resolución son complejos e intensivos en tiempo, recurso extremadamente valioso, por ejemplo, en un examen. Por otra parte, al programar y automatizar cálculos de antemano, se evita en gran medida el error humano de cálculo, fuente muy frecuente de resultados erróneos.

El progreso tecnológico ligado a los ordenadores personales, cada vez más potentes, portables y autónomos, los hace idóneos para la resolución de este tipo de problemas, ya que los procesos iterativos y de cálculo numérico son fácilmente automatizables en cualquier lenguaje de programación de alto nivel, aligerando notablemnte la carga de trabajo del estudiante o ingeniero al ahorrarle tediosos cálculos repetitivos, así como mejorando notablemente la precisión y el tiempo de cálculo.

Es por ello que se ha concebido una serie de scripts en el entorno de programación MATLAB que acompañan a este documento, y que pretenden ser herramientas de apoyo para la asimilación de la materia y la resolución de problemas, pero no sustituyen en absoluto a la necesaria comprensión de la asignatura. Cabe destacar que, hasta donde llega el conocimiento del autor, no existe aún programa que permita interpretar y plantear un

problema térmico, sino solamente resolverlo una vez identificadas las ecuaciones que rigen su

comportamiento, para lo cual es imprescindible el saber hacer de un ingeniero.

Introducción a la resolución de problemas de Transmisión de Calor mediante MATLAB Página 3 de 22

Los scripts han sido creados desde cero por el autor a partir de los contenidos de la asignatura, no estando sujetos a ningún tipo de propiedad intelectual. Por ello, pueden ser modificados o alterados en la forma que el usuario considere oportuna. Se distribuyen bajo licencia Creative Commons BY-NC-SA, lo que implica que las únicas limitación que se impone a su uso y redistribución son que se cite al autor, que no se comercialicen con fines lucrativos y que su redistr ibución sea bajo esta misma licencia.

Esta obra está licenciada bajo una Licencia Creative Commons Atribución-NoComercial-CompartirIgual 3.0 Unported. Para ver una copia de esta licencia, visita http://creativecommons.org/licenses/bync-sa/3.0/ o envía una carta a Creative Commons, 444 Castro Street, Suite 900, Mountain View, California, 94041, USA.


Cálculo del calor intercambiado por radiación entre un gas y una pared con qrad.m

Iniciaremos la explicación con un ejemplo básico de la aplicabilidad de la computación a la resolución de problemas térmicos. En el caso que nos ocupa, hemos de calcular el calor intercambiado en régimen permanente por radiación entre un determinado gas, no necesariamente gris, y una pared o recinto con la que está en contacto.

Para ello, haremos uso de la fórmula 

=

·( ·

ସ ௚

−

·

ସ ௣)

Que establece una diferencia entre la cantidad de calor absorbido por el gas (emitido por la pared a calor emitido por el mismo gas a

௚

௣

) y el

. En caso de que el gas no sea gris, la absortividad se debe evaluar a la

temperatura de la pared, de donde proviene la radiación a absorber.

Además, calcularemos también el coeficiente de película radiativo, especialmente útil en problemas en los que la radiación y la convección están acopladas de manera no lineal y se hace necesario calcular un coeficiente de película de convección-radiación,

procedimiento que se verá en apartados siguientes. Dicho coeficiente se

define como  

ℎ =

(

௚

−

௣)

Para efectuar el cálculo, utilizaremos el primero de los scripts, llamado qrad.m. Su código se presenta a continuación: function [qr,hr] = qrad(epsilon, alpha, Tg, Tp) %Calcula el calor intercambiado en R.P. por radiación entre %un gas no necesariamente gris y una 'pared' negra. %Introducir temperaturas en K %Jaime López, 2012 sigma = 5.67*10^-8; qr = sigma*(epsilon*Tg^4 - alpha*Tp^4); hr = qr/(Tg-Tp); end


Como es obvio, los parámetros de llamada a la función son la absortividad, la emisividad, y las temperaturas de pared y gas; y los parámetros de salida son, por este orden, el calor intercambiado por radiación y el ya mencionado coeficiente de película radiativo.

Este script, en concreto, no es especialmente útil por sí solo, pero sí lo es en combinación con otras rutinas que se detallarán más adelante. Por otra parte, es un ejemplo didáctico en cuanto a la estructura de programación y la sintaxis de MATLAB. Veamos ahora brevemente cómo utilizarlo.

Primeramente, debemos asegurarnos de que todos y cada uno de los programas a ejecutar se encuentran en el Path de MATLAB. Para ello, es conveniente buscar la carpeta que los contiene en la ventana Current Folder ,

y una vez localizados, indicarle a MATLAB que pretendemos hacer uso de ella mediante la opción

“ Add to path > Selected folders and subfolders” . A partir de ahora, el programa reconocerá las llamadas que hagamos a dichos programas ya que ya sabe dónde se encuentran.

૙ ૞

E JEMPLO

= .

Calcular el calor intercambiado por radiación por un gas no-gris ( caliente (

=

૞૙૙

) que entra en un recinto más frío (

=

૜૙૙

≠

૙ ૜૞

= .

)

).

Para efectuar el cálculo, pasaremos a llamar a nuestra función en el Command Window . Para ello, tecleamos >> [qr,hr] = qrad(0.35, 0.5, 500, 300)

Y obtenemos por pantalla qr = 1.0107e+003

hr = 5.0534

Lo cual quiere decir que los resultados que buscábamos son



= 1010.7

y ℎ = 5.5034


/

ଶ

.

Cálculo de las propiedades convectivas para agua en ebullición en película a presión atmosférica con filmboil.m Uno de los temas de la asignatura está dedicado al estudio detallado de la ebullición en película, tomándose como sustancia de referencia para dicho estudio el agua a presión atmosférica. Uno de los problemas que se abordan consiste en el cálculo de las propiedades convectivas de dicha sustancia, en particular de su numero de Nusselt

y posteriormente de su coeficiente de película, parámetros de vital importancia para estudiar la

transmisión de calor por convección entre el agua y cualesquiera cuerpos o fuentes caloríficas que estén provocando dicha ebullición en película.

Para llevar a cabo este cálculo, se utilizan correlaciones empíricas que utilizan las propiedades termofísicas del agua en las condiciones del problema para dar un valor del número de Nusselt, que caracterizará la transmisión convectiva en dichas condiciones. En dos problemas distintos, a igualdad de números de Nusselt, podemos afirmar que existe semejanza y por tanto la transferencia convectiva será cualitativamente equivalente entre ambos.

Una de las correlaciones de mayor aplicación para el cálculo del nñumero de Nusselt es la siguiente:

í í

í

= 0.62 · é



ቈሺ ቈሺ

ሻ ሺ ሻ ሺ −

௟

௩ ௩

= 0.67 ·

௟

௩

௦ −

−

௩ ௩

ᇱ ℎ௙௚

௩

௦

௦௔௧

ᇱ ℎ௙௚

−

ሻ቉ ሻ቉

ଷ ଵ/ସ

௦௔௧

ଷ ଵ/ସ

Y una vez disponemos de dicho número adimensional, el coeficiente de película convectivo se calcula fácilmente como ℎ஼௏ =

·

௩

La función filmboil.m pretende simplificar este engorroso cálculo, frecuente fuente de errores numéricos cuando se abordan problemas de este tipo. Por ello, se ha estructurado la entrada de datos de manera similar a las tablas de propiedades termofísicas del agua que presentan Incropera et al. en las distintas


ediciones de su libro Fundamentals of Heat and Mass Transfer . De este modo, se reduce la posibilidad de equivocarse en el orden de magnitud o las unidades de los datos.

El código de la función en cuestión es el siguiente:

function [Nu,hconv] = filmboil(D,vv,kv,Cpv,Ts,geom) %Calcula el hconv para ebullición en película sólo para AGUA a Patm % %Entradas: (tal y como vienen en tablas Incropera) %%Recordar buscar propiedades a Tmdvap = (Tp+Tg)/2 ! % D en m % vv*10^6 en m^2/s % kv*10^3 en W/mK % Cpv en kJ/kgK % Ts en ºC %Jaime López, 2012 format bank Tsat = 100; hfgp=2257+.8*Cpv*(Ts-Tsat); if geom == 'cil' Nu = .62*(((9.81*10^6*hfgp*D^3)/(10^-9*vv*kv*(Ts-Tsat))))^.25; hconv = Nu*kv*10^-3/D; elseif geom == 'esf' Nu = .67*(((9.81*10^6*hfgp*D^3)/(vv*kv*(Ts-Tsat))))^.25; hconv = Nu*kv*10^-3/D; end


Cálculo del coeficiente de película combinado de convección-radiación para ebullición en película con htotal.m Después de haber visto en detalle dos rutinas que permiten calcular parámetros necesarios para la resolución de problemas sencillos de ebullición en película, pasamos ahora a abordar el caso completo en el cual debemos considerar que la radiación y la convección no son independientes y por tanto se hace necesario calcular el coeficiente de película combinado de convección-radiación ℎ para el problema.

Este acoplamiento no es en absoluto lineal, siendo una de las expresiones más aceptadas para su cálculo la que propone Bromley , que consiste en una ecuación trascendente en ℎ cuya resolución nos da el valor buscado: ℎସ/ଷ = (ℎ௖௩ )ସ/ଷ + ℎ · ℎଵ/ଷ



Por otra parte, el mismo autor propone que si ℎ ≪ ℎ௖௩ , la solución de la ecuación se puede aproximar por 3 ℎ = ℎ௖௩ + ℎ 4

Dado que la resolución manual de la ecuación trascendente requiere de una aproximación iterativa, podemos optimizar tiempo y esfuerzo mediante el cálculo computacional. Para ello se ha creado el programa htotal.m, que hará uso de los dos programas anteriores para el cálculo del coeficiente de película buscado.

Primeramente deberemos utilizar filmboil.m para obtener el coeficiente de película convectivo a partir de los datos del problema, dato que necesitaremos a la entrada de htotal.m. Se ha preferido desacoplar estos dos programas debido a que, de lo contrario, htotal.m necesitaría un número excesivo de parámetros de entrada (todas las propiedades termofísicas del agua más las propiedades radiativas), lo que resultaría incómodo. Una vez dispongamos de ℎ௖௩ lo introduciremos como parámetro de htotal.m, cuyo código se muestra a continuación:


function [h]=htotal(hconv, epsilon, alpha, Ts, Tsat) %Calcula el h combinado de convección-radiación %Introducir temperaturas en ºC!!!! %Jaime López, 2012 format bank Ts=Ts+273; Tsat=Tsat+273; %Cálculo del h de sólo radiación mediante qrad.m [~,hrad] = qrad(epsilon, alpha, Tsat, Ts); %Valor semilla de h para el inicio del cálculo hprev = hconv+.75*hrad; %Aproximación a la resolución de la ecuación trascendente para casos con %radiación poco significativa if hrad<.5*hconv h = hprev; else %Definición de la ecuación trascendente como Function Handle s= @(h)h.^1.3333 - hconv^1.3333 + hrad*h.^.3333; %Resolución de la ecuación trascendente en el intervalo de 0%-150% %del valor semilla (realmente, cálculo de ceros de la función s) [x,fval] = fzero(s,[0 1.5*hprev]); %Presentación de posibles soluciones h=x; end

Como queda patente en este caso, es posible (y frecuente, y muy recomendable) efectuar llamadas a una función dentro de otra. Este recurso constituye el elemento constructivo básico de la programación estructurada. De hecho, es posible llamar a una función y obviar alguno de sus argumentos de entrada o, más frecuentemente, de salida. Es el caso de la llamada a qrad.m, donde hemos sustituido con una tilde (~) el lugar de la variable de salida qrad que no nos interesa en este momento, por lo que el programa la obviará en ejecución y no la calculará explícitamente. Además, es interesante conocer la manera en que MATLAB trata a las funciones de una o varias variables, en su forma Function Handle. Es muy conveniente definir así las funciones de cara a operar numéricamente sobre ellas (derivación, integración, cálculo de ceros) en vez de simbólicamente o mediante vectores, opciones que nos permiten un rango de maniobra mucho más limitado. En la documentación de MATLAB existe ayuda extensa sobre la sintaxis de los Function Handles y su utilización en distitnas operaciones, la cual es muy recomendable consultar.


ଵ

En nuestro caso, hemos implementado la condición de radiación poco significativa como ℎ < ଶ ℎ௖௩ , caso en el cual daremos como solución la aproximada propuesta por Bromley . Para solucionar ecuaciones trascendentes en MATLAB, el procedimiento más adecuado consiste en despejar a un miembro todas las variables e igualarlas a una función auxiliar (que en nuestro caso hemos llamado ) y cuyos ceros buscaremos mediante el comando fzero en un intrervalo de nuestra elección. Para determinar este intervalo, aprovecharemos como valor semilla el valor aproximado de ℎ dado por Bromley para el caso de radiación poco significativa, que ya hemos tenido que calcular previamente. Parece razonable que si la radiación sí contribuye de manera importante a nuestro problema, el ℎ final se encontrará previsiblemente entre 0 y un 150% del valor dado por la aproximación, con lo cual determinaremos el rango de búsqueda de ceros de la función auxiliar .

E JEMPLO

(P ROBLEMA 1 ATC F EBRERO 2010)

Sea un cilindro de 50cm de diámetro y 10m de longitud introducido verticalmente a 400ºC en un gran tanque de agua saturada a presión atmosférica. El cilindro mantiene constante su temperatura superficial en todo momento por medio de una fuente térmica en su interior. Suponiendo que el coeficiente de ebullición en película vale / ૛ , calcular el calor total que pierde el cilindro por convección y radiación combinadas y su coeficiente de película de convección-radiación.

ૡ૙

>> [h]=htotal(80, 1, 1, 400, 100) hrad = 35.11

h = 106.34

Lo cual es un resultado incluso más preciso que el proporcionado en la resolución oficial, no habiendo sido necesaria la tediosa iteración manual.


Cálculo de absortividades de gases participativos con absort.m

Otra parte importante de la asignatura se dedica al estudio de la transferencia radiante en la que intervienen gases participativos, esto es, que no son transparentes a la radiación y por tanto participan del intercambio tanto absorbiendo como emitiendo radiación en función de su temperatura, la de las superficies con las que están en contacto y sus características inherentes, que dependerán del gas a estudiar en cada caso.

Principalmente, el estudio se centra en los dos gases con los que más frecuentemente se topará un ingeniero en el ejercicio de su profesión, el ଶ

ଶ

(como producto de combustión de hidrocarburos o carbones) y el

en su forma de vapor, resultante de la vaporización de agua líquida. Dado que el estudio de estos gases

tiene un fuerte carácter empírico, se recurre a correlaciones y gráficos para multitud de cálculos a lo largo del capítulo. Uno de los casos particulares es el del cálculo de la absortividad de dichos gases, el cual se abordará en este apartado.

Si actuamos según el procedimiento propuesto por J.P. Holman, para calcular la absortividad de un gas en unas determinadas condiciones, recurriremos a los gráficos de H.C. Hottel (eligiendo el gráfico convenientemente en función de la sustancia a estudio) para primero calcular su emisividad. Entraremos con la temperatura y el factor de un

௚′

௖

·

௘

y de ahí leeremos el

௚

. Ahora, el procedimiento requiere el cálculo

que se obtiene también del gráfico, pero utilizando como datos de entrada la temperatura y, en

esta ocasión, un factor

௖

Una vez disponemos de

·

௘

௚

y

·(

ᇱ ௚

௣/ ௚ )

.

, estamos en condiciones de aplicar las ecuaciones (8.55) y (8.56) de

Holman para calcular la absortividad, dependiendo del gas con el que tratemos, de la manera que sigue

ቆቇ ቆቇ

଴.଺ହ

஼

=

ᇱ ஼

௚

௣

଴.ସହ

ௐ

=

ᇱ ௐ

௚

௣


Para facilitar el cálculo y posibilitar su incorporación en programas de mayor entidad, se ha elaborado la función absort.m, cuyo sencillo código es el siguiente:

function [alpha]=absort(Tg,Tp,epsprima,subs) %Calcula la absortividad de gases participativos (CO2, H2O) %Introducir temperaturas en ºC!!!! %Si no, comentar líneas 10 y 11 %e' se evalúa a (Tp,pc·Le·(Tp/Tg)) %Jaime López, 2012 format bank % Paso a K Tg = Tg+273; Tp = Tp+273; if subs=='CO2' alpha = epsprima*((Tg/Tp)^.65); elseif subs=='W' alpha = epsprima*((Tg/Tp)^.45); end end


Cálculo de funciones de radiación del cuerpo negro con planck.m

Uno de los cálculos más frecuentes cuando se abordan en el estudio de la radiación cuando se trata con cuerpos reales, o al menos no grises, es el de las funciones de radiación del cuerpo negro tabuladas como función de

଴→ఒ

. Éstas están

entre cero y uno, pero es frecuente desconocer alguno de estos datos (por

ejemplo, la temperatura del cuerpo) y tener que definir un proceso de iteración para poder hallar el facto de banda que satisface el balance energético a resolver, lo cual nos proporcionará la temperatura solución del problema.

Dado a que esto es tremendamente trabajoso de manera manual, se ha programado una calculadora de funciones de radiación para simplificar el trabajo y posibilitar la integración del cálculo de factores de

banda en

problemas más complejos, a fin de automatizar la iteración cuando se requiere para resolver el problema.

Para ello, se ha partido de la definición de las funciones de radiación en la ecuación (12.28) de Incropera el al., ఒ ଴ →ఒ

≡

ఒ

଴ ஶ

ఒ,௕

଴

ఒ,௕

=

଴

ఒ ,௕

=

ସ

න

ఒ்

ఒ,௕

ሺ ሻ න ሺ ሻሺ ሻ ఒ்

=

ହ

଴

଴

Que a su vez requiere definir la potencia emisiva total del cuerpo negro

௕

con arreglo a la sustitución de la

ecuación (12.24) en la (12.9) del mismo texto,

න ቂ ቀቁ ቃ ஶ

௕

=

ଵ

ହ

଴

ଶ

exp

−1

Para implementar el cálculo de la integral que da el valor de la función de radiación con una precisión aceptable no es suficiente el uso de los valores de las constantes de radiación primera y segunda propuestos en el texto, ya que contienen redondeos y conducen a errores numéricos. Para evitar esto, ha sido necesario definir las constantes de radiación a partir de su relación con las constantes fundamentales ℎ, , dándole a estas últimas un valor lo más exacto posible. ଵ

=2 ℎ

ଶ ଴

  ;

ଶ

=

ℎ

଴


଴

,

Presentamos a continuación el código de planck.m, sobre el que efectuaremos después algunas aclaraciones:

function F = planck(lambda1, T, lambda2) %Este script calcula las funciones de banda. %Introducir T en K, lambda en um %Jaime López, 2012 format long c = 299792458*10^6; h = 6.6260695729*10^-34; k = 1.380648813*10^-23; sigma = 5.67037321*10^-8; C1 = 10^12*2*pi*h*c^2; C2 = h*c/k; LT1 = lambda1*T; %Integración trapezoidal % LT=1:1:LT1; % Eb = C1./(sigma.*(LT.^5).*((exp(C2./LT))-1)); % F1 = trapz(LT,Eb); %Integración por Cuadratura de Lobatto Eb = @(LT) C1./(sigma.*(LT.^5).*((exp(C2./LT))-1)); F1 = quadl(Eb,0,LT1,10^-8); if nargin>2 LT2 = lambda2*T; F2 = quadl(Eb,0,LT2,10^-8); F = (F2-F1); else F = F1; end end


El código del programa planck.m permite elegir entre dos tipos de integración: t rapezoidal (más rápida e imprecisa) o por cuadratura de Lobatto (más lenta y precisa). Comentando o descomentando las líneas de selección a voluntad, el usuario puede decicir cuál utilizar en cada momento.

Además, cabe reseñar que hemos utilizado de nuevo los Function Handles como vehículo para definir como variable de trabajo, en vez de sólo o T por separado. Así, la función que realmente integra el programa es la que hemos dado en llamar Eb = @(LT) C1./(sigma.*(LT.^5).*((exp(C2./LT))-1));

Que no es más que la expresión Function Handle del integrando

ሺ ሻ ሺ ሻቂ ቀ ቁ ቃ ଵ

=

ହ

ଶ

exp

−1

Por último, indicar que existe una condicional al final del programa if nargin>2

Que en caso de que la llamada al programa se haya realizado con sólo dos argumentos ( ଵ , ), nos dará directamente el valor de la función de radiación

଴→ఒభ

como el valor de la integral calculada. Sin embargo,

si se ha introducido un tercer argumento ( ଶ ), calculará también la integral que da lugar a resultado del problema será el factor de banda ఒభ →ఒమ

ሺ ሻ ଵ,

ଶ,

=

଴→ఒమ

ሺ ሻ ଶ

−

଴→ఒభ ( ଵ

଴→ఒమ

y el

)

Lo cual resultará muy útil en problemas más complejos donde debemos considerar múltiples frecuencias de corte. E JEMPLO

ૠ૙૙

Calcular la fracción de radiación que el Sol ( )

~

૞ૡ૙૙

) emite en el visible (

>> F = planck(0.4, 5800, 0.8) F = 0.46109541343799

(Resultado que podemos comprobar utilizando las tablas de funciones de radiación)


૝૙૙ ૃ ≤

≤

Cálculo de emisividades/absortividades medias con

emismedia.m

Una de las razones que llevó al desarrollo de planck.m fue crear la posibilidad de calcular automáticamente emisividades y absortividades espectrales totales a partir de datos espectrales y de temperatura, lo cual puede ser muy útil a la hora de obtener una primera aproximación a la solución de problemas complejos.

Con este programa podemos evitarnos el tener que computar por separado las integrales de factores de banda o buscarlas en tablas, ya que Planck.m se encarga de evaluarlas numéricamente por nosotros. function e = emismedia(T, e1, lambda1, e2, lambda2, e3, lambda3) %Este script calcula emisividades y absortividades medias a partir de datos %de emisividad/absortividad espectral por bandas, usando el script de las %funciones de banda planck.m %ei son las emisividades de cada banda, y lambdai las freq de corte %SUPERIORES de cada banda %Jaime López, 2012 p1 = planck(lambda1, T); if nargin <=4 %Si sólo hay dos bandas, limitadas por 0 e inf e = e1*p1 + e2*(1-p1); else %Si hay tres bandas p2 = planck(lambda1, T, lambda2); p3 = planck(lambda2, T, lambda3); e = e1*p1 + e2*p2 + e3*p3; end format short end

E JEMPLO

(I NCROPERA ET AL . 12.5)

Calcular la emisividad hemisfñerica total de una superficie difusa a 1600K con las propiedades espectrales del gráfico >> e = emismedia(1600, 0.4, 2, 0.8, 5, 0, 9999) e = 0.5578

Lo cual coincide con el resultado del texto.


Resolución paramétrica de un problema completo con ayuda de MATLAB Como aplicación ilustrativa de la utilidad de lo hasta ahora mostrado en este breve manual, pasaremos ahora a resolver un problema completo de la asignatura de Ampliación de Transmisión de Calor propuesto en examen en la convocatoria de Julio 2012.

Todos los problemas, como se indicó en la introducción, requieren indudablemente de una componente humana de planteamiento, imprescindible para su correcta resolución. Una vez planteadas las ecuaciones que darían solución al problema, utilizaremos algunos de los programas hasta ahora citados para resolver paramétricamente, esto es, para hallar la temperatura solución en para todas las combinaciones de frecuencias de corte inferiores y superiores posibles en un rango razonable.


Acorde con la solución oficial publicada, denominaremos a las superficies como sigue

Es claro que la placa sólo participará en el intercambio radiante en la banda donde su reflectividad no es del 100%, según el enunciado para banda como

ଵ

൏ ൏

comprendida entre

ଵ

ଶ

y

. Sin embargo, nosotros parametrizaremos ahora esta

, y operaremos en función de ello. Mantendremos, en la banda no rerradiante

ଶ

distribución y el valor de

5൏ ൏8 ൌ ൌ 0.5 ൌ

,

ఒ

ఒ

ఒ

ఒ

. Del mismo modo, también mantendremos la

dados en el enunciado original.

Siguiendo la resolución oficial, calcularemos el calor absorbido por la placa como ஶ

ൌන ൑ ൑ 9

ఒమ

ൌන

஺஻

ఒ

ఒ

଴

Siempre y cuando

ଵ

ଶ

ఒ

ൌ

ఒ

ఒభ

ఒ

ఒ

ሺ െ ሻ ଶ

ଵ

, ya que de lo contrario al menos parte de la radiación emitida no sería

absorbida por la placa en la banda de absortividad no nula.

Efectuamos ahora un balance a la placa,

ൌ ൅ ൅ ቊ ൌ ൌ ൉ ഥ൉ ൉ ൉ሺሺ െെ ሻ ሻ ஺஻ ଵ

஼௩ ଵ௔

஼௩ ଵ௔

ଵ

 ଵଶ

ଵ

஼ௗ ଵ

 ଵଶ

ଵ

௔

ସ ଵ

ଵ

ସ ଷ

Donde la emisividad media se calcula a partir de los factores de banda

ഥ ൌ ൉ൣ ଵ

ఒ

ఒభ →ఒమ

ሺ , , ሻ൧ൌ  ൉ ൣ ሺ ሻ െ ଵ

ଶ

ఒ

଴→ఒమ

ଶ

଴→ఒభ

ሺ ሻ൧ ଵ

Por tanto, el balance quedará ఒ

ఒ

ሺ െ ሻ ൌ ൉ ൉ ሺ െ ሻ ൅ ൉ ൉ൣ ଶ

ଵ

ଵ

ଵ

௔

ଵ

ఒ

ఒభ →ఒమ

ሺ , , ሻ൧൉ ሺ െ ሻ ଵ

ଶ

ଵ


ସ ଵ

ସ ଷ

Ahora, necesitamos plantear un proceso iterativo para calcular la

ଵ

que hace que se cumpla el balance, la

cual será nuestra solución buscada. Para hacer el problema paramétrico, crearemos una serie de vectores de iteración, hipótesis y solución, así como varios bucles anidados que irán evaluando la ecuación del balance a lo largo del vector de temperaturas hipótesis hasta encontrar aquella (o aquellas) que hace que la igualdad se cumpla. Como en este tipo de problemas numéricos lo que realmente se busca es una solución aproximada, daremos un cierto margen al código para que, en vez de la temperatura exacta, nos proporcione aquellas que generan un residuo (diferencia entre el calor de cada miembro de la ecuación) pequeño, para luego hacer una media entre ellas. Las variables de control de los bucles serán la frecuencia de corte inferior i y la frecuencia de corte superior j

(0 ≤ λ ଵ = ; λ ଶ =

≤9

). El vector de hipótesis vendrá representado por

=

(500, 501, … , 1999, 2000), y de dimensión 1 × , y la ecuación del balance será evaluada en sus

componentes siempre y cuando > , o lo que es lo mismo, cuando

ଶ

>

ଵ

, ya que de lo contrario la

superficie de la placa sería rerradiante en todo el espectro, y en régimen permanente estará en equilibrio con el fluido que la baña, esto es, a

ଵ

= 500 .

Como el comando find nos devuelve como resultado las posiciones del vector objetivo donde se cumple la condición dada (y no sus valores), debemos sumarle 500 a ese resultado, ya que nuestro vector de hipótesis comienza en 500K. Una vez obtenidas esas soluciones, las almacenaremos en la matriz

para

posteriormente poder representarlas de forma visual. Presentamos en la página siguiente el código completo del programa, en el que se ha hecho uso de lo explicado a lo largo de este documento, incluyendo explícitamente las llamadas al programa planck.m para el cálculo de factores de banda.


%Resolución paramétrica del Problema 2 Ex ATC Julio 2012 %Jaime López, 2012 %Vector de hipótesis T = 500:1:2000; %Determinación del número de iteraciones n = length(T); p = 90; %Preasignación de matrices de resultados para mejora de rendimiento E = zeros(1,n); B = zeros(p+1,p+1); %Valores de la FCI en um for i=0:p %Valores de la FCS en um for j=0:p %Condición de no reflexión total if j>i %Evaluación de hipótesis for k=1:n %Ecuación de balance E(k)=(.5.*5.67.*10^8.*planck(i+.001,T(k),j+.001).*((T(k)^4)- 500^4) 3.17.*(T(k)-500) -(.5*5000*(j-i))); end

+

%Búsqueda, selección y almacenamiento matricial de soluciones c = find (E>-100&E<100); [sol] = mean(c(:))+500; B(i+1,j+1) = sol; else B(i+1,j+1)=500; end iter=[i,j] end end

Se presenta en la página siguiente, por último, la representación gráfica del resultado, calculado con un mallado de grosor 0.1

para longitud de onda, y de grosor 1

para la temperatura.



Introducción a la resolución de problemas de Transmisión de Calor con MATLAB (2).pdf

Recommend Documents