ー 解の精度と計算の労力のトレードオフ ー

講義予定(案)

1. ( 9/ 2) 数値シミュレーションの手続き (テキスト第１章)

2. ( 9/ 9) 偏微分方程式と解析解 (テキスト第２章)

3. ( 9/16) 休講

4. ( 9/30) 差分方程式とそのスキーム (テキスト第３章) +変換 (テキスト第４章)

5. (10/ 7) 計算 (テキスト第５章)+連立一次方程式の解法(テキスト第６章)

6. (10/21) 流れ関数‐ポテンシャルによる解法(テキスト第７章)

7. (10/28) 流速‐圧力を用いた解法 (テキスト第７章)

8. (11/ 4) 熱流体解析と多相流解析

9. (11/11) 乱流の数値解析 by 金子暁子先生

10. (11/18) 数値解析の実際 by 渡辺正先生(JAEA)

11. (11/25) (予備日)

第５章 計算

ー 解の精度と計算の労力のトレードオフ ー

物理現象 (熱流体挙動)

計算モデル (差分方程式)

厳密解

代数方程式 数学モデル

(偏微分方程式)

差分解 数値解

工学への応用

(開発、設計、性能評価) 仮説

(定式化)

変換

(差分近似) (離散化)

数値計算

適切性 安定性

収束性 適合性 Laxの同等定理

解の精度と計算の労力のトレードオフ

一次元Laplace方程式

(0 1)

2 0

2 = ≤ x ≤

dx u d

を境界条件

( )

u

( )

u

のもとで解く。

線形定常熱伝導問題

差分方程式

差分方程式

( ) ( )

1 2

, 1 ,

, 2 , 1 1

+

+ Δ =

= Δ

−

=

m m

i i x i I x I

x L

変域は格子間隔ΔxでI_m+1個に区切られているとする。すなわち 2 0

2

1

1 =

Δ

+

− ₋

+

x

u u

u_{i i i あるいは} − u_i+1 + 2u_i − u_i−1 = 0

格子点 は内点であり は境界点で

ある。境界点は と指定する。

Im

i = 2,3,L, i = 1, i = I_m+2 1

,

0 2

1 = ₊ =

Im

u u

差分方程式の具体的表示

差分方程式を、境界条件を組み込んで具体的に書くと、

0 1

2

0 2

0 2

0 2

0

!

1 1

4 3

2

3 2

=

− +

−

=

− +

−

=

− +

−

=

− +

−

+

+

−

m m

m m

m

I I

I I

I

u u

u u

u

u u

u

u u

L L

L L

差分方程式の行列表示

連立一次方程式を変数ベクトルuと定数ベクトルbを用いて行 列の形に表すと

( ) ( )

⎟⎟

⎟⎟

⎟⎟

⎟⎟

⎠

⎞

⎜⎜

⎜⎜

⎜⎜

⎜⎜

⎝

⎛

−

−

−

−

−

−

−

−

=

=

=

=

2 1

1 2

1 1 2

1

1 2

1

1 2

1 0

0

Im ,

2 1

L L

L L

A

b u

u u

u

b Au

T T

解法の手順

1290 '*computation

1300 '*boundary condition

1320 '*initial guess for nonlinear equation 1340 '*iteration

1360 '*specification of matrix A & vector b 1380 '*Gauss eliminator

1400 '*data shift 1420 '*plotting

1440 '*data transfer 1460 '*exact solution

1480 '*judgement of convergence 1500 '*return to matrix-specification

Program: SHEAT

非線形定常熱伝導問題

熱伝導率がα(u)のように温度の関数である場合の熱流束

( )

q_x

∂

− ∂

= α

熱伝導率α(u)

( )

(

)

α u u

定常非線形熱伝導方程式

(

)

(

)

⎦

⎢ ⎤

⎣

⎡ + x

dx u du dx

d σ

( )

( )

u

差分方程式

差分方程式

( ) ( )

[ ]

1 ,

0

1 2 1

1 1 0

1 2

2 1 1

1 12

1 12

=

=

+ +

+

=

⎟ =

⎠

⎜ ⎞

⎝

⎛

Δ

− − Δ

− Δ

+

± ±

−

