平行平板間流れの層流熱伝達（理論解析と数値計算） - 名古屋工業大学

CAE 懇話会中部地区
「流体伝熱基礎講座 (第 6 回)」講義資料（その 2）
平行平板間流れの層流熱伝達
（理論解析と数値計算）
名古屋工業大学田川正人
はじめに
対流は自然対流と強制対流に通常は分類される．自然対流では速度場と温度場が連成す
るのに対して，強制対流では速度場は温度場とは独立に形成される（温度は受動スカラと
みなされる）．本講習の解析対象である強制対流は，熱交換器をはじめいろいろな工業製
品で広く利用されており，その性能向上を図るために活発な研究開発が展開されている．
強制対流は層流／乱流といった流れの状態や，円管や矩形ダクトなど装置の幾何学形状
で大別して議論されるのが普通である．強制対流熱伝達の数値解析では，解析技法が層流
と乱流で相当に異なる．なかでも，乱流熱伝達の数値計算では，k − ε モデルに代表され
る乱流モデルや LES (Large Eddy Simulation)，DNS (Direct Numerical Simulation) など独
特の手法が適用される．これに対して，層流熱伝達の数値解析は比較的単純であり，特殊
な条件下でなら理論解析で熱伝達率を予測できる場合もある．
本稿では，最も基本的な流れの一つである平行平板間の層流を対象として，そこでの強
制対流熱伝達を理論および数値計算で解析する．理論解析は熱伝達という現象の理解を深
める一助となるだけでなく，数値計算の信頼性を検証する際の基準を与える．理論解析を
適用できるのは特殊な場合に限られ，解析対象や境界条件が強く制約される．一方，数値
シミュレーションにはそのような制約が少なく，様々な境界条件の下で解を得ることがで
きる点で優れている．ただし，数値計算の精度および信頼性を適切に評価することが重要
である．
1. 平行平板間の流れと熱伝達
本節では，2 枚の平行平板の間に形成される層流を対象として，固体壁から流体への熱
伝達に関する基礎的事項をまとめておく．
1.1 プラントル数の影響
図 1 に示すように，一様な流れが平板面と平行に流れているとする．この状況で平板を
加熱すると，温度境界層が形成されて，平板先端付近では図 1 のように発達する．温度境
界層の挙動はプラントル数 Pr によって支配され，Pr が 1 より大きい場合（たとえば水や
油）には速度境界層より薄く，逆に 1 より小さい場合（たとえば液体金属）には厚くなる．
一方，空気 (Pr = 0.7) のように Pr が 1 に近い流体中では，温度境界層は速度境界層と同
じように発達する．このように，Pr は速度場と温度場の相似性（アナロジー）の指標と
なる非常に重要な無次元量である．
1
T u
Pr < 1
温度境界層
主流
速度境界層
Pr ~ 1
Pr > 1
Tw
図 1: 速度境界層と温度境界層の関係（プラントル数 Pr の影響）
加熱または冷却（等熱流束壁，等温壁）
y
Tmi
x
速度境界層
壁面摩擦
係数
Tmo
温度境界層 (Pr > 1)
熱伝達率
HFD
0
TFD
x
図 2: 平行平板間（円管内）入口部での壁面摩擦係数と熱伝達率
1.2 平行平板間入口部での流れと熱伝達
平行平板間や円管内部に速度境界層や温度境界層が形成されると，それらは最終的に合
体する．その様子を図 2 に示す．速度境界層や温度境界層が十分発達するまでの領域を
「助走区間」とよぶ．平行平板間の層流では，速度場が十分に発達 (Hydrodynamically Fully
Developed = HFD) すると，速度は流れ方向に変化しなくなり，放物線分布をなす．一方，
温度場の場合には，それが十分発達（Thermally Fully Developed = TFD）した場合でも，
分布形状は下流向かって変化する．その変化の様子は 3 節の理論解析が示すように，壁面
の熱的境界条件に強く依存する．
平行平板間層流における速度境界層と温度境界層の発達過程を図 3 に示す．壁面の熱的
境界条件としては，等熱流束壁 (qw = −λ(∂T /∂y)|y=0 = const.) と等温壁 (Tw = const.)
が代表的である．平行平板間の入口部では，壁面摩擦係数と熱伝達率が非常に大きな値を
とり，下流に発達するにしたがって一定値に漸近する．温度場の境界条件に着目すれば，
図 3 のように入口直後から加熱（あるいは冷却）する場合ばかりでなく，加熱開始点を速
度場が十分発達した領域に設定する場合（途中加熱）や，片側の板のみ加熱する場合（片
側加熱）などが考えられる．このように，単純な流れであっても，熱的境界条件の設定は
多様である．
2
u (y)
速度場が十分に
発達 (HFD)
速度助走区間
Tmi
Tw
Tw
T (x, y)
0
温度場が十分
に発達 (TFD)
温度助走区間
図 3: 平行平板間流や円管流における速度場と温度場の発達過程
2. 理論解析
ここでは，最も基本的な平行平板間の層流熱伝達を理論的に解析し，熱伝達率を求めて
みる．以下では，簡単のために，平行平板間に形成される層流を「チャネル流」とよぶ．
2.1 基礎式
理論解析では非圧縮性流体の定常 2 次元流れを仮定する．これを記述するために，次に
示す連続の式 (1) と運動方程式 (2) を用いる．解析座標系として，流れ方向に x，壁面垂直
方向に y をとり，各方向の速度成分を u, v とする．
∂u ∂v
+
=0
∂x ∂y
(1)
∂u
1 ∂p
∂u
∂2u ∂2u
+v
=−
+ν
+
u
∂x
∂y
ρ ∂x
∂x2 ∂y 2
∂v
∂v
1 ∂p
∂2v ∂2v
u
+v
=−
+ν
+
∂x
∂y
ρ ∂y
∂x2 ∂y 2













