Demo: A 2-D incompressible CFD problem using triangles and a projection time integration method

  

The main program

We present an program that solves the 2-D time-dependent incompressible Navier-Stokes equations. For this demo, we solve these equations by a time integration projection method. Space discretisation is handled by a triangular finite element method.

  • We start by including the classes for fluid dynamics: 

       #include "OFELI.h"
   #include "Fluid.h"

  • Like for the previous example, we define an instance of class IPF to handle project data. We also define as data the Reynolds number Re. 

       IPF proj(argv[1]);
   Verbosity = proj.getVerbose();
   Mesh mesh(proj.getMeshFile());
   theTimeStep = proj.getTimeStep();
   theFinalTime = proj.getMaxTime();
   int plot_flag = proj.getPlot();
   double Re = proj.getDouble("Reynolds");

  • To store the solution we define two files, for velocity and pressure by using the class IOField. The resulting velocity and pressure vectors will be stored in corresponding files for postprocessing. 

       IOField vff(proj.getMeshFile(),proj.getString("v_file"),mesh,IOField::OUT);
   IOField pff(proj.getMeshFile(),proj.getString("p_file"),mesh,IOField::OUT);

  • We now define vectors that will contain the solution (u and p for velocity and pressure respectively) and then we define an instance of class TINS2D3TS that solves 2-D incompressible Navier-Stokes equations using the projection method with P1 elements for velocity and pressure. We also trasmit two data to this instance. 

       Vect<double> u(mesh,"Velocity",0), p(mesh,"Pressure",0,1);
   TINS2DT3S eq(mesh);
   eq.Reynolds(Re);
   eq.setTolerance(proj.getTolerance());

    We have next to declare vectors for initial conditions, boundary conditions, ... To initialize these vectors, the class Prescription can be used. We note that these conditions are read in the same file as the project file. For instance to assign initial conditions, we can retrieve this information in this file, then transmit it to the equation by the member function setInput: 

       Vect<double> bc(mesh), bf(mesh);
   Prescription pr(mesh,proj.getDataFile());
   pr.get(INITIAL_FIELD,u);
   eq.setInput(INITIAL_FIELD,u);

  • We are now ready to start computations, i.e. looping over time steps. We have to assign necessary data at each time step and then run the step.

    For each datum, we use the function get to obtain from the Prescription instance the vector, and then the function set to transmit to the TINS2DT3S instance the vector. We may also note that at each time step the vectors u and p contain respectively the velocity and pressure for the current time step, so we can save them in appropriate files. 

       TimeLoop {
      eq.setInput(BOUNDARY_CONDITION,pr.get(BOUNDARY_CONDITION,theTime));
      eq.setInput(BODY_FORCE,pr.get(BODY_FORCE,theTime));
      eq.runOneTimeStep();
      if (plot_flag>0 && theStep%plot_flag==0) {
         p.setTime(theTime); u.setTime(theTime);
         vff.put(u); pff.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.01"/>
      <max_time value="0.5"/>
      <verbose value="1"/>
      <plot value="10"/>
      <mesh_file value="cavity-30.m"/>
      <parameter label="Reynolds" value="1"/>
      <parameter label="v_file" value="cavity.v"/>
      <parameter label="p_file" value="cavity.p"/>
   </Project>
   <Prescription>
      <BoundaryCondition code="2" dof="1">1</BoundaryCondition>
   </Prescription>
</OFELI_File>

Note that we have also a section describing prescription data where a boundary condition is given: We impose the value 1.0 to the horizontal velocity dof="1" for nodes with a code equal to 2. All other data that are not given are by default set to 0.