Modelagem computacional do fluxo unidimensional de água em meio não saturado do solo

Computational model of one-dimensional flow of water in an unsaturated soil

ARTIGOS CIENTÍFICOS

ENGENHARIA DE ÁGUA E SOLO

Edson Wendland^I; Maria de L. P. Pizarro^II

^IProf. Associado, Departamento de Hidráulica e Saneamento, Escola de Engenharia de São Carlos, Universidade de São Paulo, São Carlos - SP, [e-mail]

^IIProfa. Associada, Academia da Força Aérea, Estrada de Aguaí s/n, Pirassununga - SP, [e-mail]

RESUMO

O estudo do fluxo de água em zonas não saturadas do solo é de grande importância para pesquisas relacionadas à disponibilidade hídrica para o desenvolvimento das plantas. Devido ao alto custo, ao tempo demandado e ao esforço humano nas investigações de campo, os modelos matemáticos, aliados às técnicas numéricas e avanços computacionais, constituem-se em uma ferramenta importante na previsão desses estudos. No presente trabalho, objetivou-se solucionar a equação diferencial parcial não linear de Richards mediante a aplicação do Método de Elementos Finitos. Na aproximação espacial, foi empregada a adaptatividade com refinamento "h" na malha de elementos finitos e, na derivada temporal, foi aplicado o esquema de Euler Explícito. A função interpolação polinomial utilizada foi de grau 2, e a que garantiu a conservação de massa da estratégia de adaptação. Para a validação do modelo, foram utilizados dados disponíveis em literatura. A utilização da função interpolação polinomial de grau 2 e o refinamento "h", com considerável redução do tempo de execução da rotina computacional, permitiram uma boa concordância do modelo em comparação a soluções disponíveis na literatura.

Palavras-chave: zona não saturada, equação de Richards, método de elementos finitos.

ABSTRACT

Study of water flow in the unsaturated soil zone is of great importance for research related to the water availability for crop development. Due to the high cost, the time required and the human effort in the field investigations, mathematical models combined with numerical techniques and computational advances are important tools in the prediction of these studies. This work aimed to solve the Richards's non-linear partial differential equation by applying the Finite Element Method. Adaptability with "h" refinement of the finite element mesh was used in the spatial approximation, while Explicit Euler scheme was applied for the time derivative. The polynomial interpolation function used was of degree two, and ensured the mass conservation of the adaptation strategy. To validate the model, data available in the literature were used. Use of the polynomial interpolation function with degree two and the "h" refinement, with considerable reduction of the computational runtime allowed good agreement in comparison to solutions available in the literature.

Keywords: unsaturated zone, Richards's equation, finite element method.

INTRODUÇÃO

O escoamento de fluidos em um meio poroso não saturado é um caso especial de escoamento simultâneo de fluidos não miscíveis (BEAR, 1972). No caso considerado, o fluido em questão é a água, e o fluido não molhante é o ar. O escoamento não saturado ocorre quando a água se movimenta através do perfil de solo cuja umidade se encontra inferior àquela correspondente à saturação, ou seja, quando parte dos espaços vazios (poros) está ocupada pelo ar, o qual é considerado estagnado, isto é, não está em movimento. No ciclo hidrológico, o escoamento em meios porosos não saturados ocorre na fase de infiltração e redistribuição ao longo da chamada zona de aeração.

A primeira experiência que quantificou o fluxo num meio poroso saturado foi publicada por DARCY (1856). A Lei de Darcy permitiu expressar a densidade de fluxo de água através de um meio poroso, por unidade de superfície, em função da condutividade hidráulica do material e do gradiente hidráulico. Essa lei constituiu-se, portanto, na base dos métodos de avaliação quantitativa de recursos hídricos subterrâneos. A Lei de Darcy foi adaptada mais tarde para solos não saturados (REICHARDT & TIMM, 2004; ALFARO SOTO et al., 2007), com a denominação de Lei de Buckingham-Darcy ou Lei de Darcy-Buckingham. Richards formulou as bases teóricas, em 1928, para descrever a percolação de água em um meio poroso não saturado (RICHARDS, 1931). Em consequência, a equação geral do movimento de água em um meio poroso não saturado, isotermicamente e isobaricamente, baseada no Princípio da Conservação da Massa, e na equação de Buckingham-Darcy, recebeu o nome de Equação de Richards. Essa equação pode ser escrita, unidimensionalmente, em três formas equivalentes: a forma mista, a baseada na umidade volumétrica θ do solo e a com base no potencial matricial ψ (PIZARRO, 2009).

