Problema N° 1

Referencia: Problema 19.1 McCabe (2002)

Un método de separar acetona de acetato de celulosa consiste en pasar una corriente de aire sobre las fibras de acetato de celulosa. Para conocer las propiedades de las mezclas aire-acetona, el departamento de control de proceso necesita una gráfica de humedad para el sistema aire-acetona. Después de una investigación, se encontró que sería satisfactorio un intervalo de humedad de 0 a 6.0 kg/kg y de temperatura de 5 °C a 55 °C. Trace las siguientes regiones de una gráfica de humedad para el sistema aire-acetona a la presión total de 101,325 kPa.

[2 pto] Líneas de humedad porcentual para 50 % y 100 %
[2 pto] Volumen saturado contra temperatura
[1 pto] Calor latente de la acetona contra temperatura
[1 pto] Calor húmedo contra humedad
[3 pto] Líneas de enfriamiento adiabático para temperaturas de saturación adiabática de 20 °C y 40 °C
[3 pto] Líneas de temperatura de bulbo húmedo de 20 °C y 40 °C.

Los datos necesarios son:

Para el vapor de acetona, Cp = 1.47 kJ/(kg K) y h/(MBky) = 1.7 kJ/(kg K)

html('<h2>Solución</h2>')

########## Datos de equilibrio
datos_t1     =  [ 10 ,   20  ,  30  ,  40  ,  50  ,  56.1,  60  ,   70  ,   80  ]   # °C
datos_pvap   =  [115.6, 179.6, 281.0, 420.1, 620.9, 760.0, 860.5, 1189.4, 1611.0]   # mmHg
datos_t2     =  [  0,  20,  40, 56.1,  60,  80]        # °C
datos_lambda =  [564, 552, 536, 521 , 517, 495]        # kJ/kg

########## Datos del enunciado
Yp_min       =    0         # kg/kg
Yp_max       =    6         # kg/kg
t_min        =    5         # °C
t_max        =   55         # °C
p_t          =  101.325     # kPa

C_pAG        =    1.47      # kJ/(kg K)
h_div_MBky   =    1.7       # kJ/(kg K)

########## Datos adicionales
M_A          =   58.08      # acetona
M_B          =   28.97      # aire
R            =    8.314     # [kJ/(kmol K)]
t_0          =    0.01      # [°C]

C_pBG(t)      =  R/M_B * (3.355 + 0.575*10^-3*(t+273.15) - 0.016*10^ 5/(t+273.15)^2) # [kJ/(kmol K)]


########## Ajuste de datos
# Conversión de unidades
datos_pvap   =  [p_mmHg * 101.325/760 for p_mmHg in datos_pvap]   # kPa

# Unir los valores en pares
datos_tpv    =  zip(datos_t1, datos_pvap)
datos_tlam   =  zip(datos_t2, datos_lambda)

# Ajuste a un modelo
%var  a0, a1, a2, a3, a4
modelo(t)    =  a0 + a1*t + a2*t^a3
ajuste       =  find_fit (datos_tpv, modelo, solution_dict = True)
p_A_sat(t)   =  modelo(a0 = ajuste[a0], a1 = ajuste[a1], a2 = ajuste[a2], a3 = ajuste[a3])

# Los spline son útiles para interpolar pero no permiten extrapolar
p_vap        =  spline(datos_tpv)
lambda_      =  spline(datos_tlam)
lambda_0     =  lambda_(t_0)

########## Funciones
%var   t, Yp
Y_p  (t, pt)     =  p_A_sat(t)/(p_t-p_A_sat(t))*M_A/M_B                # [kg/kg]
Y_s  (t, pt)     =  M_A*p_A_sat(t)/(M_B*(pt-p_A_sat(t)))               # [kg/kg]
C_s  (t, Yp)     =  C_pBG(t) + C_pAG*Yp                                # [kJ/(kg K)]
H    (t, Yp)     =  C_s(t, Yp) * (t-t_0) + lambda_0 * Yp               # [kJ/kg]
H_s  (t, pt)     =  C_s(t, Y_s(t,pt)) * (t-t_0) + lambda_0 * Y_s(t,pt) # [kJ/kg]
v_H  (t, Yp, pt) =  R*(t+273.15)/pt * (1/M_B+Yp/M_A)                   # [m³/kg]


# ************************** SOLUCIÓN **************************
COLOR=sorted(colors)

