Example problem: Flow in a 2D channel with an oscillating wall

In this example we consider our first time-dependent Navier-Stokes problem and demonstrate how to apply periodic boundary conditions.

The example problem

We consider finite-Reynolds-number flow in a 2D channel that is driven by the oscillatory tangential motion of the "upper" wall:

Unsteady flow in a 2D channel with an oscillating wall.

Here is a sketch of the problem:

Sketch of the problem.

The flow is governed by the 2D unsteady 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( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) \ \ \ \ \ \ \ \ \ \ (1) \]

and the continuity equation

\[ \frac{\partial u_i}{\partial x_i} = 0 \ \ \ \ \ \ \ \ \ \ (2) \]

in the domain

\[ D = \bigg\{(x_1,x_2) \ | \ x_1 \in [0,1], \ x_2 \in [0,1] \bigg\}. \]

We apply the Dirichlet (no-slip) boundary condition

\[ \left. \mathbf{u}\right|_{\partial D_{lower}}=(0,0), \ \ \ \ \ \ \ \ \ \ (3) \]

on the lower, stationary wall, $ \partial D_{lower} = \big\{(x_1,x_2) \ | \ x_2=0 \big\} $ , apply the Dirichlet (no-slip) conditions

\[ \left. \mathbf{u}\right|_{\partial D{upper}}=(\sin\left(\omega t\right),0), \ \ \ \ \ \ \ \ \ \ (4) \]

on the upper, moving wall , $ \partial D_{upper} = \big\{(x_1,x_2) \ | \ x_2=1 \big\} $ , and apply periodic boundary condition on the "left" and "right" boundaries:

\[ \left.\mathbf{u}\right|_{x_1=0} = \left.\mathbf{u}\right|_{x_1=1}. \ \ \ \ \ \ \ \ \ \ (5) \]

Initial conditions for the velocities are given by

\[ {\bf u}(x_1,x_2,t=0) = {\bf u}_{IC}(x_1,x_2), \]

where $ {\bf u}_{IC}(x_1,x_2)$ is given.

The exact solution

The above problem has an exact, time-periodic parallel flow solution of the form

\[ {\bf u}_{exact}(x_1,x_2,t) = U(x_2,t) \ {\bf e}_1 \mbox{ \ \ \ \ \ and \ \ \ \ } p_{exact}(x_1,x_2,t) = 0, \]

where $ U(x_2,t) $ is governed by

\[ Re\, St\frac{\partial U}{\partial t} = \frac{\partial^2 U}{\partial x_2^2}, \]

subject to the boundary conditions $ U(x_2=0,t)=0$ and $ U(x_2=1,t)=\sin(\omega t)$ . The solution is given by

\[ U(x_2,t) = Re\left\{\frac{e^{i\omega t}}{e^{i\lambda}-e^{-i\lambda}} \left(e^{i\lambda y}-e^{-i\lambda y}\right)\right\}, \]


\[ \lambda = i\sqrt{i\omega Re\, St}. \]


The two animations below show the computed solutions obtained from a spatial discretisation with Taylor-Hood and Crouzeix-Raviart elements, respectively. In both cases we set $ \omega=2\pi, \ Re = ReSt = 10 $ and specified the exact, time-periodic solution as the initial condition, i.e. $ {\bf u}_{IC}(x_1,x_2) = {\bf u}_{exact}(x_1,x_2,t=0)$ . The computed solutions agree extremely well with the exact solution throughout the simulation.

Plot of the velocity field computed with 2D Crouzeix-Raviart elements, starting form the exact, time-periodic solution.
Plot of the velocity field computed with 2D Taylor-Hood elements, starting form the exact, time-periodic solution.

If the simulation is started from other initial conditions, i.e. $ {\bf u}_{IC}(x_1,x_2) \ne {\bf u}_{exact}(x_1,x_2,t=0)$ , the velocity field initially differs noticeably from the time-periodic solution $ {\bf u}_{exact}(x_1,x_2,t)$ but following the decay of initial transients we have

