Como a Matemática Pode Salvar 20 Milhões de Toneladas de Grãos por Ano

Dentro do modelo que usa física de percolação, CFD e cinética de Arrhenius para prever a deterioração de grãos armazenados — semanas antes que qualquer termômetro perceba.

Por Giselle Falcão · 18 min de leitura

Existe um problema que custa ao Brasil mais de R$ 7 bilhões por ano e que quase ninguém fora do agronegócio conhece: a deterioração de grãos dentro de silos.

Não estou falando de pragas óbvias ou de ratos roendo sacas. Estou falando de um processo invisível, silencioso, governado por equações diferenciais parciais acopladas, que transforma milhões de toneladas de milho, soja e trigo em massa contaminada por micotoxinas — tudo por causa de um gradiente de temperatura de 5°C entre o centro e a parede do silo.

Nos últimos meses, mergulhei nesse problema. O resultado é um modelo matemático-computacional que integra três componentes que nunca haviam sido acoplados antes: convecção natural em meio poroso, teoria da percolação topológica e liberação distribuída de CO₂.

Neste artigo, vou destrinchar tudo — da física ao código, das equações às limitações — para que você entenda não apenas o que o modelo faz, mas por que funciona e onde ele falha.

O problema: o que acontece dentro de um silo

Um silo de grãos não é um depósito estático. É um reator bioquímico passivo. Os grãos respiram. Os fungos respiram. E essa respiração produz calor, CO₂ e água metabólica — criando um ciclo de retroalimentação positiva que, se não detectado a tempo, é irreversível.

Mas o gatilho de tudo isso não é biológico. É termodinâmico.

O ciclo funciona assim:

Gradiente Térmico (ΔT)
        ↓
Convecção Natural do Ar (ar quente sobe, frio desce)
        ↓
Migração de Umidade (vapor viaja com o ar)
        ↓
Condensação + Hotspot (vapor condensa nas zonas frias)
        ↓
Respiração → CO₂ + Calor + Água metabólica
        ↓
Fungos + Micotoxinas
        ↓
  ↩ RETROALIMENTAÇÃO → mais calor → mais gradiente → mais migração

A parede metálica do silo esquenta durante o dia e esfria à noite. Essa diferença de temperatura gera correntes de convecção natural dentro da massa de grãos — o ar quente sobe pelo centro e desce pelas paredes. Esse ar carrega vapor d’água, que condensa nas regiões frias (geralmente o topo e as paredes), criando hotspots de umidade. E é nesses hotspots que os fungos prosperam.

A termometria convencional — o termômetro pendurado dentro do silo — só detecta o problema quando a temperatura já subiu vários graus. Nesse ponto, o dano já está feito. O CO₂ é um indicador 3 a 5 semanas mais rápido. — White et al., Canadian Agricultural Engineering, 1982

A solução: três modelos em um

A maioria dos modelos existentes na literatura trata a migração de umidade como um problema de difusão pura — a equação de Richards, basicamente. Outros, mais sofisticados, adicionam convecção natural via CFD. Mas nenhum, até onde pesquisamos, faz as três coisas ao mesmo tempo.

O modelo que desenvolvemos acopla:

✅ Os 3 componentes do modelo acoplado

1. Convecção natural — via Lei de Darcy + Boussinesq, capturando as correntes de ar que transportam vapor
2. Percolação topológica — via teoria de percolação de sítios, modelando como grãos úmidos formam caminhos conectados
3. CO₂ distribuído — tratando a concentração como campo espacial (não como um número médio no topo do silo)

Por que a convecção importa?

A análise do número de Péclet responde a essa pergunta de forma definitiva. Para velocidades de convecção de ~0,05 m/s, comprimento do silo L = 10 m e difusividade D ≈ 7×10⁻⁸ m²/s:

Pe = u·L / Deff ≈ (0,05 × 10) / (7 × 10⁻⁸) ≈ 7 × 10⁶ >> 1

Um Péclet de 7 milhões. Isso significa que ignorar a convecção — como fazem os modelos de Richards — é como tentar modelar um rio usando apenas a lei de Fick. A física simplesmente não fecha.

Por que a percolação importa?

Aqui está a parte que eu considero mais elegante do modelo. A teoria da percolação, originária da física estatística, descreve como caminhos conectados se formam em redes aleatórias.

Imagine cada grão no silo como um nó em uma rede. Se a umidade desse grão está acima do limiar crítico (14% para milho), ele é “ativo”. Quando a fração de grãos ativos ultrapassa um limiar — chamado de limiar de percolação pꟲ — forma-se um caminho contínuo de alta umidade que atravessa todo o silo.

┌─────────────────┐   ┌─────────────────┐   ┌─────────────────┐
│   p = 0,40      │   │  p = 0,59 (pꟲ)  │   │   p = 0,70      │
│                 │   │                 │   │                 │
│ Clusters        │──▶│ TRANSIÇÃO DE    │──▶│ Migração        │
│ isolados        │   │ FASE: cluster   │   │ pervasiva       │
│                 │   │ percolante      │   │                 │
│ 🟢 Reversível   │   │ 🟡Ponto crítico│   │ 🔴 Irreversível  │
└─────────────────┘   └─────────────────┘   └─────────────────┘

