print "=== Problema N° 1 ==="
# A continuación se muestran los datos de la isoterma de adsorción de glucosa (en agua) utilizando alúmina.

c = [0.004, 0.0087, 0.019, 0.024, 0.094, 0.195] # g/cm³, kg/L
q = [0.026, 0.053,  0.075, 0.082, 0.123, 0.129] # g/g, kg/kg
print "c en kg/L y q en kg/kg"

# q = m_glucosa/m_alúmina
#
# a. Ajuste a el modelo de Langmuir y Freundlich y justifique cuál de los dos modelos es mejor
# b. Suponga que va a realizar la purificación de 15 L de disolución 1,39 mol/L de glucosa utilizando una etapa con 25 kg de alúmina pura. ¿Cuál es la concentración final de la disolución?
# c. ¿Cuál es la concentración final de la disolución si utiliza 2 etapas a flujo cruzado, en ambas utilizando 25 kg de alumina pura?
# d. ¿Cuál es la cantidad de alúmina en cada etapa que minimiza la cantidad de carbón de la configuración anterior y cuánta alúmina se ahorraría usando el mínimo?

datos_eq = zip(c,q)            # Combina los datos de equilibrio
%var c_A, q_A_max, K_Ac, n

L       =   15         # [L]
c0      =    1.39      # [mol/L]
M       =  180.1559    # [kg/kmol]
S1      =   25         # [kg]
q0      =    0         # [kg/kg]
S2      =   25         # [kg]

c0      =  c0*M/1000   # [kg/L]

print "parte a) Ajustar a el modelo de Langmuir y Freundlich"

show("===== Isoterma de Langmuir =====")
qA_L(c_A) = q_A_max*K_Ac*c_A / (1+K_Ac*c_A)
ajuste    = find_fit(datos_eq, qA_L, solution_dict = True)
qA_L      = qA_L(q_A_max= numerical_approx(ajuste[q_A_max], digits = 4), K_Ac= numerical_approx(ajuste[K_Ac], digits=4))

show(r"$K_{A,c} = $ ", numerical_approx(ajuste[K_Ac], digits=4), r" L/kg; $q_{A,máx} = $",  numerical_approx(ajuste[q_A_max], digits=4), " kg/kg")
show("$q_{A_L} = $", qA_L)

show("===== Isoterma de Freundlich =====")
qA_F(c_A) = K_Ac*c_A^n
ajuste    = find_fit(datos_eq, qA_F, solution_dict = True)
qA_F      = qA_F(n= numerical_approx(ajuste[n], digits=4), K_Ac=numerical_approx(ajuste[K_Ac], digits=4))

show(r"$K_{A,c} = $ ", numerical_approx(ajuste[K_Ac], digits=4), r" L/kg; $n = $",  numerical_approx(ajuste[n], digits=4))
show("$q_{A_F} = $", qA_F)

# Gráfica
gr_datos  =  scatter_plot(datos_eq, axes_labels = ['$c_A$', '$q_A$'])
gr_Langm  =  plot(qA_L, xmin = 0, xmax = max(max(c),c0), legend_label = 'Langmuir', color = 'blue')
gr_Frlch  =  plot(qA_F, xmin = 0, xmax = max(c), legend_label = 'Freundlich', color = 'green')

show(gr_datos+gr_Langm+gr_Frlch)


print "parte b) Una etapa"
def Calcula_etapa(W_L, W_S, c_entra, q_entra, q_isot):
    R.<c_sale> = QQ[]
    q_sale     = q_isot(c_sale)
    f          = c_sale - (W_S/W_L * (q_entra - q_sale) + c_entra)
    csale      = find_root(f, 0, c_entra)
    return [csale, q_sale(csale)]

c1, q1 = Calcula_etapa(L, S1, c0, q0, qA_L)
show (r'$c_1 = %s$ kg/L; $q_1 = %s$ kg/kg' %(numerical_approx(c1, digits = 3), numerical_approx(q1, digits=3)))

