1 小林研朝ゼミ 桑名分 第 - PowerPoint Presentation

Slide1

1

小林研朝ゼミ

桑名分 第

12

回（

2018/06/01

）

流体現象の数値

シミュレーション

キャビティ

流れ

2

当初の予定

次回以降の予定

6

月

5

日（火）　出張

6

月

15

日（金）出張

次回：

6

月

19

日（火）

3

流体に興味がある人

→ふつうに聞いてください

流体

に

は興味ないけど、プログラミングに興味がある人

→言語やテーマは何でも、プログラムを動かしてみると

　いざプログラミングしよう　と思ったとき楽だと思います。

　（ゼロから始めるのと、一から始めるのではだいぶ違う）

流体にもプログラミングにも興味がない人→できるだけ、実生活に馴染みのあるネタを紹介します。　面白いな、と思ったところだけ聞いてください。

モチベーション

4

プログラミングでいう「

Hello world!

」

Cavity

：空洞、くぼみ、凹み

キャビティ流れ　とは

くぼみ上部の流れに引きずられて

くぼみ内部にも流れができる

一断面

新幹線

車両

連結部の騒音の原因

空力騒音は速度

の約

6

乗に比例して

増大

5

具体例に対して、プログラムを動かして可視化までやってみる

細かい

式の導出などは、後日。

C

言語のサンプルプログラムの紹介。

可視化ソフトウェア（無料）の紹介：

ParaView

今日やりたかったこと

具体例に対して、プログラムを動かして可視化までやってみる→手順の紹介のみ。C言語のサンプルプログラムの紹介。

→概要のみ。

使い慣れて

いる可視化ソフトウェア「

Clef3D

」の紹介

→研究室内・ゼミ等で使えるライセンスを持っています。

今日

やること

6

手順１（こんな式を解きます）

連続の式

・・・

x

軸方向

運動量保存を表す式（

Navier

–Stokes

方程式）

・・・

y

軸方向

：密度

：粘性係数

：拡散係数

：時間

独立変数　　　従属変数　　　　　　定数

：座標

1

：座標

2

：

x

軸方向の速度

：

y

軸方向の速度

：圧力

＊色々な液体を同じ式で計算できる

7

手順２（無次元化する場合が多い）

連続の式

・・・

x

軸方向

運動量保存を表す式（

Navier

–Stokes

方程式）

・・・

y

軸方向

：時間

独立変数　　　従属変数　　　　　　定数

：座標

1

：座標

2

：

x

軸方向の速度

：

y

軸方向の速度

：圧力

＊定数を一つにできる

：レイノルズ数

油：

Re

小さい

水：

Re

大きい

8

手順３（時間発展の形に変形）

フラクショナル・ステップ法（他にもいろいろ）

T

ime evolution

時間発展

を計算する必要がある

前頁の形だと計算できない

間違えていました。

9

手順４（差分化・差分近似）

微分を差分（加減乗除）で計算できるようにする

・・・前進差分

・・・後退差分

・・・中心差分

例

この項

を中心差分で近似した例

10

手順５（格子作成、諸条件の設定）

どんな形の領域で、どんな条件で計算するか決める

ごく短い区間に区切って

、各地点において

ご近所の地点での値を使って計算する。

ごく短い時間に区切って

、各瞬間で計算する。

現在の値と、ご近所の地点での値を使って、

一瞬先の未来の値を計算する。

初期条件：

　計算を始める前の状態

境界条件：壁など

11

手順６（計算結果の可視化）

膨大なデータを、人間が見やすい形に表示する

12

プログラムのながれ

Fortran

を

C

に移植中

いい加減な

フォローチャート

朝ゼミ第

6

回

ラプラス

方程式

反復

解法

Initial condition

Time evolution

Boundary condition

Calculate u*, v*

（

A

）

Calculate RHS

（

B

）

Poisson Equation

（

C

）

Boundary

condition

for Pressure

err<EPS

Calculate new u, v

（

D

）

Yes

No

NMAX

