Startpage >> Main >> Bingham

Bingham

Solution to Bingham fluid

info pdf file

Thanks to Ali Roustaei arousta@alumni.ubc.ca

download example: Bingham.edp or return to 2D examples

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);

return to 2D examples

Page last modified on June 27, 2018, at 03:31 PM