Ejemplo Pde Tool

Centro de Estudios Energticos
Departamento de Ingeniera Mecnica

Direccin de Investigacin y Postgrado
UNEXPO Vicerrectorado de Puerto Ordaz
Uso del PDETOOL de MatLab.

Para mostrar el uso del pdetool de MatLab
considere la seccin transversal de la figura,
que corresponde a un horno largo. Los datos y
propiedades son las siguientes:
k1=0,3 W/mK; Cp1=960 J/kgK; 1= 2325 kg/m3;
k2=1,1 W/mK; Cp2=955 J/kgK; 2=2645 kg/m3;
k3=80,0 W/mK; Cp3=447 J/kgK; 1=7850 kg/m3;
L = 2 m; H = 1,5 m; e1 = 0,15 m; e2 = 0,2 m;
s = 0,5 m; x = 1 m; a = 0,05 m; D = 0,06 m;
hi = 800 W/m2K; Ti = 600 C;
he = 20 W/m2K; Te = 40 C;
Te, he
x/2
k2
e2
k1
e1
s
s
e2
e1
Ti, hi
Resistencias,
D, k3
L
Se pide: determinar la distribucin de temperatura a travs de las paredes externa del

horno para condiciones de rgimen estacionario, la temperatura mxima en la superficie
exterior y la perdida de calor por unidad de longitud del horno; si el calor generado en
cada resistencia (por unidad de longitud) es de 3x105 W/m2.


Antes de inicial la solucin del problema se debe simplificar el
dominio, aprovechando la simetra que ste tiene. El dominio
simplificado queda de la siguiente manera (ver Figura)
El procedimiento a seguir para resolver el problemas, usando la
herramienta pdetool del MatLab es el siguiente:
1. En la ventana de comando del matlab se ejecuta pdetool.
2. En la ventana grfica del pdetool se ejecuta especifica el
tipo de anlisis.
Te, he
x/2
k2
k1
3. Se activa y configura el grid.
e1
Ti, hi
L/2
e2


4. Se definen los limites del grid.
4. Se dibuja la geometra.
4.1. Se construye el rectngulo mayor, de dimensiones 1,35x1,85.
1
3
2
4
Nota: Para llamar al Object Dialog, se hace doble
click sobre el rectngulo


4.2. Se dibuja el rectngulo de dimensiones 1,15x1,65.
1
2
3
1.15
1.65
4


4.3. Se dibuja el circulo que representa las resistencia de dimetro 0,06 m y centro
(1,05;0,3).
1
0.5
0.03
0.03
4
Nota: Para llamar al Object Dialog, se hace
doble click sobre la elipse (o circulo)


4.4. Se dibujan los crculos
restantes de dimetro 0,06 m y
cuyos centros son:
a) (1,05;1).
b) (1,55;1).


4.5. Se dibuja el rectngulo interior de dimensiones 1x1,5
1

4.6. En el Set formula se cambia el signo (+)

detrs de la regin R3 a (-), es decir la
operacin debe quedar: (R1+R2+E1+E2+E3)-R3.


5. Se definen las propiedades de los dominios. Para esto se pasa al modo PDE.
1
2
k2
k1
5.1. Se hace click sobre

el dominio de la pared
externa y luego doble
click en el mismo, para
abrir la ventana de
propiedades
(PDE
Specification)
doble
click
k3
1
5.2. En la ventana PDE Specification se selecciona
Eliptic y se fijan los valores:
k = 1.1 (Conductividad); Q = 0 (Generacin de calor);
h = 0 (Coef. Conveccin) y Text = 0 (Temp. del medio)


5.3. Para la pared interna y cada resistencia se debe proceder de forma similar a los
pasos 5.1 y 5.2; pero para los siguientes valores:
- Pared Interna:
k = 0,3; Q = 0; h = 0 y Text = 0
- Resistencias:
k = 80; Q = 1,2e6 ; h = 0 y Text = 0
Pared
interna
Resistencias


