Efeitos numéricos na simulação de
escoamentos gás-sólido em leito �uidizado
borbulhante utilizando a teoria cinética

dos escoamentos granulares

Meire Pereira de Souza

Bauru
2009



Efeitos numéricos na simulação de
escoamentos gás-sólido em leito �uidizado
borbulhante utilizando a teoria cinética

dos escoamentos granulares

Meire Pereira de Souza

Orientador: Hélio Aparecido Navarro

Dissertação apresentada à faculdade de Engenharia da
Universidade Estadual Paulista �Júlio de Mesquista Fi-
lho�, Campus de Bauru, como parte dos requisitos para
obtenção do título de mestre em Engenharia Mecânica.

Bauru
2009



Dedico este trabalho aos meus amados
pais Luiza e Renato (in memoriam).



Agradecimentos

A Deus, pela dádiva da vida, por ter me concedido sabedoria para a realização deste
trabalho, por me abençoar e me confortar nas horas difíceis.

Aos meus pais Luiza e Renato (in memoriam) pelo amor incondicional, pelos ensina-
mentos, e pelo apoio em todos os momentos da minha vida.

Aos meus irmãos Marisa, Amauri, Mauricéia e Renato Fillho que sempre me apoiaram
em tudo, que são e sempre serão meus melhores amigos.

Ao Francisco, pelo apoio, amor e carinho.
Ao Professor Dr. Hélio Aparecido Navarro, pela amizada, orientação e apoio no

desenvolvimento do trabalho.
Ao Professor Dr. Luben Cabezas-Gómez, pelo auxílio com material bibliográ�co e

transmissão de conhecimentos.
Aos Professores do Departamento de Engenharia Mecânica da Unesp, Campus de

Bauru, que contribuíram no desenvolvimento deste trabalho.
Aos amigos, que estiveram presentes e contribuíram direta ou indiretamente na rea-

lização do trabalho.
Aos funcionários da Pós Graduação da FEB, pelo atendimento e colaboração.
Ao Núcleo de Engenharia Térmica e Fluidos do Departamento de Engenharia Mecâ-

nica da EESC/USP, pela acolhida e facilidades oferecidas.
A CAPES, pelo auxílio �nanceiro concedido para a realização do trabalho.
A todos que eu não tenha mencionado, mas que contribuíram para realização deste

trabalho, meu muito obrigada.



�Aprender é a única coisa de que a
mente nunca se cansa, nunca tem
medo e nunca se arrepende�
Leonardo da Vinci



Sumário

Simbologia Principal vi

Resumo ix

Abstract x

1 Introdução 1
1.1 Aplicações dos escoamentos gás-sólido em leitos �uidizados . . . . . . . . . 1
1.2 A Presença da difusão numérica na simulação do escoamento gás-sólido em

leitos �uidizados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Modelagem numérica aplicada a �uidização . . . . . . . . . . . . . . . . . . . . 3
1.4 Objetivos e apresentação do trabalho . . . . . . . . . . . . . . . . . . . . . . . 5

2 Teoria Hidrodinâmica 7
2.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 Formulação Euleriana-Euleriana clássica . . . . . . . . . . . . . . . . . . . . . 9
2.3 Teoria Cinética dos Escoamentos Granulares . . . . . . . . . . . . . . . . . . . 9
2.4 Formulação do modelo das duas fases separadas . . . . . . . . . . . . . . . . . 12
2.5 Modelo hidrodinâmico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.5.1 Leis de fechamento . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

3 Metodologia de Solução Numérica 22
3.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.2 Estudo sobre a discretização de termos de convecção e difusão . . . . . . . . . 23

3.2.1 Esquemas de primeira ordem . . . . . . . . . . . . . . . . . . . . . . . 23
3.2.2 Esquemas de alta ordem . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2.3 Variável normalizada para o limitador de �uxo . . . . . . . . . . . . . . 27
3.2.4 Fator de ponderação �Downwind� . . . . . . . . . . . . . . . . . . . . . 30

i



SUMÁRIO ii

3.2.5 Equação de transporte escalar . . . . . . . . . . . . . . . . . . . . . . . 34
3.3 Solução numérica do modelo hidrodinâmico para escoamento gás-sólido . . . . 38

3.3.1 Equação da conservação da quantidade de movimento . . . . . . . . . . 40
3.3.2 Eliminação parcial do acoplamento na interface . . . . . . . . . . . . . 44
3.3.3 Equação para correção da pressão do �uido . . . . . . . . . . . . . . . . 46
3.3.4 Equação para a correção da fração de volume do sólido . . . . . . . . . 50

3.4 Metodologia computacional . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

4 Resultados 58
4.1 Introdução . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.2 Características do leito para as simulações numéricas . . . . . . . . . . . . . . 58
4.3 Resultados das simulações . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

4.3.1 Comparação entre os dois esquemas de discretização para os termos
convectivos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

4.3.2 Considerações sobre a in�uência da malha computacional nos resultados
de simulação . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

5 Conclusões e recomendações 78
5.1 Considerações gerais . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
5.2 Conclusões . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.3 Recomendações para futuros trabalhos . . . . . . . . . . . . . . . . . . . . . . 80

A Teoremas 89
A.1 Teorema de Leibniz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
A.2 Teorema de Gauss-Ostragradskii . . . . . . . . . . . . . . . . . . . . . . . . . . 89

B "Convection Boundedness Criterion"(CBC) 90
B.1 Critério de Harten . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
B.2 Exemplos de limitadores . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93



Lista de Figuras

1.1 Diferentes regimes de �uidização quando aumenta-se a velocidade do gás. . . . . . . 2
1.2 Leito �uidizado borbulhante com jato central. . . . . . . . . . . . . . . . . . . . . 5

2.1 Método das médias para modelagem de escoamento bifásico Ishii (1975) e Ishii e

Mishima (1984) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 Procedimento geral simpli�cado para formulação da TCEG, baseado no trabalho de

Therdthianwong (1994). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3 Procedimento geral para formulação do modelo de duas fases apresentado por Enwald,

Peirano e Almstedt, (1996) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.4 Regime plástico e viscoso para escoamentos granulares . . . . . . . . . . . . . . . . 16
2.5 Função de distribuição radial g0 Syamlal, Rogers e O'Brien (1993). . . . . . . . . . 18
2.6 Viscosidade dinâmica do sólido Syamlal, Rogers e O'Brien (1993). . . . . . . . . . . 19

3.1 Volume de controle típico em uma malha unidimensional . . . . . . . . . . . . . . 24
3.2 Localização dos nós baseados na direção do escoamento. . . . . . . . . . . . . . . . 28
3.3 Restrições do Limitador Universal no Diagrama de Variável Normalizada. . . . . . . 29
3.4 Restrições do Limitador Universal em termos do fator DWF . . . . . . . . . . . . . 31
3.5 Fatores Downwind em função de

∼
φC . . . . . . . . . . . . . . . . . . . . . . . . . . 32

3.6 Localização dos nós. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.7 Volume de controle para a componente x da equação de quantidade de movimento. . 41
3.8 Condição de Deslizamento livre e de Não Deslizamento na parede leste. . . . . . . . 43
3.9 Condições de contorno de correção da pressão. . . . . . . . . . . . . . . . . . . . . 50
3.10 Interface do usuário MFIX . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

4.1 Geometria e condições iniciais e de contorno usadas na simulação considerando o

sistema de coordenadas cartesianas . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.2 Contorno da fração de vazio (εg), para o esquema FOUP . . . . . . . . . . . . . . 62
4.3 Contorno da fração de vazio (εg) (Bouillard, Gidaspow e Lyczkowski (1991)) . . . . 63

iii



LISTA DE FIGURAS iv

4.4 Contorno da fração de vazio (εg), utilizando o esquema FOUP (Guenther e Syamlal

(2001)) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.5 Contorno da fração de vazio (εg), para o esquema Superbee . . . . . . . . . . . . . 64
4.6 Contorno da fração de vazio (εg), utilizando o esquema Superbee (Guenther e Syamlal

(2001)) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.7 Per�s radiais médio no tempo da velocidade, para a fase sólida 0, 15 metros acima da

entrada do leito . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.8 Velocidade axial média no tempo, 0, 05 metros acima da entrada do leito para os

esquemas FOUP e Superbee . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.9 Velocidade axial média no tempo, 0, 20 metros acima da entrada do leito para os

esquemas FOUP e Superbee . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
4.10 Velocidade axial média no tempo, 0, 29 metros acima da entrada do leito para os

esquemas FOUP e Superbee . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
4.11 Contorno da temperatura granular (θ), em t=0, 36s. . . . . . . . . . . . . . . . . . 69
4.12 Fração de vazio (εg) utilizando o esquema FOUP . . . . . . . . . . . . . . . . . . 70
4.13 Fração de vazio (εg) utilizando o esquema Superbee . . . . . . . . . . . . . . . . . 70
4.14 Fração de vazio (εg), em t = 0, 36s . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.15 Fração de vazio (εg), em t = 0, 59s . . . . . . . . . . . . . . . . . . . . . . . . . . 72
4.16 Fração de vazio (εg), instantes antes de atingir à superfície do leito . . . . . . . . . 73
4.17 Dados experimentais de Kuipers et al. (1993) . . . . . . . . . . . . . . . . . . . . 74
4.18 Velocidade axial média no tempo, para a fase sólida, em várias alturas do leito, para

as quatro malhas distintas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75
4.19 Velocidade axial da fase sólida (Vs) . . . . . . . . . . . . . . . . . . . . . . . . . . 76

B.1 Representação do diagrama CBC em variáveis normalizadas. . . . . . . . . . . . . 91



Lista de Tabelas

3.1 Fórmulas de Discretização em termos do fator DWF. . . . . . . . . . . . . . . . . 31

4.1 Dados numéricos da simulação . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
4.2 Modelo hidrodinâmico utilizado nas simulações . . . . . . . . . . . . . . . . . . . 61
4.3 Tempo de processamento (CPU) utilizado nas simulações . . . . . . . . . . . . . . 77

v



Simbologia Principal

A Área da face do volume de controle, [m2].
CDs Função de arrasto para uma partícula num meio in�nito.
dp Diâmetro de partícula, [m].
=

D Tensor taxa de deformação, [s−1].
e Coe�ciente de restituição para a colisão entre as partículas sólidas.
−→
f Força de iteração entre as fases, [N/m3].
−→g Campo gravitacional, [m/s2].
g0 Função de distribuição radial.
=

I Tensor unitário.
I2D Segundo invariante do tensor taxa de deformação.
pg Pressão da fase gasosa, [Pa].
ps Pressão da fase sólida, [Pa].
pk Pressão da fase k, [Pa].
Res Número adimensional de Reynolds para fase sólida.
tr Traço da matriz.
−→vg Velocidade média local da fase gasosa, [m/s].
−→vs Velocidade média local da fase sólida, [m/s].
um, vm, wm Velocidades nas direções x, y, z da fase gaosa/sólida

no volume de controle, [m/s].
Vrs Correlação para a velocidade terminal para a fase sólida, [m/s].

vi



SIMBOLOGIA PRINCIPAL vii

Símbolos Gregos

β Coe�ciente de arrasto da interface.
Γ Coe�ciente de difusão.
εg Fração de vazio.
εs Fração volumétrica da fase sólida.
θ Temperatura granular, [m2/s2].
λg Viscosidade volumétrica da fase gasosa, [kg/ms].
λs Viscosidade volumétrica da fase sólida, [kg/ms].
µg Viscosidade dinâmica do gás, [kg/ms].
µs Viscosidade dinâmica do sólido, [kg/ms].
ρg Densidade do gás, [kg/m3].
ρs Densidade do sólido, [kg/m3].
ρ
′
m Densidade efetiva macroscópica da fase m (gás ou sólido), [kg/m2].

=
σg Tensor das tensões para fase gasosa, [Pa].
=
σs Tensor das tensões para fase sólida, [Pa].
=
σk Tensor das tensões da fase k, [Pa].
=
τg Tensor das tensões viscosas do gás, [Pa].
=
τs Tensor das tensões viscosas do sólido, [Pa].
φ Propriedade transportada.

Sobrescritos

p Regime plástico.
v Regime viscoso.



SIMBOLOGIA PRINCIPAL viii

Subscritos

g Fase gasosa.
s Fase sólida.
w Face oeste do volume de controle.
e Face leste do volume de controle.
n Face norte do volume de controle.
s Face sul do volume de controle.
b Face da base do volume de controle.
t Face do topo do volume de controle.

Abreviaturas e siglas

EDP Equação diferencial parcial.
IIT Illinois Institute of Technology.
LFB Leito �uidizado borbulhante.
LFC Leito �uidizado circulante.
MFIX Multiphase Flow with Interphase eXchanges.
NETL National Energy Technology Laboratory.
SIMPLE Semi Implicit Method for Pressure Linked Equations.
TCEG Teoria Cinética dos Escoamentos Granulares.



Resumo

Souza, M. P. (2008) Efeitos numéricos na simulação de escoamento gás-sólido em leito
�uidizado borbulhante utilizando a teoria cinética dos escoamentos granulares. Bauru,
2008. 93p. Dissertação (Mestrado) - Faculdade de Engenharia de Bauru, Universidade
Estadual Paulista.

