Uzawa scheme for stokes bingham flow

Alice · May 9, 2024, 11:10am

Hi all,
I have written uzawa scheme for Bingham stokes flow on freefem. It is supposed to change when the bn number changes. However when I compile with different bn numbers it does not change. Could anyone help me?
Here is the code
real Bn = 10.;real Mu=1.;

real r = 1.; real r2=5e-6;
real L = 1; //Pipe length
real D = 1.; //Pipe height
real L1=.1;
real L2=.9;
real rho=1.; real alpha=1.;
macro Gradient(u) [dx(u), dy(u)] //
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //
macro UgradV(ux,uy,vx,vy) [ [ux,uy]'[dx(vx),dy(vx)] , [ux,uy]'[dx(vy),dy(vy)] ]//
macro StrainRate(ux,uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))] //

func f1=1.;
func f2=1.;
border b1(t=0., 1.){x=t; y=0.;};
border b2(t=0., 1.){x=L; y=Dt; };
border b3(t=L, L2){x=t; y=D;};
border b4(t=L2, L1){x=t; y=D;};
border b5(t=L1, 0.){x=t; y=D;};
border b6(t=0., 1.){x=0.; y=D-Dt;};
int nnL = 50;
int nnD = 50;
mesh Th=buildmesh(b1(nnL)+b2(nnD)+b3(5)+b4(5)+b5(5)+b6(nnD));
//Th = adaptmesh(Th, 1./128., IsMetric=1, nbvx=10000);
fespace Wh(Th,P2);
fespace Qh(Th,P1);
fespace Xh(Th,P2);

Wh vx,vy,ux1,uy1,dux,duy;
Wh ux,uy;
Qh p, q, dp,p1;

fespace Lh(Th,P1);
Lh TmGxx, TmGyy, TmGxy;
Lh Txx, Tyy, Txy;
Lh Gdotxx, Gdotyy, Gdotxy;
Lh dTxx=0., dTyy=0., dTxy=0.;
Lh TmGxx1, TmGyy1, TmGxy1;
Lh TmGxx0, TmGyy0, TmGxy0;
func u5= 1.;
func u4=1.;
func u3=1.;

TmGxx=0.;
TmGyy=0.;

TmGxy=0.;

for(int j=0; j<100; ++j)
{

solve bing([ux,uy],[vx,vy])
=int2d(Th)((Mu)*(Gradient(ux)‘*Gradient(vx)
+ Gradient(uy)’*Gradient(vy)
)
)

+int2d(Th)( Bn*[TmGxx,TmGyy,TmGxy]'* StrainRate(vx,vy)
)
 -int2d(Th)([f1 , f2]' *[vx ,vy]
	)
	
	+ on(b1, ux=0., uy=0.)
			+ on(b2, ux=0., uy=0.)
            + on(b3, ux=u3, uy=0.)
			+ on(b4, ux=u4, uy=0.)
			+ on(b5, ux=u5, uy=0.);

Gdotxx = 2.*dx(ux);
Gdotyy = 2.*dy(uy);
Gdotxy = dx(uy)+dy(ux);

for(int k=0; k<Lh.ndof; ++k){
	real TaGxx=TmGxx[][k]+r2*Gdotxx[][k],
		 TaGyy=TmGyy[][k]+r2*Gdotyy[][k],
		 TaGxy=TmGxy[][k]+r2*Gdotxy[][k];
	real TaGnorm = sqrt( .5*(TaGxx^2+TaGyy^2+2.*TaGxy^2) );
	real res = TaGnorm;
	if( 1.<res ){
		
		real q = 1./res;
		
		TmGxx[][k] = TaGxx*q;
		TmGyy[][k] = TaGyy*q;
		TmGxy[][k] = TaGxy*q;
	}
	else {
		
		TmGxx[][k] = TaGxx;
		TmGyy[][k] = TaGyy;
		TmGxy[][k] = TaGxy;
	}
}ux0[]=ux[];
uy0[]=uy[];

}
fespace Rh(Th, P2);
Rh dyg,dyg2;
solve streamlines (dyg, dyg2,solver=UMFPACK)
= int2d(Th)(
dx(dyg)*dx(dyg2)
+ dy(dyg)dy(dyg2)
)
+ int2d(Th)(
- dyg2(dy(ux) - dx(uy))
)
+ on(1,2,3,4,5,6, dyg=0);
plot(dyg,ps=“ldsre1500.eps”,nbiso=10);
Here is the uzawa bingham scheme
laplace u^n+1+divergence lamda^n=f
lamdaa^n+1=P(lamda^n+D(n^n+1))
P(q)=q if |q|<1 p(q)=q/|q| if |q|>1

