Demo problem: Flow in a 2D channel with an oscillating leaflet

In this example we consider the flow in a 2D channel which is partially obstructed by an oscillating leaflet. We consider the case where the motion of the leaflet is prescribed – this is a "warm-up exercise" for the corresponding FSI problem in which the leaflet is an elastic structure.

The Problem

The figure below shows a sketch of the problem: A 2D channel of height $ H^*_{tot}$ and length $ L^*_{left} + L^*_{right} $ is partially occluded by a (zero-thickness) leaflet of height $ H^*_{leaflet} $. The leaflet is parametrised by a Lagrangian coordinate $ \xi^* $ so that the position vector to a material point on the leaflet is given by $ {\bf R}^*_w(\xi^*,t^*) $, and we assume that the leaflet performs time-periodic oscillations with period $ T^*. $ Steady Poiseuille flow with average velocity $ U^* $ is imposed at the left end of the channel while we assume that the outflow is parallel and axially traction-free.

Sketch of the problem in dimensional

We non-dimensionalise all length and coordinates on the channel width, $ H^*_{tot} $ , time on the natural timescale of the flow, $ H^*_{tot}/U^* $, the velocities on the mean velocity, $ U^* $, and the pressure on the viscous scale. The problem is then governed by the non-dimensional Navier-Stokes equations

\[ Re \left( St \frac{\partial u_i}{\partial t} + u_j \frac{\partial u_i}{\partial x_j} \right) = - \frac{\partial p}{\partial x_i} + \frac{\partial }{\partial x_j} \left[ \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \right], \]

where $ Re = \rho U^* H_0^* / \mu $ and $ St = 1 $, and

\[ \frac{\partial u_i}{\partial x_i} = 0, \]

subject to parabolic inflow

\[ {\bf u} = 6 x_2 (1-x_2) {\bf e}_1 \]

at the inflow cross-section; parallel, axially-traction-free outflow at the outlet; and no-slip on the stationary channel walls, $ {\bf u} = {\bf 0} $. The no-slip condition on the leaflet is

\[ {\bf u} = \frac{\partial {\bf R}_w(\xi,t)}{\partial t} \]

and the leaflet performs oscillations with non-dimensional period $ T = T^* U^* / H^*_{tot} $. Here is a sketch of the non-dimensional version of the problem:

Sketch of the problem in dimensionless variables.

An interesting feature of this problem is that even though the leaflet is assumed to have negligible thickness its presence will generate a pressure jump between its two faces. (The velocities are continuous because the no-slip condition imposes the same velocity on both faces.) When discretising the problem with QTaylorHood elements (for which the pressure varies continuously across the element boundaries), the mesh must therefore be "opened up" with a cut along the position of the leaflet. This is done in the constructor of the ChannelWithLeafletMesh which forms the basis of RefineableAlgebraicChannelWithLeafletMesh used to discretise this problem. (See Further comments and exercises for a more detailed discussion of the mesh.)


The figure below shows a snapshot of the flow field (pressure contours and instantaneous streamlines) for a Reynolds number of $ Re=20 $ and an oscillation period of $ T=20 $, as well as the corresponding fluid mesh. Note how oomph-lib's automatic mesh adaptation has refined the mesh near the tip of the leaflet where the pressure has a singularity.

Mesh (top) and flow field (bottom).

The corresponding animation illustrates the algebraic node update strategy (implemented with an AlgebraicMesh, discussed in more detail in another tutorial) and the evolution of the flow field. Note that the instantaneous streamlines intersect the (impermeable) leaflet because the leaflet is not stationary.

The global parameters

As usual we use a namespace to define the (single) global parameter, the Reynolds number.

/// Global parameters
namespace Global_Physical_Variables
/// Reynolds number
double Re=20.0;
} // end_of_namespace

Specification of the leaflet geometry

We specify the leaflet geometry and its time-dependent motion by representing it as a GeomObject. The GeomObject has one Lagrangian and two Eulerian coordinates, and its geometry is characterised by its length, the x-coordinate of its origin, $ X_0 $, and the period and amplitude of the horizontal and vertical tip deflection.

