数学とか語学とか楽しいよね

フランス語、ドイツ語、ロシア語、アラビア語、オランダ語、英語、スペイン語、ラテン語とか数学とか数値計算(有限要素法、有限体積法、差分法、格子ボルツマン法、数理最適化、C++コード付き)とか勉強したことをまとめます。右のカテゴリーから興味のある記事を探してください。最近はクラシックの名演も紹介しています。Amazonアソシエイトを使用しています。

【Navier-Stokes方程式】SIMPLE法によるNavier-Stokes方程式の離散化

SIMPLE法(Semi-Implicit Method for Pressure. Linked Equations)によるNavier-Stokes方程式の離散化について説明していきます。やっとこさSIMPLE法が理解出来ました。いきなり空間離散化して項を消したりするのでわかりづらいのです。実はこれHSMAC法とほぼ同じなんです!『数値流体工学』のpp.158-162を読んでいてはっとしました。

数値流体工学

数値流体工学


この手法は流速と圧力を直接未知数として計算する手法です。SIMPLE法では圧力は完全に陰的に、流速は半陰的に取り扱われます。また、昔はよく定常状態を計算するために用いられていました。


以下の内容は『流体解析の基礎』のpp.215-220を参考にしています。

流体解析の基礎

流体解析の基礎


簡単のため二次元直交座標系で考えます。まず保存型のNavier-Stokes方程式

 \displaystyle
\begin{eqnarray} 
\frac{\partial u}{\partial t} + \frac{\partial (u^2)}{\partial x} + \frac{\partial (uv)}{\partial y} &=& -\frac{\partial p}{\partial x} + \frac{1}{Re} \left( \frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2}\right) \\
\frac{\partial v}{\partial t} + \frac{\partial (uv)}{\partial x} + \frac{\partial (v)^2}{\partial y} &=& -\frac{\partial p}{\partial y} + \frac{1}{Re} \left( \frac{\partial^2 v}{\partial x^2} + \frac{\partial^2 v}{\partial y^2}\right)
\end{eqnarray}

です。連続式は

 \displaystyle \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0

です。Navier-Stokes方程式と連続式を以下のように時間に対して離散化します。

 \displaystyle
\begin{eqnarray} 
\frac{u^{n+1}-u^n}{\Delta t} &+& \frac{\partial (u^n u^{n+1})}{\partial x} + \frac{\partial (v^n u^{n+1})}{\partial y} \\ &=& -\frac{\partial p^{n+1}}{\partial x} + \frac{1}{Re} \left( \frac{\partial^2 u^{n+1}}{\partial x^2} + \frac{\partial^2 u^{n+1}}{\partial y^2}\right) \\
\frac{v^{n+1}-v^n}{\Delta t} &+& \frac{\partial (u^n v^{n+1})}{\partial x} + \frac{\partial (v^n v^{n+1})}{\partial y} \\&=& -\frac{\partial p^{n+1}}{\partial y} + \frac{1}{Re} \left( \frac{\partial^2 v^{n+1}}{\partial x^2} + \frac{\partial^2 v^{n+1}}{\partial y^2}\right) \\
&&\frac{\partial u^{n+1}}{\partial x} + \frac{\partial v^{n+1}}{\partial y} = 0
\end{eqnarray}

圧力は完全に陰的に(n+1ステップの値を用いる)、流速は半陰的に(nステップとn+1ステップの値を用いる)取り扱われます。半陰的に取り扱うことにより非線型項を線型化し、非線型代数方程式を解くことを回避しています。

さて、上式の運動量方程式をある点Pで離散化(中心差分でも一次の風上差分でもok)すると、以下のように記述できます。つまりP(中心)の周囲の点W(西)、E(東)、N(北)、S(南)の5点の情報を使って離散化できます。

 \displaystyle
\begin{eqnarray} 
a_P u_P &=& a_W u_W + a_E u_E + a_N u_N + a_S u_S + \Delta y (p_E - p_W) + b\\
c_P v_P &=& c_W v_W + c_E v_E + c_N v_N + c_S v_S + \Delta x (p_N - p_S) + d
\end{eqnarray}

ここで、 u、v は未知数(n+1ステップでの値)で、 a、b、c、d は係数です。具体的な係数は別の記事でまとめます。今回はそれよりも「計算の流れ」に注目しましょう。

次に、推測した圧力を  p^* が与えられたとして(例えば  p^n)、上式の離散化式を用いて計算した流速を  u^*,v^* とします。真の値  p^{n+1}、u^{n+1}、v^{n+1} とのずれ(補正量)をそれぞれ  p'、u'、v' とすると

 \displaystyle
