Simple channel flow using the Incompressible Solver

What's new

• running a simulation with the incompressible Navier-Stokes solver
• how to use the immersed boundary method
• post-processing, i.e. plotting and checking physical values

Prerequisites

• the BoSSS framework
• a visualization tool, e.g Paraview or VisIt
• the knowledge of how to setup a database for BoSSS

This tutorial will explain the basic features of the incompressible Navier-Stokes solver in the BoSSS framework.

First, the simple testcase of a 2D channel flow will be explained. After that, there will be a short part about the immersed boundary feature of our incompressible flow solver.

Therefore the flow around a cylinder will be investigated using the immersed boundary method.

Note that BoSSS, at the present time contains no stand-alone single-phase solver that is fully recomended - although there are some legacy solvers, e.g. SIMPLE. Instead, the two-phase-solver with immersed boundary is used, where the two-phase option ist deactivated.

Problem statement

The flow is described by the unsteady Navier-Stokes equations in the fluid region

$$ \rho_f\left(\frac{\partial \vec{u}}{\partial t}+ \vec{u} \cdot \nabla \vec{u}\right) +\nabla p - \mu_f \Delta \vec{u} = \vec{f} $$

and the continuity equation

$$ \nabla \cdot \vec{u} = 0 \quad \forall\ t \in (0,T)\quad \textrm{in}\ \Omega $$

In the equations above

• $\vec{u}$ is the velocity vector
• $p$ the pressure.
• The fluid density is denoted by $\rho_f$
• $\mu_f=\rho_f \cdot \nu_f$ is the dynamic viscosity of the fluid.

Channel

First, we initialize the new worksheet; Note:

1. This tutorial can be found in the source code repository as as channel.ipynb. One can directly load this into Jupyter to interactively work with the following code examples.
2. In the following line, the reference to BoSSSpad.dll is required. You must either set #r "BoSSSpad.dll" to something which is appropirate for your computer (e.g. C:\Program Files (x86)\FDY\BoSSS\bin\Release\net5.0\BoSSSpad.dll if you installed the binary distribution), or, if you are working with the source code, you must compile BoSSSpad and put it side-by-side to this worksheet file (from the original location in the repository, you can use the scripts getbossspad.sh, resp. getbossspad.bat).

#r "BoSSSpad.dll"
using System;
using System.Collections.Generic;
using System.Linq;
using System.Diagnostics;
using ilPSP;
using ilPSP.Utils;
using BoSSS.Platform;
using BoSSS.Platform.LinAlg;
using BoSSS.Foundation;
using BoSSS.Foundation.XDG;
using BoSSS.Foundation.Grid;
using BoSSS.Foundation.Grid.Classic;
using BoSSS.Foundation.Grid.RefElements;
using BoSSS.Foundation.IO;
using BoSSS.Solution;
using BoSSS.Solution.Control;
using BoSSS.Solution.GridImport;
using BoSSS.Solution.Statistic;
using BoSSS.Solution.Utils;
using BoSSS.Solution.AdvancedSolvers;
using BoSSS.Solution.Gnuplot;
using BoSSS.Application.BoSSSpad;
using BoSSS.Application.XNSE_Solver;
using static BoSSS.Application.BoSSSpad.BoSSSshell;
Init();

The below script needs to be able to find the current output cell; this is an easy method to get it.

Now, a new database has to be created. In this worksheet, we use a temporary database which will be deleted after the worksheet has been executed. For your calculation, you might consider some non-temporary alternative,

cf. OpenOrCreateDatabase or OpenOrCreateDefaultDatabase:

var myDb = CreateTempDatabase();

Creating database 'C:\Users\jenkinsci\AppData\Local\Temp\2016208662'.

Create a new control object for setting up the simulation:

var c = new XNSE_Control();

How to define/change input data

First, the DG polynomial degree is set: (degree 2 for velocity and 1 for pressure).

c.SetDGdegree(2);

Domain and Grid variables are set (i.e. we get a channel with length 22 and height 4.1)

double xMin        = -2; 
double xMax        = 20;
double yMin        = -2; 
double yMax        = 2.1;
int numberOfCellsX = 44; 
int numberOfCellsY = 8;

Basic database options

c.SetDatabase(myDb);
c.savetodb   = true;
c.saveperiod = 1;

Setting some variables for database saving. Here it is also possible to define tags which can be helpful for finding a particular simulation in the BoSSS database

