// potential2d-v2.edp // http://nalab.mind.meiji.ac.jp/~mk/ana/potential2d-v2.edp // 2次元非圧縮ポテンシャル流 // 速度ポテンシャル,速度を求め、等ポテンシャル線, 速度場を描く border Gamma1(t=0,pi/6) { x = cos(t); y = sin(t); } border Gamma2(t=pi/6,pi) { x = cos(t); y = sin(t); } border Gamma3(t=pi,4*pi/3) { x = cos(t); y = sin(t); } border Gamma4(t=4*pi/3,2*pi) { x = cos(t); y = sin(t); } int m=20; mesh Th=buildmesh(Gamma1(m)+Gamma2(m)+Gamma3(m)+Gamma4(m)); plot(Th, wait=1, ps="Th.eps"); // 次の2行は区分1次多項式を使うという意味 fespace Vh(Th,P1); Vh phi, psi, v, v1, v2; // 境界条件の設定 func Vn1=-2.0; // Gamma1 では勢いよく注水 func Vn3=1.0; // Gamma3 からゆるく排水 func Vn=(x>0)*(-2.0)+(x<0)*1.0; // xの符号でGamma1, Gamma3 を区別してまとめる // 速度ポテンシャルφを求める solve Laplace(phi,v) = int2d(Th)(dx(phi)*dx(v)+dy(phi)*dy(v)) -int1d(Th,Gamma1)(Vn1*v)-int1d(Th,Gamma3)(Vn3*v); // -int1d(Th,Gamma1,Gamma3)(Vn*v); // 平均を0にする (解の一意性がないため、時々生じる大きなズレを消す real mean=int2d(Th)(phi)/int2d(Th)(1); phi=phi-mean; // φの等高線 (等ポテンシャル線) を描く plot(phi,ps="contourpotential.eps",wait=1); // ベクトル場 (v1,v2)=∇φ を描く (ちょっと雑なやり方) v1=dx(phi); v2=dy(phi); plot([v1,v2],ps="vectorfield.eps",wait=1); // 等ポテンシャル線とベクトル場を同時に描く plot([v1,v2],phi,ps="both.eps", wait=1); // 流れ関数を求める func angle=(atan2(y,x)>=0)*atan2(y,x)+(atan2(y,x)<0)*(atan2(y,x)+2*pi); solve LaplaceD(psi,v) = int2d(Th)(dx(psi)*dx(v)+dy(psi)*dy(v)) +on(Gamma1,psi=-2*angle(x,y)) +on(Gamma2,psi=-pi/3) +on(Gamma3,psi=angle(x,y)-4*pi/3) +on(Gamma4,psi=0); plot(psi,ps="streamline.eps",wait=1); plot([v1,v2],phi,psi,ps="all.eps", wait=1);