// resolution de NS avec THERMIQUE dans une cavite differentielle // PYL oct 2005 / fev 06 / dec 17 // inspire de l'exemple NSP1P2.edp //file //http://www.lmm.jussieu.fr/~lagree/COURS/ENSTA/web/CLib/thermiKVIT/srcrun/thermiquekimarcheKviT.edp //freefem++ thermiquekimarcheKviT.edp //http://www.lmm.jussieu.fr/~lagree/COURS/ENSTA/web/CLib/thermiKVIT/index.html // voir le cours de convection naturelle // exec("echo KviT kimarche"); real s0=clock(); real h0=1; real L0=1; int n=20; // longueur L0, epaisseur h0 // definition des cotes Maillage du carré border b(t=0,1) { x= t*L0; y = 0 ; }; border d(t=0,1) { x= L0; y = h0 * t ; } ; border h(t=1,0) { x= L0*(t); y = h0 ;}; border g(t=0,1) { x= 0; y = h0 * (1-t); }; mesh Th= buildmesh(b(n)+d(n)+h(n)+g(n)); plot(Th,wait=0,ps="cuve.eps"); fespace Vh2(Th,P2); fespace Vh(Th,P1); Vh2 u2,v2,up1,up2; Vh2 u1,v1; Vh2 T,TT,Told; exec("pwd"); real reynolds=500.; u1=0 ; u2=0 ; up1=u1; up2=u2; Told=0; T=0; Vh p=0,q; real J1=0,J2=-1; real dt = 0.1; real alpha=1/dt; real nu=1./reynolds; int i=0,iter=0; //{ ifstream f("JJ1.txt"); f>>J1;} //{ ifstream f("JJ2.txt"); f>>J2;} // dans JJ1.txt il y a 0.00 // dans JJ2.txt il y a -1.00 //{ ifstream f("reynolds.txt"); f>>reynolds;} //dans reynolds.txt il y a 1000 problem NS ([u1,u2,p],[v1,v2,q]) = int2d(Th)( alpha*( u1*v1 + u2*v2) + nu * ( dx(u1)*dx(v1) + dy(u1)*dy(v1) + dx(u2)*dx(v2) + dy(u2)*dy(v2) ) - p*q*(0.000001) - p*dx(v1) - p*dy(v2) - dx(u1)*q - dy(u2)*q ) + int2d(Th)(J2*Told*v2) + int2d(Th)( -alpha*convect([up1,up2],-dt,up1)*v1 -alpha*convect([up1,up2],-dt,up2)*v2 ) + on(b,u1=0,u2=0) + on(h,u1=0,u2=0) + on(g,u1=0,u2=0) + on(d,u1=0,u2=0) ; problem NST (T,TT,solver=Crout,init=i) = int2d(Th)( T*TT/dt + nu * ( dx(T)*dx(TT) + dy(T)*dy(TT)) ) + int2d(Th) ( -(1./dt)*convect([u1,u2],-dt,Told)*TT ) + on(d,T=1) + on(g,T=-1) //+ on(h,T=0) //+ on(b,T=1) ; Vh2 psi,phi; problem streamlines(psi,phi) = int2d(Th)( dx(psi)*dx(phi) + dy(psi)*dy(phi)) + int2d(Th)( -phi*(dy(u1)-dx(u2))) + on(b,psi=0) + on(g,psi=0) + on(h,psi=0) + on(g,psi=0); //------------------------------------------------------- int nbiter = 100; real coefdt = 0.25^(1./nbiter); real coefcut = 0.25^(1./nbiter) , cut=0.01; real tol=0.5,coeftol = 0.25^(1./nbiter); real t=0; for (iter=0;iter<=nbiter;iter++) { real errl2 = sqrt(int2d(Th)( (u1-up1)^2 + (u2-up2)^2 +(T -Told)^2) ); cout << "t = "<< t <<" ---------------- "<< iter <<" +++++++++++++++ errl2 " << errl2<< endl; alpha=1./dt; nu=1./reynolds; for (i=0;i<=10;i++) { up1=u1; up2=u2; NS; NST; Told=T; t=t+dt; if ( !(i % 4)) { plot(coef=.2,cmm=" [u1,u2] et T ",T,[u1,u2],ps="u1u2T.eps"); plot( T, ps="ture.eps"); real x,y; n=100; { ofstream gnu("NplotT.gp"); gnu << "# vals" << endl; for (int i=0;i<=n;i++) { x =(i*L0)/n*1.0; y=h0/2; gnu<< x <<" "<< T(x,y) << " "<< T(y,x) << endl; } } { ofstream gnu("NplotU.gp"); gnu << "# vals" << endl; for (int i=0;i<=n;i++) { x =(i*L0)/n*1.0; y=h0/2; gnu<< x <<" "<< u1(x,y)<<" "<< u2(x,y)<< endl; } } } cout << "CPU " << clock()-s0 << "s "<< endl; cout << "~~~~~~~~~~~~~~~~~~ J2=" << J2 << " R=" << reynolds << " " <>J1;} // ifstream f("JJ2.txt"); f>>J2;} //{ ifstream f("reynolds.txt"); f>>reynolds;} // dans reynolds.txt il y a 1000 // sauvegarde gnuplot if ((t>15) &&((errl2 ) < (5e-5))) break; //if (iter>= nbiter) break; } streamlines; plot(wait=0,psi,ps="psi.eps"); cout << "t = "<< t << endl; cout << "CPU " << clock()-s0 << "s " << endl; // // // //