\[ \lim_{t\to \infty} {\bf u}(x_1,x_2,t) = {\bf u}_{exact}(x_1,x_2,t). \]

This is illustrated in the following plot which shows the evolution of the L2-"error" between the computed and the time-periodic solutions for two different initial conditions. The red line was obtained from a simulation in which $ {\bf u}_{IC}(x_1,x_2) = {\bf u}_{exact}(x_1,x_2,t=0)$ ; the blue line was obtained from a computation in which the simulation was started by an "impulsive start", $ {\bf u}_{IC}(x_1,x_2) = {\bf 0}$ .

Plot of the L2-'errors' between the computed and time-periodic solution for two different initial conditions.

The animations of the simulations for the "impulsive start" (for Taylor-Hood and Crouzeix-Raviart elements) show how the velocity profile approaches the time-periodic solution as the simulation progresses.

The global parameters

As usual, we use a namespace to define the problem parameters, the Reynolds number, $ Re$ , and the Womersley number, $ Re\, St$ . We also provide two flags that indicate the length of the run (to allow a short run to be performed when the code is run as a self-test), and the initial condition (allowing a start from $ {\bf u}_{exact}$ or an impulsive start in which the fluid is initially at rest).

/// Namespace for global parameters
/// Reynolds number
double Re;
/// Womersley = Reynolds times Strouhal
double ReSt;
/// Flag for long/short run: Default = perform long run
unsigned Long_run_flag=1;
/// Flag for impulsive start: Default = start from exact
/// time-periodic solution.
} // end of namespace
Namespace for global parameters.
unsigned Long_run_flag
Flag for long/short run: Default = perform long run.
double ReSt
Womersley = Reynolds times Strouhal.
double Re
Reynolds number.
unsigned Impulsive_start_flag
Flag for impulsive start: Default = start from exact time-periodic solution.

The exact solution

We use another namespace to define the exact, time-periodic parallel-flow solution:

/// Namespace for exact solution
namespace ExactSoln
/// Exact solution of the problem as a vector
void get_exact_u(const double& t, const Vector<double>& x, Vector<double>& u)
double y=x[1];
// I=sqrt(-1)
complex<double> I(0.0,1.0);
// omega
double omega=2.0*MathematicalConstants::Pi;
// lambda
complex<double> lambda(0.0,omega*Global_Parameters::ReSt);
lambda = I*sqrt(lambda);
// Solution
complex<double> sol(
exp(complex<double>(0.0,omega*t)) *
/(exp(I*lambda)-exp(-I*lambda)) );
// Assign real solution
/// Exact solution of the problem as a scalar
void get_exact_u(const double& t, const double& y,double& u)
// I=sqrt(-1)
complex<double> I(0.0,1.0);
// omega
double omega=2.0*MathematicalConstants::Pi;
// lambda
complex<double> lambda(0.0,omega*Global_Parameters::ReSt);
lambda = I*sqrt(lambda);
// Solution
complex<double> sol(
exp(complex<double>(0.0,omega*t)) *
/(exp(I*lambda)-exp(-I*lambda)) );
// Assign real solution
} // end of exact_solution
Namespace for exact solution.
void get_exact_u(const double &t, const Vector< double > &x, Vector< double > &u)
Exact solution of the problem as a vector.

The driver code

We use optional command line arguments to specify which mode the code is run in: Either as a short or a long run (indicated by the first command line argument being 0 or 1, respectively), and with initial conditions corresponding to an impulsive start or a start from the time-periodic exact solution (indicated by the second command line argument being 1 or 0, respectively). If no command line arguments are specified the code is run in the default mode, specified by the parameter values assigned in the namespace Global_Parameters.