A contribuição central do nosso modelo é traduzir essa transição de fase em uma equação que modifica a difusividade efetiva:

Deff(p) = D₀ [1 + αₚ · P∞(p)]

Eq. 10 — Acoplamento percolação–difusão: a difusividade salta quando o cluster percolante se forma.

Na prática, isso significa que a partir do dia ~12 (quando a fração de grãos úmidos atinge pꟲ), a difusividade do meio sofre um salto não-linear. Sem este acoplamento, o modelo subestima o CO₂ em ~900 ppm aos 30 dias.

As equações que governam o silo

Para quem quer ir fundo, aqui estão as 5 equações de conservação que formam o coração do modelo. Todas em coordenadas cilíndricas axissimétricas (r, z):

Eq. 1 — Conservação de massa do ar

∂(ερₐ)/∂t + ∇·(ρₐu) = 0

Eq. 2 — Momentum (Darcy)

u = –(K/μ)(∇p – ρₐgẑ)

Permeabilidade estimada por Ergun: K = dₚ²ε³/[150(1–ε)²] ≈ 1,8 × 10⁻⁹ m².

Eq. 3 — Conservação de energia

(ρCₚ)eff ∂T/∂t + (ρCₚ)ₐ u·∇T = ∇·(keff∇T) + Qresp

O termo advectivo (ρCₚ)ₐ u·∇T é o que diferencia o modelo de formulações puramente difusivas.

Eq. 4 — Conservação de umidade

ρb ∂M/∂t + ∇·(ρᵥu) = ∇·(Deff∇M) + DT∇²T

Aqui entra o acoplamento com a Eq. 10: Deff não é constante, ele depende do estado da rede de percolação.

Eq. 5 — Conservação de CO₂ (distribuído)

∂C_CO₂/∂t + u·∇C_CO₂ = ∇·(D_CO₂∇C_CO₂) + R_CO₂

Com D_CO₂ = Dmol/τ (tortuosidade τ ≈ 1,5). O CO₂ é tratado como campo distribuído, não como valor médio no plenum.

Eq. 6–7 — Cinética de Arrhenius

R_CO₂ = A · exp(–Eₐ/RT) · f(w)
f(w) = 1 / [1 + exp(–κ(w – wc))]

Eₐ ≈ 50 kJ/mol, wc ≈ 14% base úmida para milho.

Eq. 8 — Isoterma de Henderson-Thompson

Me = [–ln(1 – aw) / (k₁(T + k₂))]^(1/k₃)

O que descobrimos: três modelos em confronto direto

Para demonstrar que cada componente do modelo agrega valor real — e não apenas complexidade — comparamos três modelos sob condições idênticas:

Modelo RMSE Umidade RMSE Temp. RMSE CO₂ Erro Hotspot A (difusivo puro) 1,82 %b.u. 2,65 °C 234 ppm 48% B (difusão + convecção) 0,62 %b.u. 1,68 °C 112 ppm 22% C (proposto) 0,29 %b.u. 1,37 °C 68 ppm 8% Redução C vs A −84% −48% −71% −83%

Os números são claros: o modelo puramente difusivo erra a localização de hotspots em quase metade dos casos. O modelo proposto reduz esse erro para 8%.

O que a análise de Sobol revelou

Usando 10.000 amostras quasi-aleatórias, identificamos o ranking de importância dos parâmetros:

Parâmetro Sᴛ (índice total) Papel 1 Umidade do grão (w) 0,45 Dominante 2 Temperatura (T) 0,31 Forte 3 Difusividade efetiva (Deff) 0,18 Moderado 4 Ventilação (kᵥ) 0,14 Moderado (alta incerteza) 5 Energia de ativação (Eₐ) 0,12 Moderado 6 Acoplamento percolativo (αₚ) 0,11 Opera via interações

A soma ΣSᴛ = 1,55 > 1 confirma interações fortes entre parâmetros — esperado em modelos acoplados.

Validação cruzada temporal (5-fold)

Para mitigar o risco de overfitting:

• R² médio de calibração: 0,910
• R² médio de validação: 0,884
• Gap: < 3% → evidência consistente de ausência de overfitting significativo

Onde o modelo falha — e eu sei disso

Nenhum modelo é perfeito, e esconder as limitações seria intelectualmente desonesto. Aqui estão as principais fragilidades que reconheço explicitamente:

🔴 Limitações conhecidas