/// GeomObject representing a vertical leaflet that performs
/// bending and stretching oscillations.
class Leaflet : public GeomObject
/// \short Constructor: Pass length (in Lagrangian coordinates),
/// the amplitude of the horizontal and vertical deflection of the tip,
/// the x-coordinate of the origin and the period of the oscillation.
/// Passes the number of Lagrangian and Eulerian coordinates to the
/// constructor of the GeomObject base class.
Leaflet(const double& length, const double& d_x,const double& d_y,
const double& x_0, const double& period, Time* time_pt)
: GeomObject(1,2), Length(length), D_x(d_x), D_y(d_y), X_0(x_0),
T(period),Time_pt(time_pt) {}
/// Destructor -- emtpy
virtual ~Leaflet(){}
/// \short Position vector, r, to the point identified by
/// its 1D Lagrangian coordinate, xi (passed as a 1D Vector) at discrete time
/// level t (t=0: present; t>0: previous).
void position(const unsigned& t, const Vector<double>& xi,
Vector<double>& r) const
using namespace MathematicalConstants;
r[0] = X_0 + D_x*xi[0]*xi[0]/Length/Length*sin(2.0*Pi*Time_pt->time(t)/T);
r[1] = xi[0]*(1.0+D_y/Length*0.5*(1.0-cos(4.0*Pi*Time_pt->time(t)/T)));
/// Steady version: Get current shape
void position(const Vector<double>& xi, Vector<double>& r) const
/// Number of geometric Data in GeomObject: None.
unsigned ngeom_data() const {return 0;}
/// Length of the leaflet
double length() { return Length; }
/// Amplitude of horizontal tip displacement
double& d_x() {return D_x;}
/// Amplitude of vertical tip displacement
double d_y() {return D_y;}
/// x-coordinate of leaflet origin
double x_0() {return X_0;}
private :
/// Length in terms of Lagrangian coordinates
double Length;
///\short Horizontal displacement of tip
double D_x;
///\short Vertical displacement of tip
double D_y;
///\short Origin
double X_0;
///Period of the oscillations
double T;
///Pointer to the global time object
Time* Time_pt;
}; //end_of_the_GeomObject

The driver code

We store the command line arguments, create a DocInfo object, and assign the parameters that specify the domain and the leaflet geometry:

/// Driver code -- pass a command line argument if you want to run
/// the code in validation mode where it only performs a few steps
int main(int argc, char* argv[])
// Store command line arguments
// Set up doc info
DocInfo doc_info;
// Parameters for the leaflet
// Height
double h_leaflet = 0.5;
// Tip deflection
double d_x = 0.25;
double d_y = -0.05;
// x-positon of root
double x_0 = 3.0;
// Period of the oscillation on the natural timescale of the flow
double period = 20.0;
//Parameters for the domain
// Length of the mesh to right and left of the leaflet
double l_left =2.0;
double l_right= 3.0;
// Total height of domain (unity because lengths have been scaled on it)
double h_tot=1.0;
// Initial number of element rows/columns in various mesh regions
unsigned nleft=8;
unsigned nright=12;
unsigned ny1=2;
unsigned ny2=2;

Next we build the problem and assign the time-stepping parameters (as usual, fewer timesteps are used during a validation run which is identified by a non-zero number of command line arguments):

//Build the problem
problem(l_left, l_right, h_leaflet,
h_tot,nleft, nright,ny1,ny2,
d_x, d_y, x_0,
// Number of timesteps per period
unsigned nsteps_per_period=40;
// Number of periods
unsigned nperiod=3;
// Number of timesteps (reduced for validation)
unsigned nstep=nsteps_per_period*nperiod;
if (CommandLineArgs::Argc>1)
double dt=period/double(nsteps_per_period);
/// Initialise timestep

We start the simulation with a steady solve, allowing up to five levels of adaptive refinement (fewer if we are performing a validation run):

/// Set max. number of adaptations (reduced for validation)
unsigned max_adapt=5;
if (CommandLineArgs::Argc>1)
// Do steady solve first -- this also sets the history values
// to those corresponding to an impulsive start from the
// steady solution
/// Output steady solution

Finally, we enter the proper timestepping loop, allowing one spatial adaptation per timestep and suppressing the re-assignment of initial conditions following an adaptation by setting the parameter first to false (see the discussion of timestepping with automatic mesh adaptation in another tutorial.)

/// Reduce the max number of adaptations for time-dependent simulation
// We don't want to re-assign the initial condition
bool first=false;
// Timestepping loop
for (unsigned istep=0;istep<nstep;istep++)
// Solve the problem
// Output the solution
// Step number
}//end of main