string sessionName   = "dt = 1E20_" + numberOfCellsX + "x" + numberOfCellsY + "_k2";
c.SessionName        = sessionName;
c.ProjectDescription = "Incompressible Solver Examples";
c.Tags.Add("numberOfCellsX_" + numberOfCellsX);
c.Tags.Add("numberOfCellsY_" + numberOfCellsY);
c.Tags.Add("k2");

The grid is generated using the previously defined parameters.

c.GridFunc       = null; 
var xNodes       = GenericBlas.Linspace(xMin, xMax , numberOfCellsX);
var yNodes       = GenericBlas.Linspace(yMin, yMax, numberOfCellsY);
GridCommons grid = Grid2D.Cartesian2DGrid(xNodes, yNodes, CellType.Square_Linear, false);

Set the geometric location of boundary conditions by edge tags; Later we will assign values depending on these tags.

Edges that get assigned "0" are "inner edges".

grid.DefineEdgeTags(delegate (double[] X) {
    if (Math.Abs(X[1] - (-2)) <= 1.0e-8)
        return "wall"; // wall at y = -2
    if (Math.Abs(X[1] - (+2.1 )) <= 1.0e-8)
        return "wall"; // wall at y = +2.1
    if (Math.Abs(X[0] - (-2)) <= 1.0e-8)
        return "Velocity_Inlet"; // velocity inlet at x = -2
    if (Math.Abs(X[0] - (+20.0)) <= 1.0e-8)
        return "Pressure_Outlet"; // pressure outlet at x = +20
    throw new ArgumentException("unexpected domain boundary"); 
});

Grid Edge Tags changed.

Save the grid in the database so that the simulation can use it

grid

{ Guid = bbe9639e-9ce3-42aa-9ddc-501591c4f006; Name = ; Cell Count = 301; Dim = 2 }

myDb.SaveGrid(ref grid);
c.SetGrid(grid);

Specification of boundary conditions with a parabolic velocity profile for the inlet

c.BoundaryValues.Clear();
c.AddBoundaryValue("Velocity_Inlet", "VelocityX", 
      (X => (4.1 * 1.5 * (X[1] + 2) * (4.1 - (X[1] + 2)) / (4.1 * 4.1))));

Fluid Properties

Note: The characteristic length and fluid density are choosen to one. Therefore, the viscosity can be defined by $\frac{1}{reynolds}$.

double reynolds            = 20;
c.PhysicalParameters.rho_A = 1;
c.PhysicalParameters.mu_A  = 1.0/reynolds;

Bool parameter whether the Navier-Stokes or Stokes equations should be solved

c.PhysicalParameters.IncludeConvection = true;

Initial Values are set to 0; Note that the following lines are only for demonstration -- if no initial value is specified, 0 is set automatically.

c.InitialValues.Clear();
c.InitialValues.Add("VelocityX", new Formula("X => 0.0", false));
c.InitialValues.Add("VelocityY", new Formula("X => 0.0", false));
c.InitialValues.Add("Pressure", new Formula("X => 0.0", false));

Timestepping properties: Most solvers in BoSSS simulate transient equations. Configuring a steady simulation confiures one very large timestep.

c.TimesteppingMode = AppControl._TimesteppingMode.Steady;

Run a simulation

The solver can be run inline (i.e. within the BoSSSpad process) by executing the Run method on the control objece c. An inline run will block BoSSSpad until the solver exits.

c.Run();

Session ID: d33f1bff-a610-4726-9661-7646852a0c16, DB path: 'C:\Users\jenkinsci\AppData\Local\Temp\2016208662'

Session directory 'C:\Users\jenkinsci\AppData\Local\Temp\2016208662\sessions\d33f1bff-a610-4726-9661-7646852a0c16'.

Grid repartitioning method: METIS

Grid repartitioning options: 

Number of cell Weights: 0

Going with agglomeration threshold: 0.1

Linearization hint: GetJacobiOperator

=============== Operator Configuration ===============

     isGravity                     :[ ]

     isVolForce                    :[ ]

     isTransport                   :[x]

     isViscous                     :[x]

     isPressureGradient            :[x]

     isInterfaceSlip               :[ ]

     isContinuity                  :[x]

     isMovingMesh                  :[ ]

     isMatInt                      :[x]

     isPInterfaceSet               :[ ]

     isImmersedBoundary            :[ ]

     withPressureDissipation       :[ ]