No presente trabalho desenvolve-se um estudo de modelagem matemática e simula-
ção numérica do escoamento bifásico gás-sólido em leito �uidizado borbulhante. Utiliza-se
o modelo Euleriano de duas fases separadas formulando o tensor das tensões da fase sólida
através da teoria cinética dos escoamentos granulares. As simulações numéricas são realizadas
através do código fonte MFIX (Multiphase Flow with Interphase eXchanges) desenvolvido no
NETL (National Energy Technology Laboratory). Os resultados de simulação numérica são
avaliados por meio da análise da in�uência dos seguintes parâmetros: malha computacional e
esquemas de discretização dos termos convectivos das equações de conservação. Com base nos
estudos teóricos e resultados obtidos durante o trabalho conclui-se que esquemas de primeira,
tais como FOUP são altamente difusivos, já os resultados obtidos utilizando o esquema de
alta ordem, Superbee, produziu resultados de melhor qualidade para as malhas testadas neste
trabalho. Além disso, os resultados mostraram-se bastante dependentes do tamanho da malha
computacional.

Palavras chaves: escoamentos bifásicos gás-sólido, modelo das duas fases separadas,
leitos �uidizados borbulhantes, teoria cinética dos escoamentos granulares, MFIX, simulação
numérica, esquemas de discretização, difusão numérica, in�uência da malha computacional.

ix



Abstract

Souza, M. P. (2008) Numerical e�ects on simulation of gas-solid �ow in the bubbling �ui-
dized bed using the Kinetic Theory of Granular Flows. Bauru, 2008. 93p. Dissertação
(Mestrado) - Faculdade de Engenharia de Bauru, Universidade Estadual Paulista.

In the present work is described a mathematical model and numerical simulation of
gas-solid two-phase �ow in a bubbling �uidized bed. It is used the Eulerian gas-solid two-�uid
model and the solid phase stress tensor is modeled considering the kinetic theory of granular
�ows. The numerical simulations were developed using the MFIX (Multiphase Flow with In-
terphase eXchanges) code developed in NETL (National Energy Technology Laboratory). The
numerical di�usion is analyzed considering a single bubbling detachment and its movement
process in a two-dimensional bubbling �uidized bed using the bubble shape as a metric for
results description. The in�uence of computational grid it is also analyzed. It is concluded
that SuperBee scheme produces the better results and analysis about estimating uncertainty
in grid re�nement should be studied.

Keywords: two-phase gas-solid �ow, two-�uid model, bubbling �uidized bed, kine-
tic theory of granular �ows, MFIX, numerical simulation, discretization schemes, numerical
di�usion, computational grid.

x



Capítulo 1

Introdução

1.1 Aplicações dos escoamentos gás-sólido em leitos
�uidizados

A �uidização em escoamentos gás-sólido constitue uma das mais importantes aplica-
ções industriais que envolvem escoamentos bifásico. Entre as variadas aplicações dos escoa-
mentos gás-sólido, duas são mais representativas e importantes economicamente: os craquea-
dores catalíticos para conversão de frações pesadas de petróleo em gasolina e os combustores
de leito �uidizado para geração de energia térmica e elétrica. Embora as reações químicas e os
processos de transferência de calor in�uenciem diretamente as taxas de conversão dos reatores
catalíticos e a e�ciência térmica de instalações energéticas, estes são altamente in�uenciados
pelos processos hidrodinâmicos, os quais determinam a distribuição espacial das fases e espé-
cies envolvidas. Como a hidrodinâmica é dominante nos processos de transporte de massa e
energia em certas escalas temporais e espaciais, é imprescindível o seu estudo e compreensão.
Ainda na atualidade a compreensão dos processos hidrodinâmicos básicos que acontecem em
escoamentos multifásicos nas instalações industriais é incompleta e insu�ciente.

Leitos �uidizados de escoamentos gás-sólido aparecem em diversas aplicações indus-
triais. Um tipo são os reatores de leito �uidizado borbulhante (LFB), e outro os reatores
de leito �uidizado circulante (LFC) que se distinguem pelo regime de �uidização. Os leitos
borbulhantes são caracterizados por altas densidades de particulado, pelo desenvolvimento de
bolhas de gás que promovem recirculação e mistura, e pelo processo de elutriação que pro-
move o arrasto de particulados mais �nos. A maior parte do leito é formada de partículas
cujas velocidades terminais são maiores que a velocidade do gás. Uma característica dos leitos
�uidizados gás-sólido é a presença de bolhas, as quais afetam o desempenho do reator. Assim,

1



CAPÍTULO 1. INTRODUÇÃO 2

é importante entender as suas características e comportamento. Comumente pode-se de�nir
bolhas como regiões com pouco ou sem sólido.

Por operarem em regime denso, os leitos �uidizados borbulhante apresentam uma alta
concentração para a fração de volume dos sólidos (εs), porém com o aumento da velocidade
super�cial do gás, a concentração média dos sólidos diminui, e diferentes regimes de �uidização
são encontrados (Figura 1.1).

Leito �xo Borbulhante Pistonado Turbulento Rápida �uidização

Figura 1.1: Diferentes regimes de �uidização quando aumenta-se a velocidade do gás.

O regime de rápida �uidização prevalecem em leito �uidizados circulantes, onde a
fração de volume dos sólidos é geralmente inferior a 10% (εs < 0.1). Tanto as unidades de
LFB quanto as de LFC possuem várias vantagens, tais como:

• Possibilidade de funcionamento contínuo;

• Condições quase isotérmica durante todo leito do reator causada pela rápida mistura
dos sólidos;

• Excelentes taxas de transferência de calor e massa entre o gás e as partículas suspensas;

• Excelentes taxas de transferência de calor entre o leito e superfícies imersas.



CAPÍTULO 1. INTRODUÇÃO 3

1.2 A Presença da difusão numérica na simulação do es-
coamento gás-sólido em leitos �uidizados

No estudo de escoamentos bifásicos podem-se descrever os leitos �uidizados como
modelos hidrodinâmicos de duas fases os quais tratam o �uido e o sólido como duas fases
separadas contínuas em que todas as partículas são consideradas idênticas caracterizadas por
um diâmetro efetivo e com propriedades de materiais idênticas. Essa idéia de descrever os leitos
�uidizados como modelo hidrodinâmico de dois �uidos existe desde a década de 60. Diversos
autores descreveram esse modelo na época, Davidson (1961) apud Guenther e Syamlal (2001),
Jackson (1963), Murray (1965), Collins (1965), Stewart (1968) apud Guenther e Syamlal
(2001) e Anderson e Jackson (1967) apud Guenther e Syamlal (2001).

O conjunto de equações proposto por esses pesquisadores é muito difícil de se resol-
ver, e soluções numéricas para predizerem bolhas apareceram mais tarde nos trabalhos de
Gidaspow e Ettehadieh (1983) apud Guenther e Syamlal (2001), Syamlal e O'Brien (1989)
apud Guenther e Syamlal (2001), Bouillard, Gidaspow e Lyczkowski (1991), Gidaspow (1994),
Sanyal e Cesmebasi, (1994), Boemer, Qi e Renz (1997), Kuipers et al. (1993). Nos estudos
realizados por esses autores foram observadas, a partir do contorno da fração de vazio, bolhas
de formas alongadas e pontiagudas ao invés de bolhas com formas arredondadas encontradas
experimentalmente e utilizadas nas análises teóricas de Collins (1965) e de Stewart (1968)
apud Guenther e Syamlal (2001). Estes estudos levaram ao questionamento se o modelo de
dois �uidos, é de alguma forma, incompleto que resultaria na predição de bolhas como formas
diferentes das reais ou se a técnica numérica utilizada era inadequada para resolver as equa-
ções corretamente. Recentemente, estudos mostraram que um dos problemas está na técnica
numérica.

1.3 Modelagem numérica aplicada a �uidização

As pesquisas em processos de �uidização gás-sólido tiveram início fundamentalmente
na segunda metade do século XX. Os primeiros trabalhos relacionavam-se principalmente a
estudos experimentais. Mais recentemente, em função da evolução nos campos das técnicas
numéricas e dos recursos computacionais, têm-se desenvolvido consideravelmente as áreas de
modelagem matemática e simulação numérica.

Existem diversas aproximações para os termos convectivos das equações de conser-
vação, entre as quais se destacam os esquemas de primeira ordem tais como �FOUP� (First



CAPÍTULO 1. INTRODUÇÃO 4

Order UPwind) de Courant, Isaacson e Reeves (1952), HYBRID de Spalding (1972) entre
outros. Embora esses esquemas �upwind� de primeira ordem tenham sido muito empregados
por serem incondicionalmente estáveis, eles introduzem uma viscosidade numérica que sua-
viza a solução. Os esquemas de alta ordem como, por exemplo, o de diferença centrada (CD),
upwind de segunda ordem �SOU� de Price, Varga e Warren (1966) e upwind de terceira ordem
�QUICK� (Quadratic Upstream Interpolation for Convective Kinematics) de Leonard (1979)
melhoram a precisão dos cálculos. Entretanto, esses esquemas apresentam problemas quanto
à limitação, ou seja, as soluções podem apresentar oscilações não-físicas nas descontinuidades.

Diversas pesquisas são desenvolvidas para construir esquemas de alta ordem e limi-
tados. Geralmente os esquemas propostos se baseiam em dois conceitos: as restrições TVD
(Total Variation Diminishing) e a abordagem em variáveis normalizadas NVD (Normalized
Variable Diagram). Este conceito normaliza as variáveis e também a variação do sentido do
escoamento. Os esquemas que satisfazem as restrições TVD apresentam soluções livres de
oscilação e convergentes. Esquemas que satisfazem as condições TVD são sempre limitados,
um exemplo é o esquema �MUSCL� (Monotonic Upstream Scheme for Conservation Laws) de
van Leer (1979) em que foi obtida precisão espacial de segunda ordem. Outro esquema de
ordem superior é o �Superbee� de Sweby (1984).

Para problemas estacionários, o esquema �SHARP� (Simple High-Accuracy Resolu-
tion Program) de Leonard, (1988), formulado no contexto de variáveis normalizadas (NVD),
é limitado e possui um bom desempenho. Outro exemplo é o esquema �SMART� (Sharp and
Monotonic Algorithm for Realistic Transport) de Gaskell e Lau (1988), que também utiliza a
formulação em variáveis normalizadas e é linear por partes. Ambos os esquemas �SHARP� e
�SMART� podem alcançar terceira ordem de precisão. Zhu e Leschziner (1988) propuseram
um esquema híbrido chamado �HLPA� (Hybrid Linear Parabolic Approximation Scheme), que
possui precisão de ordem dois, e em 1998 Varonos e Bergeles (1998) propuseram o esquema
�VONOS� (Variable-Order Non-Oscillatory Scheme) .

Utilizando os conceitos TVD e NVD, Song et al. (2000) construíram um esquema
com precisão de ordem três e limitado. O esquema foi designado como �WACEB� (Weighted-
Average Coe�cient Ensuring Boundedness). Resultados numéricos para alguns problemas de
convecção mostraram que este esquema possui a habilidade do esquema QUICK em reduzir
a difusão numérica sem introduzir oscilações. Entretanto, segundo Alves, Oliveira e Pinho
(2003) o esquema WACEB possui problemas de convergência para escoamentos viscoelásticos.
Eles propuseram um esquema para esses tipos de escoamentos denominado �CUBISTA� (Con-
vergent and Universally Bounded Interpolation Scheme for Treatment of Advection) baseado



CAPÍTULO 1. INTRODUÇÃO 5

em restrições TVD e associado ao número de Courant.
Portanto, a técnica numérica pode in�uenciar na solução do problema e para cada

tipo de aplicação um determinado esquema pode-se adequar melhor.

1.4 Objetivos e apresentação do trabalho

A presente pesquisa tem como objetivo investigar os efeitos de difusão numérica
quando esquemas de ordens diferentes são utilizados na discretização dos termos convectivos
nas equações de conservação para problemas de escoamento gás-sólido em leitos �uidizados
borbulhantes.

Para isto, utiliza-se um modelo para a simulação de um leito �uidizado borbulhante
com um jato central, conforme �gura 1.2. Jatos de ar entram na base do leito com uma
velocidade especí�ca e uniforme para manter o leito em um estado de mínima �uidização.
Através de um jato central é injetado gás com velocidade superior provocando a formação de
bolhas.

Figura 1.2: Leito �uidizado borbulhante com jato central.

Os parâmetros analisados no presente trabalho são os esquemas de discretização e
a malha computacional. A avaliação dos resultados de simulação numérica será efetuada
comparando-se parâmetros tais como, fração de vazio, velocidade do gás e do sólido, e tem-
peratura granular, em diferentes seções do leito.



CAPÍTULO 1. INTRODUÇÃO 6

Os experimentos computacionais foram realizados utilizando-se o código MFIX - �Mul-
tiphase Flow with Interphase eXchange�, desenvolvido pelo Departamento de Energia NETL
(National Energy Technology Laboratory). O MFIX é um código cujo modelo hidrodinâmico
descreve reações químicas e transferência de calor em escoamentos densos ou diluídos, escoa-
mentos estes que aparecem tipicamente em processos de conversão de energia e em reatores
químicos. Outras informações do código MFIX podem ser obtidas no trabalho de Syamlal,
Rogers e O'Brien (1993) e na web: www.m�x.org.