The Problem class

The problem class has the usual member functions to perform actions after the mesh adaptation and before every implicit timestep:

/// Problem class
template<class ELEMENT>
class ChannelWithLeafletProblem : public Problem
/// Constructor: Pass the length of the domain at the left
/// of the leaflet lleft,the length of the domain at the right of the
/// leaflet lright,the height of the leaflet hleaflet, the total height
/// of the domain htot, the number of macro-elements at the left of the
/// leaflet nleft, the number of macro-elements at the right of the
/// leaflet nright, the number of macro-elements under hleaflet ny1,
/// the number of macro-elements above hleaflet ny2,the x-displacement
/// of the leaflet d_x,the y-displacement of the leaflet d_y,the abscissa
/// of the origin of the leaflet x_0, the period of the moving leaflet.
ChannelWithLeafletProblem(const double& l_left,
const double& l_right, const double& h_leaflet,
const double& h_tot,
const unsigned& nleft, const unsigned& nright,
const unsigned& ny1, const unsigned& ny2,
const double& d_x,const double& d_y,
const double& x_0, const double& period);
/// Destructor (empty)
/// Overloaded access function to specific mesh
RefineableAlgebraicChannelWithLeafletMesh<ELEMENT>* mesh_pt()
// Upcast from pointer to the Mesh base class to the specific
// element type that we're using here.
return dynamic_cast<RefineableAlgebraicChannelWithLeafletMesh<ELEMENT>*>(
/// Update after solve (empty)
/// Update before solve (empty)
/// Actions after adaptation: Pin redundant pressure dofs
/// Update the velocity boundary condition on the moving leaflet
/// Doc the solution
void doc_solution(DocInfo& doc_info);
/// Pointer to the GeomObject
GeomObject* Leaflet_pt;

The problem constructor

We construct the timestepper and the GeomObject that represents the leaflet, and pass it pointers to them to the constructor of the RefineableAlgebraicChannelWithLeafletMesh (discussed in more detail in Further comments and exercises).

/// Constructor
template <class ELEMENT>
const double& l_left,
const double& l_right, const double& h_leaflet,
const double& h_tot,
const unsigned& nleft, const unsigned& nright,
const unsigned& ny1, const unsigned& ny2,
const double& d_x,const double& d_y,
const double& x_0, const double& period)
// Allocate the timestepper
add_time_stepper_pt(new BDF<2>);
//Create the geometric object that represents the leaflet
Leaflet_pt = new Leaflet(h_leaflet, d_x, d_y, x_0, period, time_pt()) ;
//Build the mesh
Problem::mesh_pt()=new RefineableAlgebraicChannelWithLeafletMesh<ELEMENT>(
l_left, l_right,

Next we create the spatial error estimator and loop over the elements to set the pointers to the relevant physical parameters.

// Set error estimator
Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;
// Loop over the elements to set up element-specific
// things that cannot be handled by constructor
unsigned n_element = mesh_pt()->nelement();
for(unsigned e=0;e<n_element;e++)
// Upcast from GeneralisedElement to the present element
ELEMENT* el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
//Set the Reynolds number
el_pt->re_pt() = &Global_Physical_Variables::Re;
// Set the Womersley number (product of Reynolds and Strouhal).
// We're assuming a Strouhal number of one, corresponding to
// a non-dimensionalisation of time on the flow's natural timescale.
el_pt->re_st_pt() = &Global_Physical_Variables::Re;
} // end loop over elements

The velocity is prescribed everywhere apart from the outflow boundary (boundary 1; see the sketch in Further comments and exercises for the enumeration of the mesh boundaries). Along the inflow (boundary 3) we apply a parabolic velocity profile with unit flux:

