SharedTareas / Tarea 1.2 | Cálculo de altura de relleno, unidades de transferencia en destilación.sagewsOpen in CoCalc
Soluciones a Tareas Masa 2

Problema N° 1

Valor: 18 puntos

Referencia: Adaptado de Problema 15.D4 Wankat (2008)

Una columna de destilación separa una alimentación formada por 40 % metanol y 60 % agua (fracciones molares). La alimentación es de dos fases con 60 % líquido. El producto destilado debe ser 92 % metanol y los fondos 4 % metanol. Se emplea un rehervidor PARCIAL y un condensador total. El reflujo es líquido saturado. La operación es a 101.3 kPa. Suponga flujo equimolar y emplee L/DL/D = 1.1. Bajo estas condiciones HtyH_{ty} = 0.396 m y HtxH_{tx} = 0.244 m para ambas secciones de enriquecimiento y agotamiento. Determine la altura requerida de ambas secciones empleando el método de las unidades de transferencia para la fase líquida.

# Datos de EVL para el sistema metanol-agua a 1 atm, fracciones molares (Fuente: Tabla 2-7 Wankat 2a ed.)
datos_x   =  [0,  0.02,  0.04,  0.06,  0.08,  0.10,  0.15,  0.20,  0.30,  0.40,  0.50,  0.60,  0.70,  0.80,  0.90,  0.95,  1]
datos_y   =  [0,  0.134, 0.23,  0.304, 0.365, 0.418, 0.517, 0.579, 0.665, 0.729, 0.779, 0.825, 0.870, 0.915, 0.958, 0.979, 1]

# Datos del enunciado
z_F       =   0.40         # mol/mol
q         =   0.60         # mol/mol
x_D       =   0.92         # mol/mol
x_B       =   0.04         # mol/mol
R_D       =   1.1
H_ty_a    =   0.396        # m
H_tx_a    =   0.244        # m
H_ty_e    =   0.396        # m
H_tx_e    =   0.244        # m

# ***** Número de puntos para resolver las integrales, mínimo 5 *****
n_puntos_agot = 50
n_puntos_enri = 50

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

datos_eq  =  zip(datos_x,datos_y)     # Crea una lista de pares de datos x,y

# Se definen algunas variables por utilizar (cuando se requiere resolver ecuaciones por ejemplo)
%var x, y, B, D, xI, yI, a0, a1, a2, a3, a4

# Ajuste de datos de equilibrio a un modelo
modelo(x) =  a1 * x + a2*x^2 + a3 *x^a4        # Modelo arbitrario que ajusta bien
ajuste    =  find_fit(datos_eq, modelo, solution_dict=True)
y_eq(x)   =  modelo(a1 = ajuste[a1], a2 = ajuste[a2], a3 = ajuste[a3], a4 = ajuste[a4])

# Balances de masa
F         = 1000 # kmol/h, Base de cálculo

BM        =  solve([    F ==     D +     B,  \
                    z_F*F == x_D*D + x_B*B], \
                    B, D, solution_dict = True)

B         =  numerical_approx(BM[0][B], digits=4)
D         =  numerical_approx(BM[0][D], digits=4)
R_B       =  (q*F + R_D*D - B) / B

print "F   =", F, "kmol/h, Base de cálculo"
print "B   =", B, "kmol/h"
print "D   =", D, "kmol/h"
print "R_B =", R_B


# Ecuaciones de las líneas de operación
y_q(x)    =  q/(q-1)*x - z_F/(q-1)
y_e(x)    =  R_D/(R_D+1) * x + x_D/(R_D+1)
y_a(x)    =  (R_B+1)/R_B * x - x_B/R_B

x_ali     =  find_root(y_q == y_e,  0, 1)  # Cruce línea q con enriquecimiento
x_aeq     =  find_root(y_q == y_eq, 0, 1)  # Cruce con equilibrio

# Composiciones de entrada y salida
y_a_entra =  y_eq(x_B)
y_a_sale  =  y_q(x_ali)
x_a_entra =  x_ali
x_a_sale  =  numerical_approx(solve(y_a_entra == y_a, x)[0].rhs(), digits = 4)

y_e_entra =  y_a_sale
y_e_sale  =  y_e(x_D)
x_e_entra =  x_D
x_e_sale  =  x_ali

print "\nEn zona de agotamiento"
print "----------------------"
print "x_entra = %.3g" %x_a_entra
print "x_sale  = %.3g" %x_a_sale
print "y_entra = %.3g" %y_a_entra
print "y_sale  = %.3g" %y_a_sale