(2)
ここで，ρ, ν は流体の密度 [kg/m3 ] と動粘性係数 [m2 /s] である．
温度場の解析には次のエネルギー式 (3) を用いる（エネルギー式の詳細については文献
[4] を参照されたい）．
u
∂T
∂2T
∂T
∂2T
+v
=α
+
∂x
∂y
∂x2
∂y 2
(3)
ここで，α は温度伝導率 [m2 /s] である．
以上の式を一定の境界条件の下で連立して解くことにより，速度成分 u(x, y), v(x, y) と
温度 T (x, y) を得る．強制対流熱伝達では速度場と温度場が連成しない（温度 T が式 (1),
(2) に影響を及ぼさない）ので，速度場の計算は温度場のそれとは独立である．
3
y
u
2H
x
図 4: チャネル流と解析座標
2.2 十分発達した速度場
速度場が十分発達した場合の速度分布は次のように導出される（先の講習の復習になる
ので読み飛ばして下さって結構です）．
解析座標系を図 4 のように設定する．流路中心を y = 0 とし，流路高さを 2H にとる．
すなわち，y = ±H に壁面がある．
速度場が十分発達したときには，次式の関係が成立する．
∂u
= 0, v = 0
∂x
(4)
すなわち，発達流では，流れ方向速度 u は下流（x 方向）に変化せず，また壁面垂直方向
速度 v は零となる．言い換えれば，式 (4) が速度場発達流を定義する．
式 (4) から，運動方程式 (2) は次のように非常に簡単になる．

1 ∂p
∂2u 

0=−
+ν 2 

ρ ∂x
∂y 

1 ∂p



0=−

ρ ∂y
(5)
上式から，u(x, y) → u(y), p(x, y) → p(x) となるので，
ν
d2 u
1 dp
=
dy 2
ρ dx
(6)
を得る．式 (6) を y について 2 回積分することにより，次式を得る．
u=
1
2µ
dp 2
y + C1 y + C2
dx
(7)
y = ±H で u = 0 より，式 (7) は
0=
1
2µ
dp
H 2 ± C1 H + C2
dx
(8)
となる．上式を連立して解けば，未定係数 C1 , C2 が次のように求まる．




C1 = 0
C2 = −
1
2µ
dp

H2 

dx
(9)
4
これを式 (7) に代入することで速度分布 u(y) を得る．結果は次式で与えられる．
u=
1
2µ
dp
1
dp
(y 2 − H 2 ) =
−
(H 2 − y 2 )
dx
2µ
dx
(10)
すなわち，速度分布はチャネル中心 y = 0 について対称な放物線となる．流れは圧力勾配
−(dp/dx) > 0 で駆動されるので，これと粘性係数 µ が速度分布を決定するパラメータと
なる．速度はチャネル中心 y = 0 で最大となり，その速度 umax は
umax = u|y=0 =
1
dp
−
H2
2µ
dx
(11)
で与えられる．そこで，この umax を用いて速度分布を表すと，式 (10) から
u = umax 1 −
y2
H2
(12)
となる．また，断面平均速度 um は次式で計算される．
um ≡
1
2H
H
−H
udy
umax
y3
=
y−
2H
3H 2
H
−H
2
= umax
3
(13)
上式から，チャネル中心の速度 umax は断面平均速度 um の 1.5 倍であることがわかる．な
お，円管内の層流の場合にも発達した速度分布は放物線となるが，円管中心の速度が断面
平均速度の 2 倍になる点でチャネル流とは異なる．
2.3 混合平均温度
温度場の解析に先立って，混合平均温度 Tm を定義しておく．混合平均温度とは，流れ
に温度分布 T (y) があるときに，対流で輸送されるエンタルピーの y 方向平均値に相当す
る量（温度）である．これは次式で定義される．
H
Tm =
−H
H
u(y)T (y)dy
H
−H
=
u(y)dy
−H
u(y)T (y)dy
(14)
2um H
ここで，um は断面平均速度である．また，式 (14) の右辺の分子と分母では，流体の物性
値が一定であることを仮定して ρcp が省略されている．式 (14) に示すように，混合平均温
H
度には流速分布 u(y) が直接関与する．すなわち，通常の平均温度 Tav = −H
T (y)dy とは
物理的意味が異なる．
2.4 十分発達した温度場
チャネル流において，温度場が十分発達したときに成立するいくつかの関係を導出する．
まず，無次元温度 Θ を次式で定義する．
Θ(x, y) =
Tw (x) − T (x, y)
Tw (x) − Tm (x)
(15)
5
ここで，Tw , Tm はそれぞれ壁温と混合平均温度である．十分発達した温度場では Θ(x, y) →
Θ(y) が実現し，式 (15) の無次元温度分布が流れ方向に変化しなくなると考えられる．す
なわち，y 方向について相似な温度分布が実現される．言い換えれば，温度場の完全発達
は次式で定義される．
∂Θ(x, y)
=0
∂x
(16)
そこで，式 (15) を x で偏微分し，その結果を次のように零とする．
∂Θ
1
=
∂x
(Tw − Tm )2
dTw ∂T
−
dx
∂x
(Tw − Tm ) − (Tw − T )
dTw dTm
−
dx
dx
= 0
(17)
上式から，
(Tw − Tm )
∂T
∂x
= (Tw − Tm )
= (T − Tm )
dTw
dTw
dTm
− (Tw − T )
+ (Tw − T )
dx
dx
dx
dTw
dTm
+ (Tw − T )
dx
dx
(18)
を得る．これから，発達した温度場では，温度分布の x 方向の変化率 ∂T /∂x が次式で与
えられる．
∂T
T − Tm
=
∂x
Tw − Tm
dTw
dx
+
Tw − T
Tw − Tm
dTm
dx
(19)
式 (19) に基づいて，壁面の熱的境界条件と温度分布の発達過程の関係を説明できる（2.6
節参照）．
2.5 チャネル流における熱伝達率の定義
平板境界層では，熱伝達率 h は qw ≡ h(Tw − T∞ ) で定義される．これに対して，チャネ
ル流や管内流の発達流では無限遠方の温度 T∞ が存在しないために，熱伝達率は壁温 Tw
と混合平均温度 Tm との温度差を用いて定義される．すなわち，チャネル流（あるいは管
内流）の熱伝達率 h は次式で定義される．
qw ≡ h(Tw − Tm )
(20)
qw は壁面での熱流束（以下，壁面熱流束とよぶ）であり，次式で与えられる．
∂T
qw = −λ
∂y
y=−H