/// Driver code for Rayleigh channel problem
int main(int argc, char* argv[])
/// Convert command line arguments (if any) into flags:
if (argc==1)
cout << "No command line arguments specified -- using defaults."
<< std::endl;
else if (argc==3)
cout << "Two command line arguments specified:" << std::endl;
// Flag for long run
/// Flag for impulsive start
std::string error_message =
"Wrong number of command line arguments. Specify none or two.\n";
error_message +=
"Arg1: Long_run_flag [0/1]\n";
error_message +=
"Arg2: Impulsive_start_flag [0/1]\n";
throw OomphLibError(error_message,
cout << "Long run flag: "
cout << "Impulsive start flag: "
int main(int argc, char *argv[])
Driver code for Rayleigh channel problem.

Next we set the physical and mesh parameters.

// Set physical parameters:
// Womersley number = Reynolds number (St = 1)
//Horizontal length of domain
double lx = 1.0;
//Vertical length of domain
double ly = 1.0;
// Number of elements in x-direction
unsigned nx=5;
// Number of elements in y-direction
unsigned ny=10;

Finally we set up DocInfo objects and solve for both Taylor-Hood elements and Crouzeix-Raviart elements.

// Solve with Crouzeix-Raviart elements
// Set up doc info
DocInfo doc_info;
//Set up problem
RayleighProblem<QCrouzeixRaviartElement<2>,BDF<2> > problem(nx,ny,lx,ly);
// Run the unsteady simulation
// Solve with Taylor-Hood elements
// Set up doc info
DocInfo doc_info;
//Set up problem
RayleighProblem<QTaylorHoodElement<2>,BDF<2> > problem(nx,ny,lx,ly);
// Run the unsteady simulation
} // end of main
Rayleigh-type problem: 2D channel whose upper wall oscillates periodically.

The problem class

The problem class is very similar to that used in the driven cavity example. We specify the type of the element and the type of the timestepper (assumed to be a member of the BDF family) as template parameters and pass the mesh parameters to the problem constructor.

/// Rayleigh-type problem: 2D channel whose upper
/// wall oscillates periodically.
template<class ELEMENT, class TIMESTEPPER>
class RayleighProblem : public Problem
/// Constructor: Pass number of elements in x and y directions and
/// lengths
RayleighProblem(const unsigned &nx, const unsigned &ny,
const double &lx, const double &ly);
RayleighProblem(const unsigned &nx, const unsigned &ny, const double &lx, const double &ly)
Constructor: Pass number of elements in x and y directions and lengths.

No action is needed before or after solving, but we update the time-dependent boundary conditions at the upper wall before each timestep, using Problem::actions_before_implicit_timestep(). The boundary values are obtained from the exact solution, defined in the namespace ExactSoln.

//Update before solve is empty
void actions_before_newton_solve() {}
/// Update after solve is empty
void actions_after_newton_solve() {}
//Actions before timestep: Update no slip on upper oscillating wall
void actions_before_implicit_timestep()
// No slip on upper boundary
unsigned ibound=2;
unsigned num_nod=mesh_pt()->nboundary_node(ibound);
for (unsigned inod=0;inod<num_nod;inod++)
// Get exact solution
double y=mesh_pt()->boundary_node_pt(ibound,inod)->x(1);
double time=time_pt()->time();
double soln;
// Assign exact solution to boundary
} // end of actions_before_implicit_timestep

The function unsteady_run(...), discussed below, performs the timestepping and documents the solution in the directory specified in the DocInfo object.

/// Run an unsteady simulation
void unsteady_run(DocInfo& doc_info);

We define the function doc_solution(...) which documents the results, and provide functions to set the initial conditions and to fix a pressure value. The problem's only member data contains an output stream in which we record the time-trace of the solution.

/// Doc the solution
void doc_solution(DocInfo& doc_info);
/// Set initial condition (incl previous timesteps) according
/// to specified function.
void set_initial_condition();
/// Fix pressure in element e at pressure dof pdof and set to pvalue
void fix_pressure(const unsigned &e, const unsigned &pdof,
const double &pvalue)
//Cast to proper element and fix pressure
/// Trace file
ofstream Trace_file;
}; // end of problem_class

The problem constructor

We start by building the timestepper, determining its type from the class's second template argument, and pass a pointer to it to the Problem, using the function Problem::add_time_stepper_pt(...).