print "\nEn zona de enriquecimiento"
print "--------------------------"
print "x_entra = %.3g" %x_e_entra
print "x_sale  = %.3g" %x_e_sale
print "y_entra = %.3g" %y_e_entra
print "y_sale  = %.3g" %y_e_sale


# Gráficas de puntos de datos y las diferentes líneas de operación
graf_data =  scatter_plot(datos_eq, xmin=0, xmax = 1, ymin = 0, ymax = 1, axes_labels = ["$x_{metanol}$", "$y_{metanol}$"])
graf_li45 =  plot(x   ,  x,  0    ,  1    , color ="black")
graf_equi =  plot(y_eq,  x,  0    ,  1    , color ="orange"      , legend_label = "Equilibrio")
graf_alim =  plot(y_q ,  x,  x_aeq,  z_F  , color = "red"        , legend_label = "Alim.")
graf_enri =  plot(y_e ,  x,  x_ali,  x_D  , color = "greenyellow", legend_label = "Enriq.")
graf_agot =  plot(y_a ,  x,  x_B  ,  x_ali, color = "green"      , legend_label = "Agot." )

# Gráfica de las composiciones alrededor del rehervidor
graf_rehe =  line([(x_B, x_B), (x_B, y_a_entra), (x_a_sale, y_a_entra)])

graf_prob =  graf_data + graf_li45 + graf_equi + graf_alim + graf_enri + graf_agot + graf_rehe

# Cálculo de condiciones de interfase
kx_div_ky_enri = numerical_approx(H_ty_e/H_tx_e * R_D/(R_D+1), digits = 4)
kx_div_ky_agot = numerical_approx(H_ty_a/H_tx_a * (R_B+1)/R_B, digits = 4)


show(r"$\left( \frac{k_x}{k_y} \right)_{enri} = \frac{H_{ty_e}}{H_{tx_e}} \frac{R_D}{R_D+1} =$",kx_div_ky_enri)
show(r"$\left( \frac{k_x}{k_y} \right)_{agot} = \frac{H_{ty_a}}{H_{tx_a}} \frac{R_B+1}{R_B} =$",kx_div_ky_agot)

# ***** Genera valores de datos a lo largo de las composiciones en entrada y salida de la zona de agotamiento
Delta_xa   = (x_a_entra - x_a_sale) / (n_puntos_agot - 1)
lista_xa   = [x_a_sale + Delta_xa * i for i in range(n_puntos_agot)]
lista_ya   = [y_a(xop) for xop in lista_xa]

lista_xIa = []    # Lista vacía
lista_yIa = []

# Para cada par (x, y) calcula (xI, yI)
for i in range(n_puntos_agot):
    xop = lista_xa[i]    # Lee un valor de x de la línea de operación
    yop = lista_ya[i]    # Lee el correspondiente valor de y

    # Encuentra el valor de xI que iguala la ecuación de la curva de equilibrio con la recta de pendiente -kx/ky
    f_resolver = y_eq(xI) == -kx_div_ky_agot * xI + kx_div_ky_agot * xop + yop
    raiz_xI    = find_root(f_resolver, 0, 1)

    # Calcula el valor de yI = y*(xI) y guarda los datos en las listas
    raiz_yI  = y_eq(raiz_xI)
    lista_xIa.append(raiz_xI)
    lista_yIa.append(raiz_yI)

# Genera la gráfica de líneas con pendiente -kx/ky
graf_lin = line([(lista_xa[0], lista_ya[0]), (lista_xIa[0], lista_yIa[0])], color = 'gray')
for i in range(n_puntos_agot-1):
    graf_lin = graf_lin + line([(lista_xa[i+1], lista_ya[i+1]), (lista_xIa[i+1], lista_yIa[i+1])], color = 'gray')

# ========== Cálculo del número de unidades de transferencia fase líquida (como pide el enunciado) ==========
# Genera lista con valores 1/(x-xI) y luego forma pares con los valores de x
f          = [1/(lista_xa[i]-lista_xIa[i]) for i in range(n_puntos_agot)]
datos_f    = zip(lista_xa, f)

# Ajuste de los datos a un modelo para integrar fácilmente
datos_x_xI = zip(lista_xa, lista_xIa)
modelo(x)  = a0 + a1*x + a2*x^2 + a3*x^a4
ajuste     = find_fit(datos_x_xI, modelo, solution_dict=True)
xI_f(x)    = modelo(a0=ajuste[a0], a1=ajuste[a1], a2=ajuste[a2], a3=ajuste[a3], a4=ajuste[a4])

