html('<h2>Solución</h2>')
y_entra = 0.04
y_sale = 0.005
T = 20
p_t = 1
L_s = 70.0
x_entra = 0
G_entra = 60.0
D_c = 0.747
FGa = 0.074
FLa = 0.17
datos_x_eq = [0, 0.0208, 0.0258, 0.0309, 0.0405, 0.0503, 0.0737, 0.096, 0.137, 0.175, 0.210, 0.241, 0.297]
datos_p_NH3 = [0, 12.0 , 15.0, 18.2 , 24.9 , 31.7 , 50.0 , 69.6 , 114 , 166 , 227 , 298 , 470 ]
datos_y_eq = [pA/(p_t*760) for pA in datos_p_NH3]
datos_X_eq = [x/(1-x) for x in datos_x_eq]
datos_Y_eq = [pA/(p_t*760-pA) for pA in datos_p_NH3]
Y_entra = y_entra/(1-y_entra)
Y_sale = y_sale/(1-y_sale)
X_entra = x_entra/(1-x_entra)
i = 0
n_i = len(datos_y_eq)
while i < n_i:
if datos_y_eq[i] >= y_entra:
n_i = i
i += 1
Datos_eq = zip([datos_X_eq[i] for i in range(n_i+1)],[datos_Y_eq[i] for i in range(n_i+1)])
datos_eq = zip([datos_x_eq[i] for i in range(n_i+1)],[datos_y_eq[i] for i in range(n_i+1)])
%var x, X, xI, yI, a0, a1, a2, a3
modelo(x) = a1*x + a2*x^2 + a3*x^3
ajuste = find_fit(datos_eq, modelo, solution_dict=True)
y_eq(x) = modelo(a1=ajuste[a1], a2=ajuste[a2], a3=ajuste[a3])
ajuste = find_fit(Datos_eq, modelo, solution_dict=True)
Y_eq(X) = modelo(a1=ajuste[a1], a2=ajuste[a2], a3=ajuste[a3])
A_c = pi/4 * D_c^2
show (r'$A_c = \frac{\pi}{4} D_c^2 = %s \text{ m}^2$' %(n(A_c, digits=4)))
G_s = G_entra*(1-y_entra)
G_sale = G_s*(1+Y_sale)
Y_op(X) = L_s/G_s * X -L_s/G_s * X_entra + Y_sale
y_op(x) = Y_op(x/(1-x))/(1+Y_op(x/(1-x)))
X_sale = numerical_approx(solve(Y_op(X)==Y_entra, X)[0].rhs())
x_sale = X_sale/(1+X_sale)
print "G_s = %.4g kmol/h" %G_s
print "G_sale = %.4g kmol/h" %G_sale
print "X_sale = %.4g kmol/kmol" %X_sale
print "x_sale = %.4g kmol/kmol" %x_sale
plot(Y_op(X), xmin=X_entra, ymin=Y_sale, xmax=X_sale, ymax=Y_entra, legend_label='Curva de operac.', color = 'red') \
+ plot(Y_eq(X), xmin=0, ymin=0, xmax=datos_X_eq[n_i], ymax=datos_Y_eq[n_i], title = 'Relaciones mol', axes_labels = ['$X_{NH_3}$', '$Y_{NH_3}$'], legend_label='Curva de equilibrio') \
+ scatter_plot(Datos_eq, ymin=0)
x_BM(xI,xL) = ((1-xL)-(1-xI))/ln((1-xL)/(1-xI))
y_BM(yI,yG) = ((1-yG)-(1-yI))/ln((1-yG)/(1-yI))
kpxa(xI,xL) = FLa/x_BM(xI,xL)
kpya(yI,yG) = FGa/y_BM(yI,yG)
n_seg = 4
Delta_y = (y_entra-y_sale)/n_seg
lista_y = [y_sale+Delta_y*i for i in range(n_seg+1)]
lista_x = [find_root(y_op == yop, x_entra, x_sale) for yop in lista_y]
lista_xI = []
lista_yI = []
for i in range(n_seg+1):
xop = lista_x[i]
yop = lista_y[i]
sol_I = solve([yI == -FLa/FGa * xI + FLa/FGa * xop + yop, yI == y_eq(xI)], xI,yI, solution_dict = True)
xI0 = sol_I[0][xI]
yI0 = sol_I[0][yI]
err_x = 1
err_y = 1
tol = 1e-16
it = 0
while it < 10 and (err_x > tol or err_y > tol):
it += 1
sol_I = solve([yI == -kpxa(xI0,xop)/kpya(yI0,yop) * xI + kpxa(xI0,xop)/kpya(yI0,yop) * xop + yop, yI == y_eq(xI)], xI,yI, solution_dict = True)
err_x = abs(xI0-sol_I[0][xI])
err_y = abs(yI0-sol_I[0][yI])
xI0 = sol_I[0][xI]
yI0 = sol_I[0][yI]
lista_xI.append(xI0)
lista_yI.append(yI0)
if(it==10): print "Excesivas iteraciones"
print(lista_x)
print(lista_y)
print(lista_xI)
print(lista_yI)
x_tope = x_entra
x_fondo = x_sale
xI_tope = lista_xI[0]
xI_fondo = lista_xI[n_seg]
y_tope = y_sale
y_fondo = y_entra
yI_tope = lista_yI[0]
yI_fondo = lista_yI[n_seg]
yeq_tope = y_eq(x_tope)
yeq_fondo = y_eq(x_fondo)
xeq_tope = find_root(y_eq == y_tope, 0, 1)
xeq_fondo = find_root(y_eq == y_fondo, 0, 1)
print "***** Parte (a), usando coeficientes individuales *****"
lxyI = line([(lista_x[0], lista_y[0]), (lista_xI[0], lista_yI[0])], color = 'gray')
for i in range(n_seg):
lxyI = lxyI + line([(lista_x[i+1], lista_y[i+1]), (lista_xI[i+1], lista_yI[i+1])], color = 'gray')
show(plot(y_op(x), xmin=x_entra, ymin=y_sale, xmax=x_sale, ymax=y_entra, legend_label='Curva de operac.', color = 'red') \
+ plot(y_eq(x), xmin=0, ymin=0, xmax=datos_x_eq[n_i], ymax=datos_y_eq[n_i], title = 'Fracciones mol', axes_labels = ['$x_{NH_3}$', '$y_{NH_3}$'], legend_label='Curva de equilibrio') \
+ scatter_plot(datos_eq, ymin=0) \
+ lxyI )
f = [1/(lista_y[i]-lista_yI[i]) for i in range(n_seg+1)]
datos_f = zip(lista_y, f)
f_int = spline(datos_f)
integral_f = numerical_integral(f_int, y_sale, y_entra)
corr = 1/2*ln((1-y_sale)/(1-y_entra))
N_tG = integral_f[0]+corr
show(plot(f_int, xmin=y_sale, xmax=y_entra, axes_labels = ['$y_{NH_3}$', '$\\frac{1}{y_{NH_3}-y_{NH_3, I}}$'], fill = 'axis') \
+ scatter_plot(datos_f) )
show (r'$N_{tG} = \int_{y_{sale}}^{y_{entra}}{\frac{1}{y-y_I}} + \frac{1}{2} \ln{\frac{1-y_{entra}}{1-y_{sale}}} = %s + %s = %s$' %(n(integral_f[0], digits=4),n(corr, digits=2),n(N_tG, digits=4)))
H_tG = (G_entra/(3600*kpya(x_fondo,xI_fondo)*A_c) + G_sale/(3600*kpya(x_tope,xI_tope)*A_c))/2
h_c = H_tG * N_tG
show (r'$H_{tG,prom} = %s \text{ m}$' %(n(H_tG, digits=4)))
show (r'$h_{c} = H_{tG} N_{tG} = %s \text{ m}$' %(n(h_c, digits=3)))
print "***** Parte (b), usando coeficientes globales *****"
lxyI = line([(lista_x[0], lista_y[0]), (lista_x[0], y_eq(lista_x[0]))], color = 'gray')
for i in range(n_seg):
lxyI = lxyI + line([(lista_x[i+1], lista_y[i+1]), (lista_x[i+1], y_eq(lista_x[i+1]))], color = 'gray')
show(plot(y_op, xmin=x_entra, ymin=y_sale, xmax=x_sale, ymax=y_entra, legend_label='Curva de operac.', color = 'red') \
+ plot(y_eq(x), xmin=0, ymin=0, xmax=datos_x_eq[n_i], ymax=datos_y_eq[n_i], title = 'Fracciones mol', axes_labels = ['$x_{NH_3}$', '$y_{NH_3}$'], legend_label='Curva de equilibrio') \
+ scatter_plot(datos_eq, ymin=0) \
+ lxyI )
f = [1/(lista_y[i]-y_eq(lista_x[i])) for i in range(n_seg+1)]
datos_f = zip(lista_y, f)
f_int = spline(datos_f)
integral_f = numerical_integral(f_int, y_sale, y_entra)
corr = 1/2*ln((1-y_sale)/(1-y_entra))
N_tOG = integral_f[0]+corr
show( plot(f_int, xmin=y_sale, xmax=y_entra, axes_labels = ['$y_{NH_3}$', '$\\frac{1}{y_{NH_3}-y_{NH_3}^*}$'], fill = 'axis') \
+ scatter_plot(datos_f))
show (r'$N_{tOG} = \int_{y_{sale}}^{y_{entra}}{\frac{1}{y-y^*}} + \frac{1}{2} \ln{\frac{1-y_{entra}}{1-y_{sale}}} = %s + %s = %s$' %(n(integral_f[0], digits=4),n(corr, digits=2),n(N_tOG, digits=4)))
show(r'$\frac{1}{F_{OG}} = \frac{1}{F_G} \frac{(1-y_A)_{iM}}{(1-y_A)_{*M}} + \frac{m\prime(1-x_A)_{iM}}{F_L (1-y_A)_{*M}}$')
m_tope = (yI_tope - yeq_tope) /(xI_tope - x_tope )
m_fondo = (yI_fondo - yeq_fondo)/(xI_fondo - x_fondo)
FOGa_tope = 1/(1/FGa * y_BM(yI_tope, y_tope)/ y_BM(yeq_tope ,y_tope) + m_tope /FLa *x_BM(xI_tope, x_tope) /y_BM(yeq_tope, y_tope) )
FOGa_fondo = 1/(1/FGa * y_BM(yI_fondo,y_fondo)/y_BM(yeq_fondo,y_fondo)+ m_fondo/FLa *x_BM(xI_fondo,x_fondo)/y_BM(yeq_fondo,y_fondo))
show (r'$m\prime_{tope} = %s$' %(n(m_tope, digits=4)))
show (r'$(F_{OG}a)_{tope} = %s$ kmol/(s m³)' %(n(FOGa_tope, digits=4)))
show (r'$m\prime_{fondo} = %s$' %n(m_fondo, digits=4))
show (r'$(F_{OG}a)_{fondo} = %s$ kmol/(s m³)' %n(FOGa_fondo, digits=4))
H_tOG = (G_entra/(3600*FOGa_fondo*A_c) + G_sale/(3600*FOGa_tope*A_c))/2
h_c = H_tOG * N_tOG
show (r'$H_{tOG,prom} = %s \text{ m}$' %(n(H_tOG, digits=4)))
show (r'$h_{c} = H_{tOG} N_{tOG} = %s \text{ m}$' %(n(h_c, digits=3)))