UNIVERSIDADE DE BRASÍLIA
FACULDADE DE TECNOLOGIA
DEPARTAMENTO DE ENGENHARIA CIVIL E AMBIENTAL
MODELO ADVECTIVO-DISPERSIVO DE TRANSPORTE DE
SOLUTOS EM SOLO NÃO-SATURADO UTILIZANDO OS
MÉTODOS DAS CARACTERÍSTICAS E DOS ELEMENTOS
FINITOS
RUBENS MACIEL WANDERLEY
ORIENTADOR: SERGIO KOIDE
DISSERTAÇÃO DE MESTRADO EM TECNOLOGIA
AMBIENTAL E RECURSOS HÍDRICOS
PUBLICAÇÃO: MTARH.DM-021A/2000
BRASÍLIA / DF: FEVEREIRO – 2000
UNIVERSIDADE DE BRASÍLIA
FACULDADE DE TECNOLOGIA
DEPARTAMENTO DE ENGENHARIA CIVIL E AMBIENTAL
MODELO ADVECTIVO-DISPERSIVO DE TRANSPORTE DE
SOLUTOS EM SOLO NÃO-SATURADO UTILIZANDO OS
MÉTODOS DAS CARACTERÍSTICAS E DOS ELEMENTOS
FINITOS
RUBENS MACIEL WANDERLEY
DISSERTAÇÃO DE MESTRADO SUBMETIDA AO DEPARTAMENTO DE
ENGENHARIA CIVIL E AMBIENTAL DA FACULDADE DE TECNOLOGIA DA
UNIVERSIDADE DE BRASÍLIA, COMO PARTE DOS REQUISITOS
NECESSÁRIOS PARA A OBTENÇÃO DO GRAU DE MESTRE.
APROVADA POR:
________________________________
SERGIO KOIDE, PhD (ENC-UnB)
(ORIENTADOR)
________________________________
NÉSTOR ALDO CAMPANA, DSc (ENC-UnB)
(EXAMINADOR INTERNO)
________________________________
JOSÉ ALMIR CIRILO, DSc (UFPE)
(EXAMINADOR EXTERNO)
DATA: BRASÍLIA/DF, 28 DE FEVEREIRO DE 2000
ii
FICHA CATALOGRÁFICA
WANDERLEY, RUBENS MACIEL
Modelo Advectivo-Dispersivo de Transporte de Solutos em Solo Não-Saturado
Utilizando os Métodos das Características e dos Elementos Finitos
[Distrito Federal] 2000.
xvi, 126 p., 297mm (ENC/FT/UnB, Mestre, Tecnologia Ambiental e Recursos
Hídricos, 2000).
Dissertação de Mestrado – Universidade de Brasília. Faculdade de Tecnologia.
Departamento de Engenharia Civil e Ambiental.
1. Transporte de solutos
2. Solo não-saturado
3. Método das características
4. Método dos elementos finitos
I. ENC/FT/UnB
II. Título (série)
REFERÊNCIA BIBLIOGRÁFICA
WANDERLEY, R.M. (2000). Modelo Advectivo-Dispersivo de Transporte de Solutos em
Solo Não-Saturado Utilizando os Métodos das Características e dos Elementos Finitos.
Dissertação de Mestrado, Publicação MTARH.DM-021A/2000, Departamento de Engenharia
Civil, Universidade de Brasília, Brasília, DF, 126 p.
CESSÃO DE DIREITOS
NOME DO AUTOR: Rubens Maciel Wanderley
TÍTULO DA DISSERTAÇÃO DE MESTRADO: Modelo Advectivo-Dispersivo de
Transporte de Solutos em Solo Não-Saturado Utilizando os Métodos das Características e dos
Elementos Finitos.
GRAU / ANO: Mestre / 2000
É concedida à Universidade de Brasília permissão para reproduzir cópias desta dissertação de
mestrado e para emprestar ou vender tais cópias somente para propósitos acadêmicos e
científicos. O autor reserva outros direitos de publicação e nenhuma parte desta dissertação de
mestrado pode ser reproduzida sem a autorização por escrito do autor.
________________________________
Rubens Maciel Wanderley
SHCES 609, Bloco C, Apto. 104, Cruzeiro Novo
CEP 70655-693 – Brasília/DF – Brasil
iii
DEDICATÓRIA
Ao meu pai, Ezequiel Luis Vanderlei, e à minha mãe, Eulina Maciel Wanderley, que, com coragem,
esperança, sabedoria, simplicidade, paciência e, principalmente, amor, souberam conduzir
corretamente os meus passos, ensinando-me a viver e compartilhar a vida.
A todos os meus irmãos, Eder, Claudia, Claudio, Carmen, Lany e Leda, que foram sempre um
complemento em minha vida.
À minha namorada, Marcia, que foi sempre amorosa, paciente e compreensiva em todos os momentos
que estivemos juntos.
iv
AGRADECIMENTOS
Agradeço primeiramente a Deus, por estar sempre presente em minha vida.
Ao professor Sergio Koide que, por mais de cinco anos, tem sido mais que um orientador, tem
sido um amigo paciente e dedicado, acreditando em meu potencial e dando-me o apoio
necessário em todo esse tempo.
Aos professores do MTARH: Cristina Celia Silveira Brandão, Marco Antônio de Souza,
Ricardo Silveira Bernardes, Oscar de Moraes Cordeiro Netto, Néstor Aldo Campana e Nabil
Joseph Eid, pelos conhecimentos transmitidos e pelo incentivo durante o período do curso.
Em especial, à professora Cristina Celia Silveira Brandão, pela amizade e pelo constante e
enorme incentivo.
A todos os amigos do mestrado, em especial, aos amigos Waldemir, Jamaci, Luzideth,
Patrícia, Antônio Carlos e Carramaschi, da turma de 97, por terem trilhado junto comigo esse
árduo caminho, sempre com apoio, respeito e grande amizade.
Aos funcionários “Boy”, André e João, pela alegria.
À CAPES pelo apoio financeiro e à UnB pela oportunidade oferecida.
Aos meus pais, Ezequiel Luiz Vanderlei e Eulina Maciel Wanderley, meus irmãos, Eder,
Claudia, Claudio, Carmen, Lany e Leda, meus sobrinhos, Jacqueline, Raphael, Thayná,
Lucas, Nayara, Marcos, Matheus, Milena e Larissa, e todos os outros familiares, pelo apoio,
pelo respeito, pela amizade, pelo carinho e pelo amor que nunca deixaram faltar em minha
vida.
À minha namorada, Marcia, pelo amor, pela paciência e pelo enorme apoio que me deu
durante todo o transcorrer de nosso tempo de união.
A todos aqueles que não foram citados, mas que fizeram parte deste trabalho.
Muito Obrigado a Todos!
v
RESUMO
Este trabalho apresenta o desenvolvimento de um modelo numérico bidimensional
para a simulação do transporte de solutos em solo não-saturado, com consideração dos
processos de produção e decaimento e de sorção de solutos. O modelo foi implementado em
linguagem FORTRAN 77, a partir de um modelo de simulação de fluxo não-saturado,
desenvolvido por Koide (1990) e modificado por Campos (1998).
A equação diferencial que rege o transporte de solutos é resolvida pelo modelo com a
combinação de dois métodos numéricos: o método das características e o método dos
elementos finitos. O processo geral de transporte é numericamente dividido em dois
processos: a advecção e a dispersão. Os processos de produção e decaimento e de sorção de
solutos são incluídos na dispersão. O método das características é utilizado para a simulação
da advecção enquanto a dispersão é simulada em seguida pelo método dos elementos finitos.
O modelo permite, para controle do modo de cálculo da advecção, a seleção do padrão
de distribuição de partículas pelo domínio e também do tipo de controle de reposicionamento
dessas partículas durante a simulação.
O modelo foi testado para cinco problemas apresentados na literatura, cujas soluções
analíticas são conhecidas. Nesses testes, foram verificados a precisão e o comportamento do
modelo em relação aos modos de cálculo e em relação ao regime de transporte, que pode ser
predominantemente advectivo ou dispersivo.
Os testes mostraram que o modelo tem, em geral, boa precisão, em especial, para os
casos de regime de fluxo permanente. Em relação aos modos de cálculo, não há grandes
diferenças entre os padrões de distribuição de partículas e um reposicionamento menos
freqüente é recomendável para melhorar a precisão do modelo, reduzindo-se o tempo de
processamento, podendo-se inclusive utilizar um menor número de partículas, diminuindo a
necessidade de uso de memória computacional.
vi
ABSTRACT
This work presents the development of a two-dimensional numerical model for the
simulation of solute transport in unsaturated soil, including the processes of production and
decay and of sorption of solute. The model was implemented in FORTRAN 77 structured
language, based on a model for the simulation of unsaturated flux, developed by Koide (1990)
and modified by Campos (1998).
The diferential equation of the solute transport is solved by the model with a
combination of two numerical methods: the method of characteristics and the finite element
method. The overall process of transport is divided in two processes: the advection and the
dispersion. The minor processes of production and decay and of sorption of solute are
included in the dispersion. The method of characteristics is used for the simulation of
advection while the finite element method is used in sequence for the simulation of the
dispersion and minor processes.
The model allows, for control of the mode of advection calculation, the selection of
the particle distribution patterns over the domain and also the type of reposicioning control of
these particles during simulation.
The model was tested for five problems presented in the literature, to wich analytical
solutions are known. In these tests, the precision and the behaviour of the model in relation to
the mode of advection calculation and the transport regim, that can be predominantly
advective or dispersive, were verified.
The tests showed that the model presented, in general, good precision, especially for
steady state flux. In terms of mode of advection calculation, there are no significative
differences between the particle distribution paterns and a less frequent reposicioning is
suggested to improve the precision and reduce processing time, allowing, in adiction, the use
of a smaller number of particles, reducing the computational memory requirement.
vii
ÍNDICE
1-
2-
INTRODUÇÃO ............................................................................................................1
1.1-
CONTEXTO GERAL ...................................................................................................1
1.2-
OBJETIVOS ...............................................................................................................2
1.3-
ESTRUTURA DA DISSERTAÇÃO..................................................................................2
FUNDAMENTAÇÃO TEÓRICA ................................................................................4
2.1-
DISTRIBUIÇÃO DA ÁGUA NA SUBSUPERFÍCIE .............................................................4
2.2-
O SOLO E SUAS CARACTERÍSTICAS ...........................................................................6
2.2.1- Características dos Solos Relacionadas à Água ....................................................8
2.2.1.1- Porosidade ..................................................................................................8
2.2.1.2- Fração Volumétrica de Água.......................................................................9
2.2.1.3- Ponto de Murchamento e Capacidade de Campo ........................................9
2.2.1.4- Condutividade Hidráulica ......................................................................... 10
2.3-
DESCRIÇÃO DO FLUXO EM SOLO NÃO-SATURADO .................................................. 10
2.3.1- Interações Entre a Água e o Solo ....................................................................... 10
2.3.2- Estado Energético da Água do Solo ................................................................... 12
2.3.3- Curva Característica do Solo: a Relação
.................................................... 14
2.3.3.1- Efeitos de Histerese ................................................................................... 16
2.3.4- Equacionamento do Fluxo Não-Saturado ........................................................... 18
2.3.5- A Condutividade Hidráulica do Solo Não-Saturado ...........................................20
2.4-
DESCRIÇÃO DO TRANSPORTE DE SOLUTOS EM SOLO NÃO-SATURADO ..................... 22
2.4.1- Descrição da Advecção ..................................................................................... 22
2.4.2- Descrição da Difusão......................................................................................... 23
2.4.3- Descrição da Dispersão ..................................................................................... 24
2.4.4- Dispersão Hidrodinâmica .................................................................................. 25
2.4.5- Descrição da Sorção .......................................................................................... 26
2.4.5.1- Modelos de Sorção em Equilíbrio .............................................................. 27
2.4.6- Equacionamento do Transporte de Solutos em Solo Não-Saturado .................... 31
3-
MODELAGEM ...........................................................................................................36
3.1-
MODELOS FÍSICOS ................................................................................................. 36
3.2-
MODELOS ANALÓGICOS ......................................................................................... 37
viii
3.3-
MODELOS MATEMÁTICOS ...................................................................................... 37
3.3.1- Modelos Determinísticos e Estocásticos ............................................................ 38
3.3.2- Modelos de Parâmetros Distribuídos e Agregados .............................................38
3.3.3- Modelos Conceituais e Empíricos ......................................................................39
3.3.4- Modelos Analíticos e Numéricos ....................................................................... 39
3.44-
IMPORTÂNCIA DOS MODELOS MATEMÁTICOS ......................................................... 40
A MODELAGEM MATEMÁTICA ..........................................................................42
4.1-
O MÉTODO DAS DIFERENÇAS FINITAS .................................................................... 44
4.2-
O MÉTODO DOS ELEMENTOS FINITOS ..................................................................... 48
4.2.1- Método dos Resíduos Ponderados ..................................................................... 51
4.2.2- O Método de Galerkin e as Funções de Base ..................................................... 52
4.3-
COMPARAÇÃO ENTRE OS MÉTODOS DAS DIFERENÇAS FINITAS E DOS ELEMENTOS
FINITOS ................................................................................................................. 53
4.4-
OUTROS MÉTODOS APLICÁVEIS À SIMULAÇÃO DO TRANSPORTE DE SOLUTOS ......... 54
4.4.1- O Método das Características ............................................................................ 54
4.4.2- O Método de Caminhamento Aleatório (Random Walk).................................... 56
5-
ESTUDOS RECENTES EM MODELAGEM NUMÉRICA DO TRANSPORTE DE
CONTAMINANTES EM SOLOS ............................................................................ 57
5.1-
6-
A ESCOLHA DO MÉTODO A SER EMPREGADO .......................................................... 62
FORMULAÇÃO DO MODELO................................................................................ 64
6.1-
APLICAÇÃO DO MÉTODO DAS CARACTERÍSTICAS .................................................... 64
6.1.1- Geração e Distribuição das Partículas Sobre o Domínio .................................... 65
6.1.2- Movimentação e Acompanhamento de Partículas ..............................................67
6.2-
APLICAÇÃO DO MÉTODO DOS ELEMENTOS FINITOS ................................................. 68
6.2.1- Funções de Base ................................................................................................ 72
7-
6.3-
INTEGRAÇÃO NO TEMPO......................................................................................... 73
6.4-
DETERMINAÇÃO DAS VELOCIDADES DE ESCOAMENTO ............................................ 75
6.5-
DETERMINAÇÃO DO COEFICIENTE DE DISPERSÃO.................................................... 78
6.6-
CRITÉRIO DE ESTABILIDADE ................................................................................... 80
6.7-
EVENTOS PARTICULARES ....................................................................................... 81
ESTRUTURA DO PROGRAMA E TESTES DO MODELO .................................. 83
ix
7.1-
ESTRUTURA DO PROGRAMA SSTRANS .................................................................. 83
7.1.1- Entrada de Dados .............................................................................................. 84
7.1.2- Montagem de Matrizes e Vetores Auxiliares ..................................................... 85
7.1.3- Determinação do Menor Comprimento Representativo do Domínio .................. 86
7.1.4- Determinação da Velocidade de Escoamento ..................................................... 86
7.1.5- Geração e Distribuição de Partículas ................................................................. 87
7.1.6- Movimentação das Partículas – Advecção ......................................................... 87
7.1.7- Cálculo da Dispersão ......................................................................................... 89
7.2-
TESTES DE VERIFICAÇÃO DO MODELO .................................................................... 89
7.2.1- Teste 1 – Injeção Contínua de Fluido com Soluto no Topo de uma Coluna
Vertical de Solo ............................................................................................... 90
7.2.2- Teste 2 – Concentração Constante no Topo de uma Coluna Vertical de Solo ..... 93
7.2.3- Teste 3 – Transporte de Soluto em uma Coluna Vertical de Solo com Fluxo
Transiente ........................................................................................................ 94
7.2.4- Teste 4 – Concentração Constante no Topo de um Perfil Vertical de Solo ......... 96
7.2.5- Teste 5 – Transporte com Produção/Decaimento e Sorção de Soluto em uma
Coluna Vertical de Solo ................................................................................... 98
8-
9-
RESULTADOS ......................................................................................................... 101
8.1-
RESULTADOS DO TESTE 1 ..................................................................................... 101
8.2-
RESULTADOS DO TESTE 2 ..................................................................................... 106
8.3-
RESULTADOS DO TESTE 3 ..................................................................................... 109
8.4-
RESULTADOS DO TESTE 4 ..................................................................................... 113
8.5-
RESULTADOS DO TESTE 5 ..................................................................................... 117
8.6-
COMENTÁRIOS FINAIS .......................................................................................... 119
CONCLUSÕES E RECOMENDAÇÕES ................................................................ 121
9.1-
CONCLUSÕES ....................................................................................................... 121
9.2-
RECOMENDAÇÕES PARA PESQUISAS FUTURAS ...................................................... 122
REFERÊNCIAS BIBLIOGRÁFICAS............................................................................. 124
x
ÍNDICE DE FIGURAS
Figura 2.1- Seção vertical de um terreno, mostrando as regiões da subsuperfície ....................5
Figura 2.2- Possíveis situações que um solo não-saturado pode apresentar ..............................7
Figura 2.3- Efeito das forças capilares e de adsorção ............................................................ 12
Figura 2.4- Perfil do potencial total da água do solo ( ) com a posição (z), mostrando os
potenciais gravitacional (
g)
e de pressão (
p)................................................... 14
Figura 2.5- Curva característica de um solo genérico ............................................................ 15
Figura 2.6- Histerese na relação
, mostrando as curvas intermediárias ........................... 17
Figura 2.7- Volume elementar representativo do solo, mostrando o balanço de massa da água
do solo ...............................................................................................................19
Figura 2.8- Efeito da presença de ar sobre a condutividade de um poro ................................. 21
Figura 2.9- Volume elementar representativo do solo, mostrando o balanço de massa do
soluto ................................................................................................................. 32
Figura 4.1- Malhas de elementos finitos ilustrando os tipos de grade .................................... 46
Figura 4.2- Grade de malha centrada ilustrando a numeração dos nós ................................... 47
Figura 4.3- Grade de elementos finitos com elementos triangulares ......................................49
Figura 4.4- Plano descrito por uma função de base linear para um elemento triangular ......... 50
Figura 4.5- Comportamento da função de base linear para um elemento triangular ilustrando o
plano formado .................................................................................................... 52
Figura 5.1- Efeitos da oscilação e dispersão numéricas em uma frente de passagem de
contaminante ......................................................................................................59
Figura 6.1- Padrões de distribuição de partículas .................................................................. 66
Figura 6.2- Elemento típico .................................................................................................. 72
Figura 6.3- Erro de aproximação na movimentação das partículas ........................................ 80
Figura 7.1- Estrutura do programa SSTRANS ......................................................................84
Figura 7.2- Fluxograma do subprograma MOVEM ............................................................... 88
Figura 7.3- Diagrama esquemático/hierárquico do programa SSTRANS .............................. 90
Figura 7.4- Malha de elementos finitos para o teste 1 (caso unidimensional) ......................... 92
Figura 7.5- Malha de elementos finitos para o teste 4 (caso bidimensional)........................... 98
Figura 8.1- Influência do padrão de geração e distribuição das partículas para o teste 1 ...... 102
Figura 8.2- Resultados para o teste 1, com alto valor para o número de Peclet .................... 103
Figura 8.3- Resultados para o teste 1, com baixo valor para o número de Peclet.................. 104
xi
Figura 8.4- Resultados para o teste 1, com controle de reposicionamento do tipo 1 ............. 105
Figura 8.5- Resultados para o teste 1, com controle de reposicionamento do tipo 2 ............. 105
Figura 8.6- Influência do padrão de geração e distribuição de partículas para o teste 2 ........ 106
Figura 8.7- Resultados para o teste 2, com alto valor para o número de Peclet .................... 107
Figura 8.8- Resultados para o teste 2, com baixo valor para o número de Peclet.................. 108
Figura 8.9- Resultados para o teste 2, com controle de reposicionamento do tipo 1 ............. 109
Figura 8.10- Resultados para o teste 2, com controle de reposicionamento do tipo 2 ........... 109
Figura 8.11- Resultados para o fluxo, no teste3 ................................................................... 110
Figura 8.12- Resultados para o teste 3, mostrando a influência dos padrões de geração e
distribuição de partículas, com controle de reposicionamento do tipo 1 ............ 111
Figura 8.13- Resultados para o teste 3, mostrando a influência dos padrões de geração e
distribuição de partículas, com controle de reposicionamento do tipo 2 ............ 113
Figura 8.14- Resultados para o teste 4, com controle de reposicionamento do tipo 1 ........... 115
Figura 8.15- Resultados para o teste 4, com controle de reposicionamento do tipo 2 ........... 116
Figura 8.16- Resultados para o teste 5, com controle de reposicionamento do tipo 1 ........... 118
Figura 8.17- Resultados para o teste 5, com controle de reposicionamento do tipo 2 ........... 119
xii
LISTA DE SÍMBOLOS
Ae ................................................. área do elemento ...................................................... [L2]
[Ai] ......... matriz de coeficientes do sistema de equações para o cálculo de vi
B ......................................................... largura ................................................................ [L]
{Bi} .......... vetor de constantes do sistema de equações para o cálculo de vi
Bd ........................................... massa específica do solo ............................................ [ML-3]
C ............................... concentração de massa do soluto na água ................................ [ML-3]
{C} ..................................... vetor de concentrações nodais
Ĉ .......................................... solução aproximada de C ............................................ [ML-3]
C* .............................. concentração específica de soluto sorvido ................... [adimensional]
C*max ................... máxima concentração específica de soluto sorvido ............. [adimensional]
C0 ............................................ concentração conhecida ............................................. [ML-3]
{Ci} .......... vetor de constantes do sistema de equações para o cálculo de vi
Cc .......................................... coeficiente de capacidade ............................................... [L-1]
D ........... coeficiente de dispersão hidrodinâmica ou coeficiente de dispersão ............ [L2T-1]
D̂ ............................................... aproximação de D ................................................ [L2T-1]
Dmn ........ componente do coeficiente de dispersão na direção mn (m, n = x, z) ........... [L2T-1]
D̂ mn .......................................... aproximação de Dmn ............................................... [L2T-1]
Dd ............................ coeficiente de dispersão mecânica do solo ............................... [L2T-1]
Dd,ij ... componente do coeficiente de dispersão mecânica na direção ij (i, j = x, y) ...... [L2T-1]
Dm ............................ coeficiente de difusão molecular do soluto ............................... [L2T-1]
Ds* ...................................... coeficiente de difusão do solo ........................................ [L2T-1]
DL .................................. coeficiente de dispersão longitudinal ................................... [L2T-1]
DT ................................... coeficiente de dispersão transversal .................................... [L2T-1]
F1, F2 .......................... variáveis auxiliares de solução analítica
[G] ............................................ matriz de condutância
H ........................................................... altura ................................................................. [L]

J ........................................ fluxo de massa total de soluto .................................... [ML-2T-1]
Ji ........................... fluxo de massa de soluto na direção i (i = x, y, z) ...................... [ML-2T-1]

J a .................................. fluxo de massa advectivo de soluto ................................ [ML-2T-1]

J d ................................... fluxo de massa difusivo de soluto ................................. [ML-2T-1]
xiii

J h ................. fluxo de massa hidrodinamicamente dispersivo de soluto ................ [ML-2T-1]

J m ..................... fluxo de massa mecanicamente dispersivo de soluto ................... [ML-2T-1]
K ............................................ condutividade hidráulica ............................................. [LT-1]
Ki .............. componente da condutividade hidráulica na direção i (i = x, z) ................ [LT-1]
K̂ i .............................................. aproximação de Ki .................................................. [LT-1]
Kd .......................................... coeficiente de distribuição .......................................... [L3M-1]
Kf .................................... constante de sorção de Freundlich ................................... [L3M-1]d
Kl ..................................... constante de sorção de Langmuir ..................................... [L3M-1]
Kp ............................................. coeficiente de partição ............................................. [L3M-1]
Ksat ............................ condutividade hidráulica do solo saturado ................................ [LT-1]
Lm .......................... menor comprimento representativo do domínio ................................. [L]
[M] ........................................ matriz de armazenamento
{N} ............................ vetor de produção/decaimento de massa
NN ............................................. número total de nós
O ............................................... operador diferencial
Pe ................................................. número de Peclet ...................................... [adimensional]
R ................................................ fator de retardação ..................................... [adimensional]
S .................................................. grau de saturação ...................................... [adimensional]
Vt ........................................ volume total da matriz sólida ............................................. [L3]
Vv ................................................. volume de vazios ...................................................... [L3]
Vw ............................................ volume da fase líquida .................................................. [L3]
W .................................. área do domínio (região do sistema)
ai, bi, ci .............................. coeficientes das funções de base
b .......................................................... largura ................................................................ [L]
d .............................................. expoente de Freundlich ................................. [adimensional]
dp ........................................... deslocamento da partícula ................................................. [L]
f(x) ................................................ função genérica
g ............................................. aceleração da gravidade .............................................. [LT-2]
h ............................................................ altura ................................................................. [L]
hb ..................... pressão de borbulhamento ou pressão de entrada de ar ............................ [L]
mt .............................................. massa total de soluto .................................................... [M]

q .............................................. fluxo de água do solo ........................... [L3L-2T-1] ou [LT-1]
xiv
qi .................................. fluxo de água na direção i (i = x, y, z) .............. [L3L-2T-1] ou [LT-1]

q b ..................................... fluxo dispersivo pelo contorno .................................... [ML-2T-1]
t ............................................................ tempo ................................................................ [T]

v ......................................... velocidade linear média real .......................................... [LT-1]
v .................................. módulo da velocidade linear média real .................................. [LT-1]
vi ......................... velocidade linear média real na direção i (i = x, z) .......................... [LT-1]
v̂ i ................................................ aproximação de vi .................................................. [LT-1]
{vi} ...................................... vetor de valores nodais de vi