− +

+

Im

i i i

i i

i i

i i

u u

u u

x u u

x u u

x

σ σ

χ

χ χ

これらの式を整理すると、

1 0

12 12

12 2 1

1 + ⎜⎝⎛ + ⎟⎠⎞ − =

− _{− − + − + i+}

i i i

i i

i u χ χ u χ u

χ

差分方程式の行列表示

連立一次方程式を行列の形に表すと

( ) ^T _I ^T

m

d u

u u

u d

u

A ⎟

⎠

⎜ ⎞

⎝

= ⎛

=

= , _{1 2 Im} , 0 0 χ +32

χ L L

⎟⎟

⎟⎟

⎟⎟

⎟⎟

⎠

⎞

⎜⎜

⎜⎜

⎜⎜

⎜⎜

⎝

⎛

−

−

−

−

−

−

−

−

=

+

+1 1

4 4

4

3 3

3

2 2

m m

m m

m

I I

I I

I

b a

c b

a c b

a

c b

a

c b

A L L

12 12

12

12 − + +

− = + =

= i b i i c i

a χ χ χ χ

初期推測値

１．行列 を指定するために数値解の予測値

( ) ( )⁰ ( )⁰

3 0

2 , ,

Im

u u

u L

Ax

( )⁰

Ax

( ) ( )

d u

A_x^{0 1} =

解 ^u₂^{( )1 , u}₃^{( )1 ,} _L ^{u( )}_I¹_m を求める。

５．この手順を繰り返す。

( ) ( )

d u

A_x^{k k+1} =

２．これより行列 が仮に定まる

３．行列式 を解き、

４．この第k回目の手順を、 とする。

６．第k回目の解u^(k)とu^(k+1)がほとんど同じであれば、

それを求めるべき解と決心する。

を設定する。

非線形熱伝導問題の厳密解

先の境界条件における非線形熱伝導問題の厳密解は

( )

(

)

u σ

σ σ σ

であるので、数値解と比較できる。

Program: SHEAT

高温の場合に熱伝導が良くなる。

X=1付近の高温領域で、温度勾配が小さくなる

Poisson方程式の数値解

Poisson方程式

(^{0 1})

2 2 2

2 = − ≤ ≤

∂ + ∂

∂

∂ g x y

y u x

u

境界条件は、一辺が１の正方形 の上でu=0とする。差分方程式は

j i j

i j

i j

i j

i j

i u u u u x g

u _{, 1} − _1, + 4 _, + _1, − _{, 1} = Δ ² _,

− _{− − + +}

ここで、Δx=Δyとした。

x y

1 2 3

4 5 6

7 8 9

11 12 13 14 15

16 17

18 19 21 22 23 24 25

20

10

差分方程式の行列表示(1/2)

行列表示

⎟⎟

⎟

⎠

⎞

⎜⎜

⎜

⎝

⎛

⎟ =

⎟⎟

⎠

⎞

⎜⎜

⎜

⎝

⎛

⎟⎟

⎟

⎠

⎞

⎜⎜

⎜

⎝

⎛

3 2 1

3 2 1

33 32

23 22

21

12 11

B B B

U U U

A A

A A

A

A A

( )

(^{i j i j})

A A

i A

ji ij

ii

≠

⎟ =

⎟⎟

⎠

⎞

⎜⎜

⎜

⎝

⎛

−

−

−

=

=

⎟ =

⎟⎟

⎠

⎞

⎜⎜

⎜

⎝

⎛

−

−

−

−

=

, 3 , 2 , 1 ,

1 1

1

3 , 2 , 1 4

1

1 4

1

1 4

x y

1 2 3

4 5 6

7 8 9

11 12 13 14 15

16 17

18 19 21 22 23 24 25

20

10

差分方程式の行列表示(2/2)

行列表示

と定義する。

(^{u u u} ) ^T

U₃ = _{7 8 9}

( _{4 5 6}) ^,

2

u T

u u

U =

( )

( )

( )

T

T

u u

g x u

g x u

u g

x B

u g

x g

x u

g x B

u u

g x u

g x u

u g

x B

24 20

9 2 23

8 2 22

17 7

2 1

19 6

2 5

2 16

