README.md

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. Currently, a velocity profile is set on the upper free flow boundary, which leads to a vertical flow into the porous medium:

Note that we neglect the influence of gravity and only solve for the stationary problem in this exercise.

You can also compile and run the test without changing anything yet and investigate the above mentioned velocity distribution with paraview.

• 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.

Task A: Change the flow direction

Open the file interface/freeflowsubproblem.hh and navigate to the part, 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.

Set a Dirichlet boundary condition for the pressure at the left and right side of the domain:

if(onLeftBoundary_(globalPos) || onRightBoundary_(globalPos))
    values.setDirichlet(Indices::pressureIdx);

Set a Dirichlet boundary condition for the velocities at the top:

if(onUpperBoundary_(globalPos))
{
    values.setDirichlet(Indices::velocityXIdx);
    values.setDirichlet(Indices::velocityYIdx);
}

Keep the coupling boundary condition:

if(couplingManager().isCoupledEntity(CouplingManager::stokesIdx, scvf))
{
    values.setCouplingNeumann(Indices::conti0EqIdx);
    values.setCouplingNeumann(Indices::momentumYBalanceIdx);
    values.setDirichlet(Indices::velocityXIdx); // assume no slip on interface
}

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:

values[Indices::velocityXIdx] = 0.0;
values[Indices::velocityYIdx] = 0.0;

Apply a fixed pressure difference between the inlet and outlet, e.g.:

if(onLeftBoundary_(globalPos))
    values[Indices::pressureIdx] = deltaP_;
if(onRightBoundary_(globalPos))
    values[Indices::pressureIdx] = 0.0;

This should make the flow go from left to right. If not already done, you can also delete the vertical velocity from dirichletAtPos, as Dirichlet-type boundary conditions should be set for the pressure.

For changing the flow direction, the boundary conditions for the porous medium (in interface/porousmediumsubproblem.hh) have to be changed as well.

Use Neumann no-flow boundaries everywhere and keep the coupling conditions.

values.setAllNeumann();

if (couplingManager().isCoupledEntity(CouplingManager::darcyIdx, scvf))
    values.setAllCouplingNeumann();

Recompile and rerun your code. Then check your results with the applied changes by visualizing them e.g. with paraview.

Task B: Include slip-condition

However, we are still missing one important feature: at the moment, the velocity component tangential to the interface gets a no-slip condition. In the next step we want to implement 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

.

To include this, just replace the no-slip condition at the interface

values.setDirichlet(Indices::velocityXIdx); // assume no slip on interface

with a Beavers-Joseph-Saffman (BJS) boundary condition for the respective momentum balance:

values.setBeaversJoseph(Indices::momentumXBalanceIdx);

at the position where the coupling boundary conditions are set in interface/freeflowsubproblem.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). For doing so, we uncomment the following line in main.cc in the subfolder interface.

freeflowVtkWriter.addField(freeflowProblem->getAnalyticalVelocityX(), "analyticalV_x");

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 in the input file to see how that affects the resulting velocity profile.

Note: In Paraview you can e.g. use the tool PlotOverLine to plot the velocity in

x

-direction (

v_x

) over the y-Axis (e.g. at the outlet or the middle of the domain). After choosing PlotOverLine make sure to give the Point 1 and the Point 2 to define the line or click on the button Y Axis and check the coordinates of the points. Click on Apply and in section Series Parameters choose to only visualize velocity_gas(m/s)_X and analyticalV_x.

Task C: Change the 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. Thus, for the following tasks subgrid needs to be correctly installed. Our hostgrid will be a Dune-Yasp grid (Dune::YaspGrid<2, Dune::TensorProductCoordinates<double, dim> >) and the bounds and resolution of the grid will be set in the params.inputfile under the group [Grid]. This hostgrid, along with elementSelector functions defining some spatial cut of this domain, are passed to the grid manager to create each subdomain grid.

To introduce this, open main.cc in the subfolder interface again and search for TODO: dumux-course-task 1.C. Comment out the first code block and uncomment the second. This will instantiate a host grid and define two helper lambda functions that are used to choose elements from the host grid for the respective sub grid. In the given case, the domain is split in two halves, separated by a sinusoidal interface.

auto elementSelectorFreeflow = [&](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 elementSelectorPorousMedium = [&](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 that you have uncommented the line for including the grid manager in the properties file, i.e.

#include <dumux/io/grid/gridmanager_sub.hh>

and do the changes in the respective lines for the Grid property.

As well we have to define the parameters only for the base grid if using dune-subgrid. For this uncomment the block [Grid] in interface/params.input and comment the groups [Freeflow.Grid] and [PorousMedium.Grid].

Note: The four additional parameters Baseline, Amplitude, Offset and Scaling define the sinus-shaped interface in the main-file (the block you previously uncommented).

The problem should now compile. However, a runtime 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

values.setCouplingNeumann(Indices::momentumYBalanceIdx);

with

values.setCouplingNeumann(scvf.directionIndex());

in freeflowsubproblem.hh in the subfolder interface.

The same is true for the BJS condition, however, here we need to consider the tangential direction:

values.setBeaversJoseph(1 - scvf.directionIndex());

Recompile and rerun to check if the final result looks something like this:

Extra Points: Rather than enforcing a pressure difference across the domain, an inflow velocity profile could be set. What changes to the left boundary conditions in the free-flow domain would you make to introduce this? What conditions can be enforced on the right boundary? Hint: A relation between velocity and position is used for the vertical velocity component in the original form of the dirichletAtPos method.