Hi everyone,
I write a program to simulate the cavity flow’s linear temporal growth rates and the neutral curve, but the results aren’t ideal,
I put the program of linear temporal growth rates as flowing, looking forward your help! Thanks very much!
load “msh3”
load “tetgen”
load “medit”
int N=70;
real [int,int] coordinate((N+1)(N+1),2);
for(int i=0;i<=N;i++)
{
for(int j=0;j<=N;j++)
{
coordinate(i(N+1)+j,0)=(cos(jpi/N)+1)/2;
coordinate(i(N+1)+j,1)=(cos(ipi/N)+1)/2;
cout<<coordinate(i(N+1)+j,0)<<" "<<coordinate(i*(N+1)+j,1)<<endl;
}
}
mesh Th2d=triangulate(coordinate(:,0),coordinate(:,1));
int[int] r1=[1,10],r2=[1,20];
mesh Th=change(Th2d,label=r1);
savemesh(Th,“Th2d.mesh”);
fespace Xh(Th,P2); //for velocity
fespace Mh(Th,P1); //for pressure and stress
fespace XhXhXhMh(Th,[P2,P2,P2,P1]);
fespace XhXhXh(Th,[P2,P2,P2]);
Xh Ux,Uy,Vx,Vy,dUx,dUy;
Mh P,dP,Q;
Ux=(y==1);
Uy=(x==0); Uy=-Uy; //The boundary conditions
Xh ux,uy,uz,ux0,uy0,uz0;
Xh uxold,uyold,vx,vy,vz;
Mh p,p0;
Mh dp,pp,pold,q;
// Physical parameter
real re=500, rebegin=0, reend=1000, Ha=0, k=0.9;
complex shift=0.+0i, omega;
real nu= 1./re,errold,abeigen,abeigenold;
real eps=1e-6;
real err=0;
complex[int] Bu1(XhXhXhMh.ndof),u(XhXhXhMh.ndof);
complex[int] Au1(XhXhXhMh.ndof),Au2(XhXhXhMh.ndof);
varf op([ux,uy,uz,p],[vx,vy,vz,q])=
int2d(Th)(-vx*(uxolddx(ux)+uxdx(uxold)+uyolddy(ux)+uydy(uxold))
- vy*(uxolddx(uy)+uxdx(uyold)+uyolddy(uy)+uydy(uyold))
- vz*(uxolddx(uz)+uyolddy(uz))
+ p*(dx(vx)+dy(vy)-1ikvz)
- nu*(dx(ux)dx(vx)+dy(ux)dy(vx)+dx(uy)dx(vy)+dy(uy)dy(vy)+dx(uz)dx(vz)+dy(uz)dy(vz))
- nukk(uxvx+uyvy+uzvz)
+ q*(dx(ux)+dy(uy))-1ikuzq
+ pq1e-18
)
+ int2d(Th)(1ishift*(uxvx+uyvy+uz*vz))
+ on(10,ux=0,uy=0,uz=0);
varf a([ux,uy,uz,p],[vx,vy,vz,q])=
int2d(Th)(-vx*(uxolddx(ux)+uxdx(uxold)+uyolddy(ux)+uydy(uxold))
- vy*(uxolddx(uy)+uxdx(uyold)+uyolddy(uy)+uydy(uyold))
- vz*(uxolddx(uz)+uyolddy(uz))
+ p*(dx(vx)+dy(vy)-1ikvz)
- nu*(dx(ux)dx(vx)+dy(ux)dy(vx)+dx(uy)dx(vy)+dy(uy)dy(vy)+dx(uz)dx(vz)+dy(uz)dy(vz))
- nukk(uxvx+uyvy+uzvz)
+ q*(dx(ux)+dy(uy))-1ikuzq
+ pq*1e-18
)
+ on(10,ux=0,uy=0,uz=0);
varf b([ux,uy,uz,p],[vx,vy,vz,q])=
int2d(Th)(1i*(uxvx+uyvy+uz*vz))
;
//find the first piont
for(int t=0;t< 20;t++)
{
solve Baseflow([dUx,dUy,dP],[Vx,Vy,Q])=
int2d(Th)( Vx*(dUxdx(Ux)+dUydy(Ux))+Vy*(dUxdx(Uy)+dUydy(Uy))
+Vx*(Uxdx(dUx)+Uydy(dUx))+Vy*(Uxdx(dUy)+Uydy(dUy))
+nu*(dx(dUx)dx(Vx)+dy(dUx)dy(Vx))
+nu(dx(dUy)dx(Vy)+dy(dUy)dy(Vy))
-dP(dx(Vx)+dy(Vy))
-Q(dx(dUx)+dy(dUy))
+1e-8dPQ
)
-int2d(Th)(Vx(Uxdx(Ux)+Uydy(Ux))+Vy*(Uxdx(Uy)+Uydy(Uy))
+nu*(dx(Ux)*dx(Vx)+dy(Ux)dy(Vx))
+nu(dx(Uy)dx(Vy)+dy(Uy)dy(Vy))
-P(dx(Vx)+dy(Vy))
-Q(dx(Ux)+dy(Uy))
)
+on(10,dUx=0,dUy=0)
;
cout<<“t=”<<t<<" dU.max “<<dUx[].max<<” “<<dUy[].max<<” "<<"dP.max= “<<dP[].max<<endl;
err= dUx[].linfty + dUy[].linfty + dP[].linfty;
Ux[]-=dUx[];
Uy[]-=dUy[];
P[]-=dP[];
cout << t << " err = " << err << " " << eps << " rey =” << re << endl;
if(err < eps) break; // converge
if( t>10 && err > 10.) break; // Blowup ???
}
if(err < eps)
{
uxold=Ux;
uyold=Uy;
//stability
int kkmax=0;
real erreigen;
matrix<complex> OP= op(XhXhXhMh,XhXhXhMh,solver=sparsesolver);
matrix<complex> A= a(XhXhXhMh,XhXhXhMh,solver=sparsesolver);
matrix<complex> B= b(XhXhXhMh,XhXhXhMh,solver=CG);
int nev=40;
complex[int] ev(nev); //to store nev eigen value
XhXhXhMh<complex>[int] [eux,euy,euz,ep](nev); //to store nev eigen vector
int ncv=4*nev;
int kk=EigenValue(OP,B,sigma=shift,value=ev,vector=eux, tol=1e-20,maxit=0,ncv=ncv);
kk=min(kk,nev); // some time the number of converged eigen value
// can be greater than nev;
abeigenold=-1;
for (int i=0;i<kk;i++)
{
//kkmax=(imag(ev[i])>imag(ev[i+1])?i:(i+1));
//kkmaxold=(imag(ev[kkmax])>imag(ev[kkmaxold])?kkmax:kkmaxold);
ux=eux[i];
uy=euy[i];
uz=euz[i];
p=ep[i];
u=[ux[],uy[],uz[],p[]];
complex v=ev[i];
Bu1=B*u;
Bu1=v*Bu1;
Au1=A*u;
u=Au1-Bu1;
erreigen=norm(u.sum);
abeigen=imag(ev[i]);
cout<<"re= "<<re<<" i= "<<i<<" ev[i]= "<<ev[i]<<" erreigen= "<<erreigen<<endl;
if(abeigen>abeigenold)
{
abeigenold=abeigen;kkmax=i;
cout <<" kkmax= "<<kkmax <<endl;
}
}
omega=ev[kkmax];
ux=eux[kkmax];
uy=euy[kkmax];
uz=euz[kkmax];
p=ep[kkmax];
cout<<"re= "<<re<<" omega= "<<omega<<endl;
//shift=omega;
if(imag(omega)>0&&abs(imag(omega))>1e-5)
{
k=k-0.05;
//shift=omega; cout<<"restart re"<<endl;
}
else if(imag(omega)<0&&abs(imag(omega))>1e-5)
{
k=k+0.05;
//shift=omega; cout<<"restart re"<<endl;
}
else
{cout<<"sucesse to find the first point"<<endl;
cout<<"k="<<k<<" re="<<re<<endl;}
}
else
{ // if blowup, increase nu (more simple)
cout << " restart re = " << re <<" converge wrong"<<endl;
}
real alp,alpbegin=0,alpend=5,err3;
int itmax=150, iter;
k=k+0.1;
while (k<alpend)
{
real erreigen;int kkmax=0;
shift=0+0i;
cout<<“shift=”<<shift<<endl;
matrix OP= op(XhXhXhMh,XhXhXhMh,solver=sparsesolver);
matrix A= a(XhXhXhMh,XhXhXhMh,solver=sparsesolver);
matrix B= b(XhXhXhMh,XhXhXhMh,solver=CG);
int nev=4;
complex[int] ev(nev); //to store nev eigen value
XhXhXhMh<complex>[int] [eux,euy,euz,ep](nev); //to store nev eigen vector
int ncv=4*nev;
int kk=EigenValue(OP,B,sigma=shift,value=ev,vector=eux, tol=1e-20,maxit=0,ncv=ncv);
kk=min(kk,nev); // some time the number of converged eigen value
// can be greater than nev;
abeigenold=-1;
for (int i=0;i<kk;i++)
{
//kkmax=(imag(ev[i])>imag(ev[i+1])?i:(i+1));
//kkmaxold=(imag(ev[kkmax])>imag(ev[kkmaxold])?kkmax:kkmaxold);
ux=eux[i];
uy=euy[i];
uz=euz[i];
p=ep[i];
u=[ux[],uy[],uz[],p[]];
complex v=ev[i];
Bu1=B*u;
Bu1=v*Bu1;
Au1=A*u;
u=Au1-Bu1;
erreigen=norm(u.sum);
abeigen=imag(ev[i]);
cout<<"re= "<<re<<" i= "<<i<<" ev[i]= "<<ev[i]<<" erreigen= "<<erreigen<<endl;
if((abeigen>abeigenold)&&(abs(real(v))<1e-5))
{
abeigenold=abeigen;kkmax=i;
cout <<" kkmax= "<<kkmax <<endl;
}
}
omega=ev[kkmax];
ux=eux[kkmax];
uy=euy[kkmax];
uz=euz[kkmax];
p=ep[kkmax];
cout<<"k= "<<k<<" re= "<<re<<" omega= "<<omega<<endl;
k=k+0.1;
}