> restart:#"m11_p11"

Una pared de hormigón de 0,5 m de espesor está en una atmósfera a 300 K y a partir de un cierto instante empieza a recibir una irradiación de 500 W/m2 que absorve en su totalidad. Se pide:
 1.Perfil de temperatura estacionario suponiendo que las pérdidas son por convección con h=10 W.m­2.K­1.
 2.Perfil de temperatura estacionario suponiendo que las pérdidas son por convección con h=10 W.m­2.K­1 y por radiación como cuerpo negro.
 3.Evolución del perfil de temperatura en régimen transitorio.

Datos:

> read`../therm_eq.m`:read`../therm_const.m`:read`../therm_proc.m`:with(therm_proc):

> su:="Hormigon":dat:=[L=0.5*m_,T0=300*K_,E=500*W_/m_^2,h=10*W_/(m_^2*K_)];

[L = `+`(`*`(.5, `*`(m_))), T0 = `+`(`*`(300, `*`(K_))), E = `+`(`/`(`*`(500, `*`(W_)), `*`(`^`(m_, 2)))), h = `+`(`/`(`*`(10, `*`(W_)), `*`(`^`(m_, 2), `*`(K_))))]

Image

Eqs. const.:

> sdat:=[get_sol_data(su)]:dat:=op(dat),op(subs(epsilon=epsilon_,sdat)),Const,SI2,SI1:k=subs(dat,k);

k = `+`(`/`(`*`(1.5, `*`(W_)), `*`(m_, `*`(K_))))

a) Perfil de temperatura estacionario suponiendo que las pérdidas son por convección con h=10 W.m­2.K­1.

Se plantean sendos balances energéticos en las entrefases:

> eq1:=h*(T[1]-T0)+epsilon*sigma*(T[1]^4-T0^4)=k*diff(T(x),x);eq2:=E=k*diff(T(x),x)+h*(T[2]-T0)+epsilon*sigma*(T[2]^4-T0^4);T(x):=T[1]+(T[2]-T[1])*x/L;sol1:=solve(subs(epsilon=0,{eq1,eq2}),{T[1],T[2]});sol1_:=expand(subs(dat,sol1));T[1]=TKC(subs(sol1_,T[1]));T[2]=TKC(subs(sol1_,T[2]));