\begin{eqnarray} 
p^{n+1} &=& p^* + p'\\
u^{n+1} &=& u^* + u'\\
v^{n+1} &=& v^* + v'
\end{eqnarray}

と書けます。これが圧力の補正式になります。真の値はもとのNavier-Stokes方程式を、 u^*,v^* は上式の離散化式を満たすことに注意して下さい。すなわち

 \displaystyle
\begin{eqnarray} 
a_P u_P^{n+1} &=& a_W u_W^{n+1} + a_E u_E^{n+1} + a_N u_N^{n+1} + a_S u_S^{n+1} + \Delta y (p_E^{n+1} - p_W^{n+1}) + b\\
c_P v_P^{n+1} &=& c_W v_W^{n+1} + c_E v_E^{n+1} + c_N v_N^{n+1} + c_S v_S^{n+1} + \Delta x (p_N^{n+1} - p_S^{n+1}) + d
\end{eqnarray}

 \displaystyle
\begin{eqnarray} 
a_P u_P^* &=& a_W u_W^* + a_E u_E^* + a_N u_N^* + a_S u_S^* + \Delta y (p_E^* - p_W^*) + b\\
c_P v_P^* &=& c_W v_W^* + c_E v_E^* + c_N v_N^* + c_S v_S^* + \Delta x (p_N^* - p_S^*) + d
\end{eqnarray}

です。この2組の式をそれぞれ引き算してやると

 \displaystyle
\begin{eqnarray} 
a_P u_P' &=& a_W u_W' + a_E u_E' + a_N u_N' + a_S u_S' + \Delta y (p_E' - p_W') \\
c_P v_P' &=& c_W v_W' + c_E v_E' + c_N v_N' + c_S v_S' + \Delta x (p_N' - p_S')
\end{eqnarray}

のように補正量の式が得られます。次の手順がSIMPLE法の特徴です。上式の隣接点に関する項(W、E、N、Sの項)を省略してしまいます!すなわち

 \displaystyle
\begin{eqnarray} 
u_P' &=&  \frac{\Delta y}{a_P} (p_E' - p_W') \\
v_P' &=&  \frac{\Delta x}{c_P} (p_N' - p_S')
\end{eqnarray}

としてしまいます。何故こんなことが可能なのでしょうか?補正量は値が収束すれば0になるので項を落としても大丈夫なのです!上式に仮の値を足すと真の値が求まります。すなわち

 \displaystyle
\begin{eqnarray} 
\bar{u}_P &=&  u^* + \frac{\Delta y}{a_P} (p_E' - p_W') \\
\bar{v}_P &=&  v^* + \frac{\Delta x}{c_P} (p_N' - p_S')
\end{eqnarray}

です。これが流速の補正式になります。本来ならば左辺は  u^{n+1}、v^{n+1} となりますが、隣接項を落としているので  \bar{u}、\bar{v} としておきます。このようにして求めた値を連続式に代入します。すると

 \displaystyle
\begin{eqnarray} 
g_P p_P' &=& g_W p_W' + g_E p_E' + g_N p_N' + g_S p_S' + h
\end{eqnarray}

と圧力に関する補正量の式となります。 g は係数で  u^*、v^* より計算可能です(ぐちゃぐちゃするので別の記事で具体的な式を紹介します)。この式を解くと圧力の補正値が求まります。


計算手順は以下のように1~3の繰り返しとなります。

