Numerical coupling of Navier-Stokes and Darcy flow for soil-water evaporation

Abstract

The objective of this work is to develop algorithms and provide a framework for an efficient coupling of free flow and porous-medium flow to simulate porous-medium-soil-water evaporation. The implementation must particularly be capable of simulating laminar free flows, be fast enough for applied research, and cover simulations in two and three dimensions with complex geometries.

We introduce a model for a compositional non-isothermal free flow coupled with a two-fluid-phase compositional non-isothermal porous-medium flow. The free flow is modeled with the Navier-Stokes, component and energy transport equations. The porous-medium flow is modeled with compositional two-fluid-phase Darcy and energy transport equations. As the pressure has different orders in the free-flow and porous-medium-flow subdomains, the coupling is not straightforward. Although the simulation of the coupled flows is motivated by a laboratory experiment to measure soil-water evaporation caused by wind blowing over a water-filled porous bed, we intend to also explore its use in other applications

The free flow is considered to be incompressible and laminar. We also assume that air and water follow nonlinear laws that describe their physical properties, and binary diffusion. Within the porous medium only creeping flows occur. Many quantities are averaged and used in a macroscopic sense. We use a formulation of two-phase Darcy law using the liquid saturation and the gas pressure as primary variables. The component mass fractions are calculated by Henry's law and the vapor pressure. The liquid phase may locally vanish leading to a variable switch, where the vapor mass fraction is tracked instead of the liquid saturation. We assume that a local thermodynamic equilibrium is valid everywhere within the domain, even across the interface. We follow the coupling concept proposed by Mosthaf et al. (2011), including the Beavers-Joseph-Saffman approach which has a sharp interface between the two subdomains.

We use a cell-centered finite volume method (FVM) on an axially parallel grid to discretize the partial differential equations of the compositional two-phase Darcy's law, the heat equation in both subdomains, and the component transport in the free-flow domain. For the Navier-Stokes equation, we use the marker and cell (MAC) scheme which moves the degrees of freedom for the velocities towards the edges of the grid elements, forming one secondary, staggered grid per dimension. The MAC scheme is stable and can be interpreted as a FVM. The coupling conditions are applied without additional variables along the coupling interface. They are incorporated as Dirichlet, Neumann or Robin boundary conditions resulting in interface fluxes.

For the porous-medium flow, we use the finite-volume implementation provided by DuMuX. The marker and cell scheme is implemented on top of Dune-PDELab utilizing the material laws from DuMuX. The grid is split into two subdomains and the grid elements can be graded. This is especially useful for developing smaller elements closer to the interface. We can use complex geometries in two or three dimensions. The coupling is provided by a Dune-Multidomain local coupling operator. The time integration is approximated with an implicit Euler scheme and an adaptive time stepping. The system of nonlinear equations is linearized by a Newton method. All contributions to the Jacobian are compiled in one system of linear equations.

The resulting matrices are difficult to solve. Although they are sparse, with a blocked structure of bands of nonzero entries, the matrices contain a saddle point problem and are nonsymmetric. We solve the matrices with direct methods. We also investigate iterative methods to get around the computational complexity and memory consumption of the direct methods: An Algebraic Multigrid (AMG) method, a Schur complement method, and a Generalized Minimal Residual method (GMRES) preconditioned with the reordering algorithm MC64 and an incomplete LU factorization with threshold and pivoting (ILUTP). We experience problems with AMG's error criteria leading to convergence problems. The Schur complement method is slow, as the Schur complement, which is not explicitly calculated, lacks preconditioners. GMRES with ILUTP shows similar results to a direct method, but reveals a restriction on the time step size for larger problem sizes, flawing a possible speedup compared to the direct methods.

We validate our implementation for proper operation with the simulation of a laboratory experiment for soil-water evaporation. The laboratory experiment consists of a water-filled sand box with a horizontal pipe installed on top of the box and a propeller creating a constant air flow. We use the implementation to investigate the influence of the Reynolds number on the evaporation rate. Further, we compare the two-dimensional simplification to different three-dimensional geometries with regard to the effects on the evaporation. For low Reynolds numbers, the geometry of the free-flow subdomain has a significant influence on the evaporation rate.

Another application involves a geological repository for nuclear waste. We investigate the water saturation in the concrete ceiling and the rock above a ventilation gallery. Our results conclude that within the first 200~years, only part of the concrete will dry, and the rock will remain unaffected. This confirms the same result by another group, though they observe evaporation rates up to the factor of ten higher.

Our third application is the water management within a polymer electrolyte membrane (PEM) fuel cell. Neglecting electrochemistry, we simulate the flow through the gas channels and the porous layer covering the membrane, including the transport of vapor and liquid water, the evaporation of water within the porous layer, and how energy and vapor are conveyed away. In comparison to the above applications, the gas phase flow is not horizontally parallel to the porous bed, but is forced to completely enter the porous medium and leave it through a second gas channel. We also briefly compare two different gas channel layouts.

We introduce the discretization of the coupling concept and its implementation. We conduct simulations of applications from different areas. We show the versatility of our approach and that it can be used as the basis for further research.