########## a) Líneas de humedad porcentual para 50 % y 100 %
Y_P1      = 0.50
Y_P2      = 1.00
gr_YP50   =  plot(Y_P1 * Y_s(t, p_t), t_min, t_max, legend_label="50 %" , color='blue')
gr_YP100  =  plot(Y_P2 * Y_s(t, p_t), t_min, t_max, legend_label="100 %", color='green')
gr_YP     =  gr_YP50 + gr_YP100
gr_YP.show(title='Humedad porcentual', ymin=Yp_min, ymax=Yp_max, axes_labels=['$t$/($\degree$C)', '$Y\'$/(kg/kg)'], figsize=(6,4))

########## b) Volumen saturado
v_sat(t, pt)  =  v_H (t, Y_s(t, pt), pt)
gr_vsat   =  plot(v_sat(t, p_t), xmin=t_min, xmax=t_max, axes_labels = ['$t$/($\degree$C)', '$v_s$/(m^3/kg)'], figsize=(6,4))

gr_vsat.show(title = 'Volumen saturado')

########## c) Calor latende de la acetona
gr_lambda =  plot(lambda_, t_min, t_max, axes_labels = ['$t$/($\degree$C)', '$\lambda$/(kJ/kg)'])
gr_lambda.show(title = 'Calor latende de la acetona', figsize=(6,4))

########## d) Calor húmedo
gr_calhum = plot(C_pBG(t_0) + C_pAG*Yp, xmin=Yp_min, xmax=Yp_max, axes_labels = ['$Y\'$/(kg/kg)', '$C_s$/(kJ/kg K)'], figsize=(6,4))
gr_calhum.show(title = 'Calor h\'umedo')

########## e) Líneas de enfriamiento adiabático de 20 °C y 40 °C
#N         =  11       # Número de puntos de cálculo
#lista_tsa =  [t_min + (t_max-t_min)/(N-1) * i for i in range(N)]
#lista_Ysa =  [Y_s(tsa, p_t) for tsa in lista_tsa]
t_sa1     =  20    # °C
t_sa2     =  40    # °C

Yp_1(t)    =  (lambda_(t_sa1)*Y_s(t_sa1,p_t) - C_pBG*(t-t_sa1))/(C_pAG*(t-t_sa1) + lambda_(t_sa1))
Yp_2(t)    =  (lambda_(t_sa2)*Y_s(t_sa2,p_t) - C_pBG*(t-t_sa2))/(C_pAG*(t-t_sa2) + lambda_(t_sa2))

gr_tsa20   =  plot(Yp_1, t_sa1, t_max, legend_label="$t_{sa} = $"+str(t_sa1)+" $\degree$C" , color='blue')
gr_tsa40   =  plot(Yp_2, t_sa2, t_max, legend_label="$t_{sa} = $"+str(t_sa2)+" $\degree$C" , color='red')
gr_tsa     =  gr_YP100 + gr_tsa20 + gr_tsa40
gr_tsa.show(title = 'L\'ineas de enfriamiento adiab\'atico', xmin= t_min, ymin=Yp_min, ymax=3, axes_labels=['$t$/($\degree$C)', '$Y\'$/(kg/kg)'], figsize=(6,4))

########## f) Líneas de temperatura de bulbo húmedo de 20 °C y 40 °C
t_w1       =  20    # °C
t_w2       =  40    # °C

Yw_1(t)    =  (Y_s(t_w1,p_t) - h_div_MBky/lambda_(t_w1) *(t-t_w1))
Yw_2(t)    =  (Y_s(t_w2,p_t) - h_div_MBky/lambda_(t_w2) *(t-t_w2))

gr_tw20    =  plot(Yw_1, t_w1, t_max, legend_label="$t_{w} = $"+str(t_w1)+" $\degree$C" , color='purple')
gr_tw40    =  plot(Yw_2, t_w2, t_max, legend_label="$t_{w} = $"+str(t_w2)+" $\degree$C" , color='orange')
gr_tw      =  gr_YP100 + gr_tw20 + gr_tw40
#gr_tw.show(title = 'L\'ineas de temperatura de bulbo h\'umedo', xmin= t_min, ymin=Yp_min, ymax=3, axes_labels=['$t$/($\degree$C)', '$Y\'$/(kg/kg)'], figsize=(6,4))

gr_twsa = gr_tw + gr_tsa
gr_twsa.show(title = 'L\'ineas de temperatura de bulbo h\'umedo y sat adiab', xmin= t_min, ymin=Yp_min, ymax=3, axes_labels=['$t$/($\degree$C)', '$Y\'$/(kg/kg)'], figsize=(6,4))

Solución

Problema N° 2

Referencia: Problema 19.3 McCabe (2002)