• Percolação 2D em rede regular: o silo é 3D, desordenado, anisotrópico. Usamos rede quadrada 2D (pꟲ = 0,593), mas o limiar real em 3D seria pꟲ ≈ 0,312. Isso significa que o limiar que reportamos é uma cota superior — na prática, a transição acontece mais cedo.
• Acoplamento sequencial (weak coupling): a percolação atualiza Deff a cada 24h, não em tempo real. Um acoplamento implícito (forte) exigiria reformulação do tensor de difusividade.
• Validação em silo único: 500 t de milho, Triângulo Mineiro, uma safra. Generalização para soja, trigo, arroz e outros climas requer novas campanhas.
• CO₂ medido manualmente a cada 48h: baixa resolução temporal. Sensores NDIR permanentes melhorariam significativamente a calibração.
• R² alto pode mascarar overfitting: mitigamos com validação cruzada 5-fold (gap < 3%), mas o risco existe com múltiplos parâmetros calibráveis (kᵥ, αₚ, τ).
• Aproximação axissimétrica: assume ∂/∂θ = 0, válido para silos cilíndricos com aquecimento uniforme. Silos com insolação lateral assimétrica ou geometria não-cilíndrica exigiriam 3D completo.

🟠 Direções de melhoria

Se eu recomeçasse amanhã, faria:

• Percolação contínua (continuum percolation) ou redes de Voronoi para representar o empacotamento real dos grãos
• Simulação 3D completa para silos com insolação assimétrica
• Acoplamento implícito entre percolação e difusividade, atualizando a cada passo de tempo
• Modelos híbridos CFD + ML para assimilação contínua de dados de sensores
• Modelos de crescimento fúngico acoplados ao transporte de umidade, para prever produção de micotoxinas
• Propagação formal de incertezas via Monte Carlo com amostragem estratificada (Latin Hypercube)

O pipeline computacional — na prática

┌──────────┐    ┌──────────┐    ┌──────────┐    ┌──────────┐
│ ENTRADA  │───▶│CFD SOLVER│───▶│PERCOLAÇÃO│───▶│  SAÍDA   │
│          │    │          │    │          │    │          │
│ T(r,z)   │    │ Eqs.1–5  │    │ Hoshen-  │    │ Índice   │
│ M, UR    │    │ SIMPLE   │    │ Kopelman │    │ de risco │
│ CO₂      │    │ FVM      │    │ p→P∞(p)  │    │ Previsão │
│ geometria│    │ 120×80   │    │ →Deff    │    │ 48h      │
└──────────┘    └──────────┘    └────┬─────┘    └──────────┘
                     ▲               │
                     └───────────────┘
                   atualiza a cada 24h simuladas

     CFD completo: 18–22 min        Modelo reduzido: < 10 s
    (Intel Xeon E5-2630)           (qualquer desktop)

O modelo completo (CFD com malha 120×80) executa em 18–22 minutos por simulação de 30 dias. Para uso operacional, o modelo reduzido de parâmetros concentrados preserva R² > 0,90 e executa em menos de 10 segundos.

A justificativa 2D — por que funciona (e onde não funciona)

Uma pergunta que qualquer parecerista faria: por que 2D se o silo é 3D?

Há uma distinção importante:

O domínio CFD é 2D axissimétrico (r, z) — e isso é rigoroso. Silos cilíndricos possuem simetria de revolução. Dados experimentais mostram desvios azimutais de apenas 1–2°C vs. gradientes radiais de 5–8°C (Jian et al. 2009). O erro da aproximação axissimétrica é ~10–20% do gradiente principal.

A percolação em rede 2D — essa sim é simplificação topológica. O limiar pꟲ²ᴰ = 0,593 é cota superior; o limiar 3D seria pꟲ³ᴰ ≈ 0,312 (Stauffer & Aharony, 1994). Na prática, a migração pervasiva no silo real acontece com ~31% dos grãos úmidos, não 59%.

Por que isso importa

O Brasil é o celeiro do mundo. Produzimos mais de 200 milhões de toneladas de grãos por ano. Perder 10% dessa produção por armazenamento inadequado não é apenas um desperdício econômico — é uma falha sistêmica de engenharia que tem solução matemática.

O modelo que apresentamos não é perfeito. É um primeiro passo. A percolação está simplificada, a validação é limitada, o acoplamento é sequencial. Mas ele demonstra algo fundamental: que é possível prever a deterioração semanas antes que qualquer termômetro detecte, usando física e matemática que já existem há décadas — só nunca haviam sido combinadas dessa forma.

O próximo passo é transformar essa matemática em produto. Sensores NDIR de baixo custo estão ficando mais baratos a cada ano. Modelos híbridos CFD + Machine Learning podem absorver dados em tempo real. E a teoria da percolação, quando estendida para redes 3D realistas, pode fornecer alertas topológicos que nenhum sensor isolado conseguiria gerar.

A matemática existe. A física existe. O que falta é integração — e vontade.

Giselle Couto Falcão é PhD em Modelagem Matemática e Computacional e pesquisadora no CEFET-MG, Pesquisa Avançada em Matemática do Clima na Sorbonne Université/Institut Henri Poincaré (IHP) e autora da metodologia CEOD, aplicada na plataforma SerGoias da Fundação Sagres entre 2022 e 2024 em mais de 800 escolas da rede pública estadual de Goiás. Este artigo é uma versão técnica e abreviada do trabalho acadêmico publicado em 2025.