In this example we shall demonstrate the spatially adaptive solution of the steady 3D NavierStokes equations using the problem of developing pipe flow.
and in the quartertube domain , subject to the Dirichlet boundary conditions: on the curved wall (where is the outer unit normal vector) on the symmetry boundaries and on the inflow boundary, and finally (parallel flow) on the outflow boundary Note that the axial velocity component, , is not constrained at the outflow. Implicitly, we are therefore setting the axial component of traction on the fluid to zero, Since in the outflow crosssection [see (5], and the flow is incompressible, this is equivalent to (weakly) setting the pressure at the outflow to zero,

The figure below shows the results computed with oomphlib's
3x3x3node 3D adaptive TaylorHood elements for the parameters , and . The large exponent imposes a very blunt inflow profile, which creates a thin boundary layer near the wall. Diffusion of vorticity into the centre of the tube smooths the velocity profile which ultimately approaches a parabolic Poiseuille profile. If you are viewing these results online you will be able to see how successive mesh adaptations refine the mesh – predominantly near the entry region where the large velocity gradients in the boundary layer require a fine spatial discretisation. Note also that on the coarsest mesh, even the (imposed) inflow profile is represented very poorly.
The problem only contains one global parameter, the Reynolds number, which we define in a namespace, as usual.
Since the 3D computations can take a long time, and since all demo codes are executed during oomphlib's
selftest procedures, we allow the code to operate in two modes:
The main code therefore starts by storing the command line arguments and setting the adaptation targets accordingly:
We then create a DocInfo
object to specify the labels for the output files, and solve the problem, first with TaylorHood and then with CrouzeixRaviart elements, writing the results from the two discretisations to different directories:
The problem class is very similar to the ones used in the 2D examples. We pass the DocInfo
object and the target errors to the Problem
constructor.
The function Problem::actions_after_newton_solve()
is used to document the solutions computed at various levels of mesh refinement:
The function Problem::actions_before_newton_solve()
is discussed below, and, as in all adaptive NavierStokes computations, we use the function Problem::actions_after_adapt()
to pin any redundant pressure degrees of freedom; see another tutorial for details.
Finally, we have the usual doc_solution()
function and include an access function to the mesh. The private member data Alpha
determines the bluntness of the inflow profile.
We start by building the adaptive mesh for the quarter tube domain. As for most meshes with curvilinear boundaries, the RefineablequarterTubeMesh
expects the curved boundary to be represented by a GeomObject
. We therefore create an EllipticalTube
with unit half axes, i.e. a unit cylinder and a pass a pointer to the GeomObject
to the mesh constructor. The "ends" of the curvilinear boundary (in terms of the maximum and minimum values of the Lagrangian coordinates that parametrise the shape of the GeomObject
) are such that it represents a quarter of a cylindrical tube of length .
Next, we build an error estimator and specify the target errors for the mesh adaptation:
Now we have to apply boundary conditions on the various mesh boundaries. [Reminder: If the numbering of the mesh boundaries is not apparent from its documentation (as it should be!), you can use the function Mesh::output_boundaries
(...) to output them in a tecplotreadable form.
If the mesh has N boundaries, the output file will contain N different zones, each containing the coordinates of the nodes on the boundary.]
For the RefineableQuarterTubeMesh
, the boundaries are numbered as follows:
GeomObject
that specifies the wall.GeomObject
that specifies the wall.We apply the following boundary conditions:
Now we assign the re_pt()
for each element and pin the redundant nodal pressures (see another tutorial for details).
We provide an initial guess for the velocity field by initialising all velocity components with their Poiseuille flow values.
Finally, we set the value of Alpha
, the exponent that specifies the bluntness of the inflow profile and assign the equation numbers.
We use the function Problem::actions_before_newton_solve()
to reassign the inflow boundary conditions before every solve. In the present problem this is an essential step because the blunt inflow profile (4) cannot be represented accurately on the initial coarse mesh (see the animation of the axial velocity profiles at the beginning of this document). As discussed in the example that illustrates the use of spatial adaptivity for timedependent problems, oomphlib
automatically (i) applies the correct boundary conditions for newly created nodes that are located on the mesh boundaries and (ii) assigns the nodal values at such nodes by interpolation from the previously computed solution. This procedure is adequate on boundaries where homogeneous boundary conditions are applied, e.g. on the curved wall, the symmetry and the outflow boundaries. However, on the inflow boundary, the interpolation from the FE representation of the blunt velocity profile (imposed on the coarse initial mesh) onto the refined mesh, does not yield a more accurate representation of the prescribed inflow profile. It is therefore necessary to reassign the nodal values on this boundary after every adaptation, i.e. before every solve.
This function remains exactly the same as in the 2D examples.
Problem::actions_before_newton_solve()
to confirm that this step is essential if the computation is to converge to the exact solution.Problem::Max_residuals
. Try increasing the value of this threshold in the Problem constructor. Is this sufficient to make the Newton method converge?]A pdf version of this document is available.