Example problem: Finite-Reynolds-number flow inside an oscillating ellipse

In this example we consider our first moving-boundary Navier-Stokes problem: The flow of a viscous fluid contained in an elliptical ring whose walls perform periodic oscillations.

oomph-lib's Navier-Stokes elements are based on the Arbitrary Lagrangian Eulerian (ALE) form of the Navier-Stokes equations and can therefore be used in moving domain problems. In this example we illustrate their use in Domain - based meshes (first discussed in the example demonstrating the solution of the unsteady heat equation in a moving domain) in which MacroElements are used to update the nodal positions in response to changes in the domain boundary. In subsequent examples, we will discuss alternative, sparse mesh update techniques that are useful in problems with free boundaries and in fluid-structure interaction problems.

Finite-Reynolds-number-flow driven by an oscillating ellipse

We consider the unsteady 2D flow of a Newtonian fluid that is contained in an oscillating elliptical ring whose wall shape is parametrised by the Lagrangian coordinate $ \xi $ as

\[ {\bf R}_w(\xi) = \left( \begin{array}{l} a(t) \cos(\xi) \\ a^{-1}(t) \sin(\xi) \end{array} \right) \]


\[ a(t) = A + \widehat{A} \cos\left(\frac{2\pi t}{T}\right). \]

$ A $ represents the average half-axis of the elliptical ring in the $ x_1 $-direction, and $ \widehat{A} $ is the amplitude of its periodic variation. The ring has constant cross-sectional area – consistent with the incompressibility of the fluid whose motion is governed by the ALE form of the Navier-Stokes equations,