O trabalho subdivide-se em cinco capítulos, sendo que no capítulo 1 foi apresentada
uma breve descrição do escoamento gás-sólido, com ênfase às suas principais características e
aplicações. Apresentou-se também uma motivação para o desenvolvimento do trabalho com
uma discussão sobre a difusão numérica na simulação do escoamento gás-sólido em leitos �ui-
dizados borbulhantes. No capítulo 2 é apresentada uma formulação matemática que modela
escoamento gás-sólido considerando o modelo de dois �uidos. Nesta etapa são descritas as
equações de balanço considerando a abordagem Euleriana-Euleriana para ambas as fases, e
também são de�nidas as equações constitutivas para o fechamento do modelo proposto. O
tensor das tensões para a fase sólida é avaliado com base na teoria cinética dos escoamentos
granulares, em que a temperatura granular é calculada por uma equação algébrica. No capí-
tulo 3 é feita uma discussão da metodologia de solução numérica utilizada no código MFIX,
com ênfase na obtenção das fórmulas discretas dos termos de convecção presentes na equação
de transporte de uma propriedade escalar. O capítulo 4 apresenta os resultados de simulação
numérica 2D obtidos a partir do modelo hidrodinâmico descrito no capítulo 2. Comparam-se
os resultados para os esquemas numéricos de discretização dos termos convectivos presentes
nas equações de conservação. Analisa-se a in�uência do tamanho da malha computacional
nos resultados de simulação. No capítulo 5 apresentam-se as conclusões sobre o estudo desen-
volvido e propostas para futuros trabalhos.



Capítulo 2

Teoria Hidrodinâmica

2.1 Introdução

Os fenômenos físicos que ocorrem na modelagem numérica de �uidos utilizam modelos
matemáticos que são fudamentalmente baseados em sistemas de equações diferenciais parciais
(EDP), com condições iniciais e de contorno prescritas. A obtenção deste sistema fechado de
equações segue uma lógica baseada na aplicação de leis físicas de conservação, tais como, a
conservação da massa, a segunda lei de Newton da conservação da quantidade de movimento e
a conservação de energia. A di�culdade de aplicação destas leis, em escoamentos bifásicos deve-
se ao fato de que este tipo de escoamento caracteriza-se pela presença de duas fases diferentes,
e de interfaces que separam essas fases entre si. Sendo assim, dependendo da geometria
das interfaces existem vários regimes de escoamentos da mistura bifásica caracterizados por
diferentes mecanismos de transporte (Ishii (1975)). Tais regimes de escoamentos podem
ocorrer simultaneamente em um sistema simples, di�cultando a modelagem dos escoamentos
bifásicos.

Pode-se contornar teoricamente este problema aplicando localmente as equações de
balanço em cada fase, e modelando apropriadamente as condições de contorno na interface.
Porém, uma formulação geral baseada em variáveis locais instantâneas e interfaces em movi-
mento, resulta em um problema de multifronteiras com a posição de interfaces desconhecidas,
o que torna impraticável para a maioria dos casos a obtenção de modelos matemáticos e por-
tanto de soluções. Ishii e Mishima (1984) resumem os três principais procedimentos adotados
pra desenvolver modelos matemáticos para escoamentos bifásicos ou multifásicos:

1. Modelo da Difusão;

2. Volume de Controle;

7



CAPÍTULO 2. TEORIA HIDRODINÂMICA 8

3. Método das Médias.

Tais procedimentos são tratados a partir de um ponto de vista macroscópio, elimi-
nando parcialmente os detalhes da formulação local instantânea.

Os dois primeiros procedimentos são baseados principalmente em hipóteses, intuição
física e similaridade assumida com escoamentos monofásicos. Por outro lado, o método das
médias é matematicamente mais rigoroso e requer uma larga manipulação das equações. Para
uma descrição de cada método pode-se consultar Ishii e Mishima (1984).

Atualmente o procedimento mais utilizado é o das equações de balanço médias ou
método das médias. Por sua vez, esse procedimento, é subdividido em três grupos, dependendo
do conceito físico básico usado para formular o problema que se analisa. Segundo Ishii (1975),
esses grupos são: médias de Euler, médias de Lagrange e média estatística de Boltzmann.
Esses gupos também se subdividem em vários subgrupos, dependendo da variável utilizada
para estabelecer as médias. Essa classi�cação é mostrada na �gura 2.1.

Método das Médias

Equações de balanço
médias

Procedimento de médias e subgrupos
de cada procedimento

Médias de Euler

Valor médio temporal
Valor médio espacial
Média volumétrica
Média de área
Média na linha
Valor médio estatístico
Valor médio misto

Médias de Lagrange

Valor médio temporal
Valor médio estatístico

Média Estatística
de Boltzmann

Propriedades de
transporte

Figura 2.1: Método das médias para modelagem de escoamento bifásico Ishii (1975) e Ishii e Mishima

(1984)

De acordo com a �gura 2.1, dependendo do procedimento de médias adotado, existem
várias formulações diferentes para um sistema bifásico. Neste capítulo apresenta-se duas destas
formulações: a formulação Euleriana-Euleriana clássica e a formulação de Euler/Boltzmann



CAPÍTULO 2. TEORIA HIDRODINÂMICA 9

(Teoria Cinética dos Escoamentos Granulares - TCEG).

2.2 Formulação Euleriana-Euleriana clássica

A formulação Euleriana-Euleriana clássica permite a obtenção do modelo denominado
modelo das duas fases separadas. Este modelo é geralmente obtido utilizando o procedimento
de médias de Euler, e constitue uma das principais formulações das equações de campo macros-
cópicas para um sistema bifásico. O modelo é formulado considerando cada fase em separado,
em termos de um sistema de equações de conservação de massa, quantidade de movimento,
e energia, para cada fase. Como ambas as fases interagem entre si, aparecem nas equações
de campo termos devidos a essa interação, que especi�cam os transportes de massa, quanti-
dade de movimento e energia através da interface. Análises mais rigorosas sobre a formulação
Euleriana-Euleriana clássica podem ser encontrados nos trabalhos de Delhaye e Achard (1976)
e Delhaye e Achard (1977).

2.3 Teoria Cinética dos Escoamentos Granulares

Na formulação do modelo das duas fases separadas para escoamentos gás-sólido é
necessário especi�car as leis constituivas que modelam o tensor das tensões para cada fase
analisada. Para fases contínuas, como o gases e os líquidos, geralmente este termo é mode-
lado considerando-se a hipótese de �uido Newtoniano, o que é bastante razoável. No caso de
fases dispersas (particulado) assumir essa mesma hipótese constitui uma aproximação muito
grosseira, ainda mais quando a fase dispersa é composta de partículas sólidas. A consideração
da fase sólida como um contínuo e como um �uido Newtoniano implica na procura de um
valor, coerente �sicamente, da viscosidade dinâmica da fase sólida µs, para resolver o sis-
tema de equações de EDPs que compõe o modelo tradicional das duas fases separadas. Para
contornar este problema, muitos pesquisadores na atualidade utilizam a Teoria Cinética dos
Escoamentos Granulares (TCEG). Esta teoria baseia-se na formulação de Euler/Boltzmann e
correspondentemente permite a obtenção de um modelo reológico teórico para a fase sólida,
através da modelagem das �utuações da velocidade das partículas com ajuda de equações
complementares (algébricas ou diferenciais).

A TCEG baseia-se nas similaridades entre o escoamento de um material granular, o
qual compreende uma população de partículas com ou sem gás intersticial, e as moléculas de
um gás (Peirano e Leckner (1998)).



CAPÍTULO 2. TEORIA HIDRODINÂMICA 10

A formulação da Teoria Cinética dos Escoamentos Granulares para o modelo de dois
�uidos, de acordo com Peirano e Leckner (1998), é desenvolvida através de dois procedimentos
fundamentais: o primeiro que considera um escoamento granular seco (sem in�uência do gás
intersticial, de acordo com Ding e Gidaspow (1990)) e o segundo, descrito por Ma e Ahmadi
(1998) que considera o gás intersticial entre as partículas do escoamento. No livro de Gidaspow
(1994) discute exaustivamente o primeiro procedimento. Peirano e Leckner (1998) apresentam
uma detalhada discussão do segundo procedimento aplicado a escoamentos turbulentos gás-
sólido em leitos �uidizados circulantes.

Sinclair e Jackson (1989) foram os primeiros a proporem uma descrição fundamental
das tensões na fase sólida no contexto do modelo das duas fases separadas usando a TCEG.
Estes autores descrevem de maneira sucinta o signi�cado físico da �temperatura do particu-
lado� ou �temperatura granular�, que constitui uma de�nição de temperatura para o sólido
semelhante a temperatura de um gás que normalmente apresenta-se na termodinâmica clássica
e na estatística. Em escoamentos granulares, Jenkins e Savage (1983) designam a temperatura
granular como um terço da velocidade ao quadrado média.

Segundo Sinclair e Jackson (1989), devido ao fato do gás não deslizar livremente pelas
paredes do duto, existe um per�l de velocidade do gás, que induz um per�l de velocidade para
o sólido em razão da força de arrasto que o gás exerce sobre as partículas. Como resultado
deste movimento cisalhante as partículas colidem entre si, gerando uma componente aleatória
do movimento das partículas. Este movimento aleatório ou �utuante é completamente inde-
pendente das �utuações da velocidade do gás, podendo sua magnitude exceder as �utuações
da velocidade do gás, devido a massa das partículas.

Estas �utuações da velocidade do particulado geram uma pressão efetiva na fase
sólida, que junto com uma viscosidade efetiva resistem as tensões cisalhantes do conjunto de
partículas. Tanto a pressão como a viscosidade efetiva, dependem fortemente da �temperatura
granular�, acima de�nida, que pode ser calculada através de uma equação algébrica ou pela
solução de uma EDP separada que representa o balanço da denominada energia pseudotérmica,
ou seja, a energia do movimento aleatório das partículas segundo Agrawal et al. (2001). A
energia cinética do movimento aleatório do particulado é análoga a do movimento térmico
das moléculas em um gás. Essa energia pseudotérmica é gerada pelo trabalho efetuado pelas
tensões de cisalhamento efetiva na fase particulada, dissipada pelas colisões inelásticas entre
as partículas e conduzida como resultado dos gradientes da �temperatura do particulado�.

A formulação da TCEG, pode ser desenvolvida seguindo o procedimento geral sim-
pli�cado mostrado na �gura 2.2. Este procedimento relata as etapas básicas usadas na mo-



CAPÍTULO 2. TEORIA HIDRODINÂMICA 11

delagem da TCEG. Primeiramente, formula-se a equação integro-diferencial de Boltzmann
para a freqüência de distribuição da velocidade f . Esta equação formula o escoamento desde
o ponto de vista microscópio do sistema. Numa segunda etapa, formula-se a equação de
transporte de Maxwell. Esta equação permite formular o escoamento desde o ponto de vista
macroscópio, proporcionando uma equação de transporte hidrodinâmica que possibilita o cál-
culo de uma propriedade ψ transportada (por exemplo, massa, quantidade de movimento ou
energia). Num terceiro passo formula-se o teorema de transporte para escoamentos densos
Jenkins-Savage apresentados em Jenkins e Savage (1983) e em Lun, Savage e Je�rey (1984).
Este teorema baseia-se na equação de transporte para gases densos apresentada por Chapman
e Cowling (1960) e é denominado também como equação de transporte de Maxwell-Chapman,
em Savage (1983), e em Gidaspow (1994), tal teorema generaliza o teorema de transporte
válido para gases e escoamentos de partículas diluídas. Por último formulam-se as equações
constitutivas que permitem calcular os termos não modelados que foram introduzidos nas
equações de transporte macroscópicas de acordo com Therdthianwong (1994).

Equação integro-diferencial
de Boltzmann

Equação de transporte
de Maxwell

Teorema de transporte denso
de Jenkins-Savage

Equações constitutivas

Sistemas de equações da TCEG

Figura 2.2: Procedimento geral simpli�cado para formulação da TCEG, baseado no trabalho de

Therdthianwong (1994).

Após o trabalho pioneiro de Sinclair e Jackson (1989), Ding e Gidaspow (1990) aplica-
ram pela primeira vez a TCEG na simulação bidimensional transiente de um leito borbulhante
a frio considerando um escoamento granular seco com base nos seguintes trabalhos: Chapman
e Cowling (1960), Savage (1983), Jenkins e Savage (1983) e Lun, Savage e Je�rey (1984).
Nesses trabalhos a função distribuição de velocidade do sólido foi modelada através da função



CAPÍTULO 2. TEORIA HIDRODINÂMICA 12

de distribuição de velocidade de Maxwell, o que fez válido este procedimento apenas para
escoamentos gás-sólido densos (leitos borbulhantes).

Segundo Peirano e Leckner (1998) uma das análises mais completas desenvolvida com
o modelo das duas fases separadas utilizando a TCEG tem sido realizada pelo grupo da EDF
(Electricité de France). A TCEG desenvolvida por este grupo considera a in�uência do �uido
intersticial, o que possibilita aplicar o modelo para todas as concentrações de sólido, assim
como modular a turbulência do escoamento. Esta modulação signi�ca considerar a in�uência
da presença das partículas sobre a turbulência da fase gasosa e vice-versa (Crowe (2000)).
As simulações desenvolvidas pelo grupo da EDF utilizam modelos de turbulência completos
para ambas as fases (de segunda ordem), e resolvem quatro equações de transporte acopladas
entre si.

A modelagem da turbulência em escoamentos multifásicos, especi�camente em escoa-
mentos bifásicos gás-sólido não será abordada neste trabalho. Além das di�culdades inerentes
à formulação do fenômeno da turbulência em escoamentos monofásicos, deve-se considerar
esses efeitos na interação entre fases através da interface que as separa.

2.4 Formulação do modelo das duas fases separadas

