Navier Stokes on a step

General Discussion
1

Hey, everyone! I’m still a problem with solve Navier Stokes equation on a step. When I plot, the vector come out of the mesh (what should not happen). Could someone help me?

image
image1263×360 30.3 KB

(It stays like this)

Note: I’m using Newton’s Method to linearize

2

Without code, or any further info, it is hard to help. I recommend you use paraview to visualize the solution.

3

I managed to solve this plot in FreeFEM, but when I export the VTK to Paraview it looks different.

image
image1412×204 106 KB

while ((err > tol) && (iter < maxIter)){

  problem NavierStokes([u1,u2,p], [v1,v2,q])
    = int2d(Th)(
      mu*(dx(u1)*dx(v1) + dy(u1)*dy(v1)
        + dx(u2)*dx(v2) + dy(u2)*dy(v2))
      + (u1k*dx(u1) + u2k*dy(u1))*v1
      + (u1k*dx(u2) + u2k*dy(u2))*v2
      - p*(dx(v1) + dy(v2))
      - q*(dx(u1) + dy(u2))
    )
    - int2d(Th)(fx*v1 + fy*v2)
    + on(1, u1=uin, u2=0)
    + on(2, u1=0, u2=0);
  
  NavierStokes;

  // Verify error
  err = sqrt(int2d(Th)((u1 - u1k)^2 + (u2 - u2k)^2));
  u1k = u1;
  u2k = u2;

vtk:

// Cabeçalho do arquivo VTK
  vtkfile << "# vtk DataFile Version 3.0" << endl;
  vtkfile << "Malha P2 com dados P1/P2" << endl;
  vtkfile << "ASCII" << endl;
  vtkfile << "DATASET UNSTRUCTURED_GRID" << endl;

  // Escrever os pontos (vértices) da malha
  vtkfile << "POINTS " << Th.nv << " float" << endl;
  for (int i = 0; i < Th.nv; ++i) {
      vtkfile << Th(i).x << " " << Th(i).y << " 0.0" << endl;
  }

  // Escrever as células (triângulos) da malha
  vtkfile << "CELLS " << Th.nt << " " << 4 * Th.nt << endl;
  for (int i = 0; i < Th.nt; ++i) {
      // Para uma malha triangular: 3 vértices por célula
      vtkfile << "3 " << Th[i][0] << " " << Th[i][1] << " " << Th[i][2] << endl;
  }

  // Escrever os tipos de células (todas serão triângulos)
  vtkfile << "CELL_TYPES " << Th.nt << endl;
  for (int i = 0; i < Th.nt; ++i) {
      vtkfile << "5" << endl;  // Tipo "5" é para triângulos
  }

  // Escrever os dados do campo de velocidade
  vtkfile << "POINT_DATA " << Th.nv << endl;
  vtkfile << "VECTORS velocity float" << endl;
  for (int i = 0; i < Th.nv; ++i) {
      vtkfile << u1[][i] << " " << u2[][i] << " 0.0" << endl;
  }
4

Maybe you can use “iovtk” to save .vtu file, it is more simple.
It is like

load "iovtk"
int[int] fforder = [1,1];
string Dataname = "u";
savevtk("u.vtu", Th, [u1, u2], dataname=Dataname, order=fforder);