Poiseuille Flow simulation of viscoplastic fluid

Hi all,
I’m writing a code to compute the Poiseuille flow of a voscoplastic fluid in the Stokes regime. However it’s giving me an error I can’t figure out. Could anyone help me?

the code is the following:

 load "iovtk"
//
//  Parameters:

 real u0 = 0.1;           	// initial guess on velocity field
 real pin = 1.;               	// input pressure constant wrt time

 real del = 0.20;	// ratio heigth/length of the channel, comes from the non-dimensionalization of the problem
 real Re = 0.5;          // Reynolds number
 real G = 1;             // Gamma and theta material parameters of the problem
 real th = 1;
 real k = th*Re;               // DKT number
 real Bi = G/(2*Re);           // Bingham number

 real dt = 0.0001;         // time step
 real eps = 1e-4;              // term for papanastasiou approximation
 real err=0;			// error for Newton-Rhapson
 real tol = 1e-6;

 real n = 20;           	// number of nodes

 int Wall = 1;		// pipe wall label
 int Inlet = 2;		// pipe inlet label
 int Outlet = 3;		// pipe outlet label


//  Boundary
 border b1 ( j = 0.0, 1.0 ) { x = 2*j-1; y = -1.0; label=Wall;}
 border b2 ( j = 0.0, 1.0 ) { x = 1.0; y = 2*j-1; label=Outlet;}
 border b3 ( j = 0.0, 1.0 ) { x = 1-2*j; y = 1.0; label=Wall;}
 border b4 ( j = 0.0, 1.0 ) { x = -1.0; y = 1-2*j; label=Inlet;}

//  Define the geometry of the problem and the finite element mesh
 mesh Mh = buildmesh ( b1(n) + b2(n) + b3(n) + b4(n));
  
  
// FE space Taylor Hood
 fespace Vh(Mh,P2);		//velocity field
 fespace Qh(Mh,P1);		//pressure
 Vh ux,uy,vx,vy,dux,duy,upx,upy;
 Qh p,q,dp,pp;
 Vh stress;
 Vh gamma,gd;
 Vh visc;
 Vh dvisc;
  
  
// Macros
 real del2 = del^2.;
 real del3 = del^3.;

 macro grad(ux,uy) [ dx(ux),dy(ux) , dx(uy),dy(uy) ]	// EOM
 macro UgradV(ux,uy,vx,vy) [ [ux,uy]'*[dx(vx),dy(vx)] , [ux,uy]'*[dx(vy),dy(vy)] ]	// EOM
 macro div(ux, uy) (dx(ux) + dy(uy)) // EOM
 macro strainrate(ux,uy) [del*dx(ux),(dy(ux)+del2*dx(uy))/2., (dy(ux)+del2*dx(uy))/2.,del*dy(uy)] // EOM


