Problem related to SUPG

General Discussion
1

I am trying to solve the Navier Stokes equation in a driven cavity but i am acctually having problem with adding the stabilization term related to SUPG.
Freefam can only solve weak formulation problem and in the integral related to the SUPG terms we are adding the strond formulation of the problem:

fespace Uh(Th,P2);
Uh ux,uu;
Uh uy,vv;

fespace Ph(Th,P1);
Ph p,pp;

//stabilization term related to SUPG
int2d(Th)(deltatau(uxdx(uu)+uydy(vv))(uxdx(ux)…)

The problem is related to the term: ux*dx(ux)
Can someone help me or does someone have some example to look at?

2

Your stabilization term is non linear so need a schema to solve this problem.

3

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.

Immagine WhatsApp 2024-03-09 ore 19.38.27_9f20ac4d
Immagine WhatsApp 2024-03-09 ore 19.38.27_9f20ac4d869×162 14 KB

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/(12nu)+(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) (u1dx(u3)+u2dy(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* (uxuu + uyvv)
+ 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);

4

Hello Did you solve your problem ?
and if you have a Simple reference regarding SUPG implementation, could you please share it

5

I remember there are some scripts in Stabfem project with SUPG methods to solve the steady solution where the stabilization term is linearized for the jacobian matrix. For the time stepping problem, I think maybe the stabilization term could be added to the source terms in some senses.

1 Like
6

Thank you for you reply
but I really did not understand any thing from stabfem project
any recommendations for simpler explanation

7

could you please Professor explain more about your reply pr suggest a reference

8

Hi, there is one example with SUPG to solve the Navier-Stokes equations in stabfem project. They introduce one stabilization term and linearize it in the Newton iteration. Maybe this could give you some ideas to write your own code.