# Calcula NTU
f_int(x)   = 1/(x-xI_f(x))
integral_f = numerical_integral(f_int, x_a_sale, x_a_entra)
N_tx_a     = integral_f[0]        # La función numerical_integral devuelve dos valores: la integral y el error estimado

# Cálculo de la altura del relleno de la sección de agotamiento
H_c_a      = numerical_approx(H_tx_a*N_tx_a, digits = 3)

# Gráfica
graf_NTU   = scatter_plot(datos_f, axes_labels = ['$x$', r'$\frac{1}{x-x_I}$'], facecolor = 'none') \
    + plot(f_int, xmin = x_a_sale, xmax = x_a_entra, color = 'green', fill = 'axis', fillcolor='green') \
    + text('$ABC = %s$' %n(N_tx_a, digits=3), ((x_a_entra+x_a_sale)/2, max(f)/2), color = 'black')


# ***** Igual que todo lo anterior pero para la zona de enriquecimiento
Delta_xe   = (x_e_entra - x_e_sale) / (n_puntos_enri - 1)
lista_xe   = [x_e_sale + Delta_xe * i for i in range(n_puntos_enri)]
lista_ye   = [y_e(xop) for xop in lista_xe]

lista_xIe  = []
lista_yIe  = []

for i in range(n_puntos_enri):
    xop    = lista_xe[i]
    yop    = lista_ye[i]

    f_resolver = y_eq(xI) == -kx_div_ky_enri * xI + kx_div_ky_enri * xop + yop
    raiz_xI    = find_root(f_resolver, 0, 1)
    raiz_yI    = y_eq(raiz_xI)
    lista_xIe.append(raiz_xI)
    lista_yIe.append(raiz_yI)

for i in range(n_puntos_enri):
    graf_lin = graf_lin + line([(lista_xe[i], lista_ye[i]), (lista_xIe[i], lista_yIe[i])], color = 'darkgray')

g          = [1/(lista_xe[i]-lista_xIe[i]) for i in range(n_puntos_enri)]
datos_g    = zip(lista_xe, g)

datos_x_xI = zip(lista_xe, lista_xIe)
modelo(x)  = a0 + a1*x + a2*x^2 + a3*x^a4
ajuste     = find_fit(datos_x_xI, modelo, solution_dict=True)
xI_g(x)    = modelo(a0=ajuste[a0], a1=ajuste[a1], a2=ajuste[a2], a3=ajuste[a3], a4=ajuste[a4])

g_int(x)   = 1/(x-xI_g)
integral_g = numerical_integral(g_int, x_e_sale, x_e_entra)
N_tx_e     = integral_g[0]

H_c_e      = numerical_approx(H_tx_e*N_tx_e, digits = 3)

graf_NTU = graf_NTU + scatter_plot(datos_g, axes_labels = ['$x$', r'$\frac{1}{x-x_I}$'], facecolor = 'none') \
    + plot(g_int, xmin = x_e_sale, xmax = x_e_entra, color = 'greenyellow', fill = 'axis', fillcolor='greenyellow') \
    + text('$ABC = %s$' %n(N_tx_e, digits=3), ((x_e_entra+x_e_sale)/2, max(g)/2), color = 'black')

# Muestra los resultados
show(graf_prob + graf_lin)
show(graf_NTU)

show(r"$H_{c,\text{agot}} = (H_{tL} N_{tL})_\text{agot} = %s \,\text{m} * %s = %s \,\text{m}$" %(numerical_approx(H_tx_a, digits = 3), numerical_approx(N_tx_a, digits = 3), H_c_a))
show(r"$H_{c,\text{enri}} = (H_{tL} N_{tL})_\text{enri} = %s \,\text{m} * %s = %s \,\text{m}$" %(numerical_approx(H_tx_e, digits = 3), numerical_approx(N_tx_e, digits = 3), H_c_e))

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


Solución

F = 1000 kmol/h, Base de cálculo B = 590.9 kmol/h D = 409.1 kmol/h R_B = 0.7769 En zona de agotamiento ---------------------- x_entra = 0.278 x_sale = 0.126 y_entra = 0.236 y_sale = 0.584 En zona de enriquecimiento -------------------------- x_entra = 0.92 x_sale = 0.278 y_entra = 0.584 y_sale = 0.92
(kxky)enri=HtyeHtxeRDRD+1=\left( \frac{k_x}{k_y} \right)_{enri} = \frac{H_{ty_e}}{H_{tx_e}} \frac{R_D}{R_D+1} = 0.8501\displaystyle 0.8501
(kxky)agot=HtyaHtxaRB+1RB=\left( \frac{k_x}{k_y} \right)_{agot} = \frac{H_{ty_a}}{H_{tx_a}} \frac{R_B+1}{R_B} = 3.712\displaystyle 3.712
Hc,agot=(HtLNtL)agot=0.244m5.19=1.27mH_{c,\text{agot}} = (H_{tL} N_{tL})_\text{agot} = 0.244 \,\text{m} * 5.19 = 1.27 \,\text{m}
Hc,enri=(HtLNtL)enri=0.244m13.2=3.23mH_{c,\text{enri}} = (H_{tL} N_{tL})_\text{enri} = 0.244 \,\text{m} * 13.2 = 3.23 \,\text{m}
# ========== Cálculo del número de unidades de transferencia para la fase de vapor (SOLO POR COMPARAR) ==========

