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.