\[ Re\left(St\frac{\partial u_i}{\partial t} + \left(u_j-u_j^M\right)\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, \]

where $ u_j^M $ is the mesh velocity. We exploit the symmetry of the problem and solve the equations in the quarter domain

\[ D = \left\{(x_1,x_2) \, \bigg| \, x_1\geq0, \, x_2\geq0,\, \left(\frac{x_1}{a(t)}\right)^2 + \big(x_2 \, a(t)\big)^2 \leq 1\right\}, \]

shown in this sketch (for $ A=1, \ \widehat{A}=1 $ and $T=1$),

Sketch of the computational domain.

The fluid is subject to no-slip boundary conditions on the curved wall,

\[ \left. {\bf u} \right|_{\partial D_{ellipse}} = \frac{\partial {\bf R}_w(\xi, t)}{\partial t} \]

and symmetry conditions on the symmetry boundaries,

\[ \left.u_1\right|_{x_2=0}=0, \qquad \left.u_2\right|_{x_1=0}=0. \]

The initial conditions for the velocity 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 a given, divergence-free velocity field that satisfies the velocity boundary conditions at $ t=0 $. No initial conditions are required for the pressure.

An exact solution

It is easy to show (by inspection) that the unsteady stagnation point flow

\[ u_1(x_1,x_2,t) = \frac{1}{a}\frac{da}{dt}\, x_1 = -\frac{2\pi\widehat{A}\sin\left(\frac{2\pi t}{T}\right)} {T\left(A+\widehat{A}\cos\left(\frac{2\pi t}{T}\right)\right)} \, x_1 \mbox{\ \ \ \ and \ \ \ \ \ } u_2(x_1,x_2,t) = -\frac{1}{a}\frac{da}{dt} \, x_2 = -\frac{2\pi\widehat{A}\sin\left(\frac{2\pi t}{T}\right)} {T\left(A+\widehat{A}\cos\left(\frac{2\pi t}{T}\right)\right)} \, x_2, \]

is an exact solution of the above problem as it satisfies the Navier-Stokes equations and the velocity boundary conditions. The pressure is given by

\[ p = \frac{2\pi\widehat{A}\, Re\, \left(x_1^2\, St\, \cos\left(\frac{2\pi t}{T}\right)A + x_1^2\, St\, \widehat{A} - x_1^2\widehat{A} + x_1^2\widehat{A}\cos^2\left(\frac{2\pi t}{T}\right) - x_2^2\, St\, \cos\left(\frac{2\pi t}{T}\right) - x_2^2\, St\, \widehat{A} - x_2^2\widehat{A} + x_2^2\widehat{A}\cos^2\left(\frac{2\pi t}{T}\right)\right)} {T^2\left(A^2 + 2A\widehat{A}\cos\left(\frac{2\pi t}{T}\right) + \widehat{A}^2\cos^2\left(\frac{2\pi t}{T}\right)\right)}. \]


The two figures below show two snapshots of the solution for $ Re=Re \ St=100 $ , extracted from an animations of the results computed with Taylor-Hood and Crouzeix-Raviart elements. In both cases, the exact solution was used as the initial condition for the velocities. The figures show "carpet plots" of the two velocity components and the pressure, respectively, and a contour plot of the pressure, superimposed on the moving mesh. The carpet plot of the velocities clearly shows that the flow is of stagnation-point type as the horizontal velocity, $u_1$ , is a linear function of $x_1$ while the vertical velocity, $u_2$ , is a linear function of $ - x_2$.

Plot of the velocity and pressure fields computed with 2D Crouzeix-Raviart elements, with Re=100 and St=1.
Plot of the velocity and pressure fields computed with 2D Taylor-Hood elements, with Re=100 and St=1.

The moving wall

As usual, we represent the moving wall as a GeomObject and define its shape by implementing the pure virtual function GeomObject::position(...). The arguments to the constructor specify the mean half-axis of the ellipse, $ A $, the amplitude of its variations, $ \widehat{A} $, and the period of the oscillation, $ T $ . We also pass the pointer to a Time object to the constructor and store it in a private data member, to allow the position(...) functions to access the current value of the continuous time.

/// \short Oscillating ellipse
/// \f[ x = (A + \widehat{A} \cos(2\pi t/T)) \cos(\xi) \f]
/// \f[ y = \frac{\sin(\xi)}{A + \widehat{A} \cos(2\pi t/T)} \f]
/// Note that cross-sectional area is conserved.
class MyEllipse : public GeomObject
/// \short Constructor: Pass initial x-half axis, amplitude of x-variation,
/// period of oscillation and pointer to time object.
MyEllipse(const double& a, const double& a_hat,
const double& period, Time* time_pt) :
GeomObject(1,2), A(a), A_hat(a_hat), T(period), Time_pt(time_pt) {}
/// Destructor: Empty
virtual ~MyEllipse() {}
/// \short Current position vector to material point at
/// Lagrangian coordinate xi
void position(const Vector<double>& xi, Vector<double>& r) const
// Get current time:
double time=Time_pt->time();
// Position vector
double axis=A+A_hat*cos(2.0*MathematicalConstants::Pi*time/T);
r[0] = axis*cos(xi[0]);
r[1] = (1.0/axis)*sin(xi[0]);
/// \short Parametrised position on object: r(xi). Evaluated at
/// previous time level. t=0: current time; t>0: previous
/// time level.
void position(const unsigned& t, const Vector<double>& xi,
Vector<double>& r) const
// Get current time:
double time=Time_pt->time(t);
// Position vector
double axis=A+A_hat*cos(2.0*MathematicalConstants::Pi*time/T);
r[0] = axis*cos(xi[0]);
r[1] = (1.0/axis)*sin(xi[0]);
/// x-half axis
double A;
/// Amplitude of variation in x-half axis
double A_hat;
/// Period of oscillation
double T;
/// Pointer to time object
Time* Time_pt;
}; // end of MyEllipse

The global parameters

As in most previous examples, we use a namespace to define and initialise global problem parameters such as the Reynolds and Strouhal numbers:

/// Namepspace for global parameters
namespace Global_Physical_Variables
/// Reynolds number
double Re=100.0;
/// Womersley = Reynolds times Strouhal
double ReSt=100.0;

We also define and initialise the parameters that specify the motion of the domain boundary and specify the exact solution.

/// x-Half axis length
double A=1.0;
/// x-Half axis amplitude
double A_hat=0.1;
/// Period of oscillations
double T=1.0;
/// Exact solution of the problem as a vector containing u,v,p
void get_exact_u(const double& t, const Vector<double>& x, Vector<double>& u)
using namespace MathematicalConstants;
// Strouhal number
double St = ReSt/Re;
// Half axis
double a=A+A_hat*cos(2.0*Pi*t/T);
double adot=-2.0*A_hat*Pi*sin(2.0*Pi*t/T)/T;
// Velocity solution
// Pressure solution
u[2]=(2.0*A_hat*Pi*Pi*Re*(x[0]*x[0]*St*cos(2.0*Pi*t/T)*A +
x[0]*x[0]*St*A_hat - x[0]*x[0]*A_hat +
x[0]*x[0]*A_hat*cos(2.0*Pi*t/T)*cos(2.0*Pi*t/T) -
x[1]*x[1]*St*cos(2.0*Pi*t/T)*A -
x[1]*x[1]*St*A_hat - x[1]*x[1]*A_hat +
x[1]*x[1]*A_hat*cos(2.0*Pi*t/T)*cos(2.0*Pi*t/T) ))
/(T*T*(A*A + 2.0*A*A_hat*cos(2.0*Pi*t/T) +
A_hat*A_hat*cos(2.0*Pi*t/T)*cos(2.0*Pi*t/T) ));
} // end of namespace