∂T
= λ
∂y

λ：流体の熱伝導率
(21)
y=H
2.6 壁の熱的境界条件と温度分布の関連
壁面の熱的境界条件として 2 つの代表例（等熱流束壁と等温壁）を取り上げて，熱的境
界条件と温度分布 T (x, y) の流れ方向変化の関連を以下に考察する．
6
・等熱流束壁
等熱流束壁は次式で表される．
dqw
=0
dx
(22)
発達流では h が一定となることに注意して，式 (22) の条件を式 (20) に与えると，
dTw dTm
−
=0
dx
dx
(23)
を得る．すなわち，等熱流束壁の条件で十分発達した温度場では
dTw
dTm
=
dx
dx
(24)
が成立する．この結果を式 (19) に代入すれば，温度分布の流れ方向変化率 ∂T /∂x は次の
ように x のみの関数となることがわかる．
∂T
dTw
dTm
=
=
= F unction(x)
∂x
dx
dx
(25)
よって，∂T /∂x は Tw と Tm の x 方向勾配と一致する．式 (25) の関係を図 5(a) に示す．
なお，等熱流束壁の場合には，混合平均温度 Tm は流路入口から勾配一定で変化する．そ
の理由は次のとおりである．
流路奥行き方向（紙面に垂直方向）の長さを W とするとき，流れ方向座標 x と x + dx
の間の熱エネルギーの出入について，次の関係が成り立つ（通常，x 軸方向熱伝導による
熱移動は考慮しなくてよい）．
mc
˙ p Tm (x + dx) = mc
˙ p Tm (x) + 2qw W dx
(26)
ここで，m,
˙ cp はそれぞれ質量流量 [kg/s]，定圧比熱 [J/(kg·K)] を表す．式 (26) と m
˙ =
ρum × W × 2H の関係から，
dTm
qw
=
dx
ρcp um H
(27)
を得る．したがって，図 5(a) に示すように，壁面熱流束一定（qw = const.）の条件では，
全域で dTm /dx = const. となる．
・等温壁
等温壁は次式で表される．
dTw
=0
dx
(28)
この関係を式 (19) に代入すると，
∂T
Tw − T dTm
=
= F unction(x, y)
∂x
Tw − Tm dx
(29)
となる．これから，等温壁の条件で十分発達した温度場の ∂T /∂x には，等熱流束壁の場
合のような簡単な関係は成立せず，∂T /∂x は x と y の関数となる．このため，等温壁条
件での理論解析は容易ではない．式 (29) の関係を図 5(b) に示す．
7
TFD
発達流
Tw−Tm = const.
Tw
Tm （勾配一定）
Tmi
x
0
(a) 壁面熱流束一定で加熱した場合
Tw =const.
Tw
Tm
Tmi
発達流
TFD：十分発達した温度場
(Thermally Fully Developed)
x
0
(b) 壁面温度一定で加熱した場合
図 5: 壁面の加熱条件と温度場の発達過程（等熱流束壁と等温壁の差異）
十分発達したチャネル流における温度分布の x 方向変化を，等熱流束壁と等温壁を対比
して図 6 に示す．図 6 からわかるように，等熱流束加熱では壁温 Tw と温度分布 T の差が
合同になるのに対して，等温加熱では相似になる．等温壁で加熱し続けると，壁面熱流束
は漸減していき，最終的にはチャネル内の温度分布が Tw で一様になる．その結果，壁か
らの伝熱量は零になる．
2.7 等熱流束壁の場合の熱伝達
2.6 節で述べたように，等熱流束壁の場合には式 (25) が成立するので，十分発達した場
合の温度分布を理論的に求めることができる．その導出過程を以下に示す．
温度場が十分発達していれば，エネルギー式 (3) に対して次の境界層近似が適用できる．
∂2T
∂x2
∂2T
∂y 2
(30)
十分発達した速度場では v = 0 となるので，式 (3) は次のように非常に簡単になる．
u
∂T
∂2T
=α 2
∂x
∂y
(31)
式 (25) から，∂T /∂x = dTm /dx が成立するので，式 (31) を以下のように書き換えることが
できる．
（実は，2.6 節で示したように，等熱流束条件で温度場が十分発達すれば，∂T /∂x
8
合同
Tw
y
Tw
qw
x
qw
(a) 等熱流束で加熱 (qw = const.)
相似
Tw
y
Tw
qw
x
qw
(b) 等温で加熱 (Tw = const.)
図 6: 等熱流束加熱と等温加熱における温度分布の特徴
は一定となる．すなわち，∂ 2 T /∂x2 = 0 となるので，境界層近似を導入せずとも式 (31) が
厳密に成り立つ．
）式 (31) と式 (12), (13) から
∂2T
∂y 2
=
3um
2α
dTm
dx
1−
y2
H2
(32)
を得る．式 (32) の右辺を次式で定義される係数 A（x のみの関数）
A≡
3um
2α
dTm
dx
(33)
を用いて簡単にし，y について順次積分すると，
∂T
∂y
= A y−
y3
3H 2
+ C1
y2
y4
T = A
−
2
12H 2
(34)
+ C1 y + C2
(35)
となる．ここで，C1 , C2 は未定係数であり，次の境界条件により決定される．

