Da Química à Engenharia: Equações Numéricas em Sistemas Reais

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 — não 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\),

— e cúbicas, na prática, não têm solução fechada útil. Essa não é a única equação 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 não é “como resolver” e sim “quando vale a pena”. 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,

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.

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.00291

Para 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.

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.7927

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.

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 CO₂ 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?

#

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:

Multiplicando por \(V^2\) e agrupando por potência de \(V\):

Agora é um polinômio de grau 3 com coeficientes explícitos — exatamente o formato que o numpy.roots espera.

Isso acontece 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:

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.

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:

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.479

E aí está o ponto: o 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 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, e sim 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):

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 cerca de quatro 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.

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.01856015225418919

Como 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:

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-14

O 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, 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 e normas de projeto até hoje.

Fecha-se assim um arco de três etapas que se repete em toda a engenharia:

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 mais usadas (e, por décadas, mais penosas de reproduzir manualmente) 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

#

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.