Commit 3bb596eb by Kilian Weishaupt

 # Exercise #5 (DuMuX course) The aim of this exercise is to get familiar with the set up of coupled free flow/porous medium flow problems. ## Problem set-up The model domain consists of two non-overlapping subdomains. Free flow is modeled in the upper subdomain, while the lower subdomain models a flow in a porous medium. Both single-phase flow and two-phase flow will be considered in the porous domain. ![](../extradoc/exercise5_setup.png) ### 0. Getting familiar with the code * Navigate to the directory exercises/exercise-coupling-ff-pm There are three sub folders: interface, models and turbulence. The problem-related files for this exercise are: * Three __main files__ for the three sub-tasks :ex_interface_coupling_ff-pm.cc, ex2_coupling_ff-pm.cc, ex2_coupling_ff-pm.cc, * Three __free flow problem files__: ex_interface_ffproblem.hh, ex2_ffproblem.hh, ex3_ffproblem.hh * Three __porous medium flow problem files__: ex_interface_pmproblem.hh, ex2_pmproblem.hh, ex3_pmproblem.hh * The __input files__: ... * The __spatial parameters files__: 1pspatialparams.hh, 2pspatialparams.hh In the main file, TypeTags for both submodels are defined. The same applies for types such as GridManager, FVGridGeometry, Problem, etc.. Since we use a monolithic coupling scheme, there is only one Assembler and one NewtonSolver. The problem files very much look like "regular", uncoupled ones with the exception that they hold a pointer to the CouplingManager which allows to evaluate the coupling conditions and to exchange information between the coupled models. The coupling conditions are realized technically in terms of boundary condition. For instance, in lines 178 and 179 in ex_interface_ffproblem.hh, couplingNeumann boundary conditions are set, which means that the free flow models evaluates the mass and momentum fluxes coming from the porous domain and uses these values as boundary conditions at the interface. Note the certain checks are performed when combining different models, e.g., the fluid system has to be the same for both domains. We will use a staggered grid for the free flow and a cell-centered finite volume method for the porous medium. Keep in mind that the staggered grid implementation distinguishes between face variables (velocity components) and cell center variables (all other variables), therefore in some cases either the stokesCellCenterIdx or the stokesFaceIdx is used respectively, while for the porous medium all variables can be accessed with darcyIdx. The main part of the Stokes/Darcy coupling is implemented in the folder dumux/multidomain/boundary/stokesdarcy. __Task__: Take a closer look at the Stokes/Darcy coupling files and try to answer the following questions before moving to the next part of the exercise: * Where are the elements at the coupling interface paired with their respective coupling partner elements of the other domain? * Where are the coupling conditions implemented and evaluated? * You will find four implementations of the massCouplingCondition. What are the differences between them? * Which data is stored and updated regularly such that it is accessible from the respective other domain? * How can we access volume variables like density and pressure from the respective other domain? * ... ### 1. Changing the interface In this part of the exercise, a simple coupled system consisting of a one-phase (1p) free flow and a one-phase flow in a porous medium is set up. Both subproblems have no-flow boundaries at the sides. A velocity profile is set on the upper free flow boundary, which leads to a vertical flow into the porous medium. * We will first change the flow direction such that the free flow is parallel to the porous medium. * Afterwards, the Beavers-Joseph-Saffman condition will be used as an interface condition for the tangential momentum transfer. * Last, we change the flat interface between the two domains to a wave-shaped one. __Tasks__: Open the file ex_interface_ffproblem.hh and navigate to line 155, where the types of boundary condition are set. Instead of applying a fixed velocity profile at the top of the domain, we want to use fixed pressure boundary conditions at the left and right side of the free flow domain, while the top represents an impermeable wall. __Change the flow direction__ Set a Dirichlet boundary condition for the pressure at the left and ride side of the domain:  cpp if(onLeftBoundary_(globalPos) || onRightBoundary_(globalPos)) values.setDirichlet(Indices::pressureIdx);  Set a Dirichlet boundary condition for the velocities at the top:  cpp if(onUpperBoundary_(globalPos)) { values.setDirichlet(Indices::velocityXIdx); values.setDirichlet(Indices::velocityYIdx); }  Keep the coupling boundary condition:  cpp if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf)) { values.setCouplingNeumann(Indices::conti0EqIdx); values.setCouplingNeumann(Indices::momentumYBalanceIdx); }  Having changed the types of boundary conditions, we must now assign the correct values for them. Set a no-slip, no-flow condition for the velocity at the top:  cpp values[Indices::velocityXIdx] = 0.0; values[Indices::velocityYIdx] = 0.0;  Apply a fixed pressure difference between the inlet and outlet, e.g.:  cpp if(onLeftBoundary_(globalPos)) values[Indices::pressureIdx] = deltaP_; if(onRightBoundary_(globalPos)) values[Indices::pressureIdx] = 0.0;  For changing the flow direction, the boundary conditions for the porous medium have to be changed as well. Use Neumann no-flow boundaries everywhere, keep the coupling conditions.  cpp values.setAllNeumann(); if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf)) values.setAllCouplingNeumann();  __Include slip-condition__ This should make the flow go from left to right. However, we are still missing one important feature: the Beavers-Joseph-Saffman slip condition at the interface: $\frac{\partial v_x}{\partial y} = \frac{\alpha}{\sqrt K} (v_x - q_{pm})\quad$ at $\quad y=0$ with $\quad q_{pm}=0$. The include this, we just set:  cpp values.setBJS(Indices::momentumXBalanceIdx);  at the position where the coupling boundary conditions are set in ex_interface_ffproblem.hh. To check if the simulation behaves as expected, we can compare the velocity profile $v_x(y)$ with the analytical solution provided by [Beavers and Joseph (1967)](https://doi.org/10.1017/S0022112067001375). For doing so, we uncomment line 212 cpp stokesVtkWriter.addField(stokesProblem->getAnalyticalVelocityX(), "analyticalV_x");  in ex_interface_coupling_ff-pm.cc. After re-compiling and re-running the executable, we should be able to see also the analytical solution of $v_x$ on the free flow domain. Play around with the grid resolution to see how that affects the velocity profile. __Cange shape of interface__ Now we want to include a non-flat interface between the two domains. We use dune-subgrid to construct two grids for the two domains from one common host grid. Comment out lines 93-103 in ex_interface_coupling_ff-pm.cc and comment lines 112-147 in the same file. This will instantiate a host grid and define two helper lambda functions that are used to choose elements from to host grid for the respective sub grid. In the given case, the domain is split in two haves, separated by a sinusoidal interface. cpp auto elementSelectorStokes = [&](const auto& element) { double interface = params.amplitude * std::sin(( element.geometry().center()[0] -params.offset) / params.scaling * 2.0 * M_PI) + params.baseline; return element.geometry().center()[1] > interface; }; auto elementSelectorDarcy = [&](const auto& element) { double interface = params.amplitude * std::sin(( element.geometry().center()[0] - params.offset) / params.scaling * 2.0 * M_PI) + params.baseline; return element.geometry().center()[1] < interface; };  Make sure to comment in line 30 in both problem files cpp #include  and do the changes in the respective lines for the Grid property. The problem should no compile and run. However, an error occurs due to the coupling conditions. So far, we assumed a flat interface, therefore the normal momentum coupling condition $[\sigma \cdot \mathbf{n}]^{FF} = p^{PM}$ was always set for a fixed $\mathbf{n} = (0,1)^T$. We need to account for the curvature of the interface and thus replace cpp values.setCouplingNeumann(Indices::momentumYBalanceIdx);  with cpp values.setCouplingNeumann(scvf.directionIndex());  The same if true for the BJS condition, however, here we need to consider the tangential direction: cpp values.setCouplingNeumann(1 - scvf.directionIndex());  The final result should look something like this: ![](../extradoc/ex_ff-pm-wave-interface.png) ### 3. Change the models 1p2c/1p2c -->1p2c/2p2c __Task__: ... ### 4. ...
 // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- // vi: set et ts=4 sw=4 sts=4: /***************************************************************************** * See the file COPYING for full copying permissions. * * * * This program is free software: you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation, either version 2 of the License, or * * (at your option) any later version. * * * * This program is distributed in the hope that it will be useful, * * but WITHOUT ANY WARRANTY; without even the implied warranty of * * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * * GNU General Public License for more details. * * * * You should have received a copy of the GNU General Public License * * along with this program. If not, see . * *****************************************************************************/ /*! * \file * \ingroup OnePTests * \brief The spatial parameters class for the test problem using the 1p cc model */ #ifndef DUMUX_1P_TEST_SPATIALPARAMS_HH #define DUMUX_1P_TEST_SPATIALPARAMS_HH #include namespace Dumux { /*! * \ingroup OnePModel * \ingroup ImplicitTestProblems * * \brief The spatial parameters class for the test problem using the * 1p cc model */ template class OnePSpatialParams : public FVSpatialParamsOneP> { using FVGridGeometry = typename GET_PROP_TYPE(TypeTag, FVGridGeometry); using Scalar = typename GET_PROP_TYPE(TypeTag, Scalar); using Problem = typename GET_PROP_TYPE(TypeTag, Problem); using GridView = typename GET_PROP_TYPE(TypeTag, GridView); using ParentType = FVSpatialParamsOneP>; using Element = typename GridView::template Codim<0>::Entity; using GlobalPosition = typename Element::Geometry::GlobalCoordinate; public: // export permeability type using PermeabilityType = Scalar; OnePSpatialParams(std::shared_ptr fvGridGeometry) : ParentType(fvGridGeometry) { permeability_ = getParam("Darcy.SpatialParams.Permeability"); alphaBJ_ = getParam("Darcy.SpatialParams.AlphaBeaversJoseph"); } /*! * \brief Function for defining the (intrinsic) permeability \f$[m^2]\f$. * * \param globalPos The global position * \return the intrinsic permeability */ PermeabilityType permeabilityAtPos(const GlobalPosition& globalPos) const { return permeability_; } /*! \brief Define the porosity in [-]. * * \param globalPos The global position */ Scalar porosityAtPos(const GlobalPosition& globalPos) const { return 0.4; } /*! \brief Define the Beavers-Joseph coefficient in [-]. * * \param globalPos The global position */ Scalar beaversJosephCoeffAtPos(const GlobalPosition& globalPos) const { return alphaBJ_; } private: Scalar permeability_; Scalar alphaBJ_; }; } // end namespace #endif