A formulação do modelo das duas fases tem como idéia geral em formular balanços
integrais de massa, quantidade de movimento e energia para um volume de controle �xo que
contenha ambas as fases. Posteriormente, aplicando-se os teoremas de Leibniz e de Gauss-
Ostrogradskii (Apêndice A), os balanços de integrais vão dar origem à dois tipos de equações:
as equações locais instantâneas para cada uma das fases, e as equações de salto locais ins-
tantâneas, que representam a interação entre as fases. Em seguida, aplica-se o procedimento
de médias de Euler às equações locais, de acordo com a variável escolhida para efetuar as
médias, como mostrado na �gura 2.3. Numa terceira etapa, aplica-se as leis de fechamento
com o objetivo de modelar os termos desconhecidos nas equações de balanço. Finalmente são
estabelecidas condições iniciais e de contorno completando a formulação do modelo das duas
fases. O esquema a seguir ilustra a formulação do modelo das duas fases descrito por Enwald,
Peirano e Almstedt, (1996).

A importância de se ter equações locais instântaneas está no fato de que estas cons-
tituem a base fundamental para todos os modelos bifásicos obtidos a partir de procedimentos
de médias, além do mais, as equações locais instântaneas permitem a modelagem direta de
escoamento separados, tais como escoamentos estrati�cado, em películas ou anulares, segundo



CAPÍTULO 2. TEORIA HIDRODINÂMICA 13

Teoremas de Gauss e de Leibniz

Procedimento de médias

Leis de fechamento

Condições iniciais e de contorno

Balanço integral de massa,
quantidade de movimento e energia

Equações locais instântaneas e
condições de salto

Equações médias de balanço

Sistema fechado de equações
diferenciais parciais

Modelo de duas fases

Figura 2.3: Procedimento geral para formulação do modelo de duas fases apresentado por Enwald,

Peirano e Almstedt, (1996)

Ishii (1975). Já a importância em se aplicar o processo de médias é que através dele eliminam-
se as �utuações locais instântaneas de um escoamento multifásico, possibilitando a obtenção
de valores médios dos movimentos e propriedades do escoamento. Com isso, de acordo com
Ishii (1975) modela-se os aspectos macroscópicos de escoamento, que são de maior interesse
em problemas de engenharia.

O desenvolvimento detalhado desta metodologia voltado para escoamentos gás-sólido
é apresentado por Enwald, Peirano e Almstedt, (1996), com base em inúmeros trabalhos da
literatura, tais como: Ishii (1975), Delhaye e Achard (1976), Delhaye e Achard (1977), Drew
(1983).

2.5 Modelo hidrodinâmico

Nesta seção são descritas as equações gerais médias para um escoamento bifásico
gás-sólido, bem como, as leis de fechamento que são adotadas neste trabalho. O modelo
hidrodinâmico adotado é modelo A de Gidaspow (1994). Neste modelo, o gradiente de pressão
do gás está incluído em cada equação de quantidade de movimento para cada fase.

Em condições isotérmicas, com duas fases (gás e sólido), e sem reações químicas, as
equações da conservação da massa e conservação da quantidade de movimento são expressas
a seguir.



CAPÍTULO 2. TEORIA HIDRODINÂMICA 14

Conservação da Massa para a fase gás.

∂

∂t
(εgρg) +∇.

(
εgρg

→
vg

)
= 0 (2.1)

Conservação da Massa para a fase sólida.

∂

∂t
(εsρs) +∇.

(
εsρs

→
vs

)
= 0 (2.2)

Nas equações de conservação da massa, equações (2.1) e (2.2), o primeiro termo à esquerda
representa a variação temporal de acumulo de massa por unidade de volume e o segundo termo
representa a variação do �uxo de massa convectivo.
Conservação da Quantidade de Movimento para a fase gás.

∂

∂t

(
εgρg

→
vg

)
+∇.

(
εgρg

→
vg
→
vg

)
= εg

→
∇ .

=
σg −

→
f +εgρg

→
g (2.3)

Conservação da Quantidade de Movimento para a fase sólida.

∂

∂t

(
εsρs

→
vs

)
+∇.

(
εsρs

→
vs
→
vs

)
=
→
∇ .

=
σs +εs

→
∇ .

=
σg +

→
f +εsρs

→
g (2.4)

onde, εg e εs são as frações de vazio do gás e do sólido respectivamente; →vg e
→
vs são as velocidades

médias locais das fases gasosa e sólida, respectivamente; =
σg e =

σs são os tensores das tensões
das duas fases;

→
f é a força de interação entre as duas fases por unidade de volume; e

→
g é o

campo gravitacional. Nas equações da conservação da quantidade de movimento (2.3) e (2.4),
o primeiro termo à esquerda representa a taxa líquida de aumento do momento, já o segundo
termo representa a tranferência de momento por convecção.

O modelo das duas fases será obtido através da aplicação das leis de fechamento,
condições iniciais e de contorno ao sistema fechado de equações diferenciais parciais, como
descrito na seção 2.4.

2.5.1 Leis de fechamento

Para se obter o sistema fechado geral das equações diferenciais parciais que modelam
um escoamento gás-sólido é necessário especi�car as leis de fechamento. Nesta seção serão
dicutidas apenas as leis de fechamento apropriadas para escoamentos bifásicos gás-sólidos.

As leis de fechamento são classi�cadas em três tipos: constitutivas, de tranferência e
topológicas.
Leis Constitutivas: especi�cam como as propriedades físicas das fases se relacionam, mas
não descrevem parâmetros de transporte entre as fases na interface. Modelam propriedades
tais como: tensor das tensões, viscosidade dinâmica e viscosidade volumétrica.



CAPÍTULO 2. TEORIA HIDRODINÂMICA 15

Leis de Transferência: são equações empíricas que descrevem as diferentes interações entre
as fases que ocorrem na interface. Modelam as propriedades que ocorrem na interface entre
as fases
Leis Topológicas: descrevem a distribuição espacial de uma variável especí�ca do escoamento
analisado.

A seguir serão de�nidas as leis de fechamento utilizadas na modelagem hidrodinâmica
do escoamento gás-sólido em leito �uidizado borbulhante do presente trabalho.

Tensor das tensões

O tensor das tensões, =
σk, para fase k = (g, s) é dado por:

=
σk= −pk

=

I +
=
τk (2.5)

onde, o primeiro termo da equação (2.5) representa a pressão hidrostática,
=

I é o tensor unitário
e o segundo termo denominado tensor das tensões viscosas refere-se ao movimento.

Para a fase gás, o tensor das tensões viscosas será modelado assumindo a relação
tensão/deformação para um �uido Newtoniano, considerando-se a hipótese de Stokes. Mais
detalhes sobre tratamento de �uidos Newtonianos, pode-se consultar White (1991).

Considerando o gás como um �uido Newtoniano o tensor das tensões viscosas da fase
gás será expresso por: =

τg= 2µg

=

Dg +λgtr(
=

D)
=

I , onde µg e λg são a viscosidade dinâmica e
volumétrica da fase gasosa, respectivamente.

=

D, denominado tensor taxa de deformação é
dado por:

=

Dg=
1

2

[
∇ −→
vg +

(
∇ −→
vg

)T
]

(2.6)

Assumindo a hipótese de Stokes, tem-se que: λg + 2
3
µg = 0. Assim, o tensor das

tensões viscosas da fase gasosa será expresso por:

=
τg= µg

[
∇ →
vg +

(
∇ →
vg

)T

− 2

3
∇. →vg

=

I

]
(2.7)

A formulação do tensor das tensões para a fase sólida é um pouco mais complexa
e será feita adotando as teorias apropriadas da literatura para escoamentos granulares. De
acordo com Jenkins e Cowin (1979) apud Syamlal, Rogers e O'Brien (1993) os escoamentos
granulares podem ser classi�cados em dois regimes distintos: um regime plástico ou lento,
no qual a tensão será causada devido ao atrito entre as partículas, e um regime viscoso ou
rápido, no qual a tensão é causada por colisões. A �gura 2.4 ilustra os diferentes regimes de
escoamentos granulares.



CAPÍTULO 2. TEORIA HIDRODINÂMICA 16

Regime plástico Regime viscoso

Figura 2.4: Regime plástico e viscoso para escoamentos granulares

Nos estudos de Syamlal, Rogers e O'Brien (1993), a teoria cinética dos escoamentos
granulares descrita por Jenkins e Savage (1983) usada para escoamentos dispersos (rápidos)
e a teoria do estado crítico de Schae�er (1987) usada para escoamentos densos (lentos) são
combinadas por meio de uma chave no ponto de compactação crítico para a mínima �uidiza-
ção, alternando a formulação do tensor das tensões da fase entre as duas diferentes relações
constitutivas:

=
σs=




−pp

s

=

I +
=
τs

p se εg ≤ ε∗g

−pv
s

=

I +
=
τs

v se εg > ε∗g
(2.8)

em que ε∗g é a fração de vazio para mínima �uidização, e os sobrescritos p e v indicam os
regimes plástico e viscoso, respectivamente.

Os tensores das tensões no regime plástico são usualmente descritos pela teoria da
mecânica dos solos de acordo com Tuzun et al. (1982). Eles originam-se da fricção entre as
partículas e são descritos por modelos fenomenológicos. Existem dois componentes principais
na teoria dos regimes plásticos. O primeiro é uma função de campo, que de�ne uma superfície
no espaço das tensões, na qual o material irá comportar-se elasticamente (ou irá permanecer
rígido, se a elasticidade for desprezada), e onde o ponto de tensão deverá encontrar-se du-
rante a deformação plástica. A função de campo pode mudar enquanto a variação plástica
ocorre. Quando ela se expande diz-se que há um trabalho de endurecimento, quando ela se
contrai o trabalho é de amolecimento, e quando ela permanece �xa, o material é considerado
perfeitamente plástico.

O segundo componente da teoria dos regimes plásticos é a regra de escoamento dos
materiais, a qual estabelece relações entre os componentes de tensão e taxas do tensor das
tensões. Quando o ponto de tensão alcança as tensões plásticas da superfície, ocorrem que
em cada ponto da superfície os componentes dos incrementos plásticos da tensão são �xadas
pela regra do escoamento. Mais detalhes sobre regras de escoamento podem ser encontrados
em Tuzun et al. (1982).

Similarmente as funções usadas nas teorias de escoamento plásticos no trabalho de



CAPÍTULO 2. TEORIA HIDRODINÂMICA 17

Jenike (1987), de�ne-se uma função arbitrária que permite certa quantidade de compressibli-
dade na fase sólida, e representa os termos da pressão do sólido para o regime de escoamento
plástico:

pp
s = εsp

∗ (2.9)

onde p∗ é representado por uma lei empírica:

p∗ = εsA
(
ε∗g − εg

)n (2.10)

com valores típicos das constantes de A = 1025 e n = 10.
A formulação para o tensor das tensões é baseada nos estudos de Schae�er (1987) que

propôs uma equação relacionada com a pressão do sólido:

=
τs

p
= 2µp

s

=

Ds (2.11)

onde
µp

s =
p∗sinφ

2
√
I2D

(2.12)

e φ é o ângulo da partícula no atrito interno e I2D é o segundo invariante do tensor taxa de
deformação:

I2D =
1

6

[
(Ds11 −Ds22)

2 + (Ds22 −Ds33)
2 + (Ds33 −Ds11)

2] +D2
s12 +D2

s23 +D2
s31 (2.13)

Os termos viscosos na equação do tensor das tensões da fase sólida são baseados em
uma forma modi�cada da teoria cinética de partículas esféricas, inelásticas e suaves desenvol-
vida por Lun, Savage e Je�rey (1984). A pressão da fase sólida é expressa por:

pv
s = K1ε

2
sθ (2.14)

com K1 = 2(1 + e)ρsg0 , onde e é o coe�ciente de restituição para colisão da fase sólida e g0 é
a função de distribuição radial no contato das partículas, dada por:

g0 =
1

εg
+

1, 5εs
ε2g

(2.15)

O coe�ciente de restituição e, é a medida da elasticidade de colisão entre duas par-
tículas, e refere-se a quantidade de energia cinética das partículas antes e após a colisão.
Uma colisão perfeitamente elástica possue coe�ciente de restituição igual a 1, e uma colisão
completamente plástica ou inelástica possue coe�ciente de restituição igual a 0. Em geral
experimentos mostram que o coe�ciente de restituição não é constante (Du et al. (2006)).

Goldschmidt, Kuipers e Van Swaaij, (2001), mostraram que a hidrodinâmica dos
leitos �uidizados para gases densos dependem fortemente da energia dissipada nos encontros



CAPÍTULO 2. TEORIA HIDRODINÂMICA 18

de partículas-partículas. Lu et al. (2004), mostraram que a temperatura granular aumenta
com o incremento de e. Estes e mais outros resultados atestam que as simulações de leitos
�uidizados são sensíveis ao coe�ciente de restituição e.

A função de distribuição radial g0, pode ser vista como uma medida para a probabili-
dade do contato entre as partículas segundo Lu et al. (2004). No escopo da teoria cinética, a
função de distribuição radial é introduzida para considerar o incremento da probabilidade das
colisões quando o gás torna-se denso. Isto signi�ca que num gás rarefeito g0 é igual a unidade,
mas tende ao in�nito quando as moléculas estão bem perto, de maneira que o movimento
seja quase impossível (Peirano e Leckner (1998)). Na literatura existe uma concordância
com o fato de que g0 deve aumentar com o aumento da fração volumétrica dos sólidos. Uma
ilustração deste aumento é mostrada na �gura 2.5, que representa a função de distribuição
radial de Syamlal, Rogers e O'Brien (1993), a qual será adotada no presente trabalho.

Figura 2.5: Função de distribuição radial g0 Syamlal, Rogers e O'Brien (1993).

O tensor das tensões da fase sólida para o regime viscoso é expresso por:

=
τs

v
= 2µv

s

=

Ds +λv
str

( =

Ds

) =

I (2.16)

onde
=

Ds é o tensor taxa de deformação para a fase sólida dado por:
=

Ds=
1

2

[
∇ →
vs +

(
∇ →
vs

)T
]

(2.17)

e µv
s e λv

s são as viscosidades dinâmica e volumétrica da fase sólida no regime viscoso, respec-
tivamente, expressos por:

λv
s = K2εs

√
θ (2.18)



CAPÍTULO 2. TEORIA HIDRODINÂMICA 19

µv
s = K3εs

√
θ (2.19)

a cosntante K2, é dada por:

K2 =
4

3
√
π
dpρs(1 + e)εsg0 − 2

3
K3 (2.20)

e a constante K3 é:

K3 =
dpρs

√
π

6(3− e)

[
1 +

2

5
(1 + e)(3e− 1)εsg0

]
+

8dpρsεs
10
√
π
g0(1 + e) (2.21)

Visto que a fase sólida, não é contínua como o gás ou o líquido, sua viscosidade dinâ-
mica não é um propriedade física real. Porém, de acordo Boemer et al. (1995), a viscosidade
dinâmica do sólido pode ser derivada através da TCEG. Como descrito por Gidaspow (1994),
o tensor das tensões compõe-se de duas parcelas, a que representa a contribuição devido as
colisões e a que considera a contribuição cinética. Estas duas parcelas estão presentes nas
correlações para a viscosidade dinâmica do sólido, µs. A �gura 2.6, representa a correlação de
Syamlal, Rogers e O'Brien (1993), a qual será usada neste trabalho. Nesta �gura, a viscosi-
dade µs é calculada em função da vazão volumétrica de sólidos assumindo partículas esféricas,
com diâmetro médio de particulado de 0, 52 mm, densidade de particulado de 2.620 kg/m3,
considerando a correlação para g0, conforme equação (2.15), coe�ciente de restituição de 0, 9

e temperatura granular igual a 0, 1 m2/s2.

Figura 2.6: Viscosidade dinâmica do sólido Syamlal, Rogers e O'Brien (1993).

Pode-se utilizar uma expressão algébrica para a temperatura granular, θ, obtida da
equação da conservação da energia de Lun, Savage e Je�rey (1984), assumindo que a ener-
gia granular é dissipada localmente, desprezando contribuições de difusão e de convecção, e



CAPÍTULO 2. TEORIA HIDRODINÂMICA 20

retendo somente os termos de dissipação e geração, a equação da energia granular algébrica é:

θ =



−(K1εs + ρs)tr

( =

Ds

)

2εsK4

+

+

√[
(K1εs + ρs)

2 tr2
( =

Ds

)
+ 4K4εs

[
2K3tr

(
=

Ds
2

)
+K2tr2

( =

Ds

)]]

2εsK4




(2.22)

onde K4 é dado por:
K4 =

12 (1− e2) ρgg0

dp

√
π

(2.23)

Forças de interação

Os mecanismos e a formulação de forças de interação foram revistos em detalhes por
diversos autores. Dos estudos na dinâmica de uma única partícula em um líquido, diversos
mecanismos diferentes foram identi�cados:

• Força de arrasto: causada pela diferença de velocidade entre as fases.

• Força buoyancy: causada pelo gradiente de pressão do �uido

• Força massa virtual: força requerida para a aceleração da massa da fase circundante,
quando a fase é acelerada.

• Força Sa�man: força ascencional (lift force) que atua sobre a partícula, devido a exis-
tência de gradiente de velocidade.

• Força Magnus: força ascensional (lift force) causada pela existência de gradiente de
pressão provocado por diferencial de velocidade devido a rotação da partícula.

• Força Basset: esta força surge como resultado da difusão de quantidade de movimento
através da camada limite da fase analisada. Resulta da aceleração de uma fase com
respeito à outra e avalia o efeito da aceleração já ocorrida sobre a resistência que a fase
oferece ao movimento.

Das forças em questão, no presente estudo, considera-se que a principal contibuição
das forças de interação entre as fases é a força de arrasto. Ela pode ser modelada através de dois
procedimentos gerais: o primeiro é feito a partir de correlações para o coe�ciente de arrasto
sobre uma partícula numa suspensão de partículas, CDs . O segundo, é a partir da queda de



CAPÍTULO 2. TEORIA HIDRODINÂMICA 21

pressão por unidade de comprimento numa suspensão de partículas, ∆P�L. O objetivo da
modelagem da força de arrasto é expressá-la em função das variáveis dependentes do sistema
de equações governantes. Independente do procedimento adotado, a força de arrasto é de�nida
em função da velocidade relativa entre as fases, e do coe�ciente de transferência de quantidade
de movimento na interface, denominado função de arrasto (Cabezas-Gómez (1999)). Mais
detalhes sobre os procedimentos de obtenção das funções de arrasto são apresentados em
Enwald, Peirano e Almstedt, (1996) e em Cabezas-Gómez (1999).

Portanto, a expressão para a força de arrasto entre as fases gás e sólida será dada
por:

→
f= β

(→
vg − →

vs

)
, onde β é o coe�ciente de arrasto na interface.

Neste trabalho será considerada a correlação de arrasto de Syamlal, Rogers e O'Brien
(1993). Os autores propuseram um modelo baseado em medidas de velocidades das partículas
nos leitos �uidizados na forma:

β =
3

4

CDs

V 2
rs

∣∣∣→vg − →
vs

∣∣∣
dp

εgεs (2.24)

onde Vrs é a velocidade terminal para a fase sólida.
E a expressão de Garside e Al-Dibouni (1977) é usada:

Vrs = 0, 5

(
A− 0, 06Res +

√
(0, 06Res)

2 + 0, 12Res(2B − A) + A2

)
(2.25)

onde:
A = ε4,14

g , e (2.26)

B =





0, 8ε1,28 se εg ≤ 0, 85

ε2,65
g se εg > 0, 85

(2.27)

o número de Reynolds para a fase sólida é dado por:

Res =
εgρg

∣∣∣→vg − →
vs

∣∣∣ dp

µg

(2.28)

e, CDs é uma função de arrasto para uma esfera:

CDs =

(
0, 63 + 4, 8

√
Vrs

Res

)2

(2.29)



Capítulo 3

Metodologia de Solução Numérica

3.1 Introdução

Uma grande variedade de fenômenos físicos, dentre eles o escoamento de �uidos podem
ser descritos através de equações diferenciais parciais (EDP). No entanto, quando tais equações
envolvem não linearidade em suas formulações, a obtenção de soluções analíticas torna-se
rara. Porém com a crescente disponibilidade dos computadores surgiu uma nova alternativa
para o estudo do escoamento de �uidos: as equações que descrevem esses escoamentos são
resolvidas numericamente com o auxílio de computadores. Esse novo ramo do conhecimento,
recebeu o nome de dinâmica de �uidos computacional (DFC), e tem como principal objetivo
complementar, e não substituir os estudos teóricos e experimentais sobre o movimento de
�uidos.

Para se obter uma solução aproximada, faz-se o uso de um método de discretização
que aproxima as equações diferenciais por um sistema de equações algébricas, que pode ser
resolvido no computador. As aproximações são aplicadas a pequenos domínios no espaço e/ou
tempo tal que a solução numérica forneça resultados em locais discretos no espaço e no tempo.
Muito embora, a precisão dos dados experimentais dependa da qualidade das ferramentas
utilizadas, a precisão da solução numérica é dependente da qualidade da discretização utilizada
(Ferziger e Peric (1999)).

De forma particular, a maneira pela qual as derivadas convectivas nas equações de
conservação são aproximadas requer atenção especial, porque esses termos são responsáveis
por causar di�culdades na solução numérica do sistema de equações, provocando a obtenção
de soluções incorretas ou inexatas e/ou a divergência do método de solução. Via de regra, o
processo de discretização introduz termos de segunda ordem adicionais nas equações de ba-
lanço discretizadas em relação às equações de balanço contínuas, que produzem a denominada

22



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 23

difusão numérica. Essa difusão numérica introduz uma dissipação na solução numérica muito
maior do que a produzida pelo coe�ciente de viscosidade molecular do �uido, resultando em
soluções incorretas devido à solução de um problema matemático com características com-
pletamente diferentes do problema físico. Ferziger e Peric (1999), apresentam uma discussão
extensa sobre a de�nição e signi�cado da difusão numérica.

Neste capítulo, serão apresentadas as deduções das fórmulas para a discretização dos
termos transiente, convectivos, difusivos e do termo fonte, a�m de se obter um conjunto de
equações algébricas (discretizadas) para as equações diferenciais que modelam o escoamento
bifásico gás-sólido.

3.2 Estudo sobre a discretização de termos de convecção
e difusão

3.2.1 Esquemas de primeira ordem

As equações de transporte que contém termos convectivos-difusivos são da forma:

ρ
∂φ

∂t
= ρu

∂φ

∂x
− ∂

∂x

(
Γ
∂φ

∂x

)
(3.1)

A estabilidade e precisão do esquema numérico dependem fortemente do método uti-
lizado para discretizar tais termos. É natural discretizar os termos utilizando uma expansão
em série de Taylor. Entretanto, em dinâmica dos �uidos computacional, o método de volume
de controle (VC) é usualmente mais preferido. O método de VC invoca as bases físicas da
dedução das equações de conservação e garante a conservação global de massa, quantidade de
movimento e energia (Patankar (1980)). Na resolução numérica em uma malha su�ciente-
mente �na os dois métodos resultariam na mesma precisão de solução. Em malhas grossas,
o método VC é mais adequado do que expansão em série de Taylor porque este último apre-
senta um erro de truncamento não desprezível nestes tipos de malhas. Assim, o método VC
é atraente, pois uma malha �na nem sempre é interessante devido ao custo computacional.

Para a aplicação da técnica de volumes �nitos, considere o domínio mostrado na �gura
3.1, e a equação de transporte (3.1).

A equação (3.1) é integrada no volume de controle em torno do ponto P :
∫

V C

ρ
∂φ

∂t
dV =

∫

V C

ρu
∂φ

∂x
−

∫

V C

∂

∂x

(
Γ
∂φ

∂x

)
dV (3.2)

As fronteiras laterais do VC são denominadas e e w e elas também são os limites de



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 24

i− 1

W

i− 1/2

w

i

P

i+ 1/2

e

i+ 1

E

x

∆xP ∆xE

∆xe

VC

Figura 3.1: Volume de controle típico em uma malha unidimensional

integração. Portanto:
∫ e

w

ρ
∂φ

∂t
dx =

∫ e

w

ρu
∂φ

∂x
dx−

∫ e

w

∂

∂x

(
Γ
∂φ

∂x

)
dx (3.3)

Pelo teorema fundamental do cálculo, tem-se que a primeira integral da equação (3.3)
sobre o volume de controle VC é dada por:

∫ e

w

ρ
∂φ

∂t
dx = [Ae − Aw]

[
(ρφ)P − (ρφ)old

P

]
(3.4)

onde A é a área da seção transversal do volume de controle.
O termo de convecção-difusão, também é integrado sobre o volume de controle, o que

resulta em:
∫ e

w

[
ρu
∂φ

∂x
− ∂

∂x

(
Γ
∂φ

∂x

)]
dx =

[
ρuφe −

(
Γ
∂φ

∂x

)]
Ae −

[
ρuφw −

(
Γ
∂φ

∂x

)]
Aw (3.5)

O passo seguinte é determinar, em função das incógnitas nos pontos discretos os
valores de

(
Γ∂φ

∂x

)
e
e

(
Γ∂φ

∂x

)
w
.

Os �uxos difusivos ∂φ
∂x

são aproximados por um simples esquema de interpolação de
segunda ordem de precisão. Assim, tem-se:

(
Γ
∂φ

∂x

)

e

= Γe
(φE − φP )

∆xe

e
(

Γ
∂φ

∂x

)

w

= Γw
(φP − φW )

∆xw

(3.6)

Para malhas irregulares, os valores do coe�ciente de difusão Γ nas faces e e w são
calculados utilizando uma média harmônica dos valores de�nidos no centro das células segundo
Patankar (1980). Por exemplo, Γe será dado por:

Γe =

[
1− fe

ΓP

+
fe

ΓE

]−1

=
ΓP ΓE

feΓP + (1− fe) ΓE

onde fe =
∆xe

∆xP + ∆xE

(3.7)

A discretização dos termos de convecção é uma tarefa fundamental. Em escoamentos
nos quais a convecção tem papel importante, como aqueles com altos números de Reynolds,



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 25

a adequada discretização dos termos convectivos é de extrema importância para a qualidade
da solução numérica.

A partir da equação (3.5), a discretização do termo de convecção equivale-se a deter-
minar um valor ou expressão para φ nas faces e e w do VC. Uma interpolação linear entre os
valores adjacentes de φ resulta em:

φe =
φE + φP

2
e φw =

φP + φW

2
(3.8)

denominada de média simples, tal interpolação origina resultados de segunda ordem de pre-
cisão. Entretanto, em escoamentos dominados por convecção, é o caso dos escoamentos gás-
sólido, este método introduz oscilações na solução e conduz a um esquema numérico instável
(Syamlal (1998)). Uma solução para estabilizar os cálculos é o esquema de discretização
�upwind�, conhecido na literatura como FOUP (First Order Up Winding).

φe =





φP se u ≥= 0

φE se u ≤= 0
e φw =





φW se u ≥= 0

φP se u ≤= 0
(3.9)