print "parte c) Dos etapas a flujo cruzado, idénticas cantidades de sólido"
c2, q2 = Calcula_etapa(L, S2, c1, q0, qA_L)
show (r'$c_2 = %s$ kg/L; $q_2 = %s$ kg/kg' %(numerical_approx(c2, digits = 3), numerical_approx(q2, digits=3)))

gr_etapa1 = line([(c0,q0),(c1,q1),(c1,0)], color = 'purple', legend_label = 'Etapa 1')
gr_etapa2 = line([(c1,q0),(c2,q2),(c2,0)], color = 'magenta', legend_label = 'Etapa 2')
show(gr_datos+gr_Langm+gr_etapa1+gr_etapa2)

print "parte d) Dos etapas a flujo cruzado, calcular sólido mínimo"
S1m = L*(c0-c_A)/(qA_L(c_A) + q0)    # Alúmina requerida en la etapa 1 para obtnener una concentración c_A
S2m = L*(c_A-c2)/(qA_L(c2) + q0)     # Alúmina requerida en la etapa 2 con una concentración de entrada c_A
STm = S1m+S2m     # Alúmina total
S_Tm, c1 = find_local_minimum(STm, c2, c0)
show (r"$S_{T}$ /kg $=$", STm)
show (r"$S_{T, mín}= %s$ kg; $c_{1}= %s$ kg/L" %(numerical_approx(S_Tm, digits=3), numerical_approx(c1, digits=3)))
show (r"$Ahorro = %s $ kg" %(numerical_approx(((S1+S2)-S_Tm), digits=3)))
show (r"$S_{1}= %s$ kg; $S_{2}= %s$ kg" %(numerical_approx(S1m(c1), digits=3), numerical_approx(S2m(c1), digits=3)))

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


=== Problema N° 1 === c en kg/L y q en kg/kg parte a) Ajustar a el modelo de Langmuir y Freundlich
===== Isoterma de Langmuir =====
KA,c=K_{A,c} = 61.43\displaystyle 61.43 L/kg; qA,maˊx=q_{A,máx} = 0.1413\displaystyle 0.1413 kg/kg
qAL=q_{A_L} = 8.678cA61.43cA+1\displaystyle \frac{8.678 \, c_{A}}{61.43 \, c_{A} + 1}
===== Isoterma de Freundlich =====
KA,c=K_{A,c} = 0.2282\displaystyle 0.2282 L/kg; n=n = 0.3019\displaystyle 0.3019
qAF=q_{A_F} = 0.2282cA0.3019\displaystyle 0.2282 \, c_{A}^{0.3019}
parte b) Una etapa
<string>:1: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...) See http://trac.sagemath.org/5930 for details.
c1=0.0632c_1 = 0.0632 kg/L; q1=0.112q_1 = 0.112 kg/kg
parte c) Dos etapas a flujo cruzado, idénticas cantidades de sólido
c2=0.00531c_2 = 0.00531 kg/L; q2=0.0347q_2 = 0.0347 kg/kg
parte d) Dos etapas a flujo cruzado, calcular sólido mínimo
/projects/sage/sage-7.3/local/lib/python2.7/site-packages/smc_sagews/sage_server.py:957: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...) See http://trac.sagemath.org/5930 for details. exec compile(block+'\n', '', 'single') in namespace, locals
STS_{T} /kg == 431.9cA+0.1152(61.43cA+1)(15cA+3.75625051500000)cA2.292\displaystyle 431.9 \, c_{A} + \frac{0.1152 \, {\left(61.43 \, c_{A} + 1\right)} {\left(-15 \, c_{A} + 3.75625051500000\right)}}{c_{A}} - 2.292
ST,mıˊn=46.3S_{T, mín}= 46.3 kg; c1=0.0365c_{1}= 0.0365 kg/L
Ahorro=3.68Ahorro = 3.68 kg
S1=32.9S_{1}= 32.9 kg; S2=13.5S_{2}= 13.5 kg
print "Problema N° 2"
# Referencia: Problema 17.D3 Wankat (2008)
#
# Una columna cromatográfica está rellena con una resina de intercambio iónico en su forma de calcio. La columna mide 0,75 m y la velocidad superficial es de 15 cm/min. A la columna se le alimenta una disolución de glucosa de 75 g/L, después de 2 min la concentración cambia a 100 g/L de glucosa. Encuentre el perfil de concentración de salida (idealizado) de la columna.

