Ejemplo basado en Problema 15-D4 Wankat
Considerando rehervidor parcial
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/D=0.9. Bajo esttas condiciones Hty=0.4 m y Htx=0.24 m para ambas secciones de enriquecimiento y agotamiento. Determine la altura requerida de ambas secciones.
Los datos de equilibrio se encuentran en la Tabla 2-7.
# Datos de equilibrio datos_x = [0.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.0] datos_y = [0.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.0] 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 # Se define el modelo ajuste = find_fit(datos_eq, modelo, solution_dict=True) # Ajuste al modelo y_eq(x) = modelo(a1 = ajuste[a1], a2 = ajuste[a2], a3 = ajuste[a3], a4 = ajuste[a4]) # Sustituye los parámetros ajustados @interact(layout=[['z_F', 'q', 'H_tx_a'], ['x_B', 'R_D', 'n_puntos_agot'], ['x_D','H_ty_a','n_puntos_enri']]) def _(z_F = input_box(default=n(0.40, digits=4), width = 15, label = "$z_F =$"), \ q = input_box(default=n(0.60, digits=4), width = 15, label = "$q =$" ), \ x_B = input_box(default=n(0.04, digits=4), width = 15, label = "$x_B =$"), \ x_D = input_box(default=n(0.92, digits=4), width = 15, label = "$x_D =$"), \ R_D = input_box(default=n(0.90, digits=4), width = 15, label = "$R_D =$" ), \ H_ty_a = input_box(default=n(0.40, digits=4), width = 15, label = "$H_{ty}$/(m) ="), \ H_tx_a = input_box(default=n(0.24, digits=4), width = 15, label = "$H_{tx}$/(m) ="), \ n_puntos_agot = input_box(default=5 , width = 15, label = "Puntos agotamiento"), \ n_puntos_enri = input_box(default=5 , width = 15, label = "Puntos enriquecimiento")): H_ty_e = H_ty_a H_tx_e = H_tx_a # Balances de masa F = 1000 # kmol/h, Base de cálculo D = F*(z_F-x_B)/(x_D-x_B) # Flujo de destilado B = F - D # Flujo de fondos R_B = (q*F + R_D*D - B) / B # Relación de vapor al fondo # 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) # Encuentra punto donde la línea q y la de enriquecimiento se tocan x_aeq = find_root(y_q == y_eq, 0, 1) # Punto donde línea q toca la de 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 # Gráficas de puntos de datos y las diferentes líneas de operación graf_dat = scatter_plot(datos_eq, xmin=0, xmax = 1, ymin = 0, ymax = 1, axes_labels = ["$x_{metanol}$", "$y_{metanol}$"]) graf_l45 = plot(x , x, 0 , 1 , color ="black", aspect_ratio=1) graf_equ = plot(y_eq, x, 0 , 1 , color ="orange" , legend_label = "Equilibrio") graf_ali = plot(y_q , x, x_aeq, z_F , color = "red" , legend_label = "Alim.") graf_enr = plot(y_e , x, x_ali, x_D , color = "greenyellow", legend_label = "Enriq.") graf_ago = plot(y_a , x, x_B , x_ali, color = "green" , legend_label = "Agot." ) # Gráfica de las composiciones alrededor del rehervidor graf_reh = line([(x_B, x_B), (x_B, y_a_entra), (x_a_sale, y_a_entra)]) graf_prob = graf_dat + graf_l45 + graf_equ + graf_ali + graf_enr + graf_ago + graf_reh # Cálculo de condiciones de interfase kx_div_ky_enri = H_ty_e/H_tx_e * R_D/(R_D+1) kx_div_ky_agot = H_ty_a/H_tx_a * (R_B+1)/R_B # ***** Genera valores de datos a lo largo de las composiciones en entrada y salida de la zona de agotamiento Delta_y = (y_a_sale - y_a_entra) / (n_puntos_agot - 1) lista_y = [y_a_entra + Delta_y * i for i in range(n_puntos_agot)] lista_x = [find_root(y_a == yop, x_a_sale-0.0001, x_a_entra+0.0001) for yop in lista_y] lista_xI = [] # Lista vacía lista_yI = [] # Para cada par (x, y) calcula (xI, yI) for i in range(n_puntos_agot): xop = lista_x[i] # Lee un valor de x de la línea de operación yop = lista_y[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_xI.append(raiz_xI) lista_yI.append(raiz_yI) # Genera la gráfica de líneas con pendiente -kx/ky graf_lin = line([(lista_x[0], lista_y[0]), (lista_xI[0], lista_yI[0])], color = 'gray') for i in range(n_puntos_agot-1): graf_lin = graf_lin + line([(lista_x[i+1], lista_y[i+1]), (lista_xI[i+1], lista_yI[i+1])], color = 'gray') # ***** Cálculo del número de unidades de transferencia ***** # Genera lista con valores 1/(x-xI) y luego forma pares con los valores de x f = [1/(lista_x[i]-lista_xI[i]) for i in range(n_puntos_agot)] datos_f = zip(lista_x, f) # Ajuste de los datos a un spline cúbico para integrar fácilmente f_int = spline(datos_f) # Calcula NTU 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 H_c_a = numerical_approx(H_tx_a*N_tx_a, digits = 3) graf_NTU = scatter_plot(datos_f, axes_labels = ['$y$', r'$\frac{1}{y_I-y}$'], 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_y = (y_e_sale - y_e_entra) / (n_puntos_enri - 1) lista_y = [y_e_entra + Delta_y * i for i in range(n_puntos_enri)] lista_x = [find_root(y_e == yop, x_e_sale-0.0001, x_e_entra+0.0001) for yop in lista_y] lista_xI = [] lista_yI = [] for i in range(n_puntos_enri): xop = lista_x[i] yop = lista_y[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_xI.append(raiz_xI) lista_yI.append(raiz_yI) for i in range(n_puntos_enri): graf_lin = graf_lin + line([(lista_x[i], lista_y[i]), (lista_xI[i], lista_yI[i])], color = 'darkgray') f = [1/(lista_x[i]-lista_xI[i]) for i in range(n_puntos_enri)] datos_f = zip(lista_x, f) f_int = spline(datos_f) integral_f = numerical_integral(f_int, x_e_sale, x_e_entra) N_tx_e = integral_f[0] H_c_e = numerical_approx(H_tx_e*N_tx_e, digits = 3) graf_NTU = graf_NTU + scatter_plot(datos_f, axes_labels = ['$x$', r'$\frac{1}{x-x_I}$'], facecolor = 'none') \ + plot(f_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(f)/2), color = 'black') # Muestra los resultados ga = graphics_array([graf_prob + graf_lin, graf_NTU]) ga.show(figsize=[12,5]) show(r"$H_{c,\text{agot}} = (H_{tL} N_{tL})_\text{agot} = %.3g \,\text{m} * %.3g = %.3g \,\text{m}$" %(H_ty_a, N_tx_a, H_c_a)) show(r"$H_{c,\text{enri}} = (H_{tL} N_{tL})_\text{enri} = %.3g \,\text{m} * %.3g = %.3g \,\text{m}$" %(H_ty_e, N_tx_e, H_c_e)) # ***** Final del código *****
Interact: please open in CoCalc