Comentario de la ventana PDE Specification.
Los problemas Elpticos (Eliptic), en transferencia
de calor, se refieren a los problemas de conduccin
estacionaria, por tanto la forma general de sta
ecuacin que se resuelve en cada dominio es:
kT Q hText T
O que es lo mismo:
Donde:
T T
k
k
Q hText T
x x x x
div ( k grad (T )) kT
T T
k
k
x x x x
K es la conductividad trmica
Q es la generacin de calor volumtrica para una pared de espesor unitario (Q = qL)
h es el coeficiente de conveccin a travs de la seccin transversal
Text es la temperatura del corriente libre del medio en contacto con las secc.
Transversal.


6. Se definen las condiciones de borde. En este caso se pasa al modo Boundary.
1
2
Simetra
(Aislado)
Conveccin
externa
Conveccin
interna
Aislado


6.1. Se fija la conveccin interna. Se debe seleccionar con el click+shift las lneas (Borde)
del lado interno (se tornan color negro) y luego de hace doble click para abrir la ventana
para definir las condiciones de borde (CB):
Conveccin
interna
El tipo de CB es Neumann y los valores son:

g = h.Text = 800*600; q = h = 800
Al comparar la ecuacin que define la
conveccin con la condicin de borde (CB)
de Neumann, se tiene:
T
hT To
n
o kT n hT To
n.k .grad (T ) kT n
kT n hT hTo g hTo
q h


6.2. Se fija la conveccin externa. Se
procede igual que en el caso anterior,
donde se selecciona Neumann como
tipo de CB y los valores a definir son: g
= 20*40 y q = 20.
6.2. Se define la CB de simetra (es igual que

aislado) y el borde aislado. Se procede igual
que en el caso anterior, donde se selecciona
Neumann como tipo de CB y los valores a
definir son: g = 0 y q = 0.
20*40
20
Nota 1: Cuando se define una CB del tipo Neumann las lneas se tornan de color azul y las del tipo
Dirichlet son de color rojo
Nota 2: Las CB del tipo Dirichlet es cuando se define la
temperatura en ste borde, por lo que las variables
representa los siguientes valores:
h 1
T To
r To


7. Se definen la malla de elementos finitos. En este caso se pasa al modo Boundary.
1
2
Tambin se puede
realizar el mallado
usando este botn
Para refinar la malla se
puede presionar ste
botn y mediante el
men siguiendo la ruta:
Mesh>Redefine Mesh


7.1. Al refinar la malla se obtiene el siguiente modelo de elementos finitos.
Ntese que en las
regiones mas pequeas,
el mallador automtico
refina mas sta regin.


8. Se resuelve el modelo. La ruta es Solve>Solve PDE.
1
2
Tambin se
puede ejecutar el
solve con ste
botn
Nota:
Para los problemas
estacionarios no es
necesario configurar
los parmetros de
solucin en
Parameters


La solucin grfica de la temperatura se muestra despus de terminar el solve.


9. Se configura la salida grfica para mejor visualizacin. La grfica anterior tiene un
contorno de colores poco utilizado, en ese sentido, los contornos de colores se pueden
cambiar usando el panel que sale de: Plot>Parameters (ver figura).
1
2
La ventana grfica o panel es el siguiente, donde en
Colormap se selecciona jet. El resultado grfico es
el siguiente:


Si se activa en Plot type la opcin Contour, y se presiona botn Plot, la grfica resultante
es:


El nmero de contorno automtico se
puede cambiar modificando en
Contour plot levels.
1
2
Los contornos se pueden cambiar

introduciendo un vector con los valores de
los contornos deasados en Contour plot
levels, por ejemplo, se desea mostrar las
isoteremas [80 100 200 300 400 500 600
700 800 900]. Para visualizar las
isotermas se debe desactivar color en
Plot type.


El resultado grfico para los
parmetros grficos anteriores, es
el siguiente.
Si se quiere superponer los
vectores de flujo de calor, se activa
Arrows en Plot type. Quedando el
panel de la siguiente manera:


El resultado grfico para
parmetros anteriores, es:
los
La grfica del flujo de calor se

obtiene cambiando la propiedad
(property) de temperature a heat
flux.
3
4


