convect2.edp

来自「FreeFem++可以生成高质量的有限元网格。可以用于流体力学」· EDP 代码 · 共 36 行

// This a the rotating hill problem with one turn.
// First 1/2 turn is a convection equation and second 1/2 a convection diffusion

border a(t=0, 2*pi)     {    x = cos(t);    y = sin(t);  }; // the unit circle
mesh th2 = buildmesh(a(35));                                 // triangulates the disk
mesh th= trunc(th2,1,split=2);
fespace Vh(th,P1);
fespace Vh2(th2,P2);
Vh2 v = exp(-10*((x-0.3)^2 +(y-0.3)^2));                  // initial condition
plot(v);

real dt = 0.17,t=0;                                                 // time step
Vh u1 = y, u2 = -x;                                        // rotation velocity
int i;
Vh2 vv,vo; // work  Finite element function 
for ( i=0; i< 20 ; i++) {
    t += dt;
    vo=v;
    v=convect([u1,u2],-dt,vo);                          // convect v by u1,u2, dt seconds, results in f
 // convec(u1,u2,dt,v,v) won't work
    plot(v,cmm="convection: t="+t + ", min=" + v[].min + ", max=" +  v[].max,wait=0);
};

problem  A(v,vv,solver=CG) = int2d(th2)(v*vv/dt + 0.01*(dx(v)*dx(vv)+dy(v)*dy(vv)) )
  + int2d(th2)(-vv*convect([u1,u2],-dt,vo)/dt)  + on(a,v=0);
  
for ( i=0; i< 20 ; i++) 
{ 
    t += dt;
  vo=v;
  A; // solve le problem A
  plot(v,cmm="convection& diffusive: t="+t + ", min=" + v[].min + ", max=" +  v[].max);
};  

plot(v,wait=1);