∂T


=0 
y=0 で
∂y
y = H で T = Tw
(36)



第 1 の条件式は温度分布がチャネル中心 y = 0 で対称となることを，第 2 の条件式は壁面
y = H で温度 T が壁温 Tw と一致することを表している．条件式 (36) を式 (35) に適用す
9
ると，




C1 = 0
H4
H2
−
2
12H 2
Tw = A
+ C2 = A
(37)
5H 2

+ C2 

12
となる．これから，温度分布 T (x, y) が次のとおり決定される．
y2
y4
5H 2
−
−
2
12H 2
12
T (x, y) = A
= Tw (x) −
AH 2
12
y
H
+ Tw (x)
4
−6
y
H
2
+5
(38)
この結果から，壁面熱流束 qw が次のように求まる．
qw = λ
∂T
∂y
=
y=H
2λAH
3
(39)
一方，混合平均温度 Tm (x) は，その定義式 (14) から，
Tm (x) =
=
1
2um H
1
2um H
H
−H
H
−H
u · T dy
y2
3um
1− 2
2
H
· Tw (x) −
AH 2
12
y
H
4
−6
y
H
2
+5
dy (40)
と表される．式 (40) の計算を容易にするために，η ≡ y/H で定義される変数 η を導入し
て変数変換し，dy = Hdη を用いると，式 (40) は次のように書き換えられる．
Tm (x) =
=
3
4
1
−1
1 − η2
3Tw (x)
4
= Tw (x) −
1
−1
Tw (x) −
AH 2 4
η − 6η 2 + 5
12
1 − η 2 dη +
AH 2
16
1
−1
34AH 2
105
dη
(η 2 − 1) η 4 − 6η 2 + 5 dη
(41)
したがって，式 (20), (39), (41) より，熱伝達率 h は
h=
qw
2λAH/3
35 λ
=
=
2
Tw − Tm
34AH /105
17 H
(42)
と求められる．この結果から，ヌセルト数 Nu を代表長さ 2H （平板間距離）を用いて定
義すると，
h · (2H)
λ
35λ 2H
70
=
·
=
= 4.118
17H λ
17
Nu ≡
(43)
10
10 2
速度場・温度場ともに助走区間
温度場のみ助走区間（速度場発達）
Nu
熱流束一定
発達流
(HFD, TFD)
10
4.36
壁温一定
1
10 −3
3.66
10 −2
10 −1
1
x / (D Re Pr)
図 7: 円管内層流における熱伝達（速度助走区間の有無と熱伝達率の変化）
を得る．一方，代表長さとして，無限に広い平行平板 (W → ∞) の水力直径 Dh ，
Dh =
断面積 × 4
2W H × 4
= lim
= 4H
W
→∞
周囲の長さ
2(W + 2H)
(44)
を採用すれば，Nu = 8.24 となる．
以上の理論解析から，平行平板間の層流が等熱流束壁で加熱（冷却）されて，速度場と
温度場がともに十分発達した場合には，ヌセルト数 Nu（熱伝達率）が Nu = 4.12 になる
ことが示された．
2.8 円管内層流の熱伝達
実用面では上述のチャネル流より円管流の方が重要である．円管内層流の熱伝達は，こ
れまでに様々な条件下で詳しく調べられている．結果の一例を図 7 に示す．図 7 では，速
度場と温度場が同時に発達する場合（速度と温度の助走区間），および十分発達した速度
場の中を温度場が発達していく場合（温度助走区間）について，ヌセルト数の x 方向変化
の解析結果が示されている．横軸は，円管内径を D として Gz ≡ (Re · Pr · D)/x で定義さ
れる無次元数（グレツ数 Gz とよばれる）の逆数がとられている．壁面の熱的境界条件は
等熱流束壁と等温壁である．図 7 から，熱伝達率は円管入口から下流に向って急減するも
のの，速度場と温度場がともに十分発達する領域（Gz−1 > 0.05）においては，等熱流束
壁で 4.36，等温壁で 3.66 の一定値に漸近することがわかる．
11
Inlet
Flow
y
0
2H
x
図 8: 数値解析用の座標（座標原点が理論解析の場合と異なる）
3. 数値解析
ここでは，平行平板間の層流（チャネル流）の熱伝達を数値シミュレーションで解析す
る．解析対象は速度場と温度場が同時に発達する場合であり，壁面での熱的境界条件は等
熱流束壁である．上述の理論解析によると，この条件で速度場と温度場が十分に発達した
ときには，Nu = 4.12 となるはずである．
3.1 基礎方程式の無次元化
数値解析にあたり，2.1 節の運動方程式とエネルギー式を無次元化する．まず，u, v, p,
T を次のとおり無次元化する．
u∗ =
u
,
Uin
v∗ =
v
p
(T − Tin )
, p∗ =
, θ=
2
Uin
ρUin
Tin
(45)
ここで，Uin , Tin はそれぞれチャネル入口での速度と温度であり，ともに入口断面におい
て一様とする．数値解析で使用する座標系を図 8 に示す．ここでは，チャネル中心軸につ
いての場の対称性を利用せずに，チャネル全域を解析対象とする．このため，座標原点を
理論解析の位置（図 4）から図 8 に示す下壁上に移動させる．また，空間座標 x, y および
時間 t を次のとおり無次元化する．
x∗ =
y
t
x
, y∗ =
, t∗ =
2H
2H
(2H/Uin )
(46)
以上により無次元化された運動方程式とエネルギー式は次式で表される．
∗
∗
∂u∗
∂p∗
1
∗ ∂u
∗ ∂u
+
u
+
v
=
−
+
∂t∗
∂x∗
∂y ∗
∂x∗ Re
∂ 2 u∗ ∂ 2 u∗
+ ∗2
∂x∗2
∂y
∗
∗
∂v ∗
∂p∗
1
∗ ∂v
∗ ∂v
+
u
+
v
=
−
+
∂t∗
∂x∗
∂y ∗
∂y ∗ Re
∂ 2v∗ ∂ 2v∗
+
∂x∗ 2 ∂y ∗ 2
∂θ
1
∗ ∂θ
∗ ∂θ
+
u
+
v
=
∂t∗
∂x∗
∂y ∗
RePr
∂2θ
∂2θ
+
∂x∗ 2 ∂y ∗ 2