回で終了

NPMAX

回で終了

(A)

(C)

(B)

(D)

間違えていました。

13

プログラムのながれ

// boundary condition

for(

i

=0;i<=

MX;i

++){

U[

i][MY]=1.0; //Upper UT[i][MY]=1.0; V[

i

][MY]=0.0;

VT[

i

][MY]=0.0;

U[

i

][0] =0.0; //Lower

UT[

i

][0] =0.0;

V[

i][0] =0.0; VT[i][0] =0.0; } for(j=0;j<=MY;j

++){ U[MX][j]=0.0; //Right UT[MX][j]=0.0;

V[MX][j]=0.0; VT[MX][j]=0.0;

U[0][j] =0.0; //Left UT[0][j] =0.0; V[0][j] =0.0;

VT[0][j] =0.0; }

14

プログラムのながれ

//

Calculate u*, v*

for(i=1;i

<=MX-1;i++){

for(j=1;j<=MY-1;j++){

UT[i][j]=

U[i][j]+DT*( -U[i][j]*(U[i+1][j]-U[i-1][j])/(2.0*DX) -V[i][j]*(U[i][j+1]-U[i][j-1])/(2.0*DY) +1.0/RE*( (U[i+1][j]-2.0*U[i][j]+U[i-1][j])/(DX*DX)

+(U[i][j+1]-2.0*U[i][j]+U[i][j-1])/(DY*DY)

)

);

VT[i][j]=

　　略

}

}

15

プログラムのながれ

//

Calculate the right hand side of Pressure

for(i=1;i

<=MX-1;i++){

for(j=1;j<=MY-1;j++){

RHS[i][j]=((UT[i+1][j

]-

2.0*UT[i][j]+UT[i-1][j])/(2.0

DX

*DX

)

+(VT[i][j+1]-

2.0*VT[i][j]+

VT[i][j-1

])/(

2.0

DY

*DY

))/DT;

}}

間違えていました。

16

プログラムのながれ

//

boundary condition of Pressure

for(

i

=0;i

<=

MX;i

++){

P[

i

][MY]=P[

i

][MY-1]; //Upper

P[

i

][0] =P[

i

][1]; //Lower

}

for(j=0;j

<=

MY;j

++){

P

U

[MX

][j

]=

P

U

[MX-1

][j]; //Right

P

U

[0

][j] =PU[1][j]; //Left}// Calculate the right hand side of Pressurefor(i=1;i<=MX-1;i++){ for(j=1;j<=MY-1;j++){

RHS[i][j]=((UT[i+1][j]-

2.0*UT[i][j]+UT[i-1][j])/(

2.0 DX*DX)

+(VT[i][j+1]-2.0*VT[i][j]+VT[i][j-1])/(2.0 DY

*DY))/DT; }}

間違えていました。

17

プログラムのながれ

for(np=0;np<=

NPMAX;np

++){

略

for(i=1;i

<=MX-1;i++){ for(j=1;j<=MY-1;j++){ PP=P[i][j]; P[i][j]=((P[i+1][j]+P[i-1][j])/(DX*DX) +(P[i][j+1]+P[i][j-1])/(DY*DY)

-RHS[i][j])/(2.0/(DX*DX)+2.0/(DY*DY));

err=err+fabs(PP-P[i][j]);

}

}

if(err<EPS

) break;

}

「

」の形に変形して反復

18

プログラムのながれ

//

Calculate new velocity

for(

i

=1;i

<=MX-1;i++){

for(j=1;j<=MY-1;j++){

U[i][j]=UT[i][j]-DT*(P[i+1][j]-P[i-1][j])/(2.0*DX); V[

i

][j]=VT[

i

][j]-DT*(P[

i

][j+1]-P[

i

][j-1])/(2.0*DY);

}

}

19

実行する

格子の情報が入ったファイルと

流れ場（

U, V, P

）の情報が入ったファイルができる。

可視化ソフトでそれぞれ読み込む。

20

可視化

する

21

可視化

する