On a Multigrid Method for the Coupled Stokes and Porous Media Flow Problem

. The multigrid solution of coupled porous media and Stokes ﬂow problems is considered. The Darcy equation as the saturated porous medium model is coupled to the Stokes equations by means of appropriate interface conditions. We focus on an e ﬃ cient multigrid solution technique for the coupled problem, which is discretized by ﬁnite volumes on staggered grids, giving rise to a saddle point linear system. Special treatment is required regarding the discretization at the interface. An Uzawa smoother is employed in multigrid, which is a decoupled procedure based on symmetric Gauss-Seidel smoothing for velocity components and a simple Richardson iteration for the pressure ﬁeld. Since a relaxation parameter is part of a Richardson iteration, Local Fourier Analysis (LFA) is applied to determine the optimal parameters. Highly satisfactory multigrid convergence is reported, and, moreover, the algorithm performs well for small values of the hydraulic conductivity and ﬂuid viscosity, that are relevant for applications.


INTRODUCTION
Coupling of free flow and a saturated porous medium models has received considerable attention due to its application in environmental and industrial context, such as in flood simulation, filtration, contamination, and so on.It is challenging to deal with a coupled system, since each part is based on a different model and an appropriate coupling at the interface is required.The fluid flow through a rigid and saturated porous medium Ω d is described by Darcy's law, which is an expression of conservation of momentum.The mixed formulation of the Darcy problem is natural for computations in the porous medium region since it allows for directly approximating the velocity.Therefore, in this work, we will consider this formulation which reads where u d = (u d , v d ) describes the velocity and p d the fluid pressure inside the porous medium.K is the hydraulic conductivity tensor, representing the properties of the porous medium and the fluid.Here, only the case K = KI, K > 0 is considered.Sinks and sources are described by the force term f d .
The free flow subproblem is modeled by using the Stokes equations for a viscous, incompressible, Newtonian fluid.It is a linearized form of the Navier-Stokes equations in the limit case when the nonlinear term becomes negligible.The motion of the Stokes flow in the region Ω f is described as represents a prescribed force, and the fluid stress tensor is given by σ f = −p f I + 2νD(u f ), with p f denoting the fluid pressure, ν representing the fluid viscosity and where D(u f ) = (∇u f + (∇u f ) T )/2 is the strain tensor.
The Darcy and Stokes systems must be coupled across the internal interface Γ by adequate conditions.To describe them, we fix the normal vector to the interface to be n = n f = −n d and we denote τ as the tangential unit vector at the interface Γ. Across Γ the continuity of fluxes and normal stresses must be imposed.This gives rise to the following two standard coupling conditions on Γ: (1) Mass conservation: As a third coupling condition, the so-called Beavers-Joseph-Saffman interface condition is widely used, which is supported by experimental findings and rigorous mathematical theory of homogenization.This condition relates the tangential velocity along the interface with the fluid stresses, that is, αu f •τ+τ•σ f •n = 0 on Γ , where α is a dimensionless parameter which needs to be experimentally determined and depends on the properties of the porous medium.An alternative to this third interface condition neglects the second term in the Beavers-Joseph-Saffman condition, giving rise to a no-slip interface condition, u f • τ = 0 on Γ .Note that there exists a rich study on numerical algorithms for the coupled free flow and unsaturated porous media flow model which results to be a nonstationary and strongly nonlinear system.To apply the monolithic multigrid approach proposed here for such a coupled system can be a possible extension to the present work.

NUMERICAL METHODS
The finite volume method [1] on a staggered grid is considered as the discretization scheme for the coupled Darcy/Stokes problem.By using this discretization we ensure that spurious oscillations do not appear in the numerical solution, and we obtain a mass conservative algorithm for the whole system.The discretization at the interface is of great importance and can be viewed as a relevant ingredient towards the construction of a highly efficient multigrid method.Since the coupled system is treated as a single problem, the equations of free flow, porous media flow and their complex interaction are all included in one discrete formulation.By such a discretization, the fully coupled system possesses a saddle-point structure which is suitable for monolithic multigrid.
As it is well-known, the efficiency and robustness of a multigrid method depend on the choice of components of the algorithm.The choice of smoother requires special attention due to the saddle point structure of the considered system.An Uzawa smoother, which was proposed for the Stokes problem in [3], is considered for the coupled system.It shall be interpreted as a combination of the symmetric Gauss-Seidel smoothing for velocities, together with a Richardson iteration for the Schur complement in the pressure field.The Richardson iteration involves a relaxation parameter which affects the convergence speed, and has to be carefully determined.The analysis of the smoother is based on the framework of local Fourier analysis (LFA) and it allows us to provide an analytic bound of the smoothing factor of the smoother as well as an optimal value of the relaxation parameter.A detailed study of the Uzawa smoother in the framework of LFA was already done in [3] for a family of Stokes problems.We work out the analysis for Darcy's equation.The relaxation parameter for the Richardson iteration is given by the expression ω = h 2 /5K.Parameter ω depends on the grid size h, and therefore it will be different on each grid of the hierarchy used in the multigrid method.
Due to the saddle point structure of the coupled problem, a geometric multigrid method together with an Uzawa smoother, can be applied for the whole system.For this purpose, it is important to note that to keep the structure of the matrix of the saddle point system on the whole grid hierarchy, interface Γ has to be present on each grid level.Regarding the smoothing process, all velocity unknowns are relaxed before the pressure unknowns will be updated.The relaxation parameter ω for the Richardson iteration for the Schur complement has to be chosen differently if we are updating pressure unknowns from the Darcy or the Stokes problems.For the rest of the components, the same operators can be used at every grid point since the discretization for both problems is performed with the same staggered arrangement of unknowns.The inter-grid transfer operators that act on the different unknowns are defined as follows: At velocity grid points six-point restrictions are considered, and at pressure grid points a four-point restriction is applied.For the prolongation operators, we choose the adjoints of the restrictions.In the monolithic multigrid method we do not distinguish the subproblems and the internal interface.All the unknowns play essentially the same role.Only the relaxation parameter of the smoother is different for each subproblem.For the discretization at the interface, the unknowns for both subproblems are included in one equation.We keep the same philosophy for the other components in the monolithic multigrid.For example (Step 6 in the algorithm), to restrict the unknowns at the interface, six points from both subgrids around it are employed.A suitable discretization for the unknowns at the interface of the coupled system is a key step in achieving robustness and efficiency of our approach.Note that the interface Γ needs to be present at all grids of the hierarchy, so that the correct coupling done by interface discretization is guaranteed on all grid levels.The proposed multigrid method for the coupled Darcy/Stokes problem can also be implemented as a multiblock version in which the Darcy and Stokes domains are assumed to be two different blocks.This is appealing from a practical point of view, for example when one has to solve the coupled problem by using 560023-2 two different codes.Moreover, this multiblock approach is easily parallelizable.Next, we describe in detail how this implementation can be done.