(47)
(48)
ここで，Re ≡ Uin · 2H/ν, Pr = ν/α である．エネルギー式 (48) に現れる Re と Pr の積を
ペクレ数 Pe とよぶ．
12
3.2 壁面熱流束 qw の無次元化
下壁側の壁面熱流束は次式で与えられる．
qw = −λ
∂T
∂y
= −λ
y=0
Tin
2H
∂θ
∂y ∗
(49)
y ∗ =0
ここで，λ は流体の熱伝導率 [W/(m·K)] である．式 (49) から qw は次のように無次元化さ
れる．
qw∗ = qw
2H
λTin
(50)
この結果，無次元化された壁面熱流束 qw∗ は次式で計算される．
qw∗ = −
∂θ
∂y ∗
=−
y ∗ =0
θi, j=1 − θi, j=0
∆y ∗
(51)
式 (51) から下側壁面での温度の境界条件 θ(i, 0) は
θ(i, 0) = θ(i, 1) + qw∗ ∆y ∗
(52)
で与えられる．
[問] 同様にして上側壁面の温度の境界条件 θ(i, Ny + 1) を求めてみよう．
一方，壁面熱流束 qw は熱伝達率 h を用いて次のように表される．
qw = h(Tw − Tm )
(53)
これを無次元化すれば，
qw∗ = Nu(θw − θm )
(54)
となる．ここで，Nu = h · 2H/λ は座標 x におけるヌセルト数であり，局所ヌセルト数と
よばれる（これに対して，x 方向のある区間で Nu を平均したものを平均ヌセルト数 Num
とよぶ）．
[問] 式 (54) を導出してみよう．
式 (54) から局所ヌセルト数は次式で表される．
Nu =
qw∗
θw − θm
(55)
すなわち，ヌセルト数の値を得るには，ある位置 x での壁温 θw と混合平均温度 θm を算出
する必要がある．下側壁面の温度は
θw =
θi, j=0 + θi, j=1
2
(56)
13
で与えられる．一方，混合平均温度を得るには，式 (14) にしたがって u∗ θ を数値積分す
る必要がある．ここでは，格子点 (i, j) での速度 u∗ と温度 θ の積に格子幅 ∆y ∗ を乗算し，
それらを流路断面にわたって合計することで θm を求める．
[問] 無次元された混合平均温度 θm が次式で表されることを示そう．
1
θm =
0
0
u∗ θdy ∗
1
=
∗
u dy
∗
1
0
u∗ θdy ∗
(57)
3.3 計算プログラム
本講座の講師である名古屋工業大学牛島達夫先生が作成されたプログラムを以下に掲載
する．このプログラムで得られた等熱流束壁のヌセルト数は，チャネル後端で Nu = 4.1248
となった．これは 2.7 節の理論値 4.118 と 0.16%の差で一致する．理論値との若干のずれ
を解消するには，y 方向の格子点を増やして計算精度を高めるとともに，x 方向の計算領
域を大きくすることで速度場と温度場を十分に発達させればよい．本プログラムの計算に
より得られた速度場と温度場の発達過程を図 9∼図 12 に示す．
このプログラムを少し改造すれば等温壁の場合も解析できる．そのようにして得られた
結果を円管流の場合（図 7）と比較してみるのもおもしろい．
c
c
c
c
c
c
c
c
c
c
c
c
SMAC (original version by Dr. T. Ushijima)
Simplified Marker and Cell method
Solving Heat Transfer in flow between parallel walls.
implicit double precision(a-h,o-z)
格子数 n
parameter(n=20, nx=n*10, ny=n)
p: 圧力
u, v: 速度
phi: 補正圧力
divup: 予測速度の発散
up, vp: 予測速度
psi: 流れ関数
the: スカラー (e.g. 温度，濃度)
dimension p(0:nx+1,0:ny+1)
+
,u(0:nx,0:ny+1),v(0:nx+1,0:ny)
+
,phi(0:nx+1,0:ny+1),divup(1:nx,1:ny)
+
,up(0:nx,0:ny+1),vp(0:nx+1,0:ny)
+
,psi(0:nx,0:ny+1)
+
c
c
c
c
,the(0:nx+1,0:ny+1),the0(0:nx+1,0:ny+1)
時間ステップ数 loop
loop=20000
レイノルズ数 re
re=50.d0
プラントル数 pr (シュミット数)
pr=0.7d0
壁面熱流束の値
q_w=1.d0
14
c
c
c
c
c
c
格子幅 dx(=dy)
dy=1.d0/dble(n)
dx=dy
流れ方向の格子幅は壁垂直方向の二倍
dx=dy*2.d0
時間ステップ幅 dt
dt=0.002d0
１時間ステップで流体が移流によって飛び出さないようにする
dt=min(dt,0.25*dx)
拡散の影響の考慮 (ここではプラントル数の影響は考慮していない)
dt=min(dt,0.2*re*dx*dx)
write(6,*) ’dt = ’,dt
ddx=1.d0/dx
ddy=1.d0/dy
ddx2=ddx*ddx
ddy2=ddy*ddy
ddt=1.d0/dt
c initial condition
c
icont=1 のときは，既存のデータを用いる．
c
それ以外は，初期化する．全部零！
icont=0
if (icont.eq.1) then
open(unit=9,file=’fort.21’
+
,form=’unformatted’, status=’unknown’)
read(9) u,v,p
close(9)
else
do 131 j=0,ny
do 132 i=0,nx
u(i,j)=1.d0
v(i,j)=0.d0
p(i,j)=0.d0
the(i,j)=0.d0
132
continue
131
continue
end if
do 133 i=0,nx
p(i,ny+1)=0.d0
u(i,ny+1)=0.d0
the(i,ny+1)=0.d0
133 continue
do 134 j=0,ny
p(nx+1,j)=0.d0
v(nx+1,j)=0.d0
the(nx+1,j)=0.d0
134 continue
c
境界条件の設定１
un=0.d0
uw=1.d0
us=0.d0
ue=0.d0
vn=0.d0
vw=0.d0
vs=0.d0
15
ve=0.d0
時間積分開始
do 1000 it=1,loop
境界条件の設定２
boundary condition
do 135 j=0,ny
right wall (east) or outlet
v(nx+1,j)=v(nx,j)
u(nx,j)=u(nx-1,j)
p(nx+1,j)=0.0
c
c
c
c
the(nx+1,j)=2.*the(nx,j)-the(nx-1,j)
c
left wall (west) or inlet
v(0,j)=vw
u(0,j)=uw
p(0,j)=p(1,j)
the(0,j)=0.d0
135
c
c
continue
u(0,ny/2+1)=uw
do 136 i=0,nx
lower wall (south)
u(i,0)=2.*us-u(i,1)
v(i,0)=vs
p(i,0)=p(i,1)
the(i,0)=the(i,1)+q_w*dy !heat flux constant
c
upper wall (north)
u(i,ny+1)=2.*un-u(i,ny)
v(i,ny)=vn
p(i,ny+1)=p(i,ny)
the(i,ny+1)=the(i,ny)+q_w*dy !heat flux constant
136
c
802
801
804
803
continue
温度の計算
do 801 j=1,ny
do 802 i=1,nx
cnvt=(ddx*((the(i,j)+the(i-1,j))*u(i-1,j)
+
-(the(i,j)+the(i+1,j))*u(i,j))
+
+ddy*((the(i,j)+the(i,j-1))*v(i,j-1)
+
-(the(i,j)+the(i,j+1))*v(i,j)))/2.d0
dift=(ddx2*(the(i+1,j)-2.*the(i,j)+the(i-1,j))
+
+ddy2*(the(i,j+1)-2.*the(i,j)+the(i,j-1)))/(pr*re)
the0(i,j)=the(i,j)+dt*(cnvt+dift)
continue
continue
do 803 j=1,ny
do 804 i=1,nx
the(i,j)=the0(i,j)
continue
continue
16
予測段，既知の速度・圧力から速度を予測する
predictor step
for u_ij
(up-u)/dt=-dp/dx-duu/dx-duv/dy+(nabla)ˆ2 u
write(6,*)’up’
do 125 j=1,ny
do 126 i=1,nx-1
c
vij=(v(i,j)+v(i,j-1)+v(i+1,j)+v(i+1,j+1))*0.25
c
cnvu=0.5*ddx*(u(i,j)*(u(i+1,j)-u(i-1,j))
c
+
-abs(u(i,j))*(u(i-1,j)-2.*u(i,j)+u(i+1,j)))
c
+
+0.5*ddy*(vij*(u(i,j+1)-u(i,j-1))
c
+
-abs(vij)*(u(i,j-1)-2.*u(i,j)+u(i,j+1)))
c
cnvu: 移流項の離散化，上記コメントアウトされたものは一次精度
cnvu=ddx*((u(i+1,j)+u(i,j))**2
+
-(u(i-1,j)+u(i,j))**2)/4.d0
+
+ddy*((u(i,j+1)+u(i,j))*(v(i+1,j)+v(i,j))
+
-(u(i,j)+u(i,j-1))*(v(i,j-1)+v(i+1,j-1)))/4.d0
fij=-ddx*(p(i+1,j)-p(i,j))-cnvu
+
+ddx2*(u(i+1,j)-2.d0*u(i,j)+u(i-1,j))/re
+
+ddy2*(u(i,j+1)-2.d0*u(i,j)+u(i,j-1))/re
up(i,j)=u(i,j)+dt*fij
126
continue
c
write(6,603) (up(i,j),i=1,nx-1)
125
continue
c
do 124 j=1,ny
up(0,j)=uw
up(nx,j)=up(nx-1,j)
124
continue
c
up(0,ny/2+1)=uw*1.01
do 224 i=1,nx-1
up(i,0)=2.*us-up(i,1)
up(i,ny+1)=2.*un-up(i,ny)
224
continue
c
for v_ij
c
(vp-v)/dt=-dp/dy-duv/dx-dvv/dy+(nabla)ˆ2 v
write(6,*)’vp’
do 122 j=1,ny-1
do 123 i=1,nx
c
uij=0.25*(u(i,j)+u(i+1,j)+u(i,j+1)+u(i+1,j+1))
c
cnvv=0.5*ddx*(uij*(v(i+1,j)-v(i-1,j))
c
+
-abs(uij)*(v(i-1,j)-2.*v(i,j)+v(i+1,j)))
c
+
+0.5*ddy*(v(i,j)*(v(i,j+1)-v(i,j-1))
c
+
-abs(v(i,j))*(v(i,j-1)-2.*v(i,j)+v(i,j+1)))
c
cnvv: 移流項の離散化，上記コメントアウトされたものは一次精度
cnvv=ddx*((u(i,j+1)+u(i,j))*(v(i+1,j)+v(i,j))
+
-(u(i-1,j+1)+u(i-1,j))*(v(i-1,j)+v(i,j)))/4.d0
+
+ddy*((v(i,j+1)+v(i,j))**2
+
-(v(i,j)+v(i,j-1))**2)/4.d0
gij=-ddy*(p(i,j+1)-p(i,j))-cnvv
+
+ddx2*(v(i+1,j)-2.d0*v(i,j)+v(i-1,j))/re
+
+ddy2*(v(i,j+1)-2.d0*v(i,j)+v(i,j-1))/re
vp(i,j)=v(i,j)+dt*gij
123
continue
c
write(6,603) (vp(i,j),i=1,nx)
122
continue
c
c
c
c
17
c
121
221
c
111
c
112
c
c
c
108
107
境界条件を合わせる
do 121 i=1,nx
vp(i,0)=vs
vp(i,ny)=vn
continue
do 221 j=1,ny-1
vp(0,j)=vw
vp(nx+1,j)=vp(nx,j)
continue
evaluate continuity
write(6,*) ’evaluate continuity’
ic=0
div=0.0
do 112 j=1,ny
do 111 i=1,nx
divup(i,j)=ddx*(up(i,j)-up(i-1,j))
+
+ddy*(vp(i,j)-vp(i,j-1))
div=div+divup(i,j)**2
ic=ic+1
continue
write(6,603) (divup(i,j),i=1,nx)
continue
write(6,*) sqrt(div/dble(ic))
ポアソン方程式を解く
solve the poisson equation (nabla)2 p=(nabla)up/dt by SOR
write(6,*) ’solve the poisson equation for pressure’
initialisation
do 107 i=0,nx+1
do 108 j=0,ny+1
phi(i,j)=0.d0
continue
continue
eps=1.D-6
最大反復数 maxitr
収束させるためには maxitr=nx*ny 程度必要
最大反復数は小さめにとっているので，計算の初めは
収束しないが，定常計算では影響なし．
maxitr=nx*ny/10
c
緩和係数 alpha=1.5∼1.7 程度に設定
alpha=1.7
do 100 iter=1,maxitr
error=0.d0
do 101 j=1,ny
do 102 i=1,nx
rhs=ddt*divup(i,j)
resid=ddx2*(phi(i-1,j)-2.d0*phi(i,j)+phi(i+1,j))
+
+ddy2*(phi(i,j-1)-2.d0*phi(i,j)+phi(i,j+1))
+
-rhs
den=2.d0*(ddx2+ddy2)
dphi=alpha*resid/den
error=max(abs(dphi),error)
phi(i,j)=phi(i,j)+dphi
102
continue
101
continue
do 103 j=1,ny
c
c
C
C
18
103
104
c
100
998
c
c
c
151
150
153
152
161
160
c
c
156
c
603
155
c
c
c
c
phi(0,j)=phi(1,j)
phi(nx+1,j)=0.0
continue
do 104 i=1,nx
phi(i,0)=phi(i,1)
phi(i,ny+1)=phi(i,ny)
continue
収束の判定
if (error.lt.eps) goto 998
continue
continue
write(6,*) ’iter =’, iter, it, error
pause
if (iter.ge.maxitr) write(6,*)’maximum iteration exceeded!’
修正段
corrector step
do 150 j=1,ny
do 151 i=1,nx-1
u(i,j)=up(i,j)-dt*ddx*(phi(i+1,j)-phi(i,j))
continue
continue
do 152 j=1,ny-1
do 153 i=1,nx
v(i,j)=vp(i,j)-dt*ddy*(phi(i,j+1)-phi(i,j))
continue
continue
do 160 j=1,ny
do 161 i=1,nx
p(i,j)=p(i,j)+phi(i,j)
continue
continue
連続の式の満足度チェック（特に必要なし）
check the continuity for n+1 th step
write(6,*) ’check the continuty for n+1 th step’
ic=0
div=0.0
do 155 j=1,ny
do 156 i=1,nx
divup(i,j)=ddx*(u(i,j)-u(i-1,j))+ddy*(v(i,j)-v(i,j-1))
ic=ic+1
div=div+divup(i,j)**2
continue
write(6,603)(divup(i,j),i=1,nx)
format(20(1X,E10.2))
continue
write(6,*) sqrt(div/dble(ic))
*********************************************
混合平均温度の計算（i=nx における値）
*********************************************
iN=nx
tm=0.d0
um=0.d0
do 170 j=1,ny
tm=tm+0.5d0*(u(iN,j)+u(iN-1,j))*the(iN,j)*dy
19
c
170
c
c
c
um=um+0.5d0*(u(iN,j)+u(iN-1,j))*dy
continue
tm=tm/um
tw=0.5d0*(the(iN,0)+the(iN,1))
s_Nu=q_w/(tw-tm)
write(6,*) s_Nu
1000 continue
output 結果の出力
流れ関数の計算
do 491 i=0,nx
psi(i,0)=0.d0
do 492 j=1,ny+1
psi(i,j)=psi(i,j-1)+0.5*dy*(u(i,j-1)+u(i,j))
492
continue
491 continue
c
do 501 j=1,ny
do 502 i=1,nx
x0=dx*(dble(i)-0.5)
y0=dy*(dble(j)-0.5)
u0=0.5*(u(i,j)+u(i-1,j))
v0=0.5*(v(i,j)+v(i,j-1))
p0=p(i,j)
divu0=divup(i,j)
psi0=0.5*(psi(i,j)+psi(i-1,j))
c ｘ座標，ｙ座標，速度 u,v 圧力 P，連続の式，流れ関数，温度
write(10,699) x0,y0,u0,v0,p0,divu0,psi0,the(i,j)
699
format (8(1X,E12.5))
502
501
C
504
503
c
continue
write(10,*)
continue
do 503 j=1,ny-1
do 504 i=1,nx-1
x0=dx*dble(i)
y0=dy*dble(j)
u0=0.5*(u(i,j)+u(i,j+1))
v0=0.5*(v(i,j)+v(i+1,j))
omega=dx*(v(i+1,j)-v(i,j))-dy*(u(i,j+1)-u(i,j))
psi0=0.5*(psi(i,j)+psi(i,j+1))
x 座標，y 座標，速度 u,v，渦度，流れ関数
write(11,699) x0,y0,u0,v0,omega,psi0
continue
write(11,699)
continue
次回の計算のために速度圧力データを保存
write(21) u,v,p
end
20
2.0
y/(2H)=0.5
u / Uin
1.5
1.0
0.5
y/(2H)=0.10
y/(2H)=0.05
0
0
2
4
6
8
10
x / (2H)
図 9: 速度分布の発達過程（中心速度は入口速度の 1.5 倍になる）
1.0
(T - Tin)/Tin
0.8
y/(2H)=0.05
0.6
0.4
y/(2H)=0.5
0.2
0
0
2
4
6
8
10
x / (2H)
図 10: 温度分布の発達過程（等熱流束壁では勾配一定で増大する．図 5(a) 参照）
21
1.0
1.1
y / (2H)
0.8
0.5
0.6
0.3 0.4
0.6
1
1.1
0.9
1.3
0.6
0.4 0.3
0.7 0.8
0.5
0.4 0.3
0.7 0.8
0.9 1
1.1
0.6
1.4
1.4
1.2
1.3
1.4
0.4
1.1
0.2
1.4
1.3
1.4
1.4
1.3
1.2
1.1
0.9 1
0.7 0.8
0.5 0.60.4
0.3
0
0
2
4
6
x / (2H)
8
10
図 11: 速度 u∗ の等高線分布
1.0
0.25
0.1
0.15
0.35
0.55
0.45
0.3
0.6
0.5
0.7
0.2
0.75
0.65
0.4
0.55
0.6
0.05
y / (2H)
0.8
0.4
0.1
0.2
0.2
0
0
0.25
0.35
0.25
2
4
0.4
0.5
0.45
x / (2H)
0.6
0.55
6
0.65
8
0.7
0.75
10
図 12: 温度 θ の等高線分布
おわりに
平行平板間層流（チャネル流）の熱伝達について理論解析と数値解析の手法を概説した．
本解析の対象とした「無限に広い平行平板間に形成される層流」はもちろん実在しない．
しかし，チャネル流は解析上の取り扱いが容易であることに加えて，工業的に重要な円管
流やダクト内流れと同様に下流で完全な発達流が形成される性質をもつ（これが平板上境
界層と本質的に異なる点である）ことから，基礎研究の対象としてしばしば取り上げられ
てきた．壁乱流の研究では特にその傾向が顕著である．理論解析は現象の本質の理解を助
けるだけでなく，データ整理の指針や数値シミュレーションの信頼性検証の基準データを
与えてくれる．ここではその一例を紹介させていただいた．
参考文献
1.
2.
3.
4.
庄司正弘，
「伝熱工学」, (1995), 東京大学出版会．
日本機械学会，JSME テキストシリーズ「伝熱工学」, (2005), 丸善．
White, F. M. “Heat and Mass Transfer,” (1988), Addison-Wesley.
田川正人，
「エネルギー式を巡る」, 伝熱，43 巻 178 号, (2004.1), pp. 26–31.
(http://www.htsj.or.jp/dennetsu/denpdf/2004_01.pdf)
22

Download Report