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

The main program

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

  • We start by including the classes for fluid dynamics: 

       #include "OFELI.h"
   #include "Fluid.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]);
   Verbosity = data.getVerbose();
   int save_flag = data.getSave();
   int pres_flag = data.getIntPar(1);
   theTimeStep = data.getTimeStep();
   theFinalTime = data.getMaxTime();

  • As usual we construct a Mesh instance. 

       Mesh ms(data.getMeshFile());

  • We now declare data vectors. The vectors u and p will stand respectively for the velocity and the pressure. We also define the vector bc for Dirichlet boundary conditions. Then we assign the boundary condition using the class Prescription. 

       Vect<double> u(ms), bc(ms), p(ms,NODE_DOF,1);
   Prescription pr(ms,data.getDataFile());
   pr.get(BOUNDARY_CONDITION,bc);
   u.setName("Velocity");
   p.setName("Pressure");

  • 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,IOField::OUT);
   IOField pf(data.getMeshFile(),data.getString("p_file"),ms,IOField::OUT);

  • Now we declare an instance of class NSP2DQ41 that solves the 2-D incompressible Navier-Stokes equations using the Q1/P0 finite element with a penalty formulation to take into account the incompressibility. Then we define the vector p as the pressure field and transmit boundary conditions to this instance. 

       NSP2DQ41 eq(ms,u);
   eq.setInput(BOUNDARY_CONDITION,bc);
   eq.setInput(PRESSURE_FIELD,p);

  • Once all the data are read, we start the time loop and run each time step. 

       TimeLoop {
      eq.runOneTimeStep();

  • We can save the resulting velocity and pressure vectors for plotting or any other purpose. For the sake of clarity we assign the time value to these vectors. 

          u.setTime(theTime);
      p.setTime(theTime);
      if (save_flag) {
         v_file.put(u);
         p_file.put(p);
      }
   }

      

    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.