cA    = var('cA')
L     =  75      # [cm]
u_s   =  15      # [cm/min]
t1    =   0      # [min]
c1    =  75      # [g/L]
t2    =   2      # [min]
c2    = 100      # [g/L]

# Datos adicionales
eps_p = 0
eps_e = 0.4
rho_S = 2        # [kg/L]
qA    = 0.51*cA
Kd    = 1        # supuesto
t0    = 0        # [min]
c0    = 0        # [g/L]

print "Usando modelo para onda de choque"
DqDc (c_antes, c_desp) = (qA(c_desp)-qA(c_antes))/(c_desp-c_antes)

u_inter = u_s / eps_e
uA(c_antes,c_desp) = u_inter/(1 + (1-eps_e)/eps_e * Kd + (1-eps_e)*(1-eps_p)/eps_e * rho_S * DqDc(c_antes,c_desp))

u_sh1 = uA(c0, c1)
u_sh2 = uA(c1, c2)
t_s1 = t1 + L/u_sh1
t_s2 = t2 + L/u_sh2

line([(t0, c0),(t1,c0),(t1,c1),(t2,c1),(t2,c2),(t_s2,c2)], xmax = t_s2+2, axes_labels=['$t$/min','$c_F$/(g/L)'])

show ("$u_{sh1} = $", numerical_approx(u_sh1, digits=4), " cm/min; $u_{sh2} = $", numerical_approx(u_sh2, digits=4), " cm/min")

if(t_s1 <= t_s2):
    show("Las ondas de choque no coinciden")
    gr_onda1 = line([(t0,0),(t1,0),(t_s1, L)], legend_label='Onda de choque 1')
    gr_onda2 = line([(t0,0),(t2,0),(t_s2, L)], legend_label='Onda de choque 2', axes_labels=['$t$/min','$L$/cm'], color = 'green', xmax = t_s2+2,)
    gr_perfil = line([(t0,c0),(t_s1,c0),(t_s1,c1)]) + line([(t_s1,c1),(t_s2,c1),(t_s2,c2),(t_s2+2,c2)], axes_labels=['$t/\mathrm{min}$','$c_{sale}/\mathrm{(g/L)}$'])
    show(gr_onda1+gr_onda2)
else:
    show("Las ondas de choque coinciden y se genera una tercera")
    t_inter = t2*u_sh2/(u_sh2-u_sh1)
    z_inter = u_sh1*t_inter
    u_sh3 = uA(c0, c2)
    show ("$u_{sh3} = $", numerical_approx(u_sh3, digits=4), " cm/min; $t_{inter} = $", numerical_approx(t_inter, digits=4), " min; $z_{inter} = $", numerical_approx(z_inter, digits=4), " cm")
    t_s3 = t_inter+(L-z_inter)/u_sh3
    gr_onda1 = line([(t0,0),(t1,0),(t_inter, z_inter)], legend_label='Onda de choque 1')
    gr_onda2 = line([(t0,0),(t2,0),(t_inter, z_inter)], legend_label='Onda de choque 2', color = 'green')
    gr_onda3 = line([(t_inter,z_inter),(t_s3, L)], legend_label='Onda de choque 3', color = 'red', axes_labels=['$t$/min','$L$/cm'], xmax = t_s3+2)
    gr_perfil = line([(t0,c0),(t_s3,c0),(t_s3,c1),(t_s3+2,c1)], axes_labels=['$t/\mathrm{min}$','$c_{sale}/\mathrm{(g/L)}$'])
    show(gr_onda1+gr_onda2+gr_onda3)

show(gr_perfil)
Problema N° 2 Usando modelo para onda de choque
ush1=u_{sh1} = 9.305\displaystyle 9.305 cm/min; ush2=u_{sh2} = 9.305\displaystyle 9.305 cm/min
Las ondas de choque no coinciden