The driver code

As in most previous unsteady demo codes, we allow the code to be run in a validation mode (in which we use a coarser mesh and execute fewer timesteps). This mode is selected by specifying an (arbitrary) command line argument that we store in the namespace CommandLineArgs.

/// Driver code for unsteady Navier-Stokes flow, driven by
/// oscillating ellipse. If the code is executed with command line
/// arguments, a validation run is performed.
int main(int argc, char* argv[])
// Store command line arguments

We create a DocInfo object to specify the output directory, build the problem with adaptive Crouzeix-Raviart elements and the BDF<2> timestepper and perform the unsteady simulation.

// Solve with Crouzeix-Raviart elements
// Create DocInfo object with suitable directory name for output
DocInfo doc_info;
//Set up problem
// Run the unsteady simulation

Then we repeat this process for adaptive Taylor-Hood elements.

// Solve with Taylor-Hood elements
// Create DocInfo object with suitable directory name for output
DocInfo doc_info;
//Set up problem
// Run the unsteady simulation
}; // end of main

The problem class

Most of the problem class is a straightforward combination of the problem classes employed in the simulation of the adaptive driven cavity and Rayleigh channel problems, in that the problem combines unsteadiness with spatial adaptivity (though in the present problem the adaptivity is only used to uniformly refine the very coarse base mesh; we refer to ***another example*** for the use of full spatial adaptivity in a moving-domain Navier-Stokes problem).

/// Navier-Stokes problem in an oscillating ellipse domain.
template<class ELEMENT, class TIMESTEPPER>
class OscEllipseProblem : public Problem
/// Constructor
/// Destructor (empty)
/// Update the problem specs after solve (empty)
/// \short Update problem specs before solve (empty)
/// Actions before adapt (empty)
/// Actions after adaptation, pin relevant pressures
// Unpin all pressure dofs
// Pin redundant pressure dofs
// Now set the first pressure dof in the first element to 0.0
} // end of actions_after_adapt

The key new feature in the current problem is the presence of the moving domain which requires updates of

  1. all nodal positions
  2. the prescribed velocities on the moving wall via the no-slip condition.

before every timestep. Since the nodal positions of the QuarterCircleSectorMesh are determined via its MacroElement / Domain representation (which updates the nodal position in response to changes in the geometry of the GeomObjects that define its boundaries), the former task may be accomplished by executing the Mesh::node_update() function; the update of the no-slip condition may be performed by calling the function FSI_functions::apply_no_slip_on_moving_wall(Node* node_pt), a helper function, defined in the namespace FSI_functions, which updates the velocity components $ u_1, u_2 [, u_3] $ according to the no-slip boundary condition

\[ {\bf u}_{Node} = \frac{\partial {\bf x}_{Node}}{\partial t} \]

where the time-derivative of the nodal positions is evaluated by the Node's positional timestepper. [Note: The function FSI_functions::apply_no_slip_on_moving_wall(...) assumes that the velocity components are stored in the Node's first 2 [3] values. This is consistent with the storage of the velocity component in all existing Navier-Stokes elements. If you develop your own Navier-Stokes elements and use a different storage scheme you use this function at your own risk.]

Here is the implementation of these tasks:

/// \short Update the problem specs before next timestep
void actions_before_implicit_timestep()
// Update the domain shape
// Ring boundary: No slip; this implies that the velocity needs
// to be updated in response to wall motion
unsigned ibound=1;
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
/// Update the problem specs after timestep (empty)
void actions_after_implicit_timestep(){}

The remaining functions are similar to those used in our previous Navier-Stokes examples and require no further explanation.

/// Doc the solution
void doc_solution(DocInfo& doc_info);
/// Timestepping loop
void unsteady_run(DocInfo& doc_info);
/// \short Set initial condition
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
} // end_of_fix_pressure
/// Pointer to GeomObject that specifies the domain bondary
GeomObject* Wall_pt;
}; // end of problem_class

The problem constructor

We start by creating a timestepper of the type specified by the Problem's template parameter and add (a pointer to) it to the Problem's collection of Timesteppers. Recall that this function also creates the Problem's Time object.

/// Constructor for Navier-Stokes problem on an oscillating ellipse domain.
template<class ELEMENT, class TIMESTEPPER>
//Create the timestepper and add it to the problem
add_time_stepper_pt(new TIMESTEPPER);