v c ................................... velocidade linear média corrigida ...................................... [LT-1]
vic .................... velocidade linear média corrigida na direção i (i = x, z) ...................... [LT-1]
x, y, z ............... direções, coordenadas cartesianas (z = profundidade) ............................. [L]
i ................................... intervalo na direção i (i = x, y, z, t) ................................ [L] ou [T]
V .................................... volume elementar representativo .......................................... [L3]
............................... contorno do domínio (contorno do sistema)
ijkl
................................. dispersividade geométrica do meio .......................................... [L]
L
......................................... dispersividade longitudinal ................................................ [L]
T
.......................................... dispersividade transversal ................................................. [L]
1,
2
.................. constantes de generalização dos modelos de sorção
........................................... função delta de Kronecker
ij
........................................ potencial total da água do solo .............................. [JM-1] ou [L]
i
.................................... fontes e retiradas de massa de soluto ............................... [ML-3T-1]
........................................................ porosidade ........................................... [adimensional]
............................................ permeabilidade intrínseca ............................................... [L2]
.............................................. condutividade relativa .................................. [adimensional]
r
.................... parâmetro de ponderação de tempo do esquema numérico
' .............. fator geral de primeira ordem de decaimento de massa de soluto .......... [ML-3T-1]
l
... constante de primeira ordem de decaimento de massa de soluto na fase líquida ....... [T-1]
s
... constante de primeira ordem de decaimento de massa de soluto na fase sólida ........ [T-1]
w
....................................... viscosidade dinâmica da água .................................... [ML-1T-1]
................................. fração volumétrica ou conteúdo de água .................... [adimensional]
xv
ˆ ................................................. aproximação de ...................................... [adimensional]
res
......................................... conteúdo residual de água ............................... [adimensional]
sat
..................................... conteúdo de água na saturação ............................ [adimensional]
w
........................................... massa específica da água ............................................ [ML-3]
' ................... fator geral de ordem zero de produção de massa de soluto ............... [ML-3T-1]
l
.......... termo de ordem zero de produção de massa de soluto na fase líquida ....... [ML-3T-1]
s
........... termo de ordem zero de produção de massa de soluto na fase sólida ............... [T-1]
.......................................... fator de tortuosidade do solo ............................. [adimensional]
................................. função de ponderação relativa ao nó i
i
...................................... função de base associada ao nó i
i
........................................... potencial mátrico do solo ................................. [L2T-2] ou [L]
ˆ ............................................... aproximação de
...................................... [L2T-2] ou [L]
e
............................................ potencial eletroquímico .................................. [L2T-2] ou [L]
g
............................................ potencial gravitacional ................................... [L2T-2] ou [L]
o
............................................... potencial osmótico ...................................... [L2T-2] ou [L]
p
.............................................. potencial de pressão ..................................... [L2T-2] ou [L]
t
................................................. potencial térmico ....................................... [L2T-2] ou [L]
xvi
1- INTRODUÇÃO
1.1- CONTEXTO GERAL
A água subterrânea vem se tornando uma das principais fontes de água potável para a
humanidade. A deterioração de nossos rios e lagos e a crescente demanda de água pelas
populações vem exigindo cada vez mais o seu uso, especialmente em zonas rurais. Segundo
Solley et al. (1988), apud Fetter (1993), mais da metade da população dos Estados Unidos é
abastecida com água subterrânea. No Brasil não existe confirmação em números, mas estimase que cerca de 51% da população é abastecida com água de poços subterrâneos (Foster e
Hirata, 1993).
A água subterrânea é aquela água que se armazena em formações geológicas no
subsolo, conhecidas como aqüíferos. Comumente, a água é ali encontrada com melhor
qualidade em relação às águas superficiais. Por isso, o seu aproveitamento, isto é, a sua
potabilização, apresenta custos muito menores, necessitando muitas vezes apenas de algum
processo corretivo e de algum tipo de desinfecção. A manutenção dessa qualidade exige
esforços no sentido de se evitar possíveis contaminações.
A maior parte dos incidentes de contaminação envolvem substâncias líquidas liberadas
sobre ou logo abaixo da superfície dos terrenos. Essas substâncias infiltram-se no solo,
atravessando-o, de modo que normalmente afetam primeiramente a água menos profunda,
localizada na região do solo logo acima da superfície da água acumulada, denominada zona
não-saturada, dirigindo-se então para o aqüífero. Como a água no subsolo está em constante
movimentação, essa contaminação pode se propagar por todo o aqüífero.
Dessa forma, a zona não-saturada apresenta-se como uma primeira barreira à
passagem de contaminantes da superfície para os aqüíferos, formando uma camada de
atenuação dos efeitos danosos da contaminação (Foster e Hirata, 1993). Neste contexto, é de
extrema importância que sejam conhecidos os processos de fluxo de água e de transporte e
armazenamento de contaminantes nessa zona. Para isso, faz-se necessário o uso de
ferramentas que possibilitem, em tempo hábil, o conhecimento e a previsão dos processos
envolvidos nessa zona. Uma dessas ferramentas é a modelagem matemática dos processos
físicos envolvidos. Um bom modelo matemático pode ajudar a prever situações de risco e,
muitas vezes, indicar soluções capazes de resolver problemas graves, desde que o modelo seja
bem manipulado e tenha sido bem calibrado e suficientemente testado.
1
1.2- OBJETIVO
O objetivo deste trabalho foi o desenvolvimento de um modelo matemático para a
simulação do transporte de contaminantes solúveis dentro da zona não-saturada do solo. Esse
modelo foi concebido a partir de um modelo de simulação de fluxo não-saturado já existente,
desenvolvido por Koide (1990) e modificado por Campos (1998). Para isso, numericamente,
o processo de transporte foi dividido em uma parte advectiva e uma dispersiva. O modelo
emprega o Método das Características para a simulação da parte advectiva e o Método dos
Elementos Finitos para a simulação da parte dispersiva.
Os testes de verificação do modelo foram feitos com base em problemas apresentados
na literatura, cujas soluções analíticas são conhecidas, e buscam avaliar o desempenho do
modelo sob diversas situações aplicáveis.
1.3- ESTRUTURA DA DISSERTAÇÃO
Esta dissertação foi estruturada em capítulos que tratam especificamente de assuntos
relativos ao tema.
O capítulo 2 descreve os princípios físicos envolvidos na movimentação da água no
subsolo e sua interação com os solos, descrevendo também os princípios físicos do transporte
de contaminantes solúveis. Nesse capítulo, são desenvolvidas as equações básicas que regem
o problema.
O capítulo 3 apresenta de forma resumida alguns conceitos relativos à modelagem de
sistemas físicos e a classificação dos diversos tipos de modelos.
No capítulo 4, são apresentados alguns conceitos envolvendo a modelagem
matemática e descritos alguns métodos numéricos aplicáveis.
O capítulo 5 discute alguns dos mais recentes desenvolvimentos na área da
modelagem matemática do transporte de contaminantes em solos, buscando dar uma maior
ênfase ao transporte na zona não-saturada. Ao final, é explicada a seleção do método
numérico empregado no trabalho.
No capítulo 6, é apresentada, de forma detalhada, a formulação do modelo
desenvolvido.
2
No capítulo 7, é discutida a estrutura do programa, escrito em linguagem FORTRAN
77, e são apresentados os problemas para teste do modelo. Os testes procuram cobrir diversas
situações simuláveis pelo modelo.
O capítulo 8 traz os resultados obtidos com os testes do modelo, analisando-os
separadamente para cada teste. Ao final, é feita uma análise geral.
Finalmente, as conclusões são apresentadas no capítulo 9, onde também são feitas
recomendações para pesquisas futuras.
3
2- FUNDAMENTAÇÃO TEÓRICA
Existem diversas fontes de contaminação que colocam em risco a qualidade das águas
subterrâneas. Podem ser das mais variadas formas, tanto projetadas quanto acidentais. Além
disso, existe uma grande variedade de substâncias contaminantes para a água subterrânea.
Fetter (1993) fornece uma extensa lista dessas substâncias. O comportamento desses
contaminantes no subsolo depende, em boa parte, da sua interação com a água. Cada material
dissolve-se na água em um certo grau, podendo ser completamente, parcialmente ou não se
dissolver. Normalmente, os contaminantes são classificados em solúveis e não solúveis. Os
materiais solúveis, também conhecidos como solutos, misturam-se completamente com a
água, passando a fazer parte da sua massa, sendo por isso denominados de Fase Líquida
Aquosa (APL – Aqueous Phase Liquid). Os materiais líquidos não solúveis ou parcialmente
solúveis não se dissolvem completamente e não se misturam a ela. Por isso, são denominados
de Fase Líquida Não-Aquosa (NAPL – Non-Aqueous Phase Liquid). Normalmente, os
NAPL’s possuem massa específica diferente da água e, por isso, são divididos em duas
categorias: os mais leves, cuja massa específica é menor que a da água, sendo denominados
de Fase Líquida Não-Aquosa Leve (LNAPL – Light Non-Aqueous Phase Liquid), e os
mais pesados, de massa específica maior que a da água, denominados de Fase Liquida NãoAquosa Densa (DNAPL – Dense Non-Aqueous Phase Liquid). Neste trabalho, o modelo
desenvolvido considera somente o transporte de solutos, cujo enfoque será dado nos itens
seguintes.
O maior problema da contaminação da água subterrânea é o fato de ser um processo
de longa duração (Fetter, 1993), devido à baixa velocidade com que a água se move no
subsolo. Assim, quando alguma contaminação é detectada, a sua remediação pode ser
extremamente lenta e bastante dispendiosa. Fetter (1993) lista uma série de casos que
evidenciam isso. Por outro lado, porém, o fato de ser lento beneficia o controle, se o problema
for detectado a tempo.
2.1- DISTRIBUIÇÃO DA ÁGUA NA SUBSUPERFÍCIE
A Figura 2.1 mostra esquematicamente uma seção vertical de um terreno. De acordo
com a proporção de água dentro dos vazios do solo, a subsuperfície pode ser dividida em duas
regiões: uma zona saturada, onde a água presente ocupa todos os vazios da matriz sólida, e
4
uma zona não-saturada, onde se considera que a água não preenche todos os vazios e, por
isso, há a presença de gases. No entanto, o principal fator de distinção entre as duas zonas está
relacionado com a pressão da água no interior dos poros. Na zona saturada, a água encontra-se
sob uma pressão igual ou superior à pressão atmosférica (pressão relativa nula ou positiva)
enquanto que na zona não-saturada a pressão da água é inferior à atmosférica (pressão relativa
negativa). Com isso, inclui-se nessa última zona a região onde os poros podem estar cheios de
água, porém com pressão relativa negativa, como é o caso da região conhecida como franja
capilar.
Zona das raízes
Zona não-saturada
pressão < 0
Zona intermediária
Nível freático
pressão = 0
Franja capilar
Zona saturada
pressão > 0
Figura 2.1- Seção vertical de um terreno, mostrando as regiões da subsuperfície
A zona não-saturada é geralmente subdividida em três zonas:
a) zona das raízes ou zona do solo: compreende a região do terreno até a profundidade
normal de crescimento do sistema radicular das plantas. Normalmente, é considerada
tendo espessura de 1m, muito embora não tenha profundidade definida;
b) zona intermediária: é considerada desde o limite inferior da zona das raízes até o limite
superior da franja capilar; e
c) zona ou franja capilar: compreende a região logo acima da zona saturada. Nela os poros
estão saturados com água, porém com pressão negativa devido a efeitos de
capilaridade. A sua espessura depende da estrutura e da textura do solo.
O fluxo subsuperficial é um importante componente do ciclo hidrológico. A água que
atinge a superfície do terreno através de chuvas, irrigações etc., pode infiltrar-se, espalhandose pelo terreno. Uma parte pode evaporar, outra percolar para baixo. Essa água que percola
pode atingir camadas mais profundas e menos permeáveis de solo e acumular-se, criando
assim a zona saturada, constituída principalmente por aqüíferos. A superfície contínua dessa
zona possui pressão nula (igual à atmosférica) e é denominada superfície freática. É essa
5
superfície que separa a zona saturada da não-saturada. Portanto, de um modo menos formal, a
zona não-saturada é a região do subsolo que se estende da superfície do terreno até a
superfície freática (Nielsen et al., 1986; Ward e Robinson, 1990; Fetter, 1993).
2.2- O SOLO E SUAS CARACTERÍSTICAS
O solo é um material poroso, bastante complexo, que normalmente é resultante das
ações de intemperismo sobre rochas. Comumente, apresenta-se como um sistema constituído
por três fases:
• fase sólida: também chamada de matriz sólida ou matriz porosa (Bear, 1972; Voss, 1984),
compreende as partículas sólidas do solo;
• fase líquida: compreende a água e os materiais dissolvidos que se encontram nos poros
entre as partículas. Também é denominada de solução de solo (Ward e Robinson,
1990); e
• fase gasosa: constituída pelo ar, pelo vapor de água e por outros gases que eventualmente
podem estar presentes nos poros.
Como essas três fases não estão em constante equilíbrio, as proporções de cada uma
dentro do solo varia com o tempo. Para um solo completamente saturado, a fase gasosa
inexiste por completo, de modo que os poros da matriz sólida estão preenchidos pela fase
líquida. Para um solo não-saturado, podem existir 4 possíveis situações (Santos Neto, 1990),
que podem ser visualizadas na Figura 2.2:
i) solo com baixo grau de saturação com gás contínuo e água descontínua. A fase líquida
apresenta pequeno volume relativamente à fase gasosa, de modo que se formam
meniscos em torno das partículas sólidas;
ii) solo com grau de saturação intermediário onde as fases líquida e gasosa possuem
volumes relativamente semelhantes e apresentam-se ambas em forma contínua;
iii) solo com elevado grau de saturação onde o volume da fase líquida é amplamente
superior ao da fase gasosa, de modo que a primeira apresenta-se contínua e a
segunda, descontínua, aparece como bolhas dispersas na fase líquida; e
iv) solo com grandes vazios de gás, dispersos numa matriz sólida saturada. Nesse caso, as
bolhas da fase gasosa são muito maiores que as partículas da matriz sólida que
encontram-se unidas pelas forças de adesão da fase líquida.
6
partícula sólida
gás
partícula sólida
gás
água
i) gás contínuo, água
descontínua
partícula sólida
gás
água
ii) gás contínuo, água
contínua
partícula sólida
água
água
iii) gás descontínuo,
água contínua
gás
iv) grandes bolhas de
gás dispersas numa
matriz saturada
Figura 2.2- Possíveis situações que um solo não-saturado pode apresentar
Fonte: Santos Neto (1990), modificado
Os volumes das fases líquida e gasosa armazenáveis em um solo dependem de duas
características básicas do mesmo: a sua textura e a sua estrutura. A textura de um solo é
definida a partir da sua distribuição granulométrica, isto é, a distribuição dos tamanhos dos
grãos sólidos. Esses grãos podem agrupar-se formando unidades maiores, chamadas
agregados, que estão separadas entre si por planos de fraqueza (Ward e Robinson, 1990). A
estrutura do solo resulta do arranjo desses agregados. A fase líquida presente no solo pode
ocupar tanto os vazios entre os grãos (vazios texturais) quanto entre os agregados (vazios
estruturais).
As partículas da matriz sólida podem ser de variados tamanhos, dependendo do tipo de
material primário, isto é, da rocha originária, e do tipo de intemperismo atuante, podendo ser
7
distribuídas em quatro classes predominantes: argilas, siltes, areias e pedregulhos. Não existe
um consenso geral em relação aos tamanhos dos grãos em cada classe. Bear (1972) apresenta
um conjunto de classificações normalmente usadas em estudos de solo.
Algumas das principais características do solo são determinadas pela fração de argila
presente. Geralmente, essa fração é considerada como sendo formada por partículas com
diâmetros efetivos menores que 2 m (Bear, 1972; Yong et al., 1996; Hillel, 1998) e consiste
de partículas minerais conhecidas como minerais de argila. Uma descrição detalhada desses
minerais é fornecida por Yong et al. (1996).
Os minerais de argila possuem uma alta superfície específica, isto é, alta relação área
superficial por volume unitário. Além disso, possuem uma carga elétrica negativa não
balanceada na sua superfície. Essa carga superficial pode afetar a estrutura do solo. O
principal efeito da carga superficial das argilas é a criação de um campo elétrico que atrai
cátions para as proximidades, criando a conhecida dupla camada eletrostática. Mais
detalhes podem ser encontrados em Nielsen et al. (1986) e Foster e Hirata (1993).
2.2.1- Características dos Solos Relacionadas à Água
2.2.1.1- Porosidade
A matriz sólida é composta de grãos sólidos e de poros ocupados por líquidos e gases.
A determinação do volume absoluto desses últimos poderia ser um fator conclusivo para a
caracterização do solo. Porém, sua proporção em relação ao volume total é muito mais útil.
A porosidade é definida como a relação entre o volume de vazios e o volume total da
matriz sólida
Vv
Vt
onde
(2.1)
= porosidade [adimensional],
Vv = volume de vazios [L3] e
Vt = volume total da matriz sólida [L3], isto é, o volume de vazios mais o volume
dos grãos sólidos.
8
2.2.1.2- Fração Volumétrica de Água
O volume que a fração líquida ocupa dentro da matriz sólida é caracterizado pela
fração volumétrica de água, também chamada conteúdo de umidade ou conteúdo de água.
É expresso como a relação entre o volume da fase líquida e o volume total da matriz sólida
Vw
Vt
onde
(2.2)
= conteúdo de água [adimensional] e
Vw = volume da fase líquida presente [L3].
Também pode ser definido o grau de saturação (S) [adimensional], que expressa, em
porcentagem, a relação entre o volume de água e o volume de vazios, ou seja,
S
Vw
100
Vv
(2.3)
2.2.1.3- Ponto de Murchamento e Capacidade de Campo
Como já descrito, a água num solo não-saturado encontra-se sob uma pressão relativa
negativa (sucção) que depende do conteúdo de água ( ). Assim, para uma planta conseguir
extrair água do solo, ela precisa produzir, através de suas raízes, uma sucção sobre a água
maior que a do solo. O ponto de murchamento de um solo é definido como o conteúdo
mínimo de água presente no solo em que as plantas conseguem extrair água, ou seja, quando a
sucção máxima produzida pela planta iguala-se à sucção do solo. Embora a sucção máxima
que se consegue produzir dependa do tipo de planta e do seu estado de crescimento, em geral,
diferenças pouco significativas entre diferentes plantas são obtidas em um mesmo solo (Ward
e Robinson, 1990), de modo que o valor do ponto de murchamento é usualmente considerado
constante para o solo.
Se uma amostra de solo saturado é posta em repouso, a água que inicialmente está
presente nos poros tende a movimentar-se para baixo, por efeito da gravidade, drenando o
9
solo. Após um certo período de tempo, verifica-se que uma certa quantidade de água continua
retida entre os grãos. Esse conteúdo de água remanescente é definido como a capacidade de
campo do solo. O tempo necessário para a drenagem não é estabelecido quantitativamente,
mas, para situações de campo, é normalmente tomado como sendo 48 horas após o término da
chuva que levou o solo à saturação (Ward e Robinson, 1990).
2.2.1.4- Condutividade Hidráulica
A condutividade hidráulica é uma medida da capacidade de um meio poroso de
permitir a passagem de um fluido (Fetter, 1993). Expressa a velocidade média com que o
fluido pode atravessar os poros da matriz sólida quando está sujeito a um gradiente hidráulico
unitário, podendo ter um valor diferente para cada direção. É uma função do tamanho,
continuidade e forma dos vazios, isto é, da estrutura e textura do solo, e das características do
fluido. Com isso, o seu valor em geral é mais alto para solos de granulometria grosseira e
mais baixo para solos de granulometria fina.
Na zona saturada, o valor da condutividade hidráulica, para uma dada direção, pode
ser considerado invariável, qualquer que seja o tipo de solo, mas para a zona não-saturada o
seu valor varia consideravelmente com o conteúdo de água. No Item 2.3.5, será feita uma
descrição mais detalhada do assunto.
2.3- DESCRIÇÃO DO FLUXO EM SOLO NÃO-SATURADO
A presença de ar nos poros de um solo não-saturado afeta fortemente as suas
características hidráulicas, de modo que a descrição do fluxo nessa zona torna-se amplamente
complexa.
2.3.1- Interações Entre a Água e o Solo
A água presente em um solo tende a movimentar-se para baixo por ação da força
gravitacional. No entanto, quando o seu conteúdo está abaixo da capacidade de campo do
solo, não ocorre essa drenagem e a água fica retida nos poros da matriz sólida, devido a forças
10
características da interação entre a água e o solo. As principais forças de retenção atuantes são
as de capilaridade, adsorção e osmose (Ward e Robinson, 1990).
As forças capilares são resultantes da tensão superficial na interface entre a água e o
ar. A atração entre as moléculas de água no líquido é maior que a das moléculas de vapor.
Com isso, a superfície livre da água líquida tende a se contrair formando uma película. No
entanto, quando a água encontra uma face sólida, ela pode aderir-se a essa face com uma
força maior ou menor que a atração entre suas moléculas. No caso dos solos, é comum que a
adesão seja maior, de modo que, nos contatos das faces sólidas com a água, essa sofra uma
elevação e, fora desses contatos, sofra um abaixamento. A força capilar é o resultado da
combinação dessas duas forças. Assim, a superfície livre da água fica encurvada, formando
um menisco, e a intensidade da força capilar dependerá dessa curvatura. Por isso, a água é
retida mais fortemente em poros menores, como os dos solos de textura fina, onde os
meniscos possuem uma curvatura mais acentuada.
As forças de adsorção são resultantes principalmente da atração eletrostática entre as
moléculas de água e as superfícies das partículas sólidas. Como as forças eletrostáticas só são
efetivas em pequenas distâncias, as moléculas de água formam finíssimos filmes sobre as
partículas do solo, quando sujeitas somente a essas forças. A quantidade de água adsorvida
eletrostaticamente depende fundamentalmente da área superficial das partículas, de modo que
aquelas com maior área específica retém mais água. Isto ajuda a explicar porque as argilas
retém mais a água do que outros tipos de solo (Ward e Robinson, 1990).
Quando a água possui solutos dissolvidos, surge então a força osmótica. Ela é o efeito
resultante da diferença de concentração de soluto através de uma membrana permeável. Nos
solos, essa membrana permeável pode ser uma barreira tal como uma bolha de ar que ocupe
um poro e que permite que moléculas de água, sob a forma de vapor, a atravessem no sentido
da solução menos concentrada para a mais concentrada. Com isso, há uma pequena atração
entre as soluções em ambos lados da bolha (Hillel, 1998).
A combinação das forças capilares e de adsorção resulta numa sucção denominada
sucção mátrica ou matricial, pois os seus efeitos geralmente atuam conjuntamente no
sentido de reter a água sobre as partículas da matriz sólida. O seu efeito pode ser visualizado
na Figura 2.3.
A combinação de todas as forças de retenção resulta numa sucção denominada sucção
total, que varia com a textura e a estrutura do solo, bem como com o seu conteúdo de água.
11
água capilar
água adsorvida
partículas
Figura 2.3- Efeito das forças capilares e de adsorção
Fonte: Ward e Robinson (1990), modificado
2.3.2- Estado Energético da Água do Solo
O movimento da água num solo não-saturado é dirigido principalmente pelas suas
energias potenciais. Logicamente, a água deslocando-se num solo possui velocidade e,
portanto, possui também energia cinética. No entanto, devido às baixas velocidades com que
se processa o movimento, essa energia é desprezível, de modo que consideram-se apenas as
energias potenciais. Por isso, o estado energético da água no solo é normalmente denominado
de potencial da água do solo.
Conceitualmente, esse potencial expressa a energia potencial total da água do solo em
relação a um referencial padrão. A energia potencial de qualquer corpo é resultante da ação de
forças que possibilitem ou restrinjam o seu movimento. No caso da água do solo, as principais
forças atuantes são as forças capilares, de adsorção, osmótica e eletroquímicas, além da força
gravitacional.
Portanto, o potencial total da água do solo é dado por
g
onde
p
o
= potencial total da água do solo,
g
= potencial gravitacional,
p
= potencial de pressão,
o
= potencial osmótico,
e
= potencial eletroquímico e
t
= potencial térmico.
12
e
t
(2.4)
Cada um desses potenciais pode ser expresso em unidades de energia por unidade de
massa da água [L2T-2]. Porém, é muito mais usual que sejam usadas dimensões de
comprimento [L], sendo conhecidos, dessa forma, como cargas, expressando a energia da
água por unidade de peso.
No estudo de potenciais, deve ser levado em consideração um nível de referência. Em
geral, para os potenciais de pressão, considera-se como referência a pressão atmosférica.
Dessa forma, uma carga nula equivale ao potencial produzido por uma pressão igual à pressão
atmosférica.
Assim, para solos não-saturados, o potencial de pressão, que engloba as forças
capilares e de adsorção, representa a sucção produzida pelo solo, sendo denominado de
potencial mátrico, de valores negativos. Em solos saturados, o potencial de pressão engloba
as forças devidas à pressão hidrostática, representando o potencial piezométrico, com valores
positivos. Em solos não-saturados, o potencial gravitacional aumenta com a elevação e tende
a drenar o solo, caso as forças de retenção não sejam suficientes (Ward e Robinson, 1990).
O potencial osmótico e o potencial eletroquímico são provocados pela presença de
solutos e podem ser desprezados em estudos de fluxo, por conseqüência de seu baixo valor
em comparação com os outros potenciais. O potencial térmico é importante somente quando o
problema trata de fluidos aquecidos em condições não isotérmicas (Santos, 1996), sendo que
pode ser desprezado na maioria dos casos, inclusive neste trabalho.
Assim, neste estudo, o potencial total pode ser descrito como a soma dos potenciais
gravitacional e de pressão, ou seja,
g
p
(2.5)
Para os potenciais gravitacionais, adota-se, em geral, um datum ou nível no terreno
como referência. Pode ser tomado qualquer nível. Na prática, é usual, em solos não-saturados,
usar como datum a superfície do terreno. Assim, em termos de carga, o potencial
gravitacional é tomado como a posição do ponto em consideração medida em um eixo vertical
com sentido para cima, isto é,
z
onde
= potencial mátrico [L] e
13
(2.6)
z = profundidade da água do solo [L].
A Figura 2.4 ilustra o comportamento do potencial total, em relação ao eixo vertical z.
z
Superfície do terreno
Plano de
fluxo zero
p
g
Nível freático
g
p
Figura 2.4- Perfil do potencial total da água do solo ( ) com a posição (z), mostrando os
potenciais gravitacional (
g)
e de pressão (
p)
Fonte: Ward e Robinson (1990), modificado
Fisicamente, a água move-se de um ponto de potencial total mais alto para um de
potencial total mais baixo, não importando os valores e sim a diferença entre os pontos. De
modo mais específico, o que rege o movimento da água do solo são os gradientes do potencial
total (
). Podem existir locais onde esses gradientes são nulos, de modo que a água, nesses
locais, não possui movimento. Na Figura 2.4, o plano de fluxo zero é o conjunto desses
pontos e representa um plano de máximo potencial (
= 0). Acima desse plano, o fluxo se dá
com sentido à superfície, com o contrário ocorrendo abaixo.
2.3.3- Curva Característica do Solo: a Relação
O potencial mátrico ( ), que representa a sucção do solo, está intrinsecamente
relacionado com o conteúdo de água ( ). Normalmente, diz-se que
14
é uma função de . Com
isso, essa relação pode ser expressa por um gráfico, denominado curva característica do solo
ou curva de retenção de umidade. A Figura 2.5 mostra uma curva característica para um
solo genérico.
Quando o solo está saturado, o potencial mátrico é nulo. Se então lhe é aplicada uma
sucção gradualmente crescente, verifica-se que, no início, praticamente nenhuma água é
retirada, de modo que o seu conteúdo permanece inalterado e igual ao conteúdo de água na
saturação (
sat).
A sucção aumenta até que se atinja uma sucção limite na qual os poros
começam a drenar e permitem a entrada de ar no solo. Essa sucção é conhecida como pressão
de borbulhamento ou pressão de entrada de ar (hb). A partir daí, uma certa quantidade de
água é perdida para cada aumento na sucção aplicada, fazendo com que o conteúdo de água
diminua consideravelmente. Com o aumento da sucção, chega-se a um ponto onde não mais
se consegue retirar água do solo. Daí em diante, o conteúdo de água torna-se novamente
constante para qualquer aumento na sucção. Esse último conteúdo de água é denominado
conteúdo residual de água (
res).
Esse comportamento é muito comum para a maioria dos
solos, porém alguns têm comportamento diferente, não sendo tratados aqui.
-106
(cm)
-104
Potencial mátrico
res
-102
= 0,1
hb
sat
= 0,5
-100
0
0,1
0,2
0,3
0,4
Conteúdo de água
0,5
Figura 2.5- Curva característica de um solo genérico
Fonte: Fetter (1993), modificado
15
0,6
A forma da curva característica depende, basicamente, do tipo de solo, com sua textura
e estrutura. Como a curva é uma resposta à entrada de ar nos poros, especialmente na pressão
de borbulhamento, ela tem uma curvatura mais acentuada para solos arenosos do que para
solos argilosos. Quando um solo está sendo drenado, os poros maiores, que suportam uma
menor força de retenção, permitem mais rapidamente a entrada de ar. Os poros menores
somente permitem a entrada de ar a sucções mais altas. Como, nos solos arenosos, a maioria
dos poros é relativamente grande, eles drenam quase todos ao mesmo tempo, com uma queda
do conteúdo de água mais acentuado. Por isso, o efeito da pressão de borbulhamento é mais
visível nesses tipos de solo. Nos solos argilosos, os poros são menores e a entrada de ar ocorre
lentamente, de modo que a curva característica tem curvatura bem menos acentuada.
2.3.3.1- Efeitos de Histerese
A principal dificuldade na obtenção e utilização de curvas características é
determinada pelos efeitos de histerese. De uma maneira geral, o processo de secagem de um
solo apresenta a relação do potencial mátrico com o conteúdo de água diferente da obtida no
umedecimento. Com isso, o potencial mátrico não vai depender somente do conteúdo de água
presente, mas também de o solo estar sendo umedecido ou secado. A Figura 2.6 ilustra os
efeitos da histerese.
Os poros esvaziam-se a um potencial mátrico normalmente maior do que quando são
preenchidos. Isto ocorre, principalmente, devido a três efeitos (Nielsen et al., 1986; Ward e
Robinson, 1990; Hillel, 1998):
i) efeito de gargalo: esse efeito é provocado principalmente pela geometria do poro. O poro
em um trecho mais estreito necessita de maior sucção para drenar do que em trecho
mais largo, pois a curvatura do menisco no seu interior é mais acentuada;
ii) efeito ângulo de contato: o ângulo de contato entre o fluido e as partículas da matriz
sólida é maior quando o solo está sendo umedecido do que quando ele está secando.
Isso acontece porque a água em contato com as partículas tende a manter a posição
inicial do contato, movimentando-se assim com uma defasagem em relação à água
no centro do poro. Com isso, o menisco tem maior curvatura durante a secagem do
que durante o umedecimento, o que produz uma maior sucção na secagem; e
iii) efeito do ar aprisionado: para uma mesma sucção mátrica, o conteúdo de água tende a
ser menor no umedecimento do que na secagem, pois alguma quantidade de ar pode
16
eventualmente ficar aprisionada nos poros. Esse efeito é mais evidente quando o solo
está sob sucção nula, pois a saturação não chega a ser total quando é umedecido. O ar
aprisionado, com o tempo, pode se dissolver, permitindo a saturação completa do
solo.
(cm)
-104
Potencial mátrico
-106
-102
Secagem
Umedecimento
-100
0
0,1
0,2
0,3
0,4
Conteúdo de água
Figura 2.6- Histerese na relação
0,5
0,6
, mostrando as curvas intermediárias
Nos solos argilosos, os efeitos de histerese podem ser acentuados pelos efeitos de
inchamento e encolhimento dos minerais de argila. Quando ocorre umedecimento, a matriz
porosa incha, aumentando o volume dos poros e, portanto, o seu volume total. O contrário
ocorre na secagem. Isso porque os minerais de argila podem absorver grandes quantidades de
água e íons entre suas partículas (Fetter, 1993; Hillel, 1998).
Quando o processo de secagem (ou umedecimento) não é completo e o processo
inverso é iniciado, o solo seguirá um caminho intermediário, a partir do ponto onde parou, até
atingir a curva característica novamente, como mostrado na Figura 2.6.
A dificuldade de se obter e/ou modelar precisamente a histerese em solos é muito
grande, de modo que, na prática, ela é geralmente ignorada, sendo usada somente a curva de
secagem ou a de umedecimento, pois os erros inerentes a essa aproximação em geral são
pequenos em termos globais.
17
2.3.4- Equacionamento do Fluxo Não-Saturado
Segundo Fetter (1993), o primeiro a reconhecer as leis básicas do fluxo não-saturado
foi Buckingham. No seu trabalho, publicado em 1907, ele descrevia que o potencial mátrico
( ) do solo era dependente do seu conteúdo de água ( ), da temperatura e da massa específica
do solo. Descrevia também que o fluxo de água através de uma seção era proporcional ao
gradiente do potencial mátrico do solo. A proporcionalidade entre o fluxo e o gradiente é dada
por um fator K( ) denominado de condutividade hidráulica, que é uma função do potencial
mátrico do solo. Como o potencial mátrico do solo depende do seu conteúdo de água, alguns
autores preferem expressar a condutividade hidráulica como uma função do conteúdo de água,
ou seja, K( ).
No entanto, somente em 1928, Richards formalizou as idéias de Buckingham e
estendeu o conceito ao gradiente do potencial total. A lei de fluxo de Buckingham é então
dada vetorialmente por

q
onde
K
(2.7)

q = fluxo de água do solo [LT-1],
K( ) = condutividade hidráulica do solo não-saturado para um dado potencial
mátrico [LT-1] e
= gradiente do potencial total de água do solo [adimensional].
Deve ser notada a semelhança dessa equação com a equação de Darcy para o fluxo
saturado, podendo ser considerada, então, como uma extensão.

O fluxo de água do solo ( q ) é um vetor que atua na direção do gradiente. Representa a
velocidade média com que a água atravessa uma seção de área unitária do solo, tendo por isso
dimensões de velocidade.
Para o volume elementar representativo de um solo não-saturado, mostrado na Figura
2.7, cujos lados tem comprimentos x, y e z, o volume é dado por V = x y z. A lei de
conservação de massa ou equação de continuidade para esse volume pode ser estabelecida
pelo princípio: a massa de água que entra no volume menos a que sai é igual à taxa de
variação de água nesse volume. A descarga de água através de uma face i é o produto do fluxo
qi pela área daquela face. A variação do fluxo dentro do volume elementar ao longo da
direção i é dada por ( qi/ i) i. Portanto, a equação da continuidade será
18
qx y z qy x z qz x y
qz
qx
x
x
qx
qz
z
z
y z
qy
qy
y
y
x z
(2.8)
V
x y
t
qz + qz / z z
z
qx
y
y
x
z
qy + qy / y y
qy
x
q x + qx / x x
qz
Figura 2.7- Volume elementar representativo do solo, mostrando o balanço de massa da água
do solo
Reorganizando a equação e considerando que não há variação de volume da matriz
sólida, a equação da continuidade se reduz a
qy
qx
x
t
y
qz
z
(2.9)
Em notação vetorial