Devido à dificuldade em se obter soluções analíticas, a modelagem numérica constitui-se numa tentativa aproximada de representar um fenômeno da natureza (WENDLAND, 2004; COSTA & CASTRO, 2007). Para a solução numérica da Equação de Richards, diferentes técnicas matemáticas têm sido aplicadas: Método das Diferenças Finitas (MILLER et al., 2006), Método de Elementos Finitos (HE & REN, 2009; JAVADI et al., 2008), Método dos Volumes Finitos (CAPUTO & STEPANYANTS, 2008), Transformada de Laplace (WENDLAND & VILHENA, 2001).

Neste trabalho, teve-se como objetivo a resolução numérica da Equação de Richards, cuja solução aproximada foi obtida pelo Método de Elementos Finitos (ZIENKIEWICZ, 1977), e a derivada temporal foi aproximada pelo Método de Diferenças Finitas (RÜBENKÖNIG, 2006; ROCHA et al., 2008), aplicando-se Euler Explícito (SPERANDIO et al., 2006). A referida equação foi resolvida utilizando-se de malhas uniformes inicialmente, e com a finalidade de obter simulações mais eficientes, a um custo computacional reduzido, empregou-se a adaptatividade com refinamento "h" na malha de elementos finitos (PIZARRO, 2009).

MATERIAL E MÉTODOS

Considerando o escoamento unidimensional não saturado e a coordenada vertical orientada positivamente para cima (BUNSRI et al., 2008 a), a Equação de Richards, para determinar os valores do potencial matricial de água, é dada por:

em que,

- capacidade hídrica específica do solo, L^-1;

θ - umidade volumétrica, L³L^-3;

ψ - potencial matricial, L;

K(ψ) - condutividade hidráulica do solo não saturado, LT^-1;

t - tempo, T, e

z - coordenada vertical, L.

As condições inicial e de contorno (condições de Dirichelet) utilizadas para a resolução da eq.(1) são: ψ(z,0) = ψ^inicial(z); ψ(0,t) = ψ₀ e ψ(-L,t) = ψ_L,-L< z <0.

Para resolver a Equação de Richards, relações constitutivas precisam ser especificadas. A relação entre θ e ψ é dada pela equação de GENUCHTEN (1980).

Com m = 1 - e a saturação efetiva

em que,

θ_s - umidade volumétrica do solo saturado, L³ L^-3;

θ_r - umidade volumétrica residual do solo, L³ L^-3;

α - parâmetro que depende do solo, L^-1, e

m e q - parâmetros que dependem do solo.

A relação entre K e θ é dada por:

em que,

K_s - condutividade hidráulica saturada, L T^-1.

A capacidade hídrica específica do solo é dada por:

A velocidade da água nos poros em solo não saturado ν_z (ROSSI et al., 2007), na direção z, é dada por:

Método dos resíduos ponderados

Para a solução numérica da Equação de Richards [eq.(1)], aplicou-se o Método de Elementos Finitos, cuja fundamentação matemática é descrita na sequência abaixo:

Sejam o domínio Ω = (-L, 0), o espaço funcional:

E o subespaço de funções-teste dado por:

em que,

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 a obtenção da solução, multiplica-se o resíduo da eq.(1) por uma função-teste ν ∈ V e integra-se sobre o domínio de Ω. Aplicando-se integração por partes na primeira parcela do membro direito da eq.(1) e como as funções-teste ν pertencem ao espaço V, ν (0) = ν (-L) = 0, tem-se:

A derivada temporal é aproximada por um quociente de diferença finita (SPERANDIO et al., 2006), e aplicando-se Euler Explícito (BEAR, 2007), obtém-se:

Assim, resolver a eq.(1) consiste em encontrar ψⁿ⁺¹(z,t) ∈ H¹ (Ω) que atenda às condições iniciais ψ (z,0) = ψ^inicial(z) e de contorno ψ(0,t) = ψ₀ e ψ(-L,t) = ψ_L, e que satisfaça:

para qualquer ν ∈ V.

Aproximação de elementos finitos