`+`(`*`(h, `*`(`+`(T[1], `-`(T0)))), `*`(epsilon, `*`(sigma, `*`(`+`(`*`(`^`(T[1], 4)), `-`(`*`(`^`(T0, 4)))))))) = `*`(k, `*`(diff(T(x), x)))
E = `+`(`*`(k, `*`(diff(T(x), x))), `*`(h, `*`(`+`(T[2], `-`(T0)))), `*`(epsilon, `*`(sigma, `*`(`+`(`*`(`^`(T[2], 4)), `-`(`*`(`^`(T0, 4))))))))
`+`(T[1], `/`(`*`(`+`(T[2], `-`(T[1])), `*`(x)), `*`(L)))
{T[1] = `/`(`*`(`+`(`*`(L, `*`(`^`(h, 2), `*`(T0))), `*`(k, `*`(E)), `*`(2, `*`(k, `*`(h, `*`(T0)))))), `*`(h, `*`(`+`(`*`(h, `*`(L)), `*`(2, `*`(k)))))), T[2] = `/`(`*`(`+`(`*`(h, `*`(L, `*`(E))), `*...
{T[1] = `+`(`*`(309.37500000000000000, `*`(K_))), T[2] = `+`(`*`(340.62500000000000000, `*`(K_)))}
T[1] = `+`(`*`(36.22500000000000000, `*`(`C`)))
T[2] = `+`(`*`(67.47500000000000000, `*`(`C`)))

i.e., la cara irradiada se pone a 67 ºC y la otra a 36 ºC (en un ambiente a 27 ºC).

Se puede dibujar a escala todo el perfil de temperaturas si se aproximan las capas límite (e.g. por exponenciales empíricas):

> T01:=subs(sol1_,dat,SI0,T0+(T[1]-T0)*exp(1e1*x));T12:=subs(dat,sol1_,SI0,T(x));T23:=subs(sol1_,dat,SI0,T0+(T[2]-T0)*exp(-1e1*(x-L)));plot({[x,T01,x=-0.5..0],[x,T12,x=0..0.5],[x,T23,x=0.5..1],[x,subs(dat,SI0,T0),x=-0.5..1],subs(dat,SI0,[[0,0],[0,1000]]),subs(dat,SI0,[[L,0],[L,1000]])},x=-0.5..1,'T'=250..350,color=black);

`+`(300, `*`(9.37500000000000000, `*`(exp(`+`(`*`(0.1e2, `*`(x)))))))
`+`(309.37500000000000000, `*`(62.500000000000000000, `*`(x)))
`+`(300, `*`(40.62500000000000000, `*`(exp(`+`(`-`(`*`(0.1e2, `*`(x))), 5.)))))
Plot_2d

b) Perfil de temperatura estacionario suponiendo que las pérdidas son por convección con h=10 W.m­2.K­1 y por radiación como cuerpo negro.

> eq1;eq2;sol1_:=fsolve(subs(epsilon=1,dat,SI0,{eq1,eq2}),{T[1],T[2]}):T1_:=subs(sol1_,T[1])*K_;T[1]=TKC(%);T2_:=subs(sol1_,T[2])*K_;T[2]=TKC(%);

`+`(`*`(h, `*`(`+`(T[1], `-`(T0)))), `*`(epsilon, `*`(sigma, `*`(`+`(`*`(`^`(T[1], 4)), `-`(`*`(`^`(T0, 4)))))))) = `/`(`*`(k, `*`(`+`(T[2], `-`(T[1])))), `*`(L))
E = `+`(`/`(`*`(k, `*`(`+`(T[2], `-`(T[1])))), `*`(L)), `*`(h, `*`(`+`(T[2], `-`(T0)))), `*`(epsilon, `*`(sigma, `*`(`+`(`*`(`^`(T[2], 4)), `-`(`*`(`^`(T0, 4))))))))
`+`(`*`(303.99919540932598259, `*`(K_)))
T[1] = `+`(`*`(30.84919540932598259, `*`(`C`)))
`+`(`*`(325.65769023542461392, `*`(K_)))
T[2] = `+`(`*`(52.50769023542461392, `*`(`C`)))

i.e., la cara irradiada se pone ahora a 53 ºC (en vez de a 67 ºC) y la otra a 31 ºC (en vez de a 36 ºC).

ADICIONAL. Ya puestos, podemos ver el efecto de la emisividad en todo el rango de 0 a 1.

> eqd1:=0=convert(series(subs(T[1]=T0+DT1,T[2]=T0+DT2,rhs(eq1)-lhs(eq1)),DT1=0,2),polynom);DT1_:=solve(eqd1,DT1);eqd2:=0=convert(series(subs(T[1]=T0+DT1,T[2]=T0+DT2,rhs(eq2)-lhs(eq2)),DT2=0,2),polynom);DT2_:=solve(subs(DT1=DT1_,eqd2),DT2);DT2__:=evalf(subs(dat,SI0,k=1.5,DT2_));DT2_e0:=subs(epsilon=0,DT2__)*K_;DT2_e1:=subs(epsilon=1,DT2__)*K_;plot(DT2__,epsilon=0..1,DT2=0..50,color=black);

0 = `+`(`/`(`*`(k, `*`(DT2)), `*`(L)), `*`(`+`(`-`(`/`(`*`(k), `*`(L))), `-`(h), `-`(`*`(4, `*`(epsilon, `*`(sigma, `*`(`^`(T0, 3))))))), `*`(DT1)))
`/`(`*`(k, `*`(DT2)), `*`(`+`(k, `*`(h, `*`(L)), `*`(4, `*`(epsilon, `*`(sigma, `*`(`^`(T0, 3), `*`(L))))))))
0 = `+`(`-`(`/`(`*`(k, `*`(DT1)), `*`(L))), `-`(E), `*`(`+`(h, `*`(4, `*`(epsilon, `*`(sigma, `*`(`^`(T0, 3))))), `/`(`*`(k), `*`(L))), `*`(DT2)))
`/`(`*`(E, `*`(`+`(k, `*`(h, `*`(L)), `*`(4, `*`(epsilon, `*`(sigma, `*`(`^`(T0, 3), `*`(L)))))))), `*`(`+`(`*`(2, `*`(k, `*`(h))), `*`(`^`(h, 2), `*`(L)), `*`(8, `*`(h, `*`(L, `*`(epsilon, `*`(sigma,...
`+`(`/`(`*`(500., `*`(`+`(6.5, `*`(3.061800000000, `*`(epsilon))))), `*`(`+`(80.0, `*`(79.606800000000, `*`(epsilon)), `*`(18.749238480000000000, `*`(`^`(epsilon, 2)))))))
`+`(`*`(40.625000000000000000, `*`(K_)))
`+`(`*`(26.805372224815968002, `*`(K_)))
Plot_2d

i.e. el efecto es casi lineal calentándose 41 ºC si no emite y 27 ºC si emite perfectamente.

c) Evolución del perfil de temperatura en régimen transitorio.

Hay que discretizar y calcular numéricamente.

Sea N el nº de capas (ver Fig. 1) y M el nº de intervalos de tiempo.

El intervalo espacial es Dx=L/N.

El intervalo temporal Dt lo elegiremos para que Fo<1/2; e.g. Fo=1/4.

La discretización es:

> eq11_24;eq_i:=subs(phi=0,h[lat]=0,epsilon[lat]=0,eq11_24_gen_);eq_0:=subs(phi=0,h[lat]=0,epsilon[lat]=0,h0=h,epsilon0=1,E0=0,eq11_24_0_);eq_N:=subs(phi=0,h[lat]=0,epsilon[lat]=0,hN=h,epsilonN=1,EN=E,eq11_24_N_);eqFo:=Fo=a*Dt/Dx^2;eqa:=eq11_5;eqDx:=Dx=L/N;eqDt:=Dt=solve(subs(eqFo,eqa,eqDx,0.25=rhs(eqFo)),Dt);

`/`(`*`(Dm, `*`(c, `*`(DT))), `*`(Dt)) = Qnet
T[i, `+`(j, 1)] = `+`(T[i, j], `*`(Fo, `*`(`+`(T[`+`(i, 1), j], `-`(`*`(2, `*`(T[i, j]))), T[`+`(i, `-`(1)), j]))))
T[0, `+`(j, 1)] = `+`(T[0, j], `*`(2, `*`(Fo, `*`(`+`(T[1, j], `-`(T[0, j]))))), `/`(`*`(2, `*`(Dt, `*`(`+`(`-`(`*`(h, `*`(`+`(T[0, j], `-`(T0))))), `-`(`*`(sigma, `*`(`+`(`*`(`^`(T[0, j], 4)), `-`(`*...
T[N, `+`(j, 1)] = `+`(T[N, j], `*`(2, `*`(Fo, `*`(`+`(T[`+`(N, `-`(1)), j], `-`(T[N, j]))))), `/`(`*`(2, `*`(Dt, `*`(`+`(E, `-`(`*`(h, `*`(`+`(T[N, j], `-`(TN))))), `-`(`*`(sigma, `*`(`+`(`*`(`^`(T[N,...
Fo = `/`(`*`(a, `*`(Dt)), `*`(`^`(Dx, 2)))
a = `/`(`*`(k), `*`(rho, `*`(c)))
Dx = `/`(`*`(L), `*`(N))
Dt = `+`(`/`(`*`(.25000000000000000000, `*`(rho, `*`(c, `*`(`^`(L, 2))))), `*`(k, `*`(`^`(N, 2)))))

Para los valores numéricos, además de los datos, se elige un valor de N del orden de 10 (i.e. aproximar una curva por diez puntos intermedios), de donde se deduce Dx y Dt (para que sea estable). Para llegar hasta el régimen estacionario, o se compara a cada paso la solución para ver cuándo ya no varía apreciablemente, o se estima por otro método el tiempo de relajación, o simplemente se calcula un número grande ve pasos y después se ve si se había estabilizado.

> N:=10;M:=1000;datnum0:=op(subs(SI0,[Dt=subs(dat,rhs(eqDt)),Dx=subs(dat,rhs(eqDx)),Fo=subs(Dt=rhs(eqDt),eqa,eqDx,dat,rhs(eqFo)),dat]));

10
1000
Dt = 625.79166666666666667, Dx = 0.50000000000000000000e-1, Fo = .25000000000000000000, L = .5, T0 = 300, E = 500, h = 10, T[f] = 0., rho = 2300., c = 653., k = 1.5, alpha = .6, epsilon_ = .8, c[o] = ...
Dt = 625.79166666666666667, Dx = 0.50000000000000000000e-1, Fo = .25000000000000000000, L = .5, T0 = 300, E = 500, h = 10, T[f] = 0., rho = 2300., c = 653., k = 1.5, alpha = .6, epsilon_ = .8, c[o] = ...
Dt = 625.79166666666666667, Dx = 0.50000000000000000000e-1, Fo = .25000000000000000000, L = .5, T0 = 300, E = 500, h = 10, T[f] = 0., rho = 2300., c = 653., k = 1.5, alpha = .6, epsilon_ = .8, c[o] = ...
Dt = 625.79166666666666667, Dx = 0.50000000000000000000e-1, Fo = .25000000000000000000, L = .5, T0 = 300, E = 500, h = 10, T[f] = 0., rho = 2300., c = 653., k = 1.5, alpha = .6, epsilon_ = .8, c[o] = ...
Dt = 625.79166666666666667, Dx = 0.50000000000000000000e-1, Fo = .25000000000000000000, L = .5, T0 = 300, E = 500, h = 10, T[f] = 0., rho = 2300., c = 653., k = 1.5, alpha = .6, epsilon_ = .8, c[o] = ...

Inicialización: T[i , j=0]=T0.

> j:=0:for i from 0 to N do T[i,j]:=subs(datnum0,T0):od:

Iteraciones.

> for j from 0 to M-1 do eq_0_:=subs(datnum0,eq_0);assign(%); for i from 1 to N-1 do eq_i_:=subs(datnum0,eq_i);assign(%);od; eq_N_:=subs(TN=T0,datnum0,eq_N);assign(%);od:

Resultado.

> Mdis:=99:pl:=subs(datnum0,[seq([seq([L*(i/N),T[i,jdis*trunc(M/Mdis)]],i=0..N)],jdis=0..Mdis)]):plot(pl,color=black);

Plot_2d

> plot(subs(datnum0,[seq([seq([subs(datnum0,M*Dt)*(j/M),T[i,j*trunc(M/Mdis)]],j=1..Mdis)],i=1..N)]),color=black);

Plot_2d

Visualización interactiva:

> with(plots):
Ts :=subs(datnum0,[seq([seq([L*(i/N),j*Dt,T[i,j]],i=0..N)],j=1..M)]): surfdata( Ts, axes=boxed, style=patch,labels=[`L`,`time`,`Temp`] );

Plot_2d

>