A discretização �upwind� é de primeira ordem de precisão, e em geral introduz uma forte
difusão numérica na solução. Apesar disso, por não produzir soluções oscilatórias, dispersivas
ela ainda é encontrada na literatura.

A motivação para o estudo de esquemas numéricos é encontrada na solução analítica
para um escoamento unidimensional em regime permanente:

φ− φP

φE − φP

=
exp

(
Pex
∆xe

)
− 1

exp (Pe)− 1
(3.10)

onde Pe = ρu∆xe

Γ
é o número de Péclet e indica a razão entre o �uxo convectivo e o difusivo,

na face leste do volume de controle. Quanto maior o número de Péclet, maior a intensidade
da convecção.

3.2.2 Esquemas de alta ordem

Da literatura, sabemos que os esquemas numéricos �upwind� de primeira ordem, apre-
sentam em geral uma forte difusão numérica. Com o objetivo de minimizar a introdução
excessiva de viscosidade numérica, e ao mesmo tempo derivar técnicas de alta ordem de con-
sistência, diversos esquemas �upwind�, com ordem igual ou superior a dois surgiram na década
de setenta. Dentre os mais populares, podem-se destacar os algoritmos SUDS - �Skew-Upwind
Di�erencing Schemes� de Raithby (1976) e QUICK - �Quadratic Upstream Interpolation for



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 26

Convective Kinematics� de Leonard (1979). Entretanto, nas vizinhanças de regiões onde gra-
dientes elevados aparecem a natureza dispersiva desses esquemas podem induzir oscilações
espúrias (�wiggles�) na solução, as quais são su�cientes para causar instabilidade numérica e,
em decorrência, deteriorar completamente a precisão do método numérico.

As di�culdades em superar oscilações/instabilidades na solução numérica levaram
ao desenvolvimento de esquemas localmente limitados, e muitas tentativas foram feitas para
encontrar um compromisso aceitável entre precisão, estabilidade e simplicidade de implemen-
tação. Hoje em dia, os esquemas localmente limitados consistem de uma das metodologias
mais bem aceitas em escoamentos incompressíveis, pois proporcionam alta ordem do esquema
QUICK e as boas qualidades de estabilidade do esquema FOUP, resultando em esquemas
incondicionalmente estáveis e de ordem de consistência razoável (Ferreira (2001)).

Atualmente, dois conceitos têm demonstrados ser de grande utilidade na construção
de esquemas de ordem elevada: as restrições TVD -�Total Variation Diminishing� e as apro-
ximações NVD -�Normalized Variavel Diagram�. O mecanismos comum e fundamental destes
dois modelos é o uso da natureza dissipativa da diferenciação de primeira ordem. Esquemas
TVD, isto é, esquemas que satisfazem as restrições TVD possuem características atrativas. As
soluções obtidas por eles são bem resolvidas, livres de oscilações e convergentes. Usualmente,
os esquemas TVD são aplicados à escoamentos compressíveis nos quais as variáveis convecti-
vas sofrem grandes alterações no seu gradiente ou no caso da presença de descontinuidades.
Segundo Kaibara, Ferreira e Navarro (2004) esquemas que satisfazem as condições TVD são
sempre limitados .

Esquemas TVD

Dada uma seqüência de aproximações discretas φ(t) = {φi(t)}i∈Z para um escalar, a
Variação Total (TV- �Total Variation�) em um nível t desta seqüência é de�nida por:

TV (φ(t)) =
∑

i∈Z
|φi+1(t)− φi(t)| (3.11)

Após Harten (1984), esquemas de diferenças que originam uma TV limitada são
denominados esquemas TVD. Uma propriedade desejável para uma solução aproximada com-
partilhar com uma exata é que sua TV diminua no tempo (Tadmor (1988)). TVD é uma
propriedade puramente escalar, ela garante que oscilações indesejáveis são removidas da solu-
ção numérica de uma lei de conservação não linear.

Formalmente, considere o esquema de diferenças explícito envolvendo (2k+1) pontos



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 27

da forma:
φn+1

i = H
(
φn

i−k, ..., φ
n
i+k

)
,∀n ≥ 0, i ∈ Z (3.12)

em que H : R2k+1 −→ R é uma função contínua e φn
i denota uma aproximação da solução

exata φ em uma malha uniforme (xi, tn), sendo xi = i∆x, tn = n∆t, com ∆x e ∆t os passos
espacial e temporal respectivamente. Por de�nição, temos que o esquema acima é TVD se:

TV
(
φn+1

) ≤ TV (φn) (3.13)

Se as seguintes propriedades de monotonicidade são mantidas como uma função de t:

• Nenhum extremo local em x pode ser criado;

• O valor de um mínimo local é não decrescente, o valor de um máximo local é não
crescente.

Então, diz-se que o esquema preserva monotonicidade. Esta condição estabelece que um
esquema está preservando monotonicidade se φn+1 permanece monótono quando φn é monó-
tono. Em outras palavras, per�s monótonos são preservados durante a evolução temporal das
soluções discretas e superestimativas não serão criadas.

Os esquemas TVD utilizam um limitador de �uxo para o valor de φ na vizinhança
da face de um volume de controle (VC) quando a variação local de φ é monótona, prevenindo
assim a introdução de oscilações espúrias na solução.

3.2.3 Variável normalizada para o limitador de �uxo

A avaliação de termos convectivos requerem uma expressão para φ, em cada face de
uma célula. Como dito anteriormente, um método amplamente utilizado em escoamentos gás-
sólidos é o FOUP, que embora seja estável é apenas de primeira ordem e portanto altamente
difusivo. Por outro lado, técnicas de discretização de alta ordem produzem oscilações não
físicas que podem causar instabilidade numérica e destruir a precisão do método numérico.
Uma maneira de contornar-se este problema é aplicação de um limitador de �uxo. Leonard,
(1988) e Gaskell e Lau (1988), propuseram o conceito de variável normalizada e critério CBC -
�Convection Boundedness Criterion� (apêndice B) para se obter esquemas capazes de resolver
gradientes elevados e ao mesmo tempo manter instabilidades na solução numérica. Leonard e
Mokhtari (1990), introduziram um �limitador universal� que pode ser expresso em função de
um valor normalizado de φ, em qualquer vizinhança da face do VC, de�nido como:

∼
φ=

φ− φU

φD − φU

(3.14)



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 28

onde os subescritos U e D, representam os valores à montante e à justante respectivamente,
baseados na direção do escoamento (�gura 3.2), e C é o valor do nó entre U e D.

U C D

f

Figura 3.2: Localização dos nós baseados na direção do escoamento.

Por de�nição
∼
φU= 0 e

∼
φD= 1. E também, a distribuição local de φ é monótona

quando 0 ≤
∼
φC≤ 1. Sobre as condições de monotonicidade, as restrições para o limitador

universal em
∼
φf podem ser escritas como:

1. φf está entre φC e φD. Ou seja:

∼
φC≤

∼
φf≤ 1 para 0 ≤

∼
φC≤ 1 (3.15)

isto inclue o caso especial quando φC = φD. Neste caso, φf = φC = φD, isto é,
∼
φf= 1

para
∼
φC= 1.

2. Se φC = φU , então φf = φC = φU . Isto é:

∼
φf= 0 para

∼
φC= 0 (3.16)

3. Para evitar a não unicidade quando
∼
φC→ 0 a fronteira OB tem uma pequema inclinação

positiva e �nita. Isto introduz uma restrição adicional:

∼
φf=

∼
φC

c
para 0 ≤

∼
φC≤ c (3.17)

em que c é uma constante. Para esquema com variação temporal c é a direção normal
do número de Courant

(
u∆t
∆x

)
.

4. Para comportamentos não monótonos
( ∼
φC< 0 ou

∼
φc> 1

)
, o limitador não impõe qual-

quer restrição, com exeção das interpolações que devem ser contínuas com respeito a
∼
φC , isto é, φf é a curva que deve passar por (0, 0) e (1, 1), com ∂

∼
φf

∂
∼

φC

> 0 e �nita.

As limitações acima, podem ser representadas em um �Diagrama Variável Normali-
zado� (NVD), isto é, no plano

∼
φC x

∼
φf , conformme a �gura 3.3.

O valor de φf calculado por quaisquer regimes de ordem superior deve atravessar a
região sombreada a�m de prevenir super ou subestimativa da solução numérica. Além disso,



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 29

∼
φf

∼
φC

O

B A

1

Figura 3.3: Restrições do Limitador Universal no Diagrama de Variável Normalizada.

Leonard (1987) mostrou que uma condição su�ciente e necessária para esquemas de segunda
ordem é que a representação funcional de φf deve passar pelo ponto (0, 5; 0, 75). Métodos de
ordem superiores a dois, não podem ser representados como uma única curva sobre o NVD.
Mais detalhes sobre o diagrama de variável normalizada para limitador de �uxo podem ser
encontrados em Leonard (1979) e Leonard (1991).

Em cada estágio de uma marcha pseudo-temporal ou solução iterativa, as restrições
para o limitador universal são aplicadas como segue:

(i) Para cada face do VC, identi�que a direção da componente da velocidade normal atual,
portanto, identi�cando φU , φC e φD.

(ii) Calcule explicitamente uma tentativa para o valor de face com alta ordem, multidimen-
sional, �upwind� baseado na estimativa de φf

(iii) Calcule o correspondente valor de face normalizado,
∼
φf , e o valor nodal adjacente

�upwind�,
∼
φC .

(iv) Se o ponto
( ∼
φC ,

∼
φf

)
não pertende à região triangular da �gura 3.3, proceda para a

próxima face do VC.

(v) Se não,
∼
φf está limitado e a uma vizinhança apropriada da restrição da fronteira para

um dado valor de
∼
φC .

(vi) O valor da face não normalizado é então reconstruído por: φf =
∼
φf (φD − φU) + φU .

(vii) O mesmo procedimento é utilizado para cada face do VC.



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 30

Então, em um algoritmo de marcha temporal, os �uxos convectivos (limitados) e
correspondentes difusivos (calculados utilizando uma diferenciação central de segunda ordem)
estão agora disponíveis para o passo de atualização explícita. Alternativamente, soluções itera-
tivas implícitas podem ser implemantadas pela introdução do fator ponderado de �downwind�,
descrito por Leonard e Mokhtari (1990).

3.2.4 Fator de ponderação �Downwind�

O limitador universal de Leonard e Mokhtari (1990) pode ser expresso em termos
do �Downwind Weighting Factor� (DWF), que ao invés de limitar diretamente os valores na
face do VC, é introduzido como uma variável auxiliar, gerando assim um esquema compacto
implícito conveniente para métodos de solução de sistema tridiagonais. Após calcular explici-
tamente uma estimativa para φf , de�ne-se:

DWF =
φf − φC

φD − φC

(3.18)

O que é equivalente à:

DWF =

∼
φf −

∼
φC

1−
∼
φC

(3.19)

Em termos do DWF , as restrições para o limitador universal, referentes à �gura 3.3, tornam-
se:

0 ≤ DWF ≤ 1 para 0 <
∼
φC≤ 1 (3.20)

A �gura 3.4 ilustra um diagrama das restrições do limitador universal em termos do fator
DWF. Note que o ponto A dado por (1, 1) na �gura 3.3 foi estendido numa linha vertical na
�gura 3.4

A partir da de�nição do fator DWF, claramente, o valor de φf na face do VC é dado
por:

φf = DWFφD + (1−DWF )φC (3.21)

Para vários esquemas de discretização, pode-se deduzir fórmulas explícitas para o
fator Downwind. As fórmulas DWF utilizadas no MFIX são mostradas na tabela 3.1, como
uma função de γ, onde γ ≡

∼
φC /

(
1−

∼
φC

)
. Algumas destas fórmulas são plotadas nas �guras

3.5.



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 31

DWF

∼
φCO

B A

1

Figura 3.4: Restrições do Limitador Universal em termos do fator DWF .

Tabela 3.1: Fórmulas de Discretização em termos do fator DWF.
Esquema de Discretização DWF
FOUP 0

Diferença Central 0, 5

Upwind de segunda ordem 1
2
γ

Esquemas TVD 0 se
( ∼
φc< 0 ou

∼
φC> 1

)
senão

Van Leer
∼
φC

Minmod 1
2
max[0,min(1, γ)]

MUSCL 1
2
max[0,min(2γ; 0, 5 + 0, 5γ; 2)]

UMIST 1
2
max[0,min(2γ; 0, 75 + 0, 25γ; 2)]

SMART 1
2
max[0,min(4γ; 0, 75 + 0, 25γ; 2)]

Superbee 1
2
max[0,min(1, 2γ),min(2, γ)]

ULTRA-QUICK 1
2
γ para

∼
φC< 0

(
1
c
− 1

)
γ para 0 ≤

∼
φC<

3c
8−6c

0, 375 + 0, 125γ para 3c
8−6c

≤
∼
φC<

5
6

1 para 5
6
≤

∼
φC≤ 1

0, 5 para
∼
φC> 1



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 32

∼
φC

DWF

Figura 3.5: Fatores Downwind em função de
∼
φC .

Esquemas com valores maiores (em um dado valor de φC normalizado),
∼
φc tendem a

ser mais compressivos do que esquemas com valores menores. Portanto o esquema Superbee
é mais compressivo que o MUSCL, por exemplo, por outro lado, o esquema Minmod é mais
difusivo.

Algoritmo para a utlização do fator �Downwind�

Por conveniência de programação, calcula-se um fator de convecção ponderado (ξ)

a partir do fator �Downwind�, o qual pode ser calculado uma vez e utilizado sem checagem
adicional na direção do �uxo. O fator de convecção ponderado ξ é de�nido como:

