        load "ppm2rnm"


// ========================================================
// Création du maillage
// ========================================================

// Paramètres
int nbtj=20; // Nb de triangles autour de chaque objet
int nbtbord=10; // Nb de triangles aux bords du domaine


int nx = 31;
int ny = 31;
    
// Maillage pour isoline
mesh Th = square(nx-1, ny-1, [(nx-1)*(x), (ny-1)*(1-y)]);

// Fespace
fespace Vh(Th, P1);





border C01(t=0, 1){x = t*(nx-1);y = 0;label=3;}  // Sud
border C02(t=0, 1){x = nx-1;y = t*(ny-1);label=4;}  // Est
border C03(t=0, 1){x = (1-t)*(nx-1);y = (ny-1);label=5;}  // Nord
border C04(t=0, 1){x = 0;y = (1-t)*(ny-1);label=6;}  // Ouest

      cout << "START MESH EXPORT" << endl;

      int NbVertices = Th.nv;
      int NbTriangles = Th.nt;
      int NbBoundaryElements = Th.nbe;
      
      cout << NbVertices << " " << NbTriangles << " " << NbBoundaryElements << endl;
      for (int i = 0; i < NbVertices; i++){
        cout << Th(i).x << " " << Th(i).y << " " << Th(i).label << endl;
      }
      for (int k = 0; k < NbTriangles; k++){
        cout << int(Th[k][0])+1 << " " << int(Th[k][1])+1 << " " << int(Th[k][2])+1 << " " << Th[k].label << endl;
      }
      for (int k = 0; k < NbBoundaryElements; k++){
        cout << Th.be(k)[0]+1 << " " << Th.be(k)[1]+1 << " " << Th.be(k).label << endl;
      }

      cout << "END MESH EXPORT" << endl;

      plot(Th, ps="figure.eps");



      // ========================================================
      // Résolution du problème
      // ========================================================

      // Paramètres
      // real nu = 1.56e-5;
      real nu = 1.;
      real dt = 0.1;
      real T = 2;  // Temps final
      int i;
      // Fespace
      fespace Uh(Th,P2);  // Espace pour la vélocité
      Uh u1, u2, v1, v2;
      Uh u1old, u2old;

      fespace Ph(Th,P1);  // Espace pour la pression
      Ph p, q;

      // Problème de Navier-Stokes
      // Par linéarisation
      problem NavierStokes1(u1, u2, p, v1, v2, q, solver=Crout, init = i)
          =   int2d(Th)
              (
                      (u1*v1 + u2*v2)/dt
                  +   nu*(dx(u1)*dx(v1) + dy(u1)*dy(v1)
                      +   dx(u2)*dx(v2) + dy(u2)*dy(v2))
                  - p*dx(v1) - p*dy(v2)
                  - q*dx(u1) - q*dy(u2)
                  - p*q*(0.000001)  // Pour la stabilité et permet d'utiliser le solveur Crout
              )
            + int2d(Th)
              (
                  - u1old*v1/dt - u2old*v2/dt
              )
            + int2d(Th)
              (
                  u1old*(dx(u1)*v1 + dy(u1)*v2)
                  + u2old*(dx(u2)*v1 + dy(u2)*v2)
              )
            + on(1, 2, 7, 8, 9, 10, u1=0, u2=0)  // Conditions de Dirichlet sur les avions
            + on(wall1, wall2, wall3, u1=Vx, u2=Vy)  // Conditions de Dirichlet à l'Est, au Nord et à l'Ouest
          ;

      // Méthode des caractéristiques
      problem NavierStokes2(u1, u2, p, v1, v2, q, solver = Crout, init = i)
          =   int2d(Th)
              (
                      (u1*v1 + u2*v2)/dt
                  +   nu*(dx(u1)*dx(v1) + dy(u1)*dy(v1)
                      +   dx(u2)*dx(v2) + dy(u2)*dy(v2))
                  - p*dx(v1) - p*dy(v2)
                  - q*dx(u1) - q*dy(u2)
                  - p*q*(0.000001)  // Pour la stabilité et permet d'utiliser le solveur Crout
              )
            + int2d(Th)
              (
                  -   convect([u1old, u2old], -dt, u1old)*v1/dt
                  -   convect([u1old, u2old], -dt, u2old)*v2/dt
              )
            + on(1, 2, 7, 8, 9, 10, u1=0, u2=0)  // Conditions de Dirichlet sur les avions
            + on(wall1, wall2, wall3, u1=Vx, u2=Vy)  // Conditions de Dirichlet à l'Est, au Nord et à l'Ouest
          ;


      // Boucle pour le temps
      u1 = 0;
      u2 = 0;
      p=0;

      i = 0;
      string[int] colour(21);
      colour[0] = "figure1.eps";
      colour[4] = "figure2.eps";
      colour[8] = "figure3.eps";
      colour[12] = "figure4.eps";
      colour[16] = "figure5.eps";
      colour[20] = "figure6.eps";
      int compteur = 0;

      //for (real t = 0; t <= T; t += dt)
      real t = 0;
      for (int compteur = 0; compteur < 24; compteur += 1)
      {
          if (i%10 == 0) i=0;  // On recalcule la matrice tous les 10 pas de temps
          if (choix == 1)
          {
              u1old = u1;
              u2old = u2;
              NavierStokes1;           
          }
              
          else if (choix == 2)
          {
              u1old = u1;
              u2old = u2;
              NavierStokes2;
          }
          if (compteur % 4 == 0)
          {
            plot([u1,u2],p, cmm="t = "+t, value=1, ps=colour[compteur]);
          }
          i++;
          t += dt;
          //compteur++;
      }
      plot([u1,u2],p, cmm="t = "+t, value=1);

Simulation Parameters

Fluide Velocity