q
t
(2.10)
Combinando as Equações (2.7) e (2.10), chega-se à equação geral do fluxo nãosaturado, conhecida como equação de Richards (Richards, 1931, apud Fetter, 1993).
t
K
19
(2.11)
Substituindo o valor de
t
Como
da Equação (2.6), vem
K
z
t
K
K
z
(2.12)
z = 1, a Equação (2.12) se torna a equação geral de fluxo em solos não-
saturados
t
K
K
(2.13)
A obtenção da Equação (2.13) leva em consideração várias hipóteses simplificadoras:
matriz sólida indeformável, temperatura e pressão do ar constantes, água incompressível,
densidade da água independente da concentração de solutos e invariante sobre o domínio.
Também assume que a presença de ar pode ser ignorada, exceto quando afeta o valor de K
(Fetter, 1993). Nielsen et al. (1986) descrevem as implicações práticas que essas hipóteses
podem ter.
2.3.5- A Condutividade Hidráulica do Solo Não-Saturado
A condutividade hidráulica é uma medida da capacidade de um meio poroso de
permitir a passagem de um fluido. Em um sistema solo-água-ar, quando o solo está saturado,
todos os poros estão preenchidos por água, de modo que a grande maioria consiste em
conexões que permitem a passagem dessa. A condutividade, nesse caso, é máxima e, como os
poros teoricamente não perdem volume, a capacidade desse solo de transmitir água é
praticamente invariável e pode ser considerada constante.
No entanto, em solos não-saturados, a presença de ar complica consideravelmente o
comportamento da água. Como o fluxo de água se dá por dentro de conexões, quanto menor
essa conexão menor a capacidade do poro de transmitir água. Assim, se uma bolha de ar
ocupa um poro parcialmente, a pequena conexão formada permite a passagem de uma
pequena quantidade de água. Se ocupar o poro inteiramente, a conexão é quebrada e a água
não tem condições de passar por aquele poro. Dessa forma, fica evidente a forte influência
que o conteúdo de água tem sobre o valor da condutividade hidráulica, ou seja, a
20
condutividade é realmente uma função do conteúdo de água e, conseqüentemente, do
potencial mátrico do solo. Como se vê, não é necessário apenas que haja água disponível ao
fluxo, mas também que essa água esteja interconectada entre os poros. Esse efeito pode ser
visualizado na Figura 2.8.
água
ar
água
partícula
partícula
(b)
Sem a bolha o fluxo é livre
(a)
A bolha de ar impede o fluxo
Figura 2.8- Efeito da presença de ar sobre a condutividade de um poro
Devido ao fato de ser dependente do conteúdo de água do solo, que é histerético em
relação ao potencial mátrico, a condutividade hidráulica também apresenta comportamento
histerético, o que a torna um fator complicador na resolução das equações de fluxo. Segundo
Nielsen et al. (1986), apesar de ser pouco pronunciada, a histerese de K não deve ser
desconsiderada. No modelo de fluxo utilizado neste trabalho (Koide, 1990; Campos, 1998), o
efeito da histerese é desconsiderado.
A condutividade também é influenciada pela temperatura. Constantz (1982), apud
Nielsen et al. (1986) e Fetter (1993), obteve a seguinte relação
K
r
w
g
(2.14)
w
onde
r(
) = condutividade relativa [adimensional],
= permeabilidade intrínseca [L2],
w
= massa específica da água [ML-3],
g = aceleração da gravidade [LT -2] e
w
= viscosidade dinâmica da água [ML-1T-1].
Nessa equação, a permeabilidade intrínseca ( ) é uma medida da facilidade que a água
tem de atravessar a matriz sólida, quando essa encontra-se saturada. Tecnicamente, é esse
termo que define o valor máximo da condutividade hidráulica. A condutividade relativa ( r) é
21
um número que varia de 0 a 1 e expressa a fração da permeabilidade intrínseca que permanece
quando o solo não encontra-se saturado. Por isso, é quem determina a dependência de K em
relação a . Os dois outros parâmetros, ou seja, a densidade e a viscosidade dinâmica, são
características da água. O efeito da temperatura ocorre sobre os seus valores. Porém, como a
densidade é muito pouco afetada, o efeito principal da temperatura está sobretudo na
viscosidade dinâmica.
A condutividade hidráulica também pode ser afetada pela presença de íons em
solução. Nielsen et al. (1986) descrevem com detalhes esse processo. Geralmente, um
aumento na concentração de íons produz uma redução na condutividade do solo. Esse
processo é parcialmente reversível se a solução for levada pelo fluxo e substituída por outra
menos concentrada. Yong et al. (1996) descrevem os principais efeitos da salinidade sobre a
condutividade hidráulica em solos argilosos. A maioria dos modelos de fluxo de água em
solos, desconsideram esses efeitos. Porém, Nielsen et al. (1986) recomendam que, devido à
sua significância, os efeitos da salinidade devam ser incluídos nos modelos de fluxo e
transporte, embora a dificuldade para isso seja imensa. Este trabalho desconsidera os efeitos
da temperatura e da presença de íons.
2.4- DESCRIÇÃO DO TRANSPORTE DE SOLUTOS EM SOLO NÃO-SATURADO
Os materiais dissolvidos na água do solo são transportados pelo fluxo dessa, porém,
em geral, com uma velocidade diferente. O transporte de solutos em um solo é o resultado de
três processos principalmente: advecção, dispersão e difusão. Além disso, uma substância
dissolvida pode sofrer reações químicas tanto com as partículas da matriz sólida quanto com
outras substâncias ou ainda aderir-se às partículas sólidas, em um processo conhecido como
sorção.
2.4.1- Descrição da Advecção
O transporte advectivo, ou simplesmente advecção, ocorre quando o soluto é
carreado pelo fluxo de água. Portanto, na ausência de um outro processo, ele ocorreria à
mesma velocidade do fluxo. A velocidade real com que o fluido se move por dentro dos poros

da matriz sólida é maior que a velocidade descrita pelo fluxo q . Isso porque nem toda a área
22
transversal está disponível para o fluxo, devido à presença dos grãos sólidos e de bolhas de ar.

A velocidade média real do fluxo relaciona-se com o fluxo q por

q

v
onde
(2.15)

v = velocidade linear média real da água [LT-1].
A massa de soluto transportada por advecção é dependente da sua concentração na
água e do fluxo de água. Pode ser expresso como

Ja
onde

qC

v C
(2.16)

J a = fluxo de massa advectivo de soluto [ML-2T-1] e
C = concentração de massa de soluto na água [ML -3].
A advecção também recebe o nome de convecção. No entanto, o termo advecção é
considerado mais adequado, porque o termo convecção sugere a idéia de movimentação
devida a diferenças de temperatura (Wang e Anderson, 1982).
O resultado do transporte advectivo puro é que, praticamente, não há espalhamento do
soluto na água. Assim, se a entrada do soluto no solo for em forma de um pulso, esse será
transportado como um bloco, não alterando o seu perfil. Essa característica é muito útil em
testes de modelos.
2.4.2- Descrição da Difusão
Normalmente, uma massa de soluto se move de uma região mais concentrada para
uma menos concentrada. Esse movimento é natural e espontâneo, sendo denominado difusão
molecular ou, simplesmente, difusão. Resulta dos movimentos aleatórios das moléculas do
soluto na água e independe da ocorrência de movimentação da água. Basta que exista um
gradiente de concentração. O fluxo de massa difusivo é expresso pela primeira lei de Fick

Jd
Ds
23
*
C
(2.17)
onde

J d = fluxo de massa difusivo de soluto [ML-2T-1],
Ds* = coeficiente de difusão do solo [L2T-1] e
C = gradiente de concentração de soluto na massa de água [ML-4].
O sinal negativo indica que o transporte ocorre na direção decrescente da
concentração.
O coeficiente de difusão está associado com a tortuosidade do meio poroso por (Bear,
1972; Nielsen et al., 1986)
Ds
onde
*
Dm
(2.18)
Dm = coeficiente de difusão molecular do soluto [L2T-1] e
( ) = fator de tortuosidade do solo [adimensional].
A tortuosidade do solo é uma medida do efeito da forma dos poros sobre o fluxo da
água (Fetter, 1993). No caso do solo não-saturado, o conteúdo de água terá efeito significativo
sobre o fator , pois a presença de ar aumenta a tortuosidade do solo. Este trabalho
desconsidera essa variação, sendo a tortuosidade considerada constante.
2.4.3- Descrição da Dispersão
Também chamada de dispersão mecânica, tem o mesmo efeito da difusão, ou seja, o
espalhamento do soluto de uma região de alta concentração para uma de baixa concentração.
Quando a água se move dentro dos poros, a sua velocidade varia em torno da velocidade
linear média. Pode ser maior em poros grandes e menor em poros pequenos. Além disso, no
centro do poro a velocidade é maior que junto às paredes. Assim, como a massa de soluto é
conduzida pela massa de água, essa variação microscópica da velocidade faz com que a massa
de soluto se disperse pelos poros da matriz sólida, tomando rumos diferentes. Por isso, Voss
(1984) considera a dispersão como um processo de pseudo-transporte, já que representa
apenas desvios do processo advectivo, ou seja, em cada poro ocorre um transporte advectivo
que globalmente conduz ao efeito da dispersão. A tortuosidade do solo também é
determinante nesse processo.
Quando o espalhamento da massa de soluto ocorre na direção do fluxo, o processo é
denominado dispersão longitudinal. Como a massa de soluto se dispersa também devido à
24
tortuosidade, há espalhamento em direções normais ao fluxo. Esse é denominado de
dispersão transversal.
A dispersão mecânica pode ser descrita também pela primeira lei de Fick como

Jm
onde

Dd v C
(2.19)

J m = fluxo de massa mecanicamente dispersivo de soluto [ML-2T-1] e

D d v = coeficiente de dispersão mecânica do solo [L2T-1].

Matematicamente, o coeficiente de dispersão D d v é representado por um tensor. A
determinação do valor cada termo do tensor é função da velocidade linear média e
caracterizado por um fator, denominado dispersividade empírica, que possui dois termos:
um representando a dispersão longitudinal e o outro representando a dispersão transversal.
2.4.4- Dispersão Hidrodinâmica
Macroscopicamente, a difusão molecular e a dispersão mecânica produzem o mesmo
efeito: o espalhamento da massa de soluto proporcional à diferença de concentração. Portanto,
os dois processos atuam juntos num fluxo de água. A sua combinação resulta na dispersão
hidrodinâmica.
Assim, a dispersão hidrodinâmica pode ser considerada como a soma dos dois
processos. Daí,

Jh

Jd

Jm
(2.20)
e

Jh
onde
Ds*

Dd v
C
(2.21)

J h = fluxo de massa hidrodinamicamente dispersivo de soluto [ML-2T-1] e
Define-se então o coeficiente de dispersão hidrodinâmica ou, simplesmente,
coeficiente de dispersão ( D ) [L2T-1] como
25
D
Ds*

Dd v
(2.22)
O coeficiente de dispersão ( D ) também é representado matematicamente por um
tensor. A determinação do valor de cada termo desse tensor será discutida mais adiante, no
Item 7.5.
Nielsen et al. (1986) questionam a validade da união dos processos porque, apesar de
macroscopicamente semelhantes, microscopicamente a difusão é um processo ativo que age
em resposta direta a um gradiente de concentração, indiferente ao fluxo, enquanto a dispersão
mecânica é um processo passivo que age somente quando há fluxo. No entanto, Bear (1972)
considera que a separação dos processos é artificial, embora o efeito global da difusão seja
significante somente quando o transporte ocorre em velocidades muito baixas. Este trabalho
considera a atuação conjunta dos processos no transporte dispersivo.
2.4.5- Descrição da Sorção
O soluto transportado pelo fluxo de água no solo está sujeito a inúmeros processos que
tendem a retirar ou introduzir uma certa quantidade de massa desse soluto na água. Alguns
desses processos são normalmente agrupados num processo mais geral, conhecido como
sorção.
A sorção é uma combinação, de pelo menos, quatro processos que atuam
conjuntamente: a adsorção, a absorção, a quimiosorção e a troca iônica. A sorção é,
portanto, um processo geral, onde a massa de soluto dissolvida na fase líquida é atraída pelas
partículas da matriz sólida. O principal efeito da sorção é que o soluto se move mais
lentamente que a água. Esse efeito é conhecido como retardação.
A adsorção ocorre para os solutos da mesma forma que para a água, como descrito no
Item 2.3.1. A força eletrostática resultante da atração entre a superfície carregada das
partículas sólidas e os íons dissolvidos do soluto faz com que ocorra a adesão desses últimos
sobre as primeiras. Como essa força tem significância apenas a uma pequena distância,
somente uma fina película de solutos será adsorvida. Essa força atrativa é, no entanto,
suficientemente fraca para permitir que alguns íons adsorvidos sejam trocados por outros em
solução ou simplesmente retirados pelo fluxo da água. Esse processo é denominado troca
iônica e é um processo que ocorre continuamente na fase líquida. Se os íons trocados tiverem
26
carga positiva, como no caso das argilas que possuem uma carga negativa na superfície das
suas partículas, atraindo e adsorvendo íons positivos, o processo é denominado troca
catiônica. No caso de íons negativos serem atraídos e adsorvidos, o processo é denominado
troca aniônica. A dimensão desses dois processos depende da abundância de íons em solução
e da capacidade de troca iônica das partículas da matriz sólida. Com isso, uma espécie iônica
pode tanto ter sua concentração diminuída ou aumentada na solução, de acordo com o
processo em andamento.
A quimiosorção ocorre através de alguma reação química que faz com que uma parte
da massa de soluto seja incorporada à matriz sólida. A absorção ocorre quando as partículas
da matriz sólida são suficientemente porosas a ponto de permitir que os solutos se difundam
para o seu interior. Esse é o caso de algumas partículas de argila, que possuem espaços entre
os seus minerais.
O resultado global do processo de sorção é a remoção de alguma quantidade de massa
de soluto da água. Se o processo sortivo é relativamente rápido em relação à velocidade do
fluxo, o soluto alcança uma condição de equilíbrio entre as fases líquida e sólida, sendo
denominado de modelo de sorção em equilíbrio. Caso contrário, não há equilíbrio entre as
fases e o processo é denominado de modelo de sorção cinético.
Como normalmente a velocidade de fluxo da água no solo é baixa, os processos de
sorção, num ambiente subterrâneo, conseguem alcançar uma condição de equilíbrio. Assim,
neste estudo serão considerados, para descrever o processo, apenas os modelos de sorção em
equilíbrio.
Devido à sorção ser conseqüência da adesão de soluto sobre os grãos sólidos, a
capacidade de um solo de promovê-la está diretamente associada à área superficial dos seus
grãos. À medida em que ocorre o processo, essa área vai sendo ocupada pela massa de soluto
e o solo vai perdendo o poder de sorver, até que não haja mais áreas livres. Dessa forma, vê-se
que o solo possui uma capacidade limitada para processar a sorção.
2.4.5.1- Modelos de Sorção em Equilíbrio
Matematicamente, a sorção é descrita por uma relação entre a massa de soluto sorvida
sobre a matriz sólida e a massa de soluto em solução, o que, graficamente, determina uma
curva denominada isoterma de sorção. Existem vários modelos de sorção, sendo que três se
27
destacam pela sua aplicabilidade e razoável aproximação ao problema em estudo: sorção
linear, de Freundlich e de Langmuir.
Os modelos de sorção são incorporados ao modelo de transporte baseados na
determinação de um coeficiente geral Kp, denominado coeficiente de partição, através da
relação
C*
t
em que
Kp
C
t
(2.23)
C* = concentração específica de soluto sorvido sobre a matriz sólida [MM -1] e
Kp = coeficiente de partição [L3M-1].
O termo C* é uma relação entre a massa de soluto sorvida e a massa da matriz sólida,
sendo, portanto, adimensional.
• Modelo de sorção linear
Quando a concentração do soluto é baixa (menor que 10-5 Molar ou menor que metade
da solubilidade do soluto (Johnson et al., 1989)), a relação entre o soluto sorvido e o soluto
dissolvido pode ser considerada linear. Essa relação é então descrita pela equação
C*
onde
KdC
(2.24)
Kd = coeficiente de distribuição [L3M-1].
O coeficiente de distribuição (Kd) é um valor constante e depende basicamente do tipo
de solo e do tipo de soluto, ou seja, depende da interação entre a matriz sólida e o soluto.
Para incorporar esse modelo de sorção ao modelo de transporte faz-se
C*
t
Kd
C
t
(2.25)
Assim, para o modelo de sorção linear, o coeficiente geral é Kp = Kd.
Devido à sua simplicidade, a isoterma de sorção linear é a mais utilizada em
simulações de transporte de solutos (Healy, 1990). Fetter (1993), no entanto, adverte que a
28
principal limitação desse modelo é que ele não estabelece um limite máximo para a
quantidade de soluto sorvido. Por isso, o modelo é limitado a baixas concentrações.
• Modelo de sorção de Freundlich
Esse é um modelo de sorção mais geral. É descrito pela equação
C*
onde
K f Cd
(2.26)
Kf = constante de sorção de Freundlich e
d = expoente de Freundlich [adimensional].
As dimensões da constante Kf dependem do valor do expoente d, ou seja, [L3M-1]d.
Novamente, esse modelo é utilizado quando a concentração de soluto é baixa, já que
sofre do mesmo problema da isoterma linear, ou seja, não há um limite para a quantidade de
soluto sorvido. Aliás, observa-se que o modelo linear é apenas um caso especial do modelo de
Freundlich, quando d = 1.
Para incorporar esse modelo ao modelo de transporte, faz-se
C*
t
dK f C d
1
C
t
Assim, o coeficiente de partição para o modelo de Freundlich fica K p
(2.27)
dK f C d-1 .
• Modelo de sorção de Langmuir
A isoterma de sorção de Langmuir baseia-se na idéia de que a matriz sólida possui
uma área finita disponível para a sorção. Com isso, existe um valor máximo para a massa de
soluto possível de ser sorvida (Johnson et al., 1989; Fetter, 1993).
O modelo de Langmuir é descrito pela equação
C*
C* max
29
K lC
1 K lC
(2.28)
onde
C*max = máxima concentração específica de soluto sorvido [adimensional] e
Kl = constante de sorção de Langmuir [L3M-1].
Para incorporar o modelo de Langmuir ao modelo de transporte, faz-se
C*
t
Kl
C* max
1 K lC
C
t
2
Kl
C* max
Assim, o coeficiente de partição fica K p
(2.29)
1 K lC
2
.
• Generalização dos modelos de sorção
Como visto acima, os três modelos podem ser incorporados ao modelo de transporte
de acordo com as Equações (2.23), (2.25), (2.27) e (2.29). Observa-se que essas equações são
semelhantes em termos de forma, apresentando uma ou duas constantes, de modo que podem
ser generalizadas.
Para as isotermas linear e de Freundlich, pode-se fazer então
C*
1C
C*
t
sendo
(2.30a)
2
1 2C
2
1
C
t
1
= Kf, para a isoterma de Freundlich, ou Kd, para a isoterma linear, e
2
= d, para a isoterma de Freundlich, ou 1, para a isoterma linear.
(2.30b)
Para a isoterma de Langmuir, faz-se
C*
C*
t
2C
1
1
(2.31a)
2C
1 2
1
30
2C
2
C
t
(2.31b)
sendo
1
= C*max e
2
= Kl.
Essas duas expressões são utilizadas no modelo de transporte, objeto deste estudo.
2.4.6- Equacionamento do Transporte de Solutos em Solo Não-Saturado
O fluxo de massa total é a combinação dos efeitos da advecção e da dispersão
hidrodinâmica. Pode ser expresso como a soma das Equações (2.16) e (2.21)

J
onde

v C
(2.32)
D C

J = fluxo de massa total de soluto [ML-2T-1].
A massa de soluto total (mt) em um volume de solo é a soma da sua massa dissolvida
mais a massa sorvida sobre a matriz sólida. A massa de soluto dissolvida é o produto do
conteúdo de água ( ) pela concentração de soluto (C) pelo volume total da matriz sólida (V t).
A massa de soluto sorvida é o produto da massa específica do solo (B d) [ML-3] pela massa de
soluto sorvida (C*) pelo volume Vt. Assim, a massa de soluto total num volume de solo é
mt
CVt
B d C* Vt
(2.33)
Da mesma forma que foi feito para o fluxo, a continuidade para o soluto no volume
elementar da Figura 2.9 diz que a massa de soluto que entra menos a que sai é igual à taxa de
armazenamento no volume. A descarga de massa de soluto numa direção i é o produto do
fluxo Ji pela área da face i. Assim, a continuidade será
Jx y z Jy x z Jz x y
Jz
Jz
z
z
Jx
x
x
Jx
C V
t
x y
31
y z
Jy
B d C* V
t
Jy
y
y
x z
(2.34)
Jz + Jz / z z
z
Jx
y
y
x
z
Jy + Jy / y y
Jy
x
Jx + Jx / x x
Jz
Figura 2.9- Volume elementar representativo do solo, mostrando o balanço de massa do
soluto
Reorganizando a equação e considerando que não há mudança de volume da matriz
sólida, a equação da continuidade se reduz a
C
t
B d C*
t
Jy
Jx
x
Jz
z
y
(2.35)
Vetorialmente,
C
t

J
B d C*
t
(2.36)
Substituindo a Equação (2.32) em (2.36), resulta
C
t
B d C*
t

v C
D C
(2.37)
Na equação acima podem ainda ser consideradas fontes e retiradas de massa de soluto.
Daí, obtém-se a equação fundamental do transporte de solutos para solos não-saturados
C
t
onde
i
B d C*
t

v C
D C
= fontes e retiradas de massa de soluto [ML-3T-1].
32
i
(2.38)
Essas fontes e retiradas de soluto são, talvez, os processos mais difíceis de se modelar,
devido às várias causas que podem ter. Decaimento radioativo, precipitação e dissolução
químicas, absorção de solutos por raízes de plantas, utilização e transformação de solutos por
microorganismos etc. são alguns dos processos mais modelados (Nielsen et al., 1986). Devido
a isso, muitas combinações de elementos, tais como as espécies químicas modeladas e as
várias condições de contorno, tornam a resolução da Equação (2.38) muito difícil.
Nesse estudo, esse termo será aproximado por reações de decaimento e produção de
primeira ordem e ordem zero, respectivamente. Assim, tem-se
i
onde
l
l
C
l
*
s Bd C
s Bd
(2.39)
= constante de primeira ordem de decaimento de massa de soluto na fase líquida
[T-1],
s
= constante de primeira ordem de decaimento de massa de soluto na fase sólida
[T-1],
l
= termo de ordem zero de produção de massa de soluto na fase líquida [ML -3T-1]
e
s
= termo de ordem zero de produção de massa de soluto na fase sólida [T -1].
Os valores de
l
e
l
dependem somente do tipo de soluto enquanto que
s
e
s
dependem da interação entre o soluto e a matriz sólida.
Substituindo a equação acima na Equação (2.38) e reorganizando os termos, tem-se
B d C*
t
C
t

v C
D C
l
C
s Bd C
*
l
s Bd
(2.40)
Expandindo os termos do lado esquerdo e a primeira parcela do lado direito da
Equação (2.40) tem-se
C
t
C
t
Bd
C*
t
l
C*
C
Bd
C
t
*
s Bd C
33

v
l

v
s Bd
C
D C
(2.41)
Voss (1984) mostra que o valor de
Bd/ t é desprezível, de modo que pode-se
considerar Bd constante. Também, da continuidade do fluxo de água, Equação (2.10), sabe-se
que

q
t

v
(2.42)
Substituindo e reorganizando a Equação (2.41), chega-se a
C
t
Bd
C*
t

v
C
D C
l
C
s Bd C
*
l
s Bd
(2.43)
Substituindo a Equação (2.23) para incorporar a sorção, reorganizando e dividindo por
, chega-se a
Bd K p

v
C
O termo adimensional
1
1
C
t
1
D C
Bd K p
lC
s Bd C
*
s Bd
l
(2.44)
é conhecido como fator de retardação (R).
Dividindo a equação por R tem-se
C
t

v
R
1
R
C
D C
'
R
'
R
(2.45)
sendo que
'
lC
'
l
s Bd C
s Bd
34
*
(2.46a)
(2.46b)
A Equação (2.45) mostra que a velocidade de movimentação do soluto é menor que a

velocidade do fluxo v , quando o fator de retardação é maior que 1. Então, a velocidade do

soluto é uma velocidade corrigida ( v c ) [LT-1], dada por

v
R

vc
(2.47)
e a equação do transporte chega finalmente a
C
t