Next we create the GeomObject that defines the curvilinear domain boundary and pass it to the Mesh constructor. (Since we will only use adaptivity to refine the mesh uniformly, it is not necessary to define an error estimator.)

// Setup mesh
// Build geometric object that forms the curvilinear domain boundary:
// an oscillating ellipse
// Half axes
// Variations of half axes
// Period of the oscillation
// Create GeomObject that specifies the domain bondary
Wall_pt=new MyEllipse(a,a_hat,period,Problem::time_pt());
// Start and end coordinates of curvilinear domain boundary on ellipse
double xi_lo=0.0;
double xi_hi=MathematicalConstants::Pi/2.0;
// Now create the mesh. Separating line between the two
// elements next to the curvilinear boundary is located half-way
// along the boundary.
double fract_mid=0.5;
Problem::mesh_pt() = new RefineableQuarterCircleSectorMesh<ELEMENT>(
//Z2ErrorEstimator* error_estimator_pt=new Z2ErrorEstimator;

Both velocity components on the curvilinear mesh boundary are determined by the no-slip condition and must therefore be pinned,

// Fluid boundary conditions
// Ring boundary: No slip; this also implies that the velocity needs
// to be updated in response to wall motion
unsigned ibound=1;
unsigned num_nod= mesh_pt()->nboundary_node(ibound);
for (unsigned inod=0;inod<num_nod;inod++)
// Pin both velocities
for (unsigned i=0;i<2;i++)
} // end boundary 1

whereas on the symmetry boundaries only one of the two velocity components is set to zero:

// Bottom boundary:
unsigned num_nod= mesh_pt()->nboundary_node(ibound);
for (unsigned inod=0;inod<num_nod;inod++)
// Pin vertical velocity
} // end boundary 0
// Left boundary:
unsigned num_nod= mesh_pt()->nboundary_node(ibound);
for (unsigned inod=0;inod<num_nod;inod++)
// Pin horizontal velocity
} // end boundary 2

Finally, we pass the pointers to $ Re $, $ Re \ St $ and the global Time object (automatically created by the Problem when the timestepper was passed to it at the beginning of the constructor) to the elements, pin the redundant nodal pressure degrees of freedom (see the discussion of the adaptive driven-cavity problem for more details), pin one pressure degree of freedom, and set up the equation numbering scheme.

// Complete the build of all elements so they are fully functional
// Find number of elements in mesh
unsigned n_element = mesh_pt()->nelement();
// Loop over the elements to set up element-specific
// things that cannot be handled by constructor
for(unsigned i=0;i<n_element;i++)
// Upcast from FiniteElement to the present element
ELEMENT *el_pt = dynamic_cast<ELEMENT*>(mesh_pt()->element_pt(i));
//Set the Reynolds number, etc
el_pt->re_pt() = &Global_Physical_Variables::Re;
el_pt->re_st_pt() = &Global_Physical_Variables::ReSt;
// Pin redundant pressure dofs
// Now set the first pressure dof in the first element to 0.0
// Do equation numbering
cout << "Number of equations: " << assign_eqn_numbers() << std::endl;
} // end of constructor

Assigning the initial conditions

This function assigns "history values" for the velocities and the nodal positions from the exact solution. It is implemented in exactly the same way as in the solution of the unsteady heat equation in a moving domain. Note that because the domain is moving, the nodal positions must be updated (according to the position of the domain boundary at the relevant previous timestep), before evaluating the exact solution at the nodal position.

/// \short Set initial condition: Assign previous and current values
/// from exact solution.
template<class ELEMENT,class TIMESTEPPER>
// Backup time in global timestepper
double backed_up_time=time_pt()->time();
// Past history for velocities must be established for t=time0-deltat, ...
// Then provide current values (at t=time0) which will also form
// the initial guess for first solve at t=time0+deltat
// Vector of exact solution value
Vector<double> soln(3);
Vector<double> x(2);
//Find number of nodes in mesh
unsigned num_nod = mesh_pt()->nnode();
// Get continuous times at previous timesteps
int nprev_steps=time_stepper_pt()->nprev_values();
Vector<double> prev_time(nprev_steps+1);
for (int itime=nprev_steps;itime>=0;itime--)
// Loop over current & previous timesteps (in outer loop because
// the mesh also moves!)
for (int itime=nprev_steps;itime>=0;itime--)
double time=prev_time[itime];
// Set global time (because this is how the geometric object refers
// to continous time )
cout << "setting IC at time =" << time << std::endl;
// Update the mesh for this value of the continuous time
// (The wall object reads the continous time from the
// global time object)
// Loop over the nodes to set initial guess everywhere
for (unsigned jnod=0;jnod<num_nod;jnod++)
// Get nodal coordinates
// Get exact solution (unsteady stagnation point flow)
// Assign solution
// Loop over coordinate directions
for (unsigned i=0;i<2;i++)
} // end of loop over previous timesteps
// Reset backed up time for global timestepper
} // end of set_initial_condition

