CFD Online Programming Seminar

流体力学シミュレーション

オンライン・プログラミング講座

都内で実施した、「流体シミュレーション・プログラム製作講座」のオンライン版を開設しました。既にアクセスコードを購入された受講者の方は、講座ページへ
※PayPal支払い後、パスワード発行ページへジャンプしなかった場合、下記のフォームまたはメールにて支払情報(日付など)と共に連絡してください。

目的・対象者

流体シミュレーションプログラムを製作する上で基礎となる計算手法を、エクセル及びプログラミングを含む実習を通じて習得します。数値計算プログラミングに特化した講座で、流体シミュレーション経験の有無は問いません。講座では、亜音速・超音速流体及び非反応非圧縮性流体をシミュレーションするコードを製作することを目的としています。

流体シミュレーションソフトウェアを自作したいが糸口が掴めない方、とりあえず数値計算の指針が欲しい方、一般的な流体シミュレーションソフトがどのようにして流れ場を解いているのか知りたい趣味人、学生、技術者などが主な対象です。自作の流体コードを作りましょう!

本講座の特徴

以前都内で実施した、「流体シミュレーション・プログラム製作講座」の内容をさらに発展させた、WEB資料・サンプルコードがベースのオンライン教材です。以前の講座受講者からのご意見・ご要望を反映し、受講者が各自のペースで知識・技術を習得できるような教材となっています。流体シミュレーションプログラミングを直感的に最短で習得できるように工夫しました。特に講座資料の数式通りに書かれたC言語・Fortranでのサンプルプログラムにより、理論がどのようにシミュレーション内で取り扱われているのかが分かりやすくなっています。以前の講座においても様々なバックグラウンドを持つ受講者の方々が、プログラムを完成させています。目次と講座のサンプルを以下に紹介しています。

受講者の必須予備知識・技能

  • 表計算ソフト(エクセル)の基本操作
    例:http://allabout.co.jp/gm/gc/297746/程度
  • FortranまたはC言語プログラムの基礎(サンプルコードは両言語にて提供)
    以下のコードのうちどちらか一方が何をしているか分かれば十分です。ポインタなどの比較的高度なプログラミング技術は使いません。また、講座内容と一貫性のあるサンプルプログラムを配布するので、現在はプログラムに詳しくないが先に数値計算について知りたい方でも受講可能です。
  • Fortranプログラム

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    C FORTRAN CODE SAMPLE
           INTEGER I,J,N
           PARAMETER(N=10)
           REAL(8) AR(N,N)
           DO J=1,N
             DO I=1,N
               AR(I,J)=DBLE(I*J)
             ENDDO
           ENDDO
           STOP
           END
    C言語プログラム

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    // C CODE SAMPLE
    #include
    #define n 10
    int i,
    double ar[n][n];

    int main(void){
        for(j=0;j<n;j++){
            for(i=0;i<n;i++){
                ar[i][j]=(double)((i+1)*(j+1));
            }
        }
    }

準備物

以下の物を各自準備してください。

  • パソコン
    Windows、Mac又はLinuxのOSで動き、以下をインストールしてあること。基本的に全て無料で揃います。

    1. 表計算ソフト
      例:エクセル、LibreOffice(無料)など。
    2. FortranまたはC言語プログラム開発環境
      VisualC++ 2010 Express Edition (無料、Win, C言語), G95 (無料、Win, Fortran), Xcode (無料、Mac, C言語&Fortran), Intel Compiler (有料, C言語&Fortran)、gcc, gfortran (無料、Mac&Linux; Windowsの場合はCygwinをインストールすることで利用可能)など。サンプルプログラムは両言語対応。
    3. バイナリデータ可視化ソフト
      Paraview(無料)–使い方
      Viewsquare(無料、Windowsのみ)など。

受講料・ライセンス購入

受講料に含まれるもの:
講座資料、拡散・移流拡散シミュレーション・エクセルシート、圧縮性・非圧縮性流体シミュレーション用プログラム (C言語&Fortran)へのアクセス権 (購入後3ヶ月間有効)。
価格:
・学生ライセンス/人:¥1,200
・一般ライセンス/人:¥4,800
団体価格などの質問は以下のフォームから別途お問い合わせください。

