CFD 3D：シンプルな水シミュレーター

はじめに

CFD（計算流体力学） -計算流体力学。

液体のさまざまなプロセスと、さまざまな種類の液体（たとえば、蜂蜜、オイル-これらはすべて液体）のシミュレーションに使用されます。



この投稿では、表面が開いて障害物がある通常の水の2Dシミュレーターを検討します（3Dバージョンでは、すべてが同様であり、 ソースコードが利用可能です）。

水面は水と空気を分離する境界であり、これにより波、水滴などをシミュレートできます。



いくつかの理由でシミュレーターを作成することにしました。 それらの1つは、インターネット上の主題に関する適切に記述されたソースの欠如です。 私が見つけたのは、Fortranまたは2dのみでした。



そのため、シミュレータの主な機能：

• リアルタイムモード
• シンプルで明確な計算スキーム+明示的、暗黙的
• オープンソース
• 2Dと3Dがあります
• 組み込みのシンプルなOpenGLレンダー
• xmlシーンの説明

計算の観点から、私はNavier-Stokesの方程式（およびインターネットで見つかったソース）から始めました。

これは、粘性のある非圧縮性流体の運動を記述する偏微分方程式のシステムです。 定義上、すべての液体は非圧縮性であるため、ここで非圧縮性という言葉はダニに似ています。 粘度は、個々の粒子間の摩擦を示し、割り当てられたボリューム全体を埋めるときに、ガスの場合のように、粒子が異なる方向に散乱するのを防ぎます。 液体粒子は、できるだけ一緒にくっつきます。



私のシミュレーターでは、わずかに修正された方程式のセットを使用しましたが、以下でさらに詳しく説明します...

流体力学方程式

したがって、一般的な（ベクトル）形式の基本的なナビエストークス方程式は次のとおりです。



ここに：

• uは、速度フィールドの従来の指定です（2Dでは、速度はx方向とy方向の2つの成分で構成され、それぞれuとvで示されます）
  
  
• pは圧力です
• tは時間です
• rhoシンボル-流体密度
• 記号nu-動粘度
• Fext-重力（g）など、流体に作用する外力

最初の方程式は運動方程式、2番目は連続方程式です。

運動方程式は、右にあるニュートン方程式ma = Fに似ています-流体に作用する力の合計-これは圧力、拡散、重力です...



2Dバリアントの場合、方程式は次の形式を取ります（ここでは、これらは既に無次元形式で提示されています）。







運動方程式は、2つの速度成分-uおよびvの2つの方程式で表されます。

新しい係数に関して、Reはレイノルズ数（無次元変数に切り替えたときに判明）であり、2つの前の係数（係数）をすぐに置き換えます。 密度と粘度。 それは水がどのように見えるかによって異なります-蜂蜜や普通の水のように Reが多ければ多いほど、液体は通常の水に似ています。

括弧内は、液体の粘性挙動（または拡散）を定義する表現です。 右側の残りの用語（重力を除く）は対流です。 左側-速度と圧力。



次に、これらの方程式を少し単純化しました。







ご覧のように、対流は消滅し、液体の挙動がわずかに悪化しますが、全体的には視覚的には非常に正常に見えます。さらに、対流を適切に離散化するために、特別なスキームが使用されます-これはかなり面倒に見えます。



方程式を解くとき、Splittingと呼ばれる特別なスキームが使用されます。これについて詳しく説明したいと思いますが、これは別の投稿の可能性が高くなります。

Bridson R. Fluid Simulation for Computer Graphics.pdfの「数値シミュレーションの概要」。

そして、流体力学の数値シミュレーションGriebel M Dornseifer T Neunhoeffer T SIAM 1998



分割方程式の解法は次のとおりです。







時間ステップでは、オイラー法が使用されます。最初のステップでも、圧力を下げて、速度FとGの中間計算を取得します。これらの速度は、圧力を考慮しないことに加えて、連続方程式も満たしません。 トリッキーな操作により、圧力を計算する方程式が見つかります。これにより、連続方程式の問題も解決されます。 すでに計算されたFとGからこの方程式を解くと、新しい圧力Pが見つかります。

次に、3番目のステップで、計算された圧力を考慮して速度が修正されます。

これにより、1つの完全なタイムステップが完了し、次のステップですべてが新たに繰り返されます。

方程式の離散化

方程式を離散化するために、有限差分法が使用されます。 時間ステップでは、標準のオイラー法が使用されます-概略的には次のようになります：u（n + 1）= u（n）+ dt * f、

ここで、fは、時間差を含まない方程式の全体として理解されます。



一般的に、新しい瞬間（u（n + 1））の速度値は、前の瞬間（u（n））の値から求められます。