vc
C
1
R
D C
'
R
'
R
(2.48)
Essa é a equação final a ser implementada no modelo de transporte de solutos, objeto
deste estudo.
35
3- MODELAGEM
No estudo do transporte de solutos pelas águas do subsolo, o principal objetivo é a
observação das concentrações, a movimentação e o destino desses solutos. Para isso, o
procedimento mais confiável seria a realização de testes de campo para observação do
comportamento do sistema solo-água para uma dada modificação no ambiente, por exemplo,
a introdução de um aterro sanitário. No entanto, esses testes são normalmente extremamente
onerosos e muito lentos, devido à baixa velocidade com que a água se move no solo,
fornecendo, dessa forma, pouca informação sobre o comportamento futuro do sistema. Assim,
em geral, é inviável a condução de testes de campo para a avaliação de situações diversas e a
previsão de seus efeitos. Por isso, a forma mais utilizada para essas verificações é a tentativa
de simular os efeitos de fenômenos naturais com o auxílio de modelos.
Um modelo é uma representação simplificada dos sistemas ou processos reais (Wang e
Anderson, 1982; Keely, 1989). Segundo Mercer e Faust (1980a) e Keely (1989), os modelos
para simulação de água subterrânea podem ser divididos em três categorias: físicos,
analógicos e matemáticos. Cada categoria possui vantagens e desvantagens em relação às
outras.
3.1- MODELOS FÍSICOS
Os modelos físicos são aqueles que representam fisicamente o sistema estudado a
partir da utilização de um protótipo, construído em escala reduzida e podendo ser do mesmo
material do sistema, reproduzindo tanto quanto possível as mesmas condições do caso real.
Possui a grande vantagem de possibilitar um conhecimento intuitivo do comportamento do
sistema. No entanto, exigem um grande espaço e bastante tempo para a geração de um
conjunto de dados significante. Quando o protótipo é constituído do mesmo material do
sistema, exigem também um enorme cuidado na obtenção e na manipulação das amostras de
campo, para que não ocorram distúrbios consideráveis sobre a sua condição natural (Keely,
1989). Além disso, o efeito de escala deve ser sempre considerado na interpretação dos
resultados, pois os fenômenos medidos em escala de laboratório podem ser significantemente
diferentes dos ocorridos nas condições de campo. Normalmente, esse tipo de modelo é
construído para um conjunto de situações específicas e essa falta de flexibilidade o torna
36
pouco utilizado na prática. Para estudo de água subterrânea, os principais modelos físicos são
os tanques de areia e as colunas de solo.
3.2- MODELOS ANALÓGICOS
Os modelos analógicos também representam o sistema fisicamente, mas seu princípio
de funcionamento é a semelhança física ou analogia. Para estudar um fenômeno complexo,
estuda-se o comportamento de um outro mais simples, análogo em termos físicos, ambos
relacionados pelas suas equações básicas. Este tipo de modelo é geralmente construído para
uma situação específica, não permitindo mudanças experimentais. Além disso, a obtenção de
dados é, de certa forma, lenta. Por isso, é pouco utilizado atualmente. Para o estudo de água
subterrânea em solo saturado, o fluxo de água é simulado por circuitos elétricos, numa
analogia entre a lei de Darcy e a lei de Ohm. Para o solo não-saturado, não se usa este tipo de
modelo, devido à dificuldade de se obter sistemas análogos (Campos, 1998).
3.3- MODELOS MATEMÁTICOS
Os modelos matemáticos são modelos não físicos que procuram representar os
fenômenos naturais por meio de equações matemáticas, quantificando as relações entre
parâmetros específicos e as variáveis em estudo. Normalmente, consistem de um conjunto de
equações diferenciais com apropriadas condições iniciais e de contorno que descrevem o
processo estudado. Quando o processo é muito complexo para ser simulado com exatidão,
algumas hipóteses simplificadoras podem ser adotadas. A utilidade de um modelo matemático
para determinada situação depende da quantidade e das características das simplificações
adotadas. No estudo de água subterrânea, esses são os modelos mais aplicados atualmente.
A sua principal vantagem é poder simular diferentes cenários rapidamente e ser
flexível, possibilitando modificações, tornando possível o estudo de diferentes problemas
locais com o mesmo modelo. Isso permite que sejam feitas previsões mais rápidas, objetivo
da modelagem. Outra vantagem desses modelos é que podem facilitar a compreensão
funcional entre causa e efeito dos processos naturais sem, no entanto, fornecer uma visão
intuitiva do processo (Keely, 1989).
Suas principais desvantagens são a necessidade de uma enorme quantidade de dados
de boa qualidade, muitas vezes descrevendo a história do sistema, e a dificuldade de
37
formulação teórica das equações de alguns processos. Além disso, ficam limitados, em termos
de cálculos, à precisão do simulador utilizado, que na maioria dos casos é um computador.
Os modelos matemáticos podem ser subdivididos em determinísticos ou estocásticos,
podendo ser analíticos ou numéricos. Podem, ainda, ser classificados como de parâmetros
distribuídos ou agregados e como conceituais ou empíricos. Cada categoria será brevemente
discutida aqui. O relatório EPA/625/R-94/001 (EPA, 1994) e Ambrose et al. (1996) fornecem
detalhes dos diversos tipos de modelos já desenvolvidos, bem como a forma de aquisição dos
mesmos.
3.3.1- Modelos Determinísticos e Estocásticos
Modelos determinísticos são os que trabalham com relações causa-efeito exatas, isto
é, definem uma resposta única a uma determinada excitação. Normalmente, assumem que o
sistema atua de forma que um dado conjunto de eventos conduz a uma situação única,
indiferente às incertezas contidas nos dados. Os modelos estocásticos atuam com base na
incerteza inerente aos parâmetros e aos processos atuantes no sistema. Por isso, não fornecem
uma resposta única mas um conjunto de respostas possíveis. Esses modelos calculam a
probabilidade, dentro de um intervalo de confiança desejado, do valor da variável em um
determinado ponto.
Os modelos determinísticos são muito usados nos problemas de previsão, pois
permitem que uma maior classe de problemas seja considerada com um mesmo modelo. Os
modelos estocásticos são muito úteis na classificação de dados e na descrição de sistemas
fracamente compreendidos (Mercer e Faust, 1980a).
3.3.2- Modelos de Parâmetros Distribuídos e Agregados
As características espaciais dos parâmetros de um modelo podem ser distribuídas ou
agregadas. Os modelos de parâmetros distribuídos funcionam com base na caracterização
espacial do sistema, de modo que, em cada ponto ou área, os parâmetros assumem valores
espeíficos. Os modelos de parâmetros agregados podem ser usados quando o sistema total
pode ser representado por características em um ponto, não levando em conta a variabilidade
espacial do sistema.
38
3.3.3- Modelos Conceituais e Empíricos
O tipo de equação ou função utilizada por um modelo matemático pode defini-lo como
conceitual ou empírico. Os modelos conceituais, também conhecidos como modelos de base
física, são aqueles cujas funções utilizadas são baseadas na física do sistema, isto é, os
processos físicos são conhecidos e representados por equações que os descrevem. Os modelos
empíricos utilizam funções empíricas para representar o sistema, procurando ajustar os
valores calculados aos dados observados, sem necessidade de conhecimento da física do
sistema.
Considerando que os modelos conceituais em grande parte utilizam-se de funções
empíricas, mas que estão relacionadas com os processos físicos, alguns autores ainda
subdividem os modelos conceituais em físicos e semiconceituais. Os semiconceituais
caracterizam os processos mas os parâmetros das equações utilizadas são definidos
basicamente de forma empírica, como no caso da equação de Richards (valores de K( ))
(Tucci, 1998).
3.3.4- Modelos Analíticos e Numéricos
Os modelos analíticos são aqueles que usam soluções exatas das equações
governantes do sistema. Estas soluções são geralmente apresentadas como expressões
matemáticas de forma fechada, com as variáveis contínuas no espaço e no tempo. No entanto,
são resultado de um grande número de simplificações, incluindo a geometria do sistema e
suas condições de contorno. Por isso, só podem ser aplicadas a situações muito restritas e,
normalmente, são usadas na verificação dos modelos numéricos.
Os modelos numéricos apresentam soluções aproximadas para as equações
governantes. Essas soluções são pontuais, ou seja, as variáveis não são contínuas no espaço e
no tempo, mas discretizadas sobre a região do sistema. No entanto, as hipóteses
simplificadoras são em número muito menor e, por isso, menos restritivas, permitindo a
aplicação do modelo a várias situações.
Existem ainda modelos semi-analíticos que empregam técnicas numéricas para
aproximar soluções analíticas complexas. Com isso, produzem valores aproximados não
contínuos das variáveis procuradas. Estes modelos apresentam uma maior flexibilidade na
39
incorporação de situações diversas, se comparado com modelos analíticos, sem aumento
significativo da necessidade de dados (EPA, 1994).
3.4- IMPORTÂNCIA DOS MODELOS MATEMÁTICOS
Uma das principais funções de um modelo é fornecer previsões sobre o estado e a
condição de um sistema no futuro. Com isso, os modelos podem ser usados para o teste de
hipóteses e verificação de seus efeitos. Às vezes, modelos preditivos podem ser combinados
com rotinas de otimização, podendo ser utilizados no gerenciamento de recursos (Keely,
1989). Outra função dos modelos é a determinação de esquemas de monitoramento,
especialmente em recursos hídricos. Com resultados de previsões, pontos de monitoramento
podem ser locados em posições estratégicas, fornecendo informações importantes para os
trabalhos de campo. Outro uso importante é a identificação de parâmetros, também chamado
de problema inverso ao problema de previsão. Nesse caso, os modelos são usados para
estimar características importantes do sistema estudado.
Quando se aplica um modelo a uma dada situação, deve-se levar em conta a validade
dos resultados obtidos. Essa validade dependerá de quanto o modelo se aproxima do
comportamento global real do sistema. Para isso, é necessário um conjunto de dados de boa
qualidade, com boa distribuição espacial. Devido à incerteza inerente a esses dados, previsões
absolutas são praticamente impossíveis de serem obtidas. Segundo O’Connell e Todini (1996)
e Sun (1996), a quantificação das incertezas nos dados, na estrutura dos modelos e,
conseqüentemente, nas previsões, constitui um importante problema na modelagem numérica.
A confiança nos resultados de um modelo depende, em boa parte, de três propriedades
do mesmo: a consistência, a convergência e a estabilidade. A consistência é conseguida
quando o modelo representa de forma discreta e com boa aproximação a solução exata das
equações do sistema contínuo. No geral, uma melhor consistência é obtida com o aumento do
número de pontos discretos no sistema. A convergência é conseguida quando, diminuindo-se
os intervalos de tempo e espaço, os erros na aproximação diminuem e a solução do modelo
tende a alcançar a solução exata da equação do sistema. A estabilidade está relacionada
basicamente com a característica temporal do modelo, sendo, portanto, própria de modelos
transientes. É conseguida quando os erros produzidos pelo modelo diminuem á medida em
que o processo de cálculo avança no tempo e, portanto, qualquer perturbação tende a
desaparecer com a evolução do processo.
40
Assim, os modelos numéricos podem ser ferramentas poderosas quando usadas
apropriadamente. As principais situações possíveis podem ser analisadas e os estados
resultantes podem ser observados, permitindo que, no caso de gerenciamento, a melhor
decisão possa ser tomada.
41
4- A MODELAGEM MATEMÁTICA
Segundo definição adotada no relatório EPA/625/R-94/001 (EPA, 1994), um modelo
matemático pode ser considerado como um conjunto de equações que expressam o sistema
físico e incluem hipóteses simplificadoras ou ainda a representação de um sistema físico por
meio de expressões matemáticas com as quais o comportamento do sistema pode ser deduzido
com conhecida precisão.
Normalmente, um modelo matemático para a hidrologia subterrânea, além das
equações de descrição do sistema (também chamadas equações governantes), inclui também
especificações em relação ao domínio do problema. Essas especificações definem a geometria
da região do sistema, as condições do contorno dessa região e, para os problemas transientes,
as condições iniciais.
No caso do transporte de solutos, as condições de contorno podem ser de três tipos:
a) de Dirichlet, onde os valores da concentração do soluto no sistema são conhecidos;
b) de Neumann, onde a taxa de entrada de massa de soluto é conhecida, sendo um valor
constante ou uma função do tempo; e
c) de Cauchy, onde a taxa de entrada do soluto é dependente do valor da sua concentração.
Um detalhamento matemático desses termos pode ser encontrado em Mercer e Faust
(1980b) e Sun (1996).
O tratamento matemático dessas condições de contorno é relativamente simples. No
entanto, o tratamento físico necessita de uma atenção especial. As condições de campo para
uma simulação necessitam de um conjunto de observações considerável. Mesmo assim,
podem conter erros de modo a conduzir a uma interpretação errônea da realidade. Além disso,
os dados podem não descrever a situação completamente. Assim, as condições de contorno no
campo podem ser apenas aproximadas e a confiabilidade disso depende muito da experiência
de quem está realizando o estudo, que deve fazer o julgamento adequado.
As condições iniciais devem descrever a situação do sistema no início da simulação,
em relação à variável estudada. No caso do transporte de solutos, devem definir, para cada
ponto, o valor da concentração do soluto. Essas condições são usadas somente em problemas
transientes, isto é, que variam com o tempo. O seu tratamento físico também é complicado
pelo fato de que os dados de campo devem ser suficientemente abrangentes para fornecer uma
cobertura total da região do sistema.
42
Os modelos matemáticos têm sido desenvolvidos há muito tempo. Inicialmente,
devido à dificuldade dos cálculos, os modelos eram praticamente todos analíticos e, por isso,
muito simplificados e restritos. Após a década de 60, com o avanço tecnológico dos
computadores digitais e das técnicas de programação, os modelos numéricos passaram a ser
amplamente utilizados (Wang e Anderson, 1982).
Existem atualmente numerosos métodos numéricos para a resolução da equação
governante do transporte de solutos. Eles podem ser divididos em três grandes categorias:
Eulerianos, Lagrangeanos e Lagrangeano-Eulerianos. Os métodos Eulerianos procuram
resolver a equação governante para um conjunto de pontos fixos sobre um sistema de
referência fixo no espaço. São mais simples no seu princípio e no tratamento da equação
governante. Os métodos Lagrangeanos tomam como referência o movimento natural da massa
do sistema, utilizando uma grade de pontos deformável ou uma grade fixa de pontos sobre um
sistema de coordenadas móveis, que acompanha o movimento do sistema. Os métodos mistos
Lagrangeano-Eulerianos utilizam uma grade fixa juntamente com um conjunto de pontos
móveis para a resolução da equação governante, sendo os cálculos realizados em duas fases.
Na primeira, é utilizada uma visão Lagrangeana do processo para a resolução da advecção,
promovendo a movimentação de pontos ou partículas que representam a massa do sistema. A
segunda fase é utilizada para verificar o posicionamento desses pontos, ligá-los à grade fixa e
resolver a dispersão por algum método Euleriano (Cheng et al., 1984; Yeh, 1990; Sun, 1996).
Segundo o relatório EPA/625/R-94/001 (EPA, 1994), as técnicas de solução numéricas
mais utilizadas na hidrologia subterrânea para o transporte de solutos são:
1) método das diferenças finitas;
2) método das diferenças finitas integradas;
3) método dos elementos finitos, de Galerkin e variacional;
4) método da colocação;
5) método dos elementos de contorno;
6) métodos de acompanhamento de partículas (particle tracking); e
7) método das características.
Os métodos das diferenças finitas e dos elementos finitos são os mais utilizados na
prática. Podem ser caracterizados como métodos Eulerianos, pois utilizam uma grade de
pontos fixa sobre um sistema de coordenadas também fixo. Atualmente, os métodos
Lagrangeanos e mistos Lagrangeano-Eulerianos têm sido freqüentemente utilizados,
apresentando um grande avanço em termos teóricos, como, por exemplo, o método das
43
características. A discussão a seguir enfoca os aspectos básicos de alguns dos métodos
mencionados acima.
4.1- O MÉTODO DAS DIFERENÇAS FINITAS
O método das diferenças finitas é o resultado da aplicação da aproximação por séries
de Taylor a uma função. Esse método é normalmente usado quando as equações governantes
são equações diferenciais de, no máximo, segunda ordem. A série de Taylor consiste em
expandir uma função em uma série infinita de parcelas para a determinação do valor
aproximado dessa função na vizinhança de um ponto de valor conhecido.
Assim, considerando x0 como o ponto conhecido e x0+ x o ponto para o qual se deseja
conhecer o valor da função f(x) e realizando a expansão de Taylor tem-se
f x0
x
f x0
x
x 2 d 2f x 0
2!
dx 2
df x 0
dx
x 3 d 3f x 0
3!
dx 3

(4.1)
Considerando os termos de segunda ordem e superior como desprezíveis e rearranjando
convenientemente a Equação (4.1) tem-se
df x 0
dx
f x0
x f x0
x
(4.2)
donde se tira uma expressão aproximada para a derivada de primeira ordem. Essa expressão é
conhecida como diferença progressiva.
Da mesma forma, se x0- x é o ponto da vizinhança para o qual deseja-se conhecer o
valor da função f(x), a série de Taylor fica
f x0
x
f x0
x 2 d 2f x 0
2! dx 2
df x 0
x
dx
x 3 d 3f x 0
3!
dx 3

Novamente, desprezam-se os termos de segunda ordem e superior e chega-se a
44
(4.3)
df x 0
dx
f x0
f x0
x
x
(4.4)
Essa expressão é conhecida como diferença regressiva.
A subtração da Equação (4.1) pela (4.3) resulta em
f x0
x
f x0
x
2 x
df x 0
dx
2
x 3 d 3f x 0
3!
dx 3
2
x 5 d 5f x 0
5!
dx 5

(4.5)
onde desconsideram-se os termos de terceira ordem e superior e chega-se a
df x 0
dx
f x0
x f x0
2 x
x
(4.6)
expressão conhecida como diferença central.
Somando-se as Equações (4.1) e (4.3) tem-se
f x0
x
f x0
x
2f x 0
2
x 2 d 2f x 0
2! dx 2
2
x 4 d 4f x 0
4! dx 4

(4.7)
e, desconsiderando-se os termos de quarta ordem e superior, chega-se a
d 2f x 0
dx
f x0
x
2
2f x 0
x
f x0
x
(4.8)
2
com a qual obtém-se uma aproximação para a segunda derivada da função.
As derivadas cruzadas podem ser aproximadas por (Sun, 1996)
2
f x 0 , y0
f x0
x, y 0
y f x0
x, y 0
x y
4 x y
f x0
x, y 0
y f x0
x, y 0
y
4 x y
45
y
(4.9)
No procedimento acima e da teoria da aproximação por séries de Taylor, as expressões
são tão mais precisas quanto menor for o espaçamento
x e, também, quanto mais linear a
função. Essas expressões podem ser obtidas para qualquer variação na direção e no tempo.
Para a aplicação do método das diferenças finitas, deve-se considerar a região do
sistema como um conjunto de pontos discretos, formando uma malha. Com isso, uma grade é
definida. Existem dois tipos de grade, mostrados na Figura 4.1. Na grade de malha centrada,
os pontos do sistema, chamados de pontos nodais, são locados na interseção das linhas de
grade, enquanto que, na malha de bloco centrado, os pontos nodais são locados nos
centróides da região delimitada pelas linhas de grade. Na prática, não há grande diferença
entre o uso das duas, exceto em casos específicos, que dependem das condições de contorno
(Faust e Mercer, 1980a).
pontos
nodais
y
y
x
x
(a) malha centrada
(b) bloco centrado
Figura 4.1- Malhas de elementos finitos ilustrando os tipos de grade
A grande maioria dos modelos que utilizam o método das diferenças finitas trabalham
com grades retangulares de espaçamentos regulares, como as da Figura 4.1. No entanto, essa
não é uma condição necessária, muito embora facilite a aplicação do método. Grades com
espaçamentos irregulares podem ser usadas em problemas que possuam áreas que necessitem
de um maior detalhamento, tais como as vizinhanças da superfície freática e proximidades de
poços.
Os pontos nodais, ou nós, da malha são geralmente designados por um par ordenado
de números inteiros, na forma genérica (i, j). A Figura 4.2 ilustra parte de uma grade de malha
centrada.
46
y
(i-1, j+1)
(i, j+1)
(i+1, j+1)
(i-1, j)
(i, j)
(i+1, j)
(i-1, j-1)
(i, j-1)
(i+1, j-1)
x
Figura 4.2- Grade de malha centrada ilustrando a numeração dos nós
A título de exemplo de aplicação, considere a Equação (2.37) no plano vertical com

