はじめに

今回の記事では障害物があった場合でのポテンシャル線がどのように表示されるかを検証したことについてまとめる．結果はこんな感じ↓(紫：速度)

f:id:takaishi78:20200130123429p:plain — output01

場の初期化

検証するフィールドの広さや、障害物の大きさを決定する．また今回のプログラムではresは離散化した際のブロック数に当たる．（なので格子点の数は+1になるので注意）

obs.size = ofVec2f(150.0, 300.0);
obs.mid = ofVec2f(300.0, 0.0);
obs.left = obs.mid - ofVec2f(obs.size.x/2.0, 0.0);
float eps = 0.000001;
    
rect.size = ofVec2f(800.0, 400.0);
rect.res = ofVec2f(100.0, 50.0);
rect.left = -1.0 * rect.size / 2.0;
rect.delta = rect.size / rect.res;

それぞれの格子点を操作してその点がどのような状況（境界であるのか、それとも障害物内なのかなど）の点であるのかを確認し初期化を行う．

for(int i = 0; i <= rect.res.y; i++){
        for(int j = 0; j <= rect.res.x; j++){
            int id = i * (rect.res.x + 1) + j;
            phi.push_back(0.0);
            psi.push_back(0.0);
            vel.push_back(ofVec2f(0.0, 0.0));
            type.push_back(Type::INSIDE);
            if(i == 0) type[id] = Type::BOTTOM;
            if(i == rect.res.y) type[id] = Type::TOP;
            if(j == 0) type[id] = Type::INLET;
            if(j == rect.res.x) type[id] = Type::OUTLET;
            if(j == nMeshX1 && i < nObsWY) type[id] = Type::OBS_LEFT;
            if(j >= nMeshX1 && j <= nMeshX2 && i == nObsWY) type[id] = Type::OBS_TOP;
            if(j == nMeshX2 && i < nObsWY) type[id]= Type::OBS_RIGHT;
            if(j > nMeshX1 && j < nMeshX2 && i < nObsWY) type[id] = Type::OBSTACLE;
        }
    }
lineNum = 20;

ポテンシャル計算

境界について

http://i.riken.jp/wp-content/uploads/2015/06/secure_4650_080121_toda.pdf