この投稿では方程式の離散的な類似体については説明しませんが、これは流体力学の数値シミュレーション（Griebel M Dornseifer T Neunhoeffer T SIAM 1998）で非常によくできています。



離散方程式を解くために使用されるのは、解のマトリックス法ではなく、反復ガウスザイデル法です。 行列の予備的な組み立てを必要とせず、一般に中間配列を必要としません。計算スキームを簡単に変更できます。おおよその解法は1回の反復後に既にあり、シミュレーション全体を大幅に加速します。



この投稿では、2Dケースを検討します。主な重点は境界条件の説明です。 それらは方程式を解くのに最大の困難を引き起こします。 シミュレートされた領域全体が、水平方向のimaxポイントと垂直方向のjmaxポイントに分割されます。 imax * jmaxセルのグリッドになります。

4つの境界線から境界点が追加されます。 合計で、サイズ（imax + 2）*（jmax + 2）の配列が取得されます。



各セルには、速度と圧力の独自の値があります-計算グリッド上のベクトル速度場とスカラー圧力場について話すのが慣習です。

Uはxのセル内の粒子速度です

Vは、yのセル内の粒子速度です。

P-圧力



通常、計算された変数（U、V、P）をセルの中心に配置するのが一般的ですが、流体モデリングの場合、これは常にソリューションに問題を引き起こします-完全に正しくなく振動していることが判明しました。 そのため、CFDはスタッガードグリッドを使用します-リープフロッグとも呼ばれます。







図から、速度はセル自体ではなく、セルの面、u-セルの右境界線、v-上部境界線にあることがわかります。

壁の境界条件

差分法で計算するには、計算グリッドの境界に値（u、v、p）を設定する必要があります。これらは左から下、右上の4つの壁です。境界条件は、滑り止めタイプとフリースリープタイプです。 他のタイプもありますが、これらは一般的な計画ではなく、より専門的な条件です。



フリースリープ -これは、液体が壁に沿って自由に滑ることを意味します。まるで摩擦がなく、液体がそれに沿って動くのを妨げるものがないかのように。

この投稿では、このタイプの境界条件は考慮しません。



滑り止めは粘着状態です。つまり、液体が壁にぶつかると速度が低下します。

これは、流体速度が壁速度と一致することを意味します（つまり、この例では、固定境界でゼロに等しい）。



たとえば、右側の境界のみを考えます。速度= 0のuコンポーネント。 それは壁に垂直であり、水は境界を貫通してはいけません。

滑りのない境界の場合のvコンポーネントも0ですが、スタッガードメッシュの場合は修正する必要があります。vコンポーネントの場合、vが直接境界にないため、式をわずかに修正する必要があります。

壁のvは、最後の2つのセル間の平均に等しくなります。 v_g =（v [imax + 1] [j] + v [imax] [j]）/ 2

v_gはゼロ（（v [imax + 1] [j] + v [imax] [j]）/ 2 = 0）であり、ここからプログラムで設定する必要がある値v [imax + 1] [j]を見つけます。

v [imax + 1] [j] =-v [imax] [j];

上限のuコンポーネントにも同じことを行う必要があります。







配列では、4つの境界線が次の座標で表されます。



左壁



u [0] [j]、jは壁全体を走ります

v [0] [j]



コードでは次のようになります。

for (j = 0; j <= jmax + 1; j++) { U[0][j] = 0.0; V[0][j] = -V[1][j]; }
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

底壁



u [i] [0]

v [i] [0]

 for (i = 0; i <= imax + 1; i++) { U[i][0] = -U[i][1]; V[i][0] = 0.0; }
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

右壁

間隔を空けたグリッドがあるため、右の壁には、最後の列ではなく最後から2番目の列に沿って通過するuの値の境界があります。したがって、セルのUの値をU [imax] [j]に設定します。



u [imax] [j]

v [imax + 1] [j]

 for (j = 0; j <= jmax + 1; j++) { U[imax][j] = 0.0; V[imax + 1][j] = -V[imax][j]; }
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

上壁

ここで、vの値の場合、境界は最後から2番目の線jmaxに沿って走るので、セルのVの値を設定します-V [i] [jmax]



u [i] [jmax + 1]

v [i] [jmax]

 for (i = 0; i <= imax + 1; i++) { U[i][jmax + 1] = -U[i][jmax]; V[i][jmax] = 0.0; }
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

メインソルバーをテストするために、すべての境界値= 0を設定できます。



