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 as
where
represents the average half-axis of the elliptical ring in the -direction, and 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,
and the continuity equation
where is the mesh velocity. We exploit the symmetry of the problem and solve the equations in the quarter domain
shown in this sketch (for and ),
Sketch of the computational domain.
The fluid is subject to no-slip boundary conditions on the curved wall,
and symmetry conditions on the symmetry boundaries,
The initial conditions for the velocity are given by
where is a given, divergence-free velocity field that satisfies the velocity boundary conditions at . No initial conditions are required for the pressure.
An exact solution
It is easy to show (by inspection) that the unsteady stagnation point flow
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
Results
The two figures below show two snapshots of the solution for , 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, , is a linear function of while the vertical velocity, , is a linear function of .
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, , the amplitude of its variations, , and the period of the oscillation, . 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.
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.
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
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).
The key new feature in the current problem is the presence of the moving domain which requires updates of
all nodal positions
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 according to the no-slip boundary condition
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:
/// Update the problem specs before next timestep
void actions_before_implicit_timestep()
{
// Update the domain shape
mesh_pt()->node_update();
// Ring boundary: No slip; this implies that the velocity needs
/// 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'sTime object.
//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:
Finally, we pass the pointers to , 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
// Now set the first pressure dof in the first element to 0.0
fix_pressure(0,0,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.
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.
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.
Next we set the number of timesteps for a normal run.
// Alternative initial conditions: impulsive start; see exercise.
//assign_initial_values_impulsive();
// 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)
{
nstep=2;
refine_uniformly();
cout << "validation run" << std::endl;
}
otherwise we refine the mesh three times and output the initial conditions
else
{
// Refine the mesh three times, to resolve the pressure distribution
// (the velocities could be represented accurately on a much coarser mesh).
refine_uniformly();
refine_uniformly();
refine_uniformly();
}
// Output solution initial
doc_solution(doc_info);
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
unsteady_newton_solve(dt);
//Output solution
doc_solution(doc_info);
}
} // end of unsteady_run
Comments and Exercises
Compare the results of the numerical simulation in which is given by the exact solution (an unsteady stagnation point flow) to that obtained from an "impulsive start" where (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
The source files for this tutorial are located in the directory:
![\[
{\bf R}_w(\xi) = \left(
\begin{array}{l}
a(t) \cos(\xi) \\
a^{-1}(t) \sin(\xi)
\end{array}
\right)
\]](form_1_dark.png)
![\[
a(t) = A + \widehat{A} \cos\left(\frac{2\pi t}{T}\right).
\]](form_2_dark.png)
![\[
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)
\]](form_6_dark.png)
![\[
\frac{\partial u_i}{\partial x_i} = 0,
\]](form_7_dark.png)
![\[
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\},
\]](form_9_dark.png)
![\[
\left. {\bf u} \right|_{\partial D_{ellipse}} = \frac{\partial {\bf
R}_w(\xi, t)}{\partial t}
\]](form_12_dark.png)
![\[
\left.u_1\right|_{x_2=0}=0, \qquad \left.u_2\right|_{x_1=0}=0.
\]](form_13_dark.png)
![\[
{\bf u}(x_1,x_2,t=0) = {\bf u}_{IC}(x_1,x_2),
\]](form_14_dark.png)
![\[
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,
\]](form_17_dark.png)
![\[
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)}.
\]](form_18_dark.png)
![$ u_1, u_2 [, u_3] $](form_25_dark.png)
![\[
{\bf u}_{Node} = \frac{\partial {\bf x}_{Node}}{\partial t}
\]](form_26_dark.png)
