Research in Fluid Mechanics

This web page is an evolving progress report on a joint research with Ana Maria Soane on computational fluid dynamics.

Development of computational software

We have developed finite element solver for the time-dependent Navier-Stokes equations \begin{align} &\rho\big(u_t + (\nabla u)u\big) = \mu \Delta u - \nabla p + f \\ &\div u = 0 \end{align} in two dimensional polygonal domains implementing the Argyris $C^1$-compatible finite elements on unstructured triangulated meshes. The software, written entirely in standard C, implements the algorithm described in:

Jian-Guo Liu, Jie Liu, Robert L. Pego, Stability and convergence of efficient Navier-Stokes solvers via a commutator estimate, Communications in Pure and Applied Mathematics, 60 (2007), pp. 1443–1487.

The algorithm calls for finite element approximation of the solutions in $H^2$ rather than the more common $H^1$ spaces. The implementation of $H^2$ finite elements is complicated by the fact that the shape functions must match as $C^1$ functions across element boundaries. On a triangular mesh this may be accomplished with fifth degree polynomials in two variables, that is, various linear combinations of the 21 monomials:

¹
x          y
x²          xy          y²
x³          x²y          xy²          y³
x⁴          x³y          x²y²          xy³          y⁴
x⁵          x⁴y          x³y²          x²y³          xy⁴          y⁵

Shown below are seven of the 21 Argyris shape function on an equilateral triangle. The other 14 are obtained by rotating these by 120 and 240 degrees.

Technical details

Our solver relies on several auxiliary utilities/libraries, including:

Triangulation:

The domain is meshed using Triangle, a triangulation software developed by Jonathan Shewchuk.

Integration:

Integration on triangles is performed using the algorithm described in:

Mark A. Taylor, Beth A. Wingate, Len P. Bos, Several new quadrature formulas for polynomial integration in the triangle, http://arxiv.org/abs/math.NA/0501496.

Sparse system solver:

The resulting linear systems are solved using the sparse system solver UMFPACK developed by Tim Davis.

Visualization:

The solution is written to a file in a format suitable for reading into Maple for post-processing and plotting.

Corner singularities

In domains with reentrant corners, such as in an L-shaped domain, the solution of the Navier-Stokes equations need not be in $H^2$ therefore the algorithm of Liu, Liu and Pego needs to be modified. We have adjusted the algorithm by recasting it in weighted Sobolev spaces where appropriate weights at reentrant corners compensate for the solution's singularities. The numerical results obtained this way show good agreement with those obtained by the usual $H^1$ methods. Currently we are investigating the connection between such corner singularities and weighted Sobolev spaces.

Gallery

In the following sections we give representative samples of solutions computed using our software.

Testing against an exact solution

The velocity field $u=(u_1,u_2)$ and the pressure field $p$ given by: \begin{align} u_1 &= \cos t \cos^2\frac{\pi x}{2} \cos\frac{\pi y}{2} \sin\frac{\pi y}{2}, \\ u_2 &= - \cos t \cos\frac{\pi x}{2} \sin\frac{\pi x}{2} \cos^2\frac{\pi y}{2}, \\ p &= \cos t \cos\frac{\pi x}{2} \sin\frac{\pi y}{2} \end{align} solve the Navier-Stokes equations in the domain $\Omega=(-1,1)\times(-1,1)$ with zero velocity on the boundary conditions. The body force $f$ may be computed by plugging these into the equations. The figures below show the absolute errors in the solution obtained using our finite elements solver.

Top row: Graphs of absolute errors in the x and y components of velocities. Bottom row, left: The graph of the absolute error in the pressure. Bottom row, right: The graph of the divergence of the velocity which should be zero in an ideal solver. Note that the vertical scale is quite small in all cases.

The driven cavity

We solve the Navier-Stokes equations for steady-state flow in a square two-dimensional cavity where the top edge (lid) moves at a constant horizontal velocity. No exact solution is available therefore we compare our solution against that obtained using the commercial Comsol Multiphysics (Femlab) software.

Left column: Graphs of the x and y components of velocities, the pressure and the divergence of velocity, all produced using our H2 solver. Right column: The corresponding graphs produced by Femlab's H1 solver.

Flow over a backstep

The diagram below depicts a flow region where the fluid enters from the left edge, flows over a step, and leaves from the right edge.

Flow over a backstep: Fluid enters from the left and leaves from the right. Streamlines computed and plotted in Femlab.

No exact solution is available for flow over a backstep, therefore we compare our solution against that produced by Femlab.

Left column: Graphs of the $x$ and $y$ components of velocities, the pressure and the divergence of velocity, all produced using our $H^2$ solver. Right column: The corresponding graphs produced by Femlab's H1 solver.

Publications

1. Ana Maria Soane and Rouben Rostamian, Free boundary problems in fluid mechanics, in Susan Friedlander, Barbara Keyfitz, Irene Gamba, and Krystyna Kuperberg, editors, Women in Mathematics: The Legacy of Ladyzhenskaya and Oleinik. Proceedings of an MSRI conference, May 2006, Berkeley.
2. Ana Maria Soane and Rouben Rostamian, Variational problems in weighted Sobolev spaces on non-smooth domains (draft of 2008), to appear in Quarterly of Applied Mathematics.
3. Ana Maria Soane and Rouben Rostamian, The optimal convergence rate of a C¹ finite element method for non-smooth domains (draft of 2009), Journal of Computational and Applied Mathematics (under review).

This web page was created on 2007–03–29 and was revised last on 2009–08–21.