4 2 2

18 13

3 2 12

2 2 15

11 1

2 1

+ +

Δ +

Δ +

+ Δ

=

+ Δ

Δ +

Δ

=

+ +

Δ +

Δ +

+ Δ

=

(^{u u u} ) ^T

U₁ = _{1 2 3}

Poisson方程式の厳密解

Poisson方程式の厳密解は展開項数をMEX,NEXとすると

( ) ^∑ ( ) ( ) ( ) ( )

= ⎟⎟ − −

⎠

⎜⎜ ⎞

⎝

⎛

⎟⎟ −

⎠

⎜⎜ ⎞

⎝

⎛

= ^{MEX NEX} −

n m

y n

x n m

m y g

x

u ^,

1 2 ,

0 sin 2 1 cos 2 1

1 2

2 1

2 2

, 4 π π

π π

π

となる。

Program: POILAP

Laplace方程式の厳密解

Laplace方程式の厳密解は、

( ) [ ( ) ]

( ) ( ) ( )

∑^∞

=

−

= −

1

cos sinh sinh

1 1

, 2

m

m

x m y

m m y m

x

u π π

π π

となる。

Program: POILAP

疎な行列(Sparse Matrix )

変数の数：格子分割の数：n×n=nn分割とすると²

係数行列の全要素数： n² × n² = n⁴ 非零の要素数の割合： (5n-4)/ n³

非零の要素数： (n+2(n-1))×n+n×(n-1)×2=(5n-4)n

n=1 100 % n=2 75 % n=5 16.8 % n=10 4.6 % n=20 1.2 % n=50 0.1968 % n=100 0.0496 %

⎟⎟

⎟

⎠

⎞

⎜⎜

⎜

⎝

⎛

⎟ =

⎟⎟

⎠

⎞

⎜⎜

⎜

⎝

⎛

⎟⎟

⎟

⎠

⎞

⎜⎜

⎜

⎝

⎛

3 2 1

3 2 1

33 32

23 22

21

12 11

B B B U

U U A

A

A A

A

A A

( )

(i j i j)

A A

i A

ji ij

ii

≠

⎟ =

⎟⎟

⎠

⎞

⎜⎜

⎜

⎝

⎛

−

−

−

=

=

⎟ =

⎟⎟

⎠

⎞

⎜⎜

⎜

⎝

⎛

−

−

−

−

=

, 3 , 2 , 1 , 1 1 1

3 , 2 , 1 4

1 1 4

1 1 4

非定常の問題： 波動方程式

波動方程式

2 0

2 2

2 =

∂

− ∂

∂

∂

x c u t

u

( ) ( )

( ) ( ) ( )

0 0 , ,

1 2

1 1

2

2 1 0

0 2

, =

∂

∂

⎩⎨

⎧

≤

≤

−

≤

= ≤

t x u x

x

x x x

u

( ) ( )

( ) ([ ) ] [( ) ]

∑^∞

= + +

+

= −

0 2 sin 2 1 cos 2 1

1 2

8 , 1

m

m

ct m

x m m

t x

u π π

π

初期条件

( ) ( )

u

境界条件 厳密解

波動方程式の数値解法

波動方程式に対するCFLスキーム

( )

n i n

i n

i n

i u u u u

u ⁺¹ = − ⁻¹ + λ² ₊₁ + 2 1− λ² + λ² ₋₁

内点を ^{i = 2,3,}_L^{, I}_m ⁻¹ とし、境界点を i = 1, i = I_m とする。

境界条件^u₁^{n = u}_I^{n = 0}

(

)

m

Program: WAVE

電信方程式

電信方程式

1 0

2 2 2

2

2 =

∂

− ∂

∂ + ∂

∂

∂

x u t

u t

u

c σ

減衰項の中心差分近似

電信方程式の差分スキーム

t u u

t

u _iⁿ _iⁿ Δ

= −

∂

∂ ⁺

2 σ 1

σ

(

)

(

)

(

)

Program: WAVE

非定常熱伝導問題の陽解法

非定常熱伝導方程式

(

)

2

2 ≤ ≤

∂

= ∂

∂

∂ x

x u t

u

初期条件 境界条件

( )

u

( )