※講座資料へのアクセスは購入後少なくとも3ヶ月間保証されています。それ以降のアクセスは不可能となる場合がありますので、手元に資料をダウンロードすることをお勧めします。

下の”今すぐ購入“ボタンからクレジットカード又はPayPalアカウントを使ってお支払いください。支払完了後、講座へのアクセスコードを記載したページへジャンプしますので、しばらくお待ち下さい。ジャンプしない場合、以下のフォームからお問い合わせください。

ライセンスの種類(税込み価格)

学校名(学生ライセンス購入時のみ)



問い合わせ

本講座に関するご意見・ご質問は随時受け付けています。以下の問い合わせフォーム又は flowsquare@gmail.com からお問い合わせください。必要であれば返答は通常数日~1週間以内に行いますが、繁忙期は長くかかる場合もあります。

複雑な事案の場合、リアルタイムでやり取りのできるSkype(ビデオ電話)での相談も受け付けています(有料、予約前払い制)。予約は講座ページ内の問い合わせフォームを使って行ってください。





 送信前にチェックを入れてください(必須)

Back to Top


目次と講座のサンプル(サンプル項目の目次のみリンク表示)

  1. はじめに
    • 1.1 背景
    • 1.2 サンプルプログラムによる数値計算例
  2. 流体方程式計算の基礎
    • 2.1 流体問題への2大アプローチ
    • 2.2 空間分割化(離散化)と格子点の表記
    • 2.3 流れ場を支配する方程式群
    • 2.4 スカラーの輸送方程式
  3. インクの拡散過程のシミュレーション
    • 3.1 インクの拡散過程とそれを支配する方程式
    • 3.2 シミュレーションにおける1階微分
    • 3.3 シミュレーションにおける2階微分
    • 3.4 シミュレーションにおける時間積分(時間発展)
    • 3.5 拡散方程式:数値計算プログラミングに向けた式変形 ★★★
    • 3.6 境界上の格子点の取り扱い(境界条件)
    • 3.7 インクの拡散過程シミュレーションのプログラミング(エクセル)
  4. インクの移流拡散過程のシミュレーション
    • 4.1 インクの移流拡散過程とそれを支配する方程式
    • 4.2 境界上の格子点の取り扱い(境界条件)
    • 4.3 移流拡散方程式:数値計算プログラミングに向けた式変形 ★★★
    • 4.4 インクの移流拡散過程シミュレーションのプログラミング
  5. 2次元圧縮性流体シミュレーション
    • 5.1 2次元圧縮性流体支配方程式
    • 5.2 支配方程式の空間離散化
    • 5.3 オイラー法(1次陽解法)による時間積分
    • 5.4 境界の取り扱い(境界条件)
    • 5.5 補助方程式・初期条件
    • 5.6 空間フィルター(~人口粘性;簡易的)
    • 5.7 時間ステップ間隔の制限
    • 5.8 流れ場の境界条件(例)
    • 5.9 2次元圧縮性流体シミュレーション・プログラムの構成例 ★★★
    • 5.10 2次元圧縮性流体シミュレーション・サンプルプログラム
  6. 2次元非圧縮性流体シミュレーション
    • 6.1 非圧縮性流体の特徴
    • 6.2 非圧縮性流体方程式の数値計算手法
    • 6.3 初期条件、計算条件の設定
    • 6.4 圧力項を無視した運動量方程式と近似速度
    • 6.5 境界条件
    • 6.6 圧力ポアソン方程式とその数値計算手法
    • 6.7 圧力効果による速度修正
    • 6.8 非圧縮性流体シミュレーションプログラムの構成例 ★★★
    • 6.9 2次元費圧縮性流体シミュレーション・サンプルプログラム
  7. 【補足】シミュレーション結果(バイナリデータ)の可視化方法
    • A1. Viewsquareのダウンロードと使い方
    • A2. Paraviewの紹介
  8. 1.はじめに

    1.1 背景

    我々の身の回りには様々な流体現象が存在しています。天気や河川などのような自然現象のみならず、自動車や航空機など流体現象を応用した工業製品なども数多く存在し、これらの流れ場を数値解析により可視化する必要性は日増しに強くなっています。市場に出回る汎用ソフトを用いることで、ある程度流れ場の数値解析を行うことが可能ですが、用いるモデル、計算スキーム、境界条件など、ユーザーの要求全てを満たすには、利用者独自でユーザー定義関数などを通じて汎用ソフトに新たな機能を組み込む必要性が出てくるでしょう。

    また、多くの流体解析ソフトは、ユーザーが指定するどのような入力パラメータに対しても、カラフルで滑らかな結果を出力するために、表には表れない”マジック”を使っており、この計算結果への効果は開発者でないユーザーが見積もることは困難です。流体解析シミュレーションコードを自ら製作することは、これらの問題点を理解するためにひと役買うでしょう。

    本講座では、流体シミュレーションプログラムを製作する上で基礎となる計算手法をエクセル・Fortran/C言語計算プログラミングを含む実習を通じて習得します。数値計算プログラミングに特化した講座で、流体シミュレーションソフトウェアを自作したいが全く糸口が掴めない方、とりあえず数値計算の指針が欲しい方、一般的な流体シミュレーションソフトがどのようにして流れ場を解いているのか知りたい趣味人、学生、技術者などが主な対象です。

    Flowsquareを用いた流体シミュレーション結果の例。左上:多摩川の流速ベクトル、左中:カルマン渦列、左下:自動車周りの流れ、右上:混合層(ケルビン・ヘルムホルツ不安定性)、右下:ジェットエンジン・コンプレッサー、中下:超音速ノズル

    図1.1 Flowsquareを用いた流体シミュレーション結果の例。左上:多摩川の流速ベクトル、左中:カルマン渦列、左下:自動車周りの流れ、右上:混合層(ケルビン・ヘルムホルツ不安定性)、右下:ジェットエンジン・コンプレッサー、中下:超音速ノズル

    Back to Top

    1.2 サンプルプログラムによる数値計算例

    本講座では、初めにエクセルを活用した1次元移流拡散シミュレーションを通じて、流体数値シミュレーションの基礎を学び、その後、Fortran又はC言語による2次元流体シミュレーションプログラムの製作を目指します。具体的には、第5節では圧縮性流体のシミュレーションプログラム、そして第6節では非圧縮性流体のシミュレーションプログラムの基礎とコーディングについて学びます。両節においては、流体の基礎方程式(ナビエ・ストークス方程式)通りに記述されたサンプルコードを参考にすることで、数値流体力学や数値計算プログラミングの知識を深めることができます。

    サンプルコード実行による標準出力と流体数値計算途中結果の例(左:圧縮性流体計算、右:非圧縮性流体計算)

    サンプルコード実行による標準出力と流体数値計算途中結果の例(左:圧縮性流体計算、右:非圧縮性流体計算)。可視化ツールはViewsquareを利用。

    Back to Top

    2.流体方程式計算の基礎

    2.1 流体問題への2大アプローチ

    流体に限らず連続体と呼ばれる物体のとらえ方には二つの考え方があります。一つは、連続体を数多くの微小粒子の集まりととらえる方法、もう一つは連続体を連続体としてとらえる方法です。前者をラグランジュ的、後者をオイラー的な考え方と呼びます。

    • ラングランジュ的
      連続体を微小粒子一つ一つの集まりと見る。それぞれの粒子は自然法則に支配されて動く。例えば、この方法では、連続体を高校物理で習うような放物線を描くボールの集まりと捉えると考えると分かり易い。ラグランジュ的な考えに基づくシミュレーションでは、計算領域に無数の粒子を配置し、それぞれの粒子の軌跡を計算する。以下の図のように、一人の観測者が無数のボール一つ一つの動きを追っている風に。
      ラグランジュ的方法

      図2.1 ラグランジュ的方法

    • オイラー的
      連続体をその名の通り、連続する“場”として見る。例えば、この方法では、対象領域の異なる位置(観測点)に非常に多くの観測者を配置し、それぞれの観測者は自分の分担領域のみの状況を報告する、というような考え方。例えば、以下の例では、各観測者がボールの有る無しを逐一報告することで、各人がボールの軌跡を最初から最後まで追わずともボールの位置が分かる。流体では、各観測者が自分の分担領域の状況(流速など)を報告することで、対象領域全体の流速分布が分かる。これは、連続体をばらばらの粒子の集まりと捉えるラグランジュ法と対象に、“場”の概念を用いている。
      オイラー的方法

      図2.2 オイラー的方法

    本講座では、連続体を場と捉えるオイラー的な概念を用います。現在販売されている多くの流体解析ソフトはオイラー的な概念に基づいた支配方程式を解いています。ちなみに、十分精度のよい数値計算手法を用いて支配方程式を解くと、ラグランジュ的・オイラー的にかかわらず同一のシミュレーション結果が得られます。

    Back to Top

    2.2 空間分割化(離散化)と格子点の表記

    先に述べたように、オイラー的な考え方に基づいて流体場を考えるとき、対象となる(シミュレーション)領域に数多くの観測点を配置する必要があります。この観測点を流体シミュレーションでは格子点と言います。ではどのように格子点を配置すればよいでしょう?格子点の配置には数多くの手法がありますが、ここでは手軽で最も精度のよいデカルト座標系に基づいた格子を用いましょう。

    皆さんの2次元流体シミュレーション領域がx方向(横)長さLx (m)、y方向(縦)長さLy (m)の長方形であったとします。この領域に、x方向とy方向それぞれ等間隔に、Nx個とNy個の格子点を配置すると下の図のような格子点配置となります。

    2次元領域上に配置したデカルト格子

    図2.3 2次元領域上に配置したデカルト格子

    ここで、領域中の総格子点数はNxNyの積となります。シミュレーションの条件としてNxNyを大きくすればするほど、上述のオイラー的概念における観測者の数が増えることになり、より詳細な観測(シミュレーション)が可能となります。上の図の赤い部分を拡大したものが下の図になります。x方向とy方向の格子点の間隔は、それぞれ00042_delta_x00044_delta_yと表記します。一般的に左からi番目、下からj番目の格子点の位置を(i, j)と書くことがあります。また、ある物理量(例:速度など)00046_Qに対して、格子点(i, j)で観測される値を00048_Qijと表記します。この表記方法は本講座において今後使っていくので、ここでしっかりと覚えてください。

    2次元領域上に配置したデカルト格子の拡大図

    図2.4 2次元領域上に配置したデカルト格子の拡大図

    それでは、格子点に関して説明を進めるために、上記の2次元領域を単純化した1次元の領域を考えましょう。上と同じ要領ですが、1次元では領域はx方向にLxの長さを持ち、その領域はNx点の格子点で分割されています。格子点で領域を分割することを専門用語では離散化と呼びます。下の図に、この1次元領域のi点周りを拡大した部分の模式図を示します。i点のすぐ左の点をi-1点、すぐ右の点をi+1点と表記します。この表記は後々頻出するので、慣れておいてください。例えば、i点から2点離れた点は、i-2点とi+2点、3点離れた点は、i-3点とi+3点と表記します。

    1次元領域上に配置した格子

    図2.5 1次元領域上に配置した格子

    具体例を出しましょう。下の図では、Lxの長さを持つ1次元領域を9点用いて分割することを考えます。従って、下の例では総格子点数Nxは9点です。両端の点を入れて9点なので、Lxの長さは8つの00042_delta_xと同じ長さになることはお分かりになるでしょうか?(実際に数えてみてください)

    1次元領域上に配置したデカルト格子と格子間隔の模式図

    図2.6 1次元領域上に配置したデカルト格子と格子間隔の模式図

    これを一般的に考えると、領域Lxの中には、Nx-1個の00042_delta_xが存在することになり、従って、00042_delta_xLxNx-1で割ったものと等価となります。この関係式は下の左式に対応します。また、任意のi番目の格子点の領域左端を距離0としたときの左端からの距離xは下の右式で表すことができます。この2つの式に、上図の値(Nx=9, i=1, 2, 3, …, 9)を代入してみて実際に成り立つか確かめてみてください。

    格子点間隔と領域長さ、格子点数の関係式

    式2.1 格子点間隔と領域長さ、格子点数の関係式

    このように配置した格子点の各点上で、これから紹介する流体を支配する方程式群を解き、各格子点上での物性値(流体速度、圧力、流体密度)などの時間変化を求めます。

    Back to Top

    2.3 流れ場を支配する方程式群

    上記で説明した各格子点において数値的に解く流体力学に関する方程式群は、以下のような偏微分方程式の形を持っています。これらの方程式を一般的に、保存方程式や輸送方程式と呼びます。ここでは、以下の方程式の意味については軽く触れる程度に説明します。

    流体の数値シミュレーション中に求める必要のある物理量は、速度(x方向速度成分のuとy方向速度成分のv)、圧力p、圧縮性流体では、密度ρ(ギリシャ文字のrho; ロー)や全エネルギーE、反応性流体や温度変動を伴う流れでは温度Tなどです。ここで、本小節内での添え字(ik)は上記で説明した格子点位置に関する添え字ではなく、ベクトルの成分を示すことに注意してください。後にすべて書き下すので、分からなければ、添え字に関してはここでは特に意識しなくても結構です。

    下式の質量保存の式は、質量(ρ)の増減に関する方程式です。

    質量保存の式

    式2.2 質量保存の式

    下式の運動量方程式は、x及びy方向の運動量密度(それぞれρuρv)の増減に関する方程式です。

    運動量保存の式

    式2.3 運動量保存の式

    上式2つは、全ての流体シミュレーションにおいて何らかの形で解く必要があります。圧縮性(音速)の流体では、上2式に加え、下のエネルギー(E)の保存に関する偏微分方程式も解く必要があります。

    エネルギー保存の式

    式2.4 エネルギー保存の式

    最後に、上式で支配される流体中を漂うインクの色や匂いの強さをシミュレーションする場合には、その“強さ”に対応するスカラー(SまたはρS)の保存に関する偏微分方程式を解く必要があります。一般的なインクの色や匂いの強さの場合、保存方程式は以下のようになります。

    スカラーの保存式

    式2.5 スカラーの保存式

    繰り返しになりますが、これらの式を流体シミュレーション中に解く必要がありますが、式の意味を完璧に理解する必要は本講座では必要ありません。

    これらの式をじっくり見てみると、各式の各項が似たような形を持っていることにお気づきになるでしょうか?例えば、これらの支配方程式群の各項を以下のように色分けすると分かり易いかもしれません。どの保存方程式においても、赤枠で囲まれた項は、d/dtの形を有しています。この項を一般的に時間項と呼びます。また、青色で囲まれた項は、どの保存方程式においても、1階の空間微分d/dxの形を有しています。この項を移流項と呼びます。紫色の点線で囲まれた項は、少し対応が難しいかもしれませんが、基本的に2階の空間微分d2/dxdxの形を有し、一般的に拡散項や粘性項と呼ばれます。

    質量保存の式

    式2.2b 質量保存の式

    運動量保存の式

    式2.3b 運動量保存の式

    エネルギー保存の式

    式2.4b エネルギー保存の式

    スカラーの保存式

    式2.5b スカラーの保存式

    このように、これらの方程式群はある程度、共通の形を有しています。従って、上式(2.2)–(2.5)のうち、一つの偏微分方程式の解き方が分かれば、残りの偏微分方程式も解けるということになります。そこで、次の節では、上式のうち式(2.5)のスカラーの保存方程式に着目して、この偏微分方程式を数値的に解く方法を学びます。

    Back to Top

    2.4 スカラーの輸送方程式

    スカラーはインクの色や匂いの強さを表す量

     

    先ほど軽く触れたように、スカラーとは流体中を漂うインクの色や匂いの強さを表す量だと考えることができます。そして、前節で説明したように、スカラーの保存方程式もまた、他の方程式群と同様に時間項、移流項、拡散項(粘性項)で構成されています。ここで、式中に現れる物理量は、

    00175_rhoは密度 (kg/m^3)、00176_ukは、各方向の速度成分を示しています。2次元の場合、添え字のkには1または2が入り、00177_u1はx方向速度成分のu00178_u2はy方向速度成分vに対応します。また、式中のDは拡散係数です。次節では、このスカラー輸送方程式を基礎として、流体シミュレーションプログラミングの要となる、偏微分方程式の数値計算手法について学びます。

    スカラーの輸送方程式

    式2.5b スカラーの輸送方程式とその構成

    Back to Top

    【補足】シミュレーション結果(バイナリデータ)の可視化方法

    A1. Viewsquareのダウンロードと使い方

    Viewsquareは、本講座のために開発された無料で使える2次元バイナリデータ可視化ソフトです。講座で用いる格子点数nx×ny格子点を有する出力データを読み込み、速度ベクトル表示及び物理量のカラー分布表示を行います。一般的なWindowsがインストールされたマシンで動作しますが、全てのマシン上での動作確認は致しておりません。ご了承ください。

    Viewsquareのダウンロードはこちらから

    ダウンロードしたViewsquareフォルダ内には、実行ファイル (Viewsquare.exe)と入力ファイル (grid.txt)が入っています。この実行ファイルと同じ階層に、本講座にて実施するシミュレーションの中間生成ファイル(サンプルコードでは、u.datu.rawなど)を出力します。入力ファイルには、いくつか項目があり、これらを適切な値に設定した後、Viewsquare.exeを起動すると出力結果が表示されます。中間生成ファイルが刻々と上書きされるケースでは、起動したViewsquareの画面を右クリックすると、新しいファイルがリロードされます。

    本講座では、サンプルコードをコンパイルすることにより生成される実行ファイルとViewsquare.exegrid.txtを同一フォルダ内に置くことを推奨しています。同名で次々に上書きされるシミュレーション中間結果を右クリックのみで簡単に表示できます。

    grid.txtに入力する項目として、q, u, v, nx, ny, skipbyte, flipがあります。qはカラー表示する物理量の中間生成ファイル名、u及びvは、速度uvが格納されている中間生成ファイルのファイル名、nx及びnyは計算領域のx方向及びy方向格子点数、skipbyteは、中間生成ファイルのヘッダーのバイト数、flipは、x方向及びy方向を反転するか否かのフラグとなっています。以下は、本講座で公開しているFortranサンプルプログラムとC言語サンプルプログラムに適した入力内容です。どちらのケースもカラー表示する物理量はuとなっています。

    Fortranサンプルプログラムで生成される中間生成ファイル表示用入力ファイル

    1
    2
    3
    4
    5
    6
    7
    q        u.dat   // filename
    u        u.dat   // filename
    v        v.dat   // filename
    nx        256   // nx
    ny        128   // ny
    skipbyte    4   // 0 (C), 4 or 8 (Fortran)
    flip        2   // Flip (0:off, 1:x-x, 2:y-y, 3:x-y)

    _

    C言語サンプルプログラムで生成される中間生成ファイル表示用入力ファイル

    1
    2
    3
    4
    5
    6
    7
    q        u.raw   // filename
    u        u.raw   // filename
    v        v.raw   // filename
    nx         256   // nx
    ny         128   // ny
    skipbyte     0   // 0 (C), 4 or 8 (Fortran)
    flip         3   // Flip (0:off, 1:x-x, 2:y-y, 3:x-y)

    Back to Top

    A2. Paraviewの紹介

    Paraviewは、オープン(無料)の多次元バイナリデータ可視化ソフトです。このツールを用いて可視化する場合、Fortranコードの以下の部分、

    1
    OPEN(10,FILE='xx.dat',FORM='UNFORMATTED')

    1
    OPEN(10,FILE='xx.raw',FORM='UNFORMATTED', ACCESS='STREAM')

    と変更してください。変更点はファイル拡張子とACCESSオプションです。C言語コードでは変更の必要はありません。

    出力した中間生成ファイルの可視化方法は、まず、(1) ファイルオープンのアイコンをクリックし、(2) Raw (binary) ファイルを選択し、(3) OKボタンをクリックします。そして、下図を参考に、(4) データのファイルサイズなどを適切に指定し、(5) 表示方法のSliceを選択した後、(6) Applyボタンを押すと、選択した中間生成ファイルがカラー表示されます。

    Paraviewの使い方(1)―(3)

    Paraviewの使い方(1)―(3)

    Paraviewの使い方(4)―(6)

    Paraviewの使い方(4)―(6)

    Back to Top