e v homogêneos e constantes, com D homogêneo e anisotrópico e sem o termo de massa
sorvida C*. Com isso, a equação se reduz a
C
t
onde
C
x
vx
vz
C
z
2
2
C
x2
D xx
2
C
x z
2D xz
C
z2
D zz
(4.10)
vm = componente do vetor velocidade na direção m [LT-1] e
Dmn = componente do tensor do coeficiente de dispersão na direção mn [L2T-1].
Para a grade da Figura 4.2, pode-se aplicar, à Equação (4.10), diferença progressiva
para a derivada no tempo e diferença central para as derivadas de primeira ordem no espaço.
Obtém-se, então,
C im, j 1
D xx
C im1, j
C im, j
t
2C im, j
x
vm
x ,i , j
C im1, j
2
D zz
C im1, j
C im1, j
2 x
C im1, j
2D xz
C im, j
1
1
vm
z ,i , j
C im, j
C im1, j
1
1
2 z
C im1, j
4 x z
2C im, j
z
C im, j
C im, j
1
1
C im1, j
1
(4.11)
1
2
onde o índice m indica a posição da variável no tempo. No esquema acima vê-se que somente
uma concentração é calculada para um passo de tempo futuro m+1 enquanto todas as outras
47
são conhecidas no passo atual m. Existe, portanto, somente uma incógnita e o esquema é
considerado na forma explícita. A vantagem desse esquema é a facilidade de programação,
mas sua solução só é estável para incrementos de tempo suficientemente pequenos (Wang e
Anderson, 1982; Sun, 1996).
Se todas as concentrações do lado direito da Equação (4.11) forem tomadas no passo
de tempo m+1, haverá então várias incógnitas para cada passo de tempo e o esquema assim
formado é considerado na forma implícita. Assim, para cada passo de tempo, ocorre a
montagem de um sistema de equações para o cálculo das concentrações. Esse esquema é mais
estável que o anterior, de mais difícil solução.
Pode-se melhorar a aproximação usando um esquema misto, que empregue uma
mistura dos esquemas implícito e explícito. Para isso, as derivadas de espaço são avaliadas em
um passo de tempo intermediário, isto é, entre m e m+1. Isso é feito utilizando-se um
parâmetro de ponderação para cada um dos esquemas, por exemplo, com
para o esquema
explícito, no passo de tempo m, e 1- para o implícito, no passo de tempo m+1. Esse esquema
também gera um sistema de equações que deve ser resolvido a cada passo de tempo. Se for
tomado
= 1/2, a aproximação é feita na metade do intervalo de tempo e o esquema é
conhecido como método de Crank-Nicholson. A grande vantagem dos esquemas implícitos
e mistos sobre o explícito é a sua estabilidade incondicional e um número menor de iterações
necessário para a resolução do sistema de equações gerado.
4.2- O MÉTODO DOS ELEMENTOS FINITOS
O método dos elementos finitos consiste, basicamente, na transformação da equação
diferencial governante num conjunto de equações integrais. Para a sua aplicação, a área do
domínio do sistema é dividida em áreas menores, denominadas elementos, que discretizam o
sistema. A forma e o posicionamento desses elementos são definidos por um conjunto de
pontos nodais, para os quais serão encontrados os valores das variáveis procuradas. Esses
pontos nodais podem ser localizados nos vértices ou nas laterais desses elementos. A posição
dos vértices e a área dos elementos não precisam ter nenhuma distribuição rígida, de modo
que o método é bastante flexível na representação da geometria do sistema. No entanto,
quanto maior o número de nós e elementos, maior será a precisão dos cálculos.
A forma dos elementos e o número de nós de cada elemento podem ser variados de
acordo com a formulação do método. Normalmente, os mais empregados são os elementos
48
triangulares com três nós ou quadrangulares com quatro nós. Porém, podem ser utilizados
elementos triangulares com seis nós, quadrangulares com oito nós etc. (Cirilo e Cabral, 1989).
Para trabalhos no plano bidimensional, os elementos triangulares são preferidos, por
possibilitarem uma melhor representação da geometria do sistema e um maior refinamento da
malha e por possuírem formulação mais simples. A Figura 4.3 ilustra uma grade de nós para o
método dos elementos finitos.
contorno
elemento
Figura 4.3- Grade de elementos finitos com elementos triangulares
O primeiro passo na formulação do método é definir uma solução aproximada para a
variável incógnita do problema. No caso do transporte de solutos, a variável desejada é a
concentração C. Essa solução é, normalmente, uma combinação linear de C, expressa como o
somatório do produto do valor da concentração em cada nó e de uma função de base
associada a cada nó, na forma
NN
Ĉ x, z, t
i
x, z C i t
(4.12)
i 1
onde
Ĉ (x, z, t) = solução aproximada de C(x, z, t),
i(x,
z) = função de base associada ao nó i,
Ci(t) = valor da concentração no nó i no tempo t e
NN = número total de nós.
Uma exigência imposta pelo método é que as funções de base sejam contínuas e
tenham derivadas contínuas em todo o domínio. Normalmente, também devem satisfazer a
todas as condições de contorno do problema.
As funções de base ( i) podem ser de qualquer ordem (linear, quadrática, cúbica etc.) e
devem ser linearmente independentes. De um modo geral, quanto maior for a ordem de
49
i,
melhor será a aproximação Ĉ (Cirilo e Cabral, 1989). As funções de base lineares são as mais
usadas em formulações com elementos triangulares, por serem mais facilmente integradas e
diferenciadas (Faust e Mercer, 1980a) e a discussão a seguir enfoca, justamente, a aplicação
do método para esses elementos.
Como o domínio é subdividido em elementos, é usual que a aproximação seja feita
sobre eles, ou seja,
Ĉ e x, z, t
3
e
(4.13)
x, z C i t
i
i 1
Utilizando-se funções de base lineares, a Equação (4.13) conduz à formação de um
plano, que representa o comportamento da variável no interior do elemento, interligando os
seus valores nos nós, como mostrado na Figura 4.4.
C
z
Cj
Ck
Ci
x
Figura 4.4- Plano descrito por uma função de base linear para um elemento triangular
O segundo passo é a obtenção de uma representação integral da equação governante.
Os dois métodos mais utilizados para isso são o método variacional e o método dos resíduos
ponderados. A principal diferença entre os dois é que, no método dos resíduos ponderados,
trabalha-se diretamente com a equação diferencial enquanto que, no método variacional, usase um funcional relacionado à equação diferencial e às condições de contorno (Faust e
Mercer, 1980a). O método variacional não será discutido aqui. Wang e Anderson (1982)
abordam o método.
50
4.2.1- Método dos Resíduos Ponderados
Nesse método, a equação governante é designada por um operador diferencial O(u),
sendo u a variável do problema. No caso da Equação (4.10), o operador diferencial seria
OC
C
t
vx
C
x
vz
2
C
z
C
x2
D xx
2
2D xz
C
x z
2
D zz
C
z2
0
(4.14)
Como Ĉ é um valor aproximado, a aplicação do operador diferencial resultará num
valor não nulo, denominado de resíduo.
O Ĉ
Re síduo
0
(4.15)
Naturalmente, quanto melhor for a aproximação Ĉ , menor será o resíduo.
A idéia do método é diminuir os resíduos tanto quanto possível. Para isso, os erros são
distribuídos dentro da região do domínio, o que é feito ponderando-se os resíduos sobre cada
nó do sistema. Daí então, impõe-se que o erro médio ao longo de todo o domínio seja nulo,
fazendo
W
onde
i
O Ĉ
i dW
0
(4.16)
= função de ponderação relativa ao nó i e
W = área do domínio (região do sistema).
A Equação (4.16) é uma expressão do erro médio, ponderado pelas funções de
ponderação ( i). Essas funções devem ser contínuas e linearmente independentes, de modo a
gerar um sistema de equações algébricas passível de solução.
Existem vários métodos de resíduos ponderados, cada um relacionado com o tipo de
função de ponderação adotado, dentre os quais merece destaque especial, pela sua
aplicabilidade, o método de Galerkin. Outros métodos podem ser encontrados em Wrobel
(1989) e Sun (1996).
51
4.2.2- O Método de Galerkin e as Funções de Base
Esse método pode ser considerado um caso especial do método dos resíduos
ponderados, no qual as funções de ponderação são as próprias funções de base, isto é,
i
(4.17)
i
O método de Galerkin é o mais comum na aplicação de elementos finitos. Também
muito comum é a utilização de funções de base lineares no método.
As funções de base são usualmente relacionadas com as áreas dos elementos e as
coordenadas dos seus nós. Por isso, são definidas separadamente para cada um deles. No caso
de funções de base lineares, em geral, adotam-se as seguintes características dentro do seu
elemento:
a)
i
é igual a 1 no nó i e 0 nos outros nós;
b)
i
varia linearmente com a distância ao longo dos lados; e
c)
i
é igual a zero ao longo do lado oposto ao nó i, isto é, do lado não conectado ao nó i.
No caso de elementos triangulares, outra característica importante é que, no centróide
do elemento,
i
possui valor de 1/3. A Figura 4.5 mostra o seu comportamento para um
elemento triangular.
z
z
1
1
i
i
x
x
(a) para um elemento conectado ao nó i
(b) para vários elementos conectados ao nó i
Figura 4.5- Comportamento da função de base linear para um elemento triangular ilustrando o
plano formado
52
4.3- COMPARAÇÃO ENTRE OS MÉTODOS DAS DIFERENÇAS FINITAS E DOS ELEMENTOS
FINITOS
De acordo com o teorema fundamental do cálculo, a diferenciação e a integração são
operações essencialmente inversas. Por isso, os dois métodos descritos anteriormente podem
ser considerados como equivalentes. O método das diferenças finitas aproxima a solução da
equação diferencial governante ao usar uma diferenciação aproximada pela série de Taylor
enquanto que o método dos elementos finitos resolve a equação diferencial governante usando
uma aproximação integral (Faust e Mercer, 1980a; EPA, 1994).
Vários estudos já foram feitos avaliando os prós e contras de cada método. Faust e
Mercer (1980a) apresentam uma breve discussão dessa comparação, que é sumarizada na
Tabela 4.1.
Tabela 4.1- Vantagens e desvantagens dos métodos das diferenças finitas e dos elementos
finitos
Vantagens
Desvantagens
MÉTODO DAS DIFERENÇAS FINITAS
Base intuitiva
Baixa precisão para alguns problemas
Fácil entrada de dados
Grades regulares
Eficientes técnicas de matriz
Fáceis mudanças no programa
MÉTODO DOS ELEMENTOS FINITOS
Geometria flexível
Base matemática avançada
Melhor precisão
Difícil entrada de dados
Avalia melhor termos de produto cruzado
Programação difícil
Fonte: Faust e Mercer (1980a), adaptado
Uma desvantagem comumente associada a esses dois métodos é a ocorrência de
problemas numéricos de oscilação e dispersão quando o problema inclui a passagem de
frentes de concentração de solutos, ou seja, grandes variações no valor da concentração num
curto intervalo de tempo. Isso ocorre principalmente se o problema for predominantemente
advectivo. Para tentar contornar esses problemas, outros métodos foram propostos.
53
4.4- OUTROS MÉTODOS APLICÁVEIS À SIMULAÇÃO DO TRANSPORTE DE SOLUTOS
Atualmente, a maior parte dos métodos numéricos aplicados para a simulação do
transporte de solutos enquadram-se na categoria dos métodos mistos Lagrangeano-Eulerianos.
Esses métodos têm como característica principal o isolamento dos processos de transporte,
simulando a advecção separadamente, o que implica num melhor tratamento dos problemas
numéricos de oscilação e dispersão. Dentre esses métodos, merecem destaque especial os
métodos das Características e do Caminhamento Aleatório.
4.4.1- O Método das Características
O método das características é uma técnica bastante utilizada na resolução da equação
do transporte de solutos para evitar os problemas numéricos de dispersão e oscilação. É um
tipo de método de acompanhamento de partículas e o seu conceito é razoavelmente simples.
Um conjunto de partículas de fluído é transportada pelo campo de fluxo e mudanças na
concentração do soluto nessas partículas é observada (Konikow e Bredehoeft, 1978).
Considere a Equação (4.10) rearranjada para
C
t
vx
C
x
vz
C
z
2
D xx
C
x2
2
2D xz
C
x z
2
D zz
C
z2
(4.18)
No lado esquerdo dessa equação estão os dois termos advectivos, com v x e vz, e no
lado direito estão os termos dispersivos.
Se for observado o movimento de uma partícula ao longo da sua linha de fluxo, a sua
posição num tempo t qualquer é dada por suas coordenadas xp(t) e zp(t) e a sua velocidade tem
componentes
dx
dt
vx
(4.19a)
dz
dt
vz
(4.19b)
54
Uma curva que satisfaça às Equações (4.19) é denominada característica, que é uma
linha de fluxo. Para uma partícula que percorre uma característica, a derivada total da
concentração em relação ao tempo é dada por
dC
dt
C
t
C dx
x dt
C dz
z dt
C
t
vx
C
x
vz
C
z
(4.20)
Com isso, a Equação (4.18) pode ser rescrita como
dC
dt
2
D xx
C
x2
2
2D xz
C
x z
2
D zz
C
z2
(4.21)
Essa equação representa o transporte de solutos ao longo de uma característica. Notase aqui que os termos de advecção não aparecem e, com isso, a resolução da equação pode ser
realizada pelas técnicas usuais dos elementos finitos e das diferenças finitas.
A Equação (4.21) pode ser explicada da seguinte forma: quando se acompanha o
movimento de uma partícula ao longo de uma característica, isto é, usando a partícula como
ponto de referência, qualquer mudança da concentração da partícula é resultado somente do
processo de dispersão. A sua resolução deve ser feita simultaneamente à das Equações (4.19),
isto é, no mesmo passo de tempo. Assim, o método consiste de duas partes básicas: uma para
o processo de advecção e outra para o de dispersão. No processo de advecção, a resolução das
Equações (4.19) é feita pela movimentação de partículas sobre as características. De um modo
geral, essas partículas são geradas sobre todo o domínio ou sobre uma região de interesse e
movidas de acordo com o campo de velocidades, que descrevem as características. A
movimentação é acompanhada sobre uma malha de pontos nodais fixos. Cada partícula
carrega uma concentração de soluto, relativa à concentração de um ponto nodal. Após
finalizada a movimentação, as concentrações dos pontos nodais mudam devido às novas
posições das partículas, sendo determinadas de acordo com a quantidade e as concentrações
das partículas próximas ao ponto. Finalizada a advecção, o processo de dispersão é finalmente
simulado pela resolução da Equação (4.21) por algum método de malha fixa, como o método
das diferenças finitas. Detalhes do método das Características podem ser encontrados em
Konikow e Bredehoeft (1978) e Sun (1996).
Apesar de eliminar os erros provocados pela dispersão e oscilação numéricas, o
método das características tem a grande desvantagem de exigir uma enorme quantidade de
55
memória computacional, em especial para um grande número de partículas, e, sob certas
circunstâncias, produzir elevados erros no balanço de massa (Zheng, 1993).
4.4.2- O Método de Caminhamento Aleatório (Random Walk)
O método de caminhamento aleatório é um método de acompanhamento de partículas,
semelhante ao método das características, que utiliza o movimento de partículas para
descrever o transporte de solutos no solo. Porém, aqui o transporte dispersivo é simulado
como um movimento aleatório das partículas, de modo que a nova posição é determinada sem
a necessidade de resolver a equação de dispersão, o que facilita bastante os cálculos (Prickett
et al., 1981).
A base do método é a característica estatística do fenômeno da dispersão, que produz
um espalhamento do soluto como uma distribuição normal ao redor de uma posição média.
No método, cada partícula de fluido se move ao longo de sua linha de fluxo, isto é, a sua
característica, simulando o processo de advecção, numa distância igual a v t. O movimento
aleatório das partículas é então calculado como duas parcelas de dispersão, uma longitudinal e
uma transversal, ambas relacionadas com o deslocamento advectivo. Com isso, a nova
posição de uma partícula é numericamente igual à soma da posição anterior, do deslocamento
provocado pela advecção, do deslocamento provocado pela dispersão longitudinal e do
deslocamento provocado pela dispersão transversal.
O método elimina os problemas de dispersão e oscilação numéricas, porém os seus
resultados mostram somente qual a concentração mais provável e não a absoluta, de modo que
é possível que duas simulações de um mesmo problema conduzam a resultados próximos, não
iguais. Além disso, a sua maior desvantagem é a necessidade de um enorme número de
partículas para garantir resultados mais precisos, o que aumenta a área de armazenamento de
dados computacionais.
Prickett et al. (1981) desenvolveram o método inicialmente para o caso
unidimensional. Walton (1991) e Sun (1996) fornecem mais detalhes sobre o assunto,
apresentando a formulação para o caso bidimensional.
56
5- ESTUDOS RECENTES EM MODELAGEM NUMÉRICA DO TRANSPORTE DE
CONTAMINANTES EM SOLOS
A modelagem numérica do fluxo e do transporte de contaminantes na subsuperfície é
um estudo relativamente recente, dentro da hidrologia. Somente nas três últimas décadas é
que sofreu um avanço considerável, acompanhando o avanço dos computadores digitais, que
permitiram, em especial, a evolução da modelagem numérica. A situação é mais complicada
ainda quando o estudo se refere à zona não-saturada do solo. Até a década de 70, os estudos
dentro dessa zona limitavam-se aos primeiros metros de profundidade do solo, abrangendo
somente a zona das raízes. Tais estudos eram realizados com finalidades principalmente
agrícolas. Somente nos últimos 20 anos, os hidrólogos passaram a preocupar-se com a zona
não-saturada como um todo, devido à sua grande importância em problemas de contaminação
de aqüíferos (Yeh et al., 1993), e vários modelos foram desenvolvidos, entre os quais SUTRA
(Voss, 1984), CHEMFLO (Nofziger et al., 1989), VS2DT (Healy, 1990) e 3DLEWASTE
(Yeh et al., 1992). No entanto, poucos desses modelos simulam o transporte bi e
tridimensional (SUTRA, VS2DT, 3DLEWASTE).
Na aplicação de modelos para o transporte bidimensional, Faust e Mercer (1980b)
apontam dois processos que representam grandes problemas para a modelagem: a dispersão e
a sorção. A dispersão surge como um problema devido à sua característica de tensor. Em
geral, os tensores funcionam com produtos cruzados (D xz e Dzx) que exigem um pouco mais
de esforço no tratamento matemático das equações. Com isso, a aproximação por diferenças
finitas convencional não trabalha muito bem com problemas que envolvem dispersão,
principalmente devido à fraca ligação diagonal entre os nós adjacentes. A utilização de um
modelo modificado de diferenças finitas, que inclua tais ligações, pode resolver esse
problema, muito embora exija um maior esforço no tratamento matemático (Faust e Mercer,
1980a). Por essa razão, poucos trabalhos têm sido desenvolvidos nesta linha. Como no
método dos elementos finitos a ligação diagonal é muito comum, este problema é
minimizado.
O outro grande problema é a sorção. Este processo, apesar de já ter sido muito
estudado, é pouco conhecido com relação aos mecanismos atuantes. Além disso, existem
vários tipos de contaminantes e cada um desses possui um comportamento diferente no
subsolo. No caso de solo parcialmente saturado, outra complicação que surge é o fato de que,
devido à presença de bolhas de ar, nem toda a superfície dos grãos da matriz sólida é capaz de
57
realizar a sorção. Com isso, ela se torna um processo dependente do conteúdo de água, o que
se reflete no coeficiente de partição (Kp), da Equação (2.23). Para manter esse coeficiente
inalterado, Yeh e Ward (1979), apud Faust e Mercer (1980b), sugerem a substituição da
massa específica do solo (Bd) por uma massa específica “efetiva” que é obtida simplesmente
dividindo o valor de Bd por
no fator de retardação. Essa dependência é normalmente
desprezada nos estudos, por acreditar-se ter pouca influência global, sendo o coeficiente de
partição considerado independente de .
Na resolução da equação de advecção-dispersão, surge um grande problema que é a
ocorrência de problemas numéricos de dispersão e oscilação, próximo a frentes de passagem
de contaminantes. Uma frente é uma grande variação de valores da variável em um pequeno
intervalo de espaço. A dispersão numérica provoca um espalhamento indesejado dessa frente,
alterando a precisão dos cálculos, e a oscilação numérica produz uma seqüência de valores
que oscilam em torno do valor correto da frente. Esses efeitos podem ser visualizados na
Figura 5.1. Sun (1996) explica que a dispersão numérica ocorre mais para problemas
dominados por advecção, quando a velocidade do fluxo é relativamente grande. Isso é
refletido pelo número de Peclet (Pe), que é dado por
Pe
onde
v x
DL
(5.1)
Pe = número de Peclet [adimensional],
v=
vx
2
vz
2
= módulo da velocidade linear média [LT -1] e
DL = coeficiente de dispersão longitudinal [L2T-1].
Quando este número é alto, o problema é basicamente advectivo e, quando é baixo, o
problema passa a ser basicamente dispersivo.
Sun (1996) analisa o efeito da dispersão numérica no método das diferenças finitas,
mostrando que é provocado pelo truncamento da expansão de Taylor. Yeh et al. (1993)
fornecem uma visão mais acurada matematicamente, afirmando que, quando o problema se
torna predominantemente advectivo, a equação do transporte passa de parabólica para quase
hiperbólica e que os métodos numéricos comuns geralmente produzem, por isso, dispersão
numérica e oscilação perto de frentes.
Segundo Rao e Hathaway (1989), a dispersão numérica tem efeito similar à dispersão
física. Por isso, apresentam um modelo de transporte utilizando o conceito de célula de
58
mistura onde a dispersão física é substituída pela dispersão numérica (também denominada de
dispersão artificial). Nesse modelo, o sistema é dividido em células e cada uma atua como um
reator de mistura completa, baseado no princípio de conservação de massa. O modelo atua
relativamente bem quando o problema não é predominantemente dispersivo.
Oscilação numérica
C
Espalhamento provocado
pela dispersão numérica
x
Solução exata
Solução numérica
Figura 5.1- Efeitos da oscilação e dispersão numéricas em uma frente de passagem de
contaminante
Yu e Singh (1995) apontam que a maioria dos modelos de transporte de contaminantes
baseados no método de Galerkin estão sujeitos à oscilação e dispersão numéricas devido a um
erro na continuidade do fluxo de contaminantes nos contornos dominados por advecção, pois
a formulação usual não aplica o teorema de Green para os termos advectivos da equação.
Além disso, durante a aplicação dos modelos, outros fatores são preponderantemente
causadores desses efeitos. Tais fatores incluem principalmente uma má escolha da malha de
elementos finitos e do tamanho do passo de tempo adotado. Por isso, sugerem cinco
modificações na formulação do método. Os efeitos dessas modificações foram testados e
mostraram resultados melhores que os do método convencional, em relação aos problemas
numéricos de dispersão e oscilação.
Hossain e Yonge (1997) estudaram o método dos elementos finitos por três técnicas
diferentes e procuraram avaliar os efeitos de dispersão e oscilação. Trabalhando com a
equação do transporte para um meio saturado, as técnicas utilizadas foram o método de
Galerkin comum, o método de Galerkin-Petrov e o método de Galerkin característico. O
59
método de Galerkin-Petrov apresentou resultados livres de oscilações, mas com maior
dispersão numérica em todas as simulações realizadas. Os métodos de Galerkin comum e
característico apresentaram dispersão e oscilação numéricas semelhantes entre si, mas, em
relação ao método anterior, a dispersão numérica foi menor e as oscilações maiores, quando
as simulações eram advectivas. No entanto, quando o coeficiente de dispersão usado era
grande, maior que um dado valor, os resultados destes dois métodos foram praticamente livres
de oscilações.
As formulações baseadas nos métodos de acompanhamento de partículas são as que
apresentam melhores resultados em termos de diminuição da dispersão numérica. Por isso, é
comum a combinação de métodos usuais com os de acompanhamento de partículas.
Yeh et al. (1993) combinaram o método dos elementos finitos de Galerkin com o
método das características modificado para tentar obter resultados numéricos livres destes
erros numéricos. O método das características modificado foi utilizado para resolver a parte
advectiva da equação de transporte enquanto que o método dos elementos finitos foi utilizado
para a parte dispersiva. Comparando os resultados do modelo desenvolvido com soluções
analíticas, foi verificado que o modelo apresentou bons resultados e que a oscilação numérica
praticamente foi eliminada, embora o problema de dispersão numérica ainda tenha se
mantido, mesmo em baixo nível. Além disso, erros no balanço de massa foram encontrados,
atingindo um máximo de cerca de 5% para uma dada simulação. Essa é uma formulação
interessante que pode incorporar alguns outros elementos de modo a eliminar a dispersão e o
erro no balanço de massa.
Li (1996) desenvolveu um modelo para solo não-saturado unidimensional a partir de
uma formulação que denominou “esquema de características de quarta ordem de seis pontos”
(6P4O), baseado no método das diferenças finitas e no método das características com o uso
de uma combinação ótima de três polinômios cúbicos para a solução da parte advectiva. Essa
formulação, apesar de ter boa precisão, provoca problemas de oscilação numérica. Tais
problemas foram resolvidos usando um limitador de fluxo. O modelo apresentou, por isso,
resultados livres de oscilações mas com alguma dispersão numérica.
Outros métodos de acompanhamento de partículas foram desenvolvidos recentemente.
Ewen (1996b) desenvolveu um modelo para o transporte de contaminantes em solos nãosaturados baseado numa nova formulação, descrita em Ewen (1996a). O modelo, denominado
SAMP (Subsystems and Moving Packets), simula o movimento unidimensional de água em
uma coluna vertical de solo. A massa de água é considerada como pequenos pacotes que
60
podem conter alguma concentração de contaminantes. A formulação do modelo é bastante
complexa, mas pode servir de base para o desenvolvimento de uma formulação eficaz.
Algumas modificações dos métodos usuais são encontradas na literatura. Delay et al.
(1993) modificaram o método de caminhamento aleatório basicamente trocando a forma de
gerenciar as partículas móveis. No método usual, as partículas são tratadas individualmente
enquanto que Delay et al. (1993) propuseram tratá-las na forma de variáveis inteiras
representando a massa total de contaminante contida em cada célula, num espaço discretizado
unidimensional. Posteriormente, Delay et al. (1996) estenderam o modelo para o caso
bidimensional.
Outra modificação do método de caminhamento aleatório é apresentada por Banton et
al. (1997). Considerando um sistema unidimensional, a modificação consiste na mudança do
domínio do problema. Enquanto os métodos usuais trabalham com o domínio espacial, o
método proposto por Banton et al. (1997) trabalha com o domínio temporal, de modo que as
derivadas espaciais da equação de transporte são transformadas para derivadas no tempo,
avaliando o avanço da concentração através do tempo de um ponto a outro.
Zheng (1993) propõe algumas mudanças no método das características para tentar
solucionar alguns problemas inerentes ao método, citados no Item 4.4.1. As mudanças
consistem basicamente em alocar dinamicamente as partículas, a fim de diminuir a
necessidade de memória computacional, e em um esquema de interpolação das velocidades
eficiente, a fim de reduzir os erros no balanço de massa. Com as mudanças, o método é
estendido para o caso tridimensional.
Considerando que a dispersão numérica pode ser reduzida com o emprego de uma
malha mais detalhada, com espaçamentos menores entre os pontos nodais, e que essa
dispersão tem como causa principal as frentes de contaminantes, Yeh (1990) propôs um
método de acompanhamento de partículas que utiliza uma técnica de adaptação automática da
malha. Em cada elemento, existe um certo número de “nós escondidos”. Quando a frente é
detectada sobre algum elemento, os seus nós escondidos aparecem e se tornam nós comuns,
subdividindo o elemento e detalhando melhor a malha naquela região. Após a passagem da
frente, os nós voltam a se “esconder”. Com isso, não é necessário utilizar uma malha
detalhada em todo o domínio, pois o detalhamento é automático e ocorre somente nas regiões
de interesse. O método, denominado por Yeh (1990) como LEZOOM (Lagrangian-Eulerian
Method with Zoomable Hidden Fine-Grid), reduz significativamente a necessidade de
memória computacional e praticamente elimina a dispersão e a oscilação numéricas.
61
Poucos estudos de transporte têm sido realizados com técnicas de diferenças finitas.
No entanto, Amokrane e Villeneuve (1996) desenvolveram um procedimento de diferenças
finitas, baseado no método de MacCormack, para a resolução da equação de fluxo de água em
solo não-saturado unidimensional, que pode ser estendido para a resolução da equação de
transporte. O método de MacCormack é um esquema explícito previsor-corretor que utiliza
alternância entre diferenças progressivas e regressivas. No estudo de Amokrane e Villeneuve
(1996), no passo previsor, é usada diferença progressiva para o cálculo das derivadas de
primeira ordem no espaço e, no passo corretor, é usada diferença regressiva. As derivadas de
segunda ordem são aproximadas normalmente. Apesar de não levar em conta a histerese da
condutividade hidráulica, o método tem precisão de segunda ordem e estabilidade
condicionada.
Hayworth e Burris (1996) também utilizam o método de diferenças finitas, porém,
com esquema de Crank-Nicholson, para estudar o transporte e o comportamento de
contaminantes num solo saturado homogêneo unidimensional. O modelo considera o processo
de sorção para atenuar a concentração de um surfactante catiônico. Os resultados apresentados
pelo modelo não foram comparados com resultados analíticos pela falta destes, porém
mostram algumas vantagens no uso de alguns processos de sorção em relação a outros.
5.1- A ESCOLHA DO MÉTODO A SER EMPREGADO
O modelo de transporte de solutos a ser implementado neste estudo será desenvolvido
com base num modelo de simulação do fluxo em meio não-saturado já existente. Tal modelo
foi desenvolvido por Koide (1990) e modificado por Campos (1998). O modelo utiliza o
método dos elementos finitos com elementos triangulares para a discretização do domínio e
diferenças finitas para o tempo. O trabalho de Campos (1998) analisou a eficiência de alguns
métodos de resolução da equação de fluxo, apresentando uma discussão sobre as
potencialidades e as deficiências de cada um. A melhor combinação dos métodos, para cada
caso a ser modelado, é apresentada. Apesar de terem sido feitos poucos testes (Campos,
1998), o modelo se mostrou bastante eficiente no cálculo do fluxo, o que encorajou a sua
utilização neste estudo.
Para a escolha do método a empregar na implementação do modelo de transporte,
objeto deste estudo, o fato do modelo de fluxo ser em elementos finitos foi levado em
consideração. Devido a isto, o uso do método das diferenças finitas ficou, portanto,
62
prejudicado, pois haveria a necessidade de forçar que a malha de pontos do sistema fosse
retangular ou então de realizar um processo de interpolação dos valores das variáveis e
parâmetros relevantes à solução da equação de transporte. O método dos elementos finitos,
neste caso, torna-se uma continuação do processo de modelagem, pois já permitiria o uso da
malha utilizada para a solução do fluxo, sem maiores problemas. No entanto, tanto no método
das diferenças finitas como no método dos elementos finitos, ocorrem os problemas de
dispersão e oscilação numéricos. Os métodos de acompanhamento de partículas, tais como os
métodos das Características e de Caminhamento Aleatório, minimizam esses problemas e
podem ser acoplados a malhas diversas.
A solução encontrada foi então utilizar um misto de método dos elementos finitos e
algum método de acompanhamento de partículas. Com isso, uma proposta interessante é a
utilização do método das características em conjunção com o método dos elementos finitos,
baseado principalmente nas idéias apresentadas nos trabalhos de Yeh et al. (1993) e Zheng
(1993), seguindo as formulações apresentadas no Item 4.
63
6- FORMULAÇÃO DO MODELO
A equação de transporte de solutos nesse modelo, Equação (2.48), será resolvida com
a combinação de dois métodos: o método das características e o método dos elementos
finitos. Para isso, o processo geral de transporte é dividido em dois processos: a advecção e a
dispersão, que inclui também os processos de produção, decaimento e sorção de solutos. O
método das características é utilizado para a resolução do processo de advecção. O método
dos elementos finitos é utilizado para a resolução da dispersão. Neste estudo, são utilizados
elementos triangulares com funções de base lineares. A formulação adotada tem por base os
trabalhos apresentados por Yeh et al. (1993) e Zheng (1993).
6.1- APLICAÇÃO DO MÉTODO DAS CARACTERÍSTICAS
Rearranjando-se a Equação (2.48), tem-se:
C
t

vc C
1
R
D C
R
R
(6.1)
Como foi mostrado no Item 4.4.1, o lado esquerdo dessa equação representa a
derivada total da concentração de soluto em relação ao tempo, ao longo de uma linha
característica, conforme a Equação (4.21). Com isso, a equação do transporte é rescrita como
dC
dt
1
R
D C
R
R
(6.2)
que agora representa o processo de dispersão e produção/decaimento de solutos ao longo de
uma característica, que é uma linha de fluxo.
A resolução da Equação (6.1) é então realizada em duas partes: uma simulando o
processo advectivo e a outra simulando o processo dispersivo. Isso é feito com a resolução

simultânea das Equações (4.19), aplicadas a v c , que simulam a advecção, e da Equação (6.2),
que simula a dispersão.
Para a aplicação do Método das Características na simulação do processo advectivo,
de um modo geral, várias partículas são geradas e distribuídas sobre o domínio e movidas de
64

acordo com o campo de velocidades corrigidas ( v c ). Essa movimentação segue as linhas
características e é acompanhada sobre uma malha de pontos nodais fixos, a malha de
elementos finitos, resolvendo, dessa forma, as Equações (4.19). Com isso, as partículas são
conduzidas de um tempo k para um tempo k*. Nesse processo, cada partícula carrega uma
determinada concentração, relacionada com a concentração Cik de um dos nós. Quando todas
as partículas são movimentadas e as suas posições finais conhecidas, as suas concentrações
são repassadas para os nós, determinando, assim, as concentrações C ik* dos nós, finalizando a
advecção.
O passo de tempo k* não é um tempo físico, pois representa um passo de tempo
intermediário entre os passos k e k+1, relacionado com o fim do processo de advecção e o
início do processo de dispersão. Dessa forma, o tempo k* torna-se um limite numérico entre
os processos e a advecção é simulada antes da dispersão.
Com a advecção finalizada, a dispersão pode então ser simulada, o que é feito pela
resolução da Equação (6.2). Para isso, a derivada total da Equação (6.2), que é definida sobre
as características, é passada para os nós da malha de elementos finitos, aplicando, para cada
nó i, uma aproximação na forma de diferenças finitas
dC i
dt
onde
Ci
t
Ci
k 1
Ci
k*
t
(6.3)
Cik+1 = concentração do nó i num passo de tempo k+1 e
Cik* = concentração do nó i num passo de tempo k*.
A resolução da Equação (6.2) é então resolvida pelo método dos elementos finitos,
simulando o processo de dispersão.
6.1.1- Geração e Distribuição das Partículas Sobre o Domínio
O primeiro passo a ser tomado no método das características é a geração e a
distribuição de partículas sobre o domínio. Essa distribuição pode ser feita de duas formas:
uniformemente sobre o domínio ou uniformemente sobre o elemento.
A primeira forma gera uma densidade uniforme de partículas sobre o domínio. Sua
principal vantagem é a criação de um padrão de resolução, em termos de partículas, sobre a
região simulada. Porém, regiões de maior interesse dentro do domínio, que necessitem um
65
maior detalhamento, como as frentes de contaminação e as proximidades do lençol freático,
onde normalmente existe um número maior de nós, passam a ter o mesmo padrão de
resolução que regiões de menor interesse. Já a segunda forma gera uma distribuição irregular
de partículas sobre o domínio, mas, no entanto, produz um padrão geométrico de distribuição
em torno dos nós, que faz com que regiões de maior interesse, que possuem um número maior
de nós e elementos, possuam, também, um número maior de partículas. No entanto, as duas
formas não resolvem o problema do detalhamento das frentes.
Assim, para promover um melhor detalhamento da malha, o modelo desenvolvido
neste trabalho incorpora a distribuição uniforme sobre o elemento.
Para o modelo, também foram implementados dois padrões de distribuição: com 3 ou
6 partículas por elemento, podendo o usuário utilizar um deles. Com 3 partículas, cada uma é
associada a um nó e posicionada sobre a mediana que parte daquele nó, a 1/3 do seu
comprimento, como mostra a Figura 6.1.a. Com 6 partículas, cada nó passa a ter associadas
duas partículas. O elemento é subdividido por suas medianas em seis subelementos e cada
partícula é colocada no centróide de um desses subelementos, como mostrado na Figura 6.1.b.
É importante notar que, para elementos muito pequenos, não haverá grande diferença, em
termos físicos, entre os dois padrões, podendo haver um aumento desnecessário do número de
partículas para a simulação.
b) com 6 partículas (2 por nó)
a) com 3 partículas (1 por nó)
Figura 6.1- Padrões de distribuição de partículas
Como as partículas são utilizadas para o carregamento de concentrações, pode ser
interessante gerar partículas somente sobre os nós contaminados, isto é, que já tenham alguma
massa de soluto, ou que sejam fonte de massa de soluto, o que reduz bastante o número de
partículas necessárias à simulação. Por isso, o modelo também incorpora a possibilidade de
geração de partículas sobre todo o domínio ou somente sobre os nós contaminados, também a
critério do usuário.
66
6.1.2- Movimentação e Acompanhamento de Partículas
Cada partícula no domínio está associada a três valores – as suas coordenadas xp(t) e
zp(t) e sua concentração Cp(t) – e a duas entidades do domínio discretizado – o nó e o
elemento. O acompanhamento das partículas é feito pela determinação dessas cinco
características a cada passo de tempo.
A movimentação de cada partícula depende da velocidade do fluxo no ponto onde ela
se encontra. Sendo conhecida a posição da partícula num passo de tempo k, o seu
deslocamento e, conseqüentemente, a sua posição no passo k+1 são proporcionais ao

intervalo de tempo t entre esses passos e à velocidade v c da posição onde a partícula se
encontra inicialmente. Assim, a movimentação é simulada como um movimento retilíneo
uniforme em cada passo, isto é, as coordenadas da nova posição são calculadas por uma
aproximação por diferenças finitas das Equações (4.19), na forma:
xp
k 1
zpk
em que
k
v xc,p t
(6.4a)
zpk
v zc,p t
(6.4b)
xp
1
(xpk, zpk) e (xpk+1, zpk+1) = coordenadas da partícula p nos passos de tempo k e k+1,
respectivamente, e
vxc,p e vzc,p = componentes da velocidade corrigida na coordenada (xpk, zpk).

A velocidade v c,p da partícula é calculada por uma interpolação linear das velocidades
dos nós, como a função de aproximação por elementos descrita no Item 4.2. Por essa

interpolação, os valores da velocidade v c,p em cada nó formarão um plano dentro do
elemento, da mesma forma que aquele apresentado na Figura 4.4, e as velocidades na posição
da partícula serão dadas por
3
v xc,p
e
( x p k , z p k ) v xc,i ( t )
(6.5a)
e
( x p k , z p k ) v zc,i ( t )
(6.5b)
i
i 1
3
v zc,p
i
i 1
67
Antes da movimentação, cada partícula está associada a um nó e a um elemento. O
elemento associado é aquele em cuja área a partícula se encontra e o nó associado é o mais
próximo da partícula, dentre os três nós do elemento. A concentração que a partícula carrega
está relacionada com a concentração desse nó. Após a movimentação, com a mudança de
posição, essas associações são refeitas e a concentração da partícula é passada para o novo nó
associado.
A concentração Cik* para cada nó i é a média aritmética das concentrações de todas as
partículas associadas ao nó. Dessa forma, a advecção é resolvida e a dispersão pode então ser
simulada.
A Figura 6.2 apresenta os procedimentos básicos na aplicação do Método das
Características.
6.2- APLICAÇÃO DO MÉTODO DOS ELEMENTOS FINITOS
Neste modelo, são utilizados elementos triangulares com funções de base lineares. Os
vértices dos triângulos são os nós da malha. Assim, para a aplicação do método dos elementos
finitos, uma solução aproximada para a concentração, na forma da Equação (4.13).
Da mesma forma, também é adotada uma aproximação similar para os componentes
do tensor de dispersão
D̂ mn e ( x, z, t )
3
e
i
(6.6)
( x, z) D mn,i ( t )
i 1
onde
m e n podem representar tanto x quanto z,
D̂ mn (x, z, t ) = solução aproximada de Dmn(x, z, t) e
Dmn,j(t) = componente do tensor de dispersão do nó j no tempo t.
A Equação (6.2) pode ser convenientemente rearranjada, fornecendo o operador
diferencial na forma
O(C)
1
D C
68
'
' R
dC
dt
0
(6.7)
tempo t = 0
tempo t = 0
a) no início da simulação, várias partículas são
geradas e distribuídas sobre o domínio
b) as partículas são associadas aos nós e aos
elementos
tempo t = k
tempo t = k
c) num tempo k qualquer, as concentrações Ck,
conhecidas nos nós, são passadas para as partículas
d) determinado o campo de velocidades, as
características são conhecidas
tempo t = k*
tempo t = k*
e) as partículas são movidas, seguindo as
características, e as suas posições finais, no tempo
k*, são conhecidas
f) as associações são refeitas
tempo t = k*
g) as concentrações das partículas são passadas para
os nós associados e as concentrações Ck* são
calculadas para os nós (média aritmética)
Figura 6.2- Procedimento básico do Método das Características
69
Pelo método dos resíduos ponderados, descrito no Item 4.2.1, esse operador, quando
aplicado às aproximações de C e D mn, gera um resíduo que é ponderado sobre cada nó no
domínio, na forma integral da Equação (4.16). Utilizando a técnica de Galerkin, essa integral
fica
1
D̂ Ĉ
W
'
dĈ
dt
' R
j dW
0
(6.8)
para j = 1 até NN.
A Equação (6.8) pode ser escrita na forma de diferenciais
1
W
Ĉ
x
D̂ xx
x
Ĉ
z
D̂ xz
z
D̂ zx
Ĉ
x
D̂ zz
Ĉ
z
'
' R
dĈ
dt
j dW
0 (6.9)
Além disso, como o domínio é discretizado em elementos e em função da forma como
as funções de base são definidas, a Equação (6.9) pode ser integrada sobre as áreas desses
elementos. Com isso, a integração total é substituída por uma soma de integrações
1
x
e
e
D̂ xx
Ĉ e
x
e
D̂ xz
'
Ĉ e
z
e
D̂ zx
z
dĈ e
' R
dt
e
j
dxdz
e
Ĉ e
x
D̂ zz
e
Ĉ e
z
(6.10)
0
Integrando-se as derivadas de segunda ordem por partes, desenvolvendo e
reorganizando, chega-se a
e
e
D̂ xx
e
Ĉ e
x
D̂ xz
e
Ĉ e
z
e
j
x
z
D̂ zx
e
Ĉ e
x
D̂ zz
Ĉ e
z
e
e
j
x
dxdz
(6.11)
e
R
dĈ
dt
e
e
j

qb
dxdz
e
70
e
j
d
e
'
'
e
j
dxdz
onde

q b = fluxo dispersivo de soluto no contorno
[ML-2T-1].
Neste estudo, somente será considerado fluxo advectivo através do contorno, de modo

