Thanks a lot for the answer, i’m am new to freefem and i didn’t understand that non linear terms are a problem for freefem. Can i ask you one more thing?

After running a first simulation (without stabilization) with lower reynolds i tried to run a second one with higher Reynolds, so i can use the first one as initial condition for the second. I’m having problem with the SUPG stabilization terms. This term is non linear, but i’am trying some escamotage to solve the problem without linearizing. (driven cavity problem).

Do you have any advice or raccomandation (or see any errors) because the code does not converge.

fespace Uh(Th,P2)

Uh ux, uu;

Uh uy,vv;

fespace Ph(Th,P1)

Ph p,pp;

unastro=1; //reference velocity

cout<<“Reynolds =”<<Re<<endl;

//SUPG PARAMETERS

real delta=1; //parametro libero

Uh U=sqrt(ux^2+uy^2); //modulo velocità

Uh Reu=U* hTriangle/(2*nu); //ReynoldsU

Uh tau; //parametro di stabilizzazione

tau=(Reu<3)* hTriangle* hTriangle/(12*nu)+(Reu>=3)* hTriangle/(2*U);

Uh pold;

macro div(u1,u2) (dx(u1)+dy(u2)) //

macro laplaciano(u1) (dxx(u1)+dyy(u1)) //

macro convettivo(u1,u2,u3) (u1*dx(u3)+u2*dy(u3)) //

macro prodint(u1,u2,v1,v2) (dx(u1)*dx(v1)+dy(u2)*dy(v2)+0.5*(dx(u2)*dx(v2)+

dx(u2)*dy(v1)+dy(u1)*dx(v2)+dy(u1)*dy(v1))) //

problem NavierStokesSUPG ([ux,uy,p],[uu,vv,pp],init=i) =

int2d(Th)(alpha* (ux*uu + uy*vv)

+ 2* invRe* prodint(ux,uy,uu,vv)

- p* pp* (0.000001) - p* div(uu,vv)

+ pp* div(ux,uy))

- int2d(Th)(alpha* (uoldx* uu + uoldy* vv))

+ int2d(Th)(convettivo(uoldx,uoldy,ux)* uu //termini convettivi trattati

+ convettivo(uoldx,uoldy,uy)* vv) //in modo semi-implicito

// Stabilizzazione SUPG

+ int2d(Th)(delta* tau* convettivo(uoldx,uoldy,uu)* alpha* ux)

+ int2d(Th)(delta* tau* convettivo(uoldx,uoldy,vv)* alpha* uy)

- int2d(Th)(delta* tau* convettivo(ux,uy,uu)*(alpha* uoldx+dx(pold)-invRe* laplaciano(uoldx)))

- int2d(Th)(delta* tau* convettivo(ux,uy,vv)*(alpha* uoldy+dy(pold)-invRe* laplaciano(uoldy)))

+ int2d(Th)(delta* tau* convettivo(uoldx,uoldy,uu)* convettivo(uoldx,uoldy,ux))

+ int2d(Th)(delta* tau* convettivo(uoldx,uoldy,vv)* convettivo(uoldx,uoldy,uy))

+ on(1,2,4,ux=0,uy=0) + on(3,ux=unastro,uy=0);

pold=p;

uoldx=ux; //simulazione precedente presa come condizione iniziale

uoldy=uy;

NavierStokesSUPG;

Th=adaptmesh(Th,[ux,uy]);

plot(Th, wait=1);

plot([ux,uy], wait=1);

for(i=0; i<100; i++){

uoldx=ux;

uoldy=uy;

NavierStokesSUPG;

if (!(i % 2)) // plot every 2 iteration

plot(p,[ux,uy]);}

plot(p, wait=1);

plot([ux,uy], wait=1, value=true);