1. nステップでの圧力  \boldsymbol p^n \boldsymbol{(p^*)} を用いて仮の流速  \boldsymbol{(u^*, v^*)} を得る
2. 仮の流速  \boldsymbol{(u^*, v^*)} を用いて圧力の補正量の式を解き、圧力の補正量  \boldsymbol{p'} を求める
3. 圧力の補正量  \boldsymbol{p'} を用いて、流速を補正(流速の補正式)し、圧力を補正(圧力の補正式)し、 \boldsymbol{(\bar{u}, \bar{v}, \bar{p})} を求める
4. 補正量が十分に小さければ  \boldsymbol{(u^{n+1} = \bar{u}, v^{n+1} = \bar{v}, p^{n+1} = \bar{p})} として次のタイムステップの値を得る、補正量が十分に小さくなければ  \boldsymbol{p^* = \bar{p}, u^* = \bar{u}, v^* = \bar{v})} として1に戻り反復を繰り返す


実は圧力は不足緩和しないといけませんし、流速の更新式にも緩和係数を入れるようです。これも別記事で説明します。

反復一回につき、仮の流速の計算と圧力の補正量の計算の際に、合計二回連立方程式を解く必要があります。

空間方向の離散化には、保存型をそのまま用いれば有限体積法を、非保存型を用いれば有限要素法や差分法を用いることができます。普通は有限体積法が用いられています。

【数か国語対照】風立ちぬ

Paul Valéryの詩"Le Cimetière marin"の一節のフランス語原文、英語訳、ドイツ語訳、日本語訳を紹介します。この一節は宮崎駿の『風立ちぬ』で印象的に登場します。


フランス語版。34秒あたりから。


フランス語
Le vent se lève !… Il faut tenter de vivre !

英語
The wind is rising! . . . We must try to live!

ドイツ語
Der Wind hebt an ! ... Man muß zu leben wagen !

日本語
風が立つ。生きねば。


ちなみに堀辰雄の『風立ちぬ』には「風立ちぬ、いざ生きめやも」という詩句が登場しますが、これはポール・ヴァレリーの詩句の直訳ではないと私は考えています。

参考

https://www.jstor.org/stable/40615872?seq=1#page_scan_tab_contents

【有限体積法】2次元非構造有限体積法のコードを作ってみました 三次元図付き

2次元非構造有限体積法のコードを作ってみました。cell-centeredと呼ばれる方法を使っています。解く方程式は2次元の定常移流拡散方程式です。流れ場はx方向が  u=x、y方向が  v=-y と与えます。領域は  (0,1)^2 で、境界条件

 x=0 においてDirichlet条件  \phi=1-y
 y=1 においてDirichlet条件  \phi=0
 y=0 においてNeumann条件  \frac{\partial \phi}{\partial y} = 0
 x=1 においてNeumann条件  \frac{\partial \phi}{\partial x} = 0

です。ファーツィガーの『コンピュータによる流体力学』4.7の計算例と同じ問題です。

コンピュータによる流体力学

コンピュータによる流体力学


コードを作るうえで一番参考になった本は、やはりVersteegの『数値流体力学』です。非構造格子の話は第11章にあります。

数値流体力学 第2版

数値流体力学 第2版

  • 作者: H.K.Versteeg,W.Malalasekera,松下洋介,齋藤泰洋,青木秀之,三浦隆
  • 出版社/メーカー: 森北出版
  • 発売日: 2011/05/31
  • メディア: 単行本(ソフトカバー)
  • クリック: 4回
  • この商品を含むブログを見る

洋書も含めて非構造格子における有限体積法が詳しく説明されている本は本当に少なく、私も作るのに非常に苦労しました。Versteegの本は高いですが値段に見合う情報が載っています。しかも日本語で読めるのは有難いです。

メッシュはFreeFem++で作成したDelaunay(ドロネーまたはデローニー)三角形メッシュを拝借しています。FreeFem++は非常に扱いやすく、簡単に数値実験ができるのでおすすめです。試してみてください。

FreeFem++の紹介 - Qiita


f:id:mutsumunemitsutan:20190418234707p:plain:w500
計算に用いたメッシュ(非構造格子)

まずメッシュですが、節点数が291、セル数が518です。左下に多めに節点を集めています。


次に計算結果です。出力の都合上、隣接点を繋いでいます。また、境界上の点(境界条件)は出力していません。

f:id:mutsumunemitsutan:20190418182612p:plain:w500
非構造格子による計算結果


比較のために構造格子で解いた場合の結果も載せておきます。17×17=289節点、16×16=256セルです。

f:id:mutsumunemitsutan:20190418182536p:plain:w500
構造格子による計算結果

コードはもう少し整理してから公開します。

【gnuplot】コマンドまとめ

gnuplotのコマンドと参考リンクをまとめます。だんだん増えます。


gnuplotでx軸とy軸のスケールを一致させたいとき

set size ratio -1

gnuplotでx軸とy軸のスケールを一致させたいとき - Qiita

【フランス語単語】se voir+過去分詞

se voir+過去分詞:~される
Elle veut réaliser un souhait qu'elle a oublié et est aidée pour cela de Jintan, qui se voit ainsi confiée la tâche très difficile de réunir au complet leur groupe d'amis d'enfance, dissous après l'accident.

【数値計算】LU分解(Doolittle法)C++コード付き

今回は連立一次方程式に対する直接法である、LU分解(Doolittle法)のC++コードを公開します。

今回解く連立一次方程式は

 Ax=b

とあらわされます。ここで、  A 4 \times 4 の行列、 x b 4 次元のベクトルです。 A b が既知で、 x が未知です。具体的には

 \left( \begin{array}{ccc} 8 & 16 & 24 & 32 \\ 2 & 7 & 12 & 17 \\ 6 & 17 & 32 & 59 \\ 7 & 22 & 46 & 105 \end{array} \right) \left( \begin{array}{c} x_1 \\ x_2 \\ x_3 \\ x_4 \end{array} \right) = \left( \begin{array}{c} 80 \\ 38 \\ 114 \\ 180 \end{array} \right)

であり、これを解くと  x_1=1, x_2=1, x_3=1, x_4=1 となります。プログラムでは  A b はテキストファイルから読み込ませています。また、LU分解の結果は一つの配列(行列)LUに格納しています。ではコードを以下に示します。

#include <iostream>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <stdio.h>
#include <time.h>
#include <stdlib.h> 

using namespace std;

//matix vector multiplication Ax=b
inline void mvm(double A[],double x[],double b[], int size)
{
	int i,j;
	
	for(i=0;i<size;i++)
	{
		b[i]=0.;
		for(j=0;j<size;j++)
		{
			b[i]+=A[i*size+j]*x[j];
		}
	}
}

//inner product//
inline double ip(double a[], double b[], int size)
{
	int i;
	
	double value=0.;
	for(i=0;i<size;i++) value+=a[i]*b[i];
	return value;
}

//input//
inline void input(double A[], double b[], int size)
{
	int i,j;
	
	ifstream fo;
	fo.open("A.txt");
	for(i=0;i<size;i++)
	{
		for(j=0;j<size;j++)
		{
			fo>>A[i*size+j];
		}
	}
	fo.close();
	
	fo.open("b.txt");
	for(i=0;i<size;i++)
	{
		fo>>b[i];
	}
	fo.close();
}

//matrix output//
inline void matout(double A[], int size)
{
	int i,j;
	
	for(i=0;i<size;i++)
	{
		for(j=0;j<size;j++)
		{
			cout<<A[i*size+j]<<" ";
		}
		cout<<endl;
	}
	cout<<endl;
}

//vector output//
inline void vecout(double b[], int size)
{
	int i;
	
	for(i=0;i<size;i++)
	{
		cout<<b[i]<<endl;
	}
	cout<<endl;
}

//LU decomposition//
inline void LUd(double A[], double LU[], int size)
{
	int i,j,k;
	double sum;
	
	for(i=0;i<size;i++)
	{
		for(j=0;j<size;j++)
		{
			if(i>j)//l_ij
			{
				sum=0.0;
				for(k=0;k<j;k++) sum+=LU[i*size+k]*LU[k*size+j];
				LU[i*size+j]=(A[i*size+j]-sum)/LU[j*size+j];
			}
			else//u_ij
			{
				sum=0.0;
				for(k=0;k<i;k++) sum+=LU[i*size+k]*LU[k*size+j];
				LU[i*size+j]=A[i*size+j]-sum;
			}
		}
	}
}

//forward substitution and backward substitution//
void inline fb(double LU[], double x[], double y[], double b[], int size)
{
	int i,j;
	double sum;
	
	y[0]=b[0];
	
	//forward//
	for(i=1;i<size;i++)
	{
		sum=0.0;
		for(j=0;j<i;j++)
		{
			sum+=LU[i*size+j]*y[j];
		}
		y[i]=b[i]-sum;
	}
	
	//backward//
	x[size-1]=y[size-1]/LU[(size-1)*size+size-1];
	
	for(i=size-2;i>=0;i--)
	{
		sum=1.0*x[i];
		for(j=i+1;j<size;j++)
		{
			sum+=LU[i*size+j]*x[j];
		}
		x[i]=(y[i]-sum)/LU[i*size+i];
	}
}

main()
{
	int size=4;
	
	double *b=new double[size];
	double *x=new double[size];
	double *A=new double[size*size];
	double *LU=new double[size*size];
	double *y=new double [size];
	
	//input//
	input(A,b,size);
	
	//LU decomposition//
	LUd(A,LU,size);
	
	//forward substitution and backward substitution//
	fb(LU,x,y,b,size);
	
	//output//
	matout(LU,size);
	vecout(x,size);
	
	delete [] b,x,A,LU,y;

    return 0;
}


参考