ξ =





DWF se u ≥ 0

1−DWF se u < 0
(3.22)

Tal método é ilustrado através do cálculo de ξ na face leste, para as outras faces
o procedimento é análogo. As localizações dos nós são mostradas na �gura 3.6. Refere-se
também a �gura 3.2 para as de�nições dos nós C, D e U .



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 33

W P

e

E EE

Figura 3.6: Localização dos nós.

1. Algoritmo 3.1 (Cálculo de
∼
φC):

Se ue ≥ 0(P ≡ C;E ≡ D;W ≡ U)

∼
φC=

φP − φW

φE − φW

(3.23)

Senão (E ≡ C;P ≡ D;EE ≡ U)

∼
φC=

φE − φEE

φP − φEE

(3.24)

2. Utilize
∼
φC na fórmula da tabela 3.1 para calcular o fator �Downwind�.

3. Aplicando a equação (3.21), obtém-se:

φe = DWFeφD + (1−DWFe)φC (3.25)

calcula-se ξ conforme algoritmo 3.2:

Algoritmo 3.2 (Cálculo de ξ):
Se ue ≥ 0(P ≡ C;E ≡ D)

φf = DWFeφE + (1−DWFe)φP (3.26)

ξe = DWFe (3.27)

Senão (E ≡ C;P ≡ D)

φf = DWFeφP + (1−DWFe)φE (3.28)

ξe = 1−DWFe (3.29)

O valor de φ na face leste, por exemplo, é escrito como:

φe = ξeφE+
−
ξe φP (3.30)

onde
−
ξe= 1− ξe.



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 34

3.2.5 Equação de transporte escalar

Equação de transporte discreta

Nas seções anteriores, fórmulas para discretização dos termos de convecção-difusão fo-
ram estudadas. Nesta seção aplicando as fórmulas discutidas obtém-se uma equação algébrica
(discretizada) para um equação de transporte para uma propriedade escalar φ:

∂

∂t
(εmρmφ) +

∂

∂xi

(εmρmvmφ) =
∂

∂xi

(
Γφ

∂φ

∂xi

)
+Rφ (3.31)

A equação (3.31) é a equação de transporte para um propriedade escalar φ da fase gás
m = g ou fase sólida m = s. Os termos na equação (3.31) representam: a taxa de variação da
quantidade φ, os efeitos de convecção e difusão, e um termo fonte Rφ para φ. Em escoamento
multifásico o termo fonte contabiliza a transferência de massa e de quantidade de movimento
na interface.

Integrando a equação (3.31) sobre um volume de controle (�gura 3.1) e escrevendo
termo a termo, da esquerda para a direita, tem-se:

• Termo Transiente
∫

V C

∂

∂t
(εmρmφ) dV ≈

[
(εmρmφ)P − (εmρmφ)old

P

] ∆V

∆t
(3.32)

• Termo Convectivo
∫

∂

∂xi

(εmρmvmφ) dV ≈
{
ξe (εmρmφ)E +

−
ξe (εmρmφ)P

}
(um)eAe

−
{
ξw (εmρmφ)P +

−
ξw (εmρmφ)W

}
(um)w Aw

+

{
ξn (εmρmφ)N +

−
ξn (εmρmφ)P

}
(vm)nAn (3.33)

−
{
ξs (εmρmφ)P +

−
ξs (εmρmφ)S

}
(vm)sAs

+

{
ξt (εmρmφ)T +

−
ξt (εmρmφ)P

}
(wm)tAt

−
{
ξb (εmρmφ)P +

−
ξb (εmρmφ)B

}
(wm)bAb

• Termo Difusivo
∫

∂

∂xi

(
Γφ

∂φ

∂xi

)
dV ≈

(
Γφ
∂φ

∂x

)

e

Ae −
(

Γφ
∂φ

∂x

)

w

Aw

+

(
Γφ
∂φ

∂y

)

n

An −
(

Γφ
∂φ

∂y

)

s

As (3.34)

+

(
Γφ
∂φ

∂z

)

t

At −
(

Γφ
∂φ

∂z

)

b

Ab



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 35

Os �uxos difusivos são aproximados utilizando a equação (3.6). Por exemplo, o �uxo
difusivo da face leste será dado por:

(
Γφ
∂φ

∂x

)

e

≈ (Γφ)e

φE − φP

∆xe

(3.35)

Os valores do coe�ciente de difusão são calculados de acordo com a equação (3.7). Assim:

(Γφ)e =

[
1− fe

(Γφ)P

+
fe

(Γφ)E

]−1

=
(Γφ)P (Γφ)E

fe (Γφ)P + (1− fe) (Γφ)E

onde fe =
∆xE

∆xP + ∆xE

(3.36)

• Termo Fonte
O termo fonte é geralmente não linear e é primeiramente linearizado como segue:

Rφ ≈ Rφ −R
′
φφP (3.37)

Para a estabilidade do esquema numérico, é essencial que R′
φ ≥ 0 (Patankar (1980)).

Então a integral do termo fonte sobre um volume de controle resulta em:
∫

V C

RφdV ≈ Rφ∆V −R
′
φφP ∆V (3.38)

Combinando as equações deduzidas acima, obtém-se:
((

ρ
′
mφ

)
P
− (

ρ
′
mφ

)0

∆t

)
∆V

+

{
ξe (εmρmφ)E +

−
ξe (εmρmφ)P

}
(um)eAe −

{
ξw (εmρmφ)P +

−
ξw (εmρmφ)W

}
(um)w Aw

+

{
ξn (εmρmφ)N +

−
ξn (εmρmφ)P

}
(vm)nAn −

{
ξs (εmρmφ)P +

−
ξs (εmρmφ)S

}
(vm)sAs

+

{
ξt (εmρmφ)T +

−
ξt (εmρmφ)P

}
(wm)tAt −

{
ξb (εmρmφ)P +

−
ξb (εmρmφ)B

}
(wm)bAb

=

{
(Γφ)e

(φE − φP )

∆xe

Ae − (Γφ)w

(φP − φW )

∆xw

Aw

}

+

{
(Γφ)n

(φN − φP )

∆yn

An − (Γφ)s

(φP − φS)

∆ys

AS

}

+

{
(Γφ)t

(φT − φP )

∆zt

At − (Γφ)b

(φP − φB)

∆zb

Ab

}
+

(
Rφ −R

′
φφP

)
∆V (3.39)

Onde se de�ne as densidades macroscópicas como: ρ′m = εmρm.
A equação (3.39) pode ser rearranjada de forma a obter o seguinte sistema de equações

lineares para φ:
aPφP =

∑

nb

anbφnb + SP (3.40)



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 36

onde SP é a contribuição dos termos transientes e a integral do termo fonte linearizado, o
subscrito nb representa a contribuição das faces E,W , N , S, T e B da célula. Antes de utilizar
a equação acima para determinação de φ, é recomendado que a equação da continuidade
discretizada multiplicada por φ seja subtraída de (3.40). A razão para tal manipulação é
discutida em detalhes por Patankar (1980).

A parte homogênea da equação diferencial parcial para φ (Eq.3.31) admite in�nitas
soluções da forma (φ+ c), onde c é uma constante arbitrária. A equação de diferenças �nitas
para φ deve admitir o mesmo número de soluções. Do contrário, uma pequena �utuação no
balanço de massa durante as iterações pode produzir grandes �utuações nos valores de φ,
afetando de maneira adversa a convergência. Pode-se mostrar que as equações de diferenças
�nitas terão a propriedade desejada se:

aP =
∑

nb

anb (3.41)

quando as constribuições das instabilidades e do termo fonte à aP são descartadas. Patankar
(1980) denomina esta restrição de Regra 4. Uma equação da forma (3.40) deduzida da equação
(3.39) não satisfaz a Regra 4.

A forma discretizada da equação da continuidade pode ser obtida da equação (3.39),
assumindo φ = 1 e substituindo o termo fonte por

∑
l Rl. Então, subtraindo a equação da

continuidade discretizada multiplicada por φ da equação (3.39), obtém-se uma equação linear
da forma (3.40), na qual os coe�cientes são de�nidos como segue:

aE = De − ξe (εmρm)E (um)eAe (3.42)

aW = Dw−
−
ξw (εmρm)W (um)w Aw (3.43)

aN = Dn − ξn (εmρm)N (vm)nAn (3.44)

aS = Ds−
−
ξs (εmρm)S (vm)sAs (3.45)

aT = Dt − ξt (εmρm)T (wm)tAt (3.46)

aB = Db−
−
ξb (εmρm)B (wm)bAb (3.47)

aP =
∑

nb

anb + a0
P +R

′
φ∆V + [[

∑
Rl]]∆V (3.48)

SP = a0
Pφ

0
P +

−
Rφ ∆V + φP [[−

∑
Rl]]∆V (3.49)

a0
P =

(εmρm)0

∆t
∆V (3.50)

Dk =
(Γφ)k Ak

∆xk

, k = e, w, n, s, t, b (3.51)



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 37

Observações

Diferentemente do escoamento monofásico, a equação da continuidade para escoa-
mento multifásico apresenta um termo fonte (

∑
l Rl) que contabiliza a transferência de massa

na interface. Desde que a equação da continuidade multiplicada por φ é subtraída da equação
(3.40), o termo aparece na equação de transporte discretizada. A inclusão deste termo fonte
diminuiria a velocidade de convergência, e a sua inclusão no coe�ciente central desestabili-
zaria as iterações quando

∑
l Rl < 0. Portanto, o termo é manipulado de forma que, sua

contribuição para aP é não negativa. Utilizando a de�nição da função de duplo colchete:

[[R]] =





0, R ≤ 0

R, R > 0
(3.52)

Da de�nição acima, segue que:
R = [[R]]− [[−R]] (3.53)

Então o termo de transferência de massa na interface pode ser escrito como:

φP

∑

l

Rl = φP

[[∑

l

Rl

]]
− φP

[[
−

∑

l

Rl

]]
(3.54)

O primeiro termo do lado direito da equção acima contribui para o termo fonte e o segundo
contribui para o coe�ciente central.

Se uma discretização em lei de potências é procurada, toma-se o fator de convecção
como zero, isto é, ξ = 0 e

−
ξ= 1 e substitui-se D por: DA(|P |), onde P é o número de

Péclet para a face. Por exemplo, substituindo-se De na equação(3.42) por DeA (|Pe|), onde
A (|Pe|) =

[[
(1− 0, 1|Pe|)5]].

Há dois pontos a serem notados a respeito do uso de esquemas de discretização de
ordem elevada. Nos esquemas com segunda ordem ou superior os fatores ξ na equação da
continuidade discretizada podem diferir dos correspondentes fatores na equação de transporte
da propriedade φ. Então a equação discretizada para φ obtida pela subtração de φmultiplicada
pela equação da continuidade não deverá satisfazer a Regra 4. Portanto, se admite a hipótese
que os fatores de convecção na equação da continuidade discretizada são os mesmos daqueles
na equação de transporte da propriedade φ, e assim satisfazer a Regra 4.

O uso de métodos de ordem elevada pode resultar na violação, em algumas regiões da
Regra 2 de Patankar (1980). A Regra 2 estabelece que todos os coe�cientes anb na equação
(3.40) devem possuir o mesmo sinal, dito positivo. As bases físicas para esta regra é que o
aumento (ou diminuição) no valor de φ em uma célula vizinha deveria causar um aumento
(ou diminuição) nos valores de φP , não o contrário. Esta regra quando combinada com a



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 38

Regra 4 também garante que a discretização produz um sistema de equações diagonalmente
dominante. A regra é estritamente satisfeita pela equação acima somente quando utilizamos
o esquema FOUP. Quando esquemas de ordem superior são utilizados, alguns coe�cientes
podem tornar-se negativos quando o comportamento local de φ é monótona. Tais violações
da Regra 2 não são motivos para inquietações, desde que o limitador utiliza considerações
mais elaboradas para garantir que a solução é limitada e �sicamente realística.

3.3 Solução numérica do modelo hidrodinâmico para es-
coamento gás-sólido

Uma extensão do algoritmo SIMPLE de Patankar (1980) é utilizada para a resolução
das equações discretizadas. Vários pontos importantes precisam ser discutidos quando este
algoritmo, desenvolvido para escoamentos monofásicos, é estendido para resolver as equações
de um escoamento multifásico. Spalding (1980) lista três pontos a seguir e os classi�ca como:
�o primeiro é óbvio, o segundo um tanto menos, e o terceiro pode facilmente passar sem ser
notado�.

(i) Há mais campos de variáveis e, portanto mais equações comparadas com escoamentos
monofásicos. Isto desacelera a computação, mas não torna o algoritmo mais complexo.

(ii) A pressão aparece nas componentes da equação de quantidade de movimento para es-
coamento monofásico, mas não há uma equação conveniente para resolver o campo de
pressão. O fundamental do algoritmo SIMPLE é a dedução de uma equação para a pres-
são - a equação de correção para a pressão. A correção da pressão propicia a correção
para o campo da velocidade tal que a equação da continuidade é satisfeita exatamente
(para a precisão da máquina). Não há uma única maneira para a dedução de tal equação
para escoamento multifásico, desde que há mais uma equação para a continuidade no
escoamento multifásico.

(iii) As equações de quantidade de movimento para o escoamento multifásico são fortemente
acopladas por meio do termo de intercâmbio de quantidade de movimento. Fazer este
termo totalemte implícito é essencial para o sucesso do esquema numérico. Esta é a
principal idéia no �Implicit Multi�eld Field� (IMF), técnica de Harlow e Amsden (1975),
que está implementada no K-FIX �Kachina-Fully Implicit Exchange�. No MFIX, as
equações de quantidade de movimento são resolvidas em todo o domínio computacional.



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 39

