> restart:#"m11_p36"

Considérese la difusión plana de energía térmica y de una especie química a través de un medio de espesor L, manteniendo en régimen estacionario los extremos en valores constantes. Se pide:

a) Ecuación y condiciones de contorno en cada caso.

b) Solución analítica.

c) Aplicación numérica para L=1 cm, siendo el medio aire y la especie vapor de agua evaporándose desde una superficie líquida.

Datos:

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

> su1:="Aire":su2:="CH3OH":su2:="H2O":dat:=[L=0.01*m_];

[L = `+`(`*`(0.1e-1, `*`(m_)))]

Eqs. const.:

> Adat:=get_gas_data(su1):Adat:=op(subs(M=M[a],[Adat]));Ygdat:=get_gas_data(su2):Yldat:=get_liq_data(su2):Ydat:=op(subs(M=M[v],[Ygdat,Yldat])):get_pv_data(su2):dat:=op(dat),Adat,Const,SI1,SI2:

M[a] = `+`(`/`(`*`(0.29e-1, `*`(kg_)), `*`(mol_))), T[b] = `+`(`*`(82., `*`(K_))), T[cr] = `+`(`*`(132., `*`(K_))), p[cr] = `+`(`*`(0.3750e7, `*`(Pa_))), c[p] = `+`(`/`(`*`(1004., `*`(J_)), `*`(kg_, `...
M[a] = `+`(`/`(`*`(0.29e-1, `*`(kg_)), `*`(mol_))), T[b] = `+`(`*`(82., `*`(K_))), T[cr] = `+`(`*`(132., `*`(K_))), p[cr] = `+`(`*`(0.3750e7, `*`(Pa_))), c[p] = `+`(`/`(`*`(1004., `*`(J_)), `*`(kg_, `...

a) Ecuación y condiciones de contorno en cada caso.

Eqs. balance:

> eqBM:=diff(rho(x,t),t)=-'diff(rho(x,t)*v(x,t),x)';eqBM_:=v(x,t)=V;eqBMi:='diff(rhoi(x,t),t)=Di*diff(rhoi(x,t),x,x)-diff(rhoi(x,t)*v(x,t),x)';eqBMi_:=0=Di*diff(yi(x),x,x)-diff(yi(x)*V,x);eqBE:=subs(phi=0,eq11_6_11);eqBE_:=0=a*diff(T(x),x,x);

diff(rho(x, t), t) = `+`(`-`(diff(`*`(rho(x, t), `*`(v(x, t))), x)))
v(x, t) = V
diff(rhoi(x, t), t) = `+`(`*`(Di, `*`(diff(rhoi(x, t), x, x))), `-`(diff(`*`(rhoi(x, t), `*`(v(x, t))), x)))
0 = `+`(`*`(Di, `*`(diff(diff(yi(x), x), x))), `-`(`*`(diff(yi(x), x), `*`(V))))
diff(T(x, t), t) = `*`(a, `*`(diff(diff(T(x, t), x), x)))
0 = `*`(a, `*`(diff(diff(T(x), x), x)))

b) Solución analítica.

Nótese que el balance energético da una recta pero el másico sólo la da si V<<Di/L.

> eqBC1:=T(0)=T0;eqBC2:=T(L)=T1;dsolT:=dsolve({eqBE_,eqBC1,eqBC2},T(x));eqBC1:=yi(0)=yi0;eqBC2:=yi(L)=yi1;dsolyi:=dsolve({eqBMi_,eqBC1,eqBC2},yi(x));dsolyilin:=yi(x)=collect(simplify(convert(series(rhs(dsolyi),V=0,3),polynom)),x);plot([seq(subs(yi1=0,yi0=1,L=1,V=1,Di=10^i,rhs(dsolyi)),i=-1..1)],x=0..1,color=black);

T(0) = T0
T(L) = T1
T(x) = `+`(`-`(`/`(`*`(`+`(`-`(T1), T0), `*`(x)), `*`(L))), T0)
yi(0) = yi0
yi(L) = yi1
yi(x) = `+`(`/`(`*`(`+`(`*`(yi0, `*`(exp(`/`(`*`(V, `*`(L)), `*`(Di))))), `-`(yi1))), `*`(`+`(`-`(1), exp(`/`(`*`(V, `*`(L)), `*`(Di)))))), `-`(`/`(`*`(`+`(`-`(yi1), yi0), `*`(exp(`/`(`*`(V, `*`(x)), ...
yi(x) = `+`(`/`(`*`(`+`(yi1, `-`(yi0)), `*`(x)), `*`(L)), yi0)
Plot_2d

Pero la V se obtiene de la ecuación de Fick:

> eqFick:=rho*V=eval(-rho*Di*subs(x=0,diff(subs(dsolyi,yi(x)),x)));eqFick:=rho*V=convert(series(rhs(%),V=0,2),polynom);

`*`(rho, `*`(V)) = `/`(`*`(rho, `*`(`+`(`-`(yi1), yi0), `*`(V))), `*`(`+`(`-`(1), exp(`/`(`*`(V, `*`(L)), `*`(Di))))))
`*`(rho, `*`(V)) = `/`(`*`(rho, `*`(`+`(`-`(yi1), yi0), `*`(Di))), `*`(L))

y V<<Di/L sólo si yi<<1.

c) Aplicación numérica para L=1 cm, siendo el medio aire y la especie vapor de agua evaporándose desde una superficie líquida.

> eqBC0:=yi0=x[v]*M[v]/(x[v]*M[v]+x[a]*M[a]);eqBC0:=yi0=x[v]*M[v]/M[a];eqBC0:=yi0='(M[v]/M[a])*pv(T)/p0';eqBC0_:=evalf(subs(T=T0,dat,Ydat,dat,eqBC0));

yi0 = `/`(`*`(x[v], `*`(M[v])), `*`(`+`(`*`(x[v], `*`(M[v])), `*`(x[a], `*`(M[a])))))
yi0 = `/`(`*`(x[v], `*`(M[v])), `*`(M[a]))
yi0 = `/`(`*`(M[v], `*`(pv(T))), `*`(M[a], `*`(p0)))
yi0 = 0.1065e-1

Suponiendo una humedad ambiente del 50%:

> eqBC1:=yi1=yi0*phi;eqBC1_:=subs(phi=0.5,eqBC0_,yi1=yi0*phi);

yi1 = `*`(yi0, `*`(phi))
yi1 = 0.5325e-2

Para el aire, tomando Di=a=k/(rho*cp):

> eqD:=a=k/(rho*c[p]);eqD:=a=k*R*T0/(p0*c[p]);eqD_:=subs(dat,eqD);eqD__:=subs(dat,nu=mu*R*T0/p0);eqD:=Di=a;V_:=subs(eqD,eqD_,eqBC0_,eqBC1_,dat,solve(eqFick,V));

a = `/`(`*`(k), `*`(rho, `*`(c[p])))
a = `/`(`*`(k, `*`(R, `*`(T0))), `*`(p0, `*`(c[p])))
a = `+`(`/`(`*`(0.1973e-4, `*`(`^`(m_, 2))), `*`(s_)))
nu = `+`(`/`(`*`(0.1486e-4, `*`(`^`(m_, 2))), `*`(s_)))
Di = a
`+`(`/`(`*`(0.1051e-4, `*`(m_)), `*`(s_)))

Del balance másico en la interfase líquida:

> eqBM0:=-(p0/(R*T0))*Di*(yi1-yi0)/L=-rholiq*dx0/dt;eqBM0_:=dx0_dt_=solve(subs(eqD,eqD_,rholiq=rho,Ydat,eqBC0_,eqBC1_,dat,eqBM0),dx0)/dt;'dx0_dt_'=rhs(%)*(86400*s_/dia_)*(1000*mm_/m_);

`+`(`-`(`/`(`*`(p0, `*`(Di, `*`(`+`(yi1, `-`(yi0))))), `*`(R, `*`(T0, `*`(L)))))) = `+`(`-`(`/`(`*`(rholiq, `*`(dx0)), `*`(dt))))
dx0_dt_ = `+`(`-`(`/`(`*`(0.7915e-8, `*`(m_)), `*`(s_))))
dx0_dt_ = `+`(`-`(`/`(`*`(.6839, `*`(mm_)), `*`(dia_))))

es decir, se evapora menos de 1 mm/dia (sólo válido sin convección, p.e. en un tubo de ensayo).

Para la difusión térmica pura, sin difusión de masa, sería:

> eqFourier:=-k*(T1-T0)/L;

`+`(`-`(`/`(`*`(k, `*`(`+`(T1, `-`(T0)))), `*`(L))))

que para 1 cm de aire, por unidad de área y de salto térmico sería:

> eqFourier_:=subs(T1=1*K_+T0,dat,SI1,eqFourier)/K_;

`+`(`-`(`/`(`*`(2.4, `*`(W_)), `*`(`^`(m_, 2), `*`(K_)))))

Pero la evaporación másica lleva asociado un flujo energético, que puede provenir axialmente:

> eqBE0:=k*(T1-T0)/L=-rholiq*(dx0/dt)*h[lv]+(p0/(R*T0))*V*c[p]*(T1-T0);eqBE0:=k*(T1-T0)/L=-rholiq*(dx0/dt)*h[lv0];DT01_:=simplify(subs(dx0=dx0_dt_*dt,eqBM0_,dat,rholiq=rho,Ydat,dat,solve(eqBE0,T1)-T0)):DT01=evalf(%,2);

`/`(`*`(k, `*`(`+`(T1, `-`(T0)))), `*`(L)) = `+`(`-`(`/`(`*`(rholiq, `*`(dx0, `*`(h[lv]))), `*`(dt))), `/`(`*`(p0, `*`(V, `*`(c[p], `*`(`+`(T1, `-`(T0)))))), `*`(R, `*`(T0))))
`/`(`*`(k, `*`(`+`(T1, `-`(T0)))), `*`(L)) = `+`(`-`(`/`(`*`(rholiq, `*`(dx0, `*`(h[lv0]))), `*`(dt))))
DT01 = `+`(`*`(7.4, `*`(K_)))

Si el ambiente está a T0=15 ºC, el agua del tubo se enfriaría hasta T0-DT=7,5 ºC si las paredes del tubo fuesen aislantes; si no, el balance energético de agua sería:

> eqBE0:=U*A[lateral]*DT+A[cross]*k*DT/L=-A[cross]*rholiq*(dx0/dt)*h[lv0];DT01__:=simplify(subs(U=5*W_/(m_^2*K_),A[lateral]=5*A[cross],dx0=dx0_dt_*dt,eqBM0_,dat,rholiq=rho,Ydat,dat,solve(eqBE0,DT))):DT01=evalf(%,2);

`+`(`*`(U, `*`(A[lateral], `*`(DT))), `/`(`*`(A[cross], `*`(k, `*`(DT))), `*`(L))) = `+`(`-`(`/`(`*`(A[cross], `*`(rholiq, `*`(dx0, `*`(h[lv0])))), `*`(dt))))
DT01 = `+`(`*`(.65, `*`(K_)))

p.e. para U=5 W/(m2.K) y Alat/Across=5, lo que indica que normalmente la transmitancia térmica lateral es la dominante.

>