A aproximação de Galerkin (BUNSRI et al., 2008b) consiste em aproximar as funções variável de estado ψ ∈ H¹(Ω) e teste ν ∈ V por funções aproximadas e , com Π¹ um subespaço finito de H¹ e . Neste trabalho, o subespaço Π¹ foi 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: e . Substituindo-se em (10), considerando-se:

e agrupando-se os termos, obtém-se:

A eq.(12) deve ser satisfeita para qualquer . Logo, tomando-se uma função φ_i por vez, pode-se escrever a eq.(12) em forma matricial:

Assim, o Método de Elementos Finitos foi aplicado para encontrar os coeficientes multiplicadores {α}ⁿ⁺¹ que satisfazem o problema algébrico definido pela eq.(13). A solução desse problema algébrico é obtida através da decomposição de Cholesky (PRESS, 2007).

Parâmetros das aplicações

A primeira aplicação (CELIA et al., 1990) considera os seguintes parâmetros do solo: α = 3,35 m^-1; q = 2; K_s= 9,92 x 10^-5 m s^-1; θ_s = 0,368 m³ m^-3, e θ_r = 0,102 m³ m^-3. O comprimento da amostra de solo é de 1 m, e as condições inicial e de contorno são ψ(z,0) = -10 m, -1 < z < 0; ψ(0,t) = -0,75 m, t > 0 e ψ(-1,t) = -10 m, t > 0.

O domínio Ωσ = [-1, 0] foi subdividido uniformemente em 2⁶ = 64 elementos, e cada elemento Ω_e possui tamanho Δz = 1/64 m. 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 1). O tempo de execução para a obtenção da solução é de 1 min 23,169s.

Uma segunda aplicação da literatura (MIRANDA et al., 2005) tem como parâmetros de simulação θ₀ = 0,216 m³ m^-3; θ_i = 0,216 m³ m^-3; θ_r = 0; θ_s = 0,443 m³ m^-3; K_S = 1,515 x 10^-5 m s^-1; α = 4,49 m^-1; q = 3,6732, e m = 0,727758. As condições inicial e de contorno são dadas por ψ(z,0) = -0,686524 m, -0,70 < z < 0; ψ(0,t) = -0,062476 m, t > 0; ψ(0,70,t) = -0,686524 m, t > 0.

O domínio Ω = [-0,70, 0] foi particionado, na malha uniforme, em 128 elementos de igual tamanho. O tempo total de simulação é t = 1,75 h, em 450.000 passos de tempo com Δt = 1,75/450.000 horas = 0,014 segundos.

RESULTADOS E DISCUSSÃO

Com o objetivo de validar o código computacional desenvolvido, são apresentadas aplicações numéricas do modelo proposto. A primeira aplicação (CELIA et al., 1990) considera a infiltração de água em uma coluna de solo homogêneo, inicialmente seco. Os parâmetros do solo considerados estão citados no Material e Métodos.

Posteriormente, a malha uniforme fina foi adaptada com o objetivo de diminuir a quantidade de graus de liberdade do problema, reduzindo o tempo de execução do programa. Pode-se observar, pela Figura 2, que os resultados da simulação considerando-se funções polinomiais com p = 1 e com parâmetro de adaptação ε (|du_En| < ε max|du|, em que, |du_En| é gradiente em um elemento, max|du| é o máximo dos gradientes e 0 < ε < 1) no valor de ε = 0,01 não correspondem aos resultados obtidos com a malha uniforme fina.

Observou-se que a frente de molhamento 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.

A metodologia para evitar a perda de massa, observada no processo de refino ou desrefino da malha, foi a utilização de ordem de interpolação maior do que 1 (PIZARRO, 2009). Utilizando-se da ordem quadrática (p = 2), foram obtidas as soluções da Equação de Richards, em que o parâmetro de adaptação ε assume os valores ε = 0,01 e ε = 0,001. As soluções com ε = 0,01 e ε = 0,001 são bastante satisfatórias, em relação aos resultados da malha uniforme fina com 2⁶ elementos e p = 1, sem o problema de perda de massa da estratégia de adaptação, bastando observar a Figura 1 e os nós da malha adaptada com p = 2 e ε = 0,01, e a malha adaptada com p = 2 e ε = 0,001 (Figuras 3A e 3B, respectivamente).

