Thank you ,Professor.
Sorry, about the incompressible version I forgot one line to determine the additive constant in the pressure
+int2d(Th)(-1.e-10*p*q)
With that the scheme converges rather nicely! Here is the full code (with other minor modifications).
I think that Mu=sqrt(2.) corresponds to the test of Olshanskii.
verbosity=0;
real Bn = 1.; real Mu=sqrt(2.);
real L = 1; //Pipe length
real D = 1.; //Pipe height
real L1=.1;
real L2=.9;
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=0.;
func f2=0.;
border b1(t=0., 1.){x=t*L; 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 = 64;
mesh Th=buildmesh(b1(nn)+b2(nn)+b3(nn)+b4(nn));
//Th = adaptmesh(Th, 1./64., IsMetric=1, nbvx=10000);
plot(Th,wait=1);
fespace Wh(Th,P2);
fespace Qh(Th,P1);
Wh ux,uy,ux0,uy0,vx,vy;
Qh p,q;
Qh normsigma,normsigmaplot;
fespace Lh(Th,P1dc);
Lh TmGxx, TmGyy, TmGxy;
Lh Gdotxx, Gdotyy, Gdotxy;
Lh TmGxx0, TmGyy0, TmGxy0;
Lh Tauxx,Tauyy,Tauxy,Taunorm;
func u3=1.;
TmGxx=0.;
TmGyy=0.;
TmGxy=0.;
TmGxx0=TmGxx;TmGyy0=TmGyy;TmGxy0=TmGxy;
ux0=0.; uy0=0.;
real errold=1.e10;
real errstop = 0.1/8.*L/nn;
for(int j=0; j<500; ++j){
solve bing([ux,uy,p],[vx,vy,q],solver=UMFPACK)
=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])
-int2d(Th)(Mu*p*Divergence(vx, vy)+Mu*q*Divergence(ux, uy))
+int2d(Th)(-1.e-10*Mu*p*q)
+ 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));
Tauxx=TmGxx+r2*Gdotxx;
Tauyy=TmGyy+r2*Gdotyy;
Tauxy=TmGxy+r2*Gdotxy;
Taunorm=sqrt(.5*(Tauxx^2+Tauyy^2+Tauxy^2));
TmGxx=Tauxx/max(1.,Taunorm);
TmGyy=Tauyy/max(1.,Taunorm);
TmGxy=Tauxy/max(1.,Taunorm);
real normsigest=sqrt(int2d(Th)((TmGxx+2.*Mu/Bn*dx(ux))^2+(TmGyy+2.*Mu/Bn*dy(uy))^2
+(TmGxy+2.*Mu/Bn*(dx(uy)+dy(ux))/sqrt(2.))^2));
real err = sqrt(int2d(Th)(square(ux - ux0)+square(uy - uy0)));
real errsigma = sqrt(int2d(Th)(square(TmGxx-TmGxx0)+square(TmGyy-TmGyy0)+square(TmGxy-TmGxy0)));
cout << err << " " << errsigma << " --iteration error u,sigma /iteration step --" << j<< endl;
if (r2*max(err,errold)<normsigest*errstop) break;
ux0=ux;
uy0=uy;
TmGxx0=TmGxx;TmGyy0=TmGyy;TmGxy0=TmGxy;
errold=err;
}// end of iteration loop on j
normsigma=sqrt((TmGxx+2.*Mu/Bn*dx(ux))^2+(TmGyy+2.*Mu/Bn*dy(uy))^2
+(TmGxy+2.*Mu/Bn*(dx(uy)+dy(ux))/sqrt(2.))^2);
plot(ux,wait=1,fill=1,value=true);
plot([ux,uy],wait=1,nbiso=10,value=true);
normsigmaplot=max(normsigma,1.37);
real[int] valiso(15);
for (int i=0; i<valiso.n;i++)
valiso[i]=1.37+i*0.016/3.;
real[int] colorhsv=[
4./6., 1 , 0.5,
4./6., 1 , 1,
5./6., 1 , 1,
1, 1. , 1,
1, 0.5 , 1
];
plot(normsigmaplot, ps = "normsigma_solution.eps",viso=valiso(0:valiso.n-1),value=1,fill=1,wait =1,hsv=colorhsv);
load "Element_P3"
fespace Rh(Th, P3);
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)
;
real errdiv=int2d(Th)(abs(Divergence(ux, uy)));
cout << "errdiv = " << errdiv << endl;
real errstream=sqrt(int2d(Th)((dx(dyg)-uy)^2+(dy(dyg)+ux)^2));
cout << "errstream = " << errstream << endl;
plot(dyg,ps="ldsre1500.eps",nbiso=10,value=true);
Dear Professor,
Thank you so much for this code. Right now I am trying to run the code but it stops at the 3rd iteration like this
It stops to compile.
How can I fix?
It looks like the iteration is just finished, just type âenterâ in the plot window to see the diagnostics!
For convenience you can add a line
cout << "end of computations, type return in the plot window to plot " << endl;
just after the end of the iteration loop.
Dear Professor,
It looks great. Thank you so much.
Dear Professor,
When I plot the streamlines ,I also want to draw rigid zones as in the examples;gray shaded areas. I couldnât figure out how to include rigid zones into the streamlines plot. Could you give me any ideas please?
I think itâs not possible to draw two different type of plots on the same window in FF++. One thing you could do is post treatment with an image freeware in order to merge two output files.
But a try is to rescale one function to be in a range close to the other. For the version of the code I have posted, the plotting section becomes
normsigmaplot=max(normsigma,1.37);
real[int] valiso(15);
for (int i=0; i<valiso.n;i++)
valiso[i]=1.37-0.0746666+2*i*0.016/3.;// doubled range
real[int] colorhsv=[
4./6., 1 , 0.5,
4./6., 1 , 1,
5./6., 1 , 1,
1, 1. , 1,
1, 0.5 , 1
];
plot(normsigmaplot, ps = "normsigma_solution.eps",viso=valiso(0:valiso.n-1),value=1,fill=1,wait =1,hsv=colorhsv);
load "Element_P3"
fespace Rh(Th, P3);
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)
;
real errdiv=int2d(Th)(abs(Divergence(ux, uy)));
cout << "errdiv = " << errdiv << endl;
real errstream=sqrt(int2d(Th)((dx(dyg)-uy)^2+(dy(dyg)+ux)^2));
cout << "errstream = " << errstream << endl;
plot(dyg,ps="ldsre1500.eps",nbiso=10,value=true,wait=1);
dyg=min(dyg,0.07);// to avoid larger values
dyg=1.37-0.0746666+(dyg+0.005)/0.093*14*0.016/3.;//rescaling
plot(dyg,normsigmaplot,viso=valiso(0:valiso.n-1),value=1,fill=1,wait =1,hsv=colorhsv);
Here the values <1.37 are for the stream function, and the values >1.37 are for \|\sigma\| (values between 1.37 and 1.44 because we focus around the value \sqrt{2}).
1 Like
Dear Professor ,
Thank you so much.
Your result looks correct. You can take for plot commands
plot(p, cmm="Pressure",ps="pressurepoiseuille.eps",wait=1,value=true);
plot([ux, uy], cmm="Velocity",ps="velpoiseuille.eps",value=true);
Please open a new topic for a new subject (your question is not related to the subject âUzawa scheme for stokes Bingham flowâ). To open a new topic press the â+â in the bottom right of the FF window.
1 Like
Thank you for your kind reply sir.