/// Problem constructor
template<class ELEMENT,class TIMESTEPPER>
(const unsigned &nx, const unsigned &ny,
const double &lx, const double& ly)
//Allocate the timestepper
add_time_stepper_pt(new TIMESTEPPER);

Next we build the mesh and pass an additional boolean flag to the constructor to indicate that periodic boundary conditions are to be applied in the $ x_1$ -direction. We will discuss the implementation of this feature in more detail in below.

//Now create the mesh with periodic boundary conditions in x direction
bool periodic_in_x=true;
Problem::mesh_pt() =
new RectangularQuadMesh<ELEMENT>(nx,ny,lx,ly,periodic_in_x,

We pin both velocity components on the top and bottom boundaries (i.e. at $ x_2=0$ and $ x_2=1$ , respectively), and pin the vertical velocity on the left and right boundaries (i.e. at $ x_1=0 $ and $ x_1=1 $ , respectively) to enforce horizontal outflow through the periodic boundaries.

// Set the boundary conditions for this problem: All nodes are
// free by default -- just pin the ones that have Dirichlet conditions
// here
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++)
// No slip on top and bottom
if ((ibound==0)||(ibound==2))
// Horizontal outflow on the left (and right -- right bc not
// strictly necessary because of symmetry)
else if ((ibound==1)||(ibound==3))
} // end loop over boundaries

Finally we pass the pointers to the Reynolds and Strouhal numbers, $ Re $ and $ Re\, St $ , to the elements. Since no traction boundary conditions are applied anywhere, the pressure is only determined up to an arbitrary constant. To ensure a unique solution we pin a single pressure value before setting up the equation numbering scheme.

//Complete the problem setup to make the elements fully functional
//Loop over the elements
unsigned n_el = mesh_pt()->nelement();
for(unsigned e=0;e<n_el;e++)
//Cast to a fluid element
ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(e));
//Set the Reynolds number, etc
el_pt->re_pt() = &Global_Parameters::Re;
el_pt->re_st_pt() = &Global_Parameters::ReSt;
// Now pin the pressure in first element at value 0 to 0.0
//Assgn equation numbers
cout << assign_eqn_numbers() << std::endl;
} // end of constructor

Initial conditions

The application of initial conditions for vector-valued problems is performed by the same procedure that we described for scalar problems, except that we now have to assign "history values" for multiple nodal values. For timesteppers from the BDF family, the "history values" represent the solution at previous timesteps. We check that the timestepper is of the appropriate type,
loop over previous time levels, determine the velocity at those times and assign them to the "history values" of the velocities. No initial conditions are required for the pressure. Note that we also have to assign "history values" for the nodal positions since oomph-lib's Navier-Stokes elements discretise the momentum equations in their ALE form. This aspect was explained in more detail in our discussion of the solution of the unsteady heat equation.

/// Set initial condition: Assign previous and current values
/// from exact solution.
template<class ELEMENT,class TIMESTEPPER>
// Check that timestepper is from the BDF family
if (time_stepper_pt()->type()!="BDF")
std::ostringstream error_stream;
error_stream << "Timestepper has to be from the BDF family!\n"
<< "You have specified a timestepper from the "
<< time_stepper_pt()->type() << " family" << std::endl;
throw OomphLibError(error_stream.str(),
// Backup time in global Time object
double backed_up_time=time_pt()->time();
// Past history needs to be established for t=time0-deltat, ...
// Then provide current values (at t=time0) which will also form
// the initial guess for the first solve at t=time0+deltat
// Vector of exact solution value
Vector<double> soln(2);
Vector<double> x(2);
//Find number of nodes in mesh
unsigned num_nod = mesh_pt()->nnode();
// Set continuous times at previous timesteps:
// How many previous timesteps does the timestepper use?
int nprev_steps=time_stepper_pt()->nprev_values();
Vector<double> prev_time(nprev_steps+1);
for (int t=nprev_steps;t>=0;t--)
// Loop over current & previous timesteps
for (int t=nprev_steps;t>=0;t--)
// Continuous time
double time=prev_time[t];
cout << "setting IC at time =" << time << std::endl;
// Loop over the nodes to set initial guess everywhere
for (unsigned n=0;n<num_nod;n++)
// Get nodal coordinates
// Get exact solution at previous time
// Assign solution
// Loop over coordinate directions: Mesh doesn't move, so
// previous position = present position
for (unsigned i=0;i<2;i++)
// Reset backed up time for global timestepper
} // end of set_initial_condition
void set_initial_condition()
Set initial condition (incl previous timesteps) according to specified function.