que q b será assumido como nulo.
Substituindo (4.13) e (6.6) em (6.11), gera-se um sistema de equações lineares que
pode ser colocado sob a forma matricial
G C
M
dC
dt
N
(6.12)
onde os termos das matrizes e vetores são dados por:
e
e
3
i
G ij
e
k 1
R
e
Nj
'
e
k
dxdz
z
e
e
i
e
'
j
e
e
e
k
z
e
D xz,k
e
k 1
k
j
dxdz
x
e
3
i
dxdz
e
3
i
e
D zx ,k
x
e
k 1
3
i
M ij
D xx,k
x
e
e
D zz,k
k 1
e
e
k
dxdz
(6.13a)
j
z
dxdz
(6.13b)
e
(6.13c)
j
dxdz
A matriz [M] da Equação (6.13b) pode ser diagonalizada se, ao invés de se substituir a
aproximação Ĉ e na Equação (6.11), for feita a substituição do valor de Ci por nó. Com isso,
tem-se que
M jj
R
e
e
e
j
dxdz
(6.14)
Essa forma de apresentação da matriz [M] é conhecida como matriz aglutinada e, quando se
trata de fluxo de água, o seu emprego resulta numa convergência mais rápida, como apontado
por Neuman (1973) e Campos (1998).
71
6.2.1- Funções de Base
Neste modelo, são utilizadas funções de base lineares. As equações dessas funções são
descritas com maiores detalhes em Zienkiewicz (1977), Wang e Anderson (1982) e Cirilo e
Cabral (1989).
No elemento triangular da Figura 6.3, os nós são numerados globalmente, em relação
ao domínio, por r, s e t. Localmente, dentro do elemento, são numerados por 1, 2 e 3. Assim,
para um elemento, as funções de base podem ser definidas por
e
i
com i = 1, 2, 3 e (x, z)
1
ai
2A e
( x, z)
(6.15)
bi x ci z
Ae e onde
a1
x 2z3
x 3z 2
b1
z2
z3
c1
x3
x2
(6.16a)
a2
x 3 z1
x 1z 3
b2
z3
z1
c2
x1
x3
(6.16b)
a3
x 1z 2
x 2 z1
b3
z1
z2
c3
x2
x1
(6.16c)
A área do elemento é dada por
1
b1c 2
2
Ae
(6.17)
b 2 c1
r
1
e
z
3
t
2
s
x
Figura 6.3- Elemento típico
Para essas funções, obtém-se as seguintes derivadas:
72
e
e
bi
i
x
2A
e
e
i
ci
z
2A e
(6.18)
Também, pode-se obter a forma geral para integração:
e l
e
e m
i
j
e n
k
dxdz
2A e
l!m!n!
l m n 2!
(6.19)
que define para as integrais das Equações (6.13) e (6.14),
e
i
e
dxdz
Ae
3
(6.20)
Assim, os termos das matrizes e vetores do sistema de equações podem ser rescritos
como
G ij
e
M ij
1
b i D xx
4A e
R
Ae
3
'
'
e
Nj
e
onde
c i D xz b j
b i D zx
c i D zz c j
(6.21a)
(6.21b)
Ae
3
(6.21c)
D mn = média dos valores de Dmn calculados nos nós do elemento.
6.3- INTEGRAÇÃO NO TEMPO
De acordo com as Equações (6.2) e (6.3), a variação temporal da concentração devida
à dispersão pode ser calculada fazendo
dC i
dt
Ci
t
73
(6.22)
onde
Ci = variação na concentração no nó i devida à dispersão.
Segundo Konikow e Bredehoeft (1978), devido aos processos de advecção e dispersão
ocorrerem simultaneamente, sendo simulados por métodos distintos, e devido às mudanças de
concentração ocorridas em um nó-fonte serem proporcionais ao gradiente de concentração
entre o fluido da fonte e o nó no domínio, o valor de Ci para um passo de tempo k+1 deve
ser calculado como um valor médio da dispersão de C nos passos de tempo k e k*. Com isso,
o sistema de equações passa a ser colocado na forma
1
G
2
k
C
k
G
k*
C
1
M
2
k*
k
M
C
t
k*
k 1
1
N
2
k
N
k*
(6.23)
gerando um esquema explícito para o cálculo de Cik+1
M
k
M
k*
C
k 1
t
G
k
C
k
k*
G
C
k*
N
k
N
k*
(6.24)
Com o valor de Cik+1 calculado, a concentração no passo de tempo k+1 para cada nó
é finalmente calculada pela aplicação de Equação (6.3)
Ci k
1
C i k*
Ci k
1
(6.25)
Agora, a variação Cik+1 deve ser passada também para as partículas para refletir sobre
essas as mudanças ocorridas devido à dispersão. Se essa variação for positiva, ela é
diretamente acrescentada à partícula, isto é,
Cp
k 1
Cp
k
Ci
k 1
(6.26)
Por outro lado, se for negativa, é acrescentada à partícula como um valor percentual em
relação ao valor de Cik*, ou seja,
Cp k
1
Cpk 1
74
Ci k
C i k*
1
(6.27)
Isso é feito com o intuito de evitar que partículas com baixa concentração passem a carregar
valores negativos de C (Konikow e Bredehoeft, 1978; Zheng, 1993).
Note-se que, nesta formulação, os valores de , R e ' foram considerados uniformes
dentro de cada elemento e que, apesar de R e ' dependerem direta ou indiretamente do valor
de C, são tomados como valores médios dos valores nos nós. Os erros numéricos decorrentes
dessas aproximações não foram verificados, mas acredita-se que tenham pouca influência
sobre o resultado global.
6.4- DETERMINAÇÃO DAS VELOCIDADES DE ESCOAMENTO
Para realizar a simulação do transporte de solutos, é necessário o conhecimento das
velocidades de fluxo da água sobre todo o domínio. É esse campo de velocidades que
determina a extensão e a direção do movimentos dos solutos na água. A velocidade é o
principal mecanismo atuante no processo de advecção e tem importância fundamental no
processo de dispersão, já que o valor do coeficiente de dispersão depende basicamente da sua
magnitude. Com isso, a precisão do modelo de transporte é diretamente afetada pela precisão
nos cálculos da velocidade.
Como a região do fluxo é discretizada em elementos finitos, a velocidade poderia ser
calculada com base na declividade do plano formado pelo potencial mátrico calculado no
elemento. Dessa forma simples, cada elemento possuirá uma velocidade única em todos os
pontos de sua área. No entanto, o campo de velocidades resultante dessa aproximação sofre
uma descontinuidade nos pontos nodais e nas fronteiras entre elementos, já que podem
pertencer a mais de um elemento, o que leva a uma violação na conservação da massa (Yeh,
1981). Para evitar esse problema, uma formulação baseada no método dos elementos finitos é
adotada (Yeh, 1981; Yeh et al., 1993).
De acordo com a lei de fluxo de Buckingham, Equação (2.7), o fluxo de água, ou
descarga específica, na direção do eixo x, é dado por:
qx
Kx
x
(6.28)
Esse fluxo é, no entanto, a velocidade aparente ou macroscópica, equivalente à
velocidade de Darcy. A velocidade real no poro, Equação (2.15), para o eixo x, é dada por:
75
qx
vx
qx
vx
Daí,
vx
Kx
(6.29)
x
Para a aplicação do método dos elementos finitos, adotam-se as aproximações
3
e
v̂ x ( x, z, t )
( x, z) v x ,i ( t )
(6.30a)
( x , z)
(6.30b)
i
i 1
3
ˆ ( x, z, t )
e
m
m (t)
m 1
3
e
K̂ x ( x, z, t )
n
(6.30c)
( x, z) K x ,n ( t )
n 1
3
ˆ ( x, z, t )
e
k
( x, z)
(6.30d)
k (t)
k 1
Aplicando a técnica de Galerkin e integrando por elementos, tem-se
3
e
e
e
i
i 1
3
v x ,i
e
m
e
m
j
3
dxdz
m 1
e
e
n
e
3
K x ,n
n 1
x
e
k
e
k
j
dxdz
(6.31)
k 1
que resulta num sistema de equações lineares que matricialmente é escrito como
Ax vx
Bx
(6.32)
onde os elementos das matrizes são dados por
3
A ij x
e
m 1
e
m e i
e
m
76
e
j
dxdz
(6.33a)
Bj
e
3
x
3
k
k
e
x
k 1
e
K x ,n
e
k 1
n
e
j
dxdz
(6.33b)
Os valores das integrais são dados pela Equação (6.19).
Na direção z, o fluxo de água é dado por
qz
z
Kz
Kz
z
Kz
z
(6.34)
e a velocidade no poro
vz
Kz
(6.35)
Kz
z
Similarmente, a equação final para a direção z fica
Az vz
Bz
Cz
(6.36)
dxdz
(6.37a)
onde
3
A ij z
e
Bj
m 1
e
e
e
k 1
K z ,n
z
n 1
3
C jz
K z ,n
e
n 1
j
3
k
k
e
m
3
z
e
m e i
e
e
n
e
j
e
e
n
e
j
dxdz
dxdz
(6.37b)
(6.37c)
Como a velocidade de escoamento depende fortemente dos valores do potencial
mátrico calculados pelo modelo de fluxo, os seus resultados estão intimamente ligados à
precisão do cálculo dos potenciais. Então, para que o campo de velocidades gerado seja
coerente, é necessário que o campo de potenciais seja coerente. Com isso, qualquer falha no
modelo de fluxo reflete-se automaticamente no modelo de transporte.
77
6.5- DETERMINAÇÃO DO COEFICIENTE DE DISPERSÃO
Segundo Sun (1996), muitos estudos estatísticos foram feitos tentando descrever
características e valores para o coeficiente de dispersão mecânica da Equação (2.19), porém,
no geral, são incapazes de dar uma formulação ou uma explicação generalizada para esse
coeficiente. Dessa forma, para se trabalhar com as equações de transporte de solutos, modelos
teóricos tem sido utilizados.
Bear (1972) estabelece que os componentes do coeficiente de dispersão mecânica,
utilizando soma implícita, podem ser expressos por
D d ,ij
onde
ijkl
vk vl
v
(6.38)
i, j, k e l podem representar tanto x quanto z,
Dd,ij = componente do tensor do coeficiente de dispersão mecânica na direção ij
[L2T-1] e
ijkl
= dispersividade (geométrica) do meio [L].
A dispersividade do meio (
ijkl)
representa a influência da geometria e da distribuição
dos poros sobre o transporte dispersivo. Portanto, está diretamente relacionada à tortuosidade
do meio (Bear, 1972), o que implica que, para um solo não-saturado, sofre influência do
conteúdo de água ( ). Neste estudo, porém, essa influência será desconsiderada, sendo que a
dispersividade terá valor constante.
A dispersividade é um tensor de quarta ordem e, para um domínio bidimensional,
possui 16 componentes, dos quais apenas aqueles com quatro índices iguais ou com dois
pares de índices iguais são não nulos (Bear, 1972). Por isso, somente 8 desses componentes
são não nulos.
Para um meio isotrópico, os componentes não nulos da dispersividade estão
relacionados com duas constantes: a dispersividade longitudinal (
transversal (
T),
L)
e a dispersividade
ambas com dimensões de comprimento [L]. Scheidegger (1961), apud Bear
(1972), considerando as propriedades de simetria do tensor de dispersividade para um meio
isotrópico, propõe a seguinte expressão:
78
L
ijkl
onde
mn
T ij kl
ik
jl
il
jk
(6.39)
é a função delta de Kronecker, definida por
0 para m
1 para m
mn
onde
T
2
n
n
(6.40)
m e n podem representar i, j, k ou l.
Substituindo (6.39) em (6.38), chega-se a
D d,ij
vi v j
T v ij
L
T
v
(6.41)
que determina uma forma geral para a determinação dos componentes do coeficiente de
dispersão mecânica.
Mas a dispersão hidrodinâmica, como descrito no Item 2.4.4, é a união dos
processos de dispersão mecânica e difusão molecular. Assim, o coeficiente de dispersão,
Equação (2.22), tem seus componentes determinados por
D ij
vi v j
T v ij
L
T
v
Dm
(6.42)
Com isso, para o caso bidimensional, a Equação (6.42) determina quatro componentes
para o coeficiente de dispersão como
D xx
D zz
D xz
Lvx
2
T vz
2
v
Tvx
2
L vz
T
L
v
79
(6.43a)
Dm
(6.43b)
2
v
D zx
Dm
vx vz
Dm
(6.43c)
onde
Dxx e Dzz = componentes do tensor do coeficiente de dispersão nas direções dos
eixos x e z, respectivamente, [L2T-1] e
Dxz e Dzx = componentes do tensor do coeficiente de dispersão nas direções
cruzadas xz e zx, respectivamente [L2T-1].
6.6- CRITÉRIO DE ESTABILIDADE
A solução explícita montada pela Equação (6.24) está sujeita a problemas de
instabilidade, devidos, principalmente, ao esquema de integração no tempo. Para garantir a
estabilidade, alguns critérios devem ser adotados. Normalmente, esses critérios dependem de
relações entre entidades geométricas do domínio e o intervalo de tempo. Como o domínio é
fixo, na maioria das vezes, a forma mais fácil de garantir a estabilidade requer que sejam
feitas mudanças no passo de tempo, subdividindo o intervalo t em intervalos menores.
Para este trabalho, a formulação não foi formalmente analisada para se verificar o
melhor critério de estabilidade. O único critério analisado faz referência principalmente ao
movimento de solutos por advecção. Com isso, basicamente, analisa-se o movimento das
partículas no domínio.
O movimento das partículas é simulado por um movimento retilíneo uniforme entre
dois passos de tempo, descrito pela Equação (6.4). Logicamente, em regiões onde as linhas de
fluxo são curvas, haverá erro nessa aproximação, proporcional ao intervalo de tempo
t,
como mostrado na Figura 6.4. Esse erro produz algum erro na solução numérica, sendo,
portanto, razoável procurar limitar a movimentação das partículas.
Posição inicial
Posição real
Posição calculada
Linha de fluxo
Figura 6.4- Erro de aproximação na movimentação das partículas
80
Essa limitação está relacionada com a geometria da malha. Não é conveniente que
uma partícula em movimento passe sobre um elemento sem que esse possa “senti-la” em sua
área. Por isso, o domínio deve ter um comprimento característico relativo aos seus elementos,
o qual seria considerado o menor comprimento representativo do domínio (Lm).
Considerando então a partícula de maior velocidade em movimento retilíneo, essa limitação
traz, da Equação (6.4) em forma vetorial,
dp
onde
v t
Lm
(6.44)
dp = deslocamento da partícula [L] e
Lm = menor comprimento representativo do domínio [L].
Isolando-se t, tem-se o critério de estabilidade adotado
t
Lm
v
(6.45)
Esse critério é conhecido como condição de Courant e é muito utilizado em estudos
de advecção pura e propagação de ondas (Abbott e Basco, 1989; Sun, 1996).
Esse critério é aplicado para a velocidade máxima do domínio. Neste estudo, como a
velocidade das partículas é sempre calculada como uma interpolação linear entre as
velocidades dos nós da malha, a máxima velocidade do domínio sempre ocorre sobre um nó.
Com isso, é mais conveniente e mais prático utilizar a máxima velocidade dos nós do que das
partículas. Dessa forma, o critério adotado torna-se mais restritivo que o critério de Courant,
já que o nó de velocidade máxima nem sempre faz parte do elemento correspondente ao
comprimento Lm.
6.7- EVENTOS PARTICULARES
Durante a movimentação das partículas, podem ocorrer alguns eventos particulares
que conduzam a problemas na resolução das equações do modelo, afetando a sua precisão.
São dois problemas básicos, ambos relacionados com a discretização adotada na formulação
do modelo.
81
Sabe-se que em contornos onde o domínio é fechado para o fluxo, água e soluto não
atravessam os limites. Porém, devido ao movimento retilíneo adotado para as partículas, pode
ocorrer o caso de o contorno ser atravessado e a partícula cair fora do domínio, trazendo o
primeiro problema para o modelo. Então, torna-se necessário alocá-la numa posição próxima
da que a sua trajetória a levaria. Isso é feito por meio de um rebatimento da posição da
partícula em relação ao contorno atravessado, resolvendo assim o problema.
O segundo problema surge quando algum elemento que originalmente gerou partículas
fica sem nenhuma partícula em sua área, o que ocorre em regiões de fluxo divergente, onde as
linhas de fluxo se afastam, e as partículas conduzidas não abrangem toda a região, não
penetrando no elemento ali colocado. Isso pode fazer com que o elemento passe a não receber
contribuição da advecção, gerando um maior erro para a dispersão. Esse problema é resolvido
ao se realizar o reposicionamento das partículas, de modo que o elemento volte a participar
novamente da advecção, reiniciando, assim, o processo. Para isso, foram implementados, no
modelo, dois tipos de controle para esse reposicionamento:
• tipo 1 – sempre que o número de partículas que tiveram mudança no seu nó associado
ultrapassa 10% do total de partículas dentro do domínio, todas as partículas são
retiradas e há um novo reposicionamento; e
• tipo 2 – sempre que pelo menos um elemento que originalmente gerou partículas fica vazio,
há um novo reposicionamento.
O controle do tipo 1 está relacionado, unicamente, com o número de partículas
presentes no domínio, no tempo em que for feita a verificação, enquanto que o tipo 2
relaciona-se com a situação dos elementos no domínio. O percentual escolhido para o controle
do tipo 1 (10%) foi selecionado arbitrariamente. Logicamente, nesse tipo existe um
reposicionamento de partículas mais rápido que para o tipo 2. Os efeitos desses controles
serão testados no modelo.
82
7- ESTRUTURA DO PROGRAMA E TESTES DO MODELO
O programa SSTRANS foi desenvolvido tomando como base o programa SSFLO,
desenvolvido por Koide (1990) e modificado por Campos (1998), o qual é utilizado para o
cálculo do fluxo em solos não-saturados, ou seja, para a determinação dos valores de
e ,
em regime transiente. O programa SSTRANS utiliza então esses valores para a resolução da
equação do transporte de solutos, Equação (2.48), empregando para isso as técnicas numéricas
dos métodos dos elementos finitos, com elementos triangulares e funções de base lineares, e
das características, ambos descritos no Item 6.
O programa é escrito em linguagem de programação estruturada FORTRAN 77.
O programa SSTRANS não exigiu mudanças na estrutura do programa SSFLO, a
qual foi mantida, modificando-se, apenas quando necessário, as funções de solo C( ), K( ) e
( ).
Para a análise do modelo, foram verificados alguns modos de cálculo baseados no
padrão de geração e distribuição de partículas e no tipo de reposicionamento. Foram
implementados quatro padrões de geração e distribuição de partículas e dois tipos de
reposicionamento, descritos nos Itens 6.1.1 e 6.7, respectivamente. Os padrões de geração e
distribuição de partículas são controlados por duas variáveis no programa – NPPD, que
controla a distribuição pelo domínio, e NPPN, que controla o número de partículas por
elemento – que podem assumir os valores 0 e 1. Por isso, foram denominados de opções,
sendo assim chamados nos capítulos posteriores:
• opção 00 – partículas geradas sobre todo o domínio, com 3 por elemento;
• opção 01 – partículas geradas sobre todo o domínio, com 6 por elemento;
• opção 10 – partículas geradas somente nos nós contaminados, com 3 por elemento; e
• opção 11 – partículas geradas somente nos nós contaminados, com 6 por elemento.
7.1- ESTRUTURA DO PROGRAMA SSTRANS
O programa SSTRANS é estruturado em subprogramas que realizam funções gerais
ou específicas. O propósito deste item é fornecer uma visão global da estrutura em que foi
montado o programa, com seus subprogramas. A Figura 7.1 representa essa estrutura, onde é
indicada, entre parênteses, a atuação do programa SSFLO.
83
INÍCIO
ENTRADA DE DADOS (SSFLO)
CÁLCULOS INICIAIS (SSFLO)
MONTAGEM DE MATRIZES E VETORES AUXILIARES
CÁLCULO DO MENOR COMPRIMENTO REPRESENTATIVO DO DOMÍNIO
DETERMINAÇÃO DO FLUXO -
,
(SSFLO)
GERAÇÃO E DISTRIBUIÇÃO DE PARTÍCULAS
MOVIMENTAÇÃO DAS PARTÍCULAS - ADVECÇÃO
CÁLCULO DA DISPERSÃO
FIM DO
TRANSPORTE?
NÃO
SIM
FIM DA
SIMULAÇÃO?
NÃO
SIM
FIM
Figura 7.1- Estrutura do programa SSTRANS
7.1.1- Entrada de Dados
Os dados para a simulação de um problema no programa SSFLO/SSTRANS são
entrados diretamente pela rotina de leitura do programa SSFLO, através de um arquivo
específico. Esses dados podem ser de dois tipos:
a) dados de descrição do problema: caracterizam o problema físico a ser simulado. São
dados da geometria do domínio, posição e número dos nós, elementos e suas
84
conectividades com os nós, fluxos de água e soluto no contorno e nós onde ocorrem,
nós de carga e/ou concentração conhecidas, condições iniciais, características dos
solos e solutos etc.; e
b) dados de descrição da simulação numérica: caracterizam a forma de resolução do
problema e a precisão desejada. São dados sobre os modos da resolução das
equações, a duração dos passos de tempo, as tolerâncias etc.
Outras características importantes de descrição do problema fazem parte do próprio
corpo do programa SSFLO, inseridas nele como subprogramas. São as funções de solo, que,
por serem específicas para cada problema, devem ser modificadas. Os subprogramas que as
representam e calculam são:
• FNCPT – calcula o conteúdo de água em função do potencial mátrico, ( );
• FNCPK – calcula a condutividade hidráulica em função do potencial mátrico, K( ); e
• FNCPC – calcula o coeficiente de capacidade, descrito em Campos (1998), em função do
potencial mátrico, Cc( ).
Além disso, como o programa trabalha por alocação de memória computacional, o
número máximo de nós, elementos, partículas, nós de carga conhecida etc. também devem ser
modificados, se necessário, para garantir a simulação. Essas modificações exigem sempre
uma recompilação do programa.
7.1.2- Montagem de Matrizes e Vetores Auxiliares
Para realizar a movimentação das partículas de uma forma correta, são montadas
algumas matrizes e vetores auxiliares que armazenam informações sobre o domínio e suas
entidades, os nós e elementos. Isso é realizado pelo subprograma INITLC2, que analisa os
dados do domínio e monta as seguintes matrizes e vetores:
• matriz de conectividade elemento-elemento (IELEL): informa, para cada elemento, quais
outros elementos possuem pelo menos um nó comum, determinando a área à sua
volta;
• matriz de conectividade nó-elemento (NOEL): informa, para cada nó, quais elementos a
ele conectados, determinando também a área à sua volta;
• vetor de nós do contorno (NOCONT): informa os pares de nós que fazem parte do
contorno; e
85
• vetores das coordenadas dos centróides dos elementos (XC e ZC): armazena, para cada
elemento, as coordenadas dos seus centróides.
7.1.3- Determinação do Menor Comprimento Representativo do Domínio
De acordo com o Item 6.6, o critério de estabilidade de Courant, Equação (6.45), é
utilizado para determinar o menor passo de tempo necessário à estabilidade da solução
numérica, baseado na relação entre a velocidade e a distância percorrida pelo soluto. Para
garantir que esse critério fosse relacionado à malha de elementos de uma forma geral, essa
distância foi escolhida como sendo o menor comprimento representativo do domínio (L m).
Esse comprimento é calculado pelo subprograma DETHMIN.
O domínio é composto por elementos triangulares. Então, para cada elemento, um
comprimento representativo é o que ocorre entre um vértice e o seu lado oposto, obtido pelo
dobro da relação entre a sua área e o comprimento do lado oposto em questão, sendo adotado,
por segurança, a sua metade. O menor comprimento representativo do domínio é então o
mínimo entre todos aqueles encontrados para os elementos.
7.1.4- Determinação da Velocidade de Escoamento
Após cada passo de tempo para o cálculo do fluxo, um conjunto de valores de
é
fornecido e, com isso, pode-se calcular a velocidade de escoamento em cada nó, para esse
passo de tempo. O cálculo da velocidade é realizado pelo subprograma VELOC, conforme
formulação descrita no Item 6.4. Para a obtenção dos valores da ( ) e K( ), são feitas
chamadas aos subprogramas FNCPT e FNCPK. Os sistemas de equações lineares para vx e vz
são montados e então resolvidos pelo subprograma SONAG2.
Calculadas as velocidades nos nós, o subprograma VELOC faz ainda uma varredura
nos valores das velocidades (v) e seleciona o maior valor. Essa velocidade máxima passa a ser
a velocidade representativa do domínio, a qual é utilizada juntamente com o menor
comprimento representativo do domínio para calcular o valor de t que satisfaça a condição
de Courant, mesmo que ocorram em nós distintos. Isso, a princípio, garante que essa condição
seja satisfeita para todo o domínio.
86
7.1.5- Geração e Distribuição de Partículas
Antes do início da simulação, no arquivo de dados, são pré-definidos os modos de
geração e distribuição de partículas sobre o domínio. Com base nessa escolha, o subprograma
POSIC gera e posiciona as partículas.
No primeiro passo de tempo, as partículas são geradas e posicionadas nos seus locais
originais. A partir daí, a cada vez que o subprograma é chamado, é verificada a necessidade
de se reposicionar as partículas. Se for necessário, as antigas partículas são desativadas e
novas partículas são geradas e distribuídas.
7.1.6- Movimentação das Partículas – Advecção
Para o cálculo do processo de advecção, as partículas devem ser movimentadas de
acordo com o campo de velocidades gerado. A movimentação de cada uma é feita pelo
subprograma MOVEM, seguindo o descrito no Item 6.1.2. No entanto, devido aos problemas
específicos que ocorrem com as partículas, esse é o mais longo e mais complexo subprograma
do programa, pois inclui o tratamento das partículas nos vários contornos, o gerenciamento
eficiente das partículas, etc. Como a velocidade adquirida pelas partículas é a velocidade