# Se toman los mismos datos de composiciones en la interfase ya calculados anteriormente
# Genera lista con valores 1/(yI-y) y luego forma pares con los valores de y
f          = [1/(lista_yIa[i]-lista_ya[i]) for i in range(n_puntos_agot)]
datos_f    = zip(lista_ya, f)

# Ajuste de los datos a un modelo para integrar fácilmente
datos_y_yI = zip(lista_ya, lista_yIa)
modelo(y)  = a0 + a1*y + a2*y^2 + a3*y^a4
ajuste     = find_fit(datos_y_yI, modelo, solution_dict=True)
yI_f(y)    = modelo(a0=ajuste[a0], a1=ajuste[a1], a2=ajuste[a2], a3=ajuste[a3], a4=ajuste[a4])

# Calcula NTU
f_int(y)   = 1/(yI_f(y)-y)
integral_f = numerical_integral(f_int, y_a_entra, y_a_sale)
N_ty_a     = integral_f[0]        # La función numerical_integral devuelve dos valores: la integral y el error estimado

# Cálculo de la altura del relleno
H_c_a      = numerical_approx(H_ty_a*N_ty_a, digits = 3)

# Gráfica
graf_NTU   = scatter_plot(datos_f, axes_labels = ['$y$', r'$\frac{1}{y_I-y}$'], facecolor = 'none') \
    + plot(f_int, xmin = y_a_entra, xmax = y_a_sale, color = 'green', fill = 'axis', fillcolor='green') \
    + text('$ABC = %s$' %n(N_ty_a, digits=3), ((y_a_entra+y_a_sale)/2, max(f)/2), color = 'black')

# ***** Igual que todo lo anterior pero para la zona de enriquecimiento
f          = [1/(lista_yIe[i]-lista_ye[i]) for i in range(n_puntos_enri)]
datos_f    = zip(lista_ye, f)

datos_y_yI = zip(lista_ye, lista_yIe)
modelo(y)  = a0 + a1*y + a2*y^2 + a3*y^a4
ajuste     = find_fit(datos_y_yI, modelo, solution_dict=True)
yI_f(y)    = modelo(a0=ajuste[a0], a1=ajuste[a1], a2=ajuste[a2], a3=ajuste[a3], a4=ajuste[a4])

f_int(y)   = 1/(yI_f-y)
integral_f = numerical_integral(f_int, y_e_entra, y_e_sale)
N_ty_e     = integral_f[0]

H_c_e      = numerical_approx(H_ty_e*N_ty_e, digits = 3)

graf_NTU = graf_NTU + scatter_plot(datos_f, axes_labels = ['$y$', r'$\frac{1}{y_I-y}$'], facecolor = 'none') \
    + plot(f_int, xmin = y_e_entra, xmax = y_e_sale, color = 'greenyellow', fill = 'axis', fillcolor='greenyellow') \
    + text('$ABC = %s$' %n(N_ty_e, digits=3), ((y_e_entra+y_e_sale)/2, max(f)/2), color = 'black')

# Muestra los resultados
show(graf_NTU)

show(r"$H_{c,\text{agot}} = (H_{tG} N_{tG})_\text{agot} = %s \,\text{m} * %s = %s \,\text{m}$" %(numerical_approx(H_ty_a, digits = 3), numerical_approx(N_ty_a, digits = 3), H_c_a))
show(r"$H_{c,\text{enri}} = (H_{tG} N_{tG})_\text{enri} = %s \,\text{m} * %s = %s \,\text{m}$" %(numerical_approx(H_ty_e, digits = 3), numerical_approx(N_ty_e, digits = 3), H_c_e))


Hc,agot=(HtGNtG)agot=0.396m4.60=1.82mH_{c,\text{agot}} = (H_{tG} N_{tG})_\text{agot} = 0.396 \,\text{m} * 4.60 = 1.82 \,\text{m}
Hc,enri=(HtGNtG)enri=0.396m10.1=4.01mH_{c,\text{enri}} = (H_{tG} N_{tG})_\text{enri} = 0.396 \,\text{m} * 10.1 = 4.01 \,\text{m}