Post processing

The function doc_solution(...) is similar to those used in the unsteady heat examples. We plot the computed solution, the time-periodic exact solution and the difference between the two, and record various parameters in the trace file. The plot of the computed solution contains tecplot instructions that generate a blue line in the top-left corner of the plot to indicate how time progresses during the simulation. The trace file contains a record of

  • the value of the continuous time, $ t $ ,
  • the coordinates of a control node, $ (x_1^{[c]},x_2^{[c]})$ ,
  • the computed velocity at the control node, $( u_1^{[c]},u_2^{[c]})$ ,
  • the time-periodic solution, evaluated at the control node, $ (u_1^{[c,exact]},u_2^{[c,exact]})$ ,
  • the difference between the computed velocities and the time-periodic solution at the control node,
  • the L2 norm of the "error" between the computed and time-periodic solution for the velocity, and
  • the L2 norm of the time-periodic solution for the velocity.
/// Doc the solution
template<class ELEMENT,class TIMESTEPPER>
ofstream some_file;
char filename[100];
// Number of plot points
unsigned npts=5;
// Output solution
// Write file as a tecplot text object
some_file << "TEXT X=2.5,Y=93.6,F=HELV,HU=POINT,C=BLUE,H=26,T=\"time = "
<< time_pt()->time() << "\"";
// ...and draw a horizontal line whose length is proportional
// to the elapsed time
some_file << "GEOMETRY X=2.5,Y=98,T=LINE,C=BLUE,LT=0.4" << std::endl;
some_file << "1" << std::endl;
some_file << "2" << std::endl;
some_file << " 0 0" << std::endl;
some_file << time_pt()->time()*20.0 << " 0" << std::endl;
// Output exact solution
// Doc error
double error,norm;
// Doc solution and error
cout << "error: " << error << std::endl;
cout << "norm : " << norm << std::endl << std::endl;
// Get time, position and exact soln at control node
unsigned n_control=37;
Vector<double> x(2), u(2);
double time=time_pt()->time();
Node* node_pt=
x[0] = node_pt->x(0);
x[1] = node_pt->x(1);
// Write trace file
Trace_file << time << " "
<< x[0] << " "
<< x[1] << " "
<< node_pt->value(0) << " "
<< node_pt->value(1) << " "
<< u[0] << " "
<< u[1] << " "
<< abs(u[0]-node_pt->value(0)) << " "
<< abs(u[1]-node_pt->value(1)) << " "
<< error << " "
<< norm << " "
<< std::endl;
} // end_of_doc_solution
void doc_solution(DocInfo &doc_info)
Doc the solution.

The timestepping loop

The function unsteady_run(...) is used to perform the timestepping procedure. We start by opening the trace file and write a suitable header for the visualisation with tecplot.

/// Unsteady run...
template<class ELEMENT,class TIMESTEPPER>
// Open trace file
char filename[100];
// Write tecplot header for trace file
Trace_file << "time" << ", "
<< "x" << ", "
<< "y" << ", "
<< "u_1" << ", "
<< "u_2" << ", "
<< "u_exact_1" << ", "
<< "u_exact_2" << ", "
<< "error_1" << ", "
<< "error_2" << ", "
<< "L2 error" << ", "
<< "L2 norm" << ", " << std::endl;
void unsteady_run(DocInfo &doc_info)
Run an unsteady simulation.

Next, we choose a value for the timestep and set up the initial conditions, either for an impulsive start...

