Fluid Dynamics Demo 1: A 2-D incompressible CFD problem using quadrilateral elements

The main program

Let us present an OFELI program that solves the 2-D time dependent incompressible Navier-Stokes equations.

We start by including the classes for fluid dynamics. We also use here a user defined class for prescibing data:
#include "OFELI.h" #include "Fluid.h" #include "User.h"

We define an instance of class IPF to deal with project data. We extract useful information from this instance. Note that all these parameters were described in the previous examples.

   IPF data("tiff2 - 1.0",argv[1]);
   int output_flag = data.getOutput();
   int save_flag = data.getSave();
   int pres_flag = data.getIntPar(1);
   theTimeStep = data.getTimeStep();
   theFinalTime = data.getMaxTime();

As usual we construct a Mesh instance. In addition, we present here another method for prescribing problem data. This one consists in constructing a class that inherits from abstract class UserData. This class is presented after the main program.
Mesh ms(data.getMeshFile()); User ud(ms);
The linear system uses a symmetric matrix. For this, we we use an instance of class SkSMatrix<double>. Then we declare a vector b that will contain the right-hand side and then the solution. We also need a vector u that contains the solution at the previous time step:
SkSMatrix<double> A(ms); Vect<double> b(ms), u(ms);
We now use the instance ud to prescribe initial solution, Dirichlet boundary condition, body force and source force (tractions): SkSMatrix<double>. Then we declare a vector that will contain the right-hand side and then the solution:
ud.setInitialData(u); Vect<double> bc(ms); ud.setDBC(bc); Vect<double> body_f(ms); ud.setBodyForce(body_f); Vect<double> bound_f(ms); ud.setSurfaceForce(bound_f);
In order to store the solution we define 2 instances of class IOField to store respectively the velocity and the pressure. The names of the storage files were primarily obtained from the project file.
IOField vf(data.getMeshFile(),data.getString("v_file"),ms,OUT); IOField pf(data.getMeshFile(),data.getString("p_file"),ms,OUT);
Once all the data are read, we start the time loop. At first, we set to 0 the right-hand side that will be assembled.
TimeLoop { b = 0;

The global matrix and right-hand side are, as usual, filled by looping over elements. For each element, an instance of class NSP2DQ41 is considered. This class deals with quadrilateral elements. In this formulation, a penalty term takes into account the divergence free constraint. We add for this the function Penal that adds the contribution of penalty term with a penalty parameter equal here to 10⁷. We also add the (lumped) mass matrix and the viscosity term with viscosity equal here to 0.1. We then add the body forces using the User instance. Assembly of element matrix is performed only at the first time step, the matrix being constant. Assembly of the element right-hand side is also made.

      MeshElements(ms) {
         NSP2DQ41 eq(theElement,u,theTime);
         eq.LMass(1./theTimeStep);
         eq.Penal(1.e07);
         eq.Viscous(0.1);
         eq.RHS_Convection();
         eq.BodyRHS(ud);
         if (theStep==1)
            eq.ElementAssembly(A);
         eq.ElementAssembly(b);
      }

Since we want to take account for boundary forces, we have to loop over sides (here edges) and construct element loads. This has also been considered in previous lessons.
MeshSides(ms) { NSP2DQ41 eq(theSide); eq.BoundaryRHS(ud); eq.SideAssembly(b); }
We can now prescribe boundary conditions by penalty method and solve the resulting linear system by factorization and backsubstitution. We finally assign the resulting solution to the vector u.
A.Prescribe(ms,b,bc,theStep-1); if (theStep == 1) A.Factor(); A.Solve(b); u = b;
We can save the resulting velocity vector for plotting or any other purpose. For this, the class IOField is used.
Note that we need to convert to class NodeVect<double> in order to save the vector with some mesh information.
NodeVect<double> uf(ms,u,2,"Velocity",theTime); if (save_flag) vf.put(uf);

If the flag pres_flag is true, then we will compute a posteriori the pressure field. This is easy when the penalty formulation is used: We start by computing a piecewise constant pressure field, then use a smoothing procedure to retrieve a continuous piecewise linear field.

      if (pres_flag) {
         double pres = 0.;
         Vect<double> pm(ms.getNbNodes());
         ElementVect<double> ep(ms,1);
         Reconstruction rr(ms);
         NodeVect<double> p(ms,1,"Pressure",theTime);
         MeshElements(ms) {
            NSP2DQ41 eq(theElement,u,theTime);
            pres += eq.Pressure(1.e07);
         }
         ep(theElementLabel) += pres;
         rr.P0toP1(ep,p);
         if (save_flag)
            pf.put(p);
      }
   }

The class 'User'

We define here the class User that defined the Dirichlet boundary conditions. All other data (body, surface forces and initial conditions) are assumed null, so they don't need be defined. This class inherits from the abstract class UserData. The constructor using mesh is already implemented in UserData. We need here to define the function BoundaryCondition in a very simple way.

class User : public UserData
{
 public:

   User(Mesh &mesh) : UserData(mesh) { }

   double BoundaryCondition(const Point &x, int code, double time=0., size_t dof=1)
   {
      time = 0;
      dof = 1;
      double ret = 0.0;
      if (code==1)
         ret = 0.0;
      if (code==2)
         ret = 1.;
      return ret;
   }
};

An example

Let us run this program with the data presented in the following project file:

<?xml version="1.0" encoding="ISO-8859-1" ?>
<OFELI_File>
<info>
   <title></title>
   <date></date>
   <author></author>
</info>
<Project name="cavity">
   <time_step value="0.1"/>
   <max_time value="1.0"/>
   <verbose value="1"/>
   <output value="1"/>
   <mesh_file value="cavity.m"/>
   <parameter label="v_file" value="cavity.v"/>
   <parameter label="p_file" value="cavity.p"/>
</Project>
</OFELI_File>

Note that only execution and time integration parameters are given. All other data were given through the class User.