corrigida ( v c ), é necessário sempre calcular o fator de retardação (R), o que é feito pelo
subprograma CALCD, que chama os subprogramas FNCPT e FNCDK para o cálculo de
( ) e do coeficiente de partição (Kp), respectivamente.
Cada partícula movimentada tem sua nova posição analisada, verificando se pertence a
algum elemento. Se pertencer, é verificado a qual nó está associada. Se este nó for um nó de
descarga, a partícula é desativada. Também é verificado se houve mudança do nó associado.
Se houve, é verificado se o nó anterior é um nó fonte, pois nesse caso outra partícula deve ser
gerada e posicionada na posição original da anterior. Se a nova posição da partícula não
pertencer a nenhum elemento, então a partícula saiu do domínio e deve ser reposicionada para
dentro, por meio de rebatimento. Para isso, é determinado o contorno que ela atravessou.
Todos esses procedimentos utilizam as matrizes auxiliares, construídas com o subprograma
INITL2.
Para ilustrar melhor o procedimento empregado pelo subprograma MOVEM, a Figura
7.2 apresenta o seu fluxograma.
87
INÍCIO - MOVEM
FIM - MOVEM
CALCULA O INTERVALO DE
TEMPO PARA O TRANSPORTE
SIM
REPASSA A VARIAÇÃO DE
CONCENTRAÇÃO NOS NÓS, DEVIDA À
DISPERSÃO, PARA AS PARTÍCULAS
INICIA A MOVIMENTAÇÃO
NÃO
ÚLTIMO
MOVIMENTO?
SELECIONA A PRÓXIMA PARTÍCULA
A PARTÍCULA
ESTÁ ATIVA?
NÃO
SIM
RECONHECE O NÓ E O
ELEMENTO ASSOCIADOS
CALCULA A VELOCIDADE DA PARTÍCULA
POR INTERPOLAÇÃO ENTRE OS NÓS DO
ELEMENTO ASSOCIADO
MOVE A PARTÍCULA;
CALCULA AS NOVAS COORDENADAS
A PARTÍCULA
ESTÁ DENTRO
DO DOMÍNIO?
CALCULA A
DISPERSÃO
NÃO
REPOSICIONA A PARTÍCULA POR
REBATIMENTO EM RELAÇÃO AO
CONTORNO ATRAVESSADO
DETERMINA A
CONCENTRAÇÃO
DOS NÓS, DEVIDA À
ADVECÇÃO
NÃO
SIM
ÚLTIMA
PARTÍCULA?
DETERMINA O ELEMENTO E
O NÓ ASSOCIADOS À
PARTÍCULA
PARA O NÓ ASSOCIADO, SOMA A
CONCENTRAÇÃO DA PARTÍCULA E
CONTA MAIS UMA;
PARA O ELEMENTO ASSOCIADO,
CONTA MAIS UMA PARTÍCULA
GERA UMA NOVA PARTÍCULA OU
REATIVA ALGUMA DESATIVADA NA
POSIÇÃO ORIGINAL DA ANTERIOR
SIM
O NÓ ASSOCIADO
DA PARTÍCULA
MUDOU?
NÃO
A PARTÍCULA FOI
CRIADA NESSE
NÓ?
NÃO
SIM
SIM
O NOVO NÓ
ASSOCIADO É UM
NÓ DE DESCARGA?
NÃO
O NÓ ASSOCIADO
ANTERIOR É UMA
FONTE DE SOLUTO?
NÃO
SIM
DESATIVA A
PARTÍCULA
Figura 7.2- Fluxograma do subprograma MOVEM
88
SIM
É nesse subprograma que é aplicada a condição de Courant. O valor de t é calculado.
Esse valor é então comparado com a metade do intervalo de tempo utilizado no cálculo do
fluxo, não podendo ser maior. Com isso, o menor intervalo é utilizado. Isso é feito para
garantir maior segurança quanto à estabilidade, forçando ainda que a cada iteração no tempo
do fluxo sejam feitas, no mínimo, duas iterações para o transporte.
7.1.7- Cálculo da Dispersão
Após a advecção, com as concentrações Ck* determinadas, é chamado o subprograma
DISPERSION, que calcula o processo da dispersão. O objetivo desse subprograma é,
basicamente, calcular as variações de C devidas à dispersão. Para isso, é necessária a
montagem das matrizes globais [M]k, [M]k*, [G]k e [G]k* e dos vetores globais {N}k e {N}k*,
o que é realizado por uma chamada ao subprograma ASSEMBCON. Esse utiliza o
subprograma CALCD para o cálculo dos termos do tensor de dispersão, aplicando-os ao
subprograma MATLOC para determinar a contribuição de cada elemento às matrizes globais,
montadas finalmente pelo subprograma MATGLB.
Assim, as variações de C devidas à dispersão são calculadas explicitamente e
repassadas para os nós.
A Figura 7.3 esquematiza a organização dos subprogramas que compõem o programa
SSTRANS.
7.2- TESTES DE VERIFICAÇÃO DO MODELO
Para verificar a capacidade do modelo e a sua precisão e eficiência, foram realizados
cinco testes utilizando problemas apresentados na literatura, cujas soluções são conhecidas.
Cada teste tem um propósito específico, visando analisar os principais aspectos do modelo
desenvolvido.
89
SSTRANS
INITLC2
DETHMIN
TRANSP
VELOC
FNCPT
FNCPK
SONAG2
POSIC
MOVEM
CALCD
FNCPT
FNCDK
DISPERSION
ASSEMBCON
CALCD
FNCPT
FNCDK
MATLOC
MATGLB
Figura 7.3- Diagrama esquemático/hierárquico do programa SSTRANS
7.2.1- Teste 1 – Injeção Contínua de Fluido com Soluto no Topo de uma Coluna
Vertical de Solo
Este primeiro teste foi realizado com o propósito de verificar o comportamento do
modelo para uma situação de fluxo permanente e uniforme em uma coluna vertical de solo
com transporte de soluto transiente, obedecendo uma condição de contorno do tipo de
Cauchy, isto é, com taxa de entrada de massa de soluto conhecida e dependente da
concentração. Neste teste, os resultados do modelo são comparados com resultados analíticos.
90
O teste envolve uma coluna vertical de solo para a qual existe, no seu topo, uma
injeção contínua de água com uma determinada concentração de soluto. O solo é homogêneo
e isotrópico e o fluxo de água em seu interior é mantido constante e unidimensional por meio
da vazão de água constante no topo e o estabelecimento de uma carga constante na base da
coluna. A vazão é tal que o solo permanece não-saturado em toda a sua extensão, por toda a
simulação, garantindo uma distribuição permanente e uniforme do conteúdo de água e,
conseqüentemente, do potencial mátrico. No início da simulação, a coluna está livre de
qualquer concentração de soluto e a única fonte de soluto é a água injetada no topo, a qual é
mantida com uma concentração constante durante toda a simulação. Não são considerados os
processos de produção e decaimento e de sorção de solutos.
A solução analítica utilizada para esse teste foi apresentada por van Genuchten (1981)
e trata do transporte de solutos unidimensional com fluxo permanente e uniforme em uma
coluna vertical semi-infinita de solo. A equação diferencial que rege o problema é dada por
R
onde 0
z
C
t
2
DL
C
z2
vz
C
z
(7.1)
.
Originalmente desenvolvida para considerar os termos de produção de soluto de
ordem zero e a sorção, a Equação (7.1) foi adaptada ao teste ao se fazer R = 1 e
= 0.
As condições iniciais e de contorno do sistema são
C z,0
0
DL
C
z
C
z z
para z
vzC
(7.2a)
0
v z C0
para
t
0
(7.2b)
z 0
finito
para
t
0
(7.2c)
A condição de contorno (7.2b) estabelece que qualquer entrada de massa de soluto no sistema
pode ser considerada como uma advecção pelo contorno, estando de acordo com a
consideração feita no Item 6.2 de que o fluxo dispersivo pelo contorno pode ser desprezado. A
condição de contorno (7.2c) estabelece que quando z tende para o infinito, o gradiente de
91
concentração permanece com um valor finito em qualquer tempo, o que, segundo van
Genuchten (1981), descreve mais corretamente o sistema semi-infinito.
Com isso, a solução do problema descrito é dada por:
C z, t
C0
1
1
2
vz2t
DL
z vz t
1
erfc
2
2 DLt
vzz
DL
1
2
exp
z
vz t
4D L t
2
(7.3)
2
vz t
v z
z vzt
exp z erfc
DL
DL
2 DLt
onde o valor de z representa a profundidade na coluna e erfc é a função erro complementar.
Como o modelo é bidimensional, para a realização do teste, que é unidimensional, fo i
adotada uma malha regular e simétrica, com a configuração apresentada na Figura 7.4,
suficientemente longa, de modo a suprimir ao máximo possível os efeitos transversais.
h
Contorno com fluxo de água com
soluto conhecido
h
Contorno de carga constante
b b
Figura 7.4- Malha de elementos finitos para o teste 1 (caso unidimensional)
92
O teste foi realizado para uma coluna de solo com altura total de 200cm em que todos
os elementos são idênticos em forma, possuindo altura h = 2cm e largura b = 1cm, sendo no
total 303 nós com 400 elementos.
Para procurar assemelhar o teste com uma situação real, as características hidráulicas
do solo adotadas foram as apresentadas por Vauclin et al. (1979), que realizou experimentos
com areia fina de granulometria uniforme. Tais características variam com o potencial
mátrico, sendo dadas por
40000
sat
onde
sat
2, 9
(7.4)
0
= 0,30 e
K
onde
40000
K sat
2990000
2990000
5, 0
0
(7.5)
Ksat = 35cm/h.
7.2.2- Teste 2 – Concentração Constante no Topo de uma Coluna Vertical de Solo
Este teste foi realizado com o objetivo de verificar o comportamento do modelo para
uma situação de fluxo permanente e uniforme em uma coluna vertical de solo com transporte
de soluto transiente, obedecendo uma condição de contorno do tipo de Dirichlet, isto é,
concentração conhecida no contorno. Novamente, os resultados do modelo são comparados
com resultados analíticos.
A situação é semelhante à do teste 1, exceto que a única fonte de massa de soluto para
o sistema é a presença de uma faixa de solo com concentração constante no topo, ao invés do
fluxo conhecido.
Este teste é um dos cinco problemas de referência adotados pelo Fórum de AdvecçãoDifusão em conjunção com a 7ª Conferência Internacional em Recursos Hídricos, realizado
em 1988 no Massachussets Institute of Technology (Yeh, 1990), somente com algumas
modificações nas condições iniciais.
93
O problema é descrito pela equação
C
t
com 0
z
2
DL
C
z2
vz
C
z
(7.6)
, que rege o transporte unidimensional de solutos com fluxo permanente e
uniforme e sem consideração dos processos de produção/decaimento e de sorção de soluto,
em uma coluna vertical semi-infinita.
As condições iniciais e de contorno do sistema são
C z,0
0
para z
C 0, t
C0
para
C ,t
0
para
t
t
(7.7a)
0
0
0
(7.7b)
(7.7c)
A solução analítica para esse problema foi obtida por Ogata e Banks (1961), apud
Fetter (1993), sendo dada por
C(z, t )
C0
z vzt
erfc
2
2 DL t
exp
vzz
z vzt
erfc
DL
2 DL t
(7.8)
A obtenção dessa solução é desenvolvida por Sun (1996), considerando o decaimento de
soluto.
Como a situação é semelhante à do teste 1, a mesma malha de elementos finitos foi
adotada, tendo sido modificada somente nos nós do topo da coluna, que passaram ser nós de
concentração constante. Também foram utilizadas as mesmas características hidráulicas para
o solo.
7.2.3- Teste 3 – Transporte de Soluto em uma Coluna Vertical de Solo com Fluxo
Transiente
Este terceiro teste foi realizado com o intuito de verificar o comportamento do modelo
para uma situação onde o fluxo de água é transiente e em forma de frente de umidade. Os
94
resultados do teste são comparados com resultados numéricos (van Genuchten, 1982; Voss,
1984; Healy, 1990), obtidos com base em um experimento de campo (Warrick et al., 1971).
Uma coluna vertical de solo, inicialmente não-saturada, passa a receber água pelo seu
topo, o qual é mantido saturado, com carga constante durante toda a simulação. A água se
infiltra no solo produzindo uma frente de umidade que se desloca, com o tempo, para a base
da coluna, que é mantida não-saturada por toda a simulação, devido ao estabelecimento de
uma carga constante, inferior à carga relativa à saturação. A simulação é conduzida em duas
partes. No início da primeira, a distribuição da carga pela coluna é conhecida e a coluna está
livre de qualquer concentração de soluto. A fonte de soluto do sistema é o topo, mantido
também com uma concentração constante durante essa parte, produzindo assim também uma
frente de contaminação. Na segunda parte, a água que flui pelo topo da coluna tem a sua
concentração anulada, de modo que a água sem soluto passa agora a alimentar o sistema. O
soluto é conservativo e, por isso, não são considerados os processos de produção e
decaimento e de sorção.
Devido à semelhança, a mesma malha de elementos utilizada nos testes 1 e 2 é
utilizada aqui, sendo que a única modificação são os nós no topo da coluna que agora
possuem carga constante.
As propriedades hidráulicas do solo variam com a carga e são dadas por
K
K sat 51164
3, 4095
K sat 13,672
0,97814
para
0,6829 0,09524 ln
0,4531 0,02732 ln
29,484cm
(7.9)
para
29,484
14,495cm
para
para
29,484cm
29,484
14,495cm
Os outros parâmetros hidráulicos do solo são Ksat = 1,574cm/h e
L
(7.10)
= 1,026cm.
As condições iniciais do sistema são
z,0
0
269,169 exp
z
z
114,288
159,23cm
e as condições de contorno são
95
para 0
z
para
60cm
z
60cm
(7.11)
0, t
200, t
C 0, t
(7.12a)
14,495cm
(7.12b)
159,23cm
C0
209meq / l
0
para t
para t
2,8h
2,8h
(7.12c)
A simulação é realizada para um período total de 9 horas e, como mostra a Equação
(7.12c), a primeira parte compreende as 2,8 horas iniciais, quando há entrada de soluto no
sistema, e a segunda compreende as 6,2 horas finais.
7.2.4- Teste 4 – Concentração Constante no Topo de um Perfil Vertical de Solo
O objetivo deste teste foi verificar o comportamento do modelo quando o domínio de
fluxo é bidimensional, onde tem efeito o processo de dispersão transversal, com os resultados
sendo comparados com uma solução analítica.
Um perfil vertical de solo recebe, no seu topo, uma injeção contínua de água. Essa
injeção é tal que mantém o fluxo de água constante e unidirecional, mantida pelo
estabelecimento de uma carga constante no topo e no pé do perfil, permanecendo o solo nãosaturado por toda a simulação. Com isso, o conteúdo de água e o potencial mátrico nos poros
são mantidos com uma distribuição uniforme e permanente no perfil. A concentração de
soluto em todo o perfil é nula no início e a fonte de soluto se encontra somente na água
injetada na metade esquerda do topo, mantida com uma concentração constante por toda a
simulação. Os processos de produção/decaimento e de sorção de soluto são desprezados.
A equação diferencial que descreve o problema é dada por
C
t
2
DL
2
C
z2
DT
O domínio de fluxo é um semi-plano com z
C
x2
vz
(7.13)
0, sendo o eixo orientado para baixo, e os
outros contornos no infinito.
As condições iniciais e de contorno do sistema são
96
C
z
C x, z,0
C
z z
0
para z
0
C
x
e
C0
C0
2
0
C x,0, t
0 e
(7.14a)
x
(7.14b)
0
x
para
x
0
para
x
0
para
x
0
(7.14c)
A Equação (7.14c) determina que só há entrada de soluto pela metade esquerda do
topo do perfil.
A solução analítica para esse problema foi encontrada por Leij e Dane (1990) e é dada
por
C x, z, t
z C0
4
DL
t
0
3
2
2
erfc
x
2 DT
exp
z vz
2 DL
d
(7.15)
A integral em (7.15) não tem solução analítica e precisa ser calculada numericamente
(Leij e Dane, 1990). Para este estudo, foi utilizada a integração numérica pelo método do
trapézio, com pequenos intervalos.
Para esse teste foi adotada a malha mostrada na Figura 7.5 que apresenta simetria em
relação ao eixo central e entre cada duas colunas. Nas colunas centrais foi adotada uma
largura menor que nas outras, buscando com isso aumentar a densidade de nós naquela região,
por ser uma região de interesse, pois é o limite físico entre a região que recebe a massa de
soluto e a que não recebe, onde terão influência significativa os efeitos da dispersão
transversal.
O perfil de solo foi simulado com altura e largura totais de 20cm cada. Assim, cada
elemento teve altura h = 2cm em toda a malha e a largura foi de b = 1cm para os elementos
centrais e, conseqüentemente, de 2cm para os outros, num total de 143 nós e 240 elementos.
As características hidráulicas do solo são as mesmas utilizadas para os testes 1 e 2,
porém adotando-se
sat
= 0,5 e Ksat = 72,71cm/h.
97
h
h
h
h
h
h
h
h
h
h
2b
2b
2b
2b
b b b b
2b
2b
2b
2b
Contorno de carga e concentração (nula) constantes
Contorno de carga e concentração (não nula) constantes – fonte de soluto
Contorno de carga constante
Figura 7.5- Malha de elementos finitos para o teste 4 (caso bidimensional)
7.2.5- Teste 5 – Transporte com Produção/Decaimento e Sorção de Soluto em uma
Coluna Vertical de Solo
Este último teste foi realizado para verificar o comportamento do modelo quando são
considerados os processos de produção e decaimento e de sorção de soluto. Novamente, os
resultados obtidos com o modelo são comparados com uma solução analítica.
A situação é semelhante à do teste 1, acrescentando somente a consideração dos
processos de produção/decaimento e de sorção do soluto. No caso, foi utilizado o modelo de
sorção linear.
Nessas condições, a equação que rege o problema é
98
C
R
t
2
C
z2
DL
vz
C
z
'C
'
(7.16)
sendo
s Bd K d
'
l
'
l
(7.17a)
s Bd
(7.17b)
que trata do transporte de solutos em regime de fluxo unidimensional, uniforme e permanente.
As condições de contorno e iniciais do sistema são
C z,0
0
DL
C
z
C
z z
0
(7.18a)
v zC0
(7.18b)
para z
vzC
z 0
0
para
t
0
(7.18c)
A solução analítica desse problema foi apresentada por van Genuchten (1981) e é dada
por
C z, t
'
'
C0
F1 z, t
F2 z, t
(7.19)
onde
vz
F1 z, t
vz
vz
vz
u
2
u
exp
exp
vz u z
2D L
vz u z
2D L
vz
v z
exp z
2 ' DL
DL
ut
R
erfc
erfc
erfc
Rz ut
2 D L Rt
Rz ut
2 D L Rt
Rz
vz t
2 D L Rt
99
(7.20a)
F2 z, t
'
exp
'
vzz
1
1
2
DL
u
vz 1
4 ' DL
't
R
Rz v z t
1
erfc
2
2 D L Rt
vz2t
v z
exp z
DLR
DL
erfc
Rz
2
vz t
RD L
1
2
exp
Rz v z t
4D L Rt
2
(7.20b)
vz t
2 D L Rt
'
'
'
exp
'
't
R
1
2
(7.20c)
vz2
Devido à semelhança, a mesma malha de elementos utilizada nos testes 1 e 2 é
utilizada aqui, sendo que a única modificação são os nós no topo da coluna que agora
possuem carga constante.
As características hidráulicas do solo são as mesmas adotadas para o teste 4.
100
8- RESULTADOS
Para avaliar o desempenho do modelo quanto à sua precisão e eficiência, foram
realizados cinco testes baseados em problemas com soluções conhecidas, descritos no Item
7.2. Essa avaliação foi feita principalmente em relação aos modos de cálculo, baseados nos
padrões de geração e distribuição de partículas e nos tipos de controle utilizados para o
reposicionamento das partículas.
Os quatro tipos de padrão de geração e distribuição de partículas são:
• opção 00 – partículas geradas sobre todo o domínio, com 3 por elemento;
• opção 01 – partículas geradas sobre todo o domínio, com 6 por elemento;
• opção 10 – partículas geradas somente nos nós contaminados, com 3 por elemento; e
• opção 11 – partículas geradas somente nos nós contaminados, com 6 por elemento.
Os dois tipos de controle de reposicionamento são:
• tipo 1 – quando o número de partículas que trocam o nó associado ultrapassa 10% do total; e
• tipo 2 – quando pelo menos um elemento que originalmente gerou partículas fica vazio.
Também foi verificado o comportamento do modelo quando se altera o número de
Peclet, modificando-se o valor do coeficiente de dispersão. Com isso, altera-se o tipo de
regime de transporte entre predominante advectivo e predominante dispersivo.
Para cada teste realizado, os resultados são analisados separadamente. Ao final, é feita
uma análise geral dos resultados.
8.1- RESULTADOS DO TESTE 1
A configuração geral do teste é descrita no Item 7.2.1.
Foram utilizados os seguintes valores na simulação:
altura da coluna: H = 2m;
taxa de entrada de água: qz = 15cm/h; e
concentração de soluto na entrada: C0 = 200mg/l.
Para a simulação, foram utilizados como parâmetros numéricos:
tempo total de simulação: t = 5h;
incremento de tempo padrão: t = 0,05h; e
opção de cálculo para o fluxo: RA Distribuída (Campos, 1998).
101
O tempo total de simulação adotado foi suficiente para que toda a coluna recebesse
massa de soluto da fonte, adquirindo a concentração máxima de 200mg/l. Como o regime de
fluxo é permanente, a opção adotada para este cálculo não teve influência nos valores obtidos.
Os principais resultados para o fluxo obtidos com os valores acima foram:
conteúdo de água:
= 0,257, em toda a coluna; e
velocidade de escoamento: vz = 58,42cm/h, em toda a coluna.
Esses resultados são empregados na solução analítica, Equação (7.3).
Para verificar a influência dos padrões de geração e distribuição de partículas, o
problema foi simulado para as quatro opções. Foi utilizada uma dispersividade
L
= 0,3cm,
com o controle de reposicionamento do tipo 1. Os resultados são mostrados na Figura 8.1.
Opção 00
Opção 01
Opção 10
Opção 11
200
180
Concentração (mg/l)
160
140
1 hora
120
2 horas
3 horas
100
80
60
40
20
0
0
20
40
60
80
100
120
Profundidade (cm)
140
160
180
200
Figura 8.1- Influência do padrão de geração e distribuição das partículas para o teste 1
Como pode-se observar, existe uma grande semelhança entre as curvas obtidas,
mostrando que a mudança de opção não teve influência significativa para este teste. O mesmo
teste foi aplicado para o controle de reposicionamento do tipo 2 e, novamente, não houve
diferença significativa entre as curvas obtidas. Devido a isso, os próximos resultados são
apresentados apenas para opção 00.
Para avaliar a precisão do modelo, os resultados obtidos para o fluxo no teste foram
aplicados à solução analítica do problema, Equação (7.3). Assim, os resultados analíticos
foram comparados com os do modelo. Essa comparação foi feita para dois valores de
102
L,
sendo possível, assim, avaliar o efeito do número de Peclet e, conseqüentemente, do tipo de
regime de transporte sobre o comportamento do modelo, observando a ocorrência de
oscilação e dispersão numéricas. Os resultados foram obtidos com controle de
reposicionamento do tipo 1.
Aplicando um valor de
L
= 0,1cm, o valor do coeficiente de dispersão calculado pelo
modelo foi de DL = 5,84cm2/h, fornecendo um número de Peclet Pe = 20. Com esse valor, o
regime de transporte é predominantemente advectivo e os resultados obtidos são mostrados na
Figura 8.2. O erro relativo máximo obtido foi de 31,4%, calculado em relação à concentração
analítica máxima.
Numérico
Analítico
200
180
Concentração (mg/l)
160
140
1,5 hora
120
100
80
60
40
20
0
0
20
40
60
80
100
120
Profundidade (cm)
140
160
180
200
Figura 8.2- Resultados para o teste 1, com alto valor para o número de Peclet
Com o valor de
L
= 2cm, o coeficiente de dispersão calculado pelo modelo foi
DL = 116,84cm2/h, sendo o número de Peclet Pe = 1, caracterizando um regime de transporte
predominantemente dispersivo. Para essa situação, os resultados são mostrados na Figura 8.3
e o erro relativo máximo obtido foi de 9,0%.
Como se vê, em ambos os casos, inclusive nos resultados da Figura 8.1, não houve o
aparecimento de nenhuma oscilação numérica próximo à frente de contaminação e os
resultados numéricos se aproximam suficientemente bem dos resultados analíticos. Isso
mostra que, a princípio, o modelo não sofre influência sensível do tipo de regime de
transporte. No entanto, percebe-se a ocorrência de um pequeno avanço da solução numérica
103
em relação à analítica e uma pequena dispersão numérica, mais perceptível nos pontos de
baixa concentração, o que não inviabiliza os resultados.
Numérico
Analítico
200
180
Concentração (mg/l)
160
140
1,5 hora
120
100
80
60
40
20
0
0
20
40
60
80
100
120
Profundidade (cm)
140
160
180
200
Figura 8.3- Resultados para o teste 1, com baixo valor para o número de Peclet
Foi avaliado também o efeito do tipo de controle de reposicionamento das partículas
sobre a precisão do modelo. Para isso, o problema foi simulado para os dois tipos de controle
e os resultados do modelo foram comparados aos da solução analítica. Foi utilizado um valor
de
L
= 0,3cm, o que fornece DL = 17,53cm2/h. Para os dois tipos de controle, foram obtidos
os resultados mostrados nas Figuras 8.4 e 8.5. Nos tempos de 1, 2 e 3 horas, os erros relativos
máximos, obtidos para o controle do tipo 1, foram, respectivamente, de 18,8%, 20,4% e
22,2% e, para o controle do tipo 2, foram, respectivamente, de 8,1%, 6,2% e 4,6%.
Os resultados mostram uma diferença sensível entre os dois tipos controle de
reposicionamento. Percebe-se que o controle do tipo 1 provoca um pequeno adiantamento da
solução numérica em relação à analítica e uma também pequena, mas sensível, dispersão
numérica. Para o tipo 2, no entanto, ocorre uma concordância quase perfeita da solução
numérica com a analítica, não ocorrendo efetivamente dispersão numérica. Isso mostra que o
tipo de controle de reposicionamento das partículas tem efeito significativo na precisão do
modelo. A causa provável disso é que, quando se usa o controle do tipo1, o freqüente
reposicionamento das partículas pode fazer com que, a cada reposicionamento, a massa de
soluto seja passada diretamente de uma partícula a jusante de um nó para outra posicionada a
104
montante do mesmo. Isso não acontece no controle do tipo 2, pois o reposicionamento das
partículas não é tão freqüente.
Numérico
Analítico
200
180
140
1 hora
120
2 horas
3 horas
100
80
60
40
20
0
0
20
40
60
80
100
120
Profundidade (cm)
140
160
180
200
Figura 8.4- Resultados para o teste 1, com controle de reposicionamento do tipo 1
Numérico
Analítico
200
180
160
Concentração (mg/l)
Concentração (mg/l)
160
140
1 hora
120
2 horas
3 horas
100
80
60
40
20
0
0
20
40
60
80
100
120
Profundidade (cm)
140
160
180
200
Figura 8.5- Resultados para o teste 1, com controle de reposicionamento do tipo 2
105
8.2- RESULTADOS DO TESTE 2
A configuração geral do teste é descrita no Item 7.2.2.
Os valores utilizados na simulação foram:
altura da coluna: H = 2m;
taxa de entrada de água: qz = 15cm/h; e
concentração constante na fonte: C0 = 200mg/l.
Os parâmetros numéricos utilizados foram os mesmos do teste 1.
Os principais resultados obtidos para o fluxo foram:
conteúdo de água:
= 0,257, em toda a coluna; e
velocidade de escoamento: vz = 58,42cm/h, em toda a coluna.
O problema foi simulado para as quatro opções de geração e distribuição, com
L
= 0,3cm e controle de reposicionamento do tipo 1. Os resultados obtidos estão mostrados
na Figura 8.6.
Opção 00
Opção 01
Opção 10
Opção 11
200
180
Concentração (mg/l)
160
140
1 hora
120
2 horas
3 horas
100
80
60
40
20
0
0
20
40
60
80
100
120
Profundidade (cm)
140
160
180
200
Figura 8.6- Influência do padrão de geração e distribuição de partículas para o teste 2
Assim como no teste 1, observa-se que não há influência significativa sobre o
comportamento do modelo, proporcionada pela opção adotada. A semelhança nas curvas
denota isso. Assim, considerando também os resultados do teste 1, pode-se dizer que, a
106
princípio, não há grandes vantagens em utilizar um número elevado de partículas e a
utilização da opção 10 gera uma economia de tempo, em termos de execução do programa, e
de utilização de memória computacional, já que o número de partículas pode ser
extremamente reduzido, sem perda de precisão dos resultados.
A solução do modelo foi, novamente, comparada com a solução analítica, para dois
valores do número de Peclet, para verificar a sua precisão em regimes de transporte
diferentes. Para um valor de
L
= 0,1cm, o coeficiente de dispersão calculado foi
DL = 5,84cm2/h e o número de Peclet Pe = 20. Para esse regime de transporte
predominantemente advectivo, os resultados estão mostrados na Figura 8.7. O erro relativo
máximo nesse caso foi de 31,6%.
Numérico
Analítico
200
180
Concentração (mg/l)
160
140
1,5 hora
120
100
80
60
40
20
0
0
20
40
60
80
100
120
Profundidade (cm)
140
160
180
200
Figura 8.7- Resultados para o teste 2, com alto valor para o número de Peclet
Para
L
= 2cm, os resultados são mostrados na Figura 8.8. Nessa situação, o
coeficiente de dispersão foi DL = 116,84cm2/h, com número de Peclet Pe = 1, em regime de
transporte predominantemente dispersivo. O erro relativo máximo obtido foi de 6,4%.
Pelas Figuras 8.7 e 8.8, percebe-se novamente que não ocorreu o problema da
oscilação numérica próximo à frente da contaminação. O que se vê é apenas um pequeno
adiantamento da frente calculada pelo modelo e uma pequena dispersão numérica em pontos
de baixa concentração, resultado provavelmente do tipo de reposicionamento adotado, o tipo
1. No entanto, os resultados numéricos se aproximam bastante dos analíticos, mostrando uma
107
boa precisão do modelo para o teste, o que evidencia novamente que o modelo não sofre
muito a influência do regime de transporte.
Numérico
Analítico
200
180
Concentração (mg/l)
160
140
1,5 hora
120
100
80
60
40
20
0
0
20
40
60
80
100
120
Profundidade (cm)
140
160
180
200
Figura 8.8- Resultados para o teste 2, com baixo valor para o número de Peclet
Foi avaliado o efeito do tipo de controle de reposicionamento das partículas sobre a
precisão do modelo, simulando o problema para os dois tipos de controle, com
L
= 0,3cm
(DL = 17,53cm2/h). Para os dois tipos de controle, foram obtidos os resultados das Figuras 8.9
e 8.10. Nos tempos de 1, 2 e 3 horas, os erros relativos máximos obtidos para o controle do
tipo 1, foram, respectivamente, de 17,8%, 19,7% e 21,6% e, para o controle do tipo 2, foram,
respectivamente, de 6,9%, 5,3% e 3,9%.
Novamente, existe uma diferença sensível entre os tipos de controle adotados,
mostrando que o reposicionamento das partículas é bastante influente sobre a precisão do
modelo. No controle do tipo 1, existe um avanço da frente de contaminação calculada pelo
modelo e uma dispersão numérica, ambos pequenos, mas sensíveis. No controle do tipo 2,
essa dispersão numérica diminui bastante e os resultados numéricos praticamente se igualam
aos analíticos. Esses resultados permitem afirmar, em primeira análise, que o modelo é
bastante preciso.
108
Numérico
Analítico
200
180
Concentração (mg/l)
160
140
1 hora
120
2 horas
3 horas
100
80
60
40
20
0
0
20
40
60
80
100
120
Profundidade (cm)
140
160
180
200
Figura 8.9- Resultados para o teste 2, com controle de reposicionamento do tipo 1
Numérico
Analítico
200
180
Concentração (mg/l)
160
140
1 hora
120
2 horas
3 horas
100
80
60
40
20
0
0
20
40
60
80
100
120
Profundidade (cm)
140
160
180
200
Figura 8.10- Resultados para o teste 2, com controle de reposicionamento do tipo 2
8.3- RESULTADOS DO TESTE 3
A configuração geral do teste é descrita no Item 7.2.3.
109
Os valores utilizados na simulação foram apresentados no Item 7.2.3.
Os parâmetros numéricos utilizados foram:
tempo total de simulação: t = 9h, sendo 2,8h para a primeira parte e 6,2h para a segunda; e
incremento de tempo padrão: t = 0,025h.
Segundo Campos (1998), a melhor alternativa para o cálculo do fluxo em solos
bastante secos é a opção Neuman Aglutinada (QG). No entanto, para encontrar a melhor
alternativa, procurando minimizar os efeitos do cálculo do fluxo sobre os do transporte, foram
analisadas todas as opções de cálculo, descritas em Campos (1998). As melhores em termos
de precisão em relação à frente de umidade real foram Solução Mista Aglutinada e Solução
Mista Distribuída, apresentando praticamente os mesmos resultados. Optou-se então por
utilizar para este teste a Solução Mista Aglutinada, cujos resultados para o fluxo são
mostrados na Figura 8.11.
Numérico
Real
Conteúdo de água
0,4
0,38
0,36
0,34
0,32
0,3
0,28
0,26
0,24
9 horas
2 horas
0,22
0,2
0,18
0,16
0,14
0
20
40
60
Profundidade (cm)
80
100
120
Figura 8.11- Resultados para o fluxo, no teste3
Nessa figura pode-se verificar a boa proximidade entre as curvas numérica e real, nos
tempos de 2 e 9 horas, sendo que os erros relativos máximos foram de 3,4%, para 2 horas, e
1,3%, para 9 horas, em relação ao conteúdo de água máximo. Isso indica que, aparentemente,
o cálculo do fluxo foi realizado com boa precisão durante toda a simulação. No entanto, o
comportamento do fluxo em tempos intermediários não pôde ser aferido, já que não se dispõe
110
dos dados reais, de forma que qualquer problema surgido entre esses tempos não pôde ser
verificado.
Para testar o transporte, foram, então, verificadas as quatro opções de cálculo, para os
dois tipos de controle de reposicionamento de partículas. Para o controle do tipo 1, foram
obtidos os resultados mostrados na Figura 8.12. No tempo de 2 horas, os erros relativos
máximos foram de 9,0%, para a opção 00, de 4,3%, para a opção 01, de 4,8%, para a opção
10, e de 20,6%, para a opção 11. No tempo de 9 horas, os erros relativos máximos foram de
54,9%, para a opção 00, de 41,8%, para a opção 01, de 53,7%, para a opção 10, e de 27,4%,
para a opção 11.
Opção 00
Opção 01
Opção 10
Opção 11
Real
220
2 horas
200
180
9 horas
Concentração (mg/l)
160
140
120
100
80
60
40
20
0
0
20
40
60
Profundidade (cm)
80
100
120
Figura 8.12- Resultados para o teste 3, mostrando a influência dos padrões de geração e
distribuição de partículas, com controle de reposicionamento do tipo 1
Observando os resultados, percebe-se que o padrão de geração e distribuição de
partículas teve forte influência sobre os resultados. Para todas as opções, em 2 horas, os
resultados são muito semelhantes, de modo que se sobrepõem e se confundem e ocorre um
pequeno adiantamento da solução numérica e uma também pequena dispersão numérica.
Para o tempo de 9 horas, no entanto, as opções diferenciam-se bastante. A opção 00
provoca um adiantamento da solução com o pico de concentração bastante semelhante ao real.
111
Para essa opção, foi verificada a menor dispersão numérica, com o formato da solução
numérica em boa concordância com a real. Para a opção 01, verifica-se que, além do avanço
ocorrido, o pico de concentração é menor que o real, com uma dispersão maior que a anterior.
A opção 10 provoca um atraso considerável e uma maior dispersão numérica, inclusive
mantendo ainda uma pequena concentração próxima ao topo da coluna. O seu pico de
concentração, no entanto, praticamente iguala-se ao real. A opção 11 é a que apresentou os
melhores resultados gerais, com a solução numérica alcançando a real, em termos de
profundidade, e com o pico de concentração quase com mesmo valor. A dispersão, no
entanto, é considerável, fazendo com que, também, ainda permaneça uma pequena
concentração próxima ao topo da coluna.
O motivo dessa grande diferença nos resultados pode ser o regime transiente de fluxo.
Com as condições de umidade sendo modificadas constantemente, o cálculo das velocidades
de fluxo pode não ter uma precisão muito boa. Como esse procedimento cálculo não foi
verificado, essa é uma hipótese a se considerar. Outro motivo para essa diferença pode estar
no próprio cálculo do fluxo, que não pôde ser comparado em tempos intermediários. No
entanto, como os cálculos realizados são os mesmos, era de se esperar que os mesmos
problemas ocorressem em todas as opções. Porém, a forma de distribuição e, naturalmente, a
movimentação das partículas sofrem forte influência do fluxo, tanto que as opções que usam
um número menor de partículas, mantiveram uma pequena concentração próxima do topo e as
duas com número maior de partículas levaram a um avanço da solução numérica.
Para o controle do tipo 2, os resultados obtidos são mostrados na Figura 8.13. Nesse
caso, no tempo de 2 horas, os erros relativos máximos foram de 6,4%, para a opção 00, de
6,2%, para a opção 01, de 23,7%, para a opção 10, e de 23,3%, para a opção 11. No tempo de
9 horas, os erros relativos máximos foram de 16,2%, para a opção 00, de 15,2%, para a opção
01, de 96,4%, para a opção 10, e de 94,7%, para a opção 11.
Novamente, os padrões de geração e distribuição de partículas tiveram efeitos
significativos nos resultados calculados pelo modelo. Na Figura 8.13, pode-se distinguir que
as opções com distribuição sobre todo o domínio foram praticamente idênticas, o mesmo
ocorrendo para as opções com distribuição somente nos nós contaminados. Para o tempo de 2
horas, em todas as opções, a solução numérica teve dispersão inferior à real, de modo que a
frente de contaminação numérica foi mais acentuada.
Para o tempo de 9 horas, verifica-se que as opções 00 e 01 geraram uma solução
numérica que praticamente alcançou a real, em termos de profundidade, com um pico mais
acentuado e um espalhamento menor, mostrando que, apesar do alcance, o processo numérico
112
atenua em parte o processo dispersivo. Para as opções 10 e 11, ocorre o mesmo problema do
controle do tipo 1, porém de forma mais acentuada, que é a concentração mantida nas
proximidades do topo.
Opção 00
Opção 01
Opção 10
Opção 11
Real
220
2 horas
200
180
9 horas
Concentração (mg/l)
160
140
120
100
80
60
40
20
0
0
20
40
60
Profundidade (cm)
80
100
120
Figura 8.13- Resultados para o teste 3, mostrando a influência dos padrões de geração e
distribuição de partículas, com controle de reposicionamento do tipo 2
Este teste mostrou que, para situações de fluxo transiente, o modelo não tem boa
precisão, de modo que apresenta resultados razoáveis qualitativamente mas não muito bons
quantitativamente, necessitando uma melhor análise e modificação do programa.
8.4- RESULTADOS DO TESTE 4
A configuração geral do teste é descrita no Item 7.2.4.
Os valores utilizados na simulação foram:
largura total do perfil: B = 20cm;
altura total do perfil: H = 20cm;
carga constante no topo e no pé do perfil:
= -23,95cm;
concentração de soluto na entrada pelo lado direito (x > 0): C 0 = 0;
113
concentração de soluto na entrada pelo lado esquerdo (x < 0): C0 = 1mg/l; e
concentração de soluto na entrada no centro (x = 0): C0 = 0,5mg/l.
Os parâmetros numéricos utilizados foram:
tempo total de simulação: t = 1h;
incremento de tempo padrão: t = 0,025h; e
opção de cálculo para o fluxo: RA Distribuída.
Com esses valores e esses parâmetros, foram obtidos para o fluxo os resultados:
conteúdo de água:
= 0,4, em todo o perfil; e
velocidade de escoamento: vz = 50cm/h em todo o perfil.
Com uma dispersividade longitudinal
dispersividade transversal
T
L
= 1cm, que fornece DL = 50cm2/h, e uma
= 0,5cm, que fornece DT = 25cm2/h, o problema foi simulado
para os quatro padrões de geração e distribuição de partículas, gerando basicamente os
mesmos resultados, assim como nos testes 1 e 2. Por isso, os próximos resultados são
apresentados apenas para a opção 00.
Para avaliar o efeito do tipo de controle de reposicionamento sobre a precisão do
modelo, o problema foi então simulado para os dois tipos de controle e os resultados foram
comparados com a solução analítica da Equação (7.15).
Os resultados são mostrados nas Figuras 8.14 e 8.15. Nos tempos de 0,25 e 0,5 hora,
os erros relativos máximos, obtidos para o controle do tipo 1, foram, respectivamente, de
13,3% e 14,8% e, para o controle do tipo 2, foram, respectivamente, de 9,5% e 13,0%.
Os resultados mostram que os dois tipos de controle tiveram precisões razoáveis e
similares, com boa aproximação dos resultados numéricos aos analíticos, ao contrário dos
testes anteriores. Observando as figuras, fica evidente que o modelo sofreu influência
consideravelmente forte da forma da malha de elementos. Essa influência ocorreu no tipo 1
com uma intensidade um pouco maior que no tipo 2. Um dos motivos da influência da malha
pode ser o esquema numérico adotado, e o fato de ser mais forte sobre o tipo 1 pode ser por
causa do número maior de reposicionamentos. Além disso, o cálculo do fluxo pode ter
resultado num campo de velocidades convergentes em algumas localidades, provocando o
avanço irregular. Nesse caso, o fluxo também estaria sofrendo o efeito da forma da malha.
Esse efeito, no entanto, não foi verificado, levantando-se, apenas, essa hipótese. Fica evidente,
portanto, a necessidade de uma melhor análise do problema com uma malha mais refinada.
114
x (cm)
-10.00 -8.00 -6.00 -4.00 -2.00
0.00
0.00
2.00
4.00
6.00
8.00
10.00
2.00
4.00
6.00
8.00
10.00
-2.00
-4.00
-6.00
z (cm)
-8.00
-10.00
-12.00
-14.00
-16.00
-18.00
-20.00
a) para 0,25 hora
x (cm)
-10.00 -8.00 -6.00 -4.00 -2.00
0.00
0.00
-2.00
-4.00
-6.00
z (cm)
-8.00
-10.00
-12.00
-14.00
-16.00
-18.00
-20.00
Analítico
Numérico
b) para 0,5 hora
Figura 8.14- Resultados para o teste 4, com controle de reposicionamento do tipo 1
115
x (cm)
-10.00 -8.00 -6.00 -4.00 -2.00
0.00
0.00
2.00
4.00
6.00
8.00
10.00
2.00
4.00
6.00
8.00
10.00
-2.00
-4.00
-6.00
z (cm)
-8.00
-10.00
-12.00
-14.00
-16.00
-18.00
-20.00
a) para 0,25 hora
x (cm)
-10.00 -8.00 -6.00 -4.00 -2.00
0.00
0.00
-2.00
-4.00
-6.00
z (cm)
-8.00
-10.00
-12.00
-14.00
-16.00
-18.00
-20.00
Analítico
Numérico
b) para 0,5 hora
Figura 8.15- Resultados para o teste 4, com controle de reposicionamento do tipo 2
116
Quanto aos efeitos da dispersão longitudinal, as Figuras 8.14 e 8.15 mostram que o
modelo simula muito bem esses efeitos, o que pode ser comprovado pelo espalhamento lateral
da contaminação e principalmente pela Figura 8.15a, onde o espalhamento calculado pelo
modelo tem concordância muito boa com a solução analítica.
8.5- RESULTADOS DO TESTE 5
A configuração geral do teste é descrita no Item 7.2.5.
Foram utilizados os seguintes valores na simulação:
altura da coluna: H = 2m;
taxa de entrada de água: qz = 20cm/h; e
concentração de soluto na entrada: C0 = 1000mg/l.
Para a simulação, foram utilizados como parâmetros numéricos:
tempo total de simulação: t = 6h;
incremento de tempo padrão: t = 0,05h; e
opção de cálculo para o fluxo: RA Distribuída.
Os principais resultados para o fluxo obtidos com os valores acima foram:
conteúdo de água:
= 0,4, em toda a coluna; e
velocidade de escoamento: vz = 50cm/h, em toda a coluna.
Considerando que, em todos os testes anteriores, realizados para uma condição
permanente de fluxo, o padrão de geração e distribuição de partículas não teve influência
significativa sobre os resultados, este teste foi realizado somente para a opção 00.
O efeito do tipo de controle de reposicionamento das partículas sobre a precisão do
modelo foi então avaliado. Para isso, foram utilizados os seguintes valores para considerar os
processos de produção e decaimento e de sorção de soluto:
constante de primeira ordem de decaimento de soluto na fase líquida:
constante de primeira ordem de decaimento de soluto na fase sólida:
termo de ordem zero de produção de soluto na fase líquida:
termo de ordem zero de produção de soluto na fase sólida:
massa específica do solo: Bd = 2500mg/cm3 ;
117
l
s
l
s
= 0,15h-1 ;
= 0,05h-1 ;
= 0,1mg/cm3 h;
= 0,003h-1;
coeficiente de distribuição: Kd = 0,0003cm3/mg; e
dispersividade longitudinal:
L
= 1,5cm (DL = 75cm2/h).
Com esses valores, foram obtidos os resultados mostrados na Figura 8.16, para o
controle de reposicionamento do tipo 1. Nos tempos de 2, 4 e 6 horas, os erros relativos
máximos obtidos para esse caso foram, respectivamente, de 24,1%, 24,5% e 24,0%.
Esses resultados mostram que o modelo se adequa muito bem à situação, calculando
de forma efetiva os processos considerados de produção, decaimento e sorção de soluto, pois
as formas das curvas analíticas se apresentam semelhantes à distribuição dos pontos
numéricos. Nota-se que nos pontos mais profundos, que não foram atingidos pelo transporte,
ocorre uma produção de soluto maior que o decaimento, com um aumento da concentração.
Nos pontos atingidos, por outro lado, verifica-se o decaimento maior. Ainda assim, ocorre um
avanço considerável da solução numérica, provocado pelo tipo de controle adotado.
Para o controle do tipo 2, os resultados obtidos são mostrados na Figura 8.17. Os erros
relativos máximos obtidos, nos tempos de 2, 4 e 6 horas, foram, respectivamente, de 3,5%,
2,3% e 2,0%.
1000
Analítico
Numérico
900
Concentração (mg/l)
800
700
600
500
2 horas
4 horas
6 horas
400
300
200
100
0
0
20
40
60
80
100
120
Profundidade (cm)
140
160
180
Figura 8.16- Resultados para o teste 5, com controle de reposicionamento do tipo 1
118
200
1000
Analítico
Numérico
900
Concentração (mg/l)
800
700
600
500
2 horas
4 horas
6 horas
400
300
200
100
0
0
20
40
60
80
100
120
Profundidade (cm)
140
160
180
200
Figura 8.17- Resultados para o teste 5, com controle de reposicionamento do tipo 2
Novamente, os resultados do tipo 2 mostram uma precisão apreciável do modelo, com
os resultados numéricos em concordância quase exata com os analíticos.
8.6- COMENTÁRIOS FINAIS
Os problemas simulados nos cinco testes apresentados permitiram avaliar a precisão e
a eficiência do modelo para diversas situações. O desempenho geral do modelo foi avaliado
pelos diversos modos de cálculo propostos para o programa.
Quanto aos padrões de geração e distribuição de partículas, os resultados apontam que,
para um regime de fluxo permanente, a opção escolhida é indiferente para a precisão, já que
não há uma variação efetiva entre os resultados. No entanto, para fluxo transiente, ocorrem
divergências apreciáveis entre os resultados, não tendo sido possível, portanto, fazer uma
análise que permitisse selecionar qual a melhor opção para o cálculo do transporte.
Quanto ao tipo de controle de reposicionamento das partículas, nos casos de regime de
fluxo permanente, os testes mostraram que há uma melhora significativa quando se usa o
controle do tipo 2, pois o controle do tipo 1 provoca sempre um avanço indesejável da frente
de contaminação numérica. Para o caso de fluxo transiente, os resultados são completamente
inconsistentes, impedindo também que seja apontado o melhor tipo de controle.
119
Quanto ao efeito do número de Peclet, os testes realizados mostram que a mudança no
regime de transporte, caracterizada pela grande variação no valor de Peclet, tem pouca ou
nenhuma influência sobre os resultados, que não sofreram muito os problemas de oscilação e
dispersão numéricas. No entanto, essa análise só foi feita nos casos de regime de fluxo
permanente. Mas, devido ao comportamento do modelo, acredita-se que o mesmo ocorra nos
casos de fluxo transiente.
Analisando de uma forma geral, para os casos de fluxo permanente, o modelo
mostrou-se bastante efetivo e com uma boa precisão no cálculo de frentes de passagem de
soluto, especialmente para o controle de reposicionamento do tipo 2. Os resultados são
amplamente satisfatórios, possibilitando o uso do modelo para diversas situações. Já para o
caso de fluxo transiente, os resultados não foram satisfatórios, sendo inclusive bastante
divergentes, impedindo uma melhor análise e necessitando mais testes.
120
9- CONCLUSÕES E RECOMENDAÇÕES
O objetivo deste trabalho foi desenvolver um modelo matemático para a simulação do
transporte de solutos em solos não-saturados. Foi então desenvolvido um modelo numérico,
denominado SSTRANS, baseado numa combinação de dois métodos numéricos: o método
dos elementos finitos e o método das características. O modelo foi implementado a partir de
um outro modelo, de simulação do fluxo em zona não-saturada, já existente, o modelo
SSFLO, desenvolvido por Koide (1990) e modificado por Campos (1998), podendo ser
considerado como uma extensão.
O desempenho do modelo foi verificado quanto à sua precisão e eficiência, utilizando
como base os padrões de geração e distribuição de partículas e os tipos de controle de
reposicionamento das mesmas. Essa verificação foi feita a partir de cinco testes realizados,
que descrevem cinco problemas apresentados na literatura, cujas soluções analíticas são
conhecidas. A precisão do modelo foi avaliada ao se comparar os resultados numéricos
produzidos pelo modelo com as soluções analíticas dos problemas. O modelo também foi
testado para avaliar o seu comportamento em relação aos problemas de oscilação e dispersão
numéricas, o que foi feito provocando-se uma variação no regime de transporte entre
predominantemente advectivo e dispersivo.
9.1- CONCLUSÕES
Embora os cinco casos testados na análise do modelo sejam problemas relativamente
simples, eles envolveram diversas situações simuláveis pelo modelo que permitiram tecer
importantes conclusões. Tais conclusões podem ser discutidas sob um aspecto considerado
preponderante, relativo ao regime de fluxo.
Para o fluxo em regime permanente, o modelo apresentou em geral uma boa precisão.
O uso dos diferentes padrões de geração e distribuição de partículas não apresentou diferenças
significativas. Com isso, pode-se afirmar que o uso de um padrão com um número menor de
partículas traz um ganho de tempo de processamento sem perda de precisão nos resultados. O
reposicionamento de partículas durante a simulação, no entanto, mostrou que um
reposicionamento freqüente é prejudicial para os resultados do modelo, provocando, na
maioria dos casos testados, um avanço considerável da frente de contaminação. Com um
reposicionamento menos freqüente, a precisão do modelo foi muito melhorada, com os
121
resultados tornando-se quase exatos. Com isso, pode-se reduzir ainda mais o tempo de
processamento ao promover um número menor de reposicionamentos, inclusive melhorando
consideravelmente a precisão nos resultados. O problema de oscilação numérica não foi
detectado e a dispersão praticamente não ocorreu, independentemente do tipo de regime de
transporte que estiver ocorrendo.
Para o regime de fluxo transiente, no entanto, o modelo não teve um comportamento
satisfatório. Os resultados para os padrões de geração e distribuição de partículas, associados
com os tipos de controle de reposicionamento, foram bastante inconsistentes, de forma que
não se pode tirar uma conclusão definitiva acerca do melhor modo de cálculo. Por isso, esse
tipo de situação exige um número maior de testes e verificações e, eventualmente,
modificações no programa do modelo. Apesar disso, mostraram que o problema da oscilação
numérica não ocorreu e que os resultados mantiveram-se razoavelmente próximos dos
esperados.
Portanto, o modelo desenvolvido produziu resultados qualitativamente bons, mas que,
quantitativamente, dependem do tipo de regime de fluxo em análise. A precisão é melhor para
um número menor de reposicionamentos das partículas, reduzindo-se o tempo de
processamento, e podendo-se inclusive utilizar o menor número possível de partículas.
9.2- RECOMENDAÇÕES PARA PESQUISAS FUTURAS
Para dar continuidade ao trabalho aqui desenvolvido ou apenas como preceitos a
serem seguidos em outros trabalhos, algumas sugestões são colocadas a seguir.
Em nenhum momento, foi verificado o balanço de massa de soluto para o transporte.
Como esse é um critério normalmente útil para estabelecer a eficiência de um modelo, alguns
testes podem ser feitos com o intuito de se verificar esse balanço. Contudo, deve ser
implementada uma rotina no programa para esse fim, seguindo, por exemplo, as formulações
apresentadas por Konikow e Bredehoeft (1978).
Somente um teste para verificar o comportamento do modelo para problemas de fluxo
transiente não foi suficiente para permitir uma análise conclusiva. Por isso, apesar das
dificuldades de se encontrar na literatura dados disponíveis, sugere-se que outros testes sejam
realizados, fornecendo uma maior base para análise. Se necessário, mudanças no programa
deverão ser feitas.
122
Deve ser ainda avaliada a precisão no cálculo do campo das velocidades de fluxo.
Como essa não foi avaliada, não foi possível determinar se os problemas ocorridos com o
fluxo transiente são resultado do cálculo do fluxo ou do cálculo das velocidades.
Com relação ao controle de reposicionamento de partículas do tipo 1, é importante que
sejam realizados testes com outros percentuais, já que o valor adotado (10%), aparentemente,
leva a reposicionamentos excessivamente freqüentes, provocando erros na precisão dos
cálculos.
Outras formas de cálculo para o transporte podem ser tentadas. A resolução da
equação de transporte pelo método explícito pode estar sujeita a erros não controláveis pelo
critério de estabilidade adotado. Pode-se experimentar o método implícito, cuja estabilidade é
incondicional ou adotar mais algum critério de estabilidade.
Alguns autores sugerem que o cálculo da movimentação das partículas tem
considerável importância na precisão dos resultados. Por isso, assim como em Yeh et al.
(1993) e Zheng (1993), essa movimentação pode ser calculada utilizando o algoritmo de
Runge-Kutta de quarta ordem. Como esse algoritmo exige um maior esforço computacional, é
mais conveniente que seja utilizado somente em regiões de onde as linhas de fluxo são curvas.
Por fim, há ainda a necessidade de se realizar testes do modelo com problemas reais,
observando o que ocorre em tais situações. Como os dados reais disponíveis na literatura não
são muitos, devem ser realizados mais experimentos em escala real para a coleta desses
dados, melhorando, com isso, a base de dados para análise deste e de outros modelos.
123
REFERÊNCIAS BIBLIOGRÁFICAS
Abbott, M.B. e Basco, D.R. (1989). Computational Fluid Dynamics: An Introduction for
Engineers. Longman Scientific & Technical, London, Inglaterra.
Ambrose, R.B., Barnwell, T.O., McCutcheon, S.C. e Williams, J.R. (1996). “Computer
models for water-quality analysis.” In: Mays, L.W. (ed.) Water Resources Handbook.
McGraw-Hill, New York, EUA, 14.1-14.28.
Amokrane, H. e Villeneuve, J.P. (1996). “A numerical method for solving the water flow
equation in unsaturated porous media.” Ground Water, 34(4), 666-674.
Banton, O., Delay, F. e Porel, G. (1997). “A new time domain Random Walk method for
solute transport in 1-D heterogeneous media.” Ground Water, 35(6), 1008-1013.
Bear, J. (1972). Dynamics of Fluids in Porous Media. American Elsevier Publishing
Company, New York, EUA.
Campos, R.F. (1998). Eficiência de Diferentes Algoritmos em Elementos Finitos para Fluxo
em Meio Saturado/Não-Saturado. Dissertação de Mestrado, publicação
MTARH.DM/009A, Departamento de Engenharia Civil, Universidade de Brasília,
Brasília, DF.
Cheng, R.T., Casulli, V. e Milford, S.N. (1984). “Eulerian-Lagrangian solution of the
convection-dispersion equation in natural coordinates.” Water Resources Research,
20(7), 944-952.
Cirilo, J.A. e Cabral, J.P. (1989). “Modelos de água subterrânea.” In: Silva, R.C.V. (ed.)
Métodos Numéricos em Recursos Hídricos. ABRH, volume 1, Rio de Janeiro, 302-380.
Constantz, J. (1982). “Temperature dependence of unsaturated hydraulic conductivity of two
soils.” Soil Science Society of America Journal, 46(3), 466-470.
Delay, F., Dzikowski, M. e de Marsily, G. (1993). “A new algorithm for representing
transport in porous media in one dimension, including convection, dispersion, and
interaction with the immobile phase with first-order kinetics.” Mathematical Geology,
25(6), 689-712.
Delay, F., Housset-Resche, H., Porel, G. e de Marsily, G. (1996). “Transport in a 2-D
saturated porous medium: a new method for particle tracking.” Mathematical Geology,
28(1), 45-71.
EPA (1994). Ground Water and Wellhead Protection. Report EPA/625/R-94/001, EPA.
Ewen, J. (1996a). “‘SAMP’ model for water and solute movement in unsaturated porous
media involving thermodynamic subsystems and moving packets: 1. Theory.” Journal
of Hydrology, 182, 175-194.
Ewen, J. (1996b). “‘SAMP’ model for water and solute movement in unsaturated porous
media involving thermodynamic subsystems and moving packets: 2. Design and
application.” Journal of Hydrology, 182, 195-207.
Faust, C.R. e Mercer, J.W. (1980a). “Ground-water modeling: numerical models.” Ground
Water, 18(4), 395-409.
Faust, C.R. e Mercer, J.W. (1980b). “Ground-water modeling: recent developments.” Ground
Water, 18(6), 569-577.
Fetter, C.W. (1993). Contaminant Hydrogeology. Prentice-Hall, New Jersey, EUA.
Foster, S. e Hirata, R. (1993). Determinação do Risco de Contaminação das Águas
Subterrâneas: Um Método Baseado em Dados Existentes. Governo do Estado de São
Paulo, Instituto Geológico, boletim n . 10, São Paulo.
Hayworth, J.S. e Burris, D.R. (1996). “Modeling cationic surfactant transport in porous
media.” Ground Water, 34(2), 274-282.
Healy, R.W. (1990). Simulation of Solute Transport in Variably Saturated Porous Media with
Supplemental Information on Modifications to the U.S. Geological Survey’s Computer
124
Program VS2D. Report 90-4025, U.S. Geological Survey, Water-Resources
Investigations.
Hillel, D. (1998). Environmental Soil Physics. Academic Press, New York, EUA.
Hossain, M.A. e Yonge, D.R. (1997). “Linear finite-element modeling of contaminant
transport in ground water.” Journal of Environmental Engineering, 123(11), 1126-1135.
Johnson, R.L., Palmer, C.D. e Fish, W. (1989).“Subsurface chemical processes.” In: EPA,
Transport and Fate of Contaminants in the Subsurface. Report EPA/625/4-89/019,
EPA, 41-56.
Keely, J.F. (1989). “Modeling subsurface contaminant transport and fate.” In: EPA, Transport
and Fate of Contaminants in the Subsurface. Report EPA/625/4-89/019, EPA, 101-131.
Koide, S. (1990). Hillslope Subsurface Flow Study by Finite Element Method. Tese de
Doutorado, University of London, London, Inglaterra.
Konikow, L.F. e Bredehoeft, J.D. (1978). Computer Model of Two-Dimensional Solute
Transport and Dispersion in Ground Water. U.S. Geological Survey, Techniques of
Water-Resources Investigations, book 7, chapter C2.
Leij, F. e Dane, J.H. (1990).“Analytical solutions of the one-dimensional advection equation
and two- or three-dimensional dispersion equation.” Water Resources Research, 26(7),
1475-1482.
Li, C.W. (1996). “Modelling variably saturated flow and solute transport into sandy soil.”
Journal of Hydrology, 186, 315-325.
Mercer, J.W. e Faust, C.R. (1980a). “Ground-water modeling: an overview.” Ground Water,
18(2), 108-115.
Mercer, J.W. e Faust, C.R. (1980b). “Ground-water modeling: mathematical models.” Ground
Water, 18(3), 212-227.
Neuman, S.P. (1973) “Saturated-unsaturated seepage by finite elements.” Journal of
Hydraulic Division of American Society of Civil Engineers, 99, 2233-2269.
Nielsen, D.R., van Genuchten, M.Th. e Biggar, J.W. (1986). “Water flow and solute transport
processes in the unsaturated zone.” Water Resources Research, 22(9), 89S-108S.
Nofziger, D.L., Rajender, K., Nayudu, S.K. e Su, P.-Y. (1989). CHEMFLO: OneDimensional Water and Chemical Movement in Unsaturated Soils. Report EPA/600/889/076, EPA.
O’Connell, P.E. e Todini, E. (1996). “Modelling of rainfall, flow and mass transport in
hydrological systems: an overview.” Journal of Hydrology, 175, 3-16.
Ogata, A. e Banks, R.B. (1961) A Solution of the Differential Equation of Longitudinal
Dispersion in Porous Media. Professional Paper 411-A, U.S. Geological Survey.
Prickett, T.A., Naymik, T.G. e Lonnquist, C.G. (1981). A ‘Random-Walk’ Solute Transport
Model for Selected Groundwater Quality Evaluations. ISWS/BUL-65/81, Bulletin 65,
Illinois Department of Energy and Natural Resources.
Rao, B.K. e Hathaway, D.L. (1989). “A three-dimensional mixing cell solute transport model
and its application.” Ground Water, 27(4), 509-516.
Richards, L.A. (1931). “Capillary conduction of liquids through porous mediums.” Physics, 1,
318-333.
Santos Neto, P.M. (1990). Compressibilidade de Solos Não-Saturados com Bolhas de Ar
Oclusas. Tese de Doutorado, COPPE/UFRJ, Rio de Janeiro.
Santos, P.C.V. (1996). Estudo da Contaminação de Água Subterrânea por Percolado de
Aterro de Resíduos Sólidos – Caso Jockey Club – DF. Dissertação de Mestrado,
publicação G.DM-32A/96, Departamento de Engenharia Civil, Universidade de
Brasília, Brasília.
Scheidegger, A.E. (1961).“General theory of dispersion in porous media.” Journal of
Geophysical Research, 66(10), 3273-3278.
125
Solley, W.B., Merk, C.F. e Pierce, R.R. (1988). Estimated Use of Water in the United States
in 1985. U.S. Geological Survey, circular 1004.
Sun, N.-Z. (1996). Mathematical Modeling of Groundwater Pollution. Springer-Verlag New
York Inc., New York, EUA.
Tucci, C.E.M. (1998).Modelos Hidrológicos. UFRGS/ABRH, Porto Alegre.
van Genuchten, M.Th. (1981).“Analytical solutions for chemical transport with simultaneous
adsorption, zero-order production and first-order decay.” Journal of Hydrology, 49,
213-233.
van Genuchten, M.Th. (1982). “A comparison of numerical solutions of the one-dimensional
unsaturated-saturated flow and mass transport equations.” Advances in Water
Resources, 5(1), 47-55.
Vauclin, M., Khanji, D. e Vachaud, G. (1979). “Experimental and numerical study of a
transient, two-dimensional unsaturated-saturated water table recharge problem.” Water
Resources Research, 15(5), 1089-1101.
Voss, C.I. (1984). A Finite-Element Simulation Model for Saturated-Unsaturated, FluidDensity-Dependent Ground-Water Flow with Energy Transport or Chemically-Reactive
Single-Species Solute Transport. Report 84-4369, U.S. Geological Survey, WaterResources Investigations.
Walton, W.C. (1991). Principles of Groundwater Engineering. Lewis Publishers, Michigan,
EUA.
Wang, H.F. e Anderson, M.P. (1982). Introduction to Groundwater Modeling. W.H. Freeman
and Company, San Francisco, EUA.
Ward, R.C. e Robinson, M. (1990). Principles of Hydrology. McGraw-Hill, 3a. ed., London,
Inglaterra.
Warrick, A.W., Biggar, J.W. e Nielsen, D.R. (1971).“Simultaneous solute and water transfer
for an unsaturated soil.” Water Resources Research, 7(5), 1216-1225.
Wrobel, L.C. (1989). “Introdução aos métodos numéricos.” In: Silva, R.C.V. (ed.) Métodos
Numéricos em Recursos Hídricos. ABRH, volume 1, Rio de Janeiro, 1-83.
Yeh, G.T. (1981).“On the computation of Darcian velocity and mass balance in the finite
element modeling of groundwater flow.” Water Resources Research, 17(5), 1529-1534.
Yeh, G.T. (1990).“A Lagrangian-Eulerian method with zoomable hidden fine-mesh approach
to solving advection-dispersion equations.” Water Resources Research, 26(6), 11331144.
Yeh, G.T. e Ward, D.S. (1979). FEMWASTE: A Finite-Element Model of Waste Transport
Through Porous Media. ORNL-5601, Oak Ridge National Laboratory.
Yeh, G.T., Sharp-Hansen, S., Lester, B., Strobl, R. e Scarbrough, J. (1992).
3DFEMWATER/3DLEWASTE: Numerical Codes for Delineating Wellhead Protection
Areas in Agricultural Regions Based on the Assimilative Capacity Criterion. Report
EPA/600/R-92/223, EPA.
Yeh, T.-C.J., Srivastava, R., Guzman, A. e Harter, T. (1993). “A numerical model for water
flow and chemical transport in variably saturated porous media.” Ground Water, 31(4),
634-644.
Yong, R.N., Mohamed, A.M.O. e Warkentin, B.P. (1996). Principles of Contaminant
Transport in Soils. Elsevier Science B.V., Amsterdam, Holanda.
Yu, F.X. e Singh, V.P. (1995). “Improved finite-element method for solute transport.”
Journal of Hydraulic Engineering, 121(2), 145-158.
Zheng, C. (1993). “Extension of the method of characteristics for simulation of solute
transport in three dimensions.” Ground Water, 31(3), 456-465.
Zienkiewics, O.C. (1977).The Finite Element Method. McGraw-Hill, London, Inglaterra.
126