//Set value of dt
double dt = 0.025;
// Initialise all history values for an impulsive start
cout << "IC = impulsive start" << std::endl;

...or for a "smooth" start from the time-periodic exact solution:

// Initialise timestep
// Set initial conditions.
cout << "IC = exact solution" << std::endl;

We choose the number of timesteps to be computed and document the initial conditions.

//Now do many timesteps
unsigned ntsteps=80;
// If validation run only do 5 timesteps
cout << "validation run" << std::endl;
// Doc initial condition
// increment counter

Finally, perform the actual timestepping and document the solution after every timestep.

//Loop over the timesteps
for(unsigned t=1;t<=ntsteps;t++)
cout << "TIMESTEP " << t << std::endl;
//Take one fixed timestep
//Output the time
cout << "Time is now " << time_pt()->time() << std::endl;
// Doc solution
// increment counter
} // end of unsteady run

Comments and Exercises

Periodic boundaries

A key feature of the current problem is the presence of periodic boundary conditions. The application of the periodic boundary condition is performed "inside" the mesh constructor and details of the implementation were therefore "hidden". We will now discuss the steps required to apply periodic boundary conditions and explain why it is easier to apply periodic boundary conditions in the mesh constructor rather than in the "driver code".

Periodic boundary conditions arise in problems that are periodic in one of their coordinate directions. It is important to realise that, even though the solution at the corresponding nodes on the two periodic domain boundaries (the left and the right boundary in the above example) are the same, one of their nodal coordinates differs. For instance, in the above example, each of the nodes on the left boundary has the same velocity values and the same $ x_2$ -coordinate as its corresponding (periodic) node on the right boundary. However, the $ x_1 $ -coordinate of the nodes on the left boundary is $ x_1=0$ , whereas that of the (periodic) nodes on the right boundary is $ x_1=1$ . It is therefore not possible to regard the nodes as identical.

In oomph-lib we create periodic nodes by allowing two (or more) nodes to access some of the same internal data. One of the nodes should be regarded as the original and the other(s) are set to be its "periodic counterpart(s)" and hence access its internal data. The "periodic counterpart(s)" are created by calling the member function

BoundaryNode::make_periodic(Node *const& node_pt)

where the pointer to the original node is passed as the argument. Note that the required functionality imposes a slight storage overhead and so in oomph-lib we only allow BoundaryNodes to be made periodic.

Here is a sketch of a 2D rectangular quad mesh. If this mesh is to be used in a problem with periodic boundary conditions in the horizontal direction (as in the above example), the pointer to node 3 on boundary 3 would have to be used when node 3 on boundary 1 is made periodic, etc. The appropriate commands are

//Get pointers to the two nodes that are to be "connected" via
//a periodic boundary
Node* original_node_pt = mesh_pt()->boundary_node_pt(3,3);
Node* periodic_node_pt = mesh_pt()->boundary_node_pt(1,3);
//Make the periodic_node_pt periodic the data from the original_node_pt
Figure of a mesh that is periodic in the horizontal direction.

Although it is possible to make nodes periodic at any time, it is usually easier to determine which nodes should be "connected" during mesh construction. We therefore strongly recommend to implement periodic boundary conditions inside the mesh constructor. The source code for the constructor of the


that we used in the above problem, illustrates a possible implementation.

Periodic boundaries in spatially adaptive computations

We note that the application of periodic boundary conditions in spatially adaptive computations is slightly more complicated because of the possible presence of hanging nodes on the periodic boundaries. We refer to another tutorial for a discussion of this aspect.


  1. Show that in the present problem the time-periodic solution can also be obtained without applying periodic boundary conditions. Show this mathematically and "by trial and error" (i.e. by changing the boolean flag that is passed to the mesh constructor). Explain why the number of unknowns increases when no periodic boundary conditions are applied.
  2. Confirm that the assignment of "history values" for the nodal positions in set_initial_conditions() is essential.

Source files for this tutorial

PDF file

A pdf version of this document is available.