potential.edp

FreeFem++可以生成高质量的有限元网格。可以用于流体力学

// file potential.edp

real S=99;
border C(t=0,2*pi) {  x=3*cos(t);  y=3*sin(t);}
border Splus(t=0,1){  x = t; y = 0.17735*sqrt(t)-0.075597*t
        - 0.212836*(t^2)+0.17363*(t^3)-0.06254*(t^4); label=S;}
border Sminus(t=1,0){  x =t; y= -(0.17735*sqrt(t)-0.075597*t
        -0.212836*(t^2)+0.17363*(t^3)-0.06254*(t^4)); label=S;}
mesh Th= buildmesh(C(50)+Splus(70)+Sminus(70));
fespace Vh(Th,P2); Vh psi,w;

solve potential(psi,w)=int2d(Th)(dx(psi)*dx(w)+dy(psi)*dy(w))+
  on(C,psi = y)+ on(S,psi=0);

plot(psi,wait=1);
border D(t=0,2){x=1+t;y=0;}
mesh Sh = buildmesh(C(25)+Splus(-90)+Sminus(-90)+D(200));
fespace Wh(Sh,P1); Wh v,vv;
int steel=Sh(0.5,0).region, air=Sh(-1,0).region;
fespace W0(Sh,P0);
W0 k=0.01*(region==air)+0.1*(region==steel);
W0 u1=dy(psi)*(region==air), u2=-dx(psi)*(region==air);
Wh vold = 120*(region==steel);
real dt=0.05, nbT=50;
int i;
problem thermic(v,vv,init=i,solver=LU)= int2d(Sh)(v*vv/dt + k*(dx(v) * dx(vv) + dy(v) * dy(vv))
                + 10*(u1*dx(v)+u2*dy(v))*vv)- int2d(Sh)(vold*vv/dt);
for(i=0;i<nbT;i++){
    v=vold; thermic;
 plot(v);
}
plot(v,wait=1);