Una torre de enfriamiento a contracorriente con relleno celular de 76 cm fue diseñada para una temperatura de aproximación de 6 °C y un intervalo de enfriamiento de 9 °C cuando la temperatura de bulbo húmedo es de 24 °C. La sección transversal de la torre es de 11 m * 11 m, y las velocidades de los flujos de aire y agua son 14 810 m³/min a 32 °C y 22.7 m³/min.

[8 pto] ¿Cuál es la entalpía del aire a la entrada y la variación de entalpía a través de la torre?
[4 pto] Calcule NTU y HTU para estas condiciones.
[8 pto] ¿Cuánto relleno más sería necesario si la temperatura del bulbo húmedo fuera de 26 °C y la temperatura de aproximación fuera de 4 °C?


html('<h2>Solución</h2>')

########## Datos del enunciado
H_c             =   76     # cm
D_t_aprox       =    6     # °C
Delta_tL        =    9     # °C
t_w             =   24     # °C
A_c             =  11*11   # m²
q_G_entra       =  14810   # m³/min
t_G_entra       =   32     # °C
q_L_entra       =   22.7   # m³/min

########## Datos adicionales
p_t             =  101.325     # kPa
R               =    8.314     # kJ/(kmol K)
t_0             =    0.01      # °C
N_Le            =    0.855
M_A             =   18.01   # agua
M_B             =   28.84   # aire

lambda_A(t)     = 2501.6 + 2.2908456564*t-0.0015914753*t^2    # [kJ/kg]

C_pAL(t)        = R/M_A * (8.712 + 1.25 *10^-3*(t+273.15) - 0.18 *10^-6*(t+273.15)^2) # [kJ/(kg K)]
C_pAG(t)        = R/M_A * (3.470 + 1.450*10^-3*(t+273.15) + 0.121*10^ 5/(t+273.15)^2) # [kJ/(kg K)]
C_pBG(t)        = R/M_B * (3.355 + 0.575*10^-3*(t+273.15) - 0.016*10^ 5/(t+273.15)^2) # [kJ/(kg K)]

# Constantes de Antoine para el agua, [0.01 °C a 100 °C], log10 P_sat/kPa = A - B / (t/°C + C)
A               =    7.110969267248
B               = 1681.106698364
C               =  229.2939885035
p_A_sat(t)      =   10^(A-B/(t+C)) # [kPa]

########### Funciones
C_s(t, Y)       =  C_pBG(t) + C_pAG(t)*Y
Y_s(t, p)       =  M_A*p_A_sat(t)/(M_B*(p-p_A_sat(t)))
H(t,Y)          =  C_s(t, Y) * (t-t_0) + lambda_A(t_0) * Y
H_s(t,p)        =  C_s(t, Y_s(t,p)) * (t-t_0) + lambda_A(t_0) * Y_s(t, p)
v_H(t,Y,p)      =  R*(t+273.15)/p * (1/M_B+Y/M_A)

########## Condiciones del líquido
rho_L           =  997    # kg/m³
t_L_sale        =  t_w + D_t_aprox       # °C
t_L_entra       =  t_L_sale + Delta_tL   # °C
L_entra         =  q_L_entra*rho_L/A_c         # kg/(m² min)

########## Condiciones del aire
Y_w             =  Y_s(t_w, p_t)
Y_entra         =  Y_w - (t_G_entra - t_w)*N_Le^(2/3)*C_pBG(t_G_entra)/lambda_A(t_w)
C_s_entra       =  C_s(t_G_entra, Y_entra)
H_entra         =  H(t_G_entra,Y_entra)
v_H_entra       =  v_H(t_G_entra, Y_entra, p_t)
G_s             =  q_G_entra/(v_H_entra * A_c)

# Línea de operación
H_G(tL)         =  L_entra*C_pAL(tL)/G_s * (tL-t_L_sale) + H_entra
H_sale          =  H_G(t_L_entra)


print "***** Parte a *****"
print "tL_sale  = %8.3f" %t_L_sale, "°C"
print "tL_entra = %8.3f" %t_L_entra, "°C"
print "L'_entra = %8.3f" %L_entra, "kg/(m² min)"
print "Y'_w     = %8.3f" %Y_w, "kg/kg"
print "Y'_entra = %8.3f" %Y_entra, "kg/kg"
print "Cs_entra = %8.3f" %C_s_entra, "kJ/(kg K)"
print "H'_entra = %8.3f" %H_entra, "kJ/kg"
print "vH_entra = %8.3f" %v_H_entra, "m³/kg"
print "G's      = %8.3f" %G_s, "kg/(m² min)"
print "H'_sale  = %8.3f" %H_sale, "kJ/kg"
print "Delta H  = %8.3f" %(H_sale-H_entra), "kJ/kg"