Post processing

The function doc_solution(...) is similar to that in the unsteady heat examples and the previous Navier-Stokes examples. We add dummy zones and tecplot geometries to facilitate the post-processing of the results with tecplot.

/// Doc the solution
template<class ELEMENT, class TIMESTEPPER>
ofstream some_file;
char filename[100];
// Number of plot points
unsigned npts;
// Output solution
some_file << "TEXT X=2.5,Y=93.6,F=HELV,HU=POINT,C=BLUE,H=26,T=\"time = "
<< time_pt()->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;
// Write dummy zones that force tecplot to keep the axis limits constant
// while the domain is moving.
some_file << "ZONE I=2,J=2" << std::endl;
some_file << "0.0 0.0 -0.65 -0.65 -200.0" << std::endl;
some_file << "1.15 0.0 -0.65 -0.65 -200.0" << std::endl;
some_file << "0.0 1.15 -0.65 -0.65 -200.0" << std::endl;
some_file << "1.15 1.15 -0.65 -0.65 -200.0" << std::endl;
some_file << "ZONE I=2,J=2" << std::endl;
some_file << "0.0 0.0 0.65 0.65 300.0" << std::endl;
some_file << "1.15 0.0 0.65 0.65 300.0" << std::endl;
some_file << "0.0 1.15 0.65 0.65 300.0" << std::endl;
some_file << "1.15 1.15 0.65 0.65 300.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;
// Plot wall posn
unsigned nplot=100;
for (unsigned iplot=0;iplot<nplot;iplot++)
Vector<double> xi_wall(1), r_wall(2);
some_file << r_wall[0] << " " << r_wall[1] << std::endl;
// Increment number of doc
} // end of doc_solution

The timestepping loop

The timestepping loop is extremely straightforward: We choose a timestep and the overall length of the simulation, initialise the timestepper(s) by calling Problem::initialise_dt(...) and assign the initial condition.

/// Unsteady run
template<class ELEMENT, class TIMESTEPPER>
// Specify duration of the simulation
double t_max=3.0;
// Initial timestep
double dt=0.025;
// Initialise timestep
// Set initial conditions.

Next we set the number of timesteps for a normal run.

// Alternative initial conditions: impulsive start; see exercise.
// find number of steps
unsigned nstep = unsigned(t_max/dt);

We over-write this number and perform a single uniform mesh refinement if the code is run in self-test mode (indicated by a non-zero number of command line arguments),

// If validation: Reduce number of timesteps performed and
// use coarse-ish mesh
if (CommandLineArgs::Argc>1)
cout << "validation run" << std::endl;

otherwise we refine the mesh three times and output the initial conditions

// Refine the mesh three times, to resolve the pressure distribution
// (the velocities could be represented accurately on a much coarser mesh).
// Output solution initial

Finally we execute the proper timestepping loop and document the solution after every timestep

// Timestepping loop
for (unsigned istep=0;istep<nstep;istep++)
cout << "TIMESTEP " << istep << std::endl;
cout << "Time is now " << time_pt()->time() << std::endl;
// Take timestep
//Output solution
} // end of unsteady_run

Comments and Exercises

  1. Compare the results of the numerical simulation in which $ {\bf u}_{IC} $ is given by the exact solution (an unsteady stagnation point flow) to that obtained from an "impulsive start" where $ {\bf u}_{IC} = {\bf 0}. $ (This is most easily implemented by replacing the call to set_initial_condition() with a call to Problem::assign_initial_values_impulsive().

    Why do we obtain the same velocity with both initial conditions and why does the pressure take a few timesteps (How many exactly? Compare simulations with BDF<4> and BDF<2> timesteppers.) to "catch up" with the exact solution? [Hint: The unsteady stagnation point flow is a potential flow, therefore the viscous terms in the Navier-Stokes equations disappear. See also chapter 3.19 in Volume 2 of Gresho & Sani's wonderful book "Incompressible Flow and the Finite Element Method".]

Source files for this tutorial

PDF file

A pdf version of this document is available.