Multiblock multigrid algorithm
We divide our domain into two different blocks corresponding to the Darcy and Stokes domains.In this way, the original staggered grid is split into two different sub-grids.Since in this version of the algorithm it is necessary to transfer information between both blocks, the mesh corresponding to the Stokes domain is extended by adding an overlap region of one cell length, as can be seen in Figure ?? (Le f t).Next, we explain in detail the two-grid version of the multiblock algorithm.For simplicity in the presentation of the algorithm, we use pre-smoothing but no post-smoothing.By recursion, the multigrid version follows straightforwardly.This multiblock algorithm requires only little data communication.In particular, each communication step involves transfer of information in only one way.Moreover, each stage in the algorithm can be performed in parallel since the data required for each operation is available in the same process.Finally, although this multiblock approach can be cast into the class of domain decomposition (DD) methods, we wish to emphasize that in our case the communication between both Darcy and Stokes problems is performed on each level in the hierarchy instead of only on the finest grid as usual in the DD methods.This is crucial to achieve a highly efficient solver for this coupled problem, as we will see in the numerical experiments section.

NUMERICAL RESULTS
Here we present a more complicated and realistic numerical test.This test addresses the coupling of Darcy and Stokes problem which is in a cross-flow filtration setting.The data in this test are taken from the experiment presented in [2] which is a micro-membrane filtration model.The domain of the coupled problem is shown in Figure ?? (Right).Ω f represents a channel on the top where the flow can go through, while Ω d represents a filter.Since the lengths of the free flow domain and the porous medium are not the same, the coupled domain is divided into four different blocks corresponding to the Darcy (Block1) and Stokes (Block2, Block3 and Block4) domains.The two-block multigrid algorithm described before can be straightforwardly adapted for fours blocks.The information transfer between Block1 and Block3 is the same as before.For the Stokes domain, two artificial boundaries are generated by the partitioning.As the communication between the subgrids in Ω f is necessary, an overlap region of one cell length is created for Block2 and Block4 along the artificial boundaries.The data located in the overlap region is computed and transferred from the neighboring subgrid in Block3.The unknowns at the artificial boundaries, i.e. u f , are updated in Block2 and Block4, and then sent to Block3.The communication is implemented on each level in the multigrid algorithm.The inflow entering into the domain Ω f is specified.At the interface, the Beavers-Joseph condition is imposed.At the bottom of the porous medium, the pore pressure is set as zero.There is an exit (see Figure ?? (Right) in red) at the right-vertical boundary of the free flow domain.The height of the exit is 0.00125 which is quite small compared to the inlet.The stress-free boundary condition is employed at the exit where the flow may leave the the domain freely.All the other imposed conditions are shown in Figure ?? (Right).Two values of hydraulic conductivity K = 0.1 and K = 10 −6 are considered in the numerical experiment.The fluid viscosity is chosen as 10 −6 .For a multigrid W(2, 2)−cycle, the multigrid convergence factor is around 0.2 for all cases, and the multigrid method exhibits a highly satisfactory behavior.LFA is also used to confirm the convergence obtained from the monolithic multigrid method.LFA is applied to both Darcy and Stokes subproblems separately, and it is shown that the worst of these factors results to be the global convergence of the multigrid for the coupled problem.Here we obtain 0.2 as the two-grid convergence factor for W(2, 2)−cycle which matches perfectly with the asyptotic convergence factor.In Figure 2 (a), we show the velocity vector corresponding to K = 0.1.Since the hydraulic conductivity of the porous medium is quite high, when the fluid travels tangentially cross the interface, the majority of the flow seeps into the filter.While only a small amount of fluid goes through the exit of the channel.In Figure 2 (b), the velocity vector corresponding to K = 10 −6 is represented.With such a low hydraulic conductivity of the porous medium, the minority of the flow penetrates the interface.Whereas most of the fluid flows towards the small exit of the channel.

FIGURE 2 .
FIGURE 2. Velocity vectors over the cross-flow filtration domain with different values of hydraulic conductivity.