fb77 · May 9, 2024, 11:45am

Can you please put all your line codes between:
three backquotes at the begining,
three backquotes at the end
so that it appears unaltered to us?

Alice · May 9, 2024, 12:17pm

Thank you has it changed?
Here is the code

‘’’
real Bn = 10.;real Mu=1.;

real r = 1.; real r2=5e-6;
real L = 1; //Pipe length
real D = 1.; //Pipe height
real L1=.1;
real L2=.9;
real rho=1.; real alpha=1.;
macro Gradient(u) [dx(u), dy(u)] //
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //
macro UgradV(ux,uy,vx,vy) [ [ux,uy]'[dx(vx),dy(vx)] , [ux,uy]'[dx(vy),dy(vy)] ]//
macro StrainRate(ux,uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))] //

func f1=1.;
func f2=1.;
border b1(t=0., 1.){x=t; y=0.;};
border b2(t=0., 1.){x=L; y=Dt; };
border b3(t=L, L2){x=t; y=D;};
border b4(t=L2, L1){x=t; y=D;};
border b5(t=L1, 0.){x=t; y=D;};
border b6(t=0., 1.){x=0.; y=D-Dt;};
int nnL = 50;
int nnD = 50;
mesh Th=buildmesh(b1(nnL)+b2(nnD)+b3(5)+b4(5)+b5(5)+b6(nnD));
//Th = adaptmesh(Th, 1./128., IsMetric=1, nbvx=10000);
fespace Wh(Th,P2);
fespace Qh(Th,P1);
fespace Xh(Th,P2);

Wh vx,vy,ux1,uy1,dux,duy;
Wh ux,uy;
Qh p, q, dp,p1;

fespace Lh(Th,P1);
Lh TmGxx, TmGyy, TmGxy;
Lh Txx, Tyy, Txy;
Lh Gdotxx, Gdotyy, Gdotxy;
Lh dTxx=0., dTyy=0., dTxy=0.;
Lh TmGxx1, TmGyy1, TmGxy1;
Lh TmGxx0, TmGyy0, TmGxy0;
func u5= 1.;
func u4=1.;
func u3=1.;

TmGxx=0.;
TmGyy=0.;

TmGxy=0.;

for(int j=0; j<100; ++j)
{

solve bing([ux,uy],[vx,vy])
=int2d(Th)((Mu)*(Gradient(ux)‘*Gradient(vx)
+ Gradient(uy)’*Gradient(vy)
)
)

+int2d(Th)( Bn*[TmGxx,TmGyy,TmGxy]'* StrainRate(vx,vy)
)
 -int2d(Th)([f1 , f2]' *[vx ,vy]
	)
	
	+ on(b1, ux=0., uy=0.)
			+ on(b2, ux=0., uy=0.)
            + on(b3, ux=u3, uy=0.)
			+ on(b4, ux=u4, uy=0.)
			+ on(b5, ux=u5, uy=0.);

Gdotxx = 2.*dx(ux);
Gdotyy = 2.*dy(uy);
Gdotxy = dx(uy)+dy(ux);

for(int k=0; k<Lh.ndof; ++k){
	real TaGxx=TmGxx[][k]+r2*Gdotxx[][k],
		 TaGyy=TmGyy[][k]+r2*Gdotyy[][k],
		 TaGxy=TmGxy[][k]+r2*Gdotxy[][k];
	real TaGnorm = sqrt( .5*(TaGxx^2+TaGyy^2+2.*TaGxy^2) );
	real res = TaGnorm;
	if( 1.<res ){
		
		real q = 1./res;
		
		TmGxx[][k] = TaGxx*q;
		TmGyy[][k] = TaGyy*q;
		TmGxy[][k] = TaGxy*q;
	}
	else {
		
		TmGxx[][k] = TaGxx;
		TmGyy[][k] = TaGyy;
		TmGxy[][k] = TaGxy;
	}
}ux0[]=ux[];
uy0[]=uy[];

}
fespace Rh(Th, P2);
Rh dyg,dyg2;
solve streamlines (dyg, dyg2,solver=UMFPACK)
= int2d(Th)(
dx(dyg)*dx(dyg2)
+ dy(dyg)dy(dyg2)
)
+ int2d(Th)(
- dyg2(dy(ux) - dx(uy))
)
+ on(1,2,3,4,5,6, dyg=0);
plot(dyg,ps=“ldsre1500.eps”,nbiso=10);

‘’’

Here is the uzawa bingham scheme
laplace u^n+1+divergence lamda^n=f
lamdaa^n+1=P(lamda^n+D(n^n+1))
P(q)=q if |q|<1 p(q)=q/|q| if |q|>1

Alice · May 9, 2024, 12:20pm

======Code=====================
real Bn = 10.;real Mu=1.;

real r = 1.; real r2=5e-6;
real L = 1; //Pipe length
real D = 1.; //Pipe height
real L1=.1;
real L2=.9;
real rho=1.; real alpha=1.;
macro Gradient(u) [dx(u), dy(u)] //
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //
macro UgradV(ux,uy,vx,vy) [ [ux,uy]'[dx(vx),dy(vx)] , [ux,uy]'[dx(vy),dy(vy)] ]//
macro StrainRate(ux,uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))] //

func f1=1.;
func f2=1.;
border b1(t=0., 1.){x=t; y=0.;};
border b2(t=0., 1.){x=L; y=Dt; };
border b3(t=L, L2){x=t; y=D;};
border b4(t=L2, L1){x=t; y=D;};
border b5(t=L1, 0.){x=t; y=D;};
border b6(t=0., 1.){x=0.; y=D-Dt;};
int nnL = 50;
int nnD = 50;
mesh Th=buildmesh(b1(nnL)+b2(nnD)+b3(5)+b4(5)+b5(5)+b6( nnD));
//Th = adaptmesh(Th, 1./128., IsMetric=1, nbvx=10000);
fespace Wh(Th,P2);
fespace Qh(Th,P1);
fespace Xh(Th,P2);

Wh vx,vy,ux1,uy1,dux,duy;
Wh ux,uy;
Qh p, q, dp,p1;

fespace Lh(Th,P1);
Lh TmGxx, TmGyy, TmGxy;
Lh Txx, Tyy, Txy;
Lh Gdotxx, Gdotyy, Gdotxy;
Lh dTxx=0., dTyy=0., dTxy=0.;
Lh TmGxx1, TmGyy1, TmGxy1;
Lh TmGxx0, TmGyy0, TmGxy0;
func u5= 1.;
func u4=1.;
func u3=1.;

TmGxx=0.;
TmGyy=0.;

TmGxy=0.;

for(int j=0; j<100; ++j)
{
solve bing([ux,uy],[vx,vy])
=int2d(Th)((Mu)*(Gradient(ux)'*Gradient(vx)

Gradient(uy)'*Gradient(vy)
)
)

+int2d(Th)( Bn*[TmGxx,TmGyy,TmGxy]‘* StrainRate(vx,vy)
)
-int2d(Th)([f1 , f2]’ *[vx ,vy]
)

on(b1, ux=0., uy=0.)
on(b2, ux=0., uy=0.)
on(b3, ux=u3, uy=0.)
on(b4, ux=u4, uy=0.)
on(b5, ux=u5, uy=0.);

Gdotxx = 2.*dx(ux);
Gdotyy = 2.*dy(uy);
Gdotxy = dx(uy)+dy(ux);

for(int k=0; k<Lh.ndof; ++k){
real TaGxx=TmGxx[k]+r2Gdotxx[][k],
TaGyy=TmGyy[][k]+r2Gdotyy[k],
TaGxy=TmGxy[k]+r2Gdotxy[][k];
real TaGnorm = sqrt( .5(TaGxx^2+TaGyy^2+2.*TaGxy^2) );
real res = TaGnorm;
if( 1.<res ){

real q = 1./res;

TmGxx[k] = TaGxxq;
TmGyy[x][k] = TaGyyq;
TmGxy[k] = TaGxy*q;
}
else {

TmGxx[k] = TaGxx;
TmGyy[k] = TaGyy;
TmGxy[k] = TaGxy;
}
}ux0=ux;
uy0=uy;

}
fespace Rh(Th, P2);
Rh dyg,dyg2;
solve streamlines (dyg, dyg2,solver=UMFPACK)
= int2d(Th)(
dx(dyg)*dx(dyg2)

dy(dyg)*dy(dyg2)
)
int2d(Th)(

dyg2*(dy(ux) - dx(uy))
)