El resultado grfico del flujo
de calor por unidad de rea
para los parmetros
anteriores, es el siguiente:
Para calcular la temperatura
mxima y el flujo de calor a
travs de la superficie
externa de la pared de
refractario, se debe exportar
la solucin y la malla de
elementos finitos.


La malla se exporta como sigue:
1
En la ventana grfica inicialmente

aparecen la variables: p e t, que
representan las coordenada de los
nodos, los nodos que unen dada lado y
los lados que unen los elementos
triangulares. Las variables se pueden
remplazar, en este caso se han colocado
las variables: x l t
Esta operacin ha enviado al workspace las coordenadas de los nodos en la variable x

(una matriz de tamao 2xN, donde N es el nmero de nodos), los nodos que forman
cada elemento en la variable t (matriz de 4xE, donde E el nmero de elementos y lado:


La solucin se exporta como sigue:
1
En la ventana grfica inicialmente

aparecen la variables: u, que representa
la solucin de la temperatura en cada
nodo. sta variable se pueden
remplazar, en este caso se utiliza Temp
La variable que almacena la solucin (en este caso Temp) es un vector de tamao Nx1,
donde N es el nmero de nodos de la malla de elementos finitos.
Para calcular la temperatura mxima y mnima en la pared exterior se debe determinar
los nodos que estn en sta superficie, por lo que hay dos procedimientos que se
ilustran a continuacin.

Procedimiento 1: Los nodos que estn en los borde se

pueden determinar grficamente activando la numeracin
(o etiqueta) de los nodos, y luego, mediante un zoom se
visualiza la numeracin y se copia en un vector.
Los nodos se muestran en la grficas siguientes:


Procedimiento 2: Grficamente es un poco engorroso definir los nodos en los bordes,
por lo que es ms recomendable usar un pequeo scritp que permita seleccionar stos
nodos. El script puede ser como el que sigue:
%Determinacin de lo nodos en el borde x = 1.35
cont = 0;
for i = 1:max(size(x))
if x(1,i)>1.349&x(1,i)<1.351
cont=cont+1;
ix(cont) = i; %ix almacena los nodos de ste borde
end
end
%Determinacin de lo nodos en el borde y = 1.85
cont = 0;
for i = 1:max(size(x))
if x(2,i)>1.849&x(2,i)<1.851
cont=cont+1;
iy(cont) = i; %iy almacena los nodos de ste borde
end
end
Los nodos de los borde son

(resultado del script):
ix = [1 6 54 55 56 57 58 59 60
61 62 368 369 370 371 372
373 409 410 411 412 413 414
415 416 417 418 536 537 538
539 540 541]
iy = [1 2 44 45 46 47 48 49 50
366 367 395 396 397 398 399
400 401 402 534 535]


Clculo de la Temperatura mxima y mnima: Al scrip anterior se le agregan las
siguiente linea, donde los valores de la temperatura mxima y mnima se amacenan en
Tmax y Tmin, respectivamente:
%Calculo de la mxima temperatura en la pared
Tmax = max(max(Temp(ix)),max(Temp(iy)));
%Calculo de la mnima temperatura en la pared
Tmin = min(min(Temp(ix)),min(Temp(iy)));
Los valores de salida son:

Tmax = 97.4472 C
Tmin = 42.8842 C


Clculo de la Perdida de Calor a travs de la Pared Externa: El calor por unidad de
longitud del horno se calcula por:
T
Q hT T dA k
dA
n
A
A
Como:
dA bds
Donde b es la longitud o profundidad perpendicular al papel y s la longitud de la

arista de la superficie convectiva.
El calor por unidad de longitud del horno es:
Q
T
hT T ds k
ds
b s
n
s
x
Q 1
hT ( x) T dx
Si la cara es paralela al eje x, se tiene:
b xo
y
Q 1
hT ( y ) T dy
Si la cara es paralela al eje y, se tiene:
b yo
Nota: Las integrales

se puede resolver
por el mtodo de los
trapecios, Simpson o
similar.


El script utilizado para calcular las integrales es el siguiente (se debe agregar a los script
anteriores):
%Calcular el flujo de calor por la pared x=1,35 m (Qx).
%Se debe interpolar la temperatura en esta superficie, es decir,
%se crea una funcin: Ty = Temp(y)
yi=0:0.025:1.85; %puntos de integracin
h = 20; %Coeficiente de conveccin
Te = 40; %Temperatura del medio
Tyi=interp1(x(2,ix),Temp(ix),yi); %Interpolacin de Temp
qxi= h*(Tyi-Te); %Flujo de calor por unidad de rea
Qx=trapz(yi,qxi); %Integracin numrica por metodo del trapecio
%Calcular el flujo de calor por la pared y=1,85 m (Qy).
yi=0:0.025:1.35;
Txi=interp1(x(1,iy),Temp(iy),yi);
qyi= h*(Txi-Te); %Flujo de calor por unidad de rea
Qy=trapz(yi,qyi);
%Calculo del flujo de calor total por unidad de longitud del horno
Qt = Qx + Qy;
Los resultados de ste

extracto de scrip, es:
Qx = 1568,4 W/m
Qy = 1004,9 W/m
Qt = 2573,3 W/m


Anlisis Transitorio del problema anterior:
Para comprender el uso del PDETOOL de MatLab en un anlisis trmico, se estudiar el
proceso de calentamiento del horno anterior, considerando que inicialmente se
encuentra a una temperatura uniforme de 35 C. En ste caso, solo se pide determinar
la distribucin de temperatura para varios instante de tiempo a lo largo de dos (2) hora,
y verifique si al cabo de ste tiempo el sistema ha alcanzado la condicin de estado
estacionario.
Los pasos seguidos para resolver el problema en estado estacionario son los mismo
que se debe seguir para resolver el problema transitorio, solo hay que modificar el paso
6, 8 y 9.
El paso 6 se debe modificar para activar la opcin de la ecuacin parablica (Parabolic)
y definir las otras propiedades restantes, como son, densidad y calor especfico.
El cuadro de propiedades para el dominio
de las resistencia queda como sigue:


El cuadro de propiedades para el
dominio de la pared interna es el
siguiente:
El cuadro de propiedades para

dominio de la pared externa es:
el
Una vez definidas las propiedades y el tipo de anlisis (Parabolic) se pasa al mallado de
forma similar a la descrita en el paso 7.
Antes de resolver el problema se deben
definir los parmetros del solve, para lo cual
se sigue: Solve>Parameters.
1
2


Despus de ejecutar la secuencia de comandos anteriores se abre la ventana que
permite definir los parmetros del solve, que para este caso son los siguientes
0:7200
El 0 representa el tiempo inicial y el 7200

(7200s o 2 hr) es el tiempo final. El incremento
de tiempo es 1 s. La estructura general de sta
entrada es: to:dt:tf, donde to es el tiempo
inicial, dt es el incremento de tiempo (cuando
no se especifica se toma como 1) y tf el tiempo
final.
El 35 representa la temperatura inicial, es
decir: T(t=to).
Error relativo de la solucin (se puede reducir

si desea una mayor precisin.
Error absoluto de la solucin (se puede reducir
si desea una mayor precisin.
Despus de configurar los parmetro se

procede a ejecutar el solve (puede se
usando el cono
o mediante el men).
1
2


Para presentar los resultados, primero de cambian los lmites del grid (1) y tambin se
desactiva el grid (2).
(1)
La grfica de temperatura se
presentar para los tiempo
1800s (30 min), 3600s (60 min),
5400s (90 min), 7200s (120 min
2 hr). Para Visualizar la
primera distribucin de
temperatura se configuran los
parmetros del plot (puede
presionar el cono
o por
comando: plot>Parameters
(2)
Aqu se define el vector de isotermas:

[40 50 100 150 200 300 350 400 450 500]
3
4


Para cada grfica se debe cambiar en Time for plot en los valores de tiempo
deseado. En este caso, las grficas para los cuatro tiempos planteados son:



Se puede hacer un video para visualizar la variacin de temperatura en el tiempo, para
esto de debe presional el cono
o llamar al comando: Plot>Parameters sta
ventana se debe configurar como sigue:
Se activa la opcin
Animation (1), luego se
presiona Options(2) y
en la venta de opcin de
animacin (3) se configura
la velocidad (Animation
rate) y el nmero de
repeticiones (Number of
repeats), se presiona OK
(4) y luego Plot (5) para
ver el video de solucin.
2
4
5


Para hacer un estudio mas detallado del proceso de calentamiento del horno y verificar si
se alcanza la condicin de estado estacionario, es necesario exportar la malla y la
solucin al Workspace. El procedimiento es similar al descrito anteriormente y las
ventanas de variables son las siguientes:
Variables exportadas en la malla.
xx coordenadas de los nodos.
Variables exportadas de la Solucin.

Temp es la temperatura en los nodos.
Temp es una matriz de
dimensiones NxM, donde N
es el nmero de nodos y M en
nmero de divisiones del
intervalo de tiempo +1.
1
Para verificar los nodos a estudiar se activa la
numeracin de nodos en la malla y se visualiza el
nodo. La malla de elementos finitos es la que se
muestra a continuacin.
2


De la malla se obtiene que los nodos a analizar son: 1, 5, 7, 9, 176, 876. Donde se
espera que el nodo mas fro es el 1.


Para construir las grficas y estudiar el calentamiento de los nodos anteriormente citados,
se debe desarrolla el siguiente script.
%Grfica de temperatura de nodos ms importantes
ti=0:1:7200; %divisin del intervalo de tiempo usado en el pdetool
%Los nodos a estudiar son: 1, 5, 7, 9, 176, 876
T1=Temp(1,:);
T5=Temp(5,:);
T7=Temp(7,:);
T9=Temp(9,:);
T176=Temp(176,:);
T876=Temp(876,:);
plot(ti,T1,'b-',ti,T5,'k-',ti,T7,'c-',ti,T9,'m-',ti,T176,'y-',ti,T876,'r-');
xlabel('tiempo (s)'); ylabel('Temp (C)');
legend('Nodo 1','Nodo 5','Nodo 7','Nodo 9','Nodo 176','Nodo 876','Location','Best');
%grfica del nodo ms fro, Nodo 1.
figure(2);
plot(ti,T1,'b-');
xlabel('tiempo (s)'); ylabel('Temp (C)');
legend('Nodo 1','Location','Best');
La grfica de salida es la siguiente:


600
39.5
39
500
38.5
38
Temp (C)
Temp (C)
400
300
Nodo 1
Nodo 5
Nodo 7
Nodo 9
Nodo 176
Nodo 876
200
100
37.5
37
36.5
36
35.5
35
1000
2000
3000
4000
tiempo (s)
5000
6000
7000
Nodo 1
8000
1000
2000
3000
4000
5000
tiempo (s)
6000
7000
8000
Se observa que los nodos en la superficie interna de la pared alcanzan prcticamente el

estado estacionario, pero el nodo ms fro y el del centro de la resistencia no han
alcanzado el estado estacionario, por lo que el sistema todava no ha alcanzado la
condicin de equilibrio. En general para que el sistema alcance el estado estacionario,
todos los nodos deben alcanzar el estado estacionario.

Ejemplo Pde Tool

Cargado por

Copyright:

Formatos disponibles

Ejemplo Pde Tool

Cargado por

Información del documento

Descripción original:

Derechos de autor

Formatos disponibles

Compartir este documento

Compartir o incrustar documentos

Opciones para compartir

¿Le pareció útil este documento?

¿Este contenido es inapropiado?

Copyright:

Formatos disponibles

Ejemplo Pde Tool

Cargado por

Copyright:

Formatos disponibles

Centro de Estudios Energticos

Departamento de Ingeniera Mecnica

Uso del PDETOOL de MatLab.

Se pide: determinar la distribucin de temperatura a travs de las paredes externa del

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

3. Se activa y configura el grid.

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

Nota: Para llamar al Object Dialog, se hace doble

4.6. En el Set formula se cambia el signo (+)

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

5.1. Se hace click sobre

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

El tipo de CB es Neumann y los valores son:

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

6.2. Se define la CB de simetra (es igual que

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

Los contornos se pueden cambiar

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

La grfica del flujo de calor se

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

Centro de Estudios Energticos

Uso del PDETOOL de MatLab.

En la ventana grfica inicialmente