html('<h2>Solución</h2>')
cA = var('cA')
F = 0.30/1000
Dc = 6/100
Ac = pi/4*Dc^2
Hc = 80
t0 = 0
c0 = 0
t1 = 5
c1 = 0.01
t2 = t1+10
c2 = 0.015
t3 = t2+15
c3 = 0
t4 = t3+60
u_s = F/Ac*100
print("u_s = %.3g cm/min"%u_s)
eps_e = 0.4
eps_p = 0.1
Kd = 1.0
qA(cA) = 20*cA/(1+350*cA)
rho_f = 0.78
rho_b = 1.2
rho_p = (rho_b - eps_e*rho_f)/(1-eps_e)
rho_S = (rho_p - eps_p*rho_f)/(1-eps_p)
print("Densidad de la partícula, rho_p = %.3g kg/L"%rho_p)
print("Densidad del sólido, rho_S = %.3g kg/L"%rho_S)
line([(t0, c0),(t1,c0),(t1,c1),(t2,c1),(t2,c2),(t3,c2),(t3,c3),(t4,c3)], xmin = t0, axes_labels=['$t$/min','$c_F$/(g/L)'])
print("Usando modelo para onda de choque")
DqDc (c_antes, c_desp) = (qA(cA = c_desp)-qA(cA = c_antes))/(c_desp-c_antes)
u_inter = u_s / eps_e
uA_sh(c_antes,c_desp) = u_inter/(1 + (1-eps_e)/eps_e * eps_p * Kd + (1-eps_e)*(1-eps_p)/eps_e * rho_S * DqDc(c_antes,c_desp))
u_sh1 = uA_sh(c0, c1)
u_sh2 = uA_sh(c1, c2)
t_s1 = t1 + Hc/u_sh1
t_s2 = t2 + Hc/u_sh2
t_s3 = 0
print("Tiempo que tardaría en recorrer Hc onda de choque 1, Delta t_sh1 = %.3g min" %(Hc/u_sh1))
print("Tiempo que tardaría en recorrer Hc onda de choque 2, Delta t_sh2 = %.3g min" %(Hc/u_sh2))
print("Momento en que saldría onda de choque 1, t_sh1 = %.3g min" %(t_s1))
print("Momento en que saldría onda de choque 2, t_sh2 = %.3g min" %(t_s2))
[t_s1, t_s2, t_s3]
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, Hc)], legend_label='Onda de choque 1')
gr_onda2 = line([(t0,0),(t2,0),(t_s2, Hc)], 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$/min','$c_{sale}$/(g/L)'])
show(gr_onda1+gr_onda2)
show ("$t_{esc1} = $", numerical_approx(t_s1, digits=3), " min; $t_{esc2} = $", numerical_approx(t_s2, digits=3), " min")
else:
show("Las ondas de choque coinciden y se genera una tercera")
t_inter = t1+(t2-t1)*u_sh2/(u_sh2-u_sh1)
z_inter = u_sh1*(t_inter-t1)
u_sh3 = uA_sh(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+(Hc-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, Hc)], 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,c2)], axes_labels=['$t$/min','$c_{sale}$/(g/L)'])
show(gr_onda1+gr_onda2+gr_onda3)
print("Tiempo salida onda de choque 3 = %.3g min"%t_s3)
print("Onda difusa")
conc = [c2+i*(c3-c2)/9 for i in range(10)]
u_A(cA) = u_inter/(1 + (1-eps_e)/eps_e * eps_p * Kd + (1-eps_e)*(1-eps_p)/eps_e * rho_S * derivative(qA))
t_esc = [t3+numerical_approx(Hc/u_A(c), digits=4) for c in conc]
print("Tiempos de salida en minutos:", t_esc)
t_max = t_esc[5]
lista_tc = list(zip(t_esc,conc))
lista_tc.insert(0,(max(t_s2,t_s3),c2))
gr_perfil + line(lista_tc, xmax = t4)