Observou-se que a solução obtida, considerando a malha adaptada com p = 2 e ε = 0,001, está bem próxima da solução com a malha uniforme fina (2⁶ elementos e p = 1), embora a solução obtida, considerando-se a malha adaptada (p = 2 e ε = 0,01), seja uma boa aproximação em relação à malha fina, devido à quantidade bem menor de nós em relação à adaptada (p = 2 e ε = 0,001) (redução de 41% no número de nós) (Figura 3). 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 é, respectivamente, 23 e 39 nós, e o tempo de processamento é, respectivamente, 1 min 28,390 s e 1 min 42,474 s.

O tamanho do menor elemento da malha adaptada é igual ao da uniforme, por isso foi utilizado o mesmo espaço de tempo para essas duas malhas, mantendo, assim, a condição de CFL. Na malha adaptada, são removidos os nós desnecessários, diminuindo-se o custo computacional. A distância máxima entre os nós varia a cada espaço de tempo. Os resultados obtidos na simulação do potencial matricial È, deste trabalho, estão em conformidade com o modelo de CELIA et al. (1990) (Figura 4). Os parâmetros de simulação da segunda aplicação da literatura (MIRANDA et al., 2005) estão descritos no Material e Métodos.

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 foi o mesmo dos elementos da malha uniforme, ou seja, 0,70/128 m, e o do maior elemento varia a cada espaço de tempo. O parâmetro de adaptação foi de ε = 0,01. O processo de adaptação da malha é realizado a cada 1.000 espaços de tempo. Para as duas simulações, adotou-se ordem de aproximação p = 2, constante para todos os elementos da malha.

Nas Figuras 5 e 6, mostram-se, respectivamente, as soluções obtidas do potencial matricial ψ e da umidade volumétrica θ na malha uniforme com 2⁷ elementos (A), e na malha adaptada com ε = 0,01(B), ambas com p = 2, em vários instantes de simulação.

Os resultados indicam a equivalência de solução para as duas malhas utilizadas, 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 da malha uniforme é de 21 min 38s, e o da malha adaptada é de 3 min e 5 s, redução de 85%, aproximadamente, no tempo de processamento. O consumo de memória é de 6 Mb quando se utilizou a malha uniforme e de 5,800 Mb, quando foi utilizada a adaptada.

O gráfico da umidade versus profundidade, observando-se o último passo de tempo (Figura 6), está em conformidade com os resultados apresentados por MIRANDA et al. (2005) (Figura7), que resolveu a Equação de Richards formulada em termos de umidade θ .

CONCLUSÕES

A simulação realizada pelo modelo proposto neste trabalho foi capaz de prever o perfil do potencial matricial È e, consequentemente, da umidade volumétrica θ, em que se constata um bom desempenho do modelo na zona não saturada do solo.

A metodologia empregada, com a utilização do refinamento "h", fez com que a malha uniforme fina seja adaptada, diminuindo-se os graus de liberdade e, com isso, minimizou-se o tempo de execução computacional, mas com eficiente qualidade de solução.

A solução obtida com os dados do trabalho de CELIA et al. (1990), considerando-se a malha adaptada com p = 2 e ε = 0,001, está bem próxima da solução com a malha uniforme fina com 2⁶ elementos e p = 1, embora a solução obtida, considerando-se a malha adaptada com p = 2 e ε = 0,01, seja uma boa aproximação em relação à malha fina, devido à quantidade menor de nós em relação à adaptada com p = 2 e ε = 0,001. Portanto, a escolha do parâmetro de adaptação ε é fundamental na solução do problema.

A resolução da Equação de Richards, utilizando-se do Método de Elementos Finitos na aproximação espacial, com ordem de interpolação quadrática (p = 2), evitou a perda de massa do processo de adaptação de malha implementado, mostrando-se equivalente aos modelos de CELIA et al. (1990), que utilizaram malhas uniformes e funções-base lineares por partes, e de MIRANDA et al. (2005), que resolveram a Equação de Richards de maneira numérica por um sistema de volumes finitos.

Recebido pelo Conselho Editorial em: 22-7-2009

Aprovado pelo Conselho Editorial em: 12-12-2009

