Bingham

Solution to Bingham fluid

info

Thanks to Ali Roustaei

real Bn = 4.8;
real dPdL = -9.6;
real r = 10.;

macro div(u) ( dx(u#x)+dy(u#y) ) //
macro StrainRate(u) [2.*dx(u#x), 2.*dy(u#y), dy(u#x)+dx(u#y)] //

mesh Th = square(3,20,[x,2.*y]);
fespace Wh(Th,[P2,P2,P1],periodic=[[2,y], [4,y]]);
fespace Qh(Th,P1);
fespace Xh(Th,P2);

Wh [ux,uy,p], [vx,vy,q];
fespace Lh(Th,P1dc);
Lh TmGxx=0., TmGyy=0., TmGxy=0.;
Lh Txx=0., Tyy=0., Txy=0.;
Lh Gdotxx=0., Gdotyy=0., Gdotxy=0.;

problem Poisse([ux,uy,p],[vx,vy,q],init=true) =
int2d(Th)(  r*(StrainRate(u)'*StrainRate(v))-p*div(v) -q*div(u) +1e-10*q*p )
+int2d(Th)( dPdL*vx )
+int2d(Th)( [TmGxx,TmGyy,TmGxy]'*[dx(vx),dy(vy),dy(vx)+dx(vy)] )
+on(1,3,ux=0,uy=0);

for(int iter=0; iter<100; ++iter)
{
Poisse;
Gdotxx = 2.*dx(ux);
Gdotyy = 2.*dy(uy);
Gdotxy = dx(uy)+dy(ux);

for(int k=0; k<Lh.ndof; ++k){
real TaGxx=Txx[][k]+r*Gdotxx[][k],
TaGyy=Tyy[][k]+r*Gdotyy[][k],
TaGxy=Txy[][k]+r*Gdotxy[][k];
real TaGnorm = sqrt( .5*(TaGxx^2+TaGyy^2+2.*TaGxy^2) );
real res = TaGnorm-Bn;
if( 0<res ){
real c = r/(1.+r)*res/TaGnorm;
real q = 1.-c;
Txx[][k] = TaGxx*q;
Tyy[][k] = TaGyy*q;
Txy[][k] = TaGxy*q;

q = 1.-2.*c;
TmGxx[][k] = TaGxx*q;
TmGyy[][k] = TaGyy*q;
TmGxy[][k] = TaGxy*q;
}
else {
Txx[][k] = TaGxx;
Tyy[][k] = TaGyy;
Txy[][k] = TaGxy;

TmGxx[][k] = TaGxx;
TmGyy[][k] = TaGyy;
TmGxy[][k] = TaGxy;
}
}
}
plot(ux,dim=3,fill=true,value=true);