//Pin the boundary nodes
unsigned num_bound = mesh_pt()->nboundary();
for(unsigned ibound=0;ibound<num_bound;ibound++)
unsigned num_nod= mesh_pt()->nboundary_node(ibound);
for (unsigned inod=0;inod<num_nod;inod++)
//do not pin the x velocity of the outflow
if( ibound != 1)
// Setup parabolic flow along the inflow boundary 3
unsigned ibound=3;
unsigned num_nod= mesh_pt()->nboundary_node(ibound);
for (unsigned inod=0;inod<num_nod;inod++)
double ycoord = mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
double uy = 6.0*(ycoord/h_tot)*(1.0-(ycoord/h_tot));
}// end of setup boundary condition

Finally, we pin the redundant pressure degrees of freedom (see another tutorial for details), and assign the equations numbers.

// Pin redudant pressure dofs
// Setup equation numbering scheme
cout << "Number of equations: " << assign_eqn_numbers() << std::endl;
}//end of constructor

Actions before the timestep

Before each timestep we update the nodal positions in the mesh and re-apply the no-slip condition on the nodes of the moving leaflet (boundaries 4 and 5; see the sketch in Further comments and exercises for the enumeration of the mesh boundaries).

/// Actions before implicit timestep: Update domain shape and
/// the velocity boundary conditions
template <class ELEMENT>
// Update the domain shape
// Moving leaflet: No slip; this implies that the velocity needs
// to be updated in response to leaflet motion
for( unsigned ibound=4;ibound<6;ibound++)
unsigned num_nod=mesh_pt()->nboundary_node(ibound);
for (unsigned inod=0;inod<num_nod;inod++)
// Which node are we dealing with?
Node* node_pt=mesh_pt()->boundary_node_pt(ibound,inod);
// Apply no slip
} //end_of_actions_before_implicit_timestep

Actions after the mesh adaptation

Once the mesh has been adapted, we free all pressure degrees of freedom and then (re-)pin any redundant ones (see another tutorial for details):

// Actions after adaptation: Pin redundant pressure dofs
template<class ELEMENT>
// Unpin all pressure dofs
// Pin redundant pressure dofs
} // end_of_actions_after_adapt

Note that the default interpolation of the (quadratic!) inflow velocity profile from father to son elements during the mesh adaptation already ensures that the inflow profile remains quadratic, therefore no further action is required.


The function doc_solution(...) documents the results.

/// Doc the solution
template<class ELEMENT>
ofstream some_file;
char filename[100];
// Number of plot points
unsigned npts;
// Output solution
// Output boundaries
} // end_of_doc_solution

Further comments and exercises

Further comments: The algebraic node update procedure

The figure below illustrates the algebraic node update procedure employed in the RefineableAlgebraicChannelWithLeafletMesh. The mesh employs four different node update functions, depending on which region a node is located in: Nodes in region I (or II) are located on straight lines that connect the upstream (or downstream) boundary with the leaflet; nodes in region III (or IV) are located on straight lines that connect the upstream and (or downstream) boundary with the straight line from the tip of the leaflet to upper channel wall.

Sketch illustrating the algebraic node update procedure.

The implementation of the node update functions is straightforward and can be found in the source files




which also illustrate how the mesh is constructed by inheritance from the SimpleRectangularQuadMesh (the main task being the creation of additional nodes in the interior to "cut open" the mesh along the position of the leaflet). The source files also contain other versions of the mesh in which the node update is performed with Domain/MacroElements, using the technique described in another tutorial.


With the node update strategy illustrated above, the position of all nodes in the fluid mesh has to be updated when the leaflet moves. This is not a particular problem in the current application where the node-update is only performed once per timestep. However, in the corresponding FSI problem , the approach is costly because of the of large number of shape derivatives to be computed.

As an exercise, we suggest to make the node-update procedure more efficient by sub-dividing the regions upstream and downstream of the leaflet into a central section in which the nodes move in response to the motion of the leaflet (the old regions I-IV) and two additional regions (regions V and VI) in which they remain stationary. This is easy because, as explained elsewhere, all AlgebraicNodes already have a default node update function that leaves them stationary.

A better node update strategy.


Source files for this tutorial

PDF file

A pdf version of this document is available.