print ""
print "***** Parte b *****"
########## NTU y HTU
%var t

plot(H_s(t,p_t), xmin = min(t_L_sale,t_G_entra)-5 , xmax = t_L_entra+2, legend_label = "$H\'_s(t)$", axes_labels=['$t$/($\degree$C)','$H\'_G$/(kJ/kg)']) \
    + line([(t_L_sale, H_entra),(t_L_entra, H_G(t_L_entra,L_entra, G_s))], color = 'purple', legend_label = "$H\'_G (t_L)$", figsize=(6,4))

# Paso f(H) --> g(t)
g(t)            =  derivative(H_G(t,L_entra, G_s)) / (H_s(t,p_t)-H_G(t,L_entra, G_s))
integral_g      =  numerical_integral(g, t_L_sale, t_L_entra)

# *********** Resultados, parte b) ********************************************************
plot(g, xmin = t_L_sale, xmax = t_L_entra, axes_labels = ['$t$/($\degree$C)', r'$\frac{d H\'(t)}{H^{*}-H}$'], fill = 'axis', color = 'red', figsize=(6,4))
print "Se hace cambio de variable para integrar más fácilmente en función de t en vez de H\'"

N_tOG           =  integral_g[0]
H_tOG           =  H_c/N_tOG           # cm

print ""
print "NTU = N_tOG = %8.3f" %N_tOG
print "HTU = H_tOG = %8.3f cm = %5.2f m" %(H_tOG, H_tOG/100)

print ""
print "***** Parte c *****"

t_w             =   26     # °C
D_t_aprox       =    4     # °C

########## Nuevas condiciones del líquido
t_L_sale        =  t_w + D_t_aprox       # °C
t_L_entra       =  t_L_sale + Delta_tL   # °C
L_entra         =  q_L_entra*rho_L/A_c         # kg/(m² min)

########## Nuevas condiciones del aire
Y_w             =  Y_s(t_w, p_t)
Y_entra         =  Y_w - (t_G_entra - t_w)*N_Le^(2/3)*C_pBG(t_G_entra)/lambda_A(t_w)
C_s_entra       =  C_s(t_G_entra, Y_entra)
H_entra         =  H(t_G_entra,Y_entra)
v_H_entra       =  v_H(t_G_entra, Y_entra, p_t)
G_s             =  q_G_entra/(v_H_entra * A_c)

########## Nueva línea de operación
H_G(tL)         =  L_entra*C_pAL(tL)/G_s * (tL-t_L_sale) + H_entra
H_sale          =  H_G(t_L_entra)

print "tL_sale  = %8.3f °C" %t_L_sale
print "tL_entra = %8.3f °C" %t_L_entra
print "Y'_w     = %8.3f kg/kg" %Y_w
print "Y'_entra = %8.3f kg/kg" %Y_entra
print "Cs_entra = %8.3f kJ/(kg K)" %C_s_entra
print "H'_entra = %8.3f kJ/kg" %H_entra
print "vH_entra = %8.3f m³/kg" %v_H_entra
print "G's      = %8.3f kg/(m² min)" %G_s

########## Nuevo NTU
plot(H_s(t,p_t), xmin = min(t_L_sale,t_G_entra)-5 , xmax = t_L_entra+2, legend_label = "$H\'_s(t)$", axes_labels=['$t$/($\degree$C)','$H\'_G$/(kJ/kg)']) \
    + line([(t_L_sale, H_entra),(t_L_entra, H_G(t_L_entra,L_entra, G_s))], color = 'purple', legend_label = "$H\'_G (t_L)$", figsize=(6,4))

# Paso f(H) --> g(t)
g(t)            =  derivative(H_G(t,L_entra, G_s)) / (H_s(t,p_t)-H_G(t,L_entra, G_s))
integral_g      =  numerical_integral(g, t_L_sale, t_L_entra)

plot(g, xmin = t_L_sale, xmax = t_L_entra, axes_labels = ['$t$/($\degree$C)', r'$\frac{d H\'(t)}{H^{*}-H}$'], fill = 'axis', color = 'red', figsize=(6,4))

N_tOG           =  integral_g[0]
H_c_nueva       =  N_tOG * H_tOG
Delta_Z         =  H_c_nueva - H_c

print ""
print "NTU = N_tOG = %7.4g" %N_tOG
print "Nueva H_c   = %7.4g cm" %H_c_nueva
print "Diferencia  = %7.4g cm" %Delta_Z

Solución