また、境界に圧力をかける必要があります。圧力は、速度としてエッジではなくセルの中央に設定されます。したがって、非常に簡単です。境界の圧力は、隣接するセルと同じに設定できます。隣接するセルから値をコピーするだけです。

 for (j = 1; j <= jmax; j++) { //   P[0][j] = P[1][j]; //   P[imax + 1][j] = P[imax][j]; } for (i = 1; i <= imax; i++) { P[i][0] = P[i][1]; P[i][jmax + 1] = P[i][jmax]; }
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

障害物の境界条件

障害物はフラグC_Bで表されます。 それらの境界条件は、外壁の場合と同じ原理で設定されています。圧力の設定の2つの例を示します。

 if (IsObstacle(i, j)) { //     -  -    P   if (IsFluid(i, j - 1)) { P[i][j] = P[i][j - 1]; } //   -  ... else if (IsFluid(i - 1, j)) { P[i][j] = P[i - 1][j]; } // ....... }
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

外周セルの場合、周囲の水セルの値の平均を取得し、たとえば障害物を取得すると、その左上に水があるため、次のように障害物セルの圧力を考慮します。

 if (IsFluid(i - 1, j) && IsFluid(i, j + 1)) { P[i][j] = (P[i][j + 1] + P[i - 1][j]) / 2; }
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

表面の境界条件

表面とその動きは、粒子を使用してモデル化されます。 （このセクションで説明されているように作成された、表面に境界条件を持つソースは、simpletestobstacleフォルダーにあります。）



最初に、粒子はセルごとに4個ずつセルごとに液体で配置されます（各コーナー近くに1個）。 次に、各ステップでの粒子は、古典力学からの単純なオイラー法を使用して移動します。 移動の速度は、セル全体の平均として各セルで取得されます（ただし、補間などを使用しても問題ありません）。

  x = particles[k].x; y = particles[k].y; //        i = (int)(x / delx); j = (int)(y / dely); u = U[i][j]; v = V[i][j]; //        x += delt * u; y += delt * v;
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

各ステップで、粒子のあるセルは水としてラベル付けされ、残りは空のセルであり、基本変数（u、v、p）は計算されません。



主な変数を計算するには、壁に必要なのと同じように、表面と隣接するセルに境界条件を設定する必要があります。 しかし、最初に、水とラベル付けされた細胞のうち、どの細胞が表面に属し、またどの側から空気が入っているかを判断する必要があります。 これらの目的のために、FLAGとFLAGSURFの2つのアレイが使用されます。 1つ目は、セルタイプ（水、空気（空）、障害物）のみが指定されています。 対応するフラグは次のとおりです（略語B-境界F-流体E-空）：

 public const int C_B = 0x0000; // public const int C_F = 0x0010;// public const int C_E = 0x1000;//
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

FLAGSURF配列は、表面のセルを決定するために使用され、残りのセルでは、そこの値= 0です。 この配列のフラグは、セルのタイプだけでなく、空の隣接セルのすべての組み合わせも決定します。 フラグは、結合できるように標準ビットマスクとして作成されます。

FLAGSURFの各値には、4辺（隣接セル）に対応する4ビットが含まれています。

ビットが1に設定されている場合、対応する隣接セルで-空。 0の場合-水があります。

ビット：0000 NSWO 0000 0000-文字は、N（北北）S（南南）W（西東）およびO（西）の4つの側面を示します。



フラグのリスト全体がソースにあるため、値のいくつかの例を示します。



public const int C_W = 0x0200; // 512

バイナリでは、値は0000 0010 0000 0000のようになります

ここで、サイドWに対応するフラグは1に設定されます。これは、現在のセルの左側が空であることを意味します。

同時に、残りの3ビットは0に設定されます。これは、残りの隣接セルが水で満たされることを意味します。







public const int C_SW = 0x0600; // 1536 0000 0110 0000 0000

ここでは、WとSの両側に対応するフラグが1に設定されています。したがって、現在のセルの左と下は空で、残りのセルは水です。







表面セルのタイプを決定し、FLAGSURF配列を埋めるとき、対応するビットはこの方法で1に設定されます。

 //   ,      1 if (FLAG[i-1][j] == GG.C_E) FLAGSURF[i][j] = FLAGSURF[i][j] | GG.C_W;
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

FLAGSURF配列は、主に表面の境界条件を設定するためにのみ必要です。 表面セルの種類が異なると、異なる境界条件が適用されるためです。 既に述べたように、表面セルの隣にある空のセルに境界条件を置く必要があります。また、表面セル自体にも境界条件が必要です。グリッドがずらされており、すべての表面セルがuv変数の計算に含まれないためです。



値を設定する原理は簡単です。 なぜなら 空気圧は水圧よりも1000倍低い-速度を制限せず、移動方向を変えることなく、水を無視して水面に沿って自由に移動させることができます。 もちろん、私のスキームでは、表面張力は考慮されていません。そうでなければ、すべてがはるかに複雑になります。



表面のセルとそれらに接する空のセルの速度UとVの値を、格子がずれていることを忘れずに記録します。

配置の値は、隣接する水セルから取得されますが、複数のセルが存在する可能性があるため、どの隣接セルからこれを取得するかを決定するだけです。







これは、貼り付ける必要のある境界条件を含む画面です。青い四角は水です。 黒いマークは、必要な境界条件です。 ラベルの一部は表面セルに、一部は空のセルに配置されていることに注意してください。 これは、ソルバーコードには次のような条件があるために発生します（ご覧のとおり、UVは右境界壁の前のセルで計算されません）。

  if (IsFluid(i, j) && IsFluid(i+1, j)){ F[i][j] = ... }
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

いくつかの例を見てみましょう。



フラッグケージC_SW

  case GG.C_SW: { U[i][j - 1] = U[i][j]; V[i][j - 1] = V[i][j]; U[i - 1][j] = U[i][j]; V[i - 1][j] = V[i][j]; } break;
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

ここでは、セル自体の値は必要ありません-ソルバーで計算されます。 ただし、隣接する空のセルに値が必要です-なぜなら ソルバーにはU [i] [j-1]などの項があります。



これらの空のセルに最も近い水セルは、そこからのセル[i] [j]であり、UとVの値が取得されます。

図から、V [i-1] [j]の値は、セル[i] [j]と[i-1] [j + 1]の両方から取得できることがわかりますが、一般的な場合、セル[ i-1] [j + 1]は非水である可能性があり、その中のVの値も境界である可能性があり、まだ下に置かれないため、正しいオプションは[i] [j]です。 。



フラッグケージC_W

  case GG.C_W: { U[i - 1][j] = U[i][j]; V[i - 1][j] = V[i][j]; } break;
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

ここではすべてが似ています。



フラグがC_NWのセル

  case GG.C_NW: { V[i][j] = V[i][j - 1]; U[i - 1][j] = U[i][j]; } break;
      
      

        
        
        
      

    
        
        
        
      
      

        
        
        
      

    
     

ここで、セル自体のUの値はソルバーで計算されます。 彼女の右側に水槽があります。 ただし、上のセルは空であるため、Vを下げる必要があります。また、U [i-1] [j]を設定する必要があります。 U [i-1] [j]。



前の場合と同じように、V [i] [j]からではなく、セルV [i] [j-1]から値V [i] [j]を取得します。V[i-1の値] [j]は境界であり、まだ知られていないことが判明する場合があります。



空のセルと表面セルの圧力を0に設定します。これは完全に正しいわけではありませんが、機能します。

表面セルでは、圧力の方程式を解くときにこれらのセルは境界であり、圧力計算はセル内で直接実行されないため、圧力が必要です。

粒子のみの流体運動アルゴリズム

ソースには、タイトルにトラックという名前のオプションがあります。 これは、粒子が液体全体ではなく表面自体のみを移動する自由表面を移動する方法です。 これは、VOFメソッドに多少似ていますが、表面の予備的な再構築がそこで行われ、かなり面倒です。 私の方法では、粒子がセルを離れると、空のマークが付けられ、粒子のない近くの液体セルに粒子が追加されます-表面がどこにあるかがわかります。 粒子が空のセルに入る場合、セルは液体としてマークされます。 もちろん、この方法には不正確な部分がかなりありますが、高速であり、複雑なコーディングを必要としません。

暗黙的な計算スキーム

ソースには暗黙的なソルバーのバージョンもあります-速度の方程式に適用されます。明示的なバージョンとのコードの違いは最小限です。 サンプリングするとき、UおよびV [i] [j]を含むすべての項は、明示的なように右ではなく方程式の左側に移動します。暗黙的（暗黙的）スキームを使用すると、明示的には不可能な大幅に大きな時間ステップを使用できます。

暗黙的については、 http：//math.mit.edu/cse/codes/mit18086_navierstokes.pdfで読むことができます。

3Dバージョン

3Dバージョンでは、すべてが2Dとの類推によって行われます。

F1-F8 + WASDキー、矢印、ER PgUp PgDownを制御してカメラを回転させます。

水源があるシーンの場合-Pキー-水圧をオフにします。

G-より良い表面レンダリング（立方体の代わりに球が使用されます）が、ひどく遅くなります。

Demosフォルダー-xmlファイル形式のシーンとパラメーター（サイズ、タイムステップ、重力、セルあたりのパーティクル数）があります。 HeightMap（高さのマップ）をペイントでペイントすることもできます。サイズは任意です-自動サイズ変更があります。



結論として、2Dおよび3Dのスクリーンショットを提供します。















ソースはここsourceforgeです。











PS CFDに精通している人がいる場合、ソリューションの速度と正確さの点で（特にサーフェスのモデリングの場合）このプロジェクトを一緒に改善することは興味深いでしょう。

私自身はレンダリングを多少レンダリングできますが、数学と物理学は私の主な分野ではありません。 コメントや役立つヒントを喜んでお届けします。