Ensino de Engenharia Química com Python

Cálculo do tempo e do espaço percorrido por uma esfera

$F=\sum F_i=m\times a \qquad \text{para a partícula a cair na vertical}\qquad m\frac{du}{dt}=\frac{\pi}{6}D_s^3(\rho_s-\rho)g-\frac{1}{2}f\rho u^2A_p\qquad com \qquad A_p=\frac{\pi}{4}D_s^2$

$se \qquad u=\frac{\mu}{D_s\rho}Re\qquad e \qquad \frac{du}{dt}=\frac{\mu}{D_s\rho}\frac{dRe}{dt}\qquad \qquad vem \qquad \qquad\frac{m\mu}{D_s\rho}\frac{dRe}{dt}=\frac{\pi}{6}D_s^3(\rho_s-\rho)g-\frac{\pi}{8}\frac{\mu^2}{\rho}f\ Re^2$

$\text{rearranjando }\frac{8\rho}{\pi\mu^2}\frac{m\mu}{D_s\rho}\frac{dRe}{dt}=\frac{8\rho}{\pi\mu^2}\frac{\pi}{6}D_s^3(\rho_s-\rho)g-f\ Re^2 \qquad \qquad logo \qquad \qquad \frac{8\ m}{\pi\mu\ D_s}\frac{dRe}{dt}=\frac{4}{3\mu^2}\rho\ D_s^3(\rho_s-\rho)g-f\ Re^2$

$\text{como } \phi=\frac{4}{3\mu^2}\rho\ D_s^3(\rho_s-\rho)g \qquad \qquad \text{que é }\frac{4}{3}G_a \qquad \qquad logo \qquad \qquad\frac{8\ m}{\pi\mu\ D_s}\frac{dRe}{dt}=\phi-f\ Re^2$

$$\boxed{t=\frac{8\ m}{\pi\mu\ D_s}\int_{Re_i}^{Re_f} \frac{dRe}{\phi-fRe^2}}$$

como $\qquad d e = u\ dt\qquad $ vem $\qquad de=\frac{8\ m}{\pi\rho\ D_s^2}\frac{\mu}{D_s\rho} Re\frac{dRe}{\phi-fRe^2}\qquad $ isto é:

$$\boxed{e=\frac{2\ m}{\rho\ A_p}\int_{Re_i}^{Re_f} \frac{Re}{\phi-fRe^2}dRe}$$
In [1]:
from __future__ import division
%pylab inline
from scipy.optimize import root
from scipy.integrate import odeint

Populating the interactive namespace from numpy and matplotlib

Qual o tempo e que profundidade atinge uma esfera largada em água demora a atingir a velocidade terminal?

Considere uma partícula esférica de pirite ($\rho_p\ =\ 4840\ kg/m^3$) d = 0,005 m em água à temp. de 20ºC.

In [2]:
# Dados do problema
mu=1.0216e-3
rf=998.9873
rp=4840
dp=0.005
In [3]:
def inte(fg,x0,limites):
    res=[x0]
    f=x0
    for i in range(1,len(limites)):
        dt=limites[i]-limites[i-1]
        f+=fg(limites[i],0)*dt
        res.append(f)
    return res
In [4]:
g=9.8
m=pi/6*dp**3*rp
Ap=pi/4*dp**2

#Cáculo de um
def f(Re):
    if Re>=1000:
        return 0.44
    elif Re>0.2 and Re<1000:
        return 24/Re*(1+0.15*Re**0.687)
    elif Re<=0.2:
        return 24/(Re+1e-16)

fip=4*rf*dp**3*(rp-rf)*g/3/mu**2

def gf(x):
    f(x)*x**2-fip

def v_term(fi):
    if fi>0.44*1000**2:
        return sqrt(3*(rp-rf)*g*dp/rf)
    elif fi<24*0.2:
        return (rp-rf)*g*dp**2/18/mu
    else:
        return root(gf,x0=1)*mu/dp/rf

um=v_term(fip)    
print("A velocidade terminal é um = %.3f m/s "%um)

Rem=um*rf*dp/mu
fi=f(Rem)*Rem**2
Res=linspace(0,Rem*0.999,1000)
us=[i*mu/dp/rf for i in Res]

def dt(Re,t):
    dt=1/abs(fi-f(Re)*Re**2)
    return dt

def de(Re,t):
    de=Re/abs((fi-f(Re)*Re**2))
    return de

ts=inte(dt,0,Res)
ts=[i*8*m/pi/mu/dp for i in ts]
e=inte(de,0,Res)
e=[i*2*m/rf/Ap for i in e]

plot(ts,us,label='$v\ (m/s)$')
plot(ts,e,label='$L\ (m)$')
#arrow(0,um,tm,0,width=0.001,head_width=0.01,linestyle='--')
#arrow(tm,um,0,-um+b,width=0.001,head_width=0.01,linestyle='--')
#arrow(tm,b,-tm,0,width=0.001,head_width=0.01,linestyle='--')
grid(True)
xlabel('$t\ (s)$')
ylabel('$L\ (m)\ or\ \\nu (m/s)$')
legend(bbox_to_anchor=(1.01, 1.01))

print("A partícula percorre %.2f m e demora %.2f s até atingir a sua velocidade terminal."%(e[-1],ts[-1]))

A velocidade terminal é um = 0.752 m/s 
A partícula percorre 0.25 m e demora 0.40 s até atingir a sua velocidade terminal.

Bibliografia

  • Tecnologia Química, Vol. II, Operações Unitárias, J. M. Coulson e J. F. Richardson, Tradução C. Ramalho Carlos, Fundação Calouste Gulbenkian, Lisboa, 2 Edição, 1987.