void  Environment::calc()
{
    for(int i = 0; i <= rect.res.y; i++)
    {
        for(int j = 0; j <= rect.res.x; j++)
        {
            int id = j + (rect.res.x + 1) * i;
            //境界面
            if(type[id] == BOTTOM || type[id] == OBS_LEFT || type[id] == OBS_TOP || type[id] == OBS_RIGHT){
                psi[id] = 0.0;
            }else if(type[id] == TOP){
                //psi[id] = rect.size.y;//UW p81
                psi[id] = 1.0;
            }else{
                float delta = 1.0 / rect.res.y;
                //psi[id] = i * rect.delta.y;//INLET, OUTLET, INSIDE
                psi[id] = i * delta;//INLET, OUTLET, INSIDE
            }
            
            if(type[id] == INLET){
                phi[id] = 0.0;
            }else if(phi[id] == OUTLET){
                //phi[id] = rect.size.x;
                phi[id] = 2.0;
            }else{
                float delta = 2.0 / rect.res.x;
                phi[id] = delta * j;
            }
        }
    }

ラプラスの方程式

連続の式

これは流出、流入が等しいことを示す．（ベクトル演算の内積は流出、流入を示す） $$ {\nabla}{\cdot}\vec{v} = 0 \tag{1} $$

$\vec{v}$は(u,v)で構成されている． $$ u = \frac{\partial{\phi}}{\partial{x}}, v = \frac{\partial{\phi}}{\partial{y}} \tag{2} $$

(1)式に(2)式を代入するとラプラス方程式を得る $$ {\nabla}^2{\cdot}{\Phi} = \frac{{\partial}^2{\Phi}} {{\partial}{x}^2} + \frac{{\partial}^2{\Phi}} {{\partial}{y}^2} = 0 \tag{3} $$

差分法

コンピュータ上で数値解析を行う場合、微分演算を差分に変換して解く必要がある．

2階微分の近似

$$ \frac {d^{2}f} {dx^{2}} = \frac{ \frac{df}{dx}_{i+1} - \frac{df}{dx}{_i}}{{\Delta}x} $$

分子の前半は前進差分、後半は後退差分を用いて表すと

$$ \frac{ \frac{df}{dx}_{i+1} - \frac{df}{dx}_i}{{\Delta}x} = \frac{\frac{f_{i+1} - f_i} {{\Delta}x} - \frac{f_i - f_{i-1}} {{\Delta}x}}{{\Delta}x} = \frac{f_{i+1} - 2f_i + f_{i-1}} {{\Delta}x^2} $$

これで差分法をもちいた二階微分の近似方法がわかった．

ラプラス方程式(差分法)

流れ関数に対応するラプラス方程式は

$$ {\nabla^2}{\Psi} = 0 \\ \frac{ {\partial^2}{\Psi} }{ {\partial} {x^2}} + \frac{ {\partial^2}{\Psi} }{ {\partial} {y^2}} = 0 \tag{4} $$

2階微分の部分に対して差分法を適応していくと(4)式は以下のように書き換えることが出来る $$ \frac{{\Psi}_{i+1, j} - 2{\Psi}_{i, j} + {\Psi}_{i-1, j}} {{\Delta}x^2} + \frac{{\Psi}_{i, j+1} - 2{\Psi}_{i, j} + {\Psi}_{i, j-1}} {{\Delta}y^2} = 0\tag{5} $$

${\Psi_{i,j}}$を左辺に移動させると

$$ {\Psi_{i,j}} = \frac{ ({\Delta}y)^2 ({\Psi_{i-1,j}} + {\Psi_{i+1,j}} ) + ({\Delta}x)^2 ({\Psi_{i,j-1}} + {\Psi_{i,j+1}} ) } { 2({ ({\Delta}x})^2 { ({\Delta}y})^2) } $$

サンプルコード

ラプラス方程式を具体的に解いたサンプル．注目している点の周囲4点を用いてポテンシャルの値を算出している．

    //FDM
    int iteration = 5000;
    float tolerance = 0.000001;
    int cnt = 0;
    float error = 0.0;
    float maxError = 0.0;
    float dx2 = rect.delta.x * rect.delta.x;
    float dy2 = rect.delta.y * rect.delta.y;
    
    float pp;
    
    while (cnt < iteration)
    {
        maxError = 0.0;
        for (int i = 1; i < rect.res.y; i++)
        {
            for (int j = 1; j < rect.res.x; j++)
            {
                int id = j + i * (rect.res.x+1);
                if(type[id] != INSIDE) continue;
                //周囲４つ
                int id1 = j     + (i-1) * (rect.res.x + 1);//y-1
                int id2 = j     + (i+1) * (rect.res.x + 1);//y+1
                int id3 = (j-1) + i     * (rect.res.x + 1);//x-1
                int id4 = (j+1) + i     * (rect.res.x + 1);//x+1
                
                //laplatian
                //psi
                pp = dy2 * (psi[id3] + psi[id4]) + dx2 *(psi[id1] + psi[id2]);
                pp /= 2.0 * (dx2 + dy2);

反復法(ガウスザイデル法)

反復法では許容誤差以下になるまで計算を繰り返し行う．

$$ |f{_i}{^n} - f{_i}{^{n-1}}| {\leq} {\epsilon} (i=1,2,N-1) $$ 流体が流れる範囲の全てのポテンシャルを計算し、前のタームとの誤差がその中でも最大となる値をmaxerroeとしている．この値が誤差の許容範囲以下に収まるまで計算を繰り返す．(任意のiで許容誤差内に収まらないといけない)

//さっきのコードの続き
//１つ前のψと比較
                error = abs(pp - psi[id]);
                if (error > maxError) maxError = error;
                psi[id] = pp;
            }
        }
        if (maxError < tolerance) break;
        cnt++;
    }