***** Parte a *****
tL_sale  =   30.000 °C
tL_entra =   39.000 °C
L'_entra =  187.040 kg/(m² min)
Y'_w     =    0.019 kg/kg
Y'_entra =    0.016 kg/kg
Cs_entra =    1.043 kJ/(kg K)
H'_entra =   73.526 kJ/kg
vH_entra =    0.891 m³/kg
G's      =  137.446 kg/(m² min)
H'_sale  =  124.889 kJ/kg
Delta H  =   51.363 kJ/kg

***** Parte b *****

Se hace cambio de variable para integrar más fácilmente en función de t en vez de H'

NTU = N_tOG =    1.829
HTU = H_tOG =   41.560 cm =  0.42 m

***** Parte c *****
tL_sale  =   30.000 °C
tL_entra =   39.000 °C
Y'_w     =    0.021 kg/kg
Y'_entra =    0.019 kg/kg
Cs_entra =    1.049 kJ/(kg K)
H'_entra =   81.736 kJ/kg
vH_entra =    0.895 m³/kg
G's      =  136.762 kg/(m² min)

NTU = N_tOG =   2.621
Nueva H_c   =   108.9 cm
Diferencia  =   32.91 cm

Problema N° 3

Referencia: Adaptado de Ejemplo 13.9 de Coulson and Richardson, Volumen 1 (1999)