( )

u

非定常熱伝導方程式に対するFTCSスキーム

(

)

n i n

i u d u u u

u ⁺¹ = + ₊₁ − 2 + ₋₁

境界条件の与え方

( )^{0, 0 0} ( ^0,1,_L)

1 = =

→

=

n u

t u

n

( )^{1,= /} ₁ ⁰ ( ^{= 0,1,}_L)

→

=

∂

∂

− n

u u

x t

u

m

m I

I

n-1 n n+1

1 2 3

x

n-1 n n+1

I_m-2 I_m-1 I_m x

Program: HEAT

Dirichlet条件： Neumann条件：

非定常熱伝導問題の陰解法

時間に後退、空間に中心差分近似するスキーム

( )

(

)

/ ,

2

1 ¹ ₁^{1 2}

1 1

=

−

=

Δ Δ

=

=

− +

+

− ₋^{+ +} ₊⁺

n I

i

x t

d u

du u

d du

m

n i n

i n

i n

i α

境界を i = 1, i = I_m + 2 とすると、点2と点I_m-1では

( )

( ) ¹_{1 1} ¹₂

1 2

1 1 2

1 3 1

2

2 1 2

1

+ +

− +

− +

−

+ +

+

+

= +

+

−

+

=

− +

n I n

I n

I n

I

n n

n n

m m

m

m d u u du

du

du u

du u

d

この場合、 である。非定常熱伝導問題も定常

熱伝導問題と同じく の形の連立一

次方程式の数値解法に帰着される。

1 0

2 1

1ⁿ⁺ = _Iⁿ⁺₊ = u m

u ^{Aun+1 = b} (^{n = 1,2,3,}_L)

境界条件：Neumann条件

格子系は の格子点からなり、境界点は とする。

1 ,

3 , 2 ,

1 −

= I_m

i L

1 ,

1 = −

= i I_m i

1ⁿ = 0 u

x=0における条件：

x=1におけるNeumann条件： ^{1 11} _{= 0}

Δ

− ⁺₋

+

x u u_I^{n n}_I

m m

1 1

1 +

− + = _Iⁿ

n

I_m u _m

u

あるいは

( )

( ) _Iⁿ _I^{n n}_I ( ) _Iⁿ _Iⁿ

n I

n n

n n

n

m m

m m

m

m d u du du d u u

du

u du

u du

u d

1 1

1 1

2 1

1 1 1

2

2 1

1 2

1 3 1

2

1 2

1 2

1

+ − + −

+ − +−

+−

+ +

+

= +

+

−

=

− +

+

−

= +

=

− +

i=2, I_m-1における差分スキーム

Program: HEATIM

表計算ソフトによる熱伝導解析

100 100 100 100 100 100 100 100 100 100

0 51.07237 64.4262 66.64282 65.19423 61.99906 57.19167 49.72089 37.43022 0 0 31.6532 39.86327 39.9896 36.95083 32.13504 25.61034 17.04673 4.261682 0 0 18.33154 22.32332 21.16744 17.95006 13.67962 8.787319 3.538683 0.884671 0 0 9.599884 11.13954 9.804685 7.44554 4.818129 2.370757 0.884245 0.221061 0 0 4.507474 4.936661 3.986152 2.683597 1.479724 0.668902 0.241968 0.060492 0 0 1.896977 1.953679 1.448326 0.878807 0.45038 0.195639 0.069671 0.017418 0 0 0.716155 0.697562 0.483257 0.27722 0.136799 0.058171 0.020536 0.005134 0 0 0.227735 0.213375 0.142508 0.079438 0.038443 0.016165 0.00568 0.00142 0

0 0 0 0 0 0 0 0 0 0

1 3 5 7 9 S1

S5 S9 0

20 40 60 80 100

80-100 60-80 40-60 20-40

0 0-20

4 _{, , 1 , 1}

, 1 ,

1 − + − − =

− u_{i+ j} u_{i− j} u_{i j} u_{i j+} u_{i j−}

(

)

, 4

1 + + − + + + −

= _{i j i j i j i j}

j

i u u u u

u

2 0

2 2

2 =

∂ + ∂

∂

∂

y u x

u

Poisson方程式：