Em 1873, o físico holandês Johannes Diderik van der Waals defendeu uma tese de doutorado que resolvia um problema incômodo: a lei dos gases ideais, \(PV=nRT\), simplesmente não descrevia gases reais perto da condensação nem explicava por que, sob pressão e temperatura suficientes, um gás vira líquido. Van der Waals propôs uma correção com dois termos extras: um para o volume que as próprias moléculas ocupam, outro para a atração entre elas. A equação resultante fez algo que a lei dos gases ideais jamais conseguiria: previu a existência de um ponto crítico, a temperatura acima da qual líquido e vapor deixam de ser fases distintas. A descoberta valeu a van der Waals o Prêmio Nobel de Física de 1910.
O preço da correção: a elegante equação linear \(PV=nRT\) vira uma cúbica em \(V\),
$$ \left(P + \frac{an^2}{V^2}\right)(V - nb) = nRT $$Só que cúbicas, na prática, não têm solução fechada útil. Essa é só uma das equações assim fora da física: da química analítica à engenharia de tubulações, o mesmo tipo de problema aparece sob disfarces diferentes.
Antes de voltar a van der Waals, porém, vamos começar por um aquecimento: uma equação quadrática, que todo mundo sabe resolver, e onde a pergunta interessante é “quando vale a pena”, e não “como resolver”. Depois vem a cúbica do gancho, e por fim uma equação que não é polinômio nenhum.
Por Que Isso Importa na Química e na Engenharia? #
- Engenharia química: equações de estado de gases reais (como a de van der Waals) são a base de qualquer projeto de processo industrial que envolva compressão, liquefação ou armazenamento de gases.
- Química analítica: calcular o pH real de uma solução (sem os atalhos simplificadores de sala de aula) importa sempre que a diluição é alta ou a precisão exigida é maior que a de um exercício de livro-texto.
- Engenharia civil e mecânica: o dimensionamento de qualquer tubulação
(água, óleo, gás) depende do fator de atrito, que por décadas exigiu
iteração manual ou leitura de ábaco antes de existir
scipy.optimize. - Indústria de processos em geral: as três equações deste artigo aparecem, das mais variadas formas, em qualquer software comercial de simulação de processos químicos.
Os métodos usados aqui, brentq para busca de raiz e numpy.roots para
polinômios, são os mesmos apresentados em Resolvendo Equações Numericamente com SciPy,
o primeiro artigo deste bloco. Como no artigo anterior, sobre física,
aqui eles resolvem problemas completos de química e engenharia, com
números reais.
1. pH de um Ácido Fraco: Quando o Atalho de Sala de Aula Basta #
Todo curso de química geral ensina a mesma receita para calcular o pH de um ácido fraco. O equilíbrio de dissociação,
$$ HA \rightleftharpoons H^+ + A^-, \qquad K_a = \frac{x^2}{C_0 - x} $$vira uma equação quadrática em \(x = [H^+]\). A receita de sala de aula assume que a dissociação é pequena (\(C_0 - x \approx C_0\)) e chega direto a \(x \approx \sqrt{K_a C_0}\), sem resolver quadrática nenhuma. É rápido e funciona na maioria dos exercícios. Mas é uma aproximação, e vale ver o quanto ela se afasta da resposta exata.
Exemplo 1: ácido acético, exato x aproximado
Enunciado: Uma solução de ácido acético (\(K_a = 1{,}8 \times 10^{-5}\), a 25 °C) tem concentração inicial \(C_0 = 0{,}10\) mol/L. Qual é o pH exato, resolvendo a quadrática completa, e qual é o pH pela aproximação de química geral?
Solução:
import numpy as np
from scipy import optimize
Ka = 1.8e-5 # ácido acético
C0 = 0.10 # mol/L
# x^2 + Ka*x - Ka*C0 = 0
roots_acid = np.roots([1, Ka, -Ka * C0])
# só as raízes reais e positivas: [H+] não pode ser negativa nem complexa
x_exact = roots_acid[np.isreal(roots_acid) & (roots_acid.real > 0)].real.max()
x_approx = np.sqrt(Ka * C0)
pH_exact = -np.log10(x_exact)
pH_approx = -np.log10(x_approx)
print(f"[H+] exato: {x_exact:.6e} mol/L -> pH = {pH_exact:.4f}")
print(f"[H+] aproximado: {x_approx:.6e} mol/L -> pH = {pH_approx:.4f}")
print(f"diferença de pH: {abs(pH_exact - pH_approx):.5f}")[H+] exato: 1.332671e-03 mol/L -> pH = 2.8753
[H+] aproximado: 1.341641e-03 mol/L -> pH = 2.8724
diferença de pH: 0.00291Para essa concentração, a diferença é minúscula: menos de \(0{,}003\) unidade de pH, irrelevante na prática. O gráfico a seguir mostra por quê: a curva é a quadrática completa, \(g(x) = x^2 + K_a x - K_a C_0\), com a raiz exata marcada onde ela cruza zero. A aproximação não resolve essa mesma equação: é a raiz de uma equação mais simples, então ela não cai exatamente sobre a curva, mas fica visivelmente perto.
import matplotlib.pyplot as plt
import numpy as np
Ka = 1.8e-5 # ácido acético
C0 = 0.10 # mol/L
roots_acid = np.roots([1, Ka, -Ka * C0])
x_exact = roots_acid[np.isreal(roots_acid) & (roots_acid.real > 0)].real.max()
x_approx = np.sqrt(Ka * C0)
x_range = np.linspace(0, x_exact * 1.5, 400)
def ph_residual(x, Ka, C0):
return x**2 + Ka * x - Ka * C0
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(x_range * 1e3, ph_residual(x_range, Ka, C0), color="C0")
ax.axhline(0, color="gray", linewidth=1, zorder=0)
ax.scatter([x_exact * 1e3], [0], color="C2", zorder=3,
label=f"raiz exata: {x_exact * 1e3:.3f} mmol/L")
ax.axvline(x_approx * 1e3, color="C1", linestyle="--",
label=f"aproximação: {x_approx * 1e3:.3f} mmol/L")
ax.legend()
ax.set_xlabel("[H+] (mmol/L)")
ax.set_ylabel("resíduo")
ax.set_title("pH de ácido fraco: raiz exata x aproximação de química geral")
Quando a Aproximação Deixa de Valer? #
A aproximação \(C_0 - x \approx C_0\) só é boa quando a dissociação é pequena perto de \(C_0\), o que falha justamente em soluções muito diluídas, onde uma fração maior do ácido se dissocia. Varrendo várias concentrações em escala logarítmica (natural em química, onde concentrações variam ordens de grandeza) e comparando pH exato com aproximado, fica visível onde a regra prática de sala de aula para de ser confiável.
import matplotlib.pyplot as plt
import numpy as np
Ka = 1.8e-5 # ácido acético
C0_range = np.logspace(-6, 0, 200) # 1 µmol/L até 1 mol/L
pH_exact_range = []
pH_approx_range = []
for C0_i in C0_range:
roots_i = np.roots([1, Ka, -Ka * C0_i])
x_exact_i = roots_i[np.isreal(roots_i) & (roots_i.real > 0)].real.max()
x_approx_i = np.sqrt(Ka * C0_i)
pH_exact_range.append(-np.log10(x_exact_i))
pH_approx_range.append(-np.log10(x_approx_i))
erro_pH = np.abs(np.array(pH_exact_range) - np.array(pH_approx_range))
fig, ax = plt.subplots(figsize=(8, 5))
ax.semilogx(C0_range, erro_pH, color="C0")
ax.axhline(0.05, color="C1", linestyle="--", label="erro de 0,05 unidade de pH")
ax.set_xlabel(r"$C_0$ (mol/L)")
ax.set_ylabel("|pH exato − pH aproximado|")
ax.set_title("Onde a aproximação de química geral deixa de valer")
ax.legend()
Abaixo de aproximadamente \(3 \times 10^{-4}\) mol/L, o erro já passa de \(0{,}05\) unidade de pH e continua crescendo. É exatamente a faixa de soluções muito diluídas (comuns em análises ambientais e traços), onde a receita de sala de aula deixa de ser confiável e a quadrática completa passa a valer a pena.
Este artigo vem chamando a solução da quadrática de “exata”, em oposição ao atalho de sala de aula. Ela é exata dentro do seu modelo, e esse modelo ignora uma coisa: a autoprotólise da água, \(2H_2O \rightleftharpoons H_3O^+ + OH^-\), que também produz \(H^+\).
O modelo completo junta o equilíbrio do ácido, o da água (\(K_w = 10^{-14}\)), o balanço de massa e o balanço de cargas \([H^+] = [A^-] + [OH^-]\). O resultado é uma cúbica:
$$[H^+]^3 + K_a[H^+]^2 - (K_aC_0 + K_w)[H^+] - K_aK_w = 0$$import numpy as np
Ka = 1.8e-5 # ácido acético
Kw = 1e-14
for C0_i in (1e-1, 1e-4, 1e-6, 1e-7):
x_quad = np.roots([1, Ka, -Ka * C0_i]).max()
r = np.roots([1, Ka, -(Ka * C0_i + Kw), -Ka * Kw])
x_completo = max(r[np.abs(r.imag) < 1e-12].real)
print(f"C0={C0_i:.0e}: pH quadrática={-np.log10(x_quad):.4f}"
f" pH com água={-np.log10(x_completo):.4f}")C0=1e-01: pH quadrática=2.8753 pH com água=2.8753
C0=1e-04: pH quadrática=4.4638 pH com água=4.4638
C0=1e-06: pH quadrática=6.0223 pH com água=6.0178
C0=1e-07: pH quadrática=7.0024 pH com água=6.7927Para o ácido acético a \(0{,}10\) mol/L, nada muda: a água é irrelevante ali. Mas repare no último caso, \(C_0 = 10^{-7}\) mol/L: a quadrática devolve pH \(7{,}0024\), ou seja, prevê que uma solução de ácido é levemente básica. É um absurdo físico, e vem justamente de ignorar a água, que, nessa diluição, contribui com mais \(H^+\) do que o próprio ácido. O modelo completo dá \(6{,}79\), corretamente ácido.
A lição vale para o artigo inteiro, e é o motivo de esta callout existir: “exato” é sempre exato em relação a um modelo. Descer um degrau de aproximação (do atalho para a quadrática) resolve um problema e revela o seguinte. O gráfico acima mede a distância entre o degrau 1 e o degrau 2; abaixo de \(10^{-6}\) mol/L, seria preciso o degrau 3.
2. Van der Waals: Voltando ao Gancho #
De volta à equação que abriu este artigo. O volume molar de um gás real
não tem solução fechada, mas brentq resolve sem dificuldade, desde que
o bracket contenha a raiz física certa.
Exemplo 2: CO2 real x CO2 ideal
Enunciado: Um mol de CO2 está a \(P = 10\) atm e \(T = 300\) K. Qual é o volume molar pela lei dos gases ideais, e qual é o volume real pela equação de van der Waals (\(a = 3{,}640\) L²·atm/mol², \(b = 0{,}04267\) L/mol)?
Solução:
from scipy import optimize
R = 0.082057366080960 # atm L / (mol K)
a_co2, b_co2 = 3.640, 0.04267 # L^2 atm / mol^2 ; L / mol
P_val, T_val, n_val = 10.0, 300.0, 1.0
def van_der_waals(V):
return (P_val + a_co2 * n_val**2 / V**2) * (V - n_val * b_co2) - n_val * R * T_val
V_ideal = n_val * R * T_val / P_val
V_vdw = optimize.brentq(van_der_waals, 0.01, 10)
print("V ideal:", V_ideal, "L")
print("V van der Waals:", V_vdw, "L")
print("desvio:", (V_vdw - V_ideal) / V_ideal * 100, "%")V ideal: 2.4617209824287998 L
V van der Waals: 2.3524663349745696 L
desvio: -4.438140968617679 %O CO2 real ocupa cerca de 4,4% menos volume que o gás ideal previria:
nessas condições, a atração entre moléculas (o termo \(a\)) supera
ligeiramente o efeito do volume próprio delas (o termo \(b\)). O gráfico
a seguir mostra as duas soluções: a linha tracejada é o volume ideal, o
ponto é a raiz real encontrada pelo brentq.
import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize
R = 0.082057366080960 # atm L / (mol K)
a_co2, b_co2 = 3.640, 0.04267 # L^2 atm / mol^2 ; L / mol
P_val, T_val, n_val = 10.0, 300.0, 1.0
def van_der_waals(V):
return (P_val + a_co2 * n_val**2 / V**2) * (V - n_val * b_co2) - n_val * R * T_val
V_ideal = n_val * R * T_val / P_val
V_vdw = optimize.brentq(van_der_waals, 0.01, 10)
V_range = np.linspace(0.05, 3, 500)
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(V_range, van_der_waals(V_range), color="C0")
ax.axhline(0, color="gray", linewidth=1, zorder=0)
ax.axvline(V_ideal, color="C1", linestyle="--", label="V ideal")
ax.scatter([V_vdw], [0], color="C2", zorder=3, label="V van der Waals")
ax.legend()
ax.set_xlabel("V (L)")
ax.set_ylabel("resíduo de van der Waals")
ax.set_title("Gás ideal x gás real (CO₂)")
E se a Cúbica Tiver Mais de Uma Raiz Real? #
brentq só encontra a raiz que você já esperava
No Exemplo 2, brentq assumiu, silenciosamente, que só havia uma raiz
física no bracket (0.01, 10). Mas a equação de van der Waals é uma
cúbica em \(V\) de verdade, e cúbicas podem ter até 3 raízes reais.
brentq nunca avisa que outras raízes existem fora do bracket fornecido;
ele simplesmente devolve a primeira que encontra dentro dele.
Para achar todas de uma vez, é preciso primeiro reconhecer a cúbica que está escondida na forma usual da equação. A manobra é multiplicar tudo por \(V^2\) e distribuir:
$$ \left(P + \frac{an^2}{V^2}\right)(V - nb) = nRT $$$$ PV - Pnb + \frac{an^2}{V} - \frac{an^3b}{V^2} = nRT $$Multiplicando por \(V^2\) e agrupando por potência de \(V\):
$$ PV^3 - (Pnb + nRT)V^2 + an^2V - an^3b = 0 $$Agora é um polinômio de grau 3 com coeficientes explícitos: exatamente o
formato que o numpy.roots espera.
As três raízes reais aparecem abaixo da temperatura crítica \(T_c\), exatamente a região
do diagrama P-V onde líquido e vapor coexistem: o “laço” nos diagramas de
fase de gases reais, e a mesma descoberta que valeu o Nobel a van der
Waals. Para achar todas as raízes de uma vez, reescrevemos a equação na
forma polinomial explícita,
\(PV^3 - (Pnb+nRT)V^2 + an^2V - an^3b = 0\), e usamos de novo o mesmo
numpy.roots da seção anterior, a ferramenta que o
próximo artigo desta série
vai explorar em profundidade.
Em vez de escolher \(T\) e \(P\) arbitrariamente até a cúbica render 3 raízes, usamos a lei dos estados correspondentes:
Definindo a temperatura e a pressão reduzidas \(T_r = T/T_c\) e \(P_r = P/P_c\) (em relação aos valores críticos \(T_c\) e \(P_c\) da substância), qualquer fluido descrito pela equação de van der Waals com \(T_r < 1\) tem uma faixa de pressões onde a cúbica admite 3 raízes reais, independentemente da substância. Isso é o que dá utilidade prática à lei.
A ressalva importa: essa é uma propriedade da equação, não um fato sobre gases reais. Gases reais só obedecem à lei dos estados correspondentes de forma aproximada, e a própria van der Waals é uma equação de estado bastante grosseira perto do ponto crítico; é por isso que existem Redlich-Kwong, Soave e Peng-Robinson, todas mais precisas e todas ainda cúbicas em \(V\), ou seja, resolvidas exatamente do mesmo jeito.
import numpy as np
R = 0.082057366080960 # atm L / (mol K)
a_co2, b_co2 = 3.640, 0.04267 # L^2 atm / mol^2 ; L / mol
T_c = 8 * a_co2 / (27 * R * b_co2)
P_c = a_co2 / (27 * b_co2**2)
print(f"CO2: T_c = {T_c:.1f} K, P_c = {P_c:.1f} atm")
def van_der_waals_coeffs(P, T, n=1.0):
c3 = P
c2 = -(P * n * b_co2 + n * R * T)
c1 = a_co2 * n**2
c0 = -a_co2 * n**3 * b_co2
return [c3, c2, c1, c0]
T_r, P_r = 0.84, 0.41
T_subcritica = T_r * T_c
P_subcritica = P_r * P_c
print(f"T_r = {T_r} -> T = {T_subcritica:.1f} K")
print(f"P_r = {P_r} -> P = {P_subcritica:.1f} atm")
raizes = np.roots(van_der_waals_coeffs(P_subcritica, T_subcritica))
print("raízes (podem incluir complexas):", raizes)
raizes_reais_positivas = np.sort(raizes[(np.abs(raizes.imag) < 1e-8) & (raizes.real > 0)].real)
print("raízes reais positivas, ordenadas:", raizes_reais_positivas)CO2: T_c = 308.0 K, P_c = 74.0 atm
T_r = 0.84 -> T = 258.7 K
P_r = 0.41 -> P = 30.4 atm
raízes (podem incluir complexas): [0.53623521 0.13527776 0.07052874]
raízes reais positivas, ordenadas: [0.07052874 0.13527776 0.53623521]Três raízes reais positivas. Fisicamente:
- a menor, \(V \approx 0{,}0705\) L/mol, é o volume molar do líquido (perto do covolume \(b\), moléculas bem compactadas);
- a maior, \(V \approx 0{,}5362\) L/mol, é o volume molar do vapor;
- a do meio, \(V \approx 0{,}1353\) L/mol, não corresponde a nenhum estado físico estável: é a parte da curva onde \(\partial P/\partial V > 0\), termodinamicamente instável.
Dizer que as três raízes marcam a região onde líquido e vapor coexistem é uma simplificação conveniente, mas não é bem isso. Ter três raízes reais significa estar dentro do laço de van der Waals, e o laço é mais largo que a coexistência de verdade.
A coexistência em equilíbrio só acontece numa pressão específica para cada temperatura: a pressão de saturação, aquela em que líquido e vapor têm o mesmo potencial químico. Ela é determinada pela construção de áreas iguais de Maxwell: traça-se a horizontal que corta o laço de modo que as duas áreas fechadas acima e abaixo dela sejam iguais.
Fora dessa pressão específica, as três raízes ainda existem, mas o que elas descrevem é outra coisa: uma das fases é estável e a outra é metaestável, isto é, líquido superaquecido ou vapor super-resfriado, estados que existem de fato (é o que permite superaquecer água no micro-ondas), mas que colapsam ao menor distúrbio. A raiz do meio, essa sim, não descreve estado nenhum em condição alguma.
Para este artigo isso não muda nada: as três raízes são as mesmas, e o ponto (que a matemática entrega opções e a física escolhe) continua de pé. Mas vale saber que o rótulo “coexistência” esconde uma construção geométrica a mais.
numpy.roots não sabe nada de termodinâmica: ele só devolve as 3 raízes
matemáticas da cúbica. A física, a estabilidade termodinâmica, é quem
decide qual delas é o volume “de verdade” em cada fase. É o mesmo tipo de
cuidado já visto no artigo sobre
polinômios com SymPy
(Poly.all_roots() x Poly.real_roots()), só que agora com um
significado físico concreto por trás de cada raiz. O gráfico a seguir
marca as três raízes sobre a curva, anotando o papel físico de cada uma.
import matplotlib.pyplot as plt
import numpy as np
R = 0.082057366080960 # atm L / (mol K)
a_co2, b_co2 = 3.640, 0.04267 # L^2 atm / mol^2 ; L / mol
T_c = 8 * a_co2 / (27 * R * b_co2)
P_c = a_co2 / (27 * b_co2**2)
T_subcritica = 0.84 * T_c
P_subcritica = 0.41 * P_c
def van_der_waals_coeffs(P, T, n=1.0):
c3 = P
c2 = -(P * n * b_co2 + n * R * T)
c1 = a_co2 * n**2
c0 = -a_co2 * n**3 * b_co2
return [c3, c2, c1, c0]
raizes = np.roots(van_der_waals_coeffs(P_subcritica, T_subcritica))
raizes_reais_positivas = np.sort(raizes[(np.abs(raizes.imag) < 1e-8) & (raizes.real > 0)].real)
V_range_sub = np.linspace(0.05, 1.5, 800)
def van_der_waals_sub(V):
return (P_subcritica + a_co2 / V**2) * (V - b_co2) - R * T_subcritica
V_liquido, V_instavel, V_vapor = raizes_reais_positivas
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(V_range_sub, van_der_waals_sub(V_range_sub), color="C0")
ax.axhline(0, color="gray", linewidth=1, zorder=0)
ax.scatter([V_liquido, V_instavel, V_vapor], [0, 0, 0], color="C2", zorder=3,
label="3 raízes reais")
ax.annotate("líquido", (V_liquido, 0), textcoords="offset points", xytext=(-55, 45),
ha="left", arrowprops=dict(arrowstyle="-", color="gray", linewidth=1))
ax.annotate("instável", (V_instavel, 0), textcoords="offset points", xytext=(40, 55),
ha="center", arrowprops=dict(arrowstyle="-", color="gray", linewidth=1))
ax.annotate("vapor", (V_vapor, 0), textcoords="offset points", xytext=(0, 12), ha="center")
ax.legend()
ax.set_xlabel("V (L/mol)")
ax.set_ylabel("resíduo de van der Waals")
ax.set_title(f"Van der Waals abaixo de T_c: 3 raízes reais (T={T_subcritica:.0f} K)")
As Isotermas, e o Que “Coexistência” Realmente Exige #
Todas as figuras desta seção mostraram o resíduo \(f(V)\), perfeito para enxergar “isto é um problema de raiz”, mas mudo sobre a física. A figura canônica do assunto é o diagrama \(P\)-\(V\) com isotermas, e vale desenhá-la, porque nela a ressalva da callout anterior deixa de ser abstrata.
Achar a pressão de saturação é, ela mesma, um problema de raiz, só que aninhado. A condição de Maxwell é que a horizontal divida o laço em áreas iguais:
$$ \int_{V_{liq}}^{V_{vap}} \left[P_{vdW}(V) - P_{sat}\right] dV = 0 $$E ela tem primitiva elementar, \(\int P dV = RT\ln(V-b) + a/V\), o que evita quadratura numérica perto da singularidade em \(V = b\):
import numpy as np
from scipy import optimize
R = 0.082057366080960 # atm L / (mol K)
a_co2, b_co2 = 3.640, 0.04267 # L^2 atm / mol^2 ; L / mol
T_c = 8 * a_co2 / (27 * R * b_co2)
P_c = a_co2 / (27 * b_co2**2)
T_subcritica = 0.84 * T_c
P_subcritica = 0.41 * P_c
def van_der_waals_coeffs(P, T, n=1.0):
c3 = P
c2 = -(P * n * b_co2 + n * R * T)
c1 = a_co2 * n**2
c0 = -a_co2 * n**3 * b_co2
return [c3, c2, c1, c0]
def P_vdw(V, T):
return R * T / (V - b_co2) - a_co2 / V**2
def dPdV_vdw(V, T):
return -R * T / (V - b_co2) ** 2 + 2 * a_co2 / V**3
def espinodais(T):
"""Os dois V onde dP/dV = 0: o fundo do vale e o topo da corcova do laço."""
return (optimize.brentq(dPdV_vdw, b_co2 * 1.001, 3 * b_co2, args=(T,)),
optimize.brentq(dPdV_vdw, 3 * b_co2, 50, args=(T,)))
def pressao_saturacao(T):
V_min, V_max = espinodais(T) # onde dP/dV = 0
P_lo, P_hi = P_vdw(V_min, T), P_vdw(V_max, T) # mínimo e máximo locais do laço
def primitiva(V):
return R * T * np.log(V - b_co2) + a_co2 / V
def desbalanco(P):
r = np.roots(van_der_waals_coeffs(P, T))
r = np.sort(r[np.abs(r.imag) < 1e-9].real)
return (primitiva(r[-1]) - primitiva(r[0])) - P * (r[-1] - r[0])
return optimize.brentq(desbalanco, max(P_lo, 1e-9) * 1.0001, P_hi * 0.9999)
P_sat = pressao_saturacao(T_subcritica)
print(f"P_sat = {P_sat:.2f} atm -> Pr = {P_sat / P_c:.3f}")P_sat = 35.44 atm -> Pr = 0.479O exemplo usou \(P_r = 0{,}41\), mas a pressão de equilíbrio naquela temperatura é \(P_r = 0{,}479\). As três raízes existem nas duas pressões, mas só numa delas líquido e vapor coexistem de verdade. Na outra, uma das fases é metaestável.
A curva de saturação de van der Waals em variáveis reduzidas é tabelada, e não depende da substância. Serve como teste do código:
| \(T_r\) | \(P_r\) calculado | \(P_r\) tabelado |
|---|---|---|
| 0,95 | 0,8119 | 0,812 |
| 0,90 | 0,6470 | 0,647 |
| 0,85 | 0,5045 | 0,504 |
| 0,80 | 0,3834 | 0,383 |
Bate nas quatro. Vale o aviso: um brentq que recebe o bracket com os
extremos invertidos não reclama e devolve um número plausível (aqui,
\(P_r = 0{,}725\) para \(T_r = 0{,}90\)). Só a comparação com a tabela
denuncia o erro.
A figura a seguir traz, à esquerda, as isotermas clássicas (acima, em e abaixo de \(T_c\)) com o ponto crítico marcado; e à direita, um zoom no laço subcrítico com a construção de Maxwell desenhada:
No painel da direita, as duas regiões sombreadas têm formatos bem diferentes (a da esquerda é alta e estreita, a da direita é baixa e larga), mas a mesma área. É essa igualdade que define \(P_{sat}\): a horizontal de equilíbrio é a única altura em que as duas se compensam.
Repare também onde o fundo do laço vai parar: abaixo de zero, até cerca de \(-2\) atm. Van der Waals prevê pressão negativa nessa faixa, o que não descreve nenhum estado de equilíbrio real; é mais um sinal de que o trecho central do laço não é física: é o preço de forçar uma única equação suave a atravessar uma região onde, de fato, coexistem duas fases. A construção de Maxwell existe justamente para substituir esse trecho pela horizontal.
Repare no painel da esquerda: acima de \(T_c\) a isoterma desce sem nenhuma ondulação. Não há laço, não há três raízes, não há transição de fase. É exatamente isso que van der Waals previu e que a lei dos gases ideais jamais poderia prever: existe uma temperatura acima da qual líquido e vapor deixam de ser coisas distintas.
3. Fator de Atrito de Colebrook: Da Régua de Cálculo ao brentq
#
Última equação, e a mais nova das três, mas a que mais aparece no dia a dia de um projeto de engenharia. Publicada por Cyril Colebrook em 1939, a partir dos experimentos que ele fez com Cedric White em 1937 (por isso também se chama equação de Colebrook-White), ela dá o fator de atrito em escoamento turbulento dentro de tubulações de forma implícita (veja o verbete Darcy friction factor formulae):
$$ \frac{1}{\sqrt{f}} = -2\log_{10}\left(\frac{\varepsilon/D}{3{,}7} + \frac{2{,}51}{Re\sqrt{f}}\right) $$O \(f\) aparece dos dois lados, e dentro de um logaritmo: não há como
isolá-lo. Cinco anos depois, em 1944, Lewis Moody publicou o ábaco que
levaria seu nome justamente para poupar os engenheiros dessa conta: com
régua e olho, lia-se o \(f\) num gráfico em vez de iterar à mão. Foi assim,
por décadas, até as calculadoras programáveis. Hoje, é uma
chamada de brentq, e, mais adiante nesta seção, reconstruímos o próprio
ábaco de Moody a partir dela.
Exemplo 3: água num tubo de aço comercial
Enunciado: Água escoa a \(2\) m/s (velocidade de projeto típica) num tubo de aço comercial de \(100\) mm de diâmetro. A rugosidade absoluta do aço comercial é \(\varepsilon \approx 0{,}045\) mm (valor tabelado de engenharia); a viscosidade cinemática da água a 20 °C é \(\nu \approx 10^{-6}\) m²/s. Qual é o fator de atrito?
Solução:
import numpy as np
from scipy import optimize
D_tubo = 0.10 # m — tubo de 100 mm
eps_aco = 0.045e-3 # m — rugosidade absoluta do aço comercial (tabelada)
eps_over_D = eps_aco / D_tubo
nu_agua = 1.0e-6 # m²/s — viscosidade cinemática da água a ~20°C
v_escoamento = 2.0 # m/s — velocidade de projeto típica em tubulações
Re = v_escoamento * D_tubo / nu_agua
print(f"eps/D = {eps_over_D}")
print(f"Re = {Re:.0f}")
def colebrook(f, eps_D, Re_):
return 1 / np.sqrt(f) + 2 * np.log10(eps_D / 3.7 + 2.51 / (Re_ * np.sqrt(f)))
f_friction = optimize.brentq(colebrook, 1e-4, 0.5, args=(eps_over_D, Re))
print("fator de atrito (brentq):", f_friction)eps/D = 0.00045
Re = 200000
fator de atrito (brentq): 0.01856015225418919Como nos exemplos anteriores, vale ver a curva da equação de Colebrook (reorganizada como resíduo) cruzando zero no fator de atrito encontrado.
import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize
D_tubo = 0.10 # m — tubo de 100 mm
eps_aco = 0.045e-3 # m — rugosidade absoluta do aço comercial (tabelada)
eps_over_D = eps_aco / D_tubo
nu_agua = 1.0e-6 # m²/s — viscosidade cinemática da água a ~20°C
v_escoamento = 2.0 # m/s — velocidade de projeto típica em tubulações
Re = v_escoamento * D_tubo / nu_agua
def colebrook(f, eps_D, Re_):
return 1 / np.sqrt(f) + 2 * np.log10(eps_D / 3.7 + 2.51 / (Re_ * np.sqrt(f)))
f_friction = optimize.brentq(colebrook, 1e-4, 0.5, args=(eps_over_D, Re))
f_range = np.linspace(0.008, 0.1, 400)
fig, ax = plt.subplots(figsize=(8, 5))
ax.plot(f_range, colebrook(f_range, eps_over_D, Re), color="C0")
ax.axhline(0, color="gray", linewidth=1, zorder=0)
ax.scatter([f_friction], [0], color="C2", zorder=3, label=f"f ≈ {f_friction:.4f}")
ax.legend()
ax.set_xlabel("f")
ax.set_ylabel("resíduo de Colebrook")
ax.set_title("Fator de atrito: Colebrook como problema de raiz")
A Iteração Manual, Por Trás da Frase “Antes de Calculadoras” #
A equação de Colebrook é implícita (\(f\) aparece dos dois lados), mas pode ser reorganizada isolando o \(f\) do lado esquerdo:
$$ f = \left[-2\log_{10}\left(\frac{\varepsilon/D}{3{,}7} + \frac{2{,}51}{Re\sqrt{f}}\right)\right]^{-2} $$Chutando um valor inicial de \(f\) e substituindo repetidamente do lado direito para o esquerdo, a iteração de ponto fixo, dá para resolver a equação sem nenhum método sofisticado. É exatamente o que engenheiros faziam à mão, ou com tabelas, antes de calculadoras e computadores.
import numpy as np
from scipy import optimize
D_tubo = 0.10 # m — tubo de 100 mm
eps_aco = 0.045e-3 # m — rugosidade absoluta do aço comercial (tabelada)
eps_over_D = eps_aco / D_tubo
nu_agua = 1.0e-6 # m²/s — viscosidade cinemática da água a ~20°C
v_escoamento = 2.0 # m/s — velocidade de projeto típica em tubulações
Re = v_escoamento * D_tubo / nu_agua
def colebrook(f, eps_D, Re_):
return 1 / np.sqrt(f) + 2 * np.log10(eps_D / 3.7 + 2.51 / (Re_ * np.sqrt(f)))
f_friction = optimize.brentq(colebrook, 1e-4, 0.5, args=(eps_over_D, Re))
def colebrook_ponto_fixo(f, eps_D, Re_):
return (-2 * np.log10(eps_D / 3.7 + 2.51 / (Re_ * np.sqrt(f)))) ** -2
f_atual = 0.02 # chute inicial razoável (ordem de grandeza típica)
historico_f = [f_atual]
for _ in range(8):
f_atual = colebrook_ponto_fixo(f_atual, eps_over_D, Re)
historico_f.append(f_atual)
for i, f_i in enumerate(historico_f):
print(f"iteração {i}: f = {f_i:.8f} erro vs. brentq = {abs(f_i - f_friction):.2e}")iteração 0: f = 0.02000000 erro vs. brentq = 1.44e-03
iteração 1: f = 0.01849038 erro vs. brentq = 6.98e-05
iteração 2: f = 0.01856372 erro vs. brentq = 3.57e-06
iteração 3: f = 0.01855997 erro vs. brentq = 1.82e-07
iteração 4: f = 0.01856016 erro vs. brentq = 9.28e-09
iteração 5: f = 0.01856015 erro vs. brentq = 4.73e-10
iteração 6: f = 0.01856015 erro vs. brentq = 2.41e-11
iteração 7: f = 0.01856015 erro vs. brentq = 1.23e-12
iteração 8: f = 0.01856015 erro vs. brentq = 6.28e-14O gráfico a seguir mostra essa convergência em escala logarítmica: o erro cai de forma aproximadamente linear a cada iteração.
import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize
D_tubo = 0.10 # m — tubo de 100 mm
eps_aco = 0.045e-3 # m — rugosidade absoluta do aço comercial (tabelada)
eps_over_D = eps_aco / D_tubo
nu_agua = 1.0e-6 # m²/s — viscosidade cinemática da água a ~20°C
v_escoamento = 2.0 # m/s — velocidade de projeto típica em tubulações
Re = v_escoamento * D_tubo / nu_agua
def colebrook(f, eps_D, Re_):
return 1 / np.sqrt(f) + 2 * np.log10(eps_D / 3.7 + 2.51 / (Re_ * np.sqrt(f)))
f_friction = optimize.brentq(colebrook, 1e-4, 0.5, args=(eps_over_D, Re))
def colebrook_ponto_fixo(f, eps_D, Re_):
return (-2 * np.log10(eps_D / 3.7 + 2.51 / (Re_ * np.sqrt(f)))) ** -2
historico_f = [0.02]
for _ in range(8):
historico_f.append(colebrook_ponto_fixo(historico_f[-1], eps_over_D, Re))
fig, ax = plt.subplots(figsize=(8, 5))
erro_ponto_fixo = [abs(f_i - f_friction) for f_i in historico_f]
ax.semilogy(range(len(erro_ponto_fixo)), erro_ponto_fixo, "o-", color="C0")
ax.set_xlabel("iteração")
ax.set_ylabel("|f − f_brentq|")
ax.set_title("Iteração de ponto fixo do Colebrook convergindo")
Convergiu para o mesmo valor do brentq em poucas iterações, sem
precisar de bracket nem de derivada, só repetir a substituição. É
menos robusto que brentq em geral (fórmulas de ponto fixo podem
divergir dependendo de como a equação é reorganizada), mas mostra que a
“mágica” por trás de uma chamada de brentq é, no fundo, a mesma ideia
manual que já existia décadas antes.
Swamee-Jain: a Aproximação Explícita Que a Indústria Usa #
Há um terceiro caminho, entre a iteração manual e o brentq, e ele completa
o arco. Em 1976, Prabhata Swamee e Akalank Jain publicaram uma fórmula
explícita que aproxima Colebrook sem iteração nenhuma, em uma linha de
conta:
import numpy as np
from scipy import optimize
def colebrook(f, eps_D, Re_):
return 1 / np.sqrt(f) + 2 * np.log10(eps_D / 3.7 + 2.51 / (Re_ * np.sqrt(f)))
def swamee_jain(eps_D, Re_):
return 0.25 / np.log10(eps_D / 3.7 + 5.74 / Re_**0.9) ** 2
print(f"{'Re':>8} {'eps/D':>9} {'brentq':>9} {'Swamee-Jain':>12} {'erro':>7}")
for Re_i in (1e4, 2e5, 1e6, 1e7):
for eps_i in (1e-7, 0.00045, 0.01):
f_exato = optimize.brentq(colebrook, 1e-4, 0.5, args=(eps_i, Re_i))
f_aprox = swamee_jain(eps_i, Re_i)
print(f"{Re_i:>8.0e} {eps_i:>9.5f} {f_exato:>9.5f} {f_aprox:>12.5f}"
f" {100*(f_aprox-f_exato)/f_exato:>+6.2f}%") Re eps/D brentq Swamee-Jain erro
1e+04 0.00000 0.03088 0.03097 +0.29%
1e+04 0.00045 0.03157 0.03175 +0.58%
1e+04 0.01000 0.04313 0.04404 +2.12%
2e+05 0.00000 0.01564 0.01553 -0.68%
2e+05 0.00045 0.01856 0.01867 +0.57%
2e+05 0.01000 0.03821 0.03836 +0.40%
1e+06 0.00000 0.01165 0.01161 -0.33%
1e+06 0.00045 0.01686 0.01695 +0.56%
1e+06 0.01000 0.03796 0.03801 +0.12%
1e+07 0.00000 0.00811 0.00815 +0.50%
1e+07 0.00045 0.01637 0.01640 +0.16%
1e+07 0.01000 0.03791 0.03792 +0.02%Erro abaixo de \(1\%\) em quase toda a faixa útil, chegando a \(2\%\) só no canto mais extremo (tubo muito rugoso em Reynolds baixo). Para dimensionar uma tubulação (onde a rugosidade tabelada já tem incerteza maior que isso), é mais que suficiente, e é por isso que Swamee-Jain aparece em planilhas de projeto até hoje.
Fecha-se assim um arco de três etapas que se repete em toda a engenharia:
| Abordagem | Custo | Quando usar |
|---|---|---|
| Implícita (Colebrook, 1939) | exige iteração | é a referência: define o valor “certo” |
| Iterativa manual (ponto fixo) | ~5 substituições à mão | quando não há computador |
| Explícita aproximada (Swamee-Jain, 1976) | uma linha, sem iteração | quando \(1\%\) basta, quase sempre, na prática |
O brentq não tornou as duas últimas obsoletas: ele tornou a primeira
barata o suficiente para ser usada direto. Mas saber que as outras existem
explica por que tanta planilha de engenharia ainda não chama solver nenhum.
Diagrama de Moody #
A equação de Colebrook só vale para escoamento turbulento
(\(Re \gtrsim 4000\)); no regime laminar o fator de atrito tem fórmula
fechada, \(f = 64/Re\). Resolvendo Colebrook via brentq para uma faixa
de Reynolds e várias rugosidades relativas, reconstruímos o
diagrama de Moody, uma
representação gráfica clássica que relaciona o fator de atrito, o número
de Reynolds e a rugosidade relativa, uma das figuras clássicas da
engenharia de fluidos.
import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize
D_tubo = 0.10 # m — tubo de 100 mm
eps_aco = 0.045e-3 # m — rugosidade absoluta do aço comercial (tabelada)
eps_over_D = eps_aco / D_tubo
nu_agua = 1.0e-6 # m²/s — viscosidade cinemática da água a ~20°C
v_escoamento = 2.0 # m/s — velocidade de projeto típica em tubulações
Re = v_escoamento * D_tubo / nu_agua
def colebrook(f, eps_D, Re_):
return 1 / np.sqrt(f) + 2 * np.log10(eps_D / 3.7 + 2.51 / (Re_ * np.sqrt(f)))
f_friction = optimize.brentq(colebrook, 1e-4, 0.5, args=(eps_over_D, Re))
Re_range = np.logspace(np.log10(4000), 8, 200) # regime turbulento
rugosidades = [0.0, 0.00001, eps_over_D, 0.01] # inclui o tubo de aço do exemplo
# nota sobre o max(..., 1e-7) que aparece abaixo: "tubo liso" é eps/D = 0, mas
# aí o argumento do log10 vira só o termo de Reynolds — o que é legítimo, e é
# de fato o limite liso. O piso de 1e-7 existe só para o caso eps/D = 0 não
# produzir um zero exato dentro do log em nenhuma variante da fórmula; é um
# valor pequeno o bastante para ser indistinguível de zero no resultado
fig, ax = plt.subplots(figsize=(9, 6))
for eps_D_i in rugosidades:
f_vals = [optimize.brentq(colebrook, 1e-4, 0.5, args=(max(eps_D_i, 1e-7), Re_i))
for Re_i in Re_range]
if eps_D_i == 0.0:
label = "tubo liso (ε/D → 0)"
elif eps_D_i == eps_over_D:
label = f"ε/D = {eps_D_i:g} (tubo do exemplo)"
else:
label = f"ε/D = {eps_D_i:g}"
ax.loglog(Re_range, f_vals, label=label)
ax.scatter([Re], [f_friction], color="black", marker="*", s=220, zorder=4,
edgecolors="white", linewidths=1.2, label="exemplo (Re, f)")
ax.set_xlabel("número de Reynolds (Re)")
ax.set_ylabel("fator de atrito (f)")
ax.set_title("Diagrama de Moody (regime turbulento, via Colebrook + brentq)")
ax.legend(fontsize=10)
O tubo de aço do exemplo (curva destacada) está claramente acima do tubo
liso: a rugosidade real do material aumenta o atrito de forma perceptível
mesmo em regime turbulento estabelecido. É esse tipo de leitura (como o
fator de atrito responde a Reynolds e à rugosidade, simultaneamente) que
o diagrama de Moody entregava antes de qualquer calculadora, e que
brentq reconstrói aqui em 800 chamadas de função (4 curvas × 200 pontos
de Reynolds cada).
Tabela-Resumo #
| Problema | Equação | Método | Resultado |
|---|---|---|---|
| pH de ácido fraco | \(x^2 + K_a x - K_a C_0 = 0\) | numpy.roots |
Exato: pH \(\approx 2{,}875\); aproximado: pH \(\approx 2{,}872\) (diferença \(0{,}003\)) |
| Onde a aproximação falha | Mesma equação, \(C_0\) variável | numpy.roots repetido |
Erro passa de \(0{,}05\) unidade de pH abaixo de \(C_0 \approx 3\times10^{-4}\) mol/L |
| Volume molar do CO2 real | \((P+an^2/V^2)(V-nb)=nRT\) | brentq |
\(V \approx 2{,}35\) L, \(4{,}4\%\) menor que o volume ideal |
| Coexistência líquido-vapor | Mesma cúbica, \(T_r=0{,}84\), \(P_r=0{,}41\) | numpy.roots |
\(T_c \approx 308{,}0\) K, \(P_c \approx 74{,}0\) atm; 3 raízes reais: líquido \(\approx 0{,}071\), instável \(\approx 0{,}135\), vapor \(\approx 0{,}536\) L/mol |
| Fator de atrito (Colebrook) | \(1/\sqrt{f} = -2\log_{10}(\ldots)\) | brentq |
\(f \approx 0{,}0186\) para tubo de aço, \(Re=200{.}000\) |
| Iteração de ponto fixo | Mesma equação, isolando \(f\) | Substituição repetida | Converge ao valor do brentq em \(\sim\)5 iterações |
| Diagrama de Moody | Colebrook para várias \(Re\) e \(\varepsilon/D\) | brentq repetido (800x) |
Reconstrói o diagrama clássico; tubo do exemplo bem acima da curva de tubo liso |
Próximos Passos: Fim do Arco de Equações Numéricas #
Este artigo fecha o arco de três partes sobre resolver equações numericamente com SciPy e NumPy. O primeiro artigo apresentou os métodos (bisseção, Brent, Newton, sistemas não lineares) em exemplos didáticos; o segundo aplicou esses métodos à física, de órbitas planetárias a projéteis com resistência do ar; este terceiro levou os mesmos métodos à química e à engenharia, sempre partindo de um cenário numérico real, nunca de um valor arbitrário.
O que os três artigos têm em comum é o tipo de problema: sempre uma
equação escalar, ou um sistema pequeno, resolvido ponto a ponto. Falta uma
peça inteira: e quando o objeto de interesse não é a raiz de uma equação
qualquer, mas um polinômio, com sua própria álgebra de derivadas,
integrais, produtos e ajuste a dados experimentais? No próximo artigo, Polinômios com NumPy,
vemos como o NumPy trata polinômios como um objeto de primeira classe, do
mesmo numpy.roots já usado neste artigo à tarefa que o SymPy não faz
bem: ajustar uma curva a dados reais com ruído.