Se va a enfriar agua desde 55 °C hasta 25 °C por contacto a contracorriente con aire que entra a 23 °C con una humedad relativa de 25 %. El aire fluye a 2000 m³/(m² h) y el agua a 950 kg/(m² h). La resistencia a la transferencia de calor y de masa puede suponerse que se encuentra en la fase gaseosa y el coeficiente volumétrico global de transferencia de masa puede obtenerse de $K_{Y'} a \rho_G = 0.2 \;\mathrm{s^{-1}}$ .

Tomando la presión como 97.0 kPa calcule:

[Puntos extra] Condiciones del aire de salida (humedad absoluta, humedad relativa, temperatura y entalpía)
Altura de relleno que se requiere
Relación entre el flujo de gas empleado con respecto al flujo mínimo requerido


html('<h2>Solución</h2>')

# ========== Datos del enunciado
t_L_entra       =   55        # °C
t_L_sale        =   25        # °C
t_G_entra       =   23        # °C
Y_R_entra       =    0.2500   # Pa/Pa
q_G_entra       = 2000        # m³/(m² h)
L_entra         =  950        # kg/(m² h)
K_Y_a_div_rho_G =    0.2*3600 # h^-1 <-- s^-1
p_t             =   97.0      # kPa

t, L            = var ('t, L')

# ========== Constantes y parámetros generales
M_A             =   18.01    # agua
M_B             =   28.84    # aire
t_0             =    0       # [°C]
R               =    8.314   # [kJ/(kmol K)]

# Kloppers, J. C., & Kröger, D. G. (2004) https://doi.org/10.1080/104077890504113
p_A_sat(t)      =  10^(10.79586*(1-273.16/(t+273.15)) + 5.02808*log(273.16/(t+273.15),10) + 1.50474*10^-4 * (1-10^(-8.29692*((t+273.15)/273.16 - 1))) + 4.2873*10^-4 * (10^(4.76955*(1-273.16/(t+273.15))) - 1) + 2.786118312) / 1000   # kPa
C_pBG(t)        =  (1045.356 - 0.3161783*(t+273.15) + 7.083814*10^-4*(t+273.15)^2 - 2.705209*10^-7*(t+273.15)^3) / 1000  # kJ/(kg K)
C_pAG(t)        =  (1360.5 + 2.31334*(t+273.15) - 2.46784*10^-10*(t+273.15)^5 + 5.91332 * 10^-13*(t+273.15)^6) / 1000    # kJ/(kg K)
C_pAL(t)        =  (8155.99 - 28.0627*(t+273.15) + 0.0511283*(t+273.15)^2 - 2.17582*10^-13 * (t+273.15)^6)/1000          # kJ/(kg K)
lambda_A(t)     =  (3483181.4 - 5862.7703*(t+273.15) + 12.139568*(t+273.15)^2 - 0.0140290431*(t+273.15)^3) / 1000        # kJ/kg
lambda_A0       =  lambda_A(t_0)

# ========== Funciones
C_s(t, Y)       =  C_pBG(t) + C_pAG(t)*Y
Y_s(t, p)       =  M_A*p_A_sat(t)/(M_B*(p-p_A_sat(t)))
H(t,Y)          =  C_s(t, Y) * (t-t_0) + lambda_A(t_0) * Y
H_s(t,p)        =  C_s(t, Y_s(t,p)) * (t-t_0) + lambda_A(t_0) * Y_s(t, p)
v_H(t,Y,p)      =  R*(t+273.15)/p * (1/M_B+Y/M_A)
rho_G(T, Y, p)  =  p*(M_B+M_A*Y)/(R*(T+273.15))

# ========== Condiciones de entrada
p_A_entra       =  Y_R_entra*p_A_sat(t_G_entra)
Y_entra         =  M_A/M_B * p_A_entra/(p_t-p_A_entra)
C_s_entra       =  C_s(t_G_entra, Y_entra)
H_entra         =  H(t_G_entra,Y_entra)
v_H_entra       =  v_H(t_G_entra, Y_entra, p_t)
rho_G_entra     =  rho_G(t_G_entra, Y_entra, p_t)

G_s             =  q_G_entra/v_H_entra

print "pA_entra = %8.3g kPa" %p_A_entra
print "Y'entra  = %8.3g kg/kg" %Y_entra
print "Cs_entra = %8.4g kJ/(kg K)" %C_s_entra
print "H'entra  = %8.3g kJ/kg" %H_entra
print "vH_entra = %8.3g m³/kg" %v_H_entra
print "G's      = %8.4g kg/(m² h)" %G_s

# ========== Línea de operación
H_G(t,L,Gs)     =  L*C_pAL(t)/Gs * (t-t_L_sale) + H_entra
H_sale          =  H_G(t_L_entra, L_entra, G_s)


# ========== Flujo limitante de gas
LdivG_lim       =  var('LdivG_lim')
Relac_lim(t)    =  solve(H_G(t,LdivG_lim,1) == H_s(t,p_t), LdivG_lim)[0].rhs()
sol_Rel_lim     =  find_local_minimum(Relac_lim(t), t_L_sale, t_L_entra)    # Devuelve [L(t_eq), t_eq]

G_s_min         =  L_entra/sol_Rel_lim[0]
t_eq            =  sol_Rel_lim[1]

# ========== Condiciones de salida
N               = 30 # Cantidad de etapas de cálculo
H_eq_menos_H(t) = H_s(t,p_t)-H_G(t,L_entra,G_s)

Delta_H         =  (H_sale-H_entra)/N
lista_H_G       =  [H_entra + Delta_H*i for i in range(N+1)]
lista_t_I       =  [find_root(H_G(t,L_entra,G_s) == H_g, t_L_sale-0.1, t_L_entra+0.1) for H_g in lista_H_G]
lista_t_eq      =  [find_root(H_s(t,p_t) == H_g, 0, t_L_entra+0.1) for H_g in lista_H_G]
lista_H_I       =  [H_s(t_I,p_t) for t_I in lista_t_I]
lista_t_G       =  vector(RR, N+1)

lista_t_G[0]    =  t_G_entra
if lista_t_I[0] == lista_t_G[0]:
    DH_DtG      =  10^100
else:
    DH_DtG    =  (lista_H_I[0]-lista_H_G[0])/(lista_t_I[0]-lista_t_G[0])

for i in range(1,N+1):
    lista_t_G[i]     = lista_t_I[i-1] - (lista_H_I[i-1] - lista_H_G[i])/DH_DtG
    if lista_t_G[i]  < lista_t_eq[i]:
        lista_t_G[i] = lista_t_eq[i]
    DH_DtG         = (lista_H_I[i]-lista_H_G[i])/(lista_t_I[i]-lista_t_G[i])

t_G_sale        =  lista_t_G[N]
Y_sale          =  numerical_approx(solve(H(t_G_sale,Y) == lista_H_G[N], Y)[0].rhs())
Yn_sale         =  Y_sale * M_B/M_A
y_sale          =  Yn_sale/(1+Yn_sale)

p_A_sale        =  p_t*y_sale
Y_R             =  p_A_sale/p_A_sat(t_G_sale)


# ========== Altura del relleno
# Paso f(H) --> g(t)
g(t)            =  derivative(H_G(t,L_entra, G_s)) / (H_s(t,p_t)-H_G(t,L_entra, G_s))
integral_g      =  numerical_integral(g, t_L_sale, t_L_entra)
graf_NTU        =  plot(g, xmin = t_L_sale, xmax = t_L_entra, axes_labels = ['$t$/($\degree$C)', r'$\frac{d H\'(t)}{H^{*}-H}$'], fill = 'axis', color = 'red')
N_tOG           =  integral_g[0]
H_tOG_entra     =  G_s*(1+Y_entra)/(K_Y_a_div_rho_G*rho_G_entra)
rho_G_sale      =  rho_G(t_G_sale, Y_sale, p_t)
H_tOG_sale      =  G_s*(1+Y_sale)/(K_Y_a_div_rho_G*rho_G_sale)
H_tOG_prom      =  mean([H_tOG_entra, H_tOG_sale])


graf_Hs     = plot(H_s(t,p_t), xmin = min(t_L_sale,t_G_entra) , xmax = t_L_entra+5, legend_label = "$H\'_s(t)$", axes_labels=['$t$','$H\'_G$'])
graf_Gmin   = line([(t_L_sale, H_entra),(t_L_entra, H_G(t_L_entra,Relac_lim(t_eq),1))], color = 'red', legend_label = "$\\frac{L C_{p,A,L}}{G\'_{s,min}}$")
graf_lop    = line([(t_L_sale, H_entra),(t_L_entra, H_G(t_L_entra,L_entra,G_s))], color = 'purple', legend_label = "$H\'_G (t_L)$") \
graf_HG_tG  = line(zip(lista_t_G,lista_H_G), color = 'green', legend_label = "$H\'_G (t_G)$")

graf_1      = graf_Hs + graf_lop + graf_HG_tG
graf_2      = graf_Hs + graf_lop + graf_Gmin


print "\n","a) Condiciones del aire de salida"
graf_1.show(figsize=(6,4))
graf_1.show(xmin = min(t_L_sale,t_G_entra), xmax = t_G_sale+0, ymax = H_s(t_G_sale+2,p_t), figsize=(6,4))
print "\n"
print "Y'sale   = %8.3g kg/kg" %Y_sale
print "YR       = %8.3g %%" %(Y_R*100)
print "H'sale   = %8.3g kJ/kg" %H_sale
print "tG_sale  = %8.3g °C (con %s deltas, ΔH = %.4g kJ/kg)" %(t_G_sale, N, Delta_H)

print "\n","b) Altura de relleno"
graf_NTU.show(figsize=(6,4))
print "\n","N_tOG = %s" %numerical_approx(N_tOG, digits=4)
print "H_tOG,entra = %s m" %numerical_approx(H_tOG_entra, digits=4)
print "H_c = %s m" %numerical_approx(N_tOG * H_tOG_entra, digits=4), "<--- Tomando solo H_tOG a la entrada"
print "H_tOG,sale = %s m" %numerical_approx(H_tOG_sale, digits=4)
print "H_tOG,prom = %s m" %numerical_approx(H_tOG_prom, digits=4)
print "H_c = %s m" %numerical_approx(N_tOG * H_tOG_prom, digits=4), "<--- Considerando H_tOG promedio"

print "\n","c) Relación Gs/Gsmin"
graf_2.show()
print "\n","G's/G's,mín = %s" %numerical_approx(G_s/G_s_min, digits=4), "<--- Respuesta"
print "G's,mín  = %s kg/(m² h)" %numerical_approx(G_s_min, digits=4)


########## Final del código ##########

Solución

pA_entra =    0.702 kPa
Y'entra  =  0.00455 kg/kg
Cs_entra =    1.015 kJ/(kg K)
H'entra  =     34.7 kJ/kg
vH_entra =    0.887 m³/kg
G's      =     2256 kg/(m² h)

a) Condiciones del aire de salida