=============== Linear Solver Configuration ===============

     Solvercode                    :Sparse direct solver PARDISO

=============== Nonlinear Solver Configuration ===============

     Solvercode                    :Newton

     Convergence Criterion         :0

     Globalization                 :Dogleg

     Minsolver Iterations          :4

     Maxsolver Iterations          :2000

======================================================

Level-Set field Phi is **exactly** zero: setting entire field to -1.

All Cells: min=301 max=301 avg=301 inb=0 tot=301

Cut Cells: min=0 max=0 avg=0 inb=0, tot=0

Starting time step 1, dt = 1.7976931348623158E+304 ...

 #Line, #Time, #Iter	L2Norm MomentumX    L2Norm MomentumY    L2Norm ContiEq      L2Norm Total        

     1,     1,     1	1.587922E+001       0.000000E+000       6.317042E+000       1.708961E+001       

     2,     1,     2	2.031857E-001       1.514689E-003       1.032251E-014       2.031914E-001       

     3,     1,     3	4.776973E-005       3.812315E-006       7.091384E-015       4.792161E-005       

     4,     1,     4	1.174523E-011       7.015082E-012       8.116256E-015       1.368071E-011       

     5,     1,     5	3.333286E-014       7.567489E-016       7.159226E-015       3.410142E-014       

     6,     1,     6	3.239210E-014       6.929550E-016       6.448747E-015       3.303506E-014       

     7,     1,     7	2.922436E-014       7.439214E-016       5.836745E-015       2.981081E-014       

Done with time step 1; solver success: True



Removing tag: NotTerminated

Postprocessing

In order to postprocess data we need to export it with the following command. This creates a folder containing the data as files (In Jupyter you need to copy the commands from the markdown into a code field to execute them)

myDb.Sessions.First().Export().Do();

You can now go to the path and Plot the data using a programm you prefer (e.g. VisIt or Paraview)

Open the ExportDirectory to view the *.plt files (does only work in BoSSSPad, but you should see the path from the commands before)

myDb.Sessions.First().OpenExportDirectory();

Some information like the console output or a log containing various physical values can be found in the session directory (does only work in BoSSSPad, but you should see the path from the commands before)

myDb.Sessions.First().OpenSessionDirectory();

Delete database

DatabaseUtils.DeleteDatabase(myDb.Path);

Immersed boundary method

It is also possible to use the immersed boundary feature of our incompressible Navier-Stokes Solver.

For this example we have to change two parts of the code: First, for a good result, we have to refine the grid at the position of the cylinder.

x-Direction (using also hyperbolic tangential distribution)

var _xNodes1 = Grid1D.TanhSpacing(-2, -1, 10, 0.5, false); 
_xNodes1     = _xNodes1.GetSubVector(0, (_xNodes1.Length - 1));
var _xNodes2 = GenericBlas.Linspace(-1, 2, 35); 
_xNodes2     = _xNodes2.GetSubVector(0, (_xNodes2.Length - 1));
var _xNodes3 = Grid1D.TanhSpacing(2, 20, 60 , 1.5, true);  
var xNodes   = ArrayTools.Cat(_xNodes1, _xNodes2, _xNodes3);

y-Direction

var _yNodes1 = Grid1D.TanhSpacing(-2, -1, 7, 0.9, false); 
_yNodes1     = _yNodes1.GetSubVector(0, (_yNodes1.Length - 1));
var _yNodes2 = GenericBlas.Linspace(-1, 1, 25); 
_yNodes2     = _yNodes2.GetSubVector(0, (_yNodes2.Length - 1));
var _yNodes3 = Grid1D.TanhSpacing(1, 2.1, 7, 1.1, true);  
var yNodes   = ArrayTools.Cat(_yNodes1, _yNodes2, _yNodes3);

Furthermore, the cylinder immersing the fluid should be described by using the zero contour of a level set function. The radius of the cylinder is set to 0.5.

c.InitialValues.Add("Phi", new Formula("X => -(X[0]).Pow2() + -(X[1]).Pow2() + 0.25", false));

Example control files for both, the channel and the flow around a cylinder can be found in the ControlExample directory. As soon as we run the simulation again we can take a look at the plots and the PhysicalData file in the session directory. There we can find for example lift and drag forces acting on the cylinder.