• ALFARO SOTO, M.A.; KUMAYAMAA, D.M.; CHANG, H.K. Calibração de um reflectômetro para estudos do fluxo de água em solo não saturado. Geociências, Rio Claro, v.26, n.4, p.357-368, 2007.
• BEAR, J. Dynamics of fluids in porous media New York: Elsevier, 1972. 764 p.
• ______. Hydraulics of groundwater. New York: Dover Publications, 2007. 569 p.
• BUNSRI, T.; SIVAKUMAR, M.; HAGARE, D. Influence of dispersion on transport of tracer through unsaturated porous media. Journal of Applied Fluid Mechanics, Iran, v.1, n.2, p.37-44, 2008.
• BUNSRI, T.; SIVAKUMAR, M.; HAGARE, D. Numerical modelling of tracer transport in unsaturated porous media. Journal of Applied Fluid Mechanics, Iran, v.1, n.1, p.62-70, 2008.
• CAPUTO, J.G.; STEPANYANTS, Y.A. Front Solutions of Richards' Equation. Transport in Porous Media, Netherlands, v.74, n.1, p.1-20, 2008.
• CELIA, M.A.; BOULOUTAS, E.T.; ZARBA, R.L. A general mass conservative numerical solution for the unsaturated flow equation. Water Resources Research, Hawaii, v.26, n.7, p.1.483-1.496, 1990.
• 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.
• DARCY, H. Les Fountaines Publiques de la Ville de Dijon: appendice D. Paris: Dalmont, 1856.
• GENUCHTEN, M. T. Van. A closed-form equation for predicting the hydraulic conductivity of unsatured soils. Soil Science Society of America Journal, Pittsburgh, v.44, n.3, p.892-898, 1980.
• HE, X.; REN, L. An adaptive multiscale finite element method for unsaturated flow problems in heterogeneous porous media. Journal of Hydrology, Reston, v.374, p.56-70, 2009.
• 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.
• MILLER, C.T.; ABHISHEK, C.; FARTHING, M.W. A spatially and temporally adaptive solution of Richards'equation. Advances in Water Resources, Southampton, v.29, n.4, p.525-545, 2006.
• 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. Engenharia Agrícola, Jaboticabal, v.25, n.3, p.677-685, 2005.
• PIZARRO, M.L.P. Simulação de fluxo e transporte de solutos na zona não saturada do solo 2009.185 f. Tese (Doutorado em Ciências da Engenharia Ambiental) - Escola de Engenharia de São Carlos, Universidade de São Paulo, São Carlos, 2009.
• PRESS, W.H. Numerical recipes: the art of scientific computing. Cambridge: Cambridge University Press, 2007. 1.235 p.
• REICHARDT, K.; TIMM, L.C. Solo, planta e atmosfera - conceitos, processos e aplicações Barueri: Editora Manole, 2004. 478 p.
• RICHARDS, L.A. Capillary conduction of liquids through porous medium. Physics, New York, v.1, p.318-333, 1931.
• ROCHA, F.A.; MARTINEZ, M.A.; 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.
• ROSSI, P.; MIRANDA, J.H.; DUARTE, S.N. Curvas de distribuição de efluentes do íon nitrato em amostras de solo deformadas e indeformadas. Engenharia Agrícola, Jaboticabal, v.27, n.3, p.675-682, 2007.
• RÜBENKÖNIG, O. The finite difference method: an introduction. Albert-Ludwigs University of Freiburg, p.1-5, 2006.
• 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. 354 p.
• WENDLAND, E. Contribuição à simulação de processos em meios porosos 2004. 270 f. Tese (Livre-Docência) - Escola de Engenharia de São Carlos, Universidade de São Paulo, São Carlos, 2004.
• WENDLAND, E.; VILHENA, M.T. Simulation of infiltration in porous medium by Laplace Transform Technique and finite difference method. Hybrid Methods in Engineering, New York, v.3, p.1-9, 2001.
• ZIENKIEWICZ, O.C. The Finite element method New York: McGraw-Hill, 1977. 787 p.

Modelagem computacional do fluxo unidimensional de água em meio não saturado do solo

Computational model of one-dimensional flow of water in an unsaturated soil

Datas de Publicação

• Publicação nesta coleção
  20 Jul 2010
• Data do Fascículo
  Jun 2010

Histórico

• Aceito
  12 Dez 2009
• Recebido
  22 Jul 2009