Y'sale   =   0.0229 kg/kg
YR       =     86.8 %
H'sale   =     87.6 kJ/kg
tG_sale  =     28.8 °C (con 30 deltas, ΔH = 1.761 kJ/kg)

b) Altura de relleno

N_tOG = 0.5450
H_tOG,entra = 2.762 m
H_c = 1.505 m <--- Tomando solo H_tOG a la entrada
H_tOG,sale = 2.835 m
H_tOG,prom = 2.799 m
H_c = 1.525 m <--- Considerando H_tOG promedio

c) Relación Gs/Gsmin

G's/G's,mín = 5.224 <--- Respuesta
G's,mín  = 431.8 kg/(m² h)

Problema N° 4

Referencia: Adaptación de problema 7.10 de Treybal

Se tiene una corriente de 1.3 m³/s de aire, con 38.0 °C de bulbo seco y 30.0 °C de bulbo húmedo, que se desea deshumidificar llevándolo hasta una temperatura de bulbo húmedo de 15.0 °C por contacto con agua a 10.0 °C fluyendo a contracorriente en una torre rellena. Se utilizará 1.4 veces el flujo mínimo de líquido.

Calcule:

Condiciones del aire de salida (humedad, entalpía, temperaturas)
Temperatura del agua de salida

html('<h2>Solución</h2>')

# ========== Datos del enunciado
q_G_entra       =    1.3      # m³/s
t_G_entra       =   38        # °C
t_w_entra       =   30        # °C
t_w_sale        =   15        # °C
t_L_entra       =   10        # °C
L_div_Lmin      =    1.4