Problema N° 2

Valor: 7 puntos

Referencia: Adaptado de Problema 15.D6 Wankat (2008)

Una columna de destilación trabaja a reflujo total, separando una mezcla binaria cuya volatilidad relativa promedio es 2,3. Se obtienen fracciones molares del componente más liviano ysaley_{sale} = 0,95 y yentray_{entra} = 0,05.
  1. Si hay 7,5 m de relleno, determine HtOGH_{tOG} promedio, para ello considere la siguiente expresión del NtOGN_{tOG}.
    • NtOG=11αln(yentra(1ysale)ysale(1yentra))+ln(1yentra1ysale)N_{tOG}= \frac{1}{1- \alpha } \ln \left( \frac{y_{entra}\left(1-y_{sale}\right)}{y_{sale}\left(1-y_{entra}\right)} \right)+ \ln \left( \frac{1-y_{entra}}{1-y_{sale}} \right)
  2. Compruebe los resultados de la parte a mediante el cálculo gráfico usando un diagrama de McCabe-Thiele.
    • NtOG=yentraysale1yydyN_{tOG}= \int_{y_{entra}}^{y_{sale}} \frac{1}{y^{*}-y} \,dy
    • Para generar la curva de equilibrio se emplea la ecuación: y=αx1+(α1)xy^{*}= \frac{ \alpha x}{1+\left( \alpha -1\right)x}

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

%var x

# Datos del enunciado
alfa      =    2.3
y_entra   =    0.05     # mol/mol
y_sale    =    0.95     # mol/mol
h_c       =    7.5       # m

y_eq(x)   =  alfa*x/(1+(alfa-1)*x)
y(x)      =  x

print "Parte a)"

N_tOG     =  1/(1-alfa) * ln((y_entra*(1-y_sale))/(y_sale*(1-y_entra)))+ln((1-y_entra)/(1-y_sale))
H_tOG     =  h_c / N_tOG

# Muestra el resultado formateado, lo que está entre símbolos de dólar se mostrará en formato matemático, para ello se escribe en LaTeX
show (r'$N_{tOG} = \frac{1}{1-\alpha } \ln \left ( \frac{y_{entra} ( 1-y_{sale} )} {y_{sale} (1-y_{entra} )} \right ) + \ln \left (\frac{1-y_{entra}}{1-y_{sale}} \right ) = %s$' %n(N_tOG, digits=4))
show (r'$H_{tOG} = h_c / N_{tOG} = %s$ m' %n(H_tOG, digits=3))

print "Parte b)"

plot(y_eq, x, 0, 1, thickness=2, axes_labels=['$x$','$y$'], title = 'Diagrama de McCabe-Thiele', fontsize=14) \
    + plot(y, x, 0, 1, color = 'black') \
    + line([(0, y_entra), (y_entra, y_entra)]) \
    + text('$y_{entra}$', (-0.05, y_entra), fontsize=16) \
    + line([(0, y_sale),(y_sale, y_sale)]) \
    + text('$y_{sale}$', (-0.05, y_sale), fontsize=16)

f        =  1/(y_eq-y)
N_tOG    =  integral_numerical(f,y_entra, y_sale)[0]
H_tOG    =  h_c / N_tOG  # m

plot(f, y_entra, y_sale, axes_labels=['$y$','$\\frac{1}{(y^*-y)}$'], fill = 'axis', fontsize=14) \
    + text('$y_{entra}$', (y_entra, -5), fontsize=16) \
    + text('$y_{sale}$', (y_sale, -5), fontsize=16) \
    + text('$N_{tOG} = %s$' %n(N_tOG, digits=4), (0.5,50), fontsize=16) \
    + text('$H_{tOG} = %s$ m' %n(H_tOG, digits=3), (0.5,40), fontsize=16)

integral_numerical(f,y_entra, y_sale)


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

Solución

Parte a)
NtOG=11αln(yentra(1ysale)ysale(1yentra))+ln(1yentra1ysale)=7.474N_{tOG} = \frac{1}{1-\alpha } \ln \left ( \frac{y_{entra} ( 1-y_{sale} )} {y_{sale} (1-y_{entra} )} \right ) + \ln \left (\frac{1-y_{entra}}{1-y_{sale}} \right ) = 7.474
HtOG=hc/NtOG=1.00H_{tOG} = h_c / N_{tOG} = 1.00 m
Parte b)
(7.474345100960964, 8.298190025081102e-14)