Para tornar o termo de intercâmbio totalmente implícito todas as equações para cada
componente de velocidade (tanto da fase gás como da fase sólida) devem ser resolvidas
juntas, o que origina uma matriz com estrutura não padrão. Uma alternativa mais barata
é utilizar o algoritmo de eliminação parcial (PEA - �Partial Elimination Algorithm�) de
Spalding (1980), que será discutido adiante.

No �uxo granular multifásico outros dois pontos críticos determinam o sucesso do
esquema numérico. Um é o tratamento das regiões com altas concentrações de sólido (regiões
próximas de empacotamento). A fração de volume de sólido varia de zero até um valor má-
ximo em torno de 0, 6 na região do empacotamento. O limite inferior é facilmente manipulável
pela formulação das equações lineares tais que valores não-negativos da fração de volume são
calculados. A restrição da fração de volume de sólido no ou abaixo do valor máximo é mais
difícil. A formação da região de empacotamento é análoga a condensação do vapor compre-
ensível em um líquido incompressível. As forças de reação que resistem à compactação do
meio granular resultam na pressão do sólido, a qual deve ser distinguida da pressão do �uido.
Esta situação foi tratada nos modelos S3 de Pritchett, Blake e Garg (1978) e pelo �Illinois
Institute of Technology� (ITT) em Gidaspow e Ettehadieh (1983) através da introdução de
uma equação de estado que relaciona a pressão do sólido com a fração de volume do sólido. A
função da pressão do sólido aumenta exponencialmente quando a fração volumétrica de sólido
aproxima-se do limite do empacotamento, e assim retarda a compactação do sólido. Este mé-
todo permite o meio granular ser suavememte compressível. O meio granular pode também
ser considerado incompressível como foi feito por Syamlal e O'Brien (1988). É também o mé-
todo utilizado no código FLUENTTM(�Fluent Users Manual�, 1996). Neste método nenhuma
equação de estado para a pressão do sólido é necessária. A fração volumétrica de sólido no
empacotamento máximo precisa ser especi�cada, que também é um parâmetro explícito ou
implícito na equação de estado utilizada para o caso de suave compressibilidade.

Uma equação semelhante à equação para a correção da pressão do �uido pode ser
desenvolvida para a pressão do sólido. Tal equação é resolvida no código FLUENTTM. O
MFIX, ao invés disso, utiliza uma equação para a correção da fração volumétrica do sólido. A
equação para a correção de pressão do sólido requer que ∂Ps

∂εs
não desapareça quando εs → 0.

A equação para a correção da fração volumétrica de sólido não possui tal restrição, mas deve
contabilizar o efeito da pressão do sólido tal que a computação seja estabilizada na região de
empacotamento.

Uma segunda questão importante é a di�culdade no cálculo dos campos das variáveis
nas interfaces nas quais uma fração volumétrica tende a zero. O campo de variáveis associado a



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 40

uma fase não está de�nido nas regiões onde a fração volumétrica é zero, e assim podem admitir
valores arbitrários. Entretanto, um algoritmo computacional não deve utilizar arbitrariedade
no conjunto de valores. Como exemplo, em seção anterior mostrou-se que o uso de uma
média harmônica para os cálculos dos valores dos coe�ciente de difusão da face de uma célula
previnirá a difusão de φ nas regiões onde a fase associada com ela é ausente. O cálculo das
componentes das velocidades em tais interfaces é mais difícil do que as quantidades escalares
devido a linearização dos termos convectivos não lineares. Em uma interface onde a fração
volumétrica da fase está próxima de zero,a componente normal da velocidade torna-se muito
grande. Desde que o produto da fração volumétrica da fase e da componente da velocidade é
ainda próxima de zero, o erro na conservação de quantidade de movimento é negligenciável.
Entretanto, valores elevados da velocidade da fase rapidamente desestabilizam os cálculos,
e o método é necessário para previnir tais desestabilizações. O MFIX utiliza um cálculo
aproximado da velocidade normal nas interfaces (de�nido por um pequeno valor inicial para
a fração volumétrica da fase).

Escoamentos gás-sólido são inerentemente instáveis. Cálculos no estado estável são
possíveis somente para poucos casos, tais como transporte pneumático (diluído) de sólidos.
Para a maioria dos escoamentos gás-sólido, uma simulação transiente é conduzida e os resul-
tados são médios-temporais. Simulações transientes divergem, se grande amplitude para o
passo temporal é admitida. Passos temporais com amplitude muito pequena torna a compu-
tação muito lenta. O MFIX ajusta automaticamente a amplitude do passo, dentro de limites
especi�cados pelo usuário, para reduzir o tempo de processamento.

3.3.1 Equação da conservação da quantidade de movimento

A discretização das equações de quantidade de movimento é semelhante àquela de
transporte escalar, exceto que os volumes de controles estão deslocados. De acordo com
Patankar (1980), se os componentes de velocidade e pressão são armazenados em um mesmo
local do grid, obtém-se uma solução aceitável. Um grid deslocado é usado para previnir
resultados não físicos para os campos de pressão. Conforme mostrado na �gura 3.7, em
relação ao volume de controle escalar, a direção x é deslocada para a metade leste da célula do
volume de controle. Similarmente as direções y e z são deslocadas para a metade das células.



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 41

w

i− 1/2

W

i

NW N

i+ 1/2

P

S

NE

i+ 1

E

i+ 3/2

e

Figura 3.7: Volume de controle para a componente x da equação de quantidade de movimento.

Equação da conservação da quantidade de movimento discreta

Para calcular a convecção na variação da quantidade de movimento, as componentes
da velocidade são requeridas nas posições E, W , N e S. Elas são calculadas a partir dos
valores vizinhos, por meio de uma média aritmética:

(um)E = fE (um)P + (1− fE) (um)e (3.55)

(vm)N = fP (vm)NW + (1− fP ) (vm)NE (3.56)

O valor da fração volumétrica requerida no cento da célula P é calculado de maneira
semelhante:

(εm)P = fP (εm)W + (1− fP ) (εm)E (3.57)

onde:
fE =

∆xe

∆xP + ∆xe

e fP =
∆xE

∆xW + ∆xE

(3.58)

Agora a componente x da equação de quantidade de movimento pode ser escrita
como:

aP (um)P =
∑

nb

anb (um)nb + bP − AP (εm)P

(
(Pg)E − (Pg)W

)

+

(∑

l

Flm (ul − um)P

)
∆Ve (3.59)

A equação (3.59) é semelhante à equação de transporte escalar discretizada, exceto
pelos dois últimos termos: o termo do gradiente de pressão é determinado baseado no campo
de pressão atual Pg e é adicionado ao termo fonte da equação linear. O termo de transferência
na interface acopla todas as equações para a mesma componente (x, y ou z).



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 42

A de�nição dos termos restantes na equação (3.59) é dada nas relações abaixo:

ae = DE − ξE (εmρm)e (um)E AE

aw = DW−
−
ξW (εmρm)w (um)W AW

an = DN − ξN (εmρm)n (um)N AN

as = DS−
−
ξS (εmρm)s (um)S AS

at = DT − ξT (εmρm)t (um)T AT

ab = DB−
−
ξB (εmρm)b (um)B AB (3.60)

aP =
∑

nb

anb + a0
P +R

′
um

∆Ve +
[[∑

Rl

]]
∆Ve + S

′

b = a0
Pu

0
m+

−
Rum ∆Ve + um

[[
−

∑
Rl

]]
∆Ve + (εmρm)e gx∆Ve+

−
S

a0
P =

(εmρm)0

∆t
∆Ve

DE =
(um)E AE

∆xE

O coe�ciente central aP e o termo fonte b contêm os termos adicionais S ′ e
−
S, os quais

contabilizam as fontes que surgem das coordenadas cilíndricas, do modelo de meio poroso e
do termo da tensão de cisalhamento.

Condições de contorno

Nas simulações realizadas no presente trabalho, utilizou-se para a fase gasosa a hipó-
tese de plena aderência na parede. Tal hipótese implica na condição de contorno denominada
condição de não deslizamento na parede. Esta estabelece que as componentes normal e tan-
gencial do vetor velocidade em uma parede estacionária são nulas, isto é, ug,W = vg,W = 0. No
caso da fase sólida, o estabelecimento das condições de contorno na parede é mais complicado
devido à di�culdade da caracterização da interação das partículas com a parede. As condições
de contorno na parede, implementadas no MFIX para a solução da equação linear são dadas a
seguir (Syamlal (1998)). A componente da velocidade na direção y na parede leste é utilizada
como um exemplo, �gura 3.8. A implementação para as outras componentes é análoga.

1. Deslizamento livre na parede

vm (i, j + 1/2, k)− vm (i− 1, j + 1/2, k) = 0 (3.61)

2. Não deslizamento na parede

vm (i, j + 1/2, k) + vm (i− 1, j + 1/2, k) = 0 (3.62)



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 43

Uma outra condição de contorno devido ao deslizamento das partículas na parede é
apresentada a seguir:

3. Deslizamento parcial na parede

∂vm

∂n
+ hv (vm − vw) = 0 (3.63)

onde ∂
∂n

denota a diferenciação ao longo da direção normal.

i− 1 i i− 1 i

Deslizamento livre Não deslizamento

Figura 3.8: Condição de Deslizamento livre e de Não Deslizamento na parede leste.

A forma discreta da equação (3.63), por exemplo, na face leste é:

vm (i, j + 1/2, k)

(
hv

2
+

1

∆xE

)
+ vm (i− 1, j + 1/2, k)

(
hv

2
− 1

∆xE

)
= hvvw (3.64)

A equação (3.64) é uma condição de deslizamento generalizada. Os coe�cientes hv e
vw podem descrever a condição de não deslizamento quando hv →∞ e vw = 0; deslizamento
livre quando hv = 0 e uma condição de deslizamento parcial quando hv 6= 0 e vw 6= 0.

Para as simulações realizadas neste trabalho, no caso da fase sólida, considerou-se a
hipótese de não deslizamento na parede.

Algoritmo para equação linear

Algoritmo 3.3 :
As equações lineares para resolução das equações de quantidade de movimento que foram
modi�cadas seguem:

1. Cálculo da velocidade média da quantidade de movimento das faces das células. (eq.(3.55
e 3.56)).

2. Cálculo do coe�ciente de convecção ξ. (alg.3.2).



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 44

3. Cálculo do coe�ciente anb. (eq.(3.60)).

4. Modi�car os coe�cientes vizinhos para ter em conta a presença de superfícies internas
(assume-se livre deslizamento nas paredes) (eq.(3.61)).

5. Cálculo do coe�ciente central e termos fontes. Para paredes impermeáveis superfícies
internas todos os coe�cientes vizinhos e a componente normal da velocidade são zero.

3.3.2 Eliminação parcial do acoplamento na interface

A presença dos termos de transferência na interface é uma característica que distingue
as equações do escoamento multifásico das do monofásico. Usualmente, os termos de transfe-
rência na interface acoplam fortemente as componentes da velocidade e temperatura de cada
fase às correspondentes variáveis da outra fase. O desacoplamento das equações por meio do
cálculo dos termos de transferência na interface a partir dos valores da iteração anterior tornará
as iterações instáveis ou forçará passos temporais pequenos. No outro extremo, a resolução
de todas as equações discretas para uma certa componente gerará uma matriz não padrão de
ordem elevada. Uma alternativa afetiva que mantém um alto grau de acoplamento entre as
equações enquanto fornece uma matriz septadiagonal padrão (3D) é o algoritmo parcial de
Spalding (1980). O algoritmo é ilustrado com a seguinte equação modelo:

εmρm
∂φ

∂t
+ εmρm (vm)i

∂φm

∂xi

=
∂

∂xi

(
Γφm

∂φm

∂xi

)
+Rφm − φm

(∑
Rlm

)
+

M∑

l=0

Flm (φl − φm)

(3.65)
(Note que Flm = Fml e Fmm = 0)

A equação correspondente discreta é:

(am)P (φm)P =
∑

nb

(am)nb (φm)nb + bm + ∆V
M∑

l=0

Flm [(φl)P − (φm)P ] (3.66)

que é semelhante em forma à equação discreta da quantidade de movimento, discutida na
seção 3.3.1.

Primeiramente, discute-se o problema com um desacoplamento direto das equações.
Por exemplo, considere o caso do �uxo de duas fases (M = 1):

(a0)P (φ0) =
∑

nb

(a0)nb (φ0)nb + b0 + ∆V F10 [(φ1)P − (φ0)P ] (3.67)

(a1)P (φ1)P =
∑

nb

(a1)nb (φ1)nb + b1 + ∆vF10 [(φ0)P − (φ1)P ] (3.68)



CAPÍTULO 3. METODOLOGIA DE SOLUÇÃO NUMÉRICA 45

quando F10 → 0 duas equações estão desacopladas e a solução para (φ0)P , por exemplo, é:

(φ0)P =

∑
nb (a0)nb (φ0)nb + b0

(a0)P

(3.69)

quando F10 →∞ as equações estão fortemente acopladas e as soluções são:

(φ0)P = (φ1)P =

∑
nb (a0)nb (φ0)nb + b0 +

∑
nb (a1)nb (φ1)nb + b1

(a0)P + (a1)P

(3.70)

Um esquema iterativo que trata o termo de transferência na interface meramente
como um termo fonte fornecerá a solução correta para valores pequenos de F10, mas falhará
no caso e