ads:
Universidade de São Paulo – Escola de Engenharia de São Carlos
1
UNIVERSIDADE DE SÃO PAULO
ESCOLA DE ENGENHARIA DE SÃO CARLOS
DEPARTAMENTO DE HIDRÁULICA E SANEAMENTO
PROGRAMA DE PÓS-GRADUAÇÃO EM CIÊNCIAS DA ENGENHARIA AMBIENTAL
MARIA DE LOURDES PIMENTEL PIZARRO
Simulação de Fluxo de Água e Transporte de Solutos na Zona
Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo
São Carlos
2009
ads:
Livros Grátis
http://www.livrosgratis.com.br
Milhares de livros grátis para download.
ads:
Universidade de São Paulo – Escola de Engenharia de São Carlos
27
MARIA DE LOURDES PIMENTEL PIZARRO
Simulação de Fluxo de Água e Transporte de Solutos na Zona
Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo
Tese apresentada à Escola de Engenharia de São
Carlos da Universidade de São Paulo como parte dos
requisitos para a obtenção do Título de Doutor em
Ciências da Engenharia Ambiental
Área de concentração: Ciências da Engenharia
Ambiental
Orientador: Prof. Assoc. Edson Cezar Wendland
São Carlos
2009
28
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Autorizo a reprodução e divulgação total ou parcial deste trabalho, por qualquer meio
convencional ou eletrônico, para fins de estudo e pesquisa, desde que citada a fonte.
Ficha catalográfica preparada pela Seção de Tratamento
da Informação do Serviço de Biblioteca – EESC/USP
Pizarro, Maria de Lourdes Pimentel
P695s Simulação de fluxo de água e transporte de solutos na
zona não-saturada do solo pelo método de elementos
finitos adaptativo / Maria de Lourdes Pimentel Pizarro ;
orientador Edson Cezar Wendland. –- São Carlos, 2009.
Tese (Doutorado-Programa de Pós-Graduação e Área de
Concentração em Ciências da Engenharia Ambiental) –-
Escola de Engenharia de São Carlos da Universidade de São
Paulo, 2009.
1. Zona não-saturada. 2. Equação de Richards.
3. Equação de advecção-dispersão. 4. Modelo numérico.
5. Chorume. 6. Aterro sanitário. 7. Método de elementos
finitos. I. Título.
Universidade de São Paulo – Escola de Engenharia de São Carlos
29
30
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Universidade de São Paulo – Escola de Engenharia de São Carlos
31
A meus pais, Maria e Vicente, a meu marido
Miguel e aos meus filhos Patrícia e Eduardo
Dedico
32
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Universidade de São Paulo – Escola de Engenharia de São Carlos
33
AGRADECIMENTOS
A
Deus.
Prof. Dr. Edson Cezar Wendland, da Escola de Engenharia de São Carlos – USP, toda
minha admiração, gratidão e reconhecimento na orientação deste trabalho.
Minha família que, com dedicação e carinho, soube compreender e incentivar durante
esta longa jornada.
Prof. Dr. Marcelo Pereira de Souza, da Escola de Engenharia de São Carlos – USP,
pelo apoio e presteza.
Prof. Dr. Valdir Schalch, da Escola de Engenharia de São Carlos – USP, pelo apoio e
colaboração.
Prof. Dr. Evaldo L. Gaeta Espíndola, da Escola de Engenharia de São Carlos – USP,
pela colaboração.
M. Sc. Eng. Tiago Luís D. Forti, doutorando da Universidade de Campinas, pela
preciosa contribuição no desenvolvimento do código computacional.
Me. Agnaldo Monteiro Farias, doutorando da Universidade de Campinas, pela valiosa
colaboração a este trabalho.
Prof. Dr. Jarbas Honório de Miranda, doutor pela Escola Superior de Agricultura “Luiz
de Queiroz” - USP, pela valiosa contribuição a este trabalho.
Me. Alessandro de Jesus Firmiano, pela amizade, companheirismo, valiosa colaboração
e discussões que ajudaram na realização da pesquisa.
34
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Prof. Dr. Laércio Aparecido Lucas, doutor pelo Instituto de Ciências Matemáticas de
São Carlos – USP, pelas valiosas sugestões dadas ao trabalho.
Profa. Dra. Sônia de Almeida, doutora pelo Instituto de Química de Araraquara -
UNESP, pelas sugestões a este trabalho.
Me. Carlos Alberto F. Bispo, doutorando da Escola de Engenharia de São Carlos –
USP, pela colaboração.
Meus amigos da Academia da Força Aérea, que sempre estiveram presentes nos
momentos difíceis.
Meus amigos do Laboratório de Hidráulica Computacional, pelo apoio e incentivo.
Seção de Pós-Graduação, em particular a Claudete e Nelson, pelo atendimento
atencioso e cordial.
Ao Sistema de Bibliotecas da USP, em especial à equipe da Biblioteca Central e da
Matemática, indispensável na realização deste trabalho.
CRHEA, pelo apoio institucional e pela formação multidisciplinar direcionada à
proteção do meio ambiente.
Departamento de Hidráulica e Saneamento, em especial a Rose e Pavi pela atenção e
presteza.
Comissão Permanente do Magistério da Aeronáutica, COPEMA, da Academia da
Força Aérea, em particular ao Brig Ar Marco Antonio Carballo Perez, ao Cel Av Pedro de
Carvalho e Silva, ao Prof. Dr. Leomarcos Formiga, ao Prof. Me. Antonio Luiz Ferrari e a
Profa. Me. Rosângela de O. Colabone, pelo apoio e pela dispensa semanal concedida para que
este trabalho pudesse ser desenvolvido.
Agradeço
Universidade de São Paulo – Escola de Engenharia de São Carlos
35
RESUMO
PIZARRO, M. L. P. Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-
Saturada do Solo pelo Método de Elementos Finitos Adaptativo. 2009. 185 f. Tese
(Doutorado) – Escola de Engenharia de São Carlos, Universidade de São Paulo, São Carlos,
2009.
Devido aos riscos de contaminação dos recursos naturais solo e água, ao alto custo, ao tempo
e ao esforço humano nas investigações de campo, os modelos matemáticos, aliados às
técnicas numéricas e aos avanços computacionais, constituem uma ferramenta importante na
previsão do deslocamento de solutos, contribuindo assim, para o controle de alterações
ambientais. No Brasil, a modelação de fluxo e transporte de solutos na zona não-saturada é
voltada, quase que exclusivamente, aos problemas relacionados às atividades agrícolas.
Entretanto, tão importante quanto a problemática dos produtos químicos nas atividades
agrícolas é a questão de poluição e contaminação do solo e da água por chorume, gerado pelos
resíduos sólidos domiciliares. Neste trabalho, é desenvolvido e validado um modelo
computacional unidimensional para simulação de fluxo e transporte de solutos na zona não-
saturada do solo. O modelo matemático é dado pela equação diferencial parcial não-linear de
Richards, que rege o movimento de água no solo, e a equação diferencial parcial linear de
advecção-dispersão, do transporte de solutos, acompanhadas das condições iniciais e de
contorno. A Equação de Richards é dada em função do potencial matricial da água e a
Equação de Transporte de Solutos estima a evolução temporal da concentração de solutos no
perfil do solo. Devido à dificuldade de se obter soluções analíticas destas equações, são
resolvidas numericamente pelo Método de Elementos Finitos. As referidas equações são
resolvidas utilizando-se malhas uniformes inicialmente. Com a finalidade de obter simulações
mais eficientes, a um custo computacional reduzido, é empregada a adaptatividade com
refinamento h na malha de elementos finitos. A função interpolação polinomial utilizada é de
grau 2 ou maior que garante a conservação de massa. Na Equação de Richards, a derivada
temporal é aproximada por um quociente de diferença finita e é aplicado o esquema de Euler
Explícito e na Equação de Advecção-Dispersão, é aproximada por um quociente de diferença
finita, aplicando-se o esquema de Euler Implícito, devido à linearidade da equação. O sistema
operacional é o Linux Ubuntu 32 bits, o ambiente de programação é o PZ, escrito em
linguagem de programação C
++
. Na validação do modelo, utilizam-se dados disponíveis na
Literatura. Os resultados são comparados, utilizando-se malhas uniformes e malhas
adaptativas com refinamento h. Usando-se as malhas uniformes para o problema de Richards
e de transporte de potássio, o tempo de execução é de 22 minutos e a memória utilizada de
6164 Kb. Com as malhas adaptadas, o tempo de execução é de 3 minutos e 27 segundos,
consumindo 5876 Kb de memória. Houve, portanto, uma redução de 84,32% no tempo de
execução, usando-se malhas adaptativas. A utilização da função interpolação polinomial de
grau 2 ou maior e o refinamento h, permitem uma boa concordância do modelo na
comparação com soluções disponíveis na Literatura.
Palavras-chave: Zona não-saturada, Equação de Richards, Equação de Advecção-dispersão,
modelo numérico, chorume, aterro sanitário, Método de Elementos Finitos.
36
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Universidade de São Paulo – Escola de Engenharia de São Carlos
37
ABSTRACT
PIZARRO, M. L. P. Simulation of Water Flow and Solute Transport in the Unsaturated
Zone of the Soil by Adaptative Finite Element Method. 2009. 185 f. Tese (Doutorado) –
Escola de Engenharia de São Carlos, Universidade de São Paulo, São Carlos, 2009.
Due to the risks of contamination of soil and water resources, the high cost, time and human
effort in the field investigations, the mathematical models, combined with numerical
techniques and computational advances, are important tools in forecasting the movement of
solutes thereby contributing to the control of environmental alteration. In Brazil, modeling of
flow and solute transport in the unsaturated zone is focused, almost exclusively, on problems
related to agricultural activities. However, as important as the problematical of chemicals
products in agricultural activities is the issue of pollution and contamination of soil and water
by leachate, generated by municipal solid wastes. In this work, an one-dimensional
computational model for simulation of flow and solute transport in the unsaturated soil has
been developed and validated. The mathematical model is given by the Richards’s non-linear
partial differential equation, which determines the movement of water in the soil, and the
advection-dispersion linear partial differential equation, of the solute transport, together with
initial and boundary conditions. The Richards equation is a function of the water pressure
head and the Solute Transport equation estimate the temporal evolution of the solutes
concentration in the soil profile. Due to the difficulty of obtaining analytical solutions of these
equations, they are solved numerically using the Finite Element Method. The governing
equations are solved using initially a uniform mesh. In order to obtain more efficient
simulations with low computational cost, adaptativity with h refinement on the finite element
mesh is implemented. The interpolation function is of degree two or higher, assuring mass
conservation. In Richards’ equation, the temporal derivative is approximated by Euler Explicit
finite difference. For the advection-dispersion equation, due to the linearity of the equation, an
implicit finite difference scheme is used. The code is written in the programming language
C++ based on the programming environment PZ using operating system Linux Ubuntu 32 bit.
Model results are validated in comparison with data available in the Literature. The results are
evaluated using uniform meshes and with h refinement adaptive mesh. Using the uniform
meshes for the problem of Richards and transport of potassium, the running time is 22
minutes and 6164 Kb of memory is used. With the adapted meshes, the execution time is 3
minutes and 27 seconds, consuming 5,876 Kb of memory. Therefore there was a reduction of
84.32% in execution time, using adaptive meshes. The interpolation function with degree two
or higher and the h refinement, with reduction of the computation time, showed a good
agreement in comparison with the Literature.
Keywords: Unsaturated Zone, Richards Equation, Advection-dispersion Equation, numerical
modeling, leachate, sanitary landfill, Finite Element Method
38
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Universidade de São Paulo – Escola de Engenharia de São Carlos
39
LISTA DE FIGURAS
Figura 1 -
Elemento de solo não-saturado.................................................................. 34
Figura 2 - Esquema do experimento montado por Darcy........................................... 35
Figura 3 - Etapas da modelação.................................................................................. 48
Figura 4- Aproximação por diferenças finitas para a derivada em relação ao tempo
do potencial hidráulico...............................................................................
56
Figura 5 - Discretização de domínios: A, uma dimensão; B, duas dimensões; C,
três dimensões............................................................................................
64
Figura 6 - Solução aproximada para elemento unidimensional com dois nós............
67
Figura 7 - Funções interpolação linear para elemento unidimensional com dois nós 68
Figura 8 - Função teste para o nó i no Método de Galerkin........................................
69
Figura 9- Níveis de refinamento da malha................................................................. 110
Figura 10 - Processo de desrefinamento. a) Dois elementos filhos; b) Elemento
grosso linear; c) Elemento grosso quadrático............................................
112
Figura 11 - a) Solução aproximada em uma malha uniforme refinada com
2
6
elementos e p = 1; b) ilustração dos nós da malha uniforme refinada
com 2
6
elementos e p = 1...........................................................................
122
Figura 12 - Comparação da solução obtida com a malha uniforme fina
(
6
2
elementos e p = 1) e a solução de Celia, Bouloutas e Zarba (1990)....
123
Figura 13 -
a) Solução aproximada em uma malha adaptada com p=1 e ε = 0,01; b)
ilustração dos nós da malha no último passo de tempo da simulação........
124
Figura 14 -
a) Solução aproximada em uma malha adaptada com p=1 e ε = 0,001; b)
ilustração dos nós da malha no último passo de tempo da simulação........
125
Figura 15 -
a) Malha uniforme fina com 2
6
elementos e p = 1 versus malha adaptada
com p = 1 e ε = 0,01; b) ilustração dos nós da malha uniforme fina com
2
6
elementos e p = 1 versus malha adaptada com p = 1 e ε = 0,01............
127
Figura 16 -
a) Malha uniforme fina com 2
6
elementos e p = 1 versus malha adaptada
com p = 1 e ε = 0,001; b) ilustração dos nós da malha uniforme fina com
2
6
elementos e p = 1 versus malha adaptada com p = 1 e ε = 0,001..........
128
40
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Figura 17 -
Comparação da solução obtida com a malha adaptada (p = 1 e ε =
0,001) e a solução de Celia, Bouloutas e Zarba (1990)..............................
129
Figura 18 -
a) Solução aproximada em uma malha adaptada com p = 2 e ε = 0,1; b)
ilustração dos nós da malha no último passo de tempo da simulação........
130
Figura 19 -
a) Solução aproximada em uma malha adaptada com p = 2 e ε = 0,01; b)
ilustração dos nós da malha no último passo de tempo da simulação........
131
Figura 20 -
a) Solução aproximada em uma malha adaptada com p = 2 e ε = 0,001;
b) ilustração dos nós da malha no último passo de tempo da simulação...
132
Figura 21 -
a) Malha uniforme fina com 2
6
elementos e p = 1 versus malha adaptada
com p = 2 e ε = 0,01; b) ilustração dos nós da malha uniforme fina com
2
6
elementos e p = 1 versus malha adaptada com p = 2 e ε = 0,01............
133
Figura 22 -
a) Malha uniforme fina com 2
6
elementos e p = 1 versus malha adaptada
com p = 2 e ε = 0,001; b) ilustração dos nós da malha uniforme fina com
2
6
elementos e p = 1 versus malha adaptada com p = 2 e ε = 0,001..........
134
Figura 23 -
Comparação da solução obtida com a malha adaptada (p = 2 e ε = 0,01)
e a solução de Celia, Bouloutas e Zarba (1990).........................................
136
Figura 24 -
Comparação da solução obtida com a malha adaptada (p = 2 e ε = 0,01)
e a solução de Prasad, Kumar e Sekhar (2001)..........................................
137
Figura 25 -
Comparação da solução obtida com a malha adaptada (p = 2 e ε = 0,01)
e as soluções de Celia, Bouloutas e Zarba (1990) e de Prasad, Kumar e
Sekhar (2001).............................................................................................
137
Figura 26 - Número de elementos da malha adaptada da Equação de Richards na
iteração 1000, 2000, 3000,..., 450000........................................................
140
Figura 27 -
Potencial Matricial ψ no instante final da simulação: passo de tempo
450000; 1,75 h; a) malha uniforme com
7
2
elementos e p = 2; b) malha
adaptada com p = 2 e
ε
= 0,01....................................................................
141
Figura 28
-
Potencial Matricial
ψ
em vários instantes da simulação; a) malha
uniforme com
7
2
elementos e p = 2; b) malha adaptada com p = 2 e
ε
=
0,01.............................................................................................................
142
Figura 29
-
Umidade
θ
no instante final da simulação (passo de tempo 450000;
1.75h); a) malha uniforme com
7
2
elementos e p = 2; b) malha adaptada
com p = 2 e
ε
= 0,01...................................................................................
144
Universidade de São Paulo – Escola de Engenharia de São Carlos
41
Figura 30
-
Umidade
θ
em vários instantes simulados; a) malha uniforme
com
7
2
elementos e p = 2; b) malha adaptada com p = 2 e
ε
= 0,01...........
145
Figura 31
-
Comparação do perfil de umidade obtido por Miranda et al. (2005)
(Observado experimentalmente e simulado) e a solução obtida, no
presente trabalho, com malha uniforme (
7
2
elementos e p = 2).................
146
Figura 32
-
Comparação do perfil de umidade obtido por Miranda et al. (2005)
(Observado experimentalmente e simulado) e a solução obtida, no
presente trabalho, com malha adaptada (p = 2 e
ε
= 0,01).........................
146
Figura 33
-
Número de elementos da malha adaptada do transporte de solutos nos
passos de tempo 1, 2, 3,...,450....................................................................
148
Figura 34
-
θ
C
no instante final da simulação (passo de tempo 450; 1,75 h); a)
malha uniforme com
6
2
elementos e p = 3; b) malha adaptada com p = 3
e
ε
= 0,01....................................................................................................
149
Figura 35
-
θ
C
em vários instantes de simulação. a) malha uniforme com
6
2
elementos e p = 3; b) malha adaptada com p = 3 e
ε
= 0,01..................
150
Figura 36
-
Concentração
C
no instante final da simulação (passo de tempo 450;
1,75 h). a) malha uniforme com
6
2
elementos e p = 3; b) malha adaptada
com p = 3 e
ε
= 0,01...................................................................................
152
Figura 37 -
Concentração
C
em vários instantes de simulação. a) malha uniforme
com
6
2
elementos e p = 3; b) malha adaptada com p = 3 e
ε
= 0,01..........
153
Figura 38 -
Comparação do perfil de concentração do íon potássio obtido por
Miranda et al. (2005) (Observado experimentalmente e simulado) e a
solução obtida, no presente trabalho, com malha uniforme (
6
2
elementos
e p = 3).......................................................................................................
154
Figura 39 -
Comparação do perfil de concentração do íon potássio obtido por
Miranda et al. (2005) (Observado experimentalmente e simulado) e a
solução obtida, no presente trabalho, com malha adaptada (p = 3 e
ε
=
0,01)...........................................................................................................
154
Figura 40
-
Curva de Regressão Linear considerando o potencial matricial na Malha
Uniforme Fina e no Método de Celia, Bouloutas e Zarba.........................
157
42
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Figura 41
-
Curva de Regressão Linear com os potenciais matriciais obtidos pelo
Modelo Proposto (Malha Adaptada) e pelo Método de Celia, Boulotas e
Zarba...........................................................................................................
159
Figura 42
-
Curva de Regressão Linear com os potenciais matriciais obtidos pelo
Modelo Proposto (Malha Adaptada) e pelo Método de Prasad, Kumar e
Sekhar.........................................................................................................
160
Figura 43
-
Curva de Regressão Linear considerando a umidade obtida pelo Modelo
Proposto (Malha Adaptada) com a umidade obtida experimentalmente
por Miranda et al........................................................................................
161
Figura 44
-
Curva de Regressão Linear considerando a umidade obtida pelo Modelo
Proposto (Malha Adaptada) com a umidade obtida por Miranda et al.
(MIDI) .......................................................................................................
162
Figura 45
-
Curva de Regressão Linear da concentração de potássio do Modelo
Proposto (Malha Adaptada) com a concentração de potássio obtida pelo
Método de Miranda et al. (Observada experimentalmente).......................
163
Figura 46
-
Curva de Regressão Linear da concentração de potássio do Modelo
Proposto (Malha Adaptada) com a concentração de potássio obtida pelo
Método de Miranda et al. (MIDI)...............................................................
164
Universidade de São Paulo – Escola de Engenharia de São Carlos
43
LISTA DE TABELAS
Tabela 1
-
Condutividade hidráulica de alguns materiais............................................. 37
Tabela 2 -
Tempo de processamento e consumo de memória na solução da Equação
de Richards...................................................................................................
139
Tabela 3 -
Tempo de processamento e consumo de memória na solução das
Equações de Richards e de Transporte de Solutos.......................................
156
Tabela 4 -
Potenciais Matriciais versus Profundidade dos Métodos a serem
comparados..................................................................................................
157
Tabela 5 -
Potenciais Matriciais versus Profundidade dos Métodos a serem
comparados..................................................................................................
158
Tabela 6 -
Umidade versus Profundidade da Malha Adaptada e de Miranda et al.
(2005) (Observada experimentalmente e MIDI)................. ........................
160
Tabela 7 -
Profundidade versus Concentração do Modelo Proposto e do Método de
Miranda et al. (Observada experimentalmente e MIDI)..............................
162
44
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Universidade de São Paulo – Escola de Engenharia de São Carlos
45
LISTA DE SIGLAS
BASIC Beginner’s All-purpose Symbolic Instructional Code
BTCs Curvas breakthrough
CCA Cobre, Cromo e Arsênio
CETESB Companhia Ambiental do Estado de São Paulo
CHAIN-2D Software para escoamento não-saturado desenvolvido por Simunek e van
Genuchten (1994)
CXTFIT stands for fitting (i.e., FIT) the C(x,t) (i.e., solute concentration as a function
of position
x
and time
t
) to measurements
EDP Equações Diferenciais Parciais
EPA U.S. Environmental Protection Agency
FEMWATER Finite Element Model of Water
FILL Flow Investigation for Landfill
FULFILL Software para escoamento não-saturado desenvolvido por Noble e Arnold
(1991)
GMRes Método do Resíduo Mínimo Generalizado
HELP Hydrologic Evaluation of Landfill Performance
HYDRUS-1D Software para escoamento não-saturado desenvolvido por Simunek, Sejna e
van Genuchten
(1998)
HYDRUS-2D Software para escoamento não-saturado desenvolvido por Simunek, Sejna e
van Genuchten
(1999)
IBGE Instituto Brasileiro de Geografia e Estatística
IMASS-C Simulação de Movimento de Água e Solutos no Solo, considerando a
presença de Cultura
IMSL International Mathematical and Statistical Libraries
LDL
t
Segunda versão de Cholesky, sendo L uma matriz diagonal inferior,
D
uma
matriz diagonal e, L
t
, a transposta de L.
LRS Leachate Recirculation Sistem
LSOR Line Sucessive Over-Relaxation
LU Decomposição LU, sendo L, uma matriz triangular inferior, e U, uma matriz
triangular superior.
46
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
MATLAB Matrix Laboratory
MEF Método de Elementos Finitos
MIDI Miscible Displacemen
MODFLOW Flow Model
MT3D Modular Mass Transport 3-Dimensions
MULTIMED Multimedia Exposure Assessment Model
PELMO 1.5 Pesticide Leaching Model 1.5
PELMO 2.0 Pesticide Leaching Model 2.0
PORFLOW Software para escoamento não-saturado desenvolvido por Runchal e Sagar
(1998)
PRZM-1 Pesticide Root Zone Model - 1
PZ Ambiente de POO (programação orientada a objeto) desenvolvido por
Devloo (1997)
SEFTRAN A Simple and Efficient Flow and Transport code
SIMASS Simulação de Movimento de Água e Solutos no Solo
SOR Successive Over-Relaxation
SSOR Symmetric Successive Over-Relaxation
SSMC State–Space Mixing Cell
SUTRA The United States Geological Survey’s Saturated-Unsaturated Flow and
Transport
TPZCompEl Classe que implementa o Elemento Computacional
TPZGeoEL Classe que implementa a Geometria do Elemento
TPZGeoMesh Classe com o conjunto dos Elementos Geométricos
TPZMaterial Classe que implementa a Formulação Integral
Universidade de São Paulo – Escola de Engenharia de São Carlos
47
LISTA DE SÍMBOLOS
PZ ambiente de programação
C
++
linguagem de programação
η
porosidade [L
3
/L
3
]
θ
umidade volumétrica [L
3
/L
3
]
S
grau de saturação [ ]
V
a
volume de água do solo [L
3
]
V
v
volume de vazios ou volume de poros do solo [L
3
]
Q
vazão [L
3
/T]
K
condutividade hidráulica L/T]
A
área da seção transversal [L
2
]
(h
1
– h
2
)
diferença de pressão entre os pontos 1 e 2 [L]
L
comprimento [L]
q
fluxo de água [L/T]
t
tempo [T]
H
potencial hidráulico ou potencial total da água no solo ou carga hidráulica [L]
grad H
gradiente do potencial hidráulico
[
]
LL
R
e
número de Reynolds
k
permeabilidade intrínseca do meio poroso [L
2
]
ρ
massa específica da água [M/L
3
]
µ
viscosidade dinâmica do fluido [M/LT]
g
gravidade [L/T
2
]
z
comprimento na direção vertical ou coordenada vertical [L]
p
ψ
componente de pressão
g
ψ
componente gravitacional
os
ψ
componente osmótica
m
ψ
componente matricial
P
0
pressão que atua sobre a água padrão
C
concentração do soluto [M/L
3
]
z
v
velocidade da água nos poros na direção z [L/T]
48
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
t
∂
∂
operador de derivada parcial em relação a t
z∂
∂
operador de derivada parcial em relação a
z
d
J
fluxo de difusão do soluto [M/L
2
T]
D
d
coeficiente de difusão [L
2
/T]
z
C
∂
∂
gradiente de concentração [M/L
4
]
ºC graus Celsius
2
2
z∂
∂
operador de derivada parcial de segunda ordem em relação a z
D*
coeficiente de difusão efetiva [L
2
/T]
ϖ
coeficiente de tortuosidade do solo [ ]
α
coeficiente de dispersividade dinâmica [L]
D
coeficiente de dispersão hidrodinâmica [L
2
/T]
D
l
coeficiente de dispersão hidrodinâmica longitudinal
[L
2
/T]
D
t
coeficiente de dispersão hidrodinâmica transversal [L
2
/T]
l
α
coeficiente de dispersividade longitudinal [L]
t
α
coeficiente de dispersividade transversal [L]
d diâmetro médio do poro [L]
l
β
valor da ordem de 1,75
t
β
valor da ordem de 0,055
S
C
quantidade de soluto [ ]
K
d
coeficiente de distribuição [L
3
/M]
R
fator de retardamento
b
ρ
densidade aparente do solo [M/L
3
]
2/1
T
meia-vida
λ
coeficiente de decaimento de 1a. ordem [T
]
1−
t
0
instante inicial [T]
t
f
instante final [T]
t
∆
intervalo de tempo [T]
φ
variável dependente
Universidade de São Paulo – Escola de Engenharia de São Carlos
49
2
∇
operador laplaciano
2
2
x
∂
∂
operador de derivada parcial de segunda ordem em relação a x
2
2
y
∂
∂
operador de derivada parcial de segunda ordem em relação a y
T
temperatura [
0
C]
T
α
coeficiente de difusividade térmica do material
x
v
velocidade positiva
y
∂
∂
operador de derivada parcial em relação a y
2
2
y∂
∂
operador de derivada parcial de segunda ordem em relação a y
y
v
velocidade de propagação de
φ
(variável dependente)
τ
coeficiente de dispersão de
φ
(variável dependente)
ψ
potencial matricial [L]
0
ψ
potencial matricial na posição z = 0 [L]
1
Γ
fronteira
0
q
fluxo inicial de água [L/T]
2
Γ
fronteira
a
constante
b
constante
c
constante
3
Γ
fronteira
[
]
C
matriz
C
•
ψ
vetor de derivadas do potencial matricial nos nós
[
]
K
matriz
K
{
}
ψ
vetor potencial matricial nos nós
{
}
F
vetor
F
F
função
conhecida
1
ψ
potencial matricial no nó 1 [L]
50
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
p
ψ
potencial matricial no nó p
ε
variável
ω
variável
{
}
t
ψ
vetor potencial matricial nos nós no instante
t
{
}
tt
∆+
ψ
vetor potencial matricial nos nós no instante
t
+
t
∆
{
}
t
F
vetor
F
no instante
t
{
}
tt
F
∆+
vetor
F
no instante
t
+
t
∆
ℑ
operador diferencial
φ
ˆ
solução aproximada
i
ϕ
funções interpolação
i
α
valores das variáveis nos nós
m
número de nós da malha
(
)
zyxR
,,
erro ou resíduo da solução aproximada
(
)
zyxv
,,
função teste
Ω
domínio do problema.
e
elemento
(
)
e
φ
ˆ
valor de
φ
ˆ
, para todo elemento e
)(e
i
ϕ
funções interpolação nos elementos (uma função interpolação por nó)
n
número de nós do elemento
i
nó
j
nó
(
)
[
]
e
ϕ
matriz formada pelas funções interpolação
{
}
α
vetor
α
cujas coordenadas são os valores das variáveis em cada nó
(
)
e
i
z
coordenadas dos nós
(
)
e
j
z
coordenadas dos nós
(
)
e
L
comprimento do elemento
(
)
(
)
(
)
(
)
e
i
e
j
e
zzL −=
(
)
Ω
1
H
espaço de funções
(
)
Ω
V
espaço de funções
(
)
ΩΠ
1
subespaço de funções de elementos finitos
Universidade de São Paulo – Escola de Engenharia de São Carlos
51
(
)
ΩΠ
1
0
subespaço de funções de elementos finitos
u
variável de estado
h
parâmetro relacionado à discretização da malha no refinamento
h
p parâmetro relacionado à ordem dos polinômios no refinamento
0
→
∆
t
símbolo de tendência: intervalo de tempo tende a zero
z
∆
tamanho da variação espacial
0
→
∆
z
variação do espaço tende a zero
Pe número de Péclet
Cr número de Courant
K(θ)
condutividade hidráulica do solo não-saturado [L/T]
q
z
fluxo vertical não-saturado [L/T]
)(
ψ
s
C
capacidade hídrica específica
[
]
L1
D(θ) difusividade hidráulica não-saturada
[
]
TL.L
v velocidade da água nos poros
A’ área da seção transversal de poros ocupados pela água
θ
i
umidade volumétrica uniforme θ
i
na profundidade infinita [L
3
/L
3
]
0
θ
umidade volumétrica na superfície do solo (umidade inicial) [L
3
/L
3
]
s
ψ
potencial matricial na saturação [L]
b parâmetro de ajuste [ ]
s
θ
umidade volumétrica no ponto de saturação [L
3
L
-3
]
B parâmetro de ajuste [ ]
K
s
condutividade hidráulica do solo saturado [LT
-1
]
a parâmetro de ajuste [ ]
γ parâmetro de ajuste [ ]
θ
*
umidade volumétrica normalizada [L
3
L
-3
]
θ
r
umidade volumétrica residual do solo [L
3
L
-3
]
α
parâmetro que depende do solo [L
-1
]
n parâmetro que depende do solo [ ]
Pb chumbo
Cu cobre
Cd cádmio
As arsênio
52
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
)(z
inicial
ψ
potencial matricial no instante t = 0 [L]
L
ψ
potencial matricial na posição L [L]
S
e
saturação efetiva
m
parâmetro que depende do solo [ ]
ϑ
elemento do espaço de funções de quadrado integrável
(
)
Ω
2
L
espaço de funções de quadrado integrável
v função teste
dz
diferencial de z
1
+
n
ψ
incógnita da Equação de Richards
nf
número de funções
ψ
função aproximante da função variável de estado
v
função aproximante da função teste
zz
D
coeficiente de dispersão mecânica [L
2
/T]
)(
zC
inicial
concentração no instante t = 0 [M/L
3
]
o
C
concentração na posição z = 0 [M/L
3
]
L
C
concentração na posição z = L [M/L
3
]
θ
C
variável de estado:
produto entre a umidade e a concentração
C
θ
função aproximante da função variável de estado
[
]
1
−
K
inversa da matriz
K
ndiv
número de vezes que se divide o elemento inicial da malha
ndiv
2
número de elementos da malha
p ordem da função polinomial por partes
1
ψ
solução do primeiro passo de tempo
grad
ε
limite pré-estabelecido para desrefinamento
ε parâmetro de adaptação
n
E
du
valor da derivada da solução aproximada de
u
em cada elemento
du
max
gradiente máximo da solução aproximada
novo
u
solução no elemento recém criado
novo
e
Ω
subdomínio
Universidade de São Paulo – Escola de Engenharia de São Carlos
53
orginal
u
solução no elemento original
∑
=
nfel
j
1
somatório
nfel
número de funções
ϕ
do elemento
novo
j
α
variáveis
novo
i
ϕ
funções interpolação
orig
1
Ω
domínio de um elemento filho
orig
2
Ω
domínio de um elemento filho
original
u
1
solução original em
orig
1
Ω
original
u
2
solução original em
orig
2
Ω
1
z
limite inferior do elemento novo
2
z
limite superior do elemento novo
54
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Universidade de São Paulo – Escola de Engenharia de São Carlos
55
SUMÁRIO
1 INTRODUÇÃO............................................................................................................ 27
1.1 Objetivo Geral............................................................................................................. 30
1.1.1 Objetivos Específicos........................................................................................
30
2 FUNDAMENTAÇÃO TEÓRICA.............................................................................. 33
2.1 Grau de Saturação – Solo Saturado e Solo Não-Saturado.......................................... 33
2.2 O Movimento da Água no Solo.................................................................................. 34
2.2.1 Lei de Darcy....................................................................................................... 34
2.2.1.1 Validade da Lei de Darcy..........................................................................
36
2.2.2 Condutividade Hidráulica.................................................................................. 37
2.2.3 Gradiente Hidráulico.......................................................................................... 38
2.3 Potencial Total............................................................................................................ 38
2.4 Processos que controlam o Transporte de Solutos...................................................... 40
2.4.1 Processos Físicos de Transferência de Massa.................................................... 40
2.4.1.1 Advecção...................................................................................................
40
2.4.1.2 Difusão...................................................................................................... 41
2.4.1.3 Dispersão...................................................................................................
42
2.4.2 Processos Químicos de Transferência de Massa...............................................
44
2.4.2.1 Sorção e Fator de Retardamento............................................................... 44
2.4.2.2 Meia – Vida (
2/1
T
)....................................................................................
46
2.4.2.3 Coeficiente de Decaimento de 1ª. Ordem –
λ
...........................................
46
2.5 Fundamentos da Modelação........................................................................................
46
2.5.1 Modelos..............................................................................................................
46
2.5.2 Pontos de partida para a Modelação.................................................................. 47
2.5.3 Processos Importantes na Modelação................................................................ 48
2.6 Problemas Transientes e Estacionários....................................................................... 49
2.7 Classificação das Equações Diferenciais Parciais.......................................................
50
2.7.1 Equações Elípticas............................................................................................. 50
2.7.2 Equações Parabólicas......................................................................................... 51
2.7.3 Equações Hiperbólicas....................................................................................... 51
56
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
2.8 Problemas de Valor de Contorno................................................................................ 52
2.9 Problemas de Valor Inicial..........................................................................................
54
2.10 Métodos Numéricos.................................................................................................. 54
2.10.1 O Método de Diferenças Finitas.................................................................... 55
2.10.1.1 Solução Aproximada da Derivada Temporal....................................... 56
2.10.1.2 Consistência, Convergência e Estabilidade......................................... 58
2.10.1.2.1 Consistência............................................................................ 59
2.10.1.2.2 Estabilidade............................................................................. 59
2.10.1.2.3 Convergência.......................................................................... 61
2.10.2 O Método de Elementos Finitos.....................................................................
61
2.10.2.1 Método dos Resíduos Ponderados e Método de Galerkin................... 63
2.10.2.2 Adaptatividade..................................................................................... 69
3 EQUAÇÕES DE RICHARDS E DE ADVECÇÃO-DISPERSÃO.......................... 71
3.1 Lei de Darcy e Equação de Richards.......................................................................... 71
3.1.1 Relação entre a Velocidade de Darcy e a Velocidade da Água nos Poros...... 74
3.2 Equação de Advecção-Dispersão................................................................................ 75
4 REVISÃO BIBLIOGRÁFICA
.......................................................................
79
4.1 Breve Histórico sobre o Fluxo de Água e o Transporte de Solutos............................ 79
4.2 Fluxo de Água nos Resíduos Sólidos Domiciliares ou na Zona Não-Saturada do
Solo
......................................................................................................................
82
4.3 Transporte de Solutos em Solos Não-Saturados......................................................... 87
4.3.1 Transporte de Chorume de Aterros Sanitários
............................................
90
5 METODOLOGIA........................................................................................................ 93
5.1 Movimento da Água no Solo: Equação de Richards.................................................. 93
5.1.1 Método dos Resíduos Ponderados................................................................... 96
5.1.2 Aproximação de Elementos Finitos................................................................. 98
5.2 Transporte de Solutos: Equação de Advecção-Dispersão...........................................
99
5.2.1 Método dos Resíduos Ponderados................................................................... 101
5.2.2 Aproximação de Elementos Finitos................................................................. 103
5.3 Implementação Computacional...................................................................................
105
5.3.1 Geometria, Espaços de Interpolação e Material...............................................
106
Universidade de São Paulo – Escola de Engenharia de São Carlos
57
5.3.2 Algébrico..........................................................................................................
106
5.4 Solução Numérica
.............................................................................................
107
5.4.1 Seqüência de Simulações................................................................................. 108
5.4.2 Refinamento do Domínio................................................................................. 109
5.4.3 Numeração das Equações.................................................................................
110
5.4.4 Adaptação da Malha.........................................................................................
110
5.4.4.1 Refinamento
h
........................................................................................ 113
5.4.4.2 Projeção da solução................................................................................
114
5.4.5 Passos de Execução do Programa.................................................................... 116
6 VALIDAÇÃO, APLICAÇÕES NUMÉRICAS E ESTUDO ESTATÍSTICO........ 119
6.1 Problema de Infiltração de Água
........................................................................
120
6.1.1 Malha Uniforme – Comparação com Resultados Publicados.......................... 121
6.1.2 Adaptatividade................................................................................................. 123
6.1.3 Utilizando Ordem de Aproximação Polinomial p = 2..................................... 129
6.1.4 Discussões entre os Resultados Obtidos.......................................................... 135
6.1.5 Comparação com Celia, Bouloutas e Zarba (1990) ........................................ 136
6.1.6 Comparação com Prasad
,
Kumar e Sekhar (2001)
.....................................
136
6.2 Problema de Transporte de Soluto
......................................................................
138
6.2.1 Movimento de Água no Solo (Equação de Richards)...................................... 138
6.2.2 Transporte de Soluto (Equação de Advecção-Dispersão)................................
147
6.2.3 Tempo de Execução e Consumo de Memória..................................................
155
6.3 Estudo Estatístico comparativo do Modelo Proposto com dados da Literatura..........
156
6.3.1 Regressão Linear dos potenciais matriciais na Malha Uniforme Fina e no
Método de Celia, Bouloutas e Zarba.................................................................................
156
6.3.2 Estudo Comparativo da solução de Richards pelo Modelo Proposto e
soluções da Literatura.......................................................................................................
158
6.3.2.1 Regressão Linear dos potenciais matriciais do Modelo Proposto e do
Método de Celia, Bouloutas e Zarba.................................................................................
158
6.3.2.2
Regressão Linear dos potenciais matriciais do Modelo Proposto e do
Método de Prasad, Kumar e Sekhar.................................................................................
159
6.3.3 Estudo Comparativo da solução da Equação de Richards usando o Método
de Miranda et al. e o Modelo Proposto.............................................................................
160
58
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
6.3.3.1
Regressão Linear das umidades do Modelo Proposto e do Método de
Miranda et al. (experimental)............................................................................................
161
6.3.3.2
Regressão Linear das umidades do Modelo Proposto e do Método de
Miranda et al. (MIDI)........................................................................................................
161
6.3.4 Estudo Comparativo das Concentrações de Potássio dadas pelo Modelo
Proposto e pelo Método de Miranda et al.........................................................................
162
6.3.4.1 Regressão Linear com as concentrações do Modelo Proposto e do
Método de Miranda et al. (experimental)..........................................................................
163
6.3.4.2 Regressão Linear das concentrações do Modelo Proposto e do
Método de Miranda et al. (MIDI).....................................................................................
163
7 CONCLUSÕES……………………………………………........................................ 165
7.1
Sugestões para pesquisas futuras................................................................................ 167
8 REFERÊNCIAS BIBLIOGRÁFICAS....................................................................... 169
Universidade de São Paulo – Escola de Engenharia de São Carlos
27
1 INTRODUÇÃO
Com o crescente avanço e desenvolvimento das indústrias e com a constante geração
de resíduos sólidos nas cidades, o problema de poluição do meio ambiente vem cada vez mais
ganhando espaço e exigindo novas alternativas de controle e de prevenção de impacto
ambiental. Os resíduos gerados pela atividade humana constituem um sério problema
ambiental e de saúde pública, causando poluição e degradação do meio ambiente, com a
contaminação do solo e dos corpos hídricos, gerando odores desagradáveis e proliferando
vetores causadores de várias doenças.
A grande preocupação com a qualidade do solo e do meio ambiente em compreender
os processos de degradação e possíveis impactos ambientais torna interessante o estudo de
problemas envolvendo fluxo e transporte de solutos em solo não-saturado.
Com o estudo do transporte de solutos no solo, é possível prever os riscos de
poluição e contaminação e os impactos que determinado soluto pode causar no sistema solo-
água, a partir do conhecimento de suas propriedades, da sua interação com o meio, de sua
movimentação e persistência no solo.
Agências Governamentais e Instituições de Pesquisa têm realizado estudos de
impacto ambiental e investigações de campo, com o objetivo de determinar as concentrações
de resíduos em vários ecossistemas, mas devido ao alto custo, ao tempo e ao esforço humano
nestes trabalhos, há necessidade de uma alternativa economicamente viável.
Frente a estes aspectos e aos avanços computacionais, os modelos matemáticos,
aliados às técnicas numéricas, constituem uma ferramenta importante na prevenção de
impactos ambientais, permitindo de maneira rápida e precisa a previsão do deslocamento de
solutos.
Na simulação computacional do transporte de solutos, há necessidade da
identificação das leis que regem os mecanismos de transporte, da seleção do modelo teórico,
da determinação dos parâmetros de transporte e da resolução das equações que regem o
problema. Não deve ser desprezada a importância das observações de campo, para rever o
modelo e fazer previsões cada vez mais eficientes (QUEIROZ, 1999).
28
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
O transporte de solutos na zona não-saturada é um processo muito complexo,
principalmente quando envolve além do transporte advectivo-dispersivo, outros fenômenos
como sorção e decaimento.
A modelagem constitui-se numa tentativa aproximada de representar um fenômeno
da natureza (WENDLAND, 2004). Portanto, antes que um modelo seja aplicado a um sistema
real, necessita de ser validado, para que seja assegurado que ele, realmente, represente o
processo físico que supostamente está sendo simulado (SPERANDIO; MENDES; SILVA,
2006). Esta tarefa é feita comparando-se o resultado do modelo com observações de campo ou
de laboratório, sendo que esta comparação não precisa ser exata, mas é preciso que tenha
evidências de que os processos físicos, químicos e bioquímicos estejam adequadamente
representados e incorporados no modelo.
Desde a década de 60, com os primeiros trabalhos realizados na área de transporte de
solutos no solo (NIELSEN; BIGGAR, 1961; BIGGAR; NIELSEN, 1962 a; BIGGAR;
NIELSEN, 1962 b), o número de pesquisas nesta área tem crescido, bem como o número de
situações experimentais aplicando a Equação de Advecção-Dispersão. Porém, estudos em
campo são ainda bastante escassos e com deficiências nos métodos experimentais disponíveis,
e o que se observa é que existe grande dificuldade em envolver todos os parâmetros
pertinentes ao ambiente, necessários à simulação do problema que se está desenvolvendo. No
Brasil, a modelação de fluxo e transporte de solutos na zona não-saturada é voltada, quase que
exclusivamente, aos problemas relacionados às atividades agrícolas (MIRANDA; DUARTE;
LIBARDI, 2004). Isto se justifica devido à agricultura intensiva estar sempre em busca de
aumento de produtividade, tendo à disposição uma grande variedade de produtos químicos.
Entretanto, tão importante quanto a problemática dos produtos químicos nas atividades
agrícolas é a questão de poluição e contaminação do solo e da água por chorume, gerado pelos
resíduos sólidos domiciliares.
A partir de meados da década de 90, houve um crescimento na quantidade de
resíduos sólidos gerados, devido ao crescimento da população, ao aumento da capacidade de
consumo, adicionado ao desenvolvimento das embalagens e dos produtos descartáveis
(BIDONE; POVINELLI, 1999; IBGE, 2006).
Dificilmente a geração de resíduos será eliminada ou reduzida à zero, pois estes são
gerados pela maioria das atividades humanas, porém a busca pela sua minimização,
obedecendo a limites legais permissíveis, é de suma importância para a sustentabilidade das
cidades (NAIME, 2001).
Universidade de São Paulo – Escola de Engenharia de São Carlos
29
A solução para a disposição final dos resíduos sólidos domiciliares é a utilização de
aterros sanitários que, com base em normas operacionais específicas (CETESB, 2004) e
fundamentos de engenharia, vão permitir a confinação dos resíduos, procurando controlar a
poluição e proteger o meio ambiente. Atualmente é uma das formas mais econômicas e
adequadas de disposição de resíduos sólidos, em vários países do mundo. Para a grande
maioria dos municípios brasileiros, não será diferente, o aterro sanitário constitui a forma
mais viável de disposição final de resíduo domiciliar.
A grande quantidade de matéria orgânica presente nos resíduos sólidos domiciliares,
a umidade própria dos resíduos, a infiltração de águas de chuvas e as atividades
microbiológicas faz com que se gerem, no interior dos aterros, os líquidos percolados
denominados chorume (SCHALCH, 1984). Os aterros devem prever a coleta e o tratamento
desse efluente líquido, normalmente de cor escura, odor desagradável e elevada DBO
(Demanda Bioquímica de Oxigênio). Portanto, chorume é um líquido produzido durante a
decomposição anaeróbica dos resíduos sólidos que tem potencial impacto ao ambiente devido
a larga variedade de substâncias orgânicas e inorgânicas contidas nesta matriz.
Geralmente, o chorume pode ser caracterizado como uma solução aquosa contendo
fósforo, cloretos, sulfatos, bicarbonatos, sódio, potássio, nitrogênio amoniacal, cálcio,
magnésio, ferro, manganês e sílica. Encontram-se presentes traços de arsênio, de cádmio, de
cromo, de cobalto, de cobre, de chumbo, de mercúrio, de níquel e de zinco (SILVA, 2002;
LIMA, 2008).
O potencial poluidor do chorume depende basicamente da sua composição, que é
conseqüência das características do resíduo depositado, depende da capacidade de atenuação
do solo e das técnicas e procedimentos de disposição adotados em cada aterro (HEIDMAN;
BRUNNER, 1973).
Diante do potencial poluidor do chorume (SCHALCH, 1984; SILVA, 2002; LIMA,
2008), o modelo computacional de fluxo de água e transporte de solutos, na zona não-saturada
do solo, pode ser aplicado à simulação de chorume, num trabalho de prevenção de
contaminação de solo e de águas subterrâneas e, em particular, na construção de aterros
sanitários, auxiliando na definição do grau de impermeabilização no fundo do aterro, que se
faz necessário.
Daí a importância deste trabalho em apresentar um modelo computacional, no
ambiente de programação PZ (DEVLOO, 1997), usado como uma biblioteca, para simular o
fluxo (RICHARDS, 1931; ANDERSON; WOESSNER, 1992; BEAR; VERRUIJIT, 1993) e o
transporte (BEAR, 1979; HINDMARSH; GRESHO; GRIFFITHS, 1984) de solutos, na zona
30
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
não-saturada do solo, através das soluções numéricas das equações diferenciais que regem
estes processos. As soluções aproximadas são obtidas através do Método de Elementos
Finitos (BREBBIA; CONNOR, 1976; ZIENKIEWICZ, 1977; ZIENKIEWICZ; MORGAN,
1983; CHEN, 2002; PORTELA; CHARAFI, 2002), por se tratar de um método mais recente,
e as derivadas temporais são aproximadas pelo Método de Diferenças Finitas (WANG;
ANDERSON, 1982; MCDONALD; HARBAUGH, 1988; SPERANDIO; MENDES; SILVA,
2006)
,
aplicando-se Euler Explícito (SPERANDIO; MENDES; SILVA, 2006) na Equação de
Richards e Euler Implícito, na Equação de Transporte. As referidas equações são resolvidas
utilizando-se malhas uniformes inicialmente e, com a finalidade de obter simulações mais
eficientes, a um custo computacional reduzido, é empregada a adaptatividade com
refinamento
h
na malha de elementos finitos.
Os resultados obtidos com o desenvolvimento deste trabalho, simulação do fluxo e
transporte de solutos na zona não-saturada do solo, estarão disponíveis às equipes
multidisciplinares, em áreas afins do conhecimento humano, buscando a preservação do solo
e dos recursos hídricos em relação à poluição e contaminação.
1.1 Objetivo Geral
O objetivo geral deste projeto é o desenvolvimento e validação de um modelo
computacional aplicado para simulação de fluxo e transporte de solutos, na zona não-saturada
do solo, por meio de soluções numéricas das equações diferenciais que descrevem esses
fenômenos.
1.1.1 Objetivos específicos
Os objetivos específicos são:
descrever o fluxo e o transporte de solutos na zona não-saturada do solo através
da Equação de Richards e da Equação de Advecção-Dispersão, respectivamente;
desenvolver o código computacional no ambiente de programação PZ, em
linguagem de programação C
++
e com a filosofia de orientação a objetos;
Universidade de São Paulo – Escola de Engenharia de São Carlos
31
validar as soluções das Equações de Richards e de Transporte de Solutos dadas
pelo código computacional, utilizando-se dados disponíveis na literatura; e
fornecer subsídios para a prevenção de poluição e contaminação de águas
subterrâneas e do solo, e, em particular, em projetos de aterros sanitários, caso
disponha dos parâmetros de transporte do chorume.
32
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Universidade de São Paulo – Escola de Engenharia de São Carlos
33
2 FUNDAMENTAÇÃO TEÓRICA
Neste capítulo, são relembrados alguns conceitos teóricos relacionados aos processos
de fluxo de água (seções 2.1, 2.2, e 2.3) e transporte de solutos (seção 2.4) na zona não-
saturada do solo, que serviram de base para o desenvolvimento desse trabalho.
Como no trabalho é desenvolvido um modelo computacional (seção 2.5), formado
por duas equações diferenciais parciais, é importante tratar sobre problemas transientes e
estacionários (seção 2.6), classificação de equações diferenciais parciais (seção 2.7) e
problemas de valor de contorno (seção 2.8) e de valor inicial (seção 2.9).
Devido à natureza das equações envolvidas nestes modelos de movimento de fluido,
os métodos numéricos surgem como alternativas para a solução deste tipo de problemas,
sendo o Método de Diferenças Finitas e o Método de Elementos Finitos muito utilizados. São
abordados também os conceitos de convergência, consistência e estabilidade de solução
(SPERANDIO; MENDES; SILVA, 2006) (seção 2.10).
2.1 Grau de Saturação – Solo Saturado e Solo Não-Saturado
O
termo solo refere-se à camada externa e agriculturável da superfície terrestre
(REICHARDT, 1996).
Além de propriedades físicas do solo como porosidade
η
e umidade
θ
, o grau de
saturação S do solo (REICHARDT; TIMM, 2004) é definido como
S =
/
η
θ
= V
a
/ V
v
(2.1)
em que
V
a
- volume de água do solo [L
3
]; e
V
v
- volume de vazios ou volume de poros do solo [L
3
].
34
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
S
indica, portanto, a fração do espaço poroso que é ocupado pela água, sendo
expresso através de um número adimensional ou em porcentagem (FRENDLUND;
RAHARDJO, 1993).
O grau de saturação será 100% quando
θ
=
η
, indicando que todo espaço poroso
V
v
está cheio de água. Um solo nestas condições é denominado solo saturado.
O solo não-saturado é um sistema multifásico, em que o grau de saturação é inferior
a um. De acordo com Fredlund e Morgenstern (1977), este sistema é constituído de quatro
fases: partículas de solo, água, ar e película contráctil (interface ar-água), conforme a Figura
1. Afirmaram ainda que, sob o ponto de vista comportamental, um solo não-saturado pode ser
definido como uma mistura de duas fases em equilíbrio (partículas de solo e película
contráctil) e duas fases que fluem (ar e água).
Nielsen, Genuchten e Biggar (1986) definem zona não-saturada do solo ou zona
vadosa como a região da terra limitada em seu topo pela superfície do solo e na parte inferior
pelo lençol freático.
2.2 O Movimento da Água no Solo
2.2.1 Lei de Darcy
O marco inicial dos estudos sobre as águas subterrâneas veio à tona com o estudo
experimental do engenheiro Henry Darcy, em Dijon, na França, em 1856. Através de um
conduto preenchido de material granular, Darcy simulou o transporte da água no meio poroso,
Figura 1 – Elemento de solo não-saturado (adaptado – Fredlund e Rahardjo, 1993).
Universidade de São Paulo – Escola de Engenharia de São Carlos
35
realizando medições do volume de água transportado neste conduto. Ele realizava, ao mesmo
tempo, medições de pressão entre dois pontos. O esquema apresentado na Figura 2 mostra
como o experimento foi montado.
Figura 2 – Esquema do experimento montado por Darcy
Através deste experimento, Darcy concluiu que a vazão por este conduto é
diretamente proporcional à área de sua seção transversal, à diferença de pressão entre dois
pontos e inversamente proporcional ao comprimento da coluna do meio poroso que a mesma
percorria (REICHARDT; TIMM, 2004), ou seja,
LhhAKQ
)12
.(.
−
=
(2.2)
em que
Q
– vazão [L
3
/T];
K
– constante de proporcionalidade ou condutividade hidráulica [L/T];
A
– área da seção transversal [L
2
];
(h
1
– h
2
)
– diferença de pressão entre os pontos 1 e 2 [L]; e
L
– comprimento [L].
O fluxo de água
q
é o volume de água
a
V
que passa por unidade de tempo
t
e pela
unidade de área da seção transversal (REICHARDT; TIMM, 2004), ou seja,
36
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
A
Q
t
A
V
q
a
==
.
[
]
TL
(2.3)
em que
t
V
Q
a
=
é a vazão
[
]
TL
3
.
Comparando (2.2) e (2.3),
obtém-se a quantificação do movimento de água em
materiais porosos saturados (REICHARDT; TIMM, 2004) dada por
gradHKq
.
−
=
[L/T] (2.4)
em que
grad H
- gradiente do potencial hidráulico
[
]
LL
; e
K
- condutividade hidráulica do solo saturado
[
]
TL
.
O sinal negativo da equação (2.4) é devido ao sentido do fluxo ser inverso ao do
gradiente (REICHARDT; TIMM, 2004).
2.2.1.1 Validade da Lei de Darcy
O escoamento laminar ocorre quando as partículas do fluido movem-se em camadas
ou lâminas segundo uma trajetória paralela. Quando as partículas movem-se em direções
diversas e aleatórias, diz-se que o escoamento é turbulento.
A Lei de Darcy é bem entendida e aceita como válida nos casos de fluxo laminar,
observado em meios microporosos. Nos casos de poros maiores, o fluxo torna-se turbulento e
a velocidade não tem uma relação linear com o gradiente hidráulico, e ainda, há uma perda de
carga adicional com o aumento da velocidade e o coeficiente de proporcionalidade diminui
quando o gradiente aumenta. O limite entre o fluxo laminar e o turbulento pode ser
determinado aproximadamente através do número de Reynolds R
e
(GILES,1976) (seção 4.2).
O fluxo laminar ocorre quando R
e
< 1 e, o turbulento é observado quando R
e
> 10 (BOUWER,
Universidade de São Paulo – Escola de Engenharia de São Carlos
37
1978). O fluxo é definido como turbulento transitório quando o R
e
está entre estes dois
valores citados anteriormente. Estes valores limites apresentados são estipulados por meio de
considerações empíricas e dependem também do tipo de arranjo das partículas e da direção do
fluxo.
2.2.2 Condutividade Hidráulica
A condutividade hidráulica
K
definida inicialmente na Lei de Darcy como o
coeficiente de proporcionalidade é um parâmetro hidrogeológico, com dimensão [L/T], que
combina as propriedades do fluido com as propriedades do meio (solo). A condutividade
hidráulica é a medida da habilidade de um meio poroso de conduzir um líquido (KUTILEK;
NIELSEN, 1994).
Vários fatores influenciam a condutividade hidráulica no solo, como por exemplo, o
tamanho e a forma das partículas, espaços vazios, composição, textura, grau de saturação,
geometria dos poros e propriedades do fluido (viscosidade e massa específica) (HILLEL,
1980).
Para um solo não-saturado,
K
é função da umidade volumétrica
θ
ou da pressão da
água nos meniscos dos poros (PREVEDELLO, 1996). Os poros ocupados pelo ar reduzem a
área efetiva do fluxo, aumentando a tortuosidade do fluxo remanescente. Assim, a
condutividade hidráulica é maior quanto maior for a umidade do solo, atingindo seu valor
máximo na saturação.
Na Tabela 1 (WENDLAND, 2004), há alguns valores de condutividade hidráulica.
Tabela 1 – Condutividade hidráulica em meio saturado de alguns materiais
Material Condutividade (cm/s)
Argila
10
-9
-10
-6
Silte, silte argiloso
10
-6
-10
-4
Areia argilosa
10
-6
-10
-4
Areia siltosa, areia fina
10
-5
-10
-3
Areia bem distribuída
10
-3
-10
-1
Cascalho bem distribuído
10
-2
-10
-1
38
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Em Wendland (2004), a condutividade hidráulica pode ser expressa em função dos
dois meios que a influenciam:
µ
ρ
gkK
..
=
(2.5)
em que
K
– condutividade hidráulica [L/T];
k
– permeabilidade intrínseca do meio poroso [L
2
];
ρ
– massa específica da água [M/L
3
];
µ
– viscosidade dinâmica do fluido [M/LT]; e
g
– gravidade [L/T
2
].
2.2.3 Gradiente Hidráulico
O gradiente hidráulico
grad H
=
zH
∂
∂
/
representa a inclinação da superfície
potenciométrica, sendo
H
,
o potencial hidráulico (seção 2.3) e
z
,
o comprimento na direção
vertical. Os valores típicos do gradiente hidráulico variam de 0,0001 a 0,05 m/m (NEWELL
et al., 1996).
2.3 Potencial Total
A água do solo, como qualquer sistema ou ponto material, contém duas formas
principais de energia, a cinética e a potencial. Como o movimento da água no solo é muito
lento, a sua energia cinética é geralmente desprezível.
O fato da água do solo estar sujeita a certo número de campos de forças, faz com que
o seu potencial seja diferente do da água pura livre.
O conceito de potencial da água do solo expressa a energia potencial específica da
água do solo em relação à da água num estado de referência padrão. Para determinar essa
Universidade de São Paulo – Escola de Engenharia de São Carlos
39
diferença, uma unidade de massa, de volume ou de peso de água deve ser levada do estado
padrão para o estado no solo.
Este estado padrão é geralmente considerado como sendo o de um reservatório de
água pura, sujeito à pressão atmosférica, à mesma temperatura que a água do solo e a uma
dada altura constante.
Como a tendência espontânea da água no solo é assumir estados de menor potencial,
conhecendo-se os potenciais da água, em diferentes pontos do solo, pode-se determinar sua
tendência de movimento.
O potencial da água é composto por cinco componentes, ou seja, térmica, de pressão,
gravitacional, osmótica e matricial, mas como os processos que ocorrem no solo são
aproximadamente isotérmicos, a componente térmica é considerada desprezível. Portanto, o
potencial total da água, no solo, ou carga hidráulica, ou potencial hidráulico
H
(
PADOIN et
al., 2006) é dado por
mosgp
H
ψψψψ
+++=
(2.6)
em que
p
ψ
- componente de pressão, que surge quando a pressão que age sobre a água do
solo é diferente e maior que a pressão P
0
,
que atua sobre a água padrão. Essa pressão é
positiva;
g
ψ
- componente gravitacional, que ocorre devido à presença do campo gravitacional
terrestre;
os
ψ
- componente osmótica, que surge porque a água no solo é uma solução de sais
minerais e outros solutos e a água padrão é pura; e
m
ψ
- componente matricial, que é a soma de todos os outros trabalhos que abrangem
a interação entre a água e a matriz sólida do solo, como por exemplo, o trabalho contra forças
de adsorção e elétricas, o trabalho capilar etc. Esses fenômenos levam a água a pressões
menores que P
0
, que age sobre a água padrão. Portanto, são pressões negativas denominadas
tensões ou sucções.
40
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
2.4 Processos que controlam o Transporte de Solutos
A interação solo/soluto é extremamente complexa, pois muitos fenômenos físicos,
químicos e biológicos podem ocorrer simultaneamente.
O movimento de solutos não depende somente do fluxo do fluido no qual estas
substâncias estão dissolvidas, mas recebe influências de muitos mecanismos aos quais estas
substâncias estão sujeitas. Os
processos físicos que controlam o transporte de solutos no solo
não-saturado são advecção, difusão e dispersão. Por outro lado, os solutos nem sempre são
inertes, ou seja, os mesmos estão submetidos a outros processos como sorção, decaimento e
biodegradação. Como conseqüência dos processos de sorção, alguns solutos se movem muito
mais devagar do que a água que os está transportando, ocorrendo, portanto, o chamado
retardamento. Ao passo que a biodegradação, o decaimento radioativo e a precipitação não
reduzem necessariamente a velocidade de movimento da pluma, mas são os responsáveis por
uma diminuição da concentração do soluto na pluma (FETTER, 1992).
2.4.1 Processos Físicos de Transferência de Massa
2.4.1.1 Advecção
A advecção é o mecanismo de transporte causado pelo fluxo de água, pois com o
deslocamento da água, os solutos presentes na mesma se deslocam na direção das linhas de
fluxo com a velocidade, em princípio, igual à velocidade média linear da água e sem alterar
sua concentração na solução. Em outras palavras, se a água se move, o soluto vai ser arrastado
pelo processo de transporte físico chamado de advecção de fluxo de massa ou fluxo
convectivo.
A equação diferencial do transporte por advecção (POLYANIN; ZAITSEV, 2004;
PINCHOVER; RUBINSTEIN, 2005), no caso unidimensional, é dada por
z
C
v
t
C
z
∂
∂
−=
∂
∂
(2.7)
Universidade de São Paulo – Escola de Engenharia de São Carlos
41
em que
C
- concentração do soluto [M/L
3
];
v
z
- velocidade da água nos poros na direção z [L/T] (seção 3.1);
t
- tempo [T]; e
z -
coordenada vertical [L].
2.4.1.2
Difusão
O transporte de soluto por difusão molecular, ou simplesmente chamado de difusão,
ocorre devido ao gradiente de concentração existente em um fluido, isto é, o soluto dissolvido
em água desloca-se de uma área de maior concentração para uma área de menor concentração,
com o objetivo de igualar a concentração em toda a massa de fluido. Isto ocorre
independentemente da velocidade do fluido, mas é acentuado pela turbulência resultante dos
mecanismos de mistura mecânica (ELBACHÁ, 1989). A difusão é proporcional ao gradiente
de concentração e, em analogia com a Primeira Lei de Fick (WENDLAND, 2004), é dada por
d
J = -D
d
z
C
∂
∂
(2.8)
em que
d
J
- fluxo de difusão do soluto [M/L
2
T];
D
d
- coeficiente de difusão [L
2
/T]; e
z
C
∂
∂
- gradiente de concentração [M/L
4
].
O sinal negativo exprime que o movimento se dá em sentido contrário ao do
gradiente. Os coeficientes de difusão geralmente variam de 1 x 10
-9
a 2 x 10
-9
m
2
/s numa
temperatura de 25 ºC (SCHNOOR, 1992). Esses valores não variam muito com a
concentração, mas, sim, com a temperatura, ou seja, chegam a reduzir-se em 50% para uma
variação de 5 ºC (ROBINSON; STOKES, 1965).
No caso em que a concentração varia com o tempo, toma-se a Segunda Lei de Fick
(FREEZE; CHERRY, 1979; FETTER, 1993), ou seja,
42
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
2
2
z
C
D
t
C
d
∂
∂
=
∂
∂
(2.9)
Há dados que mostram que em solos, a difusão é consideravelmente menor que em
uma solução livre, principalmente quando se trata de solos com granulometria fina, devido à
tortuosidade
das trajetórias de fluxo (MITCHELL, 1991). Neste caso, usa-se um coeficiente
de difusão efetiva
D*
dado por
D*
=
d
D
ϖ
(2.10)
em que
ϖ
- o coeficiente de tortuosidade do solo [ ] (BEAR, 1972).
ϖ
é sempre menor que 1,0 e é determinado por meio de ensaios de laboratório.
2.4.1.3 Dispersão
Dispersão ou transporte por mistura mecânica é o processo de diluição e redução de
concentração de soluto quando o mesmo é carregado por advecção através do meio poroso
(BEAR, 1972). É decorrente da dispersão em canais individuais, do desenvolvimento de
velocidades médias diferentes em canais diferentes devido à variação das dimensões dos
poros ao longo das linhas de fluxo e do desvio da trajetória das partículas em decorrência da
tortuosidade, reentrâncias e interligações entre os canais (BEAR, 1972; HILLEL, 1980).
O coeficiente de dispersão mecânica é definido como o produto entre a velocidade da
água nos poros na direção z,
v
z
, e o coeficiente de dispersividade dinâmica α, que é a
propriedade de um meio poroso provocar a dispersão de um traçador que nele se desloca
(BEAR, 1972).
A dispersão que ocorre na direção do fluxo é denominada de dispersão longitudinal e
a que ocorre na direção perpendicular ao fluxo é chamada de dispersão transversal.
O processo de difusão molecular não pode ser separado da dispersão mecânica no
fluxo de água pelo solo. Esses dois processos são combinados para definir um parâmetro
denominado coeficiente de dispersão hidrodinâmica
D
. Portanto, o coeficiente de dispersão
Universidade de São Paulo – Escola de Engenharia de São Carlos
43
hidrodinâmica
D
é a soma do coeficiente de difusão
d
D
com o coeficiente de dispersão
mecânica
zz
D
. Os coeficientes de dispersão hidrodinâmica longitudinal
D
l
[L
2
/T] e de
dispersão hidrodinâmica transversal
D
t
[L
2
/T] (SAFFMAN, 1959; BEAR, 1972; BUNSRI;
SIVAKUMAR; HAGARE, 2008 b), quando se considera o coeficiente de difusão efetiva, são,
respectivamente,
*
DvD
zll
+=
α
(2.11)
*
DvD
ztt
+=
α
(2.12)
em que
zl
v
α
e
zt
v
α
- coeficientes de dispersão mecânica longitudinal e transversal,
respectivamente, [L
2
/T];
l
α
e
t
α
- coeficientes de dispersividade longitudinal e transversal [L],
respectivamente;
v
z
- velocidade da água nos poros na direção z [L/T]; e
D
*
- coeficiente de difusão efetiva (D
*
=
ϖ
D
d
), e
ϖ
, o coeficiente de tortuosidade
do solo [ ].
As expressões para as dispersividades (SCHEIDEGGER, 1963) são
d
ll
β
α
=
(2.13)
d
tt
β
α
=
(2.14)
em que
d
- o diâmetro médio do poro [L]; e
l
β
e
t
β
tem valores da ordem de 1,75 e 0,055, respectivamente, e que, naturalmente,
não podem ser considerados universais.
As equações (2.11) e (2.12) se tornam (SCHEIDEGGER, 1963), portanto,
*
..75,1
DvdD
zl
+=
(2.15)
44
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
*
..055,0
DvdD
zt
+=
(2.16)
O coeficiente de difusão efetiva
D
*
está na faixa de 0,005 a 0,05 m
2
/ano.
Os ensaios de laboratório mostram que para velocidades baixas de fluxo, o
coeficiente de dispersão hidrodinâmica é igual ao coeficiente de difusão efetivo, enquanto
que, para velocidades altas, o coeficiente de dispersão hidrodinâmica aumenta linearmente em
função da velocidade (PERKINS; JOHNSTON, 1963).
Quando se trata de transporte através de um solo argiloso, a difusão irá normalmente
controlar o mecanismo de transporte, e a dispersão mecânica será desprezível (GILLHAM;
CHERRY, 1982; ROWE, 1987).
2.4.2
Processos Químicos de Transferência de Massa
Dependendo do solo e da solução contaminada, muitos processos químicos poderão
ocorrer. Geralmente essas reações causam um retardamento do fenômeno de transporte dos
solutos.
Não é o objetivo a análise de todos os processos químicos e biológicos que poderiam
ocorrer no mecanismo de transporte de solutos. Serão analisados, a seguir, somente os
processos que fizeram parte da equação diferencial que rege o transporte de solutos deste
trabalho.
2.4.2.1 Sorção e Fator de Retardamento
O retardamento que ocorre durante a migração de solutos é totalmente atribuído aos
processos de sorção (adsorção-desorção).
A sorção é determinada experimentalmente pela medida da distribuição do soluto em
um sedimento particular, solo ou rochas.
Há diversos modelos que relatam a quantidade de soluto
S
C
e são conhecidos como
isotermas de sorção. A isoterma linear (VALOCCHI, 1984) é dada por
Universidade de São Paulo – Escola de Engenharia de São Carlos
45
S
C
= K
d
C
(2.17)
em que
K
d
- coeficiente de distribuição [L
3
/M]; e
C
- concentração do soluto [M/L
3
].
Comumente a sorção é quantificada geoquimicamente pelo coeficiente de
distribuição
K
d
(FREEZE; CHERRY, 1979).
A isoterma linear é apropriada quando a sorção aumenta uniformemente com o
aumento da concentração. Este modelo tem sido considerado em casos de concentrações
baixas de soluto e para sólidos com baixo potencial de sorção (WEBER Jr.; McGINLEY;
LYNN, 1991).
Campos e Elbachá (1991) definem fator de retardamento como a capacidade de
retenção ou efeito tampão do solo, para um elemento ou composto existente em um resíduo.
Para Valocchi (1984), o fator de retardamento representa a defasagem entre a velocidade de
avanço do soluto e a velocidade de avanço da frente de molhamento da solução percolante.
Desta forma, sendo o fator de retardamento um parâmetro que, indiretamente, expressa a
capacidade do solo em reter íons, fica clara sua dependência em relação às interações entre a
fase líquida e a fase sólida, durante a percolação da solução no solo.
O fator de retardamento
R
é a razão entre a velocidade do fluido percolante e a
velocidade da frente de contaminação. Seu valor poderá ser obtido da curva característica de
transporte, obtida a partir de ensaios de coluna realizados em laboratório.
O fator de retardamento
R
(TSAI et al., 2000) é dado por
θ
ρ
bd
K
R +=
1 (2.18)
em que
b
ρ
- densidade aparente do solo [M/L
3
];
θ
- umidade volumétrica [L
3
/L
3
]; e
K
d
- coeficiente de distribuição [L
3
/M].
A habilidade do solo em reter substâncias é limitada. Se a fonte de contaminação
tiver alimentação contínua, a taxa de retenção tende a diminuir com o tempo, podendo zerar
46
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
(YONG et al., 1992). Diz-se que o solo atingiu sua capacidade de retenção (
ponto de
equilíbrio).
A quantidade da substância que permanece dissolvida na água percolante aumenta
à medida que a quantidade acumulada no solo se aproxima da sua capacidade de retenção
(YONG et al., 1992).
2.4.2.2 Meia – Vida (
2/1
T
)
A meia-vida de uma reação é o tempo necessário para que a concentração chegue à
metade da concentração inicial (CHAPRA, 1997).
2.4.2.3 Coeficiente de Decaimento de 1ª. Ordem –
λ
As reações podem ser definidas como reações de ordem zero, reações de 1ª.
ordem e
reações de 2ª. ordem. Uma reação de 1ª. ordem especifica o decaimento exponencial da
concentração no tempo, onde a curva de concentração aproxima-se de forma assíntota do zero
no tempo (CHAPRA, 1997).
2.5 Fundamentos da Modelação
2.5.1 Modelos
Modelo é um sistema que consegue reproduzir, pelo menos em parte, o
comportamento de um processo natural (WENDLAND; RÜBER, 1998).
O objetivo da modelagem é o de prever ou predizer cenários onde estão envolvidas
variáveis desconhecidas, como, por exemplo, a variação da carga hidráulica ou a distribuição
de concentrações de solutos no solo, tanto no tempo quanto no espaço (BEDIENT; RIFAI;
NEWELL, 1994).
Universidade de São Paulo – Escola de Engenharia de São Carlos
47
Os modelos podem ser físicos, analógicos ou matemáticos.
O modelo físico é a representação através de um protótipo, em escala menor, na
maioria dos casos, de um sistema. Na Hidráulica, usa-se a teoria da semelhança para
estabelecer os modelos reduzidos (TUCCI, 1998).
No modelo analógico, o processo estudado é modelado no sistema mais conveniente
através da analogia entre os diferentes fenômenos. Por exemplo, representa-se o sistema
hidráulico por um circuito elétrico, que tem menor custo, baseado na analogia entre as
equações do escoamento hidráulico e as de um circuito elétrico.
Os modelos matemáticos são ferramentas matemáticas para representar uma versão
simplificada de um problema real, ou seja, através de termos matemáticos, tenta-se
compreender os processos físicos, químicos e biológicos.
2.5.2 Pontos de partida para a Modelação
No momento de planejar o modelo, algumas questões importantes devem ser
respondidas:
1. Qual é o problema a ser modelado?
2. Qual é o objetivo? Quais respostas deverão ser obtidas?
3. É necessário um modelo para se resolver o problema?
4. Quais são os dados disponíveis (dados de entrada)?
5. Os resultados do modelo podem ser verificados através de medições?
6. Quais processos serão considerados no desenvolvimento do modelo?
- fluxo de água;
- transporte de soluto;
- transporte de chorume.
Com estas questões resolvidas, realizar a construção do modelo (WENDLAND,
2004).
48
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
2.5.3 Processos Importantes na Modelação
As várias etapas para a elaboração de um modelo para reproduzir o comportamento
de um processo natural são apresentadas na Figura 3.
Problema Físico (Suposição)
Modelo Conceitual (Simplificação)
Coleta de Dados (Incertezas)
• Geologia
• Parâmetros físicos
• Condições de contorno
Modelo Físico (Simplificação) Modelo Matemático (Simplificação)
Modelo Numérico (Aproximação)
Medição
(Precisão)
Modelo Computacional (Precisão)
Interpretação (Conhecimentos básicos)
Prognóstico (Decisões)
As definições das etapas (WENDLAND, 2004) são dadas a seguir:
Problema Físico
→
Suposições: Definição de questões importantes e suposição do
processo natural;
Modelo Conceitual
→
Simplificações: Descrição dos processos envolvidos e
descrição qualitativa do comportamento do sistema natural, para que simplificações
possam ser realizadas;
Modelo Matemático
→
Simplificações: Tradução do modelo teórico em termos
matemáticos. Descrição do processo físico usando-se as relações matemáticas,
levando-se em conta a conservação de massa e de energia. Alguns fenômenos físicos
podem ser simplificados. Nessa etapa, definem-se as condições iniciais e de
Figura 3 – Etapas da modelação (WENDLAND, 2004).
Universidade de São Paulo – Escola de Engenharia de São Carlos
49
contorno. Portanto, tem-se o conjunto de equações diferenciais associadas às
condições iniciais e de contorno e pode-se resolver o problema analiticamente, caso
seja possível;
Modelo numérico
→
Aproximações: Descrição aproximada da equação matemática
diferencial. A formulação algébrica é uma aproximação da equação diferencial, com
o objetivo de determinar as variáveis (carga hidráulica, concentração, por exemplo)
em pontos discretos do modelo;
Modelo computacional
→
Erros de arredondamento: Tradução para a linguagem
computacional do modelo numérico. Nesta fase, ocorre a resolução do sistema com
diferentes técnicas numéricas, obtendo-se resultados;
Interpretação
→
Conhecimentos físicos: As grandezas calculadas são interpretadas
baseando-se nos conhecimentos físicos específicos. Ocorre a comparação dos
resultados numéricos com valores medidos, verificando-se a necessidade ou não de
calibração ou ajuste. Nesta etapa, é feita também a validação do modelo, por meio de
dados da literatura ou experimentais. Retornar à
primeira etapa, caso a solução não
seja satisfatória; e
Prognóstico
→
Decisões: Aplicação do modelo a questões específicas (por
exemplo: transporte de chorume). Após a análise de sensibilidade e com os
resultados obtidos, tomar as decisões econômicas.
2.6 Problemas Transientes e Estacionários
Na natureza, podem-se distinguir dois tipos básicos de fenômenos físicos, os
estacionários e os transientes.
Os problemas estacionários são aqueles que estão em um estado de equilíbrio, nos
quais a propriedade de interesse não se altera com o passar do tempo.
Os problemas transientes envolvem a variável temporal das grandezas físicas de
interesse, ou seja, aqueles que evoluem no tempo. A partir dos valores iniciais dessas
grandezas em certo instante
t
0
,
calculam-se, pela solução numérica da equação diferencial
parcial, quando não se dispõe da solução analítica, seus novos valores em sucessivos
intervalos de tempo
t
∆
, até alcançar o instante final
t
f
.
50
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Os problemas transientes são regidos por equações que descrevem a variação de uma
grandeza física no tempo e no espaço e os problemas estacionários são descritos por equações
que relacionam somente a distribuição espacial da grandeza.
Estes dois processos, estacionário e transiente, geralmente aparecem juntos. Os
fenômenos estacionários são representados por equações diferenciais parciais elípticas,
enquanto que os fenômenos transientes são modelados por equações diferenciais parabólicas
ou hiperbólicas. Quando apresentam mecanismos de dissipação de energia, como, por
exemplo, na difusão de calor e no escoamento de fluidos viscosos, os fenômenos ditos
dissipativos são descritos por equações parabólicas. Quando não dissipativos, são
representados por equações hiperbólicas (seção 2.7).
2.7 Classificação das Equações Diferenciais Parciais
As equações diferenciais parciais que descrevem os fenômenos de interesse da
dinâmica de fluidos computacional podem ser classificadas em três categorias básicas, ou
seja, equações elípticas, parabólicas e hiperbólicas (POLYANIN; ZAITSEV, 2004;
PINCHOVER; RUBINSTEIN, 2005).
2.7.1 Equações Elípticas
As equações elípticas, em geral, representam os problemas de equilíbrio
(estacionários), cuja equação modelo é dada pela Equação de Laplace (POLYANIN;
ZAITSEV, 2004; PINCHOVER; RUBINSTEIN, 2005). Esta equação pode ser escrita em
coordenadas cartesianas tridimensionais como
=∇
φ
2
0
2
2
2
2
2
2
=
∂
∂
+
∂
∂
+
∂
∂
zyx
φφφ
(2.19)
em que
φ
- variável dependente; e
Universidade de São Paulo – Escola de Engenharia de São Carlos
51
2
∇
- operador laplaciano.
2.7.2
Equações Parabólicas
A equação modelo para problemas parabólicos é a equação transiente de difusão de
calor (POLYANIN; ZAITSEV, 2004; PINCHOVER; RUBINSTEIN, 2005), ou seja,
2
2
x
T
t
T
T
∂
∂
=
∂
∂
α
(2.20)
em que
T
- temperatura; e
T
α
- coeficiente de difusividade térmica do material.
É importante lembrar que os mecanismos dissipativos são expressos,
matematicamente, pelo termo difusivo
2
2
x
T
∂
∂
.
Os problemas transientes necessitam de valores para a variável dependente
T
em
t
= 0 (condições iniciais) e condições de contorno para
t
> 0, ou seja, são problemas de valor
inicial (seção 2.9).
2.7.3
Equações Hiperbólicas
Equações Hiperbólicas estão relacionadas a problemas de vibração ou de convecção,
em que os fenômenos dissipativos são mínimos ou podem ser desprezados. As situações
descritas por equações hiperbólicas necessitam tanto de condições iniciais, como, em geral, de
condições de contorno. Portanto, são também problemas de valor inicial.
52
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
A equação modelo para problemas hiperbólicos é a equação de convecção
(POLYANIN; ZAITSEV, 2004; PINCHOVER; RUBINSTEIN, 2005) que, na forma
unidimensional, é escrita como
x
v
t
x
∂
∂
−=
∂
∂
φ
φ
(2.21)
em que
φ
- variável dependente; e
x
v -
velocidade positiva.
O produto
x
v
x
∂
∂
φ
é chamado de termo advectivo.
A falta de mecanismos dissipativos faz com que quaisquer descontinuidades
presentes nas condições iniciais se propaguem para a solução em
t
> 0. Isso implica que as
equações hiperbólicas admitem soluções descontínuas, e que o método numérico utilizado
deve ser capaz de lidar com elas.
O processo de advecção que possua mecanismos dissipativos obedece à Equação
Parabólica de Advecção-Dispersão, que é também chamada de Equação de Transporte
(POLYANIN; ZAITSEV, 2004; PINCHOVER; RUBINSTEIN, 2005), ou seja,
2
2
y
y
v
t
y
∂
∂
+
∂
∂
−=
∂
∂
φ
τ
φφ
(2.22)
em que
y
v
- a velocidade de propagação de
φ
;
e
τ
- o coeficiente de dispersão de
φ
.
2.8
Problemas de Valor de Contorno
Geralmente, as Equações Diferenciais Parciais (EDP) representam problemas físicos
e apresentam uma família de soluções possíveis. Porém, é necessário definir condições
Universidade de São Paulo – Escola de Engenharia de São Carlos
53
auxiliares, de modo a caracterizar unicamente a situação modelada. Estas condições são
denominadas, dependendo do problema, de condições
iniciais
e de
contorno.
Essas condições
auxiliares deverão se apresentar somente na medida certa para se obter o que se chama de
problema
bem posto,
pois se forem prescritas em excesso, poderá haver incompatibilidade
entre elas e o problema não terá solução e se forem insuficientes, o problema será indefinido,
podendo ter infinitas soluções (CUMINATO; MENEGUETE JR., 1999).
Modelos matemáticos de fluxo em solo saturado ou não-saturado, mas estacionários,
são classificados como problemas de valor de contorno. Em problemas de valor de contorno,
o analista pode especificar o valor da quantidade desconhecida ou variável de campo, como o
potencial hidráulico ou potencial matricial, ao longo das fronteiras do domínio.
Estes valores referidos como condições de contorno quando são combinados com as
equações de fluxo ou transporte resultam em modelos matemáticos que podem ser resolvidos
para valores da variável de campo em qualquer ponto do domínio. Existem três tipos de
condições de contorno. O primeiro tipo, denominado de condição de
Dirichlet
(SPERANDIO;
MENDES; SILVA, 2006), ocorre quando os valores da grandeza, que se está procurando, são
conhecidos na fronteira e podem ser usados para o cálculo dos pontos internos:
0
ψ
ψ
=
em
1
Γ
(2.23)
No segundo tipo de condição, denominado de
Neumann
(SPERANDIO; MENDES;
SILVA, 2006), apenas o gradiente (normal ou tangencial) é conhecido na fronteira,
0
qq
=
em
2
Γ
(2.24)
logo, a variável de interesse na fronteira também é uma incógnita e deve ser determinada
como parte do processo de solução. Se o gradiente normal na fronteira for nulo, diz-se que a
condição de fronteira utilizada é natural ou homogênea.
Um terceiro tipo de condição (equação 2.25), chamada de
Robin
(ou Cauchy), é uma
combinação linear dos tipos anteriores (WENDLAND, 2004).
cbaq
=
+
ψ
em
3
Γ
(2.25)
54
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
2.9 Problemas de Valor Inicial
Modelos matemáticos de fluxo de água subterrânea ou transporte de soluto em solos
saturados ou não-saturados e transientes são classificados como problemas de valor inicial.
Em problemas de valor inicial, condições de contorno, isto é, valores especificados da
variável de campo (potencial hidráulico ou potencial matricial ou concentração de soluto) e
suas derivadas (fluxo de água subterrânea ou fluxo de soluto) são especificados da mesma
maneira que os problemas de valor de contorno. Além disso, valores das variáveis de campo
precisam ser especificados para todos os pontos do domínio em algum tempo inicial
t
0
e estes
valores especificados se referem às condições iniciais. Quando as condições de contorno e as
condições iniciais são combinadas com as equações de fluxo ou transporte, resultam em um
modelo matemático que pode ser resolvido para valores da variável de campo em qualquer
ponto do domínio em qualquer tempo
t
>
t
0
.
2.10
Métodos Numéricos
Problemas envolvendo o movimento de fluidos não possuem soluções analíticas, a
não ser que se façam algumas simplificações, daí a extrema importância dos métodos
numéricos como alternativa para a solução de problemas matematicamente complexos. O
objetivo dos modelos numéricos é dar sustentação às decisões no gerenciamento dos recursos
hídricos, conhecendo a distribuição de carga hidráulica e a concentração de solutos em função
do espaço e do tempo. Diante da grande diversidade de possíveis aplicações, esforços têm
sido realizados no desenvolvimento de métodos numéricos de solução, pretendendo a
eficiência e a flexibilidade. Neste trabalho, são apresentados o Método de Diferenças Finitas e
o Método de Elementos Finitos.
Universidade de São Paulo – Escola de Engenharia de São Carlos
55
2.10.1 O
Método de Diferenças Finitas
O Método de Diferenças Finitas
é o
método mais antigo (é o método numérico dos
anos 60) e o mais divulgado, devido à simplicidade na compreensão, aprendizado e
implementação da técnica. Foi o primeiro método a ser utilizado para a solução sistemática de
problemas de água subterrânea. Os fundamentos matemáticos foram estabelecidos por Taylor
e Lagrange, matemáticos do século XVIII (GOMES, 2000).
No Método de Diferenças Finitas, as soluções de problemas de condição inicial e de
contorno são expandidas por Séries de Taylor e por meio de truncamento adequado dessas
séries, as derivadas parciais, tanto a espacial quanto a temporal, das equações diferenciais são
aproximadas por quocientes de diferenças.
Discretiza-se a região simulada em uma malha regular e para cada nó da malha é
construída uma equação algébrica. Portanto, no Método de Diferenças Finitas, uma equação
diferencial, de natureza contínua, é substituída por uma série de equações algébricas em
pontos discretos e essas equações lineares algébricas deverão ser resolvidas (McDONALD;
HARBAUGH, 1988).
Se o domínio possuir mais de uma variável, o conceito citado acima é realizado para
cada uma das variáveis separadamente (CUMINATO; MENEGUETE JR., 1999).
Geralmente, este método de aproximação apresenta-se suficientemente robusto e representa a
equação diferencial governante de forma adequada.
A vantagem desse método está na velocidade em preparar o sistema linear de
equações algébricas e na simplicidade das operações que aparecem no decorrer do processo.
As desvantagens do Método de Diferenças Finitas são a malha de discretização
precisa ser regular, a aproximação numérica do fluxo na fronteira, ou seja, a condição de
contorno do segundo tipo não é a trivial e uma determinada
singularidade não pode ser
refinada localmente, sem causar influências no domínio. Maiores detalhes sobre o Método de
Diferenças Finitas poderão ser encontrados em Morton e Mayers (2005), Rübenkönig (2006)
e Kaw e Kalu (2008).
56
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
2.10.1.1 Solução Aproximada da Derivada Temporal
Seja o sistema de equações diferenciais ordinárias cuja solução fornece valores de
ψ
e
t
∂
∂
ψ
em cada nó na malha de elementos finitos (seção 2.10.2):
[ ] [ ]
{ } { }
FKC
=+
•
ψψ
(2.26)
com
•
ψ
∂
∂
∂
∂
=
t
t
p
ψ
ψ
M
1
e
{ }
=
p
ψ
ψ
ψ
M
1
Muitos métodos existem para resolver este sistema de equações, mas o Método das
Diferenças Finitas tem sido o escolhido para o fluxo das águas subterrâneas e transporte de
solutos.
Pelo Teorema do Valor Médio do Cálculo Diferencial e Integral (Figura 4), pode-se
calcular a derivada em relação ao tempo da função
ψ
, em algum ponto
ε
, do intervalo
t
a
tt
∆
+
, ou seja,
Figura 4 – Aproximação por diferenças finitas para a derivada em relação ao tempo do potencial hidráulico.
( )
(
)
(
)
t
ttt
t
∆
−
∆
+
=
∂
∂
ψ
ψ
ε
ψ
(2.27)
t
∆
ttt
∆
+
ε
)( tt
∆
+
ψ
)(
t
ψ
t
ttt
t
∆
−
∆
+
=
∂
∂
)()(
)(
ψ
ψ
ε
ψ
Universidade de São Paulo – Escola de Engenharia de São Carlos
57
( ) ( ) ( )
t
t
ttt
∆⋅
∂
∂
+=∆+
ε
ψ
ψψ
(2.28)
Fazendo
tt
∆
+
=
ε
,
( ) ( ) ( ) ( )
t
t
t −⋅
∂
∂
+=
εε
ψ
ψεψ
(2.29)
Definindo-se a variável
ω
como
t
t
∆
−
=
ε
ω
(2.30)
pode-se escrever
( ) ( ) ( )
∂
∂
∆⋅+=
ε
ψ
ωψεψ
t
tt
( ) ( )
(
)
(
)
t
ttt
tt
∆
−
∆
+
∆⋅+=
ψ
ψ
ωψεψ
(
)
(
)
(
)
(
)
ttt
∆
+
⋅
+
−
=
ψ
ω
ψ
ω
ε
ψ
1
(2.31)
que poderá ser estendido para o vetor potencial
{
}
ψ
e o vetor
{
}
F
, ou seja,
{
}
(
)
{
}
{
}
ttt ∆+
⋅+−=
ψωψωψ
1 (2.32)
{
}
(
)
{
}
{
}
ttt
FFF
∆+
⋅+−=
ωω
1
(2.33)
Substituindo-se as equações (2.32) e (2.33) na equação (2.26) tem-se
[
]
[
]
(
)
{
}
=
∆
+
∆+ tt
KtC
ψ
ω
[
]
(
)
[
]
(
)
{
}
(
)
{
}
{
}
(
)
tttt
FFtKtC
∆+
+
−
∆
+
∆
−
−
ω
ω
ψ
ω
11
(2.34)
(
)
(
)
(
)
(
)
tttt
ψ
ω
ψ
ω
ψ
ε
ψ
⋅
−
∆
+
⋅
+
=
58
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
A solução inicia especificando-se o valor inicial de
{
}
ψ
, isto é, os valores do
potencial no tempo 0
0
=
=
tt
,
{
}
0
t
ψ
= valores especificados.
O sistema de equações lineares (equação (2.34)) é resolvido para obter valores de
{
}
ψ
em
{
}
tt ∆+
ψ
. Tem-se, portanto,
{
}
{
}
ttt ∆+
=
0
ψ
ψ
e repete-se o processo de solução para o
próximo tempo. Dependendo da escolha de
ω
, pode-se definir:
MethodDifferenceForward
→
=
0
ω
(EULER EXPLÍCITO)
[
]
{
}
[
]
[
]
(
)
{
}
{
}
tttt
FtKtCC
∆+∆−=
∆+
ψψ
(2.35)
Central→=
2
1
ω
Difference or Crank – Nicholson Method
[ ] [ ]
{ }
[ ] [ ]
{ } { } { }
( )
tttttt
FF
t
K
t
CK
t
C
∆+∆+
+
∆
+
∆
−=
∆
+
222
ψψ
(2.36)
MethodDifferenceBackward
→
=
1
ω
(EULER IMPLÍCITO)
[
]
[
]
(
)
{
}
[
]
{
}
{
}
ttttt
FtCKtC
∆+∆+
∆+=∆+
ψψ
(2.37)
2.10.1.2
Consistência, Convergência e Estabilidade
Quando se resolve uma equação diferencial parcial numericamente, a preocupação é
garantir que a solução, obtida com o esquema numérico de aproximação, represente uma
aproximação razoável da solução exata do problema matemático. Para que isso ocorra, exige-
se que a solução numérica convirja para a solução exata, quando os incrementos espacial e
temporal tendem a zero ( 0
→
∆
t
e 0
→
∆
z
). É preciso que o sistema de equações algébricas,
resultante do processo de aproximação numérica, seja consistente com a equação diferencial
parcial do problema e o método utilizado para a solução do sistema de equações seja estável
(SPERANDIO; MENDES; SILVA, 2006).
Universidade de São Paulo – Escola de Engenharia de São Carlos
59
2.10.1.2.1 Consistência
A Consistência
está relacionada com a aproximação do sistema contínuo de equações
por um sistema discreto. Um esquema numérico é consistente quando, no limite, as
aproximações numéricas se tornarem equivalentes matematicamente às equações originais, ao
se refinar as aproximações por elementos finitos ou diferenças finitas (SPERANDIO;
MENDES; SILVA, 2006).
Para se testar a consistência, substitui-se a solução exata na equação algébrica que
resultou da aproximação e expande todos os valores nodais em Séries de Taylor, em torno de
um mesmo ponto. Para que haja consistência, é necessário que a expressão resultante seja
composta da equação diferencial parcial original acrescida de um resto que tende a zero
quando a malha é refinada. Maiores detalhes poderão ser encontrados em Rolnik (1998).
2.10.1.2.2 Estabilidade
A Estabilidade
é uma propriedade relacionada, basicamente, com o esquema de
integração no tempo. Um método numérico é instável quando uma pequena perturbação,
como, por exemplo, um erro de truncamento, tende a crescer quando o processo de cálculo
avança no tempo. Esse crescimento, na maioria das vezes, é de ordem exponencial e o erro
gerado cresce acima de limites razoáveis, depois de um número pequeno de passos no
processo computacional (SPERANDIO; MENDES; SILVA, 2006).
No caso de equação de transporte, a instabilidade pode se apresentar como
concentrações negativas, concentrações que ultrapassem o valor da condição de contorno ou
um avanço da frente de concentração.
É preciso encontrar métodos de estabilização adequados para eliminar as oscilações
numéricas e, também, métodos de solução de sistemas de equações eficientes para que se
possa minimizar o tempo de simulação desse tipo de problema.
O problema de instabilidade na solução da equação de transporte é devido à presença
do termo advectivo da equação, ou seja, é devido à dispersão numérica, que está relacionada
com a propagação de erros devido a problemas de instabilidade do sistema de equações após a
discretização espacial que representa o problema de transporte de massa, podendo causar o
60
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
aparecimento de concentrações negativas ou a determinação de valores maiores que os
observados.
Para que as soluções do problema de transporte de solutos não apresentem
instabilidade devido à dispersão numérica é importante observar duas condições:
(i) na discretização espacial, define-se o número de Péclet (Pe) que é um parâmetro
adimensional usado para verificar qual mecanismo, convecção-dispersão ou difusão, domina
o processo de transferência de solutos (ROTH, 2006).
O número de Péclet é dado por:
Pe
=
D
Lv
z
(2.38)
em que
L
-
comprimento na direção z [L];
D
- coeficiente de dispersão hidrodinâmica [L
2
/T]; e
v
z
- velocidade da água nos poros, na direção
z
[L/T].
Quando o número de Péclet é alto (por exemplo,
Pe
= 200), tem-se um regime de
advecção dominante, e, quando for baixo (por exemplo,
Pe
= 2), tem-se um regime de
dispersão dominante (PINDER; GRAY, 1977).
(ii) na discretização temporal, define-se o número de Courant, que é o quociente entre
magnitude da velocidade
v
z
e o tamanho do intervalo de tempo
t
∆
pelo tamanho da
discretização espacial
z
∆
, ou seja,
Cr
=
z
tv
z
∆
∆
(2.39)
Daus, Frind e Sudicky (1985) apresentaram um estudo utilizando o Método de
Elementos Finitos para a solução da Equação de Transporte de Solutos, determinando os
valores limites para o número de Péclet e o número de Courant. Para que na solução da
Equação da Convecção-Dispersão não ocorra instabilidade é necessário que o número de
Péclet seja menor ou igual a 2 (
Pe
≤
2).
Universidade de São Paulo – Escola de Engenharia de São Carlos
61
Pe
≤
2 e
Cr
1
≈
(número de Courant próximo de 1) é critério para atingir a máxima
estabilidade numérica (PINDER; GRAY, 1977).
2.10.1.2.3 Convergência
Um método numérico aplicado a uma equação diferencial é convergente se a solução
do esquema aproximado (solução numérica) tende para a solução exata da equação
diferencial, à medida que se diminuam os incrementos espacial e temporal, desde que não
haja nenhum erro de contorno.
Como as técnicas numéricas comuns são convergentes quando aplicadas às equações
diferenciais, não significa que, na prática, quando 0
→
∆
t
e 0
→
∆
z
, a solução numérica
sempre atinja a solução exata da equação diferencial, porque erros de modelagem e de
arredondamento poderão ocorrer em qualquer problema computacional (HORNBECK, 1975).
Para equações que descrevem o escoamento de fluidos, geralmente, é impossível
demonstrar teoricamente a convergência. Em problemas que apresentam solução analítica, o
erro da solução numérica, em relação à solução exata, deve ser computado em malhas cada
vez mais refinadas. Há necessidade, portanto, de o erro da solução numérica se reduzir a zero,
à medida que
0
→
∆
t
e
0
→
∆
z
. Como as malhas devem ser extremamente refinadas, o
processo se torna caro.
2.10.2 O Método de Elementos Finitos
O Método de Elementos Finitos
é o método numérico mais usado atualmente nos
diversos ramos da Engenharia, inicialmente utilizado em problemas relacionados à
Aeronáutica, Engenharia Estrutural e Mecânica dos Solos e será utilizado neste trabalho. Foi
desenvolvido com o objetivo de resolver problemas físicos com geometria complexa e
eliminar algumas das dificuldades do Método de Diferenças Finitas. A aplicação do Método
de Elementos Finitos a problemas de água subterrânea é relativamente recente comparada
com o Método de Diferenças Finitas e foi descrita detalhadamente por Pinder e Gray (1977),
62
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Huyakorn e Pinder (1983) e Istok (1989). A flexibilidade do método é útil na solução de
transporte de contaminantes (WANG; ANDERSON, 1982).
O Método de Elementos Finitos
troca a equação diferencial do problema por uma
formulação integral (REDDY, 1986), que é resolvida, de forma aproximada, pelo Método de
Galerkin, o mais comumente usado para resolver fluxo de água e problemas de transporte de
solutos (ZIENKIEWICZ, 1977). No Método de Diferenças Finitas, as equações diferenciais
parciais são aproximadas por quocientes de diferenças, enquanto que no
Método de
Elementos Finitos
,
as formulações integrais são aproximadas por somatórios.
O domínio do problema é dividido em subdomínios, com o objetivo de representar
domínios complexos como um conjunto de subdomínios mais simples, chamados de
Elementos Finitos. Cada elemento finito é ligado aos elementos vizinhos por meio de nós. Ao
conjunto dos elementos e dos nós dá-se o nome de malha.
É importante lembrar que no Método de Elementos Finitos a
discretização pode
ocorrer com malhas não-estruturadas, nas quais a posição dos nós na malha pode ser escolhida
aleatoriamente.
A partir dos anos 70, o Método de Elementos Finitos está sendo utilizado devido à,
(às, ao):
- possibilidade da utilização de elementos irregulares (triângulos e retângulos) na geometria
da região estudada descrevendo-a melhor;
- condições de contorno poder ser facilmente agregadas ao sistema de equações;
- facilidade de implementação em regiões onde foi feito refinamento local;
- simplicidade de um programa desenvolvido ser aplicado a diferentes geometrias e processo
físicos, sem necessidade de adequação do código a cada problema.
Uma desvantagem do Método de Elementos Finitos foi apresentada por Liggett e
Liu, (1983) que é a necessidade de escolha pelo usuário de elementos de dimensão adequada,
de forma a atender critérios de estabilidade e convergência. Outra desvantagem do método é a
dificuldade em modelar meios infinitos e a exigência pelo método de grande número de dados
necessários à discretização do domínio do problema.
Em uma simulação pelo Método de Elementos Finitos, devem-se considerar os
seguintes passos básicos (HUYAKORN; PINDER, 1983; WENDLAND; RÜBER, 1998):
(1)
Discretização
: O domínio do problema é discretizado em elementos finitos, juntando-se
um número finito de nós;
Universidade de São Paulo – Escola de Engenharia de São Carlos
63
(2
) Aproximação
: Encontrar a formulação integral, sendo o Método de Resíduos Ponderados
(BEAR, 2007), o mais usado. A forma para a função teste também deverá ser especificada,
sendo o Método de Galerkin (BEAR, 2007), (seção 2.10.2.1), o mais utilizado.
(3)
Construção das matrizes por elemento
: Uma equação matricial é formada, relacionando
as variáveis nodais do elemento entre si, com base na equação diferencial parcial que descreve
o processo físico;
(4)
Construção da matriz de coeficientes global (assembling)
: As matrizes dos elementos
são agrupadas em uma matriz de coeficientes global, estabelecendo, assim, o sistema de
equações algébricas válido para todo o domínio;
(5)
Condições de contorno
: são agregadas ao sistema de equações;
(6)
Solução do sistema de equações
: O sistema de equações algébricas é resolvido pela
Eliminação de Gauss ou pelo algoritmo de Choleski ou pelo Método CG e PCG ou pelos
métodos multigrid, obtendo-se os valores das incógnitas procuradas nos nós. No caso deste
trabalho, potencial matricial
ψ
e concentração C; e
(7)
Pós-processamento
: Nesta fase, há condições de determinar grandezas derivadas das
incógnitas primárias, por exemplo, a velocidade da água nos poros.
Outras considerações sobre o Método de Elementos Finitos poderão ser encontradas
em Pinder e Gray (1977), Huyakorn e Pinder (1983), Giralt e Raviat (1986), Olsson e Heyden
(2001), Morton e Mayers (2005), Javadi, Al- Najjar e Evans (2008) e Kaw e Kalu (2008).
2.10.2.1 Método dos Resíduos Ponderados e Método de Galerkin
Considere a equação diferencial
(
)
(
)
(
)
0,,,,
=
−
ℑ
zyxFzyx
φ
(2.40)
em que
ℑ
- operador diferencial;
φ
- variável; e
F
- função conhecida.
64
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
superfície do solo
superfície do solo
O primeiro passo para a solução do fluxo de água subterrânea ou problema de
transporte de soluto pelo Método de Elementos Finitos é discretizar o domínio do problema,
substituindo-se este domínio por uma coleção de nós e elementos conhecida de malha de
elementos finitos. Os elementos consistem de dois ou mais pontos unidos por segmentos de
reta ou arco. Há diferentes tipos de elementos para problemas de uma, duas ou três dimensões
(Figura 5).
A
B
C
Figura 5 – Discretização de domínios: A, uma dimensão; B, duas dimensões; C, três dimensões.
Os elementos podem ser de qualquer tamanho. O tamanho e a forma de cada
elemento na malha podem ser diferentes e, tipos diferentes de elementos podem ser usados
numa malha.
elemento
nós
elemento
elemento
Universidade de São Paulo – Escola de Engenharia de São Carlos
65
A malha de elementos finitos é construída através de programas computacionais,
como, por exemplo, o PZ.
Quando preparada a malha de elementos finitos, é importante lembrar que a precisão
da solução obtida e o nível de esforço computacional requerido para obter a solução serão
determinados em grande parte pelo número de nós da malha. Para determinar a precisão da
solução obtida pelo Método de Elementos Finitos, basta repetirem os cálculos com a malha
refinada para ver se a troca de resultados é significativa. Por esta razão, é melhor iniciar com
uma malha constituída de poucos nós e refiná-la depois.
O segundo passo no Método de Elementos Finitos é encontrar a formulação integral
para o fluxo de água subterrânea ou equação de transporte de soluto.
O Método de Resíduos Ponderados é, em geral, o mais usado nestes casos.
No Método de Resíduos Ponderados, uma solução aproximada
φ
ˆ
para o problema
de valor inicial ou de contorno é definida, ou seja,
( ) ( )
ii
m
i
zyxzyx
αϕφ
,,,,
ˆ
1
∑
=
=
(2.41)
em que
i
ϕ
- funções interpolação;
i
α
- valores das variáveis nos nós; e
m
- número de nós da malha.
Quando a solução aproximada é substituída na equação diferencial (2.38), a equação
diferencial não é satisfeita exatamente,
(
)
(
)
(
)
(
)
0,,,,,,
ˆ
≠=−ℑ zyxRzyxFzyx
φ
(2.42)
em que
(
)
zyxR
,,
- erro ou resíduo da solução aproximada.
O erro varia ponto-a-ponto no domínio.
No Método de Resíduos Ponderados, exige-se que
66
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
(
)
(
)
∫
Ω
=Ω 0,,,, dzyxRzyxv
(2.43)
em que
(
)
zyxv
,,
- função teste; e
Ω
- o domínio do problema.
Ω
poderá ser um comprimento em um problema unidimensional, uma área em
problemas bi-dimensionais e um volume em problemas tri-dimensionais.
Para resolver a equação
( ) ( ) ( )
(
)
[
]
∫
=−ℑ
Ω
Ωφ
0dz,y,xFz,y,x
ˆ
z,y,xv (2.44)
é preciso especificar a forma matemática da solução aproximada
φ
ˆ
e a função teste v.
O valor de
φ
ˆ
, para todo elemento e,
(
)
e
φ
ˆ
, é dado por
( )
( )
( )
∑
=
=
n
i
i
e
i
e
zyx
1
,,
ˆ
αϕφ
(2.45)
em que
)(
e
i
ϕ
- funções interpolação nos elementos (uma função interpolação por nó);
i
α
- valores das variáveis em cada nó; e
n - número de nós do elemento.
Por exemplo, a solução aproximada de um elemento unidimensional com dois nós i
e j
pode ser escrita por
(
)
(
)
(
)
(
)
(
)
(
)
j
e
ji
e
i
e
zzz
αϕαϕφ
+=
ˆ
(2.46)
ou, na forma de matriz,
(
)
(
)
(
)
[
]
{
}
αϕφ
ee
z
ˆ
=
(2.47)
Universidade de São Paulo – Escola de Engenharia de São Carlos
67
em que
( ) ( )
=
z
e
j
z
e
i
e
ϕϕϕ
(2.48)
{ }
=
j
i
α
α
α
(2.49)
Figura 6 – Solução aproximada para elemento unidimensional com dois nós.
Para o elemento da Figura 6, as funções interpolação (Figura 7) são funções lineares
de z:
( )
( )
(
)
( )
e
e
j
e
i
L
zz
z
−
=
ϕ
( )
( )
(
)
( )
e
e
i
e
j
L
zz
z
−
=
ϕ
(2.50)
em que
(
)
e
i
z e
(
)
e
j
z
- coordenadas dos nós;
(
)
e
L
- comprimento do elemento
(
)
(
)
(
)
(
)
e
i
e
j
e
zzL −= .
O valor de
(
)
e
i
ϕ
é um no nó i e decresce linearmente para zero no nó j, enquanto o
valor de
(
)
e
j
ϕ
é um no nó j e decresce linearmente para zero no nó i.
No nó
(
)
(
)
e
i
zzi
= , tem-se que
)(e
i
zz
inó
=
)(e
j
zz
jnó
=
)(e
L
j
α
i
α
)(
ˆ
)(
z
e
φ
68
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
(
)
(
)
(
)
(
)
(
)
(
)
iji
e
jii
e
ii
e
zzz
ααϕαϕφ
=+=
ˆ
(2.51)
No nó
(
)
(
)
e
j
zzj =
, tem-se que
(
)
(
)
(
)
(
)
(
)
(
)
jjj
e
jij
e
ij
e
zzz
ααϕαϕφ
=+=
ˆ
(2.52)
e no ponto médio do elemento
( ) ( )
+
=
2
e
i
e
j
zz
z ,
( )
( )
( )
( )
( )
( )
2
ˆ
ji
j
e
ji
e
i
e
zzz
α
α
αϕαϕφ
+
=+=
(2.53)
Figura 7 – Funções interpolação linear para elemento unidimensional com dois nós.
A forma para a função teste v da equação (2.44) também precisa ser especificada.
Para isso será usado o Método de Galerkin, que é o mais usado em fluxo de águas
subterrâneas e problemas de transporte de soluto. Neste método, a função teste para um nó
será idêntica à função interpolação usada para definir a solução aproximada
φ
ˆ
. Portanto, para
o elemento unidimensional com dois nós,
1/2
1/2
0
1
1
0
1
0
2/1
)(e
ϕ
)(
)(
z
e
j
ϕ
)(
)(
z
e
i
ϕ
inó jnó
z
Universidade de São Paulo – Escola de Engenharia de São Carlos
69
( )
L
zz
zv
j
i
−
= para
i
zz
≥
(2.54)
( )
L
zz
zv
i
j
−
= para
i
zz
≥
(2.55)
como aparece na Figura 8.
Figura 8 - Função teste para o nó i no Método de Galerkin.
Com as formas da solução aproximada e da função teste já especificadas, pode-se
calcular a integral da equação (2.44) para obter o sistema de equações lineares
[
]
{
}
{
}
FK
=
α
(2.56)
que poderá ser resolvido em cada nó da malha.
2.10.2.2 Adaptatividade
A aproximação de elementos finitos é definida sobre um espaço de funções
aproximadas. No capítulo 5, são definidos os espaços de funções
(
)
Ω
1
H
e
(
)
Ω
V
e os
subespaços de funções de elementos finitos
(
)
(
)
Ω⊂ΩΠ
11
H
e
(
)
ΩΠ
1
0
. Os subespaços de
funções de elementos finitos
1
Π
e
1
0
Π são construídos por funções polinomiais por partes
1
0
)(zv
i
i
v
inó jnó
z
j
v
70
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
com suporte compacto. As variáveis de estado e teste são aproximadas, respectivamente, por:
∑
=
≈
nf
j
jj
u
1
ϕα
e
∑
=
≈
nf
i
i
v
1
ϕ
, sendo
i
ϕ
e
j
ϕ
funções polinomiais.
Pela definição do Método de Elementos Finitos, os parâmetros que podem ser
alterados para modificar os espaços
(
)
(
)
Ω⊂ΩΠ
11
H
e
(
)
ΩΠ
1
0
de funções de aproximação
são:
•
h: parâmetro relacionado à discretização da malha, normalmente associado ao
"tamanho" do elemento;
•
p: parâmetro relacionado à ordem dos polinômios, de cada elemento, utilizados
como base para o espaço V.
Portanto, os espaços
(
)
(
)
Ω⊂ΩΠ
11
H
e
(
)
ΩΠ
1
0
podem ser alterados aumentando-se o
número de elementos ou aumentando a ordem de polinômios da base
i
ϕ
e
j
ϕ
.
A adaptatividade consiste no enriquecimento/melhoria do espaço de funções
(
)
(
)
Ω⊂ΩΠ
11
H
e
(
)
ΩΠ
1
0
, onde este enriquecimento pode ser feito através do refinamento
dos parâmetros h, p ou ainda a combinação desses. O refinamento
h
consiste na redução ou
aumento do tamanho dos elementos da malha. Nesse método, graus de liberdade são
adicionados através da divisão de elementos existentes, onde a malha necessita ser refinada,
ou eliminam-se graus de liberdade, através do agrupamento de elementos, onde a malha pode
ser desrefinada. O refinamento p consiste na elevação da ordem dos polinômios da base de
funções teste.
Nesse trabalho, utiliza-se o refinamento
h
que consiste no refinamento do parâmetro
h
.
Universidade de São Paulo – Escola de Engenharia de São Carlos
71
3 EQUAÇÕES DE RICHARDS E DE ADVECÇÃO-DISPERSÃO
Devido à importância das Equações de Richards e de Advecção-Dispersão no
modelo desenvolvido neste trabalho, são abordadas, neste capítulo, a demonstração da
Equação de Richards, nas três formas de apresentação, a partir da adaptação da Lei de Darcy
para solos não-saturados e do Princípio da Conservação da Massa (seção 3.1), a relação da
velocidade de Darcy com a velocidade da água nos poros (seção 3.1.1), bem como a
demonstração da Equação de Advecção-Dispersão (seção 3.2).
3.1 Lei de Darcy e Equação de Richards
No caso de movimento na fase líquida, para efeito de quantificação do movimento, o
potencial total da água não inclui a componente osmótica, porque o potencial osmótico não
causa movimento significativo de água, apenas os sais se movem até o equilíbrio ser atingido.
Define-se, então, o potencial hidráulico
H
ou potencial total, sem a componente
osmótica, ou seja,
mgp
H
ψψψ
++=
(3.1)
Como
mp
e
ψ
ψ
se referem a pressões, o primeiro, às positivas e o segundo, às
negativas, podem ser agrupadas em uma única componente
mp
ψ
ψ
ψ
+
=
. A
componente
gravitacional
g
ψ
pode ser expressa em termos da altura
z
, sendo a coordenada
z
orientada
positivamente para cima. Portanto, o potencial hidráulico (STONE et al., 2006) se apresenta
como
zH
+
=
ψ
(3.2)
72
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
A água no estado líquido se move quando diferenças de potencial hidráulico
H
existem nos diferentes pontos do sistema, no sentido do decréscimo do potencial
H
.
Em Reichardt e Timm (2004), Buckingham adaptou a Lei de Darcy para solos não-
saturados, introduzindo a dependência da umidade nessa equação, formulando a Lei de
Darcy-Buckingham
,
gradHKq
).(
θ
−
=
(3.3)
em que
K(θ)
- condutividade hidráulica do solo não-saturado [LT
1−
].
Como as equações de fluxo e de transporte deste trabalho serão tratadas
unidimensionalmente e em solos não-saturados, é interessante considerar a Lei de Darcy-
Buckingham do fluxo vertical não-saturado,
q
z
, obtida substituindo-se a equação (3.2) na
(3.3),
z
z
Kq
z
∂
+
∂
−=
)(
)(
ψ
θ
)1
)(
)(( +
∂
∂
−=
z
Kq
z
ψ
θ
z
KKq
z
∂
∂
−−=
ψ
θθ
)()( (3.4)
Combinando-se a equação (3.4) com a Equação da Continuidade (WENDLAND,
2004) dada por
t
z
q
z
∂
∂
=
∂
∂
−
θ
(3.5)
obtém-se
∂
∂
∂
∂
+
∂
∂
=
∂
∂
z
K
zz
K
t
ψθ
(3.6)
Universidade de São Paulo – Escola de Engenharia de São Carlos
73
que é a forma mista da Equação de Richards, unidimensional (HE; REN, 2009).
Sabe-se
que
t
C
tt
s
∂
∂
=
∂
∂
∂
∂
=
∂
∂
ψ
ψ
ψ
ψ
θ
θ
).(.
(3.7)
em que
)(
ψ
s
C
- capacidade hídrica específica
[
]
L1
.
Substituindo-se (3.7) em (3.6), obtém-se a Equação de Richards, formulada em
termos de potencial
ψ
, unidimensional,
t
C
s
∂
∂
ψ
ψ
).(
∂
∂
∂
∂
+
∂
∂
=
z
K
zz
K
ψ
(3.8)
Tem-se que
zCzz
s
∂
∂
=
∂
∂
∂
∂
=
∂
∂
θ
ψ
θ
θ
ψ
ψ
.
)(
1
.
(3.9)
Substituindo-se (3.9)
em (3.6),
obtém-se a Equação de Richards, formulada em
termos de umidade
θ
, unidimensional,
∂
∂
∂
∂
+
∂
∂
=
∂
∂
z
D
zz
K
t
θ
θ
θ
)(
(3.10)
em que
D(θ) =
)(
)(
θ
θ
s
C
K
- difusividade hidráulica não-saturada
[
]
TL.L
.
Portanto, o fluxo de água em meio poroso não-saturado é dado pela Equação de
Richards, que expressa o Princípio da Conservação da Massa e a Lei de Darcy-Buckingam e
pode ser escrita, unidimensionalmente, nas três formas equivalentes, equações (3.6), (3.8) e
(3.10).
Os termos do membro esquerdo das equações (3.6), (3.8) e (3.10) descrevem os
efeitos da drenagem e preenchimento dos poros. Com relação aos membros direitos, a
derivada espacial da carga hidráulica descreve a difusão da água em meio poroso.
74
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
3.1.1 Relação entre a Velocidade de Darcy e a Velocidade da Água nos Poros
O fluxo de água dado pela Lei de Darcy, embora tenha as dimensões de uma
velocidade, não representa exatamente a velocidade com que a água se move no solo.
A velocidade real v da água no solo, também chamada de velocidade da água nos
poros, é dada pelo volume de água
a
V
que passa por unidade de tempo pela área da seção
transversal de poro A’ ocupados pela água, ou seja,
v =
t
A
V
a
.
'
(3.11)
Considerando-se um solo saturado, a velocidade real v da água no solo será
v =
n
q
(3.12)
Quando o solo é não-saturado, a área disponível para o fluxo será menor ainda,
porque a água só vai se mover pelos poros que estão cheios de água e, portanto,
v =
θ
q
(3.13)
Como a Lei de Darcy-Buckingham na direção z é dada pela equação (3.4), a
velocidade da água nos poros, na direção z, é dada por
v
z
= -K (θ) (
1+
∂
∂
z
ψ
)
θ
1
(3.14)
Universidade de São Paulo – Escola de Engenharia de São Carlos
75
3.2 Equação de Advecção-Dispersão
Nesta seção, é derivada a Equação de Transporte de Soluto na zona não-saturada do
solo, na forma unidimensional, utilizada neste trabalho. O transporte de solutos ocorre por três
processos, ou seja, advecção, difusão e dispersão mecânica (seção 2.4).
O fluxo do transporte de soluto por advecção é dado por
CqJ
zm
=
(3.15)
em que
m
J - fluxo de massa do soluto no solo [M/L
2
T];
z
q
- fluxo de água na direção vertical [L/ T]; e
C - concentração do soluto [M/L
3
].
O fluxo do transporte de soluto por difusão é dado por
z
C
DJ
dd
∂
∂
−=
)(
θ
(3.16)
em que
d
J - fluxo de difusão do soluto no solo [M/L
2
T];
d
D - coeficiente de difusão [L
2
/T];
C - concentração do soluto [M/L
3
];
θ
- umidade volumétrica [L
3
/L
3
]; e
z - distância [L].
O fluxo do transporte de soluto por dispersão mecânica é dado por
z
C
DJ
zzzz
∂
∂
−=
)(
θ
(3.17)
em que
76
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
zz
J
- fluxo de dispersão mecânica do soluto no solo [M/L
2
T];
zz
D
- coeficiente de dispersão mecânica [L
2
/T];
C - concentração do soluto [M/L
3
];
θ
- umidade volumétrica [L
3
/L
3
]; e
z
- distância [L].
Pela seção 2.4.1.3,
zzd
DDD
+
=
(3.18)
em que
D - coeficiente de dispersão hidrodinâmica [L
2
/T];
d
D - coeficiente de difusão [L
2
/T]; e
zz
D
- coeficiente de dispersão mecânica [L
2
/T].
A Lei da Conservação da Massa para transporte de solutos requer que a taxa de
variação de massa no volume de controle seja igual à taxa de soluto entrando mais a taxa de
produção desse soluto no volume de controle por vários processos físicos e químicos.
Considerando a taxa de soluto que entra no volume de controle, na direção z,
como
z
J
e a taxa de soluto que sai, na direção z, do volume de controle como )(
zz
J
z
J
∂
∂
+ ,
obtém-se, na direção z,
)()()(
2
2
C
z
DCq
z
J
z
zz
θ
∂
∂
+
∂
∂
−=
∂
∂
−
(3.19)
em que
zzdmz
JJJJ
+
+
=
(3.20)
z
J
- fluxo de soluto no solo [M/L
2
T];
m
J
- fluxo de massa do soluto no solo [M/L
2
T];
d
J
- fluxo de difusão do soluto no solo [M/L
2
T]; e
Universidade de São Paulo – Escola de Engenharia de São Carlos
77
zz
J
- fluxo de dispersão mecânica do soluto no solo [M/L
2
T].
A taxa da produção de soluto devido à sorção/desorção entre um soluto e o meio
poroso no volume de controle é dada por (ISTOK, 1989)
t
C
K
t
C
dbsorção
∂
∂
−=
∂
∂
ρ
θ
)
)(
(
(3.21)
em que
ρ
b
– densidade aparente do solo [ML
-3
]; e
K
d
– coeficiente de distribuição [L
3
M
-1
].
A taxa da produção de soluto devido à decaimento radioativo é dada por (ISTOK,
1989)
)()
)(
(
CKC
t
C
dbdecaimento
ρθλ
θ
+−=
∂
∂
(3.22)
em que
λ
– coeficiente de decaimento de primeira ordem [ T
-1
];
ρ
b
– densidade aparente do solo [ML
-3
]; e
K
d
– coeficiente de distribuição [L
3
M
-1
].
Considerando-se a taxa de variação de massa como
t
C
∂
∂
)(
θ
, a Lei da Conservação da
Massa e as equações (3.19), (3.21) e (3.22), obtém-se a Equação de Advecção-Dispersão dada
por
( ) ( )
(
)
( ) ( )
CKCCK
tz
C
D
z
Cv
z
C
t
bdbdzzz
ρθλρ
θ
θθ
+−
∂
∂
−
∂
∂
∂
∂
=
∂
∂
+
∂
∂
(3.23)
C
– concentração do soluto na solução do solo [ML
-3
];
v
z
– velocidade da água nos poros, na direção
z
[LT
-1
];
78
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
zz
D
– coeficiente de dispersão mecânica desde que a contribuição da difusão
molecular para a dispersão seja desprezada [L
2
T
-1
];
λ
– coeficiente de decaimento de primeira ordem [ T
-1
];
K
d
– coeficiente de distribuição [L
3
M
-1
];
ρ
b
– densidade aparente do solo [ML
-3
]; e
θ
– umidade volumétrica [ L
3
L
-3
].
Como o
fator retardo
R
devido à sorção do solo é dado por
R
= 1 +
θ
ρ
bd
K
(3.24)
obtém-se a Equação de Advecção- Dispersão considerada neste trabalho:
( ) ( )
(
)
CR-
θλ
θ
θθ
∂
∂
∂
∂
=
∂
∂
+
∂
∂
z
C
D
z
Cv
z
CR
t
zzz
(3.25)
Universidade de São Paulo – Escola de Engenharia de São Carlos
79
4 REVISÃO BIBLIOGRÁFICA
Neste capítulo, estão descritos estudos relacionados à simulação de fluxo e transporte
de solutos, que contribuíram para o desenvolvimento deste trabalho, fornecendo
fundamentação e alicerce, daí a importância deste. Na seção 4.1, é apresentado um breve
histórico sobre o fluxo de água e o transporte de solutos. Na seção 4.2, são mostrados estudos
sobre a modelagem do fluxo de água nos resíduos sólidos domiciliares ou na zona não-
saturada do solo. Na seção 4.3, é abordado o transporte de solutos em aterros sanitários ou em
solos não-saturados e na seção 4.4, o transporte de chorume em aterros sanitários ou na zona
não-saturada do solo.
4.1 Breve Histórico sobre o Fluxo de Água e o Transporte de Solutos
A primeira experiência que quantificou o fluxo num meio poroso saturado foi
publicada por Darcy (1856). A Lei de Darcy constitui a base dos métodos de avaliação
quantitativa de recursos hídricos subterrâneos (MANOEL FILHO, 2000).
A Lei de Darcy foi adaptada mais tarde para solos não-saturados, com a
denominação de Lei de Buckingham-Darcy ou Lei de Darcy-Buckingham (RICHARDS,
1931).
Richards formulou as bases teóricas para descrever a percolação da água em um
meio poroso não-saturado, em 1928, expondo o Princípio de Buckingham. A equação geral do
movimento da água em um meio poroso não-saturado, utilizando o Princípio da Conservação
da Massa foi apresentada em 1931 (RICHARDS, 1931).
Na década de 30, a teoria do fluxo em meios porosos foi muito usada na área de
engenharia de petróleo (MUSKAT, 1937). O transporte de solutos era aplicado
principalmente em estudos de intrusão de água marinha, em aqüíferos terrestres. Geralmente,
o método de análise se resumia num cálculo advectivo, no qual se supunha que o soluto se
80
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
movia com a velocidade média da água subterrânea e não recebia influências de processos de
adsorção, reação cinética, etc.
Na década de 50, o problema de transporte foi tema de pesquisa no campo da
engenharia química, principalmente transporte em meios porosos. Ficou reconhecido que o
transporte advectivo puro não descrevia completamente a dinâmica dos solutos. Foram
desenvolvidas teorias de dispersão hidrodinâmica (DE JOSSELIN DE JONG, 1958;
SAFFMAN, 1959; SCHEIDEGGER, 1963). Expressões descrevendo o processo de adsorção,
a partir de estudos em colunas de intercâmbio iônico, foram obtidas (VERMEULEN;
HIESTER, 1952). Nesta época, também, iniciou-se o estudo sobre o destino de materiais
nucleares no subsolo, considerando o regime do fluxo subterrâneo como um sistema
hidrodinâmico e físico-químico em conjunto, onde a qualidade da água depende de uma
seqüência de processos que ocorrem na trajetória do fluxo. Nas décadas seguintes, este
conceito foi aplicado ao estudo da evolução da qualidade da água subterrânea natural e à
problemas de poluição devido a materiais não-radioativos.
Também na década de 50, surgiram as primeiras aplicações da simulação, avaliando
a circulação das águas, usando analogias com modelos elétricos, em forma de circuitos, onde
o fluxo da água era representado por correntes elétricas. Definiram-se, portanto, as bases para
o início da simulação do transporte nas décadas futuras. Algumas considerações essenciais
foram identificadas, ou seja, (1) o desenvolvimento de uma teoria de transporte de solutos; (2)
a importância da simulação como uma ferramenta analítico-preditiva; e (3) o surgimento de
um ponto de vista hidroquímico onde a qualidade da água depende de uma seqüência de
reações que ocorrem na trajetória do fluido.
Durante os anos 60, ocorreram contribuições sobre a teoria de advecção-dispersão
(BEAR, 1961). A introdução dos avanços computacionais, na década seguinte, ajudou na
solução de problemas de campo cada vez mais complexos. O desenvolvimento do Método de
Elementos Finitos e a utilização de códigos de cálculo tridimensionais tiveram grande
influência neste processo.
Na década de 70, foi proposto o primeiro modelo numérico que solucionava a
Equação de Richards, usando a técnica “Line Sucessive Over-Relaxation” (LSOR), apesar das
instabilidades numéricas e dificuldades de convergência do modelo (FREEZE, 1971).
Após esta solução, uma variedade de técnicas, utilizando o Método de Diferenças
Finitas e o Método de Elementos Finitos, tem sido usada para resolver a Equação de Richards,
formulada tanto em termos de umidade volumétrica
θ
, quanto em termos de potencial
matricial
ψ
. Vários trabalhos foram apresentados por Neuman (1973), Narasimhan e
Universidade de São Paulo – Escola de Engenharia de São Carlos
81
Witherspoon (1976), Haverkamp
et al
. (1977), Hayhoe (1978), Haverkamp e Vauclin (1979),
Cooley (1983), Hornung e Messing (1983), Huyakorn, Lester e Mercer (1983), entre outros,
mas apresentavam significativos erros no balanço de massa.
Na década de 80, foi realizado o estudo da trajetória percorrida por uma substância
orgânica e foram determinados os correspondentes parâmetros de transporte e de
biodegradação (RUBIN, 1983). Modelos geoquímicos, considerando o equilíbrio entre a água
subterrânea e seu entorno, foram desenvolvidos (PARKHURST; THORSTENSON;
PLUMMER, 1990), bem como a introdução de reações de equilíbrio inorgânicas na
simulação do transporte de solutos (GILLHAM; CHERRY, 1982). Os trabalhos de Nielsen,
van Genuchten e Biggar (1986) e Milly (1988) contêm uma revisão geral da literatura
abordando estes assuntos.
Allen e Murphy (1986) e Celia, Ahuja e Pinder (1987) apresentaram a nova
formulação da Equação de Richards, a mista, que é matematicamente equivalente às
anteriores, ou seja, a formulada em termos de umidade volumétrica
θ
e em termos de
potencial matricial
ψ
. Allen e Murphy (1986) chamaram sua aproximação de Método Quasi-
Newton, enquanto Celia, Ahuja e Pinder (1987) referiram-se ao método como Método de
Picard Modificado, sendo que foi apresentado um excelente balanço de massa nas soluções
numéricas dos dois casos. Zarba (1988) usou o método iterativo de Picard Modificado, com
aproximações com diferenças finitas no tempo e elementos finitos no espaço e também
demonstraram um perfeito balanço de massa.
Na década de 90, Celia, Bouloutas e Zarba (1990) chegaram a conclusões
importantes, em que aproximações com elementos finitos podem produzir soluções
oscilatórias, apesar da conservação da massa, para problemas de infiltração, em solos
inicialmente muito secos.
Ross e Bristow (1990) linearizaram a Equação de Richards, usando a transformação
de Kirchhoff. A técnica apresentou bons resultados apenas para solos homogêneos, além de
aumentar os esforços computacionais. Srivastava e Yeh (1991) propuseram uma solução
analítica simplificada da Equação de Richards, sendo que as funções não-lineares presentes na
equação foram aproximadas por funções exponenciais.
Pode-se ressaltar a grande quantidade de códigos de cálculo de transporte, para o
caso de fluxo saturado (PADOIN et al., 2006).
Nos últimos dez anos, os temas de pesquisa concentram-se nas áreas de transporte de
solutos em meios fraturados (WENDLAND, 2004), de fluxo em sistemas multifásicos
(WENDLAND; FLENSBERG, 2002) e de transporte de solutos em meios não-saturados
82
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
(HANCOCK et al., 2008). Alguns dos trabalhos referentes a fluxo de água e transporte de
solutos na zona não-saturada do solo serão citados nas seções 4.2, 4.3 e 4.4. O crescimento da
potência dos computadores, realizando elevado número de operações em curto tempo, e o aumento de
técnicas de processamento permitirão a resolução de problemas de grande porte, representados por
modelos sofisticados.
4.2 Fluxo de Água nos Resíduos Sólidos Domiciliares ou na Zona Não-Saturada do Solo
Philip (1957) propôs uma solução semi-analítica da Equação de Richards no estudo
da infiltração de água em um meio poroso homogêneo, de profundidade infinita, com
umidade volumétrica uniforme
θ
i
e com a superfície do solo (condição de contorno) mantida
com umidade volumétrica
0
θ
>
θ
i
.
Straub e Lynch (1982) aplicaram a teoria do fluxo não-saturado em aterros sanitários
de resíduos sólidos domiciliares. Para modelar as características não-saturadas dos resíduos
sólidos domiciliares, usaram a equação (4.1) do potencial matricial
)(
θ
ψ
e a equação (4.2) da
condutividade hidráulica
K(θ)
dadas por
b
ss
−
=
]/[)(
θθψθψ
(4.1)
em que
s
ψ
- potencial matricial na saturação [L];
b
- parâmetro de ajuste [ ];
s
θ
- umidade volumétrica no ponto de saturação [L
3
L
-3
]; e
θ
– umidade volumétrica [L
3
L
-3
].
B
ss
KK
]/[)(
θθθ
=
(4.2)
em que
K(θ)
- condutividade hidráulica do solo não-aturado [LT
-1
];
B
- parâmetro de ajuste [ ]; e
Universidade de São Paulo – Escola de Engenharia de São Carlos
83
K
s
- condutividade hidráulica do solo saturado [LT
-1
].
Os valores de
s
ψ
= 100 cm,
b
= 7 e
B
= 8 ou
B
= 9 mostraram um bom ajuste aos
resultados experimentais realizados.
Korfiatis et al. (1984) formularam e calibraram o modelo matemático para simular o
movimento vertical unidimensional da água, em resíduos sólidos domiciliares, constituído
pela Equação de Richards e pelas equações (4.1) e (4.2), propostas por Straub e Lynch (1982).
Análises de sensibilidade mostraram que uma grande variação no parâmetro de ajuste
b
do potencial matricial teve pouco efeito na simulação, enquanto que uma pequena variação
no parâmetro de ajuste
B
da condutividade afetava significativamente os resultados da
simulação.
Foram as definições dadas à
θ
s
e
K
s
que constituíram a diferença básica entre os
estudos de Korfiatis et al. (1984) e de Straub e Lynch (1982). Korfiatis et al. definiram
θ
s
como a umidade volumétrica na saturação, ao passo que Straub e Lynch definiram
θ
s
igual à
capacidade de campo dos resíduos. Analogamente, Korfiatis et al. definiram
K
s
como a
condutividade hidráulica na saturação, enquanto que Straub e Lynch a definiram igual à taxa
de aplicação diária de água.
Demetracopoulos et al. (1986) desenvolveram análise de sensibilidade no modelo de
Korfiatis et al
.
(1984). O trabalho avaliou os resultados do modelo para solos saturados e não-
saturados. As simulações com superfícies não-saturadas foram mais sensíveis às mudanças na
condutividade hidráulica e no parâmetro de ajuste
B
. As variações do tamanho da malha e de
tempo na solução numérica pouco influenciaram nos resultados do trabalho.
Celia, Bouloutas e Zarba (1990) resolveram a Equação de Richards, formulada em
termos de potencial
ψ
, pelo Método de Elementos Finitos. Na aproximação espacial, os
autores consideram as funções-base lineares por partes e na derivada temporal, é aplicado o
esquema de Euler Implícito.
Noble e Arnold (1991) simularam o fluxo de água em um aterro sanitário, utilizando
a Equação de Richards, formulada em termos de umidade volumétrica, e em um experimento
em escala de laboratório, usando o modelo unidimensional FULFILL, elaborado por eles. O
modelo desenvolvido calcula, portanto, a variação da umidade e a variação do fluxo no
tempo, com condições iniciais e de contorno, através do Método de Diferenças Finitas,
incluindo o efeito da força gravitacional.
84
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Noble e Arnold (1991) estudaram as equações de
K(θ)
e
)(
θ
ψ
usadas por Straub e
Lynch (1982) e Korfiatis et al. (1984) e formularam uma relação exponencial para
K(θ)
e
)(
θ
ψ
(equações (4.3), (4.4)):
K (θ) = K
s
e
γ (θ*-1)
(4.3)
ψ
(θ)
=
s
ψ
e
-aθ*
(4.4)
Considerando
θ
*
= (θ – θ
r
) / (θ
s
– θ
r
)
(4.5)
K(θ)
- condutividade hidráulica do solo não-saturado [LT
1−
];
ψ
(θ)
– potencial matricial do solo não-saturado [L];
a
- parâmetro de ajuste [ ];
γ - parâmetro de ajuste [ ];
θ
*
- umidade volumétrica normalizada [L
3
L
-3
];
θ
s
- umidade volumétrica do solo no ponto de saturação [L
3
L
-3
]; e
θ
r
- umidade volumétrica residual do solo [L
3
L
-3
].
Uma diferença importante entre a formulação exponencial e a de potência para
K(θ)
e
)(
θ
ψ
é que a primeira assume um valor máximo finito para
ψ
quando a umidade é zero
(
θ
= 0), ao passo que a equação de potência prevê um
ψ
infinito, quando
θ
= 0.
Ahmed
et al.
(1992) desenvolveram o modelo
Flow Investigation for Landfill
(FILL)
aperfeiçoado mais tarde por Khanbivardi, Ahmed e Gleason
(1995). Na modelagem do
movimento de água, o modelo usa uma combinação de teorias de fluxo em meio saturado e
não-saturado. Para simular o fluxo de água no sistema de drenagem, usa-se a teoria do fluxo
em meio saturado e a teoria de fluxo em meio não-saturado para simular o fluxo de água
através dos resíduos sólidos domiciliares.
Al-Yousfi, Pohland e Vasuki (1992) formularam a equação da condutividade
hidráulica em função da umidade volumétrica (equações (4.6) e (4.7)). Para isto,
desenvolveram uma análise estatística de acordo com a teoria da entropia probabilística, que
usa o conceito de que o sistema tem uma tendência natural de se aproximar e de se manter em
Universidade de São Paulo – Escola de Engenharia de São Carlos
85
seu estado mais provável, juntamente com as técnicas de maximização e minimização (teoria
dos jogos) e de observações aleatórias (teoria da informação).
Para
θ > θ
r
( ) ( )
−
−
−
+−−=
sr
rs
explnKK
θ
θ
θθ
θθθ
1
1
1 (4.6)
Para
θ
≤
θ
r
(
)
0
=
θ
K
(4.7)
K(θ)
- condutividade hidráulica do solo não-saturado [LT
1
−
].
A equação (4.6) ajustou-se muito bem às equações de Noble e Arnold (1991) e às de
Korfiatis et al. (1984).
Zeiss e Major (1993) estudaram os padrões de fluxo de umidade nos resíduos sólidos
domiciliares e determinaram as variáveis que eram influenciadas pela variação da
compactação desses resíduos, usando colunas experimentais. Essas colunas foram preenchidas
com resíduos sólidos domiciliares e foram estudados, em função da compactação, a
densidade, a porosidade, a capacidade de campo, a condutividade hidráulica e o fluxo através
de caminhos preferenciais (
channeling
).
Os resultados do experimento mostraram que o fluxo por caminhos preferenciais não
era influenciado significativamente pelo aumento da densidade e conseqüente diminuição da
porosidade. Houve, também, mudanças muito pequenas no tempo de percolação da água, na
capacidade de campo e na condutividade hidráulica não-saturada, em função da variação da
compactação.
Zeiss e Uguccioni (1997) tentaram caracterizar o regime de fluxo nos canais ou
macroporos, determinando o número de Reynolds do escoamento (equação (4.8)), ou seja,
µ
ρ
qd
R
e
=
(4.8)
86
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
em que
R
e
- número de Reynolds [ ];
ρ
- densidade ou massa específica da água [M/L
3
];
q
- taxa específica de descarga [L/T];
d
- diâmetro médio do poro [L]; e
µ
- viscosidade [M/LT].
Em cinco dos oito contêineres estudados, o número de Reynolds encontrado foi
maior que dez, indicando que o fluxo estava acima do limite do fluxo laminar ou Darciano.
Bordier, Rathle e Zimmer (1997) realizaram experimentos que consistiram em
recircular água em materiais porosos grosseiros como, por exemplo, brita e geosintéticos, com
diferentes vazões e diferentes gradientes hidráulicos entre a entrada e a saída. Mostraram que
as equações de fluxo em materiais porosos grosseiros não apresentavam uma relação linear
entre o gradiente hidráulico e a velocidade macroscópica prevista por Darcy, até mesmo nos
experimentos com gradiente hidráulico baixo. O número de Reynolds para todos os
experimentos foi acima de 10.
Venkataramani e Rao (1998) investigaram os limites entre o fluxo laminar, o fluxo
turbulento transitório e o fluxo turbulento completo e mostraram que o fluxo turbulento
transitório foi o encontrado geralmente nas velocidades pesquisadas (10
-2
a 10
-4
m/s).
Prasad, Kumar e Sekhar (2001) simularam o fluxo unidimensional de água em zona
não-saturada. A solução da Equação de Richards, formulada em termos de potencial
ψ
, foi
obtida pelo Método de Elementos Finitos. Na aproximação espacial, são consideradas as
funções-base lineares lagrangeanas (FARTHING et al., 2007) e, na derivada temporal, é
aplicado o esquema de Euler Implícito. Os resultados obtidos foram validados com os dados
de Celia, Bouloutas e Zarba (1990). Foram mostrados os efeitos da variação dos parâmetros
do solo não-saturado
α
e
n
(equações (5.4) e (5.5), seção 5.1) em relação ao potencial
matricial e a condutividade hidráulica.
Vasconcellos e Amorim (2001) realizaram a simulação numérica da infiltração de
água em meios porosos não-saturados homogêneos. Utilizaram o Método de Diferenças
Central no espaço, o esquema de Euler Explícito no tempo, o Método Iterativo de Picard para
resolver o sistema de equações algébricas e as médias aritmética, harmônica e geométrica
para o cálculo da condutividade hidráulica.
Universidade de São Paulo – Escola de Engenharia de São Carlos
87
Prevedello et al
.
(2002) desenvolveram um modelo numérico na linguagem
Beginner’s All-purpose Symbolic Instructional Code
(BASIC), que simula o processo de
infiltração vertical da água no solo homogêneo ou estratificado, com uma dada umidade
volumétrica inicial. Este modelo foi desenvolvido linearizando a Equação de Richards no
espaço e usando o Método de Newton-Raphson para resolver a equação no tempo.
Miller, Abhishek e Farthing
(2006) resolveram a Equação de Richards, em uma
dimensão, usando método adaptativo, tanto no tempo quanto no espaço. A solução dessa
equação é obtida pelo Método de Diferenças Finitas. A discretização espacial h-adaptativa é
implementada no Método de Linhas Estruturadas. A adaptatividade temporal usa ordem
variável, a aproximação com passo de tempo variável é baseada no esquema de Euler
Implícito.
Caputo e Stepanyants (2008) estudaram o fluxo de água para três modelos de
retenção de água no solo, o de Brooks-Coley, Mualen-van Genuchten e Storm-Fujita. Os
resultados foram comparados. Todos os modelos mostraram estabilidade no movimento das
frentes de molhamento. A Equação de Richards foi integrada numericamente pelo Método de
Volumes Finitos.
4.3 Transporte de Solutos em Solos Não-Saturados
Multimedia Exposure Assessment Model
(MULTIMED) foi desenvolvido pela U.S.
Environmental Protection Agency (EPA) para simular o movimento de contaminantes no
aterro (SALHOTRA et al., 1990; SHARP-HANSEN et al., 1990). O modelo consiste de
módulos que estimam liberação de contaminantes de um aterro para o ar, o solo, o lençol
freático e a água da superfície.
Ritterling e Stansbury (1998) usaram o MULTIMED, que utiliza métodos de solução
analítica e semi-analítica, para resolver equações matemáticas que descrevem o fluxo e o
transporte de contaminantes.
Klein et al. (1997) executaram modelo de simulação para predizer o fluxo de água, a
concentração de pesticida no solo e a concentração de pesticida no chorume avaliado em 39
lisímetros. Foram utilizados três modelos computacionais diferentes,
Pesticide Root Zone
Model - 1
(PRZM-1) (CARSEL et al., 1984),
Pesticide Leaching Model 1.5
(PELMO 1.5)
(KLEIN, 1994) e
Pesticide Leaching Model 2.0
(PELMO 2.0) (KLEIN, 1995). Para calcular a
88
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
evapotranspiração, o modelo PRZM-1 usa a Equação de Hamon (HAMON, 1961), enquanto
que o PELMO 1.5 e o PELMO 2.0 utilizam a Equação de Haude (SEILER; GAT, 2007). Para
calcular a sorção, usam a Equação de Freundlich. PELMO 2.0 prediz o chorume acumulado
do lisímetro com boa aproximação, enquanto que PRZM-1 e PELMO 1.5 são menos precisos.
É apropriado para medir volumes de chorume e concentrações de pesticidas no chorume.
PELMO inclui advecção e dispersão.
O modelo matemático
Simulação de Movimento de Água e Solutos no Solo
(SIMASS) (COSTA et al
.,
1999) simulou o transporte unidimensional de água e soluto no
solo, sob condições de escoamento não-permanente. As equações são resolvidas
numericamente pelo Método de Diferenças Finitas. O modelo permite, entre outras
características, obter as distribuições de umidade volumétrica e de concentração de nitrato no
solo, utilizando-se condições de contorno do tipo potencial constante e do tipo fluxo constante
e, ainda, estimar a condutividade hidráulica do solo não-saturado. O desempenho do SIMASS
foi comparado com o modelo CXTFIT, versão 1.0 (COSTA, 1998; ROSSI; MIRANDA;
DUARTE, 2007).
Tsai et al. (2000) empregaram o Método Analítico Finito para resolver as equações
de fluxo e transporte de soluto (BEAR, 1979), em duas dimensões, na zona não-saturada. As
soluções numéricas analíticas finitas obtidas foram verificadas com as soluções de elementos
finitos do código
Finite Element Model of Water
(FEMWATER) (YEH; WARD, 1980; YEH,
1987).
O modelo computacional
Miscible Displacement
(MIDI) (MIRANDA, 2001) foi
desenvolvido em linguagem de programação Visual Basic 5.0 para simulação da dinâmica de
solutos no solo. As equações foram resolvidas numericamente em um sistema de volumes
finitos. O modelo apresentou bom ajuste para os perfis de concentração de nitrato e de
umidade simulados em relação aos medidos, em condições de laboratório, em coluna vertical
de solo não-saturado.
Carey, Bidwell e Mclaren (2002) mostraram que alta concentração de cromo foi
encontrada no chorume dos lisímetros que receberam a aplicação de CCA, ou seja, soluções
de cobre, cromo e arsênio, na Nova Zelândia. Foi feita uma simulação da chuva sobre os
lisímetros, por um período de 102 dias, depois da aplicação de CCA. Em média, 26% do
cromo aplicado foram coletados no chorume após os 102 dias e 74% do cromo ficaram retidos
no perfil do solo após a lixiviação. Nem o cobre, nem o arsênio foram detectados no chorume,
indicando que estes elementos ficaram retidos no perfil do solo. Foi utilizado o modelo
State–
Space Mixing Cell
(SSMC) (BIDWEL, 1999) para simular os processos do transporte
Universidade de São Paulo – Escola de Engenharia de São Carlos
89
advectivo-dispersivo e as curvas breakthrough (BTCs), operando dentro de
Matrix
Laboratory
(MATLAB) desenvolvido por Moler (1970).
Abdou e Flury (2004) fizeram uma simulação numérica bi-dimensional num solo
arenoso, em condições de fluxo de água não-saturado. Usou-se o Método de Elementos
Finitos e o código CHAIN-2D (SIMUNEK; van GENUCHTEN, 1994). As funções )(
ψ
θ
e
K
)(
ψ
, consideradas neste trabalho, foram dadas por van Genuchten (1980). Foram feitas
comparações de fluxo e transporte com três diferentes estruturas de solo (isotrópico,
horizontal e vertical). Os resultados mostraram que a formação de empoçamento ocorre na
base do lisímetro para as três estruturas do solo e que ele ocorre mais rapidamente e foi mais
pronunciado com a estrutura vertical (efeito do fluxo preferencial). Curvas breakthrough de
um soluto conservativo (brometo) mostraram que solutos se movem mais rapidamente em
campo do que no lisímetro. Menos diferenças entre lisímetros e campo foram encontradas
com a estrutura horizontal do que com as estruturas isotrópica e vertical.
Buczko et al. (2004) simularam o transporte de cromo na zona não-saturada para
predizer a chegada do contaminante no lençol freático. O trabalho consistiu de três partes:
a) estimativa das características emitidas dos metais pesados na zona contaminada
perto da superfície do solo;
b) simulações de transporte baseadas em experimentos com lisímetro monolítico; e
c) simulação a longo-prazo da emissão do contaminante, transporte e inibição para o
lençol freático.
Para simular o fluxo da água e o transporte do cromo, foi usado o modelo de
elementos finitos HYDRUS-1D (SIMUNEK; SEJNA; VAN GENUCHTEN
,
1998).
Miranda et al. (2005), utilizando o modelo MIDI (MIRANDA, 2001), apresentaram a
simulação do deslocamento de potássio em colunas verticais de solo não-saturado. Os
parâmetros de transporte do íon potássio foram determinados no Latossolo Vermelho-
Amarelo, fase arenosa, através do modelo CXTFIT (TORIDE; LEIJ; van GENUCHTEN,
1999), desenvolvido pelo U. S. Salinity Laboratory-USDA-Riverside-CA, versão 2.1.
Concluiu-se que o modelo foi capaz de simular, com bom ajuste, a umidade
volumétrica e o deslocamento do potássio, mostrando inclusive o seu retardamento, em
relação à frente de molhamento.
Costa e Castro (2007) propuseram um método numérico-analítico, ou seja, numérico
em relação ao espaço e analítico em relação ao tempo, para problemas transientes de
90
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
transporte de contaminantes em uma dimensão. A equação unidimensional de transporte de
contaminante usada neste trabalho foi dada por Istok (1989).
Rocha et al. (2008) simularam as concentrações de nitrato e de amônio no solo,
considerando-se as transformações biológicas e o efeito da temperatura e do teor de umidade
volumétrica do solo, através de modificações no modelo de transporte de soluto no solo
denominado
Simulação de Movimento de Água e Solutos no Solo, considerando a presença
de Cultura
(SIMASS-C) (COSTA, 1998), aprimorado por Corrêa (2001). A elaboração da
rotina computacional, na linguagem Delphi 7.0, resultou no SIMASS-C modificado (ROCHA
et al., 2008). Foram considerados os processos de mineralização e nitrificação no modelo,
resultando numa melhoria na estimativa da concentração de nitrato e de amônio.
4.3.1 Transporte de Chorume de Aterros Sanitários
Lu (1996) apresentou um modelo matemático para prever a quantidade e qualidade
de chorume do incinerador de resíduo do aterro. O modelo foi baseado na equação de fluxo na
zona não-saturada (STRAUB; LYNCH, 1982), na equação de transporte de soluto (STRAUB;
LYNCH, 1982) e na equação de balanço hídrico do sistema de coleta do chorume
(DEMETRACOULOS, 1988). Usando um esquema de Diferenças Finitas Explícitas no
aterro, as equações diferenciais parciais foram discretizadas na direção vertical e
transformadas nas equações diferenciais ordinárias, que foram resolvidas usando o integrador
Adams-Moulton do IMSL (International Mathematical and Statistical Libraries) de 1987.
Chen e Wang (1997) estudaram o fenômeno de transporte de chorume no aterro
sanitário de Taichung (Taiwan). Para identificar os parâmetros e simular o modelo
matemático do transporte das matérias orgânicas do chorume foi utilizado o SEFTRAN
(NAGRA, 1994). O modelo é de transporte tri-dimensional e foi simplificado para uma
dimensão, usa o Método de Elementos Finitos e o Método de Galerkin (ISTOK, 1989) e seu
mecanismo inclui convecção, dispersão, retardamento e decaimento biológico.
Kim
et al.
(1999) utilizaram o Método de Diferenças Finitas para reduzir a
contaminação das águas subterrâneas próximas ao Aterro de Kimpo na Coréia. O potencial
matricial e a concentração de poluentes foram modelados usando-se
Flow Model
(MODFLOW) (McDONALD; HARBAUGH, 1988) e
Modular Mass transport 3-Dimensions
Universidade de São Paulo – Escola de Engenharia de São Carlos
91
(MT3D) (ZHENG, 1990). O primeiro é um modelo de fluxo tridimensional e MT3D, um
modelo de transporte.
McCreanor e Reinhart (2000) usaram o modelo
The United States Geological
Survey’s Saturated-Unsaturated Flow and Transport
(SUTRA) (VOSS, 1984) para modelar
chorume em massas de resíduo homogêneas anisotrópicas e heterogêneas. SUTRA é bi-
dimensional e utiliza os Métodos de Diferenças Finitas e de Elementos Finitos.
Fatta, Naoum e Loizidou (2002) estudaram a caracterização do chorume originário
do aterro Ano Liosia, em Atenas, Grécia, e a qualidade do aqüífero local, com os mesmos
modelos usados por Kim et al. (1999). O problema foi resolvido pelo Método de Diferenças
Finitas.
Bou-Zeid e El-Fadel (2004) usaram o modelo
Hydrologic Evaluation of Landfill
Performance
(HELP) (SCHROEDER et al., 1994) para estimar a quantidade de chorume e a
percolação na subsuperfície. Usaram também o modelo tri-dimensional PORFLOW
(RUNCHAL; SAGAR, 1998) para simular o fluxo unidimensional de água no solo e o
transporte de contaminante, em sistema de coordenadas cartesianas e cilíndricas. Foi feita
uma análise sensitiva com os parâmetros que controlam o transporte do chorume.
Haydar e khire (2005) fizeram um estudo numérico do sistema de recirculação de
chorume (LRS) consistindo de valas horizontais. Utilizou o modelo com elementos finitos
HYDRUS-2D (SIMUNEK; SEJNA; VAN GENUCHTEN, 1999), que simula a água e o
transporte de solutos na zona não-saturada.
Bunsri, Sivakumar e Hagare (2008 a) mostraram que o traçador cloreto de sódio pode
ser usado para avaliar o movimento de água em solo não-saturado. O transporte do traçador
em zona não-saturada do solo é dado pela Equação de Advecção- Dispersão. Os parâmetros
não-lineares da Equação de Richards foram obtidos pelas Equações de van Genuchten (1980).
O modelo foi resolvido usando Elementos Finitos de Galerkin.
92
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Universidade de São Paulo – Escola de Engenharia de São Carlos
93
5 METODOLOGIA
Neste capítulo, é apresentado o modelo matemático que descreve a simulação do
fluxo, o transporte de solutos e os principais métodos aplicados na solução das equações desse
modelo.
O modelo matemático para descrever o fluxo e o transporte de solutos, na zona não-
saturada do solo, consiste na solução de equações diferenciais parciais de segunda ordem, isto
é, a Equação Diferencial Parcial (equação 5.1) conhecida como Equação de Richards
(RICHARDS, 1931)
,
que rege o movimento da água no solo, com a Equação de Advecção-
Dispersão (BEAR, 1979; HINDMARSH; GRESHO; GRIFFITHS
,
1984) (equação 5.22), do
transporte de soluto (MUALEM, 1976; GENUCHTEN, 1980; SINGH; KANWAR, 1995;
OLIVEIRA, 1999), que busca a evolução temporal da concentração de solutos na zona não-
saturada do solo acompanhadas das condições iniciais e de contorno (RICHTMYER;
MORTON, 1967; SOD, 1987).
Na seção 5.1, são apresentados a Equação de Richards, com os parâmetros físico-
hídricos do solo, descritos pelas Equações de van Genuchten (1980), as condições iniciais e
de contorno, o princípio do Método dos Resíduos Ponderados e a aproximação de elementos
finitos, enquanto que na seção 5.2, a Equação de Advecção-Dispersão, com apresentação de
suas parcelas, as condições iniciais e de contorno, o princípio do Método dos Resíduos
Ponderados e a aproximação de elementos finitos. Na seção 5.3, é apresentada a
implementação computacional das Equações de Richards e de Advecção-Dispersão.
5.1 Movimento da Água no Solo: Equação de Richards
O processo do movimento de água no solo é descrito pela Equação de Richards, uma
equação de advecção-difusão, não-linear, que pode ser escrita como uma lei de conservação
para o conteúdo de água, a quantidade de água contida num dado volume de solo. O termo
94
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
advectivo é devido à gravidade, enquanto o termo difusivo vem da Lei de Darcy. A Equação
de Richards, estimando os valores do potencial matricial de água, e considerando o
escoamento unidimensional não-saturado e a coordenada vertical orientada positivamente
para cima, é dada por
( ) ( )
+
∂
∂
∂
∂
=
∂
∂
1
z
K
zt
C
s
ψ
ψ
ψ
ψ
(5.1)
em que
−=
ψ
θ
ψ
d
d
C
s
)(
solo do específica hídrica capacidade
[L
1
−
];
−
θ
umidade volumétrica [L
3
L
-3
];
−
ψ
[
]
L matricial potencial
;
−
)(
ψ
K
saturado-não solo do hidráulica adecondutivid
[LT
1
−
];
z
– coordenada vertical [L]; e
t
– tempo [T].
As condições iniciais e de contorno utilizadas na resolução da equação 5.1 são:
•
)()0,(
zz
inicial
ψψ
= Lz
≤
≤
0,
•
0
),0(
ψ
ψ
=
t
0,
>
t
•
L
tL
ψ
ψ
=
),( 0,
>
t
Para resolver a Equação de Richards, algumas relações constitutivas precisam ser
especificadas.
Considere
S
e
=
[
]
m
n
)1/(1
αψ
+
(5.2)
com
m
= 1 -
n
1
e
S
e
, a saturação efetiva, dada por
Universidade de São Paulo – Escola de Engenharia de São Carlos
95
S
e
=
rs
r
θθ
θ
θ
−
−
(5.3)
em que
θ
s
– umidade volumétrica do solo saturado [L
3
L
–3
];
θ
r
– umidade volumétrica residual do solo [L
3
L
–3
];
α
– parâmetro que depende do solo [L
-1
]; e
m
e
n
– parâmetros que dependem do solo, [ ].
Das equações (5.2) e (5.3), obtém-se a relação entre
θ
e
ψ
dada pela Equação de van
Genuchten (1980):
n
m
e
S
1
1
1
1
α
ψ
−
=
(5.4)
A relação entre
K
e
θ
é dada pela Equação de Mualen (1976) (REICHARDT, 1996):
( )
2
1
2
1
11
−
−
−−
−
−
=
m
m
rs
r
rs
r
s
KK
θθ
θθ
θθ
θθ
θ
(5.5)
em que
K(θ)
- condutividade hidráulica do solo não-saturado [LT
1
−
]; e
s
K
- a condutividade hidráulica saturada [LT
–1
].
A capacidade hídrica específica do solo é dada por
( )
(
)
( )
[ ]
1
1
.1
...
+
−
+
−
=
m
n
n
sr
n
s
nm
C
ψα
ψθθα
ψ
(5.6)
Para a solução numérica da Equação de Richards (equação 5.1), é aplicado o Método
de Elementos Finitos, cuja fundamentação matemática é descrita na seqüência.
96
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
5.1.1 Método dos Resíduos Ponderados
Sejam o domínio
(
)
L,0
=
Ω
, o espaço funcional
(
)
(
)
(
)
{
}
1,;
221
≤Ω∈∂Ω∈=Ω
αϑϑ
α
LLH
(5.7)
e o sub-espaço de funções teste dado por
(
)
(
)
(
)
{
}
0,00;
1
==Ω∈= LHV
ϑϑϑ
(5.8)
em que
(
)
Ω
2
L
- espaço de funções de quadrado integrável.
O princípio do Método dos Resíduos Ponderados consiste em minimizar o resíduo no
domínio de estudo. Para obtenção da solução da Equação de Richards, basta multiplicar a
equação (5.1) por uma função teste
V
∈
ν
e integrá-la sobre o domínio
Ω
:
( ) ( ) ( )
0
00
=
+
∂
∂
∂
∂
−
∂
∂
∫∫
LL
s
dzK
z
K
z
dz
t
C
νψ
ψ
ψν
ψ
ψ
(5.9)
ou
( ) ( )
(
)
0
000
=
∂
∂
−
∂
∂
∂
∂
−
∂
∂
∫∫∫
LLL
s
dz
z
K
dz
z
K
z
dz
t
C
ν
ψ
ν
ψ
ψν
ψ
ψ
(5.10)
Assumindo que
(
)
ψ
K é suave e aplicando o Método de Integração por partes em
( )
ν
ψ
ψ
∂
∂
∂
∂
∫
z
K
z
L
0
dz
obtém-se
Universidade de São Paulo – Escola de Engenharia de São Carlos
97
( ) ( )
(
)
=
∂
∂
−
∂
∂
+
∂
∂
∫∫∫
LLL
s
dz
z
K
dz
dz
d
z
Kdz
t
C
000
ν
ψ
νψ
ψν
ψ
ψ
( ) ( )
Lzz
z
K
z
K
==
∂
∂
+
∂
∂
−=
ν
ψ
ψν
ψ
ψ
0
(5.11)
Como as funções teste
ν
pertencem ao espaço V,
(
)
(
)
00
=
=
L
ν
ν
,
( ) ( )
(
)
0
000
=
∂
∂
−
∂
∂
+
∂
∂
∫∫∫
dz
z
K
dz
dz
d
z
Kdz
t
C
LL
s
L
ν
ψ
νψ
ψν
ψ
ψ
(5.12)
A derivada temporal
t
∂
∂
ψ
é aproximada por um quociente de diferença finita
(SPERANDIO; MENDES; SILVA, 2006):
t
t
nn
∆
ψψψ
−
≅
∂
∂
+ 1
(5.13)
e, aplicando-se Euler Explícito (BEAR, 2007) (seção2.10.1.1), obtém-se
( ) ( )
(
)
0
00
1
0
=
∂
∂
−
∂
∂
+
∆
−
∫∫∫
+
dz
z
K
dz
dz
d
z
Kdz
t
C
n
L
n
n
L
nn
n
s
L
ν
ψ
νψ
ψν
ψψ
ψ
(5.14)
Assim, resolver a equação (5.1) consiste em encontrar
(
)
Ω∈
+ 11
),( Htz
n
ψ
que atenda
as condições iniciais
)()0,(
0
zz
inicial
ψψ
=
e de contorno
0
),0(
ψ
ψ
=
t
e
L
tL
ψ
ψ
=
),(
, e que
satisfaça:
(
)
=
∆
+
∫
dzC
t
nn
s
L
νψψ
1
0
1
( ) ( )
(
)
∫∫∫
∂
∂
+
∂
∂
−
∆
=
L
n
L
n
n
L
nn
s
dz
z
K
dz
dz
d
z
KdzC
t
000
1
ν
ψ
νψ
ψνψψ
(5.15)
para qualquer
Vv
∈
. A única incógnita do problema é
1+n
ψ
que deve ser calculada
iterativamente até que se obtenha o tempo total de simulação desejado.
98
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
5.1.2
Aproximação de Elementos Finitos
É necessário especificar a forma da solução aproximada e da função teste. O Método
de Galerkin é o mais usado em fluxo de água e transporte de solutos. A aproximação de
elementos finitos ou aproximação de Galerkin (BUNSRI; SIVAKUMAR; HAGARE, 2008 a)
consiste em aproximar as funções variável de estado
(
)
Ω∈
1
H
ψ
e teste
Vv
∈
por funções
aproximadas
(
)
Ω⊂Π∈
11
H
ψ
e
Vv ⊂Π∈
1
0
, com
1
Π
um subespaço finito de
1
H
e
(
)
(
)
(
)
{
}
0,00;
11
0
==ΩΠ∈=Π L
ϑϑϑ
. Neste trabalho, o subespaço
1
Π
será construído
por funções polinomiais por partes com suporte compacto, com grau p, que é a ordem
polinomial de aproximação. Assim, tem-se que:
∑
=
=
nf
j
jj
1
ϕαψ
e
∑
=
=
nf
i
i
v
1
ϕ
. Substituindo-se em
5.15, obtém-se
( )
=
∆
∫
∑∑
==
+
L
nf
i
i
nf
j
j
n
j
n
s
dzC
t
0
11
1
1
ϕϕα
∫
∑
∫
∑
∫
∑
===
∂
∂
+
∂
∂
−
∆
=
L
nf
i
i
n
L
nf
i
i
n
n
L
nf
i
i
n
n
s
dz
z
K
dz
dz
d
z
KdzC
t
0
1
0
1
0
1
1
ϕ
ϕ
ψ
ϕψ
(5.16)
em que
( )
∑
=
=
nf
j
j
n
j
n
1
ϕαψ
,
∑
=
=
∂
∂
nf
j
j
n
j
n
dz
d
z
1
ϕ
α
ψ
,
(
)
n
s
n
s
CC
ψ
=
e
(
)
nn
KK
ψ
=
(5.17)
Agrupando-se os termos:
( )
=
∆
∑∑
∫
= =
+
nf
i
nf
j
L
ij
n
j
n
s
dzC
t
1 1
0
1
1
ϕϕα
∑
∫∫∫
=
∂
∂
+
∂
∂
−
∆
=
nf
i
L
i
n
L
i
n
n
L
i
n
n
s
dz
z
K
dz
dz
d
z
KdzC
t
1
000
1
ϕ
ϕ
ψ
ϕψ
(5.18)
Universidade de São Paulo – Escola de Engenharia de São Carlos
99
A equação (5.18) deve ser satisfeita para qualquer
1
0
Π∈v
. Logo, tomando-se uma
função
i
ϕ
por vez, pode-se escrevê-la em forma matricial:
[
]
{
}
{
}
FK
n
=
+1
α
(5.19)
em que
∫
∆
=
L
ij
n
sij
dzC
t
K
0
1
ϕϕ
(5.20)
∫∫∫
∂
∂
+
∂
∂
−
∆
=
L
i
n
L
i
n
n
L
i
n
n
si
dz
z
K
dz
dz
d
z
KdzC
t
F
000
1
ϕ
ϕ
ψ
ϕψ
(5.21)
Assim, o Método de Elementos Finitos é aplicado para encontrar os coeficientes
multiplicadores
{
}
1
+
n
α
que satisfazem o problema algébrico definido pela equação (5.19).
5.2 Transporte de Solutos: Equação de Advecção-Dispersão
Para simular o transporte de solutos, em uma dimensão, na zona não-saturada do
solo, será utilizada a Equação Linear de Advecção-Dispersão (BEAR, 1979):
( ) ( )
(
)
sCR- +
∂
∂
∂
∂
=
∂
∂
+
∂
∂
θλ
θ
θθ
z
C
D
z
Cv
z
CR
t
zzz
(5.22)
Como o fator retardo R é dado por
R = 1 +
θ
ρ
bd
K
(5.23)
tem-se
( ) ( )
(
)
( ) ( )
s ++−
∂
∂
−
∂
∂
∂
∂
=
∂
∂
+
∂
∂
CKCCK
tz
C
D
z
Cv
z
C
t
bdbdzzz
ρθλρ
θ
θθ
(5.24)
100
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
em que
C – concentração do soluto na solução do solo [ML
-3
];
R – fator retardo devido à sorção do solo [ ];
v
z
– velocidade da água no poro, na direção z [LT
-1
];
zz
D
– coeficiente de dispersão mecânica desde que a contribuição da difusão
molecular para a dispersão seja desprezada [L
2
T
-1
];
λ – coeficiente de decaimento de primeira ordem [ T
-1
];
s – fonte ou sumidouro;
K
d
– coeficiente de distribuição [L
3
M
-1
];
ρ
b
– densidade aparente do solo [ML
-3
]; e
θ – umidade volumétrica [ L
3
L
-3
].
A velocidade da água nos poros é obtida a partir da Equação de Darcy-Buckingham
(REICHARDT; TIMM, 2004)
q =
)()( zgradK
+
−
ψ
θ
(5.25)
Portanto, o fluxo de água em solo não-saturado, na direção z, é dado por
q
z
=
z
KK
∂
∂
−−
ψ
θθ
)()(
(5.26)
Como a velocidade da água nos poros em solo não-saturado (BEAR, 1979) é dada
por
v=
θ
q
(5.27)
conclui-se que a velocidade da água nos poros em solo não-saturado, na direção z, é dada por
v
z
= -K (θ) ( 1+
∂
∂
z
ψ
)
θ
1
(5.28)
O coeficiente de dispersão mecânica pode ser expresso por (ANDERSON, 1979)
D
zz
= α
l
v
z
(5.29)
Universidade de São Paulo – Escola de Engenharia de São Carlos
101
em que
α
l
– coeficiente de dispersividade longitudinal [L].
O coeficiente de decaimento de primeira ordem pode ser relacionado com o tempo de
meia vida,
2/1
T
, (BEAR, 2007) dado por
2/1
2ln
T
=
λ
(5.30)
Na resolução da equação (5.22), serão utilizadas as condições iniciais e de contorno:
•
)()0,( zCzC
inicial
=
,
Lz
≤
≤
0
•
0
),0( CtC
=
,
0
>
t
•
L
CtLC
=
),( , 0
>
t
Observe que K
d
, ρ
b
, α
l
e
2/1
T são dados de entrada, enquanto que
ψ
,
z
∂
∂
ψ
, θ e K(θ)
são obtidos da Equação de Richards (equação (5.1)). Definido o modelo matemático, o
interesse situa-se na solução, ou seja, na resolução simultânea das Equações de Richards e de
Advecção-Dispersão. A solução numérica da Equação de Advecção-Dispersão (equação
(5.22)) é obtida através do Método de Elementos Finitos, da mesma forma que a Equação de
Richards, cuja fundamentação matemática é descrita na seqüência.
5.2.1 Método dos Resíduos Ponderados
São utilizados os mesmos espaços de aproximação utilizados no problema de
Richards. Considerem-se o domínio
(
)
L,0
=
Ω
, o espaço funcional
(
)
(
)
(
)
{
}
1,;
221
≤Ω∈∂Ω∈=Ω
αϑϑ
α
LLH
(5.31)
e o sub-espaço de funções teste dado por
102
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
(
)
(
)
(
)
{
}
0,00;
1
==Ω∈= LHV
ϑϑϑ
(5.32)
em que
(
)
Ω
2
L
- espaço de funções de quadrado integrável.
Como o princípio do Método dos Resíduos Ponderados consiste em minimizar o
resíduo no domínio de estudo, para a obtenção da solução, basta multiplicar a equação (5.22)
por uma função teste V
∈
ν
e integrá-la sobre o domínio
Ω
.
(
)
(
)
(
)
−
∂
∂
∂
∂
=
∂
∂
+
∂
∂
∫∫∫
zd
z
C
D
z
zd
z
C
zd
t
CR
zz
LL
z
L
ν
θ
ν
θν
ν
θ
000
∫∫
+Ω−
LL
zdsdCR
00
ννθλ
(5.33)
Fazendo-se integrações por partes, obtém-se
(
)
[ ]
(
)
−
∂
∂
=−+
∂
∂
=
=
=
=
∫∫
Lz
z
zz
L
z
Lz
z
z
L
z
C
Dzd
dz
d
CCzd
t
CR
0
0
0
0
ν
θ
ν
θννθνν
θ
(
)
zdszdRzd
dz
d
z
C
D
LLL
zz
ννθλ
ν
θ
∫∫∫
+−
∂
∂
−
000
(5.34)
Como as funções teste
ν
pertencem ao espaço V,
(
)
(
)
,00
=
=
L
ν
ν
(
)
(
)
−
∂
∂
−=−
∂
∂
∫∫∫
zd
dz
d
z
C
Dzd
dz
d
Czd
t
CR
L
zz
L
z
L
000
ν
θ
ν
θνν
θ
zdszdCR
LL
ννθλ
∫∫
+−
00
(5.35)
A derivada temporal será aproximada por um esquema de Euler Implícito (BEAR,
2007) (seção 2.10.1.1):
Universidade de São Paulo – Escola de Engenharia de São Carlos
103
(
)
zd
dz
d
z
C
Dzd
dz
d
Czd
t
CRCR
L
n
zz
L
n
z
L
nn
∫∫∫
∂
∂
−=−
∆
−
+
+
+
0
1
0
1
0
1
ν
θ
ν
θνν
θθ
-
zdszdCR
LL
n
ννθλ
∫∫
+−
+
00
1
(5.36)
e a solução consiste em encontrar
(
)
Ω∈
+
11
HC
n
θ
que atenda as condições iniciais e de
contorno
)()0,(
0
zCzC
inicial
θθ
=
,
0
),0( CtC
θ
θ
=
e
L
CtLC
θ
θ
=
),(
, e que satisfaça:
(
)
−
∂
∂
−+
∆
=
∆
∫∫∫∫
+
+
+
zd
dz
d
z
C
Dzd
dz
d
Czd
t
CR
zd
t
CR
L
n
zz
L
n
z
L
n
L
n
0
1
0
1
00
1
ν
θ
ν
θνν
θ
ν
θ
zdszdCR
LL
n
ννθλ
∫∫
+−
+
00
1
(5.37)
para qualquer
Vv
∈
. Ressalta-se que
z
∂
∂
ψ
e demais variáveis que dependem de
ψ
são obtidas
da solução da Equação de Richards.
A variável de estado da formulação integral apresentada é o produto
C
θ
entre a
umidade volumétrica e a concentração. Para se obter a concentração do soluto
C
, divide-se
C
θ
pela umidade volumétrica
θ
calculada na solução do problema de Richards. Optou-se por
utilizar
C
θ
como variável de estado em virtude de sua solução ser bastante suave, enquanto
que
θ
e C apresentam uma frente de propagação que pode conter fortes gradientes, o que
pode induzir oscilações na solução de elementos finitos.
As condições iniciais e de contorno de
C
θ
devem ser obtidas pelo produto de
θ
e
C
. O valor de
θ
é obtido da simulação de Richards e
C
é dado de entrada do problema.
5.2.2 Aproximação de Elementos Finitos
A aproximação de elementos finitos ou aproximação de Galerkin consiste em
aproximar as funções variável de estado
(
)
Ω∈
1
HC
θ
e teste Vv
∈
por funções aproximadas
(
)
Ω⊂Π∈
11
HC
θ
e
Vv ⊂Π∈
1
0
, com
1
Π
um subespaço finito de
1
H
e
104
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
(
)
(
)
(
)
{
}
0,00;
11
0
==ΩΠ∈=Π L
ϑϑϑ
. Neste trabalho, o subespaço
1
Π
será construído
por funções polinomiais por partes com suporte compacto, com grau p, que é ordem
polinomial de aproximação. Assim, tem-se que:
∑
=
=
nf
j
jj
C
1
ϕαθ
e
∑
=
=
nf
i
i
v
1
ϕ
. Substituindo-se
em 5.37:
( ) ( ) ( )
−+
∆
=
∆
∫
∑∑
∫
∑∑∑
∫
∑
==
+
====
+
zd
dz
d
zd
t
R
zd
t
R
L
nf
i
i
nf
j
j
n
jz
L
nf
i
i
nf
j
j
n
j
nf
i
i
L
nf
j
j
n
j
0
11
1
0
111
0
1
1
ϕ
ϕανϕϕαϕϕα
( )
∫
∑∑
∫
∑
∫
∑∑
===
+
==
+
+−
−
L
nf
i
i
nf
i
i
L
nf
j
j
n
j
L
nf
i
i
nf
j
j
n
jzz
zdszdRzd
dz
d
dz
d
D
0
11
0
1
1
0
11
1
ϕϕϕαλ
ϕ
ϕ
α
(5.38)
Agrupando os termos:
( ) ( )
+
+−
∆
∫
∑∑
∫
∑∑∑
∫
∑
==
+
==
+
==
+
zd
dz
d
dz
d
Dzd
dz
d
zd
t
R
L
nf
i
i
nf
j
j
n
jzz
L
nf
i
i
nf
j
j
n
jz
nf
i
i
L
nf
j
j
n
j
0
11
1
0
11
1
1
0
1
1
ϕ
ϕ
α
ϕ
ϕανϕϕα
( ) ( )
∫
∑
∫
∑∑∑
∫
∑
=====
+
+
∆
=+
L
nf
i
i
L
nf
i
i
nf
j
j
n
j
nf
i
i
L
nf
j
j
n
j
zdszd
t
R
zdR
0
1
0
111
0
1
1
ϕϕϕαϕϕαλ
(5.39)
ou ainda:
( ) ( )
+
+−
∆
∑∑
∫∫∫
= =
+++
nf
i
nf
j
L
i
j
n
jzz
L
i
j
n
jzi
L
j
n
j
zd
dz
d
dz
d
Dzd
dz
d
zd
t
R
1 1
0
1
0
1
0
1
ϕ
ϕ
α
ϕ
ϕανϕϕα
( )
}
( )
∑
∫∫
∑
∫
= =
+
+
∆
=+
nf
i
L
i
L
i
nf
j
j
n
ji
L
j
n
j
zdszd
t
R
zdR
1
00
1
0
1
ϕϕϕαϕϕαλ
(5.40)
A equação (5.40) deve ser satisfeita para qualquer
1
0
Π∈v
. Logo, tomando-se uma
função
i
ϕ
por vez, pode-se escrevê-la em forma matricial:
[
]
{
}
{
}
FK
n
=
+
1
α
(5.41)
Universidade de São Paulo – Escola de Engenharia de São Carlos
105
em que
zdRzd
dz
d
dz
d
Dzd
dz
d
zd
t
R
K
i
L
j
L
i
j
zz
L
i
jzi
L
jij
ϕϕλ
ϕ
ϕ
ϕ
ϕνϕϕ
∫∫∫∫
++−
∆
=
0000
(5.42)
( )
∫∫
∑
+
∆
=
=
L
i
L
i
nf
j
j
n
ji
zdszd
t
R
F
00
1
ϕϕϕα
(5.43)
Assim, o Método de Elementos Finitos é aplicado para encontrar os coeficientes
multiplicadores
{
}
1
+
n
α
que satisfazem o problema algébrico definido pela equação (5.41).
Resta explicar como foi realizada a implementação computacional neste trabalho (seções 5.3 e
5.4).
5.3
Implementação Computacional
A implementação computacional das formulações de elementos finitos apresentadas
foi realizada utilizando-se um pacote de elementos finitos de domínio público denominado PZ
(DEVLOO, 1997). Esse pacote é escrito em linguagem de programação C++ e seguindo a
filosofia de orientação a objetos (LIPPMAN; LAJOIE, 1998).
O pacote ou ambiente PZ é um conjunto de classes que, interagindo, oferecem uma
determinada funcionalidade. O objetivo do ambiente PZ é implementar algoritmos avançados
de elementos finitos escondendo sua complexidade atrás de uma interface de usuário
simplificada. O ambiente PZ possui um grande tempo de desenvolvimento e uma série de
funcionalidades, tais como:
- adaptatividade hp;
- simulações unidimensionais, bidimensionais e tridimensionais;
- padrões de refinamento;
- estimativa de erros;
- métodos de resolução de sistemas lineares diretos e iterativos;
- esquemas de multigrid;
106
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
- elementos contínuos (Método de Elementos Finitos) e descontínuos (Galerkin
descontínuo);
- grande flexibilidade para implementação de novas formulações variacionais.
O ambiente PZ é constituído de módulos bastante separados, os quais combinados
originam as funcionalidades necessárias a um código de elementos finitos. Alguns módulos
podem ser destacados, devido a sua importância neste trabalho:
- geometria e espaços de interpolação;
- material; e
- algébrico.
5.3.1
Geometria, Espaços de Interpolação e Material
O ambiente PZ apresenta uma distinção muito forte entre geometria e interpolação.
Para a definição da topologia da malha, existem elementos geométricos (
TPZGeoEL)
associados a malhas geométricas (
TPZGeoMesh).
O espaço de interpolação é implementado pelos elementos computacionais,
derivados da classe base
TPZCompEl.
As classes de material implementam formulações variacionais. Sua função é calcular
a contribuição de um elemento para a matriz
[
]
K e o vetor
[
]
F do problema em um ponto de
integração. A interface é definida pela classe
TPZMaterial.
5.3.2 Algébrico
O ambiente PZ possui um módulo algébrico, onde são implementados as matrizes e
métodos de resolução de sistemas lineares. Vários modos de armazenamento de matrizes
estão disponíveis no ambiente PZ:
- matriz cheia;
- matriz em banda;
- matriz simétrica em banda;
- matriz simétrica Sky-line; e
Universidade de São Paulo – Escola de Engenharia de São Carlos
107
- matriz esparsa.
Os métodos diretos de resolução de sistemas lineares disponíveis são:
- Decomposição de Cholesky;
- Decomposição de LU;
- Decomposição de LDL
t
; e
- Método frontal.
Os métodos iterativos de resolução de sistemas lineares disponíveis são:
- Jacobi;
- SOR e SSOR;
- Métodos de Krilov: gradiente conjugado;
- GMRes: gradiente bi-conjugado.
5.4 Solução Numérica
Em problemas transientes, a escolha do passo de tempo t
∆
é uma decisão importante.
Um passo de tempo muito grande pode comprometer a exatidão da solução transiente
procurada, enquanto que um passo de tempo muito pequeno eleva o custo computacional. Em
esquemas temporais explícitos, deve-se ainda atentar que passos de tempo elevados podem
comprometer a estabilidade do método, deteriorando-se completamente a qualidade da
solução. O ideal é obter uma solução de boa qualidade com o maior passo de tempo possível,
determinando o valor do passo de tempo mais conveniente para o problema.
Em seguida, deve-se discretizar o domínio do problema, isto é, representar o domínio
em uma malha de elementos finitos.
A precisão da solução obtida e o esforço computacional requeridos são determinados
pelo grau de refinamento da malha de elementos finitos. Uma malha não–refinada tem um
número menor de nós e de graus de liberdade e dará uma precisão pior que uma malha
refinada. Entretanto, um maior número de nós na malha requer um maior esforço
computacional. Como alternativa, utilizam-se malhas adaptativas. Neste trabalho, as malhas
adaptadas adotam um maior nível de refinamento nas regiões onde a solução varia mais
108
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
fortemente e um menor refinamento nas regiões em que a solução é suave. Com isso, obtêm-
se soluções de melhor qualidade a um custo computacional reduzido.
A resolução do problema é feita iterativamente no tempo. A cada passo de tempo,
faz-se a construção do problema de elementos finitos, ou seja,
[
]
{
}
{
}
FK
=
α
e a solução do
problema
{
}
[
]
{
}
FK
1
−
=
α
. A malha é adaptada (seção 5.4.4) a cada n passos de tempo.
A solução do problema algébrico
{
}
[
]
{
}
FK
1
−
=
α
é feita através da decomposição de
Cholesky (PRESS, 2007) para o problema de Richards. A decomposição de Cholesky pode
ser aplicada em matrizes simétricas definidas positivas, sendo o método direto de inversão de
sistemas mais eficiente. O método de armazenamento de matriz utilizado é o
Sky-line
(PRESS, 2007). Para o problema do soluto, utilizou-se decomposição LU (PRESS, 2007), em
virtude da matriz K não ser simétrica. O método de armazenamento é de matriz em banda.
5.4.1 Seqüência de Simulações
Para a simulação da Equação de Richards, pelo Método dos Elementos Finitos, a
informação dos dados do problema é necessária, ou seja, os parâmetros do solo
α, n, K
s
, θ
s
e
θ
r
, o comprimento L do solo e as condições iniciais e de contorno.
A simulação do transporte de soluto requer, além dos coeficientes elencados na
equação (5.22), a solução da Equação de Richards, no passo de tempo avaliado.
O problema de transporte requer resultados obtidos na simulação da Equação de
Richards para o cálculo da velocidade da água nos poros, que é função de
θ, K (θ) e
z
∂
∂
ψ
. Por
isso, a simulação do soluto requer o cálculo anterior da simulação do problema de Richards.
As simulações da Equação de Richards e do Transporte de Soluto são realizadas no
mesmo processo. A solução de um passo de tempo do problema de Richards deve ser seguida
da solução de um passo de tempo do Transporte de Soluto. Com isso, é necessário o
armazenamento de apenas a solução do passo anterior de cada problema na memória do
computador, minimizando o uso de recursos computacionais. Os passos de tempo utilizados
podem ser diferentes para os dois problemas. O problema do transporte pode ter um passo de
tempo maior por utilizar um esquema implícito de derivada temporal. De modo a simplificar a
implementação, os passos de tempo serão adotados múltiplos um do outro, ou seja, o passo de
Universidade de São Paulo – Escola de Engenharia de São Carlos
109
tempo do problema de Transporte de Soluto será n vezes o passo de tempo do problema da
Equação de Richards. Assim, são realizados n passos de tempo do problema de Richards e em
seguida, um passo de tempo do problema de Transporte. Repete-se o processo até o número
final de iterações, atingindo-se o tempo final de simulação.
As variáveis de estado das simulações são o potencial matricial
ψ
e o produto
umidade-concentração
θ
C. A partir do potencial matricial
ψ
, obtido no problema de Richards
e dos dados do solo, é possível calcular a umidade
θ
(equação (5.4)). Dividindo-se a solução
do Transporte de Soluto
θ
C pela umidade
θ
, obtém-se a concentração C em cada tempo de
simulação.
5.4.2 Refinamento do Domínio
Neste trabalho, a malha é criada utilizando-se os recursos de refinamento do PZ. Ao
ser refinado, o elemento cria os novos nós necessários e, como conseqüência, os seus
elementos filhos. Por exemplo, considere o intervalo de
z=0 a z=L como um elemento da
malha de elementos finitos. Para criar a malha de elementos finitos, divide-se esse elemento
inicial
ndiv vezes, obtendo-se uma malha com
ndiv
2
elementos, de mesmo tamanho, ∆z =
ndiv
L
2
. Essa forma de construção da malha não é necessária, mas mostrou-se uma forma
bastante simples de geração de malha. Poder-se-ia também gerar a malha, construindo-se
todos os seus nós e elementos. A forma adotada para geração de malha, será bastante útil no
processo de adaptação de malha.
Todo elemento na malha PZ está associado a um nível. Os elementos inicialmente
possuem nível zero e quando um elemento é refinado, os novos elementos gerados (seus
filhos) possuem níveis com uma unidade a mais que seus pais (Figura 9).
110
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Figura 9 – Níveis de refinamento da malha
5.4.3 Numeração das Equações
A numeração dos nós não precisa seguir nenhum padrão. Os nós têm seus índices
dados pela ordem em que são criados e inseridos na malha. O mesmo ocorre para os
elementos. As equações do problema são associadas aos graus de liberdade dos elementos.
Elementos lineares (ordem p=1) têm dois graus de liberdade, correspondentes a suas duas
funções de base lineares, e que são compartilhados com seus vizinhos. Elementos quadráticos
(ordem p=2) possuem ainda um terceiro grau de liberdade, este não compartilhado com
nenhum vizinho.
A numeração das equações é importante, pois se reflete no tamanho da banda da
matriz
[
]
K
. Quanto menor a banda, menor o custo de armazenamento e da decomposição da
matriz.
O ambiente PZ dispõe de rotinas de renumeração de equações visando a
minimização de banda. Esse recurso foi utilizado neste trabalho.
5.4.4 Adaptação da Malha
Neste trabalho, a malha inicial da simulação é a malha uniformemente refinada, na
qual é aplicada a solução inicial
0
ψ
. Após a solução do primeiro passo de tempo (ψ
1
),
procede-se à adaptação da malha. Como critério, utilizou-se o valor do gradiente da solução
pai
filho
filho
(pai) (pai)
filho filho filho filho
nível 0
nível 1
nível 2
nível 3
Universidade de São Paulo – Escola de Engenharia de São Carlos
111
(
dz
d
ψ
para o problema de Richards e
dz
Cd
θ
para o problema do soluto). Se o gradiente no
interior de um elemento for inferior a um limite pré-estabelecido
grad
ε
, diz-se que a solução é
suave no elemento e este elemento é marcado como passível de desrefinamento. Os demais
elementos devem ser refinados.
Essa estratégia de adaptação é baseada no comportamento da solução de Richards, a
qual apresenta uma frente de propagação que pode ser mais ou menos acentuada em função
dos parâmetros do solo. Soluções que apresentam frentes de propagação acentuada requerem
um grau de refino da malha maior na região da frente de propagação, visando evitar
oscilações na solução de elementos finitos. Na frente de propagação, encontram-se os maiores
gradientes da solução.
Elementos pai unidimensionais são refinados em dois filhos unidimensionais com
um novo nó, no ponto médio do elemento pai.
No processo de refinamento, a solução dos elementos filhos é obtida por interpolação
da solução de seu elemento pai.
Os elementos marcados para refinamento são refinados até um nível limite
estabelecido, dado pela variável
ndiv na geração da malha. Elementos que já estão refinados
até o limite não são alterados. Como a solução se propaga no domínio, os elementos vizinhos
à frente de propagação também devem ser refinados, para serem adequados à solução do
próximo passo de tempo.
Observa-se que para o problema de Richards, o refinamento apenas dos elementos
vizinhos é suficiente para capturar a frente de propagação no próximo passo de tempo devido
ao esquema explícito utilizado na discretização temporal. O esquema de Euler Explícito
permite que a frente se propague apenas um elemento a cada passo de tempo. Isso é
conseqüência das restrições de estabilidade do esquema de Euler Explícito.
Já o problema de transporte de soluto foi discretizado no tempo, utilizando-se um
esquema implícito. No esquema implícito, não há limite de passo de tempo com respeito a
estabilidade do método. O limite do passo de tempo deve ser observado apenas quanto ao erro
de truncamento da derivada temporal discretizada. Dessa forma, pode ocorrer de a frente da
solução do soluto se propagar mais de um elemento por vez. Isso não foi observado nos testes
realizados neste trabalho, mas se o fosse seria necessária a adequação da estratégia de
adaptação ou a limitação do passo de tempo.
O desrefinamento só pode ocorrer, no ambiente PZ, com filhos do mesmo pai. Dois
elementos vizinhos de mesmo pai marcados para desrefinamento são convertidos em um
112
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
único elemento maior (seu pai). Um elemento marcado para desrefinamento só é eliminado e
convertido em um único elemento maior (o pai), se o elemento vizinho for do mesmo pai e se
for também marcado para ser desrefinado (Figura 9).
Se o elemento resultante do desrefinamento (elemento grosso) possui ordem de
interpolação linear (p=1), a solução no nó intermediário é desprezada, mantendo-se apenas as
soluções dos nós de extremidade. Isso pode acarretar em perda da conservação de massa, o
que pode afetar a qualidade da solução. Para minimizar a perda de massa, pode-se utilizar um
valor muito pequeno de
grad
ε
. Desse modo, apenas regiões extremamente suaves seriam
desrefinadas. O impacto dessa estratégia é uma menor capacidade de desrefinamento da
malha, o que implica em maior número de graus de liberdade e maior esforço computacional.
Uma segunda estratégia seria utilizar ordem de interpolação quadrática (p=2) ou maior. A
função de base quadrática permite a interpolação da solução dos elementos filhos para o pai
garantindo a conservação de massa. (Figura 10).
Elementos finos lineares
Elemento grosso linear (perda de
conservação de massa)
Elemento grosso quadrático (garante
conservação de massa)
Figura 1
0
–
Processo de desrefinamento
a) Dois elementos finos lineares
b) Elemento grosso linear
c) Elemento grosso quadrático
A
B
C
Universidade de São Paulo – Escola de Engenharia de São Carlos
113
5.4.4.1 Refinamento h
Com objetivo de diminuir os graus de liberdade do problema e com isso diminuir o
seu tempo de execução, a malha uniforme fina é adaptada, através do refinamento
h,
resultando em uma malha adaptada, da seguinte maneira:
a)
Calcula-se o valor da derivada da solução aproximada de u em cada elemento,
isto é,
n
E
du ;
b)
Os elementos são divididos em dois grupos: elementos com solução com forte
gradiente e elementos com solução suave. Os que têm forte gradiente são
marcados para refinamento e os de solução suave são candidatos para
desrefinamento. A divisão dos grupos é baseada no gradiente da solução no
elemento
n
E
du
comparado ao gradiente limite
grad
ε
. O valor de
grad
ε
adotado é
uma porcentagem do gradiente máximo da solução aproximada, isso é
grad
ε
=ε
dumax
.Assim, quando
n
E
du <ε
dumax
, 0 < ε < 1,
diz-se que a solução
u é suave nesta região;
c)
Desrefinam-se, segundo as possibilidades, os elementos candidatos;
d)
Os elementos marcados para refinamento são refinados até que tenham o nível
máximo
ndiv definido;
e)
Os elementos vizinhos aos elementos refinados (onde deve estar a frente de
propagação) são refinados para serem adequados à solução do próximo passo de
tempo.
A escolha adequada do parâmetro de adaptação
ε é fundamental na solução do
problema.
A ordem de aproximação p dos elementos não é modificada em todo o processo.
Todos os elementos da malha apresentam a mesma ordem de aproximação.
114
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
5.4.4.2 Projeção da solução
No processo de refinamento e desrefinamento de elementos é necessária a
transferência da solução entre a malha original e a malha adaptada.
O refinamento de elementos é implementado pelo PZ através do método Divide do
elemento computacional. O método Divide oferece, como opção, a transferência da solução
do elemento pai para seus elementos filhos. A transferência é feita por projeção
2
L
.
No desrefinamento de elementos, adotou-se o mesmo procedimento. Para a
transferência da solução dos elementos filhos para o elemento pai utilizou-se também
projeção
2
L
. O desrefinamento de elementos é implementado pelo PZ através do método
Coarsen do elemento computacional. O método Coarsen não implementa a transferência de
solução dos elementos filhos para o pai. A transferência foi implementada neste trabalho. Para
isso, o método Coarsen foi reimplementado com a opção de transferência de solução.
As Equações de Richards e de Transporte de Solutos foram implementadas neste
trabalho. O PZ permite a implementação de novas formulações externamente a seu código,
isto é, o PZ é utilizado como uma biblioteca que é instalada no computador. As equações
devem ser implementadas em classes derivadas da classe TPZMaterial. Foram implementadas
as classes TPZRichardsEquation e TPZSoluteTransport. O material do soluto aponta para a
malha de Richards para permitir a transferência de solução da malha de Richards para as
Equações de Transporte de Solutos.
O processo de refinamento e desrefinamento pode ser feito para toda a malha ou
elemento a elemento. Neste trabalho, a transferência da solução é feita elemento a elemento,
isto é, quando um elemento é refinado, seus filhos interpolam sua solução. Quando dois
elementos são agrupados, o elemento novo projeta a solução de seus filhos sobre si.
Tanto o processo de interpolação quanto o de projeção adotados garantem que a
solução dos elementos vizinhos não é modificada. Por isso, observa-se na figura 10b que a
projeção da solução para elementos com p=1 implica em perda da conservação de massa, uma
desvantagem da transferência de solução adotada. A transferência de solução entre elementos
é realizada através de projeção
2
L
. A solução no elemento recém criado
novo
u
definido no
subdomínio
(
)
21
, zz
novo
e
=Ω
é dada em função da solução no elemento original
original
u .
Universidade de São Paulo – Escola de Engenharia de São Carlos
115
Em caso de refinamento, a solução em cada um dos filhos é dada por: encontrar
∑
=
=
nfel
j
novo
j
novo
j
novo
u
1
ϕα
que satisfaça
( )
∑
∫
∑
∫
∑
= ==
=
nfel
i
z
z
novo
i
original
nfel
i
z
z
novo
i
novo
j
novo
j
nfel
j
u
1 11
2
1
2
1
ϕϕϕα
(5.44)
em que
nfel
é o número de funções
ϕ
do elemento.
Para o desrefinamento, impõe-se ainda que a solução nos limites do elemento novo
não se altere. Se o elemento novo tem limites em
z =
1
z
e z =
2
z
, os elementos filhos estão no
domínio
+
=Ω
2
,
21
11
zz
z
orig
e
+
=Ω
2
21
2
,
2
z
zz
orig
. Denominando-se
original
u
1
a solução
original em
orig
1
Ω
e
original
u
2
a solução original em
orig
2
Ω
, a projeção
2
L
da solução consiste em
encontrar
∑
=
=
nfel
j
novo
j
novo
j
novo
u
1
ϕα
tal que
)()(
111
zuzu
originalnovo
=
e
)()(
222
zuzu
originalnovo
=
que
satisfaça:
( )
∑
∫
∑
∫
∑
= ==
=
nfel
i
z
z
novo
i
original
nfel
i
z
z
novo
i
novo
j
novo
j
nfel
j
u
1 11
2
1
2
1
ϕϕϕα
(5.45)
com
Ω∈
Ω∈
=
origoriginal
origoriginal
original
zu
zu
u
22
11
,
,
e
nfel
, o número de funções
ϕ
do elemento novo.
No refinamento, a solução nos elementos filhos é idêntica a do pai. Os dois filhos
compõem um espaço de funções mais rico que o do pai. De fato, o espaço de funções do pai
está contido no espaço dos filhos, por isso a transferência é exata.
Na transferência de solução dos filhos para o pai no processo de desrefinamento, não
existe a interpolação da solução. O espaço de funções do pai não contém o espaço de funções
dos filhos e por isso, a solução transferida é diferente da original. Logo, diz-se que é feita uma
projeção.
116
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
5.4.5 Passos de Execução do Programa
Os passos de execução do programa podem ser descritos como a seguir.
1. A malha inicial da simulação de Richards é criada a partir de um único elemento. O
elemento é refinado uniformemente
ndiv vezes até que a malha tenha
ndiv
2
elementos. Aplica-
se a solução inicial
0
ψ
na malha.
2. A malha da simulação do transporte de soluto é criada refinando-se um elemento
ndiv_soluto vezes. Aplica-se a solução inicial
0
C
θ
na malha.
3. Para cada passo de tempo
t
∆
deve-se:
(a) calcular a solução no próximo passo de tempo para a simulação de Richards;
(b) Se o passo de tempo deve ser adaptado, adapta-se a malha da simulação de Richards;
(c) O passo de tempo da simulação de transporte do soluto é n
t
∆
. Se o passo de tempo for
múltiplo de n:
i. calcula-se a solução de transporte de soluto.
A. adapta-se a malha.
A adaptação de malha tem os seguintes passos. A variável do problema
u pode ser
tanto o
ψ
da simulação de Richards quando o C
θ
da simulação de transporte de soluto.
1. Calcula-se a máxima derivada de cada elemento
n
E
du . A derivada é calculada em 5 pontos
do elemento, tomando-se o maior valor encontrado.
2. A maior derivada encontrada entre todos os elementos da malha define a variável
dumax
.
3. Elementos com máxima derivada
<
n
E
du ε
dumax
, 0 < ε < 1, são marcados para
desrefinamento. Caso contrário, são marcados para refinamento.
4. Elementos marcados para desrefinamento são desrefinados.
5. Elementos marcados para refinamento são refinados.
6. Os elementos imediatamente à esquerda e à direita dos elementos marcados para
refinamento, também são refinados.
Universidade de São Paulo – Escola de Engenharia de São Carlos
117
Desrefinamento de elementos:
1. Para que um elemento seja desrefinado, seu irmão (elemento filho do mesmo pai) também
deve estar marcado para desrefinamento. Quando isso ocorre, os elementos são desrefinados.
2. Após o desrefinamento, deve-se transferir a solução dos dois elementos mais finos (filhos)
para o elemento desrefinado. A transferência é feita por projeção
2
L
da solução dos elementos
finos (filhos) para o desrefinado (pai). A projeção é feita respeitando-se o valor da solução
dos elementos vizinhos. Como o espaço de funções do elemento desrefinado é menor que o
espaço dos dois elementos finos, pode ocorrer perda de informação da solução na projeção.
Entretanto, a projeção
2
L
para elementos com p > 1 garante a conservação de massa. Isto é,
há perda de informação, mas não há perda de conservação.
Refinamento de elementos:
1. Os elementos são refinados até um nível máximo ndiv. Com isso, o tamanho do menor
elemento da malha é constante ao longo de toda a simulação, valendo
ndiv
L
2
.
2. No refinamento, um elemento é dividido em dois elementos menores.
3. A solução do elemento maior (pai) é transferida para os dois elementos refinados. A
transferência é feita por projeção
2
L
. Como o espaço de funções dos dois elementos finos
(filhos) contém o espaço de funções do elemento maior, a projeção é feita de forma exata,
sem nenhuma perda de informação.
118
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Universidade de São Paulo – Escola de Engenharia de São Carlos
119
6 VALIDAÇÃO, APLICAÇÕES NUMÉRICAS E ESTUDO
ESTATÍSTICO
Neste capítulo, são apresentadas aplicações numéricas do modelo proposto, com o
objetivo de validar o código computacional e um estudo estatístico comparativo do Modelo
Proposto com soluções disponíveis na Literatura. As simulações utilizam os recursos de
adaptatividade de malha descritos no Capítulo 5.
Na seção 6.1, são obtidas as soluções aproximadas da Equação de Richards,
utilizando-se, primeiramente, malhas uniformes e, posteriormente, malhas com refinamento
h
e, para resolver o problema da conservação da massa, observadas no processo de refino ou
desrefino da malha, utilizam-se funções de interpolação polinomial de grau 2. Na derivada
temporal, é aplicado o esquema de Euler Explícito. São elaborados os gráficos com os
resultados obtidos e comparados com os trabalhos de Celia, Bouloutas e Zarba (1990) e
Prasad, Kumar e Sekhar (2001).
As malhas utilizadas nos trabalhos de Celia, Bouloutas e Zarba (1990) e Prasad,
Kumar e Sekhar (2001) são uniformes e com funções-base lineares, enquanto que neste
trabalho, além de considerar as malhas uniformes, é empregada a adaptatividade com
refinamento
h, na aproximação espacial, e as funções-base são de grau polinomial 2.
Na seção 6.2, apresenta-se o problema de Transporte de Soluto, com os parâmetros e
condições iniciais e de contorno do trabalho de Miranda et al. (2005), para validar a Equação
de Advecção-Dispersão. Como para resolver a Equação de Advecção-Dispersão há
necessidade da solução da Equação de Richards, também foram apresentados os resultados
dessa equação através da metodologia empregada neste trabalho (seção 6.2.1). Miranda et al.
(2005) estima a concentração do soluto no solo, resolvendo numericamente em um sistema de
volumes finitos, a cada incremento de tempo
∆
t de simulação. O problema de transporte é
resolvido, neste trabalho, com malhas uniformes e considerando refinamento
h, com funções
de interpolação de grau p = 2, na solução da Equação de Richards e grau p = 3, no problema
de Transporte de Solutos. Na derivada temporal, é aplicado o esquema de Euler Implícito. São
120
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
elaborados os gráficos com os resultados obtidos e comparados com o trabalho de Miranda et
al. (2005) (seção 6.2.2).
Na seção 6.3, usando ferramentas da Estatística, foram realizadas comparações entre
os resultados obtidos no Modelo Proposto e soluções disponíveis na Literatura.
6.1 Problema de Infiltração de Água
O modelo de elementos finitos deste trabalho que resolve a Equação de Richards,
formulada em termos do potencial
ψ
, é, portanto, aplicado ao problema escolhido na
Literatura (CELIA; BOULOUTAS; ZARBA, 1990). O problema considera a infiltração da
água em uma coluna de solo homogêneo inicialmente seco. Os parâmetros do solo
considerados são
•
;35,3
1
−
= m
α
•
;2
=
n
•
smxK
s
5
1092,9
−
=
;
•
;368,0
=
s
θ
e
•
.102,0
=
r
θ
O comprimento da amostra de solo é de 1 m e as condições iniciais e de contorno são
•
mzmz 01,10)0,(
≤
≤
−
−
=
ψ
•
0,75,0),0(
>
−
=
tmt
ψ
•
0,10),1(
>
−
=
−
tmt
ψ
Na resolução numérica do problema, foi feito um estudo de adaptatividade da malha.
Foram realizados vários testes, utilizando-se, primeiramente, malhas uniformes e malhas
adaptativas.
Considerou-se primeiramente uma malha uniforme bastante refinada (chamada
malha fina), na qual é aplicada a solução inicial
0
ψ
. Após a solução do primeiro passo de
tempo (
1
ψ
), procedeu-se à adaptação da malha. Na comparação da precisão dos resultados
Universidade de São Paulo – Escola de Engenharia de São Carlos
121
obtidos, considera-se a solução aproximada obtida a partir da malha uniforme fina como
sendo a solução de referência.
6.1.1 Malha Uniforme – Comparação com Resultados Publicados
O domínio Ω = [-1, 0] foi subdividido uniformemente em 2
6
= 64 elementos e cada
elemento
Ω
e
possui tamanho ∆z = 1/64. O tempo de simulação é de 1 dia. O passo de tempo
utilizado é ∆
t = 1 segundo. Foram utilizadas funções polinomiais por partes de grau p = 1,
inicialmente, devido à simplicidade em se trabalhar com funções lineares e por se tratar de
malha uniforme (Figura 11). O tempo de processamento para obtenção da solução é de 1min
23,169s. Essa malha fina possui 65 nós.
A Figura 11a apresenta o gráfico da solução aproximada do problema, considerando-se
a malha uniforme fina, com p = 1. O gráfico da Figura 11b representa os 65 nós da malha
uniforme de elementos finitos, considerando-se p = 1.
122
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Figura 11 - a) Solução aproximada em uma malha uniforme refinada com 2
6
elementos e p = 1.
b) ilustração dos nós da malha uniforme refinada com 2
6
elementos e p = 1.
Universidade de São Paulo – Escola de Engenharia de São Carlos
123
A solução obtida, considerando-se a malha uniforme fina com
6
2
elementos e p = 1,
mostra-se equivalente ao modelo de Celia, Bouloutas e Zarba (1990) que resolveram a
Equação de Richards, formulada em termos de potencial
ψ
, também pelo Método de
Elementos Finitos. Na aproximação espacial, os autores consideram as malhas uniformes,
funções-base lineares por partes e na derivada temporal, é aplicado o esquema de Euler
Implícito (Figura 12).
Figura 12 - Comparação da solução obtida com a malha uniforme fina (
6
2
elementos e p = 1) e a solução de
Celia, Bouloutas e Zarba (1990).
6.1.2 Adaptatividade
A malha uniforme fina é adaptada com o objetivo de diminuir a quantidade de graus
de liberdade do problema, reduzindo o tempo de execução do programa. Entretanto, isso deve
ser feito sem que haja perda substancial da qualidade da solução.
Após a solução do primeiro passo de tempo (
1
ψ
), procedeu-se à adaptação da malha.
Os parâmetros de adaptação
ε (seção 5.4.4.1) utilizados foram ε = 0,01 e ε = 0,001. Os
resultados considerando
ε = 0,01 e ε = 0,001 são apresentados nas Figuras 13 e 14,
respectivamente.
124
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Figura 13 - a) Solução aproximada em uma malha adaptada com p=1 e ε = 0,01.
b) ilustração dos nós da malha no último passo de tempo da simulação.
Universidade de São Paulo – Escola de Engenharia de São Carlos
125
Figura 14 - a) Solução aproximada em uma malha adaptada com p=1 e ε = 0,001.
b) ilustração dos nós da malha no último passo de tempo da simulação.
126
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
As Figuras 15 e 16 comparam as soluções considerando os parâmetros de adaptação
ε = 0,01 e ε = 0,001 com as soluções obtidas com a malha uniforme fina (2
6
elementos e
p = 1
)
, respectivamente.
Pode-se observar na Figura 15, que os resultados da simulação considerando-se
funções polinomiais com p = 1 e
ε = 0,01 não correspondem aos resultados obtidos com a
malha uniforme fina. Observa-se que a frente de umedecimento da malha adaptada está
“atrasada” em relação à da malha uniforme. Isso se deve à perda de massa da estratégia de
adaptação adotada quando p = 1, conforme mencionado na seção 5. 4. 4.
Já a solução com a malha adaptada e
ε = 0,001 mostrou-se bastante adequada, embora
o número de nós da malha seja muito próximo do número de nós da malha uniforme, ou seja,
62 nós (Figura 16).
Universidade de São Paulo – Escola de Engenharia de São Carlos
127
Figura 15 - a) Malha uniforme fina com 2
6
elementos e p = 1 versus malha adaptada com p = 1 e ε = 0,01.
b) ilustração dos nós da malha uniforme fina com 2
6
elementos e p = 1 versus malha adaptada com
p = 1 e ε = 0,01.
128
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Figura 16 - a) Malha uniforme fina com 2
6
elementos e p = 1 versus malha adaptada com p = 1 e ε = 0,001.
b) ilustração dos nós da malha uniforme fina com 2
6
elementos e p = 1 versus malha adaptada
com p = 1 e ε = 0,001
Universidade de São Paulo – Escola de Engenharia de São Carlos
129
Como a solução com a malha adaptada e ε = 0,001 mostrou-se bastante adequada,
embora o número de nós da malha seja muito próximo do número de nós da malha uniforme,
a Figura 17 mostra a comparação com a solução de Celia, Bouloutas e Zarba (1990).
Figura 17 - Comparação da solução obtida com a malha adaptada (p = 1 e ε = 0,001) e a solução de Celia,
Bouloutas e Zarba (1990).
6.1.3 Utilizando Ordem de Aproximação Polinomial p = 2
Como mencionado na seção 5.4.4, uma abordagem para minimizar a perda da
conservação de massa é a utilização de ordem de interpolação maior do que 1. Utilizando-se
ordem quadrática (p = 2), foram obtidos os resultados das Figuras 18, 19 e 20. Nelas, o
parâmetro de adaptação
ε assume os valores ε = 0,1, ε = 0,01 e ε = 0,001, respectivamente.
Observa-se que a solução com
ε = 0,01 e ε = 0,001 são bastante satisfatórias (Figuras 19 e
20), sem o problema de perda de massa da estratégia de adaptação, já que a frente de
propagação é “igual” à da malha uniforme. A malha uniforme é a referência sempre e
qualquer adaptação que seja feita, com o objetivo de tornar a simulação mais rápida, não deve
piorar a qualidade da aproximação. Já a simulação com ε = 0,1 mostrou-se inadequada,
possivelmente por ser uma malha demasiadamente “grosseira” (Figura 18).
130
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Figura 18 - a) Solução aproximada em uma malha adaptada com p = 2 e ε = 0,1.
b) ilustração dos nós da malha no último passo de tempo da simulação.
Universidade de São Paulo – Escola de Engenharia de São Carlos
131
Figura 19 - a) Solução aproximada em uma malha adaptada com p = 2 e ε = 0,01.
b) ilustração dos nós da malha no último passo de tempo da simulação.
132
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Figura 20 - a) Solução aproximada em uma malha adaptada com p = 2 e ε = 0,001.
b) ilustração dos nós da malha no último passo de tempo da simulação.
Universidade de São Paulo – Escola de Engenharia de São Carlos
133
As Figuras 21 e 22 comparam os resultados das malhas adaptadas com ordem de
interpolação p = 2 e
ε = 0,01 e ε = 0,001, respectivamente, com os resultados da malha
uniforme fina com
6
2
elementos e p = 1.
Figura 21 - a) Malha uniforme fina com 2
6
elementos e p = 1 versus malha adaptada com p = 2 e ε = 0,01.
b) ilustração dos nós da malha uniforme fina com 2
6
elementos e p = 1 versus malha adaptada
com p = 2 e ε = 0,01.
134
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Figura 22 – a) Malha uniforme fina com 2
6
elementos e p = 1 versus malha adaptada com p = 2 e ε = 0,001.
b) ilustração dos nós da malha uniforme fina com 2
6
elementos e p = 1 versus malha adaptada
com p = 2 e ε = 0,001
Universidade de São Paulo – Escola de Engenharia de São Carlos
135
6.1.4 Discussões entre os Resultados Obtidos
Dentre os resultados obtidos, as soluções aproximadas em uma malha adaptada
(refinamento
h) com
p = 2 e
ε = 0,01 e ε = 0,001 apresentaram melhor aproximação em
relação à malha uniforme fina (Figuras 21 e 22), embora, a malha adaptada (refinamento
h)
com
p = 2 e
ε = 0,001 apresentou melhor aproximação em relação à malha uniforme fina com
2
6
elementos e p = 1 (Figura 22).
O número de nós para a obtenção da solução considerando-se a malha adaptada com
p = 2 e
ε = 0,01 e a adaptada com p = 2 e ε = 0,001, no último passo de tempo, é,
respectivamente, 23 e 39 nós, e, o tempo de processamento é, respectivamente, 1 min 28,390s
e 1 min 42,474s. (Figuras 19 e 20). A solução obtida considerando-se a malha adaptada com
p
= 2 e
ε = 0,01 é uma boa aproximação em relação à malha uniforme fina devido à quantidade
bem menor de nós em relação à malha adaptada com
p = 2 e
ε = 0,001 (redução de 41%,
aproximadamente, na quantidade de nós) (Figura 21).
O tamanho do menor elemento da adaptada é igual ao da uniforme, por isso é usado
o mesmo passo de tempo para essas duas malhas e a distância máxima entre os nós varia a
cada passo de tempo. Se na malha adaptada, o tamanho do elemento fosse menor, o passo de
tempo teria que diminuir, devido à condição de CFL. Na malha adaptada, são removidos os
nós desnecessários, diminuindo o custo computacional. Portanto, a qualidade das
aproximações é igual nas malhas uniforme e adaptada. Poder-se-ia ter uma outra abordagem,
ou seja, manter o número de nós constante na malha, o que implicaria na malha adaptada
elementos de tamanho menor que na malha uniforme, melhorando a qualidade da
aproximação com o mesmo custo computacional. Não é o que foi feito, mas nas duas
abordagens, obter-se-ia o mesmo resultado. No modelo proposto, como o custo
computacional diminui, pode-se aumentar o
ndiv e ter malhas mais finas, com um tempo
pequeno para processamento da simulação.
136
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
6.1.5 Comparação com Celia,
Bouloutas e Zarba (1990)
Como a solução obtida considerando-se a malha adaptada com
p = 2 e
ε = 0,01 é uma
boa aproximação em relação à malha uniforme fina, é interessante comparar com o modelo de
Celia, Bouloutas e Zarba (1990) (Figura 23).
Figura 23 - Comparação da solução obtida com a malha adaptada (p = 2 e ε = 0,01) e a solução de Celia,
Bouloutas e Zarba (1990).
6.1.6 Comparação com Prasad, Kumar e Sekhar (2001)
Na figura 24 é feita a comparação da solução obtida considerando-se a malha
adaptada com
p = 2 e
ε = 0,01 com o modelo de Prasad, Kumar e Sekhar (2001) que valida
seu modelo com os parâmetros do solo, condições iniciais e de contorno do trabalho de Celia,
Bouloutas e Zarba (1990), utilizando malhas uniformes com funções-base lineares
lagrangeanas.
Universidade de São Paulo – Escola de Engenharia de São Carlos
137
Figura 24 - Comparação da solução obtida com a malha adaptada (p = 2 e ε = 0,01) e a solução de Prasad,
Kumar e Sekhar (2001).
A Figura 25 faz uma comparação de resultados da metodologia empregada neste
trabalho, usando refinamento
h e polinômio de interpolação quadrático, evitando a perda de
massa, diminuindo os graus de liberdade e o custo computacional, com os dois modelos
escolhidos na Literatura.
Figura 25 - Comparação da solução obtida com a malha adaptada (p = 2 e ε = 0,01) e as soluções de Celia,
Bouloutas e Zarba (1990) e de Prasad, Kumar e Sekhar (2001).
138
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
6.2 Problema de Transporte de Soluto
Nesta seção, é apresentada uma aplicação numérica do modelo proposto no Capítulo
5, que descreve o movimento de água no solo e o transporte de solutos, na simulação do
deslocamento do íon potássio em colunas de solo não-saturado. Os dados das análises físico-
hídricas do Latossolo Vermelho-Amarelo, fase arenosa, os parâmetros da curva de retenção,
ajustada segundo o modelo de van Genuchten (1980), bem como os parâmetros de transporte
do potássio, necessários para resolução do problema, foram obtidos do trabalho de Miranda et
al. (2005).
A simulação do transporte de soluto (seção 6.2.2) requer a solução da Equação de
Richards, no passo de tempo avaliado, para o cálculo da velocidade da água nos poros, além
dos coeficientes elencados na equação (5.22). Por isso, a simulação do soluto requer o cálculo
anterior da simulação do problema de Richards (seção 6.2.1).
6.2.1 Movimento de Água no Solo (Equação de Richards)
Para a simulação da Equação de Richards, pelo Método dos Elementos Finitos, a
informação dos dados do problema é necessária, ou seja, os parâmetros do solo
α, n, K
s
, θ
s
e
θ
r
, o comprimento L do solo e as condições iniciais e de contorno.
Os parâmetros de simulação (MIRANDA et al., 2005) utilizados são
;/0216,0
33
0
cmcm=
θ
;/0216,0
33
cmcm
i
=
θ
;0
=
r
θ
;/443,0
33
cmcm
s
=
θ
;457,5
1
−
= hcmK
S
;0449,0
1
−
= cm
α
;6732,3
=
n e
.727758,0
=
m
Universidade de São Paulo – Escola de Engenharia de São Carlos
139
As condições iniciais e de contorno são dadas por
•
007,6524,68)0,(
≤
≤
−
−
=
zcmz
ψ
•
0,2476,6),0(
>
−
=
tcmt
ψ
•
0,6524,68),70(
>
−
=
−
tcmt
ψ
Na resolução numérica do problema de movimento de água no solo, utiliza-se o esquema
de Euler Explícito (BEAR, 2007) (seção 2.10.1). O domínio
Ω = [-70, 0] foi particionado, na malha
uniforme, em 128 elementos de igual tamanho. O tempo total de simulação é
t = 1,75 h, em 450000
passos de tempo com
∆t = 1,75/450000 horas = 0,014 segundos.
Foram realizadas duas simulações, uma com a malha uniforme de 128 elementos e outra
com malha adaptada. O tamanho do menor elemento da malha adaptada é o mesmo dos elementos
da malha uniforme, ou seja, 70/128 cm, e, o do maior varia a cada passo de tempo. O parâmetro de
adaptação
ε (seção 5.4.4.1) foi de ε = 0,01. O processo de adaptação da malha é realizado a cada
1000 passos de tempo. Para as duas simulações, adotou-se ordem de aproximação p = 2 constante
para todos os elementos da malha.
As Figuras 27 e 28 mostram as soluções obtidas para o potencial matricial
ψ com a
malha uniforme com
7
2
elementos, p = 2 e a malha adaptada, com p = 2 e ε = 0,01 no instante
final de simulação (Figura 27) e em vários instantes de simulação (Figura 28). Os resultados
indicam a equivalência de solução para as duas malhas utilizadas, ou seja, a malha uniforme
com
7
2
elementos, p = 2, e a malha adaptada com p = 2 e ε = 0,01, observando que a malha
adaptada tem um número menor de graus de liberdade, o que acarreta na redução do tempo de
execução computacional. O número de elementos da malha adaptada, nas diversas iterações,
varia de 9 a 26 elementos, portanto, de 10 a 27 nós. O tempo de processamento para encontrar
a solução utilizando-se a malha uniforme é de 21 min 38s e o da malha adaptada, de 3 min e 5
s, redução de 86 %, aproximadamente. O consumo de memória é 6 Mb quando se utiliza a
malha uniforme e 5,800 Mb, quando é adaptada (Tabela 2).
Tabela 2 – Tempo de processamento e consumo de memória na solução da Equação de Richards
Malhas Tempo de Processamento Memória
Uniforme 22 min 38 s 6 Mb
Adaptada 3 min 5 s 5,8 Mb
140
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
O número de elementos da malha adaptada na iteração 1000, 2000, 3000,..., até
450000 é apresentado na Figura 26.
Figura 26 – Número de elementos da malha adaptada da Equação de Richards na iteração 1000, 2000,
3000,...,450000.
Universidade de São Paulo – Escola de Engenharia de São Carlos
141
Figura 27 - Potencial Matricial ψ no instante final de simulação: passo de tempo 450000; 1,75 h.
a) malha uniforme com
7
2
elementos, com p = 2. b) malha adaptada com p = 2 e ε = 0,01.
142
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Figura 28 - Potencial Matricial ψ em vários instantes de simulação. a) malha uniforme com
7
2
elementos e p=2.
b) malha adaptada com p = 2 e ε = 0,01.
Universidade de São Paulo – Escola de Engenharia de São Carlos
143
As Figuras 29 e 30 mostram as soluções obtidas para a umidade
θ
com a malha
uniforme com
7
2
elementos, p = 2 e a malha adaptada, com p = 2 e ε = 0,01 no instante final
da simulação (Figura 29) e em vários instantes de simulação (Figura 30). Os resultados
indicam a equivalência de solução para as duas malhas utilizadas, ou seja, a malha uniforme
com
7
2
elementos, p = 2, e a malha adaptada com p = 2 e ε = 0,01.
Foram realizados vários testes, utilizando-se malhas uniformes e malhas adaptadas,
com refinamento
h.
Os resultados da umidade
θ
(Figuras 29 e 30) mostram que o umedecimento se deu
até cerca de 35 cm de profundidade após 1 hora e 45 minutos.
144
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Figura 29 - Umidade
θ
no instante final de simulação (passo de tempo 450000; 1,75 h). a) malha uniforme com
7
2
elementos, com p = 2. b) malha adaptada com p = 2 e ε = 0,01.
Universidade de São Paulo – Escola de Engenharia de São Carlos
145
Figura 30 - Umidade
θ
em vários instantes simulados. a) malha uniforme com
7
2
elementos e p = 2.
b) malha adaptada com p = 2 e ε = 0,01.
146
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
O gráfico da umidade x profundidade (Figura 29) é equivalente aos resultados
apresentados em Miranda et al. (2005) (Figuras 31 e 32), que resolveu a Equação de Richards,
formulada em termos de umidade
θ
.
Figura 31 - Comparação do perfil de umidade obtido por Miranda et al. (2005) (Observado experimentalmente e
simulado) e a solução obtida, no presente trabalho, com malha uniforme (
7
2
elementos e p = 2).
Figura 32 - Comparação do perfil de umidade obtido por Miranda et al. (2005) (Observado experimentalmente e
simulado) e a solução obtida, no presente trabalho, com malha adaptada (p = 2 e ε = 0,01).
Universidade de São Paulo – Escola de Engenharia de São Carlos
147
6.2.2 Transporte de Soluto (Equação de Advecção-Dispersão)
No problema de transporte de soluto, busca-se a evolução temporal da concentração
C do soluto no solo, nesse caso, o potássio. A equação depende de variáveis obtidas na
simulação da Equação de Richards, apresentada na seção 6.2.1. Como variável de estado do
problema de transporte de soluto, adotou-se
C
θ
(seção 5.2.1). A utilização de
C
θ
como
variável do problema ao invés de
C justifica-se devido à suavidade da solução nesta variável,
como será possível observar nos resultados do problema (Figuras 34 e 35).
Os parâmetros utilizados foram os mesmos da seção 6.2.1 acrescidos de
687,4
=
R
cm,
l
616621
=
α
;
;0
=
λ
0
=
s
; e
12
min7303,2
−
= cmD
zz
.
conforme definidos na seção 5.2.
As condições iniciais e de contorno para a concentração
C são
•
(
)
335
/50/1050, mkgcmkgzC =×=
−
para
070
≤
≤
−
z
•
(
)
335
/50/105,70 mkgcmkgtC =×=−
−
para 0
>
t
•
(
)
334
/500/105,0 mkgcmkgtC =×=
−
para
0
>
t
As condições iniciais e de contorno para a variável de estado concentração x
umidade
C
θ
são
•
(
)
3365
/08,1/1008,11050216,00, mkgcmkgzC =×=××=
−−
θ
para
070
≤
≤
−
z
•
(
)
3365
/08,1/1008,11050216,0,70 mkgcmkgtC =×=××=−
−−
θ
para 0
>
t
•
(
)
3364
/5,221/105,221105443,0,0 mkgcmkgtC =×=××=
−−
θ
para
0
>
t
Para o problema do transporte de soluto, foi possível utilizar o esquema de Euler
Implícito (BEAR, 2007) (seção 2.10.1.1), sem maiores dificuldades, devido à linearidade da
148
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
equação. O uso do esquema implícito é muito vantajoso, pois não implica em restrição ao
tamanho do passo de tempo. Desta maneira, foi utilizado um passo de tempo para o transporte
do soluto 1000 vezes maior do que para a Equação de Richards, isto é,
∆t = 14 segundos. O
número de passos de tempo executado é de 450 com um instante final de simulação de 1,75
horas.
O domínio
Ω = [-70, 0] foi particionado, na malha uniforme, em 64 elementos de
mesmo tamanho.
Foram realizadas duas simulações: uma com a malha uniforme de 64 elementos e
outra com malha adaptada. O tamanho do menor elemento da malha adaptada é o mesmo dos
elementos da malha uniforme, ou seja, 70/64 cm, e, o maior varia a cada passo de tempo. O
número de elementos nos passos de tempo do transporte de soluto é apresentado na Figura 33.
Figura 33 – Número de elemento da malha adaptada do transporte de solutos nos passo de tempo 1,2,3,...,450.
O parâmetro de adaptação
ε foi de 0,01. O processo de adaptação da malha é
realizado em todo passo de tempo. Devido à suavidade da solução
C
θ
, adotou-se, para as
duas simulações, ordem de aproximação p = 3, constante para todos os elementos da malha. A
suavidade da solução permite o emprego de maior ordem de aproximação sem receio das
oscilações numéricas presentes na solução da Equação de Richards.
As Figuras 34 e 35 mostram as soluções obtidas para variável de estado concentração
x umidade
C
θ
com a malha uniforme com
6
2
elementos, p = 3 e a malha adaptada com p = 3
e
ε = 0,01 no instante final de simulação (Figura 34) e em vários instantes de simulação
(Figura 35). Os resultados indicam também que a escolha de
C
θ
como variável de estado é
acertada (Figuras 34 e 35), devido à suavidade da solução.
Universidade de São Paulo – Escola de Engenharia de São Carlos
149
Figura 34 -
θ
C no instante final de simulação (passo de tempo 450; 1,75 h). a) malha uniforme com
6
2
elementos, com p = 3. b) malha adaptada, com p = 3 e ε = 0,01.
150
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Figura 35 -
θ
C em vários instantes de simulação. a) malha uniforme com
6
2
elementos e p = 3.
b) malha adaptada, com p = 3 e ε = 0,01.
Universidade de São Paulo – Escola de Engenharia de São Carlos
151
As Figuras 36 e 37 mostram as soluções obtidas para variável concentração C com a
malha uniforme com
6
2
elementos, p = 3 e a malha adaptada com p = 3 e ε = 0,01 no instante
final de simulação (Figura 36) e em vários instantes de simulação (Figura 37).
O umedecimento, que se deu até cerca de 35 cm de profundidade, após 1 hora e 45
minutos, também pode ser observado nos resultados da concentração do soluto
C (Figura
36). Pode-se observar que a concentração em profundidade superior a 35 cm mantém-se
inalterada, com valor igual a
3
/50 mkg
. É interessante observar que para profundidades
menores, entre 20 e 30 cm de profundidade, no instante final de simulação, a concentração
assume valores inferiores à concentração inicial (Figuras 36 e 37) e as maiores concentrações
de potássio estão nas camadas superiores. Esse resultado deve-se ao fato do potássio em
solução ter sido adsorvido pelo solo, mostrando que a frente de molhamento que avançou para
profundidades maiores, provavelmente, apresentava concentração menor desse soluto.
152
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Figura 36 - Concentração C no instante final de simulação (passo de tempo 450; 1,75 h). a) malha uniforme com
6
2
elementos com p = 3. b) malha adaptada, com p = 3 e ε = 0,01.
Universidade de São Paulo – Escola de Engenharia de São Carlos
153
Figura 37 - Concentração C em vários instantes de simulação. a) malha uniforme com
6
2
elementos e p = 3.
b) malha adaptada com p = 3 e ε = 0,01.
154
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
A solução da concentração C no instante final de simulação (passo de tempo 450;
1,75 h), tanto com a malha uniforme (
6
2
elementos com p = 3) quanto com a malha adaptada
(p = 3 e
ε = 0,01), observada nas Figuras 36 e 37, é equivalente aos resultados numéricos e
experimentais de Miranda et al. (2005) (Figuras 38 e 39).
Figura 38 - Comparação do perfil de concentração do íon potássio obtido por Miranda et al. (2005) (Observado
experimentalmente e simulado) e a solução obtida, no presente trabalho, com malha uniforme
(
6
2
elementos e p = 3).
Figura 39 - Comparação do perfil de concentração do íon potássio obtido por Miranda et al. (2005) (Observado
experimentalmente e simulado) e a solução obtida, no presente trabalho, com malha adaptada
(p = 3 e ε = 0,01).
Universidade de São Paulo – Escola de Engenharia de São Carlos
155
A solução do problema de transporte torna-se mais simples do que a solução da
Equação de Richards devido:
a)
a equação ser linear e permitir o uso do esquema implícito sem a necessidade de
resolver sistemas de equações não-lineares.
b)
o esquema implícito permitir refinar a malha tanto quanto se queira sem a
preocupação de perda de estabilidade do método.
c)
a solução
C
θ
ser suave e não induzir nenhuma oscilação numérica (Figuras 34 e
35). Devido a isso, nesse teste, optou-se por utilizar ordem de aproximação
maior do que para a Equação de Richards (p = 3).
As oscilações presentes na solução da concentração
C são devidas exclusivamente à
solução da Equação de Richards. Observa-se, portanto, que especial cuidado deve ser dado ao
problema da Equação de Richards.
6.2.3 Tempo de Execução e Consumo de Memória
O tempo de execução e o consumo de memória referem-se ao processo que executa
as simulações do problema de Richards e do Transporte de Soluto, em seqüência. Foram
realizadas duas simulações, uma com malhas uniformes e outra com malhas adaptadas.
Usando as malhas uniformes para o problema de Richards e o Transporte de Soluto,
o tempo de execução total foi de 22 minutos e a memória utilizada de 6164 Kb. Com as
malhas adaptadas, o tempo de execução total foi de 3 minutos e 27 segundos consumindo
5876 Kb de memória, redução de, aproximadamente, 84 % no tempo de execução (Tabela 3).
As simulações foram executadas em um computador com processador Intel®
Core™2 Duo T8300 de 2.40 GHz com 3 Gb de memória RAM, em sistema operacional
Linux Ubuntu 32 bits.
No consumo de memória, não estão incluídos os resultados de todos os passos de
tempo, uma vez que estes são escritos em disco. Apenas os vetores solução do passo de tempo
corrente e do passo de tempo anterior de cada um dos problemas (Equação de Richards e
Transporte de Soluto) são mantidos em memória.
156
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Tabela 3 – Tempo de processamento e consumo de memória na solução das Equações de Richards e de
Transporte de Solutos
Malhas Tempo de Processamento Memória
Uniforme 22 min 6,164 Mb
Adaptada 3 min 27 s 5,876 Mb
6.3 Estudo Estatístico comparativo do Modelo Proposto com dados da Literatura
Na seção 6.3 é apresentado um estudo estatístico comparativo da solução numérica
da Equação de Richards e da Equação de Transporte de Solutos obtidas pela adaptatividade
com refinamento
h na malha de elementos finitos com as soluções de Celia, Bouloutas e
Zarba (1990), de Prasad, Kumar e Sekhar (2001) e de Miranda et al. (2005).
No estudo da Regressão, a finalidade é determinar a dependência de uma variável
chamada variável dependente em relação a uma outra variável chamada variável
independente.
Usando aplicativos Statistica e Bioestat de Estatística, foram feitas as regressões
lineares das variáveis em questão, apresentando as curvas e equações de regressão linear,
verificando-se, pelo coeficiente de determinação, a medida da proporção da variabilidade em
uma variável que é explicada pela variabilidade da outra (MORETTIN; BUSSAB, 2004).
6.3.1 Regressão Linear dos potenciais matriciais na Malha Uniforme Fina e no Método
de Celia, Bouloutas e Zarba
O objetivo desta seção é apresentar um estudo comparativo da solução numérica da
Equação de Richards pelo Método Proposto neste trabalho usando a Malha Uniforme Fina
(2
6
elementos e p = 1) com a solução obtida pelo Método de Celia, Bouloutas e Zarba
(1990). A partir da Figura 12 é elaborada a Tabela 4, usada na construção da Curva de
Regressão Linear.
Universidade de São Paulo – Escola de Engenharia de São Carlos
157
Tabela 4 – Potenciais Matriciais versus Profundidade dos Métodos a serem comparados
Profund.
-0,05 -0,10 -0,15 -0,2 -0,3 -0,4 -0,5 -0,55 -0,55 -0,56 -0,57 -0,6 -0,65 -0,95
Celia
1
ψ
-0,75 -0,75 -0,80 -0,80 -0,76 -1 -1,3 -2,5 -3 -5,5 -7,5 -10 -10 -10
Uniforme
u
ψ
-0,75 -0,75 -0,75 -0,75 -0,75 -1 -1,5 -2,5 -3 -5,5 -7,5 -10 -10 -10
A Curva de Regressão Linear, tendo como variável dependente, o potencial matricial
u
ψ
na Malha Uniforme Fina (2
6
elementos e p = 1) e como variável independente, o
potencial matricial
1
ψ
obtido pelo Método de Celia, Bouloutas e Zarba (1990), é dada pela
Figura 40. O Coeficiente de Determinação é R
2
= 0,9998 e a equação de Regressão Linear é
u
ψ
= - 0,0087 + 0,9993
1
ψ
(6.1)
Figura 40 - Curva de Regressão Linear considerando o potencial matricial na Malha Uniforme Fina e no Método
de Celia, Bouloutas e Zarba
O Coeficiente de Determinação, com valor R
2
= 0,9998, bem próximo de 1, implica
que o potencial matricial obtido pelo Método de Celia, Boulotas e Zarba (1990) é explicativo
para o potencial matricial do Modelo Proposto usando a Malha Uniforme Fina.
158
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
6.3.2 Estudo Comparativo da solução de Richards pelo Modelo Proposto e soluções da
Literatura
O objetivo desta seção é apresentar um estudo comparativo da solução numérica da
Equação de Richards, em termos do potencial matricial, pelo Modelo Proposto, usando a
Malha Adaptada (
ε
01,0
=
e p = 2), com as soluções obtidas pelos Métodos de Celia, Boulotas
e Zarba (1990) e de Prasad, Kumar e Sekhar (2001) em determinados valores da
profundidade.
A partir da Figura 25 é elaborada a Tabela 5, usada na construção das Curvas de
Regressão Linear.
Tabela 5 – Potenciais Matriciais versus Profundidade dos Métodos a serem comparados
Profund
-0,15 -0,30 -0,50 -0,53 -0,55 -0,56 -0,575 -0,58 -0,6 -0,65 -0,85 -0,90
Celia:
1
ψ
-1 -1 -1,25 -1,6 -3 -4 -8,5 -9,25 -10 -10 -9,9 -10
Prasad:
2
ψ
-1 -1,1 -1,25 -1,65 -1,75 -1,8 -6,5 -8 -9,8 -9,8 -9,8 -9,9
Adapt:
3
ψ
-0,8 -0,9 -1,5 -1,75 -3 -5 -9,5 -9,5 -9,75 -9,8 -10 -10
6.3.2.1 Regressão Linear dos potenciais matriciais do Modelo Proposto e do Método de
Celia, Bouloutas e Zarba
A Curva de Regressão Linear, tendo como variável dependente, o potencial matricial
3
ψ
, do Modelo Proposto (Malha Adaptada (ε 01,0
=
e p = 2)) e como variável independente, o
potencial matricial
1
ψ
obtido pelo Método de Celia, Boulotas e Zarba (1990), é dada pela
Figura 41.
O Coeficiente de Determinação é R
2
= 0,9902 e a equação da Curva de Regressão
Linear é
3
ψ
= - 0,132 +1,0040
1
ψ
(6.2)
Universidade de São Paulo – Escola de Engenharia de São Carlos
159
Como R
2
= 0,9902, muito próximo de 1, significa que o potencial matricial
1
ψ
explica
muito bem o potencial matricial
3
ψ
( Figura 41).
Figura 41 - Curva de Regressão Linear com os potenciais matriciais obtidos pelo Modelo Proposto (Malha
Adaptada) e pelo Método de Celia, Boulotas e Zarba
6.3.2.2 Regressão Linear dos potenciais matriciais do Modelo Proposto e do Método de
Prasad, Kumar e Sekhar
A Curva de Regressão Linear, tendo como variável dependente, o potencial matricial
3
ψ
do Modelo Proposto (Malha Adaptada (ε 01,0
=
e p = 2)) e como variável independente, o
potencial matricial
2
ψ
obtido pelo Método de Prasad, Kumar e Sekhar (2001), é dada pela
Figura 42. O Coeficiente de Determinação é R
2
= 0,9233 e a equação da Curva de Regressão
Linear é
3
ψ
=- 0,7485 +0,9907
2
ψ
(6.3)
Como R
2
= 0,9233, muito próximo de 1, significa que o potencial matricial
2
ψ
explica
muito bem o potencial matricial
3
ψ
( Figura 42).
160
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Figura 42 - Curva de Regressão Linear com os potenciais matriciais obtidos pelo Modelo Proposto (Malha
Adaptada) e Método de Prasad, Kumar e Sekhar
6.3.3 Estudo Comparativo da solução da Equação de Richards usando o Método de
Miranda et al. e o Modelo Proposto
O objetivo desta seção é apresentar um estudo comparativo da solução numérica da
Equação de Richards, formulada em termos da umidade, considerando o Modelo Proposto,
usando a Malha Adaptada (
ε
01,0
=
e p = 2) e o perfil de umidade obtido por Miranda et al.
(2005) (Observado experimentalmente e simulado - MIDI).
A partir da Figura 32 é elaborada a Tabela 6 usada na construção da Curva de
Regressão Linear.
Tabela 6 – Umidade versus Profundidade para a Malha Adaptada e Miranda et al. (2005) (Observada
experimentalmente e MIDI)
Profundidade
-70 -50 -44 -40 -34 -30 -26 -20 -10 -8
Adaptada
adapt
θ
2 2 2 1 42 36 41 43 44 44
Observada
Obs
θ
2 2 1 6 18 30 38 40 39 41
MIDI
MIDI
θ
2 2 8 16 2 42 44 44 44 44
Universidade de São Paulo – Escola de Engenharia de São Carlos
161
6.3.3.1 Regressão Linear das umidades do Modelo Proposto e do Método de Miranda et
al. (experimental)
O Coeficiente de Determinação é R
2
= 0,8452 e a Curva de Regressão Linear, tendo
como variável dependente, a umidade do Modelo Proposto usando a Malha Adaptada
(
ε
01,0
=
e p = 2) e
como variável independente, a umidade obtida pelo Método de Miranda et
al. (2005) (Observada experimentalmente), é dada pela Figura 43. A Equação de Regressão
Linear é dada por
adapt
θ
= 2,5904 +1,0773
Obs
θ
(6.4)
Figura 43 - Curva de Regressão Linear considerando a umidade obtida pelo Modelo Proposto (Malha
Adaptada) com a umidade obtida experimentalmente por Miranda et al.
6.3.3.2 Regressão Linear das umidades do Modelo Proposto e do Método de Miranda et
al. (MIDI)
O Coeficiente de Determinação é R
2
= 0,8005 e a Equação de Regressão Linear, tendo
como variável dependente, a umidade do Modelo Proposto considerando a Malha Adaptada
(
ε 01,0
=
e p = 2) e
como variável independente, a umidade obtida pelo Método de Miranda et
al. (2005) (MIDI), é
adapt
θ
= 0,5706 +0,9649
MIDI
θ
(6.5)
162
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
A curva de Regressão Linear é dada pela Figura 44.
Figura 44 - Curva de Regressão Linear considerando a umidade obtida pelo Modelo Proposto (Malha
Adaptada) com a umidade obtida por Miranda et al. (MIDI)
6.3.4 Estudo Comparativo das Concentrações de Potássio dadas pelo Modelo Proposto e
pelo Método de Miranda et al.
São apresentadas a concentração de potássio e a Curva de Regressão, tendo como
variável dependente, a concentração de potássio obtida pelo Modelo Proposto considerando a
Malha Adaptada (
ε
01,0
=
e p = 3), e
como variáveis independentes, os perfis de concentração
de potássio obtidos pelo Método de Miranda et al. (2005) (Observado experimentalmente e
simulado (MIDI)) (Figura 39).
A partir da Figura 39, é elaborada a Tabela 7 usada na construção das Curvas de
Regressão Linear, considerando a Concentração de Potássio.
Tabela 7 – Profundidade versus Concentração do Modelo Proposto e do Método de Miranda et al. (Observada
experimentalmente e MIDI)
Profundidade
-2 -4 -8 -10 -12 -14 -18 -20 -24 -30- -40 -60 -70
Miranda Obs:
Obs
C
460 400 280 260 240 180 100 40 40 20 20 30 50
Miranda MIDI :
MIDI
C
480 400 320 260 200 122 80 40 40 20 50 50 50
Adapt:
adapt
C
480 400 300 140 80 20 10 5 5 2 50 50 50
Universidade de São Paulo – Escola de Engenharia de São Carlos
163
6.3.4.1 Regressão Linear com as concentrações do Modelo Proposto e do Método de
Miranda et al. (experimental)
A Curva de Regressão Linear, tendo a Concentração de Potássio do Modelo Proposto,
considerando a Malha Adaptada (
ε
01,0
=
e p = 3), como variável dependente e a
Concentração de Potássio obtida pelo Método de Miranda et al. (2005) (Observada
experimentalmente), como variável independente, é dada pela Figura 45. Como o Coeficiente
de Determinação é R
2
= 0,8220, próximo de 1, significa que a concentração de potássio da
Malha Adaptada é bem explicada pela Concentração de Potássio de Miranda et al. (2005)
(experimental). A equação da Curva de Regressão Linear é
adapt
C
= -3,6781 + 0,9697
Obs
C
(6.6)
Figura 45 – Curva de Regressão Linear da concentração de potássio do Modelo Proposto ( Malha Adaptada)
com a concentração de potássio obtida pelo Método de Miranda et al. (Observada experimentalmente)
6.3.4.2 Regressão Linear das concentrações do Modelo Proposto e do Método de
Miranda et al. (MIDI)
A Curva de Regressão Linear, tendo a Concentração de Potássio do Modelo
Proposto, considerando a Malha Adaptada (
ε
01,0
=
e p = 3) como variável dependente e a
Concentração de Potássio obtida pelo Método de Miranda et al. (2005) (MIDI), como variável
164
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
independente, é dada pela Figura 46. O Coeficiente de Determinação é R
2
= 0,9172 e a
equação da Curva de Regressão Linear é
adapt
C
= - 40,9400 + 1,0058
MIDI
C
(6.7)
Como o Coeficiente de Determinação é R
2
= 0,9172, próximo de 1, significa que a
concentração de potássio da Malha Adaptada é bem explicada pela Concentração de Potássio
de Miranda et al. (2005) (MIDI).
Figura 46 – Curva de Regressão Linear da concentração de potássio do Modelo Proposto (Malha Adaptada)
com a concentração de potássio obtida pelo Método de Miranda et al. (MIDI)
Universidade de São Paulo – Escola de Engenharia de São Carlos
165
7 CONCLUSÕES
A simulação realizada pelo modelo proposto neste trabalho é capaz de prever, o
perfil de potencial matricial
ψ
, e, conseqüentemente, a partir da equação (5.4), o perfil de
umidade volumétrica
θ
, bem como a evolução temporal da concentração C do soluto no solo,
em que se constata um bom desempenho do modelo na zona não-saturada do solo,
comparando-se com resultados de Celia, Bouloutas e Zarba (1990), Prasad, Kumar e Sekhar
(2001) e Miranda et al. (2005).
A estratégia de se usar malha adaptada, através do refinamento
h, onde é adotado um
maior nível de refinamento nas regiões onde a solução varia mais fortemente e um menor
refinamento nas regiões em que a solução é suave, fez com que se obtivessem soluções tão
eficientes quanto às com malhas uniformes, diminuindo os graus de liberdade, com um custo
computacional reduzido.
A solução da Equação de Richards, com os dados de Celia, Bouloutas e Zarba (1990),
considerando-se a malha adaptada com
p = 2 e
ε = 0,001 apresentou menor erro em relação à
solução com a malha uniforme fina com 2
6
elementos e p = 1 (Figura 22). A solução obtida
considerando-se a malha adaptada com
p = 2 e
ε = 0,01 é uma boa aproximação em relação à
malha fina (Figura 21), devido à quantidade menor de nós em relação à adaptada com
p = 2 e
ε = 0,001. A simulação do potencial matricial, com ordem de interpolação p = 2 e ε = 0,1
comprovou ser inadequada, possivelmente por ser uma malha demasiadamente “grosseira”, em
comparação com as simulações com ordem de interpolação p = 2 e parâmetros de adaptação
ε = 0,01 e ε = 0,001. Portanto, a escolha do parâmetro de adaptação ε é fundamental na solução do
problema.
A solução da Equação de Richards, obtida pelo Método de Elementos Finitos na
aproximação espacial, com função de interpolação quadrática, evitando a perda de massa do
processo de adaptação de malha implementado, mostrou-se equivalente aos modelos de Celia,
Bouloutas e Zarba (1990) e Prasad, Kumar e Sekhar (2001) (Figura 25), que usam malhas
uniformes e funções-base lineares lagrangeanas. Uma outra estratégia para minimizar a perda
166
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
de massa seria utilizar um valor muito pequeno de
grad
ε
. Desse modo, apenas regiões
extremamente suaves seriam desrefinadas. O impacto dessa estratégia seria uma menor
capacidade de desrefinamento da malha, o que implicaria em maior número de graus de
liberdade e maior esforço computacional.
A solução da Equação de Richards, obtida neste trabalho, mostrou-se também
equivalente ao modelo de Miranda et al. (2005), que a resolve numericamente em um sistema
de volumes finitos, formulada em termos de umidade
θ
. No problema de transporte de
potássio, adotou-se a ordem de interpolação p = 3, sem receio de oscilações numéricas, como
mostram os resultados nas Figuras 34 e 35, devido à suavidade da solução da variável de
estado
θ
C, adotada no problema.
O modelo computacional deste trabalho aplicado à solução da Equação de Richards e
na simulação de fluxo e transporte de potássio, com os dados do trabalho de Miranda et al.
(2005), mostrou a redução, de aproximadamente, 84 % no tempo de processamento, usando-
se malhas adaptadas em relação ao tempo de processamento em malhas uniformes. Já a
diferença de consumo de memória não foi tão significativa em relação às duas malhas
empregadas.
Através do deslocamento do potássio, observa-se que os locais de maiores
concentrações de potássio coincidem com os locais de maiores valores de umidade,
demonstrando seu deslocamento por fluxo de massa, podendo-se ter controle da localização
do potássio no solo em função da irrigação realizada (Figuras 30 e 37).
Para previsão de deslocamento de solutos no solo, os resultados obtidos mostram a
importância da determinação dos parâmetros do solo e de transporte exigidos pelas Equações
de Richards e de Advecção-Dispersão, respectivamente, e das condições iniciais e de
contorno.
Com a solução da Equação de Richards, obtém-se a velocidade da água nos poros
que é necessária para a resolução da Equação de Transporte de Solutos que prevê o
deslocamento destes no solo, contribuindo para a prevenção e realização de previsões de
poluição e contaminação do meio ambiente. Em particular, fornece subsídios para projetos de
aterros sanitários, caso disponha de parâmetros de transporte de chorume.
Com as ferramentas da Estatística, usando-se critérios comparativos, foram
apresentados os resultados obtidos no Modelo Proposto, tanto considerando a Malha
Uniforme Fina quanto a Malha Adaptada, em relação às soluções da Literatura. Os
coeficientes de determinação próximos de 1 (um), calculados nesta seção, indicam uma boa
Universidade de São Paulo – Escola de Engenharia de São Carlos
167
explicação entre essas variáveis. As equações de regressão e gráficos fornecem uma boa
estimativa na comparação desses dados.
7.1 Sugestões para pesquisas futuras
Outras aplicações da metodologia estudada neste trabalho poderão ser efetuadas
usando-se as duas outras formulações da Equação de Richards, a formulação em função de
umidade e a mista, efetuando-se comparações entre as soluções e estudando a influência de
cada solução na Equação de Transporte de Solutos. Todos esses aspectos poderão ser ainda
estudados, considerando-se a adaptatividade
hp. Poderá ser feita também a implementação do
código PZ, considerando-se as Equações de Richards e de Advecção-Dispersão nas
dimensões 2 e 3.
168
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
Universidade de São Paulo – Escola de Engenharia de São Carlos
169
8 REFERÊNCIAS BIBLIOGRÁFICAS
ABDOU, H. M.; FLURY, M. Simulation of water flowand solute transport in free-drainage
lysimeters and field soils with heterogeneous structures
.
European Journal of Soil Science
,
Pullman, n. 55, p. 229–241, June 2004.
AHMED, S. et al. Two dimensional leachate estimation through landfills.
Journal of
Hydraulic
Engineering
, Reston, v. 118, n. 2, p. 306-322, Fevereiro 1992.
ALLEN, M. B.; MURPHY, C. L. A finite element collocation method for variably saturated
flow in two space dimensions.
Water Resources Research
, v. 22, p. 1537-1542, 1986.
AL-YOUSFI, B. A.; POHLAND, F. G.; VASUKI, N.C. Design of landfill leachate
recirculation systems based on flow characteristic. In: PURDUE UNIVERSITY
INDUSTRIAL WASTE CONFERENCE, 47.; 1992, Indiana.
Proceedings…
Indiana: Purdue
University, West Lafayette, Maio 1992. p. 91-100.
ANDERSON, M. P. Using models in simulate the movement of contaminants through
groundwater systems.
Critical Rev. in Envir. Control
, New Jersey, v. 9, n. 2, p. 97-156,
1979.
ANDERSON, M. P.; WOESSNER, W. W.
Applied groundwater modeling:
Simulation of
flow and advective transport. San Diego: Academic Press, 1992. 381 p.
BEAR, J. Some experiments on dispersion.
Journal of Geophysical Research,
v. 66, n. 8, p.
2455-2467, 1961.
______. Dynamics of fluids in porous media
. New York: Elsevier, 1972. 764 p.
______
.
Hydraulics of groundwater
. New York: Dover Publications, INC. Mineola, 2007.
569 p.
170
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
BEAR, J.; VERRUIJIT, A.
Modeling groundwater flow and pollution
. Holland: D. Reidel
Dordrecht, 1993. 414 p.
BEDIENT, P. B.; RIFAI, H. S.; NEWELL, C. J.
Ground Water Contamination
: Transport
and Remediation. New Jersey: Prentice-Hall, 1994. 560 p.
BIDONE, F. R. A.; POVINELLI, J.
Conceitos Básicos de Resíduos Sólidos
. São Carlos:
Escola de Engenharia de São Carlos, Universidade de São Paulo, 1999. 109 p.
BIDWELL, V. J. State–space mixing cell model of unsteady solute transport in unsaturated
soil.
Environmental Modelling and Software
, Amsterdan, v. 14, n. 2, p. 161-169. 1999.
BIGGAR, J. W.; NIELSEN, D. R. Miscible displacement in soils: II Behavior of tracers.
Soil
Science Society of America Proceedings
, Pittsburg, v. 26, p. 125-128, 1962 a.
BIGGAR, J. W.; NIELSEN, D. R. A some comments on molecular diffusion and
hydrodynamic dispersion in porous media.
Journal of Geophysical Research
, Columbia, v.
67, p. 3636-3637, 1962 b.
BORDIER, C.; RATHLE, J.; ZIMMER, D. Hydraulic functioning and clogging diagnosis of
leachate collection system. In: International Landfill Symposium, 6.;
1997,
Cagliari.
Proceedings Sardinia 97…
Cagliari: CISA publisher, 1997. v. 3. p. 361-372.
BOUWER, H.
Groundwater Hydrology
. New York: MacGraw-Hill, 1978. 480 p.
BOU-ZEID, E.; EL-FADEL, M. Parametric sensitivity analysis of leachate transport
simulations at landfills.
Waste Management
, New York, n. 24, p. 681–689, 2004.
BREBBIA, C. A.; CONNOR, J. J.
Finite Elements for Fluid Flow
. London: Butterworths,
1976. 124 p.
BUCZKO, U.; HOPP, L.; BERGER, W.; DURNER, W.; PEIFFER, S.; SCHEITHAUER, M.
Simulation of chromium transport in the unsaturated zone for predicting contaminant entries
into the groundwater.
J. Plant Nutr. Soil Sci
, Weinheim, n. 167, p. 284–292, 2004.
BUNSRI, T.; SIVAKUMAR, M.; HAGARE, D. Numerical Modelling of Tracer Transport in
Unsaturated Porous Media.
Journal of Applied Fluid Mechanics
, v. 1, n. 1, p. 62-70, 2008
a.
Universidade de São Paulo – Escola de Engenharia de São Carlos
171
BUNSRI, T.; SIVAKUMAR, M.; HAGARE, D. Influence of Dispersion on Transport of
Tracer through Unsaturated Porous Media.
Journal of Applied Fluid Mechanics
, v. 1, n. 2,
p. 37-44, 2008 b.
CAMPOS, T. M.; ELBACHÁ, A. T. Avaliação do fator de retardamento por adsorção no
transporte de zinco em solos argilosos. In: SIMPÓSIO SOBRE BARRAGENS DE
REJEITOS E DISPOSIÇÃO DE RESÍDUOS, 1991, Rio de Janeiro.
Anais
… Rio de Janeiro:
REGEO, 1991. p. 271-282.
CAPUTO, J. G.; STEPANYANTS, Y. A. Front Solutions of Richards’ Equation.
Transport
in porous media
, Netherlands
,
v. 74, n. 1, p. 1-20, 2008.
CAREY, P. L.; BIDWELL, V. J.; MCLAREN, R. G.
Chromium(VI) leaching from large
undisturbed soil lysimeters following application of a simulated copper-chromium-arsenic
(CCA) timber preservative.
Australian Journal of Soil Research
, Australia, v. 40, n. 2, p.
351–365, 2002.
CARSEL, R. F.; SMITH, C. N.; MULKEY, L. A.; DEAN, L. A.; JOWISE, D. User’s manual
for the pesticide root zone model (PRZM) release1.
U. S. Envirnmental Agency
, EPA,
Athens, Georgia, USA, v. 600, n. 3, p. 84-109, 1984.
CELIA, M. A.; BOULOUTAS, E. T.; ZARBA, R. L. A general mass conservative numerical
solution for the unsaturated flow equation.
Water Resources Research,
v. 26, p. 1483-1496,
1990.
CELIA, M. A.; AHUJA, L. R.; PINDER, G. F. Orthogonal collocation and alternating-
direction procedures for unsaturated flow problems.
Adv. Water Resources,
v.
10, p. 178-
187, 1987.
CHAPRA, S. C.
Surface Water-Quality Modeling
. Texas: McGraw-Hill, 1997. 844 p.
CHEN, P. H.; WANG, C. Y. Investigation into municipal waste leachate in the
unsaturated
zone of red soil.
Environment International
, Amsterdam, v. 23, n. 2, p. 237-245, 1997.
CHEN, Y. H. Finite element analysis of contaminant transport in groundwater
.
Applied
Mathematics and Computation
, California, v. 127, Issue 1, p. 23-43, 25 mar. 2002.
COOLEY, R. L. Some new procedures for numerical solution of variably saturated flow
problems.
Water Resources Research
, v. 19, p. 1271-1285, 1983.
172
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
CORRÊA, M. M.
Desenvolvimento e teste de modelo de transporte unidimensional de
solutos no solo
. 2001. 104 p. Tese (Doutorado em Solos e Nutrição de Plantas) –
Universidade Federal de Viçosa, Viçosa, 2001.
COSTA, C. T.; CASTRO, M. A. H. Uma metodologia numérico–analítica aplicada a
problemas transparentes de transporte de contaminante.
Revista Brasileira de Recursos
Hídricos
, Porto Alegre, v. 12, n. 4, p. 183-197, 2007.
COSTA, S. N.
Desenvolvimento de um modelo numérico para simulação do transporte
de água e solutos no solo, sob condições de escoamento não-permanente na vertical
.
1998. 153 p. Tese (Doutorado em Engenharia Agrícola) – Universidade Federal de Viçosa,
Viçosa, 1998.
COSTA, S. N.; MARTINEZ, M. A.; MARTINS, J. H.; FERREIRA, P. A. SIMASS-Modelo
para simular o Transporte de água e solutos no solo I:
Desenvolvimento e teste de
sensibilidade
. Revista Brasileira de Engenharia Agrícola e Ambiental
, Campina Grande, v.
3, n. 2, p. 183-189, 1999.
CUMINATO, J. A.; MENEGUETE JUNIOR, M.
Discretização de Equações Diferenciais
Parciais:
Técnicas de Diferenças Finitas. São Carlos: ICMC/USP, 1999. 203 p.
DARCY, H.
Les Fountaines Publiques de la Ville de Dijon - Appendice D
. Paris:
Dalmont, 1856. 647 p.
DAUS, A. D.; FRIND, E. O.; SUDICKY, E. A. Comparative error analisys in finite element
formulation of advection-dispersion equation.
Advance in Water Resources
, v. 8, p. 86-95,
1985.
DE JOSSELIN DE JONG, G. Longitudinal and transverse diffusion in granular deposits.
Transactions, American Geophysical Union,
Washington, n. 39, p. 67–74, 1958.
DEMETRACOPOULOS, A. C.; KORFIATIS, G. P.; BOURODIMOS, E. L.; NAWY, E. G.
Unsaturated flow through solid waste landfills: model and sensitivity analysis.
Water
Resources Bulletin
, California, v. 22, n. 4, p. 601-609, 1986.
DEMETRACOPOULOS, A. C
. Overview of landfill bottom liner hydraulics.
Water
Resources Bulletin
, California, v. 24, p. 49–56, 1988.
Universidade de São Paulo – Escola de Engenharia de São Carlos
173
DEVLOO, P. R. B. PZ: An object oriented environment for scientific programming.
Computer Methods in applied mechanics and Engineering
, Amsterdam, v. 150, p. 133 -
153, 1997.
ELBACHÁ, A. T.
Estudo da Influência de Alguns Parâmetros no Transporte de Massa
em Solos Argilosos
. 1989. 178 p. Dissertação (Mestrado em Engenharia Civil) -
Departamento de Engenharia Civil, PUC-RIO, Rio de Janeiro, 1989.
FARTHING, M. W.; KEES, C. E.; RUSSEL, T. F.; MILLER, C. T. An ELLAM
approximation for advective-dispersive transport with nonlinear sorption.
Advances in Water
Resources
, v. 29, p. 657-675, 2006.
FATTA, D.; NAOUM, D.; LOIZIDOU, M. Integrated environmental monitoring and
simulation system for use as a management decision support tool in urban areas.
Journal of
Environmental Management
, Oxford, n. 64, p. 333-343, 2002.
FETTER, C. W.
Contaminant Hydrogeology
. New Jersey: Prentice Hall, 1992. 317 p.
FETTER, C. W.
Contaminant Hydrogeology
. New York: Macmillan Publishing Company,
1993. 458 p.
FREDLUND, D. G.; MORGENSTERN, N. R. Stress state variables for unsaturated soils.
Journal of Geotechnical Engineering Division
, ASCE, New York, v. 103, n. 5, p. 447-466,
1977.
FREDLUND, D. G.; RAHARDJO, H.
Soil Mechanics for Unsaturated Soils.
New York:
John Wiley & Sons, 1993. 517 p.
FREEZE, R. A. Three-dimensional transient, saturated-unsaturated flow in a groundwater
basin.
Water Resources Research
, v. 7, n. 2, p. 347-366, 1971.
FREEZE, R. A.; CHERRY, J. A.
Groundwater
. New Jersey: Prentice-Hall, 1979. 604 p.
GILES, R. V.
Mecânica dos Fluidos e Hidráulica
.
São Paulo: McGraw-Hill, 1976. 412 p.
GILLHAM, R. W.; CHERRY, J. A. Contaminant Migration in Saturated Unconsolidated
Geologic Deposits.
Spec. Pap. Geol. Soc. Am.,
Virginia,
n.
189, p. 31-62, 1982.
174
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
GIRAULT, V.; RAVIAT, P. A.
Finite Element Methods for Navier-Stokes Equations –
Theoria and algorithms, Springer-Verlag, 1986.
GOMES, C. C.
Método de elementos analíticos para modelagem matemática de fluxo
hídrico subterrâneo regional
. 2000. Tese (Doutorado em Engenharia Civil) Departamento
de Engenharia Hidráulica e Ambiental. Universidade Federal do Ceará, Fortaleza, 2000.
HAMON, W. R. Estimating potential evapotranspiration.
J Hydraul Div Proc Am Soc Civil
Eng,
Oxford,
n. 87, p.107–120, 1961.
HANCOCK, T. C. et al. Pesticide Fate and Transport throughout Unsaturated Zones in Five
Agricultural Settings.
J Environ Qual,
Pittsburgh, v. 37, p. 1086-1100, 2008.
HAVERKAMP, R.; VAUCLIN, M.; TOUMA, J.; WIERENGA, P.; VACHAUD, G.
Comparison of numerical simulation models for one-dimensional infiltration.
Soil Sci. Soc.
Am. J.,
Pittsburgh, v. 41, p. 285-294, 1977.
HAVERKAMP, R.; VAUCLIN, M. A note on estimating finite difference interblock
hydraulic conductivity values for transient unsaturated flow problems.
Water Resources
Research
, v. 15, p. 181-187, 1979.
HAYDAR, M. M.; KHIRE, M. V.
Leachate recirculation using horizontal trenches in
bioreactor landfills
.
Journal of Geotechnical and Geoenvironmental Engineering,
ASCE,
New York, p. 837-847, 2005.
HAYHOE, H. N. Study of the relative efficiency of finite difference and Galerkin tecniques
for modeling soil-water transfer.
Water Resour. Res.,v.
14, p. 97-102, 1978.
HE, X.; REN, L.
An adaptive multiscale finite element method for unsaturated flow problems
in heterogeneous porous media.
Journal of Hydrology
, v. 374, p. 56-70, 2009.
HEIDMAN, J. A.; BRUNNER, D. R.
Solid waste and water quality.
Water Pollution
Control Fed
., California, v. 45, n. 6, p. 1198-1202, Junho 1973.
HILLEL, D.
Fundamentals of soil physics
. New York: Academic, 1980. 413p.
HINDMARSH, A. C.; GRESHO, P. M.; GRIFFITHS, D. F. The stability of explicit Euler
time integration for certain finite difference approximations of the multidimensional
Universidade de São Paulo – Escola de Engenharia de São Carlos
175
advection-diffusion equation
. Int. Journal for numerical methods in fluids
, Livermore, v.
4, p. 853-897, 1984.
HORNBECK, R. W.
Numerical Methods
. Pittsburgh: Carnegie Mellon Univ, 1975. 326p.
HORNUNG, U.; MESSING, W. Truncation errors in the numerical solution of horizontal
diffusion in saturated/unsaturated media.
Adv. Water Resour
., v. 6, p. 165-168, 1983.
http://pt.wikipedia.org/wiki/MATLAB. Acesso em: 13 jun. 2009.
HUYAKORN, P. S.; LESTER, B. H.; MERCER, J. W. An efficient finite element technique
for modeling transport in fractured media, 1: Single-species transport.
Water Resour. Res.
, v.
19, n. 3, p. 841-854, 1983.
HUYAKORN, P. S.; PINDER, G. F.
Computational Methods in Subsurface Flow
. New
York: Academic Press, 1983. 473p.
INSTITUTO BRASILEIRO DE GEOGRAFIA E ESTATÍSTICA.
Pesquisa Nacional de
Saneamento
. Disponível em:
<http://www.ibge.gov.br/home/presidencia/notícias/27032002>. Acesso em: 17 out. 2006.
ISTOK, J. D.
Groundwater Modeling by the Finite Element Method
. Washington:
American Geophysical Union, 1989. 495 p.
JAVADI, A. A.; Al-NAJJAR, M. M.; EVANS, B. Finite element modelling of contaminant
transport through soils – Case study.
Journal of Geotechnical and Geoenvironmental
Engineering
, ASCE, New York, v. 134, n. 2, p. 214-230, 2008.
KAW, A.; KALU, E. E.
Numerical Methods with Applications
. Florida: University of
South Florida, 2008. 728 p.
KHANBIVARDI, R. M.; AHMED, S.; GLEASON, P. J. Flow Investigation for Landfill
Leachate (FILL).
Journal of Environmental Engineering
, ASCE, New York, v. 121, n. 1, p.
45-57, Janeiro 1995.
KIM, G.; KOO, J.; SHIN, J.; SHON, J.; LEE, S. A study of methods to reduce groundwater
contamination around a landfill in korea.
Journal of Environmental Hydrology
, Texas, v. 7,
n. 10, p. 1-9, 1999.
176
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
KLEIN, M. Evaluation and Comparison of Pesticide Leaching models for Registration
Purposes, Results of Simulations performed with the Pesticide Leaching Model.
Journal of
Environmental Science & Health
, New York, v. 29, n. 6, p.1197-1209, 1994.
KLEIN, M.
PELMO
: Pesticide Leaching Model, User manual version 2.0. Alemanha:
Fraunhofer-Institut für Umweltchemie und Ökotoxikogie, 1995. 103 p.
KLEIN, M.; MÜLLER, M.; DUST, M.; GÖRLITZ, G.; GOTTESBFIREN, B.; HASSINK, J.;
KLOSKOWSKI, R.; KUBIAK , R.; RESSELER, H.; SCHILLER, H.; STEIN, B.;
VEREEEKEN, H. Validation of the pesticide leaching model pelmo using lysimeter studies
performed for registration.
Cheraosphere
, Great Britain, v. 35, n. 11, p. 2563-2587, 1997.
KORFIATIS, G. P.; DEMETRACOPOULOS, A. C.; BOURODIMOS, E. L.; NAWY, E.
Moisture transport in a solid waste column.
Journal of Environmental Engineering
, ASCE,
New York, v. 110, n. 4, p. 789-796, Agosto 1984.
KUTILEK, M.; NIELSEN, D. R.
Soil hydrology
. Berlin: Verlag, 1994. 370 p.
LIGGETT, J. A.; LIU, P. L. F.
The Boundary Integral Equation Method for Porous
Media Flow
. London: Allen & Unwin, 1983. 225 p.
LIMA, E. N.
Caracterização e estudo do comportamento térmico de chorume, de
composto maturado e derivados
. 2008. 149 f. Tese (Doutorado em Química) - Instituto de
Química de Araraquara, Universidade Estadual Paulista “Júlio de Mesquita Filho”,
Araraquara, 2008.
LIPPMAN, S. B.; LAJOIE, J.
++
C
Primer
. 3rd ed. Massachusetts: Addison Wesley, 1998.
LU, C. A model of leaching behaviour from msw incinerator residue landfills
.
Waste
Management & Research
, Northampton, n. 14, p. 51–70, 1996.
MANOEL FILHO, J. Água Subterrânea: Histórico e Importância. In: FEITOSA, F. A. C.;
MANOEL FILHO, J.
Hidrogeologia
. Fortaleza: CPRM/REFO-LABHID/UFPE, 2000, 391p.
McCREANOR, P. T.; REINHART, D. R.
Mathematical modeling of leachate routing in a
leachate recirculating landfill.
Wat. Res. Elsevier Science Ltd.
, Great Britain, v. 34, n. 4, p.
1285-1295, 2000.
Universidade de São Paulo – Escola de Engenharia de São Carlos
177
McDONALD, M. G.; HARBAUGH, A. W.
MODFLOW
: A modular three-dimensional
finite-difference groundwater flow model
. Virginia: U. S. Geological Survey, 1988.
MILLER, C. T.; ABHISHEK, C.; FARTHING, M. W. A spatially and temporally adaptive
solution of Richards’equation.
Advances in Water Resources,
v. 29, p. 525–545, 2006.
MILLY, P. C. D. Advances in modeling of water in the unsaturated zone.
Transp. Porous
media
, Boston, n. 3, p. 491-514, 1988.
MIRANDA, J. H.
Modelo para simulação da dinâmica de nitrato em colunas verticais de
solo não-saturado.
2001. 79 p. Tese (Doutorado em Irrigação e Drenagem) – Escola Superior
de Agricultura “Luiz de Queiroz”, Universidade de São Paulo, Piracicaba, 2001.
MIRANDA, J. H.; DUARTE, S. N.; LIBARDI, P. L. (2004). MIDI: Modelo para simulação
do deslocamento de solutos em colunas verticais de solo não-saturado. In: SIMPÓSIO
BRASILEIRO DE SOLOS NÃO-SATURADOS, 5., 2004, São Carlos.
Anais
... São Carlos:
USP, 2004. p. 155-160.
MIRANDA, J. H.; DUARTE, S. N.; LIBARDI, P. L.; FOLEGATTI, M. V. Simulação do
deslocamento de potássio em colunas verticais de solo não-saturado.
Eng. Agrícola,
Jaboticabal, v. 25, n. 3, p. 677-685, 2005.
MITCHELL, J. K. Conduction phenomena: from theory to geotechnical pratice.
Géotechnique
, Canada, v. 41, n. 3, p. 299-340, 1991.
MOLER, C. John Murphy meets scientific computing's leading'celebrity', Cleve Moler, co-
founder of The MathWorks
. Scientific Computing World
, Cambridge,
2006.
MORETTIN, Pedro A.; BUSSAB, Wilton O.; Estatística Básica, São Paulo: Editora Atlas,
2004
MORTON, K. W.; MAYERS, D. F.
Numerical Solution of Partial Differential Equations,
An Introduction.
Cambridge: Cambridge University Press, 2005. 278 p.
MUALEM, Y. A new model for predicting the hydraulic conductivity of unsatured porous
media
.
Water Resources Research
, v. 12, n. 3, p. 513-522, 1976.
MUSKAT, M.
Flow of homogeneous fluids
. New York: McGraw-Hill Book, 1937. 468 p.
178
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
NAGRA (1994). Bericht zur Langzeitsicherheit des Endlagers SMA am Standort Wellenberg.
NAGRA Technischer bericht NTB 94-106. NAGRA, Wettingen, Switzerland.
NAIME, J. M.
Um novo método para estudos dinâmicos, in situ, da infiltração da água
na região não-saturada do solo.
2001. 146 p. Tese (Doutorado em Ciências da Engenharia
Ambiental) – Escola de Engenharia de São Carlos, Universidade de São Paulo, São Carlos,
2001.
NARASIMHAN, T. N.; WITHERSPOON, P. A. An integrated finite difference method for
analyzing fluid flow in porous media
.
Water Resour. Res
., v. 12, p. 57–64, 1976.
NEUMAN, S. P. Saturated-unsaturated seepage by finite elements.
J. Hydraul. Div. Am.
Soc. Civ. Eng.
, Oxford, n. 99, p. 2233-2250, 1973.
NIELSEN, D. R.; BIGGAR, J. W. Miscible displacement in soils: I. Experimental
Information.
Soil Science Society of America,
Pittsburgh, v. 25, n. 1, p. 1-5, 1961.
NIELSEN, D. R.; van GENUCHTEN, M. T.; BIGGAR, J. W. Water flow and solute transport
processes in the unsaturated zone.
Water Resour. Res
., v. 22, p. 89-108, 1986.
NEWELL, C. J.; MCLEOD, R. K.; GONZALEZ, J. R. BIOSCREEN User’s Manual.
National Attenuation Decision Support System. Version 1.3.
National Risk Management
Research Laboratory
, EPA/600/R-96/087, August 1996.
NOBLE, J. J.; ARNOLD, A. E. Experimental and mathematical modeling of moisture
transport in landfills
.
Chemical Engineering Communication
, London, v. 100, p. 95-111,
1991.
OLIVEIRA, L. F. C.
Modelo para transporte de solutos no solo e no escoamento
superficial
. 1999. 171 p. Tese (Doutorado em Engenharia Agrícola) – Universidade Federal
de Viçosa, Viçosa, 1999.
OLSSON, K. G.; HEYDEN, S.
Introduction to the finite element method problems.
Byggnadsmekanik, Lund: Ottosen & Petersson, Prentice Hall, 2001. 149 p.
PADOIN, E. L.; BORGES, P. A. P.; DILL, S. L.; BORGES, R. S. (2006). Análise do
deslocamento da água em solos saturados e paralelização de seus métodos computacionais
. In:
CONGRESSO DE INICIAÇÃO CIENTÍFICA E TECNOLÓGICA EM ENGENHARIA,
FEIRA DE PROTÓTIPOS, 21., 2006, cidade.
Resumos…
2006.
Universidade de São Paulo – Escola de Engenharia de São Carlos
179
PARKHURST, D. I.; THORSTENSON, C. C.; PLUMMER, L. N.
PHREEQE
: A computer
program for geochemical calculations. Reston: E. U. Geol. Survey, 1990. 195 p.
PERKINS, T. K.; JOHNSTON, O. C. A Review of Difusion and Dispersion in Porous Media.
Society of Petroleum Engineering Journal
, London, v. 3, n. 1, p. 70-84, 1963.
PHILIP, J. R. The theory of infiltration. 1. The infiltration equation and its solution.
Soil
Science
, Pittsburg, v. 83, p. 345-357, 1957.
PINCHOVER, Y.; RUBINSTEIN, J.
An Introduction to Partial Differential Equations.
Cambridge: Cambridge University Press, 2005. 400 p.
PINDER, G. F.; GRAY, W. G.
Finite Element Simulation in Surface and Subsurface
Hydrology
. New York: Academic Press, 1977. 295 p.
POLYANIN, A. D; ZAITSEV, V. F.
Handbook of Nonlinear Partial Differential
Equations
. Boca Raton- London: Chapman & Hall/CRC Press, 2004. 814 p.
PORTELA, A.; CHARAFI, A.
Finite Elements Using Maple:
A Symbolic Programming
Approach. Berlin: Springer, 2002. 325 p.
PRASAD, K. S. H.; KUMAR, M. S. M.; SEKHAR, M. Modelling flow through unsaturated
zones: Sensitivity to unsaturated soil properties
.
Sãndhanã
, v. 26, n. 6, p. 517-528, 2001.
PREVEDELLO, C. L.
Física do solo com problemas resolvidos
. Curitiba: Salesward-
Discovery, 1996. 446 p.
PREVEDELLO, C. L.; LOYOLA, J. M. T.; COSTABILE, M. S.; HORODENSKI, J. Solução
Numérica para o Processo da Infiltração da Água no Solo.
Scientia Agraria
, Paraná, v. 3, n.
1-2, p. 29-39, 2002.
PRESS, W. H.
Numerical recipes:
the art of scientific computing.
Cambridge: Cambridge
University Press, 2007. 1235p
QUEIROZ, P. I. B. Sobre a modelagem de transporte de contaminantes no solo. In:
CONGRESSO BRASILEIRO E GEOTECNIA AMBIENTAL, 4., 1999, São José dos
Campos.
Anais
…São José dos Campos: REGEO’99, 1999. p. 307-314.
180
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
REDDY, J. N.
Applied Functional Analysis and Variational Methods in Engineering
.
New York: McGraw-Hill, 1986. 546 p.
REICHARDT, K. Dinâmica da material e da energia em ecossistemas.
Escola Superior de
Agricultura “Luiz de Queiroz”, Universidade de São Paulo
, Piracicaba, p. 21-25, 1996.
REICHARDT, K.; TIMM, L. C.
Solo, Planta e Atmosfera
– Conceitos, Processos e
Aplicações. Barueri- SP: Editora Manole Ltda, 2004. 478p.
RICHARDS, L. A. Capillary conduction of liquids through porous medium.
Physics,
New
York
,
v
.
1, p. 318-333, 1931.
RICHTMYER, R. D.; MORTON, K. W.
Difference methods for initial value problems.
2.
ed. London: Wiley-Interscience. 1967. 405 p.
RITTERLING,
J. M.; STANSBURY, J. S. Alternative use of MULTIMED model for subtitle
D landfill applications
.
Practice Periodical of Hazardous, Toxic, and Radioactive Waste
Management,
Kansas, 1998.
ROBINSON, R. A.; STOKES, R. H. 1965.
Eletrolyte solutions.
2nd. ed. London:
Butterworths Scientific Publications, 1965. 571 p.
ROCHA, F. A.; MARTINEZ, M. A.; ANTÔNIO T. MATOS, A. T.; REINALDO B.
CANTARUTTI, R. B.; SILVA, J. O. DA. Modelo numérico do transporte de nitrogênio no
solo. Parte I: Desenvolvimento e teste do modelo
.
Revista Brasileira de Engenharia
Agrícola e Ambiental
, Campina Grande,
v. 12, n. 1, p. 47–53, 2008.
ROLNIK, V. P.
Obtenção Automatizada da Equação Modificada
.1998. Dissertação
(Mestrado em Ciências de Computação) - Departamento de Ciências de Computação e
Estatística, Instituto de Ciências Matemáticas e de Computação, Universidade de São Paulo,
São Carlos, 1998.
ROSS, P. J., BRISTOW, K. L. Simulating Water Movement in Layered and Gradational Soils
sing the Kirchhoff Transform.
Soil Sci. Soc. Am. J
., Pittsburg, v. 54, p. 1519-1524, 1990.
ROSSI, P.; MIRANDA, J. H.; DUARTE, S. N. Curvas de distribuição de efluentes do íon
nitrato em amostras de solo deformadas e indeformadas.
Eng. Agríc
., Jaboticabal, v. 27, n. 3,
p. 675-682, set./dez. 2007.
Universidade de São Paulo – Escola de Engenharia de São Carlos
181
ROWE, R. K. 1987. Pollutant transport through barrier. Geotechnical Practice for waste
Disposal.
Special Geotechnical Publication,
ASCE, New York, n. 13, p.159-181, 1987.
ROTH, K.
Lecture notes in soil physics
. Germany: Institute of Environmental Physics,
University of Heidelberg, 2006. 21 p.
RÜBENKÖNIG, O. The Finite Difference Method (FDM) - An introduction.
Albert
Ludwigs University of Freiburg
, Germany, p.1-5, 2006.
RUBIN, J. Transport of reacting solutes in porous media: Relation between mathematical
nature of problem formulation and chemical nature of reactions.
Water Resour. Res
., v. 19,
p. 1231-1252, 1983.
RUNCHAL, A. K.; SAGAR, B. PORFLOW, a Model for Fluid Flow, Heat and Mass
Transport in Multifluid, Multiphase Fractured or Porous Media.
Analytic and
Computational Research Inc.
, West Los Angeles, CA, 1998
SAFFMAN, P. G. A theory of dispersion in a porous medium.
Journal of Fluid Mechanics,
Cambridge, v. 6, n. 3, p. 321–349, 1959.
SALHOTRA, A. M.; MINEART, P.; SHARP-HANSEN, S; ALLISON, T. Multimedia
exposure assessment model (MULTIMED) for evaluating the land disposal of wastes: Model
theory
.
U. S. Envir. Protection Agency
, Athens, Ga, 1990.
SÃO PAULO (Estado). Secretaria do Meio Ambiente.
Companhia Ambiental do Estado de
São Paulo.
Disponível em: <http://www.cetesb.sp.gov.br/>. Acesso em: 25 mar. 2004.
Acesso em: 11 abr. 2006.
SCHALCH, V.
Produção e características do chorume em processo de decomposição de
lixo urbano
. 1984. 103 f. Dissertação (Mestrado em Hidráulica e Saneamento) – Escola de
Engenharia de São Carlos, Universidade de São Paulo, São Carlos, 1984.
SCHEIDEGGER, A.
Physics of Flow through Porous Media
. Toronto, Canada: University
of Toronto Press, 1963. 313 p.
SCHROEDER, P. R; LIOYD, C. M; ZAPPI, P. A; AZIZ, N. M.
Hydrologic evaluation of
landfill performance (help) model:
user's guide for version 3. Cincinnati: EPA, 1994. 94 p.
182
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
SCHNOOR, J. L.
Environmental Modeling:
Fate and Transport of Pollutants in Water, Air
and Soil. New York: John Wiley & Sons, Inc. 1992. 436 p.
SEILER, K. P.; GAT, J.
Groundwater recharge from run-off, infiltration and percolation
.
Portland: Water Science & Technology Library S., 2007. 241p.
SHARP-HANSEN, S.; TRAVERS, C.; HUMMEL, P.; ALLISON, T.
A subtitle D landfill
application manual for the multimedia exposure assessment model (MULTIMED).
Athens, Ga: U.S. Envier. Protection Agency, 1990. 229 p.
SILVA, A. C.
Tratamento do Percolado de Aterro Sanitário e Avaliação da Toxicidade
do Efluente Bruto e Tratado
. 2002. 79 p. Dissertação (Mestrado em Engenharia Civil) –
COPPE/UFRJ, Rio de Janeiro, 2002.
SIMUNEK, J.; van GENUCHTEN, M. T. The CHAIN-2D Code for Simulating Two-
Dimensional Movement of Water, Heat, and Multiple Solutes in Variably-Saturated Porous
Media, Version 1.1.
Research Report,
U. S. Salinity Laboratory, USDA, ARS, Riverside,
California, n. 136, 194 p, 1994.
SIMUNEK, J.; SEJNA, M.; van GENUCHTEN, M. T.
The HYDRUS-1D software package
for simulating the one-dimensional movement of water, heat, and multiple solutes in
variably-saturated media.
Riverside, California: U.S. Salinity Laboratory, 1998.
SIMUNEK, J.; SEJNA, M.; van GENUCHTEN, M. T.
The HYDRUS-2D software
package for simulating the 2D movement of water, heat, and multiple solutes in variably
saturated media:
Version 2.0. Riverside, California: U.S. Salinity Laboratory, Agriculture
Research Service, USDA, 1999.
SINGH, P.; KANWAR, R. S. Simulating NO
3
-N transport to subsurface drain flows as
affected by tillage under continuous corn using modified RZWQM.
Transactions of the
ASAE
, Cambridge, v. 38, n. 2, p. 499-506, 1995.
SOD, G.A.
Numerical methods in fluid dynamics:
Initial and initial boundary
-
value
problems. Cambridge: Cambridge University Press, 1987.
SPERANDIO, D.; MENDES, J. T.; SILVA, L. H. M.
Cálculo numérico:
Características
Matemáticas e Computacionais dos Métodos Numéricos. São Paulo: Pearson Prentice Hall,
2006. 354p.
SRIVASTAVA, R.; YEH, T. J. Analytical Solutions for One-dimensional, Transient
Infiltration toward the Water Table in Homogeneous and Layered Soils.
Water Resour. Res
.
v. 27, n. 5, p. 753-762, 1991.
Universidade de São Paulo – Escola de Engenharia de São Carlos
183
STONE, L. F.; SILVEIRA, P. M. de; MOREIRA, H. A. A.; BRAZ, A. J. B. P.
Evapotranspiração do feijoeiro irrigado em plantio direto sobre diferentes palhadas de culturas
de cobertura
.
Pesquisa Agropecuária Brasileira
, Brasília, v. 41, n. 4, 2006.
STRAUB, W. A.; LYNCH, D. R. Models of landfill leaching:moisture movement and
inorganic strength.
Journal of Environmental Engineering Division
, ASCE, New York, v.
108, p. 231–250, 1982.
TORIDE, N.; LEIJ, F. J.; van GENUCHTEN, M. T. The CXTFIT Code for Estimating
Transport Parameters from Laboratory or Field Tracer Experiments, Version 2.1.
Research
Report
, U.S. Salinity Laboratory, Riverside, California, n. 137. 119 p. 1999.
TSAI, W. S.; LEE, T. H.; CHEN, C. J.; LIANG, S. J.; KUO, C. C.
Finite analytic model for
flow and transport in unsaturated zone.
Journal of Engineering Mechanics,
Reston
,
v. 126,
n. 5,
p. 470-479, Maio 2000.
TUCCI, C. E. M.
Modelos Hidrológicos
. Porto Alegre: Editora da Universidade/UFRGS,
1998. 669 p.
VALOCCHI, A. J. Describing the transport of ion-exchanging contaminants using an
effective Kd approach.
Water Resource Research
, v .20, n. 4, p. 499-503, 1984.
van GENUCHTEN, M. T. A closed-form equation for predicting the hydraulic conductivity
of unsaturated soils.
Soil Sci. Soc. Am. J
., Pittsburgh, n. 44, p. 892-898, 1980.
VASCONCELLOS, C. A. B.; AMORIM, J. C. C. Simulação Numérica da Infiltração da
Água em Meios Porosos Não-saturados Homogêneos. In: SIMPÓSIO BRASILEIRO DE
RECURSOS HÍDRICOS, 14., 2001, Aracaju-SE.
Anais
... Aracaju-SE: ABRH, 2001.
VENKATARAMANI, E. S.; RAO, P. R. M. Darcian, Transitional, and Turbulent Flow
through Porous Media.
Journal of Hydraulic Engineering
, ASCE, New York, v. 124, n. 8,
p. 840-846, 1998.
VERMEULEN, T.; HIESTER, N. K. Ion exchange chromatography of trace elements.
Industrial Engineering Chemistry
, Korea,
n. 44, p. 636–651, 1952.
VOSS, C. I.
SUTRA
, Saturated-Unsaturated Transport, A Finite Element Simulation Model
for Saturated-Unsaturated, Fluid-Density-Dependent Ground-Water Flow with Energy
Transport or Chemically Reactive Single Species Solute Transport. Reston, VA: US
Geological Survey, National Center, 1984.
184
“Simulação de Fluxo de Água e Transporte de Solutos na Zona Não-Saturada do Solo pelo Método de Elementos Finitos Adaptativo”
WANG, H. S.; ANDERSON, M. P. Introduction to groundwater modeling, finite differences
and finite element methods. New York: Academic Press, Inc, 1982. 237 p.
WEBER JR, W. J.; McGINLEY, P. M.; LYNN, E. K. Sorption Phenomena in Subsurface
Systems: Concepts, Models and Effects on Contaminant Fate and Transport.
Water Resourse
Res
., v. 25, n. 5, p. 499-528, 1991.
WENDLAND, E.
Contribuição à simulação de processos em meios porosos
. 2004. 270 p.
Tese (Livre Docente) - Escola de Engenharia de São Carlos, Universidade de São Paulo, São
Carlos, 2004.
WENDLAND, E.; RÜBER, O.
Hydrogeologishe Modelle, Lehrstuhl für Angewandte
Hydrogeologie
. Ruhr: Universität Bochum, 1998. 141 p.
WENDLAND, E.; FLENSBERG, D. A hybrid numerical scheme for multi-phase flow in
heterogeneous porous media . In: HASSANIZADEH, S. M.; SCHOTTING, R. J.; GRAY, W.
G.; PINDER, G. F.
Computational Methods in Water Resources
. 1 ed. Amsterdam:
Elsevier, 2002. v. 1, p. 297-304.
YEH, G
. T.
FEMWATER
: a finite-element model of water flow through saturated-
unsaturated porous media. Oak Ridge, Tenn: Oak Ridge National Laboratory, 1987.
YEH, G. T.; WARD, D. S.
FEMWATER
: a finite-element model of water flow through
saturated-unsaturated porous media
. Rep. ORNL-5567. Oak Ridge, Tenn: Oak Ridge National
Laboratory, 1980.
YONG, R. N. et al.
Principles of contaminant transport in soils. Development in
geotechnical Engineering.
Netherlands: Elsevier Science Publishers, 1992. n. 73, 327 p.
ZARBA, R. L.
A numerical investigation of unsaturated flow
. 1988. M. S. Thesis - Dep.
of Civ. Eng., Mass. Inst. of Technol., Cambridge, 1988.
ZEISS, C.; MAJOR, W. Moisture flow through municipal solid waste: patterns and
characteristics.
Journal of Environmental Systems
, v. 22, n. 3, p. 211-231, 1993.
ZEISS, C.; UGUCCIONI, M. Modified flow parameters for leachate generation.
Water
Environment Research
, v. 69, n. 3, p. 276-285, may/june, 1997.
ZHENG, C.
MT3D
, a modular three-dimensional transport model. Rockville, Maryland,
USA: S. S. Papadopolus & Assoc, 1990. 173 p.
Universidade de São Paulo – Escola de Engenharia de São Carlos
185
ZIENKIEWICZ, O. C.
The Finite Element Method
. New York: McGraw-Hill, 1977. 787 p.
ZIENKIEWICZ, O. C.; MORGAN, K.
Finite Elements and Approximation
. New York:
John Wiley & So, 1983.
Livros Grátis
( http://www.livrosgratis.com.br )
Milhares de Livros para Download:
Baixar livros de Administração
Baixar livros de Agronomia
Baixar livros de Arquitetura
Baixar livros de Artes
Baixar livros de Astronomia
Baixar livros de Biologia Geral
Baixar livros de Ciência da Computação
Baixar livros de Ciência da Informação
Baixar livros de Ciência Política
Baixar livros de Ciências da Saúde
Baixar livros de Comunicação
Baixar livros do Conselho Nacional de Educação - CNE
Baixar livros de Defesa civil
Baixar livros de Direito
Baixar livros de Direitos humanos
Baixar livros de Economia
Baixar livros de Economia Doméstica
Baixar livros de Educação
Baixar livros de Educação - Trânsito
Baixar livros de Educação Física
Baixar livros de Engenharia Aeroespacial
Baixar livros de Farmácia
Baixar livros de Filosofia
Baixar livros de Física
Baixar livros de Geociências
Baixar livros de Geografia
Baixar livros de História
Baixar livros de Línguas
Baixar livros de Literatura
Baixar livros de Literatura de Cordel
Baixar livros de Literatura Infantil
Baixar livros de Matemática
Baixar livros de Medicina
Baixar livros de Medicina Veterinária
Baixar livros de Meio Ambiente
Baixar livros de Meteorologia
Baixar Monografias e TCC
Baixar livros Multidisciplinar
Baixar livros de Música
Baixar livros de Psicologia
Baixar livros de Química
Baixar livros de Saúde Coletiva
Baixar livros de Serviço Social
Baixar livros de Sociologia
Baixar livros de Teologia
Baixar livros de Trabalho
Baixar livros de Turismo
