load "MUMPS_seq" int m=20;// 1/h int withadapmesh=1; real g=-9.81, nuf=1./100., nus=10000, rhof=0.1, rhos = 1, dt=0.025, T= 100; real L=6,XM=1, r1=0.3, r2=0.05, xc=0.4, yc=L-0.3, omega=10*pi/180; real Ladp = 0.5, dhhh = 0.1/Ladp, hmax= 1./4; int lrec =1, lellipse=2; // labels of B.C. // This is a rectangle (0,1)x(0,L) border a1(t=0,XM){x=t;y=0;label=lrec;} border a2(t=0,L){x=XM; y=t;label=lrec;} border a3(t=XM,0){x=t; y=L;label=lrec;} border a4(t=L,0){x=0;y=t;label=lrec;} // This is an ellipse centered at (xc,yc), of radii r1,r2 and rotated by omega border disk(t=0,2*pi){ x=xc+r1*cos(t)*cos(omega)+r2*sin(t)*sin(omega); y=yc-r1*cos(t)*sin(omega)+r2*sin(t)*cos(omega); label=lellipse; } //plot( a1(m) + a2(3*m) +a3(m)+a4(3*m)+disk(2*m), wait=1); int ndisk= 8*(r1+r2)*m; func Boundary = a1(m*XM) + a2(m*L) +a3(m*XM)+a4(m*L)+disk( ndisk); func bool bbdisk0(real omega,real[int] bb) { real[int] xx(ndisk),yy(ndisk); real dt = 2*pi/ndisk,t=0; for( int i=0; i< ndisk; ++i) { xx[i]=r1*cos(t)*cos(omega)+r2*sin(t)*sin(omega); yy[i]=r1*cos(t)*sin(omega)+r2*sin(t)*cos(omega); t += dt; } bb[0]=xx.min; bb[1]=xx.max; bb[2]=yy.min; bb[3]=yy.max; return 1; } func real Xmaxe(){ int n=8*(r1+r2)*m; real[int] xx(n); for ( int i=0; i< n; ++i) { real t = i*2*pi/n; xx[i] = xc+r1*cos(t)*cos(omega)+r2*sin(t)*sin(omega); } return xx.max; } func Xr = (x-xc)*cos(-omega)+(y-yc)*sin(-omega); func Yr = -(x-xc)*sin(-omega)+(y-yc)*cos(-omega); real hmin = 0.03; real epss = hmin/4; real dhmin = 0.005; real[int] bb(4); //func dellipse2 = sqr(Xr/r1)+sqr(Yr/r2)-1; func hellipse = min( 1./m + (max(abs(y-yc),Ladp)-Ladp)*dhhh , hmax); //min(hmin + abs(dellipse2)*dhmin,0.3); mesh Th=buildmesh( Boundary); Th=adaptmesh(Th,hellipse,IsMetric=1); plot(Th,wait=1); fespace Wh(Th,[P1b,P1b,P1]); fespace Qh(Th,P0); Qh nu,rho, phi; Wh [u1,u2,p],[v1,v2,q], [u1old,u2old,pold]=[0,0,0]; macro div(u1,u2) (dx(u1)+dy(u2))// macro D(u1,u2) [ dx(u1), (dx(u2)+dy(u1))*0.5, (dx(u2)+dy(u1))*0.5, dy(u2)]// // Navier-Stokes equations with gravity force, integrated in the rectangle, // with varying density rho and varying viscosity but div(u)=0 in each triangle problem NStokes ([u1,u2,p],[v1,v2,q]) = int2d(Th)( rho*(u1*v1+u2*v2)/dt + 2.*nu*(D(u1,u2)'*D(v1,v2) ) - p*q*1e-8 - p*div(v1,v2) - q*div(u1,u2) ) -int2d(Th)( + rho*g*v2 + rho*convect([u1old,u2old],-dt,u1old)*v1/dt + rho*convect([u1old,u2old],-dt,u2old)*v2/dt ) + on(lrec,u1=0,u2=0); //real J=rhos*pi*r1*r2*(r1^2+r2^2)/4; verbosity=0; for(real t=0;t 0 et bb[1]+xc < 1 xc = max(xc,-bb[0]+epss); xc = min(xc,XM-bb[1]-epss); // real xmc =Xmaxe()*1.1; cout << " xmc = " << xc << " " << bb[0] << " " << bb[1] << endl; assert(xc+bb[0] > 0); assert(xc+bb[1] < XM); // xc = max(0., min(xc,XM- xmc), xmc); if ( yc < 1*(abs(sin(omega))*r1 +abs(cos(omega)*r2))) break; cout << " mov : "<< du1 << " "<< du2 << " " << domega << " / " << xc << " " << yc << " " << omega*180/pi << endl; plot(coef=0.1,cmm=" [u1,u2] et p ",p,[u1,u2], value=true,wait=false); //plot(a1(m) + a2(3*m) +a3(m)+a4(3*m)+disk(2*m),WindowIndex=1); // Mesh is rebuilt with new ellipse position //plot(Boundary); Th=buildmesh( Boundary ) ; Th=adaptmesh(Th,hellipse, IsMetric=1); [u1old,u2old,pold]=[u1,u2,p]; } plot(coef=0.1,cmm=" [u1,u2] et p ",p,[u1,u2], value=true,wait=false);