# ========== Otros datos
p_t             =  101.325    # kPa
N_Le            =    0.855

# El resto de funciones y constantes ya definidas en problema 1

# ========== Condiciones de entrada
Y_w_entra       =  Y_s(t_w_entra, p_t)           # kg/kg
Y_entra         =  Y_w_entra - (t_G_entra - t_w_entra)*N_Le^(2/3)*C_pBG(t_G_entra)/lambda_A(t_w_entra)    # kg/kg
C_s_entra       =  C_s(t_G_entra, Y_entra)       # kJ/(kg K)
H_entra         =  H(t_G_entra,Y_entra)          # kJ/kg
v_H_entra       =  v_H(t_G_entra, Y_entra, p_t)  # m³/kg
G_s             =  q_G_entra/v_H_entra           # kg/s

print "Y'w,entra = %8.3g kg/kg" %Y_w_entra
print "Y'entra  = %8.3g kg/kg" %Y_entra
print "Cs_entra = %8.4g kJ/(kg K)" %C_s_entra
print "H'entra  = %8.3g kJ/kg" %H_entra
print "vH_entra = %8.3g m³/kg" %v_H_entra
print "G's      = %8.4g kg/s" %G_s


# ========== Condiciones de salida
t_G_sale        =  t_w_sale                      # °C
Y_sale          =  Y_s(t_w_sale, p_t)            # kg/kg
C_s_sale        =  C_s(t_G_sale, Y_sale)         # kJ/(kg K)
H_sale          =  H(t_G_sale,Y_sale)            # kJ/kg
v_H_entra       =  v_H(t_G_sale, Y_sale, p_t)    # m³/kg

print ""
print "a) Condiciones de salida del aire"
print "Por el proceso de deshumidificación se sabe que el aire sale saturado"
print "t_G,sale = t_w = %.3g °C" %t_G_sale
print "Y'sale   = %8.3g kg/kg" %Y_sale
print "H'sale   = %8.3g kJ/kg" %H_sale


# ========== Línea de operación
H_G(t,L,Gs)     =  L*C_pAL(t)/Gs * (t-t_L_entra) + H_sale
t_L_max         =  find_root(H_s(t, p_t) == H_entra, t_L_entra, t_G_entra)

# ========== Flujo mínimo de líquido
LdivG_lim       =  var('LdivG_lim')
Relac_lim(t)    =  solve(H_G(t,LdivG_lim,1) == H_s(t,p_t), LdivG_lim)[0].rhs()
sol_Rel_lim     =  find_local_maximum(Relac_lim(t), t_L_entra, t_L_max)    # Devuelve [L(t_eq), t_eq]
L_min           =  sol_Rel_lim[0]*G_s
t_eq            =  sol_Rel_lim[1]

L_entra         =  L_div_Lmin * L_min
t_L_sale        =  find_root(H_G(t, L_entra, G_s) == H_entra, t_L_entra, t_L_max)


# ========== Gráfica
graf_Hs     = plot(H_s(t,p_t), [t_L_entra-1, t_L_max+1], legend_label = "$H\'_s(t)$", axes_labels=['$t$','$H\'_G$'])
graf_Lmin   = plot(H_G(t,L_min,G_s), [t_L_entra, t_L_max], color = 'red', legend_label = "$\\frac{L_{min} C_{p,A,L}}{G\'_{s}}$")
graf_lop    = plot(H_G(t,L_entra,G_s), [t_L_entra, t_L_sale], color = 'purple', legend_label = "$H\'_G (t_L)$") \
graf_P2     = graf_Hs + graf_Lmin + graf_lop

print "\n","b) Temperatura del líquido de salida"
graf_P2.show(figsize=(6,4))
print "\n"
print "L'min    = %8.3g kg/s" %L_min
print "L'       = %8.3g kg/s" %L_entra
print "t_L,sale = %8.3g °C" %t_L_sale, "<-----"


########## Final del código ##########

Solución

Y'w,entra =   0.0273 kg/kg
Y'entra  =   0.0243 kg/kg
Cs_entra =    1.054 kJ/(kg K)
H'entra  =      101 kJ/kg
vH_entra =     0.92 m³/kg
G's      =    1.413 kg/s

a) Condiciones de salida del aire
Por el proceso de deshumidificación se sabe que el aire sale saturado
t_G,sale = t_w = 15 °C
Y'sale   =   0.0107 kg/kg
H'sale   =     42.1 kJ/kg

b) Temperatura del líquido de salida

L'min    =    0.985 kg/s
L'       =     1.38 kg/s
t_L,sale =     24.4 °C <-----