on(1,2,3,4,5,6, dyg=0);
================end of the code=====================

fb77 · May 9, 2024, 1:10pm

Sorry it appears still modified by the text interpreter. You have to use only back quotes (you can copy/paste this one ` if you are not sure of the character on your keyboard).
Normally while composing your post you can see the result after the text interpreter in a separate right window, in which you can check what we will see.

Alice · May 9, 2024, 1:51pm

I hope it can be seen.
Thanks


real Bn = 10.;real Mu=1.;

real r = 1.; real r2=5e-6;
real L = 1;	//Pipe length
real D = 1.;	//Pipe height
real L1=.1;
real L2=.9;
real rho=1.; real alpha=1.;
macro Gradient(u) [dx(u), dy(u)] //
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //
macro UgradV(ux,uy,vx,vy) [ [ux,uy]'*[dx(vx),dy(vx)] , [ux,uy]'*[dx(vy),dy(vy)] ]//
macro StrainRate(ux,uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))] //

func f1=1.;
func f2=1.;
border b1(t=0., 1.){x=t; y=0.;};
border b2(t=0., 1.){x=L; y=D*t; };
border b3(t=L, L2){x=t; y=D;};
border b4(t=L2, L1){x=t; y=D;};
border b5(t=L1, 0.){x=t; y=D;};
border b6(t=0., 1.){x=0.; y=D-D*t;};
int nnL = 50;
int nnD = 50;
mesh Th=buildmesh(b1(nnL)+b2(nnD)+b3(5)+b4(5)+b5(5)+b6(nnD));
//Th = adaptmesh(Th, 1./128., IsMetric=1, nbvx=10000);
fespace Wh(Th,P2);
fespace Qh(Th,P1);
fespace Xh(Th,P2);

Wh vx,vy,ux1,uy1,dux,duy;
Wh ux,uy;
Qh p, q, dp,p1;

fespace Lh(Th,P1);
Lh TmGxx, TmGyy, TmGxy;
Lh Txx, Tyy, Txy;
Lh Gdotxx, Gdotyy, Gdotxy;
Lh dTxx=0., dTyy=0., dTxy=0.;
Lh TmGxx1, TmGyy1, TmGxy1;
Lh TmGxx0, TmGyy0, TmGxy0;
func u5= 1.;
func u4=1.;
func u3=1.;

	
TmGxx=0.; 
	TmGyy=0.;
	
	TmGxy=0.;
	
	for(int j=0; j<100; ++j)
	{  
solve bing([ux,uy],[vx,vy])
=int2d(Th)((Mu)*(Gradient(ux)'*Gradient(vx)
               +   Gradient(uy)'*Gradient(vy)
			   )
			  )
		
	+int2d(Th)( Bn*[TmGxx,TmGyy,TmGxy]'* StrainRate(vx,vy)
	)
	 -int2d(Th)([f1 , f2]' *[vx ,vy]
		)
		
		+ on(b1, ux=0., uy=0.)
				+ on(b2, ux=0., uy=0.)
                + on(b3, ux=u3, uy=0.)
				+ on(b4, ux=u4, uy=0.)
				+ on(b5, ux=u5, uy=0.);

	Gdotxx = 2.*dx(ux);
	Gdotyy = 2.*dy(uy);
	Gdotxy = dx(uy)+dy(ux);

	for(int k=0; k<Lh.ndof; ++k){
		real TaGxx=TmGxx[][k]+r2*Gdotxx[][k],
			 TaGyy=TmGyy[][k]+r2*Gdotyy[][k],
			 TaGxy=TmGxy[][k]+r2*Gdotxy[][k];
		real TaGnorm = sqrt( .5*(TaGxx^2+TaGyy^2+2.*TaGxy^2) );
		real res = TaGnorm;
		if( 1.<res ){
			
			real q = 1./res;
			
			TmGxx[][k] = TaGxx*q;
			TmGyy[][k] = TaGyy*q;
			TmGxy[][k] = TaGxy*q;
		}
		else {
			
			TmGxx[][k] = TaGxx;
			TmGyy[][k] = TaGyy;
			TmGxy[][k] = TaGxy;
		}
	}ux0[]=ux[];
	uy0[]=uy[];
	
	
	
}
fespace Rh(Th, P2);
Rh dyg,dyg2;
solve streamlines (dyg, dyg2,solver=UMFPACK)
    = int2d(Th)(
          dx(dyg)*dx(dyg2)
        + dy(dyg)*dy(dyg2)
    )
    + int2d(Th)(
        - dyg2*(dy(ux) - dx(uy))
    )
    + on(1,2,3,4,5,6, dyg=0);
	plot(dyg,ps="ldsre1500.eps",nbiso=10);

Here is the algoritm

fb77 · May 9, 2024, 2:13pm

Great, it is clear now.

I find two minor problems:

The border b4 should be discretized as b4(40) instead of b4(5)
You strain rate is not normalized correctly, it should be
macro StrainRate(ux,uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))/sqrt(2.)]
in order for the scalar product
[TmGxx,TmGyy,TmGxy]'* StrainRate(vx,vy)
to give the correct value, and accordingly
Gdotxy = 2.*(dx(uy)+dy(ux))/sqrt(2.);
real TaGnorm = sqrt( .5*(TaGxx^2+TaGyy^2+TaGxy^2) );

I think that also you should take the strain rate space Lh as P1dc instead of P1, so that Grad(u) belongs to Lh when u is in P2.

But I think that the main problem is that your parameter r2 is much too small. This implies that the convergence of u in the iteration is very bad (u appears in the strain rate update only behind the factor r2).
You should take r2 of order unity. But even in that case I am not sure that the iteration converges, in any case it could take thousands of iterations. You could try to add a relaxation term in (u^{k+1}-u^k) in the velocity equation (k is the iteration index), with a factor of order r2/h^2, see the PhD thesis of Duc Nguyen Hoai Theorem 4.2.6 p.71
https://theses.hal.science/tel-03078670/

Alice · May 9, 2024, 6:55pm

Thank you so much. I will correct the code what you have written.

Alice · May 10, 2024, 6:07am

Dear François Bouchut
I did the things you wrote but the streamlines don’t change when the Bn number changes. I even added the relaxation term.

fb77 · May 10, 2024, 10:56am

I think that the iteration is far from converged after your 100 iterations. You shoud check the difference ||u^{k+1}-u^k|| for each level k to have an idea. The method is known to be very slow. Later I will propose you a slightly different scheme with faster convergence.

Alice · May 10, 2024, 11:21am

Dear Professor I checked the difference it starts with 9e-5 ends as 1e-3 at the 100th iteration.

fb77 · May 10, 2024, 12:10pm

The first iteration is meaningless, but 1e-3 at the end looks not bad. Can you send me your code, I will look at it if I have time this afternoon.

Alice · May 10, 2024, 12:21pm

Dear Professor ,
Here is the latest code.

real Bn = 1.;real Mu=1.;

real r = 1.; 
real L = 1;	//Pipe length
real D = 1.;	//Pipe height
real L1=.1;
real L2=.9;
real rho=.1; real alpha=1.;real r2=10.;
macro Gradient(u) [dx(u), dy(u)] //
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //
macro UgradV(ux,uy,vx,vy) [ [ux,uy]'*[dx(vx),dy(vx)] , [ux,uy]'*[dx(vy),dy(vy)] ]//
macro StrainRate(ux,uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))/sqrt(2.)] //

func f1=1e-1;
func f2=1e-1;
border b1(t=0., 1.){x=t; y=0.;};
border b2(t=0., 1.){x=L; y=D*t; };
border b3(t=L, L2){x=t; y=D;};
border b4(t=L2, L1){x=t; y=D;};
border b5(t=L1, 0.){x=t; y=D;};
border b6(t=0., 1.){x=0.; y=D-D*t;};
int nnL = 32;
int nnD = 32;
mesh Th=buildmesh(b1(nnL)+b2(nnD)+b3(4)+b4(24)+b5(4)+b6(nnD));
//Th = adaptmesh(Th, 1./128., IsMetric=1, nbvx=10000);
fespace Wh(Th,P2);
fespace Qh(Th,P1);
fespace Xh(Th,P2);

Wh vx,vy,ux1,uy1,dux,duy;
Wh ux,uy,ux0,uy0;
Qh p, q, dp,p1;

fespace Lh(Th,P1dc);
Lh TmGxx, TmGyy, TmGxy;
Lh Txx, Tyy, Txy;
Lh Gdotxx, Gdotyy, Gdotxy;
Lh dTxx=0., dTyy=0., dTxy=0.;
Lh TmGxx1, TmGyy1, TmGxy1;
Lh TmGxx0, TmGyy0, TmGxy0;
func u5= 16.;
func u4=16.;
func u3=16.;

	
TmGxx=0.; 
	TmGyy=0.;
	
	TmGxy=0.;
	ux0=0.; uy0=0.;
	for(int j=0; j<100; ++j)
	{  
solve bing([ux,uy],[vx,vy])
=int2d(Th)((Mu+1./rho)*(Gradient(ux)'*Gradient(vx)
               +   Gradient(uy)'*Gradient(vy)
			   )
			  )
		-int2d(Th)((1./rho)*(Gradient(ux0)'*Gradient(vx)
               +   Gradient(uy0)'*Gradient(vy)
			   )
			  )  
	+int2d(Th)( [TmGxx,TmGyy,TmGxy]'* StrainRate(vx,vy)
	)
	 -int2d(Th)([f1 , f2]' *[vx ,vy]
		)
		
		+ on(b1, ux=0., uy=0.)
				+ on(b2, ux=0., uy=0.)
                + on(b3, ux=u3, uy=0.)
				+ on(b4, ux=u4, uy=0.)
				+ on(b5, ux=u5, uy=0.);

	Gdotxx = 2.*dx(ux);
	Gdotyy = 2.*dy(uy);
	Gdotxy = sqrt(2.)*(dx(uy)+dy(ux));

	for(int k=0; k<Lh.ndof; ++k){
		real TaGxx=TmGxx[][k]+r2*Gdotxx[][k],
			 TaGyy=TmGyy[][k]+r2*Gdotyy[][k],
			 TaGxy=TmGxy[][k]+r2*Gdotxy[][k];
		real TaGnorm = sqrt( .5*(TaGxx^2+TaGyy^2+TaGxy^2) );
		real res = TaGnorm;
		if( 1.<res ){
			
			real q = 1./res;
			
			TmGxx[][k] = TaGxx*q;
			TmGyy[][k] = TaGyy*q;
			TmGxy[][k] = TaGxy*q;
		}
		else {
			
			TmGxx[][k] = TaGxx;
			TmGyy[][k] = TaGyy;
			TmGxy[][k] = TaGxy;
		}
	}
	dux[] = ux[] - ux0[];
   duy[] = uy[] - uy0[];
   real err = sqrt(int2d(Th)(square(dux)+square(duy)));
   cout << err << "--iteration error  /iteration step  --" << j<< endl;
   ux0[]=ux[];
	uy0[]=uy[];
}
fespace Rh(Th, P2);
Rh dyg,dyg2;
solve streamlines (dyg, dyg2,solver=UMFPACK)
    = int2d(Th)(
          dx(dyg)*dx(dyg2)
        + dy(dyg)*dy(dyg2)
    )
    + int2d(Th)(
        - dyg2*(dy(ux) - dx(uy))
    )
    + on(1,2,3,4,5,6, dyg=0);
	plot(dyg,ps="ldsre1500.eps",nbiso=10);

fb77 · May 10, 2024, 6:11pm

Alice,
The results obtained with your code look correct (but take care that in the last version you sent, you forget the factor Bn in front of [TmGxx,TmGyy,TmGxy]'* StrainRate(vx,vy)). Indeed for Bn=1. or Bn=10. it converges quite fast, this is because the solution is smooth, there is no plug zone (zone where strainrate=0).
For Bn=100 it converges very slowly (and 100 iterations are not enough), this is because there is a plug zone.
I remarked that you put Neumann BC on boundary b6. Is this intentional or is this that you forgot to put ux=uy=0?

Otherwise there is a way to have faster convergence. It works if you consider in your equation 2\mathop{\rm div}(\mathop{\rm StrainRate})(u) instead of \Delta u. It follows that you get
int2d(Th)(2.*Mu*(StrainRate(ux,uy)'*StrainRate(vx,vy)))
instead of Grad'*Grad previously.
Then if you do not put the relaxation terms and if you choose the optimal coefficient r2=Mu/Bn, you get faster convergence.

Alice · May 11, 2024, 11:32am

Dear Professor,
Thank you so much. I updated the code according to what you mentioned. Here is the latest code.
The results for lid driven cavity squeezes to upward as bn increases as what is required.
However I when bn is lower such as 1 and the force changes like in this paper the result does not match with the paper “Analysis of semi-staggered finite-difference method with application
to Bingham flows ,Maxim A. Olshanskii”. Here mu=1,force=0.,u=1 when y=1.
How can I fix?

real Bn = 1.; real Mu=1.;

real r = 1.; 
real L = 1;	//Pipe length
real D = 1.;	//Pipe height
real L1=.1;
real L2=.9;
real rho=1.; real alpha=1.;real r2=Mu/Bn;
macro Gradient(u) [dx(u), dy(u)] //
macro Divergence(ux, uy) (dx(ux) + dy(uy)) //
macro UgradV(ux,uy,vx,vy) [ [ux,uy]'*[dx(vx),dy(vx)] , [ux,uy]'*[dx(vy),dy(vy)] ]//
macro StrainRate(ux,uy) [dx(ux), dy(uy), (dy(ux)+dx(uy))/sqrt(2.)] //

func f1=1.;
func f2=1.;
border b1(t=0., 1.){x=t; y=0.;};
border b2(t=0., 1.){x=L; y=D*t; };
border b3(t=L, 0.){x=t; y=D;};

border b4(t=0., 1.){x=0.; y=D-D*t;};
int nn = 50;

mesh Th=buildmesh(b1(nn)+b2(nn)+b3(nn)+b4(nn));
Th = adaptmesh(Th, 1./64., IsMetric=1, nbvx=10000);
fespace Wh(Th,P2);
fespace Qh(Th,P1);
fespace Xh(Th,P2);

Wh vx,vy,ux1,uy1,dux,duy;
Wh ux,uy,ux0,uy0;
Qh p, q, dp,p1;

fespace Lh(Th,P1dc);
Lh TmGxx, TmGyy, TmGxy;
Lh Txx, Tyy, Txy;
Lh Gdotxx, Gdotyy, Gdotxy;
Lh dTxx=0., dTyy=0., dTxy=0.;
Lh TmGxx1, TmGyy1, TmGxy1;
Lh TmGxx0, TmGyy0, TmGxy0;

func u3=1.;

	
TmGxx=0.; 
	TmGyy=0.;
	
	TmGxy=0.;
	ux0=0.; uy0=0.;
	for(int j=0; j<10; ++j)
	{  
solve bing([ux,uy],[vx,vy])
=int2d(Th)(2.*Mu*(StrainRate(ux,uy)'*StrainRate(vx,vy)))



	+int2d(Th)( Bn*[TmGxx,TmGyy,TmGxy]'* StrainRate(vx,vy)
	)
	 -int2d(Th)([f1 , f2]' *[vx ,vy]
		)
		
		+ on(b1, ux=0., uy=0.)
				+ on(b2, ux=0., uy=0.)
                + on(b3, ux=u3, uy=0.)
				+ on(b4, ux=0., uy=0.);
				

	Gdotxx = 2.*dx(ux);
	Gdotyy = 2.*dy(uy);
	Gdotxy = sqrt(2.)*(dx(uy)+dy(ux));

	for(int k=0; k<Lh.ndof; ++k){
		real TaGxx=TmGxx[][k]+r2*Gdotxx[][k],
			 TaGyy=TmGyy[][k]+r2*Gdotyy[][k],
			 TaGxy=TmGxy[][k]+r2*Gdotxy[][k];
		real TaGnorm = sqrt( .5*(TaGxx^2+TaGyy^2+TaGxy^2) );
		real res = TaGnorm;
		if( 1.<res ){
			
			real q = 1./res;
			
			TmGxx[][k] = TaGxx*q;
			TmGyy[][k] = TaGyy*q;
			TmGxy[][k] = TaGxy*q;
		}
		else {
			
			TmGxx[][k] = TaGxx;
			TmGyy[][k] = TaGyy;
			TmGxy[][k] = TaGxy;
		}
	}
	dux[] = ux[] - ux0[];
   duy[] = uy[] - uy0[];
   
   real err = sqrt(int2d(Th)(square(dux)+square(duy)));
   cout << err << "--iteration error  /iteration step  --" << j<< endl;
   ux0[]=ux[];
	uy0[]=uy[];
}
fespace Rh(Th, P2);
Rh dyg,dyg2;
solve streamlines (dyg, dyg2,solver=UMFPACK)
    = int2d(Th)(
          dx(dyg)*dx(dyg2)
        + dy(dyg)*dy(dyg2)
    )
    + int2d(Th)(
        - dyg2*(dy(ux) - dx(uy))
    )
    + on(1,2,3,4,5,6, dyg=0);
	plot(dyg,ps="ldsre1500.eps",nbiso=10);

fb77 · May 13, 2024, 1:15pm

In that paper the flow is taken incompressible, whereas in your code you use compressible. This is why the solution is different (but looks correct).
You can modify easily your code to do the incompressible case, by using the extra variables p,q as
solve bing([ux,uy,p],[vx,vy,q])
and adding
-int2d(Th)(p*Divergence(vx, vy)+q*Divergence(ux, uy))
But then it happens that the scheme converges extremely slowly (you need thousands of iterations). In order to improve it you should use the Augmented Lagrangian method as in that paper.

Another point is that your computation of the stream function could be misleading. You look for \psi satisfying
-\Delta\psi=\partial_y u_x-\partial_x u_y.
An important property that one would like that if \partial_x u_x+\partial_y u_y=0 then \partial_x\psi=u_y and -\partial_y\psi=u_x.
You see then that the Dirichlet condition on \psi is not appropriate. It is better to take
\frac{\partial\psi}{\partial n}=u_y n_x-u_x n_y on \partial\Omega.
We are thus led to the problem
-\Delta \psi=f in \Omega,\qquad\frac{\partial\psi}{\partial n}=g on \partial\Omega.
To solve this system you need the compatibility condition (to see that it is necessary, integrate the equation over \Omega)
\int_\Omega f=-\int_{\partial\Omega}g
Assuming this, the system is not invertible since \psi is defined up to a constant. You can determine it by imposing \int_{\partial\Omega}\psi=0.
In order to get the solution you can then solve
-\Delta \psi=f in \Omega,\qquad\frac{\partial\psi}{\partial n}+\varepsilon\psi=g on \partial\Omega,
for some small \varepsilon. This system is well-posed. This leads to the code

fespace Rh(Th, P2);
Rh dyg,dyg2;
solve streamlines (dyg, dyg2,solver=UMFPACK)
    = int2d(Th)(
          dx(dyg)*dx(dyg2)
        + dy(dyg)*dy(dyg2)
    )
    + int2d(Th)(
        - dyg2*(dy(ux) - dx(uy))
    )
    + int1d(Th)(
        - dyg2*(uy*N.x-ux*N.y)
    )
    + int1d(Th)(1e-8*dyg*dyg2)
    ;

Alice · May 14, 2024, 8:54am

Dear Professor,
Thank you so much for your valuable precious support. When I changed the streamfunction as you wrote at the last message,the graph looks like this.

Previously it looked like this
.What do you think that need to change?
Regards

fb77 · May 14, 2024, 10:32am

I think that the first (with non homogeneous Neumann BC) is the correct one.
But as far as div(u) is non zero, this is not much meaningful.

Alice · May 14, 2024, 11:00am

Dear Professor ,
Thank you so much. So for the incompressible case will the stream function algoritm change ?

fb77 · May 14, 2024, 11:16am

The formulation with Neumann BC is especially good for the incompressible case, since as I said it is not much meaningfull when div u is not zero.

Topic		Replies	Views
Navier Stokes with Newton method General Discussion	6	62	July 1, 2024
Sample Navier-Stokes Solver General Discussion	1	236	December 4, 2022
Newtons methods for unsteady problem General Discussion	9	523	November 9, 2022
Convect function failing General Discussion	2	378	May 12, 2022
Parallel fluid code General Discussion	10	640	June 9, 2021

Uzawa scheme for stokes bingham flow

Related Topics