// Set initial conditions and other quantities, anche questo da cambiare
  
 real fin= 1. ;// normal force inlet  
 real fout = -1.;// normal force outlet 
 
 upx = 1.-y^4.;
 upy = 0.;

 ux = 1.-y^4.;
 uy = 0.;
 
 dux = 0.;
 duy = 0.;
 
 gamma = sqrt( .5*( del2*dx(upx)^2. + (dy(upx)+del2*dx(upy))^2. + del2*dy(upy)^2. ) );
 visc = Bi/(gamma+eps) + 2.*exp(-2.*k*gamma);
 dvisc = Bi/(gamma+eps)^2. + 2.*k*exp(-2.*k*gamma);
 

 // this works
 //gd =  dvisc*sqrt(abs(strainrate(ux,uy)'*strainrate(dux,duy)));
 //plot (gd, value=1, fill=1, dim=3, wait = 1); 
 
 //cout << "dvisc: " << dvisc << "\n";
  
 int maxit=50;  //  maximum number of iterations
 for (int i=0;i<=maxit;i++)		// Newton-Rhapson DF(u)(du) = -F(u)
 {
 	solve PoiseDKT([dux,duy,dp],[vx,vy,q]) = 
	int2d(Mh)( dvisc * sqrt(strainrate(ux,uy)'*strainrate(dux,duy)) 
		   * ( (strainrate(ux,uy) + strainrate(dux,duy))'*grad(vx,vy) )
		   + div(dux,duy)*q + div(vx,vy)*dp
		   - 1e-8*dp*q 					// stabilization term
		 )
	- int2d(Mh)( visc * strainrate(ux,uy)'*grad(vx,vy)
	 	      - div(ux,uy)*q - div(vx,vy)*p
	 	   )
 	+ on(Wall, dux = 0., duy = 0.) ;
 	
 	// update u and other quantities
 	upx[] = ux[];  upy[] = uy[];
 	ux[] -= dux[]; uy[] -= duy[]; p[] -= dp[];
 	
 	gamma = sqrt( .5*( del2*dx(upx)^2. + (dy(upx)+del2*dx(upy))^2. + del2*dy(upy)^2. ) );
 	visc = Bi/(gamma+eps) + 2.*exp(-2.*k*gamma);
 	dvisc = Bi/(gamma+eps)^2. + 2.*k*exp(-2.*k*gamma);
 	
	err= dux[].linfty + duy[].linfty + dp[].linfty;
	cout << " error = " << err << "\n";
	
	if(err < tol){
		cout << " converges |" <<" iter " << n << " |error = " << err << "\n";
		break;
	}
	if( n>3 && err > 10.){
		cout << " does not converge |" <<" iter " << n << " |error = " << err << "\n"; 
		break;
	}

 }

 stress = 2*exp(-2*k*gamma(ux,uy))*gamma(ux,uy)  + Bi*gamma(ux,uy) /(gamma(ux,uy)+eps);
 
 plot (ux, value=1, fill=1, dim=3, wait = 1);

 int[int] orderOut = [1, 1, 1, 1];
 savevtk("poiseuille_explicit.vtu", Mh, ux, uy, [ux, uy, 0], stress, dataname="ux uy vector stress", order=orderOut);


 cout << "\n";
 cout << "End of execution.\n";

I’m using a Newton-Rhapson iteration. Aside from the fairness of what I wrote in the weak form inside the iteration, I’m interested to understand why the error I receive is:

...

102 :  	solve PoiseDKT([dux,duy,dp],[vx,vy,q]) = 
  103 : 	int2d(Mh)( dvisc * sqrt(abs(strainrate(upx,upy)    [del*dx(upx),(dy(upx)+del2*dx(upy))/2., (dy(upx)+del2*dx(upy))/2.,del*dy(upy)] '*strainrate(dux,duy)    [del*dx(dux),(dy(dux)+del2*dx(duy))/2., (dy(dux)+del2*dx(duy))/2.,del*dy(duy)] ) error operator (  <10LinearCombI7MGauche4C_F0E> 
 List of choices 
	 (	  <5F_KN_IddddE> :   <3KN_IdE> )
	 (	  <d> :   <St7complexIdE> )
	 (	  <l> :   <l> )
	 (	  <d> :   <d> )

 Error line number 103, in file de_kee_stationary_newton.edp, before  token )

  current line = 103
Compile error : 
	line number :103, )
error Compile error : 
	line number :103, )
 code = 1 mpirank: 0

It looks like there’s a problem multiplying dvisc * sqrt(...). I suspect the problem lies in the fact that the square root is not a linear operation and inside the square root there is my increment [dux, duy] to be computed by the solve command. But I would like another opinion on the matter.

It is difficult to understand such error message, but you are right, the problem is that you must write a variational formulation which is linear with respect to the unknown. For Newton-Raphson, the formula DF(u)(du) = -F(u) is obviously linear with respect to du, hence you must revise your solve PoiseDKT to be coherent with the Newton-Raphson algorithm!
Moreover you must separate (i.e. put in different int2d()) terms which are proportional to the unknown and terms which do not depend on the unknown.

Thank you for your answer, I have made the following alterations to the code:

//
 load "iovtk"
//
//  Parameters:
 real pin = 1.;               	// input pressure constant
 
 real Re = 0.5;          	// Reynolds number
 real G = 1;             	// Gamma and theta material parameters of the problem
 real th = 15;
 real k = th*Re;               	// DKT number
 real Bi = G/(2*Re);           	// Bingham number

 real eps = 1e-4;              	// term for regularization of viscosity
 real err = 0;			// error for Newton-Rhapson
 real tol = 1e-6;

 real n = 20;           	// number of nodes

 int Wall = 1;			// pipe wall label
 int Inlet = 2;			// pipe inlet label
 int Outlet = 3;		// pipe outlet label


//  Boundary
  border b1 ( j = 0.0, 1.0 ) { x = 2*j-1; y = -1.0; label=Wall;}
  border b2 ( j = 0.0, 1.0 ) { x = 1.0; y = 2*j-1; label=Outlet;}
  border b3 ( j = 0.0, 1.0 ) { x = 1-2*j; y = 1.0; label=Wall;}
  border b4 ( j = 0.0, 1.0 ) { x = -1.0; y = 1-2*j; label=Inlet;}

//  Define the geometry of the problem and the finite element mesh
 mesh Mh = buildmesh ( b1(n) + b2(n) + b3(n) + b4(n));
 plot(Mh, wait = 1);
  
// FE space
 fespace Vh(Mh,P1);		//velocity field
 Vh ux,uy,vx,vy,dux,duy;
 Vh stress;
 Vh gamma,gammae;
 Vh visc;
 Vh dvisc;
  
// Macros
 macro Grad(ux,uy) [ dx(ux), dy(ux) , dx(uy), dy(uy) ]	// EOM
 macro div(ux, uy) (dx(ux) + dy(uy)) // EOM
 macro strainrate(ux,uy) [ dx(ux),(dy(ux)+dx(uy))/2. , (dy(ux)+dx(uy))/2.,dy(uy) ] // EOM
 
// Set initial conditions and other quantities
 ux = 1.-y^4.;
 uy = 0.;
 
 dux = 0;
 duy = 0;
 
 gamma = sqrt( .5*( dx(ux)^2. + (dy(ux)+dx(uy))^2. + dy(uy)^2. ) );
 gammae = sqrt( .5*( dx(ux)^2. + (dy(ux)+dx(uy))^2. + dy(uy)^2. ) + eps^2.);
 visc = Bi/(gammae) + 2.*exp(-2.*k*gamma);
 dvisc = Bi/(2.*(gammae)^3.) + 2.*k*exp(-2.*k*gamma)/gammae;

 
 // first an explicit step to feed the Newton loop
 solve init([ux,uy],[vx,vy]) = 
 int2d(Mh)( visc * (strainrate(ux,uy)'*Grad(vx,vy)) - [pin, 0]'*[vx, vy] )
 	    + on(Wall, ux = 0., uy = 0.) ;
 
 gamma = sqrt( .5*( dx(ux)^2. + (dy(ux)+dx(uy))^2. + dy(uy)^2. ) );
 gammae = sqrt( .5*( dx(ux)^2. + (dy(ux)+dx(uy))^2. + dy(uy)^2. ) + eps^2. );
 visc = Bi/gammae + 2.*exp(-2.*k*gamma);
 dvisc = Bi/(2.*(gammae)^3.) + 2.*k*exp(-2.*k*gamma)/gammae;
  
 int maxit=50;  //  maximum number of iterations
 for (int i=0;i<=maxit;i++)		// Newton-Rhapson DF(u)(du) = -F(u)
 {
 	solve PoiseDKT([dux,duy],[vx,vy]) = 
	int2d(Mh)( 
		   visc * (strainrate(dux,duy)'*Grad(vx,vy)) 
		   - dvisc * (strainrate(ux,uy)'*strainrate(dux,duy)) * (strainrate(ux,uy)'*Grad(vx,vy))
		 )
	- int2d(Mh)( visc * strainrate(ux,uy)'*Grad(vx,vy) - [pin, 0]'*[vx, vy] )
 	+ on(Wall, dux = 0., duy = 0.) ;
 	
 	// update u and other quantities
 	ux[] += dux[]; uy[] += duy[]; 

 	gamma = sqrt( .5*( dx(ux)^2. + (dy(ux)+dx(uy))^2. + dy(uy)^2. ) );
 	gammae = sqrt( .5*( dx(ux)^2. + (dy(ux)+dx(uy))^2. + dy(uy)^2. ) + eps^2. );
 	visc = Bi/gammae + 2.*exp(-2.*k*gamma);
 	dvisc = Bi/(2.*(gammae)^3.) + 2.*k*exp(-2.*k*gamma)/gammae;
 	
	err= dux[].linfty + duy[].linfty;
	cout << " error = " << err << "\n";
	
	if(err < tol){
		cout << " Converges |" <<" iter " << i << " |error = " << err << "\n";
		break;
	}
	if( i=50 && err > 10.){
		cout << " Does not converge |" << "error = " << err << "\n"; 
		break;
	}

 }
 
 plot([u1,u2],cmm="u1,u2 math ",wait=1);
 plot (ux, value=1, fill=1, dim=3, wait = 1);

I’m trying to run the code but I receive this error message from the terminal:

86 :  int2d(Mh)( visc * (strainrate(ux,uy)       [ dx(ux),(dy(ux)+dx(uy))/2. , (dy(ux)+dx(uy))/2.,dy(uy) ] '*Grad(vx,vy)       	 [ dx(vx), dy(vx) , dx(vy), dy(vy) ]	) - [pin, 0]'*[vx, vy] ) error operator -  <10LinearCombISt4pairI7MGauche6MDroitE4C_F0E>, <10LinearCombI6MDroit4C_F0E> 
 List of choices 
	 (	  <l> :   <l>, <l> )
.....

Then I believe it gives a long list of choices for the elements that can be chosen for the " ’ * " operator and at the end it confirms that the error is at line 86:

 Error line number 86, in file de_kee_stokes_newton.edp, before  token )

  current line = 86
Compile error : 
	line number :86, )
error Compile error : 
	line number :86, )
 code = 1 mpirank: 0

To clarify what I’m trying to do, my idea in writing the code is based on the following calculations:

Here D(u) is the symmetric part of the velocity gradient, and u is the velocity. I believe the formulas to be linear w.r.t. the increment du, because both the operator " ’ * " and the term D(du) are linear. Could the problem be that D(du) contains the derivatives of the increment du? Or am I using the syntax incorrectly somewhere?

As I warned you, you must put in different int2d terms which are proportional to the unknown and terms which do not depend on the unknown.
Here in solve init the term strainrate(ux,uy)'*Grad(vx,vy) depends on the unknown, whereas [pin, 0]'*[vx, vy] does not. You have to separate them.

Just a small remark in Newton-Rhapson you can find directly u^{n+1} because the boundary condition on du^{n+1} must be 0 and if you decompose F un 3 part : constant, linear , non nonlinear : F= C + L + NL , when DF(u)(du) = L(du) + DNL(u)(du) and then the Newtow Rhapton become
because we have u^{n+1} = u^{n} +du^{n+1} with DF(u^{n})(du^{n+1}) = -F(u^{n})
DF(u^n)(u^{n+1}) = DF(u^n)(u^{n}) +DF(u^n)(du^{n+1}) = DF(u^n)(u^{n}) -F(u^n)
and with the decomposition you have :slight_smile:
DF(u^n)(u^{n+1}) = DNL(u^n)(u^{n}) -C - NL(u^n)

More easy to code.