熱対流場において粘性効果が形状最適化の結果に与える影響に関する考察
Considerations on the Effect of Viscosity Effects on Shape Optimization Results in Thermal Convection Field
学 青木 崇*1
正 倉橋 貴彦*1，正 片峯 英次*2
Takashi AOKI
*1, Takahiko KURAHASHI
*1and Eiji KATAMINE
*2*1 長岡技術科学大学 Nagaoka University of Technology
*2 岐阜工業高等専門学校 National Institute of Technology, Gifu College
代表著者：青木崇
[email protected]
キーワード：逆問題とデータ同化の最新展開，形状最適化，熱対流，有限要素法，力法
1. 緒 言
電子機器においては，機器の発熱による温度上昇を抑えるために熱対流を用いた冷却が行われている．
機器の一例として，電子やイオンの加速を目的としたクライストロンに用いられる「クライストロン電 源」と呼ばれるものがある．これについて，近藤ら(1)の研究においては発熱する電源をヒーターと見做 し，効率的に冷却を行うためのヒーターの位置や形状について実験的に検討を行っている．その研究にお いては，冷却効率向上のためにはヒーターを領域内の下部に広範囲に配置する必要があるという知見が得 られている．これは，領域内で全体的に対流を起こすことで熱の移動が大きくなるためであるとされてい る．一方で，クライストロン電源の冷却に用いられる流体は絶縁油であり，他の流体を冷却に用いた場合 には粘性効果の大小が異なり，最適な配置も異なることが考えられる．
本論文では，粘性効果が形状最適化の結果に与える影響を検討するために，流体の粘性係数を変更し，あ る境界の放熱量を最大化する問題について検討する．まず，非定常自然対流場の支配方程式を示し，随伴変 数法および有限要素法を用いて形状更新に必要な形状勾配関数を導出する．次に，導出した形状勾配関数に 基づいて力法(2)を適用し，設計境界を移動させることにより放熱量を最大化するヒーターの位置および形状 を決定する．
2. 随伴変数法に基づく放熱量最大化問題の定式化 2・1 目的汎関数
放熱量最大化問題の目的汎関数 𝐽 を定義する．最大化問題における目的汎関数は，マイナスの符号を用い て式（1）のように表される(3)．
𝐽 = −1
2∫ ∫ 𝑄𝑡𝑡𝑓 Γ 𝑖𝑞2
0 𝑑Γ𝑑𝑡 = −1
2∫ ∫ 𝑄𝑡𝑡𝑓 Γ 𝑖ℎ̂2(𝜃 − 𝜃𝑓)2
0 𝑑Γ𝑑𝑡 （1）
ここで，𝑄𝑖 は重み定数，𝑞 は各計算ステップにおける放熱量，ℎ̂ は熱伝達係数，𝜃 は温度，𝜃𝑓 は外気温度を
表している．𝑡0, 𝑡𝑓, Γ はそれぞれ数値解析における初期時刻，最終時刻，数値計算境界を表している．重み定
数 𝑄𝑖 は，放熱量を計算する境界では1とし，他の境界では0として与える．
2・2 支配方程式
非定常自然対流場における支配方程式として，Boussinesq近似されたNavier-Stokes方程式，連続の式，熱 伝達方程式を考える．Navier-Stokes 方程式，連続の式は非圧縮性とする．これらの支配方程式は式（2）～
（4）のように表される．
𝑢𝑖̇ + 𝑢𝑗𝑢𝑖,𝑗+1
𝜌𝑝,𝑖− 𝜈𝑢𝑖,𝑗𝑗+ 𝛽𝑔(𝜃 − 𝜃0)𝑒𝑖= 0 in 𝑡 ∈ [𝑡0, 𝑡𝑓] in Ω （2） 𝑢𝑖,𝑖= 0 in 𝑡 ∈ [𝑡0, 𝑡𝑓] in Ω （3）
𝜃̇ + 𝑢𝑖𝜃,𝑖− 𝛼𝜃,𝑖𝑖= 0 in 𝑡 ∈ [𝑡0, 𝑡𝑓] in Ω （4）
ここで，𝑢𝑖 は流速，𝑝 は圧力，𝜃 は温度である．𝜌, 𝜈, 𝛽, 𝑔, 𝜃0, 𝛼 はそれぞれ流体の密度，動粘性係数，体積 膨張率，重力加速度，基準温度，温度拡散率を表している．また，𝑒1= 𝑒3= 0, 𝑒2= −1 である．
2・3 境界条件および初期条件
非定常自然対流場における境界条件および初期条件を式（5）～（12）のように定義する．式（13）は，
領域Ω における面積制約を表している．
𝑢𝑖= 𝑢̂𝑖 at 𝑡 = 𝑡0 in Ω （5）
𝑢𝑖= 𝑢̂𝑖 in 𝑡 ∈ [𝑡0, 𝑡𝑓] on Γ1u and Γdesign （6） 𝑡𝑖= 𝑡̂ = {−𝑖 1𝜌𝑝𝛿𝑖𝑗+ 𝑢𝑖,𝑗} 𝑛𝑗 in 𝑡 ∈ [𝑡0, 𝑡𝑓] on Γ2u （7） 𝜃 = 𝜃̂ at 𝑡 = 𝑡0 in Ω （8） 𝜃 = 𝜃̂ in 𝑡 ∈ [𝑡0, 𝑡𝑓] on Γ1θ （9）
𝑏 = 𝑏̂ = −𝛼𝜃,𝑖𝑛𝑖 in 𝑡 ∈ [𝑡0, 𝑡𝑓] on Γ2θ （10）
𝑞 = ℎ̂(𝜃 − 𝜃̂𝑓) = −𝛼𝜃,𝑖𝑛𝑖 in 𝑡 ∈ [𝑡0, 𝑡𝑓] on Γ3θ （11） 𝑝 = 𝑝̂ at 𝑡 = 𝑡0 in Ω （12）
∑𝑚𝑥𝐴𝑖
𝑖=1 = 𝐴𝑖𝑛𝑖𝑡𝑖𝑎𝑙 in Ω （13）
ここで，𝛿𝑖𝑗, 𝑛𝑗, (∙)̂ はそれぞれクロネッカーのデルタ，境界における単位法線ベクトル，既知関数を表して いる．Γ1u, Γ2u, Γ1θ, Γ2θ, Γ3θ, Γdesign はそれぞれ流れ場におけるDirichlet境界，流れ場におけるNeumann境 界，温度場におけるDirichlet境界，温度場におけるNeumann境界，Robin境界，設計境界を表している．
また，式（13）において，𝑚𝑥 は全要素数，𝐴𝑖𝑛𝑖𝑡𝑖𝑎𝑙 は領域面積である．
2・4 随伴方程式および形状勾配関数の導出
任意の随伴変数 𝑢𝑖∗, 𝑝∗, 𝜃∗を用いることで，目的汎関数 𝐽 はLagrange関数 𝐽∗に拡張され，式（14）で表さ れる．
𝐽∗= 𝐽 + ∫ ∫ 𝑢𝑖∗{𝑢𝑖̇ + 𝑢𝑗𝑢𝑖,𝑗+1
𝜌𝑝,𝑖− 𝜈𝑢𝑖,𝑗𝑗+ 𝛽𝑔(𝜃 − 𝜃0)𝑒𝑖}
Ω 𝑡𝑓
𝑡0 𝑑Ω𝑑𝑡 + ∫ ∫ 𝑝𝑡𝑡𝑓 Ω ∗𝑢𝑖,𝑖
0 𝑑Ω𝑑𝑡
+ ∫ ∫ 𝜃𝑡𝑡𝑓 Ω ∗(𝜃̇ + 𝑢𝑖𝜃,𝑖− 𝛼𝜃,𝑖𝑖)
0 𝑑Ω𝑑𝑡 （14）
𝑢𝑖∗ は随伴流速，𝑝∗ は随伴圧力，𝜃∗ は随伴温度である．紙面の都合上，詳細は省略するが，Lagrange関数
の停留状態（𝛿𝐽∗= 0）を考慮すると，随伴方程式および随伴条件は次のように決定される．
−𝑢𝑖̇ + 𝑢∗ 𝑗,𝑖𝑢𝑗∗− (𝑢𝑗𝑢𝑖∗),𝑗− 𝑝,𝑖∗− 𝜈𝑢𝑖,𝑗𝑗∗ + 𝜃∗𝜃,𝑖= 0 in 𝑡 ∈ [𝑡0, 𝑡𝑓] in Ω （15）
−1
𝜌𝑢𝑖,𝑖∗ = 0 in 𝑡 ∈ [𝑡0, 𝑡𝑓] in Ω （16）
−𝜃∗̇ − (𝜃∗𝑢𝑖),𝑖− 𝛼𝜃,𝑖𝑖∗ + 𝛽𝑔𝑢𝑖∗𝑒𝑖− ∫ ∫ 𝑄ℎ̂𝑡𝑡𝑓 Γ 2(𝜃 − 𝜃𝑓)
0 𝑑Γ𝑑𝑡 = 0 in 𝑡 ∈ [𝑡0, 𝑡𝑓] in Ω （17）
𝑢𝑖∗= 0 at 𝑡 = 𝑡𝑓 in Ω （18） 𝜃∗= 0 at 𝑡 = 𝑡𝑓 in Ω （19）
𝑢𝑖∗= 0 in 𝑡 ∈ [𝑡0, 𝑡𝑓] on Γ1u and Γdesign （20）
𝜃∗= 0 in 𝑡 ∈ [𝑡0, 𝑡𝑓] on Γ1θ （21）
𝑡𝑖∗= 𝑡̂ = {𝑢𝑖∗ 𝑗𝑢𝑖∗+ 𝑝∗𝛿𝑖𝑗+ 𝜈𝑢𝑖,𝑗∗ }𝑛𝑗 in 𝑡 ∈ [𝑡0, 𝑡𝑓] on Γ2u （22） 𝑏∗= 𝑏̂ = {𝜃∗ ∗𝑢𝑖+ 𝛼𝜃,𝑖∗}𝑛𝑖 in 𝑡 ∈ [𝑡0, 𝑡𝑓] on Γ2θ （23）
さらに，設計境界の座標に対するLagrange関数の勾配 𝐺𝑖は，式（24）から導出される．
∫ ∫ 𝑡𝑡𝑡𝑓 Γ 𝑖∗𝛿𝑢𝑖
0 𝑑Γ𝑑𝑡 + ∫ ∫ 𝑏𝑡𝑡𝑓 Γ ∗𝛿𝜃
0 𝑑Γ𝑑𝑡
= ∫ ∫Γ1u+Γ2u𝑡𝑖∗𝛿𝑢𝑖 𝑡𝑓
𝑡0 𝑑Γ𝑑𝑡 + ∫ ∫Γdesign𝑡𝑖∗𝑢𝑖,𝑗𝛿𝑥𝑗 𝑡𝑓
𝑡0 𝑑Γ𝑑𝑡 + ∫ ∫Γ 𝑏∗𝛿𝜃
1θ+Γ2θ 𝑡𝑓
𝑡0 𝑑Γ𝑑𝑡 + ∫ ∫Γdesign𝑏∗𝜃,𝑖𝛿𝑥𝑖 𝑡𝑓
𝑡0 𝑑Γ𝑑𝑡
= ∫ ∫Γ 𝑡𝑖∗𝛿𝑢𝑖
1u+Γ2u 𝑡𝑓
𝑡0 𝑑Γ𝑑𝑡 + ∫ ∫Γ 𝑏∗𝛿𝜃
1θ+Γ2θ 𝑡𝑓
𝑡0 𝑑Γ𝑑𝑡 + ∫Γ (𝐺𝑗𝛿𝑥𝑗+ 𝐺𝑖𝛿𝑥𝑖)
design 𝑑Γ （24）
勾配 𝐺𝑖 は，線形弾性体の外力として考え，変位値を修正勾配 𝐺𝑖∗ として用いる手法，力法(2)を適用するこ
とで変更される．式（25）に示すように，修正された勾配 𝐺𝑖∗を用いて，目標境界上の座標は更新され る．ここで，𝑙 は形状更新ステップ回数，𝜂 は形状更新のステップ幅である．
𝑥𝑖(𝑙+1)= 𝑥𝑖(𝑙)− 𝜂𝐺𝑖∗(𝑙) on Γdesign （25）
2・5 解析手順
本研究における形状最適化は，以下のステップを繰り返すことによって行われる．
Step 1 初期形状を与える．
Step 2 状態方程式を用いて，初期時刻 𝑡0 から最終時刻 𝑡𝑓 の方向へ，流速分布 𝑢𝑖 ，圧力分布 𝑝 ，温度分布
𝜃 を解析する．
Step 3目的汎関数の収束判定を行う．収束判定条件式（ |𝐽(𝑘)− 𝐽(𝑘−1)|⁄𝐽(0)< 𝜀 ）を満たした場合，解析を終 了する．
Step 4 上記のStep 2で得られた流速分布 𝑢𝑖 ，圧力分布 𝑝 ，温度分布 𝜃 および随伴方程式を用いて，最終時
刻 𝑡𝑓 から初期時刻 𝑡0 の方向へ，随伴流速分布 𝑢𝑖∗ ，随伴圧力分布 𝑝∗ ，随伴温度分布 𝜃∗ を解析する．
Step 5 上記のStep 2およびStep 4で得られた結果を用いて，形状勾配関数 𝐺𝑖 を導出する．
Step 6 力法に基づいて形状更新を行い，Step 2に戻る．
3. 解析例
3・1 解析モデルおよび解析条件
解析モデル，解析モデルの有限要素メッシュおよび解析条件をそれぞれ図1，図2，表1に示す．図1に おいて赤線に示す円を初期形状とし，粘性効果が小さい場合と大きい場合の解析条件として 𝜇 = 0.5 [Pa ∙ s]
と 𝜇 = 1 [Pa ∙ s] の2パターンの最適化解析を行う．目的汎関数である放熱量は，Robin境界 Γ3θ において計 算する．また，温度既知境界 Γ1θ を設計境界とし，すべての境界に対してノンスリップ境界条件 𝑢𝑖= 0 [m/s]
を与える．さらに，Neumann境界 Γ2θ においては断熱境界として 𝑞 = 0 [W/m2] を与える．図1において，
𝐿 = 1 [m] , 𝑅 = 0.2 [m] である．
Fig. 1 Computational model. Fig. 2 Finite element mesh for initial shape.
Table 1 Computational conditions.
Number of nodes 2518 Time steps 𝑡𝑓 [s] 3000
Number of elements 4774 Time increment ∆𝑡 [s] 3
Convergence criterion 𝜀 [−] 10−5 Density 𝜌 [kg/m3] 998.2
Initial condition
Temperature 𝜃𝑖𝑛𝑖 [K] 293.15 Kinematic viscosity 𝜈 [m2/s] 𝜇/𝜌 Flow velocity 𝑢𝑖𝑛𝑖 [m/s] 0 Thermal expansion coefficient 𝛽 [1/K] 2.06 × 10−4
Gravity acceleration 𝑔 [m/s2] 9.8
Boundary condition
Temperature 𝜃̂ [K] 298.15 Reference temperature 𝜃0 [K] 293.15 Heat flux 𝑞̂ [W/m2] 0 Outside air temperature 𝜃𝑓 [K] 293.15 Flow velocity 𝑢𝑖 [m/s] 0 Temperature diffusivity 𝛼 [m2/s] 1.43 × 10−4
Heat transfer coefficient ℎ̂ [W/m2∙ K] 1
3・2 解析結果
図3，図4にそれぞれ粘性係数 𝜇 = 0.5 [Pa ∙ s], 𝜇 = 1 [Pa ∙ s] の場合の解析結果を示す．解析結果として，
改善形状における (a)有限要素メッシュ，(b)温度分布，(c)流速・圧力分布を示す．また，図5にそれぞれの 粘性係数における，初期値を100 [%] とした目的汎関数と面積の変化を示す．目的汎関数は，式（1）で示 した目的汎関数の絶対値を示している．形状更新の際には，FreeFem++上の関数「adaptmesh」を用いてリメ ッシュを行った．
図3, 4より，粘性効果の大小にかかわらず𝑥2 方向に広がるように変化していることが確認できる．一般 的に自然対流による流れは下から上へと発生するため，𝑥2 方向に広がることで流れを促進させるような効 果が得られていると考えられる．熱源の中心部においては，粘性効果が小さい場合は𝑥1 方向への変化があ まり見られなかった．一方で，粘性効果が大きい場合は𝑥1方向に括れるような傾向となった．このことか ら，粘性効果が大きい流体の場合は熱源をI字型のように配置することで，領域の下から上へと移動する 流れを妨げないようにすることが望ましいと考えられる．粘性効果が小さい流体の場合は，対流を促進す るために渦の数を増やすような形状の変化になったと考えられる．また，図5において目的汎関数は増加 途中であるが，メッシュがつぶれたため解析を途中で終了している．
4. 結 言
本研究では，放熱量を最大化する問題において，粘性効果が形状最適化の結果に与える影響について検 討を行った．流体の粘性係数のみを変更して最適化を行った結果，増加した放熱量の値はほぼ同じであっ たが異なる形状が得られた．結果より，粘性効果が対流の大きさや渦の数に大きく影響を与えると考えら れるため，冷却に用いる流体に合わせて機器を最適に配置することが重要であると言える．一方，設計境 界の長さを制限していないことによって放熱量が増加していることが考えられるため，設計境界の周長を 考慮して最適化を行う必要がある．
(a) Finite element mesh. (b) Distribution of temperature. (c) Distributions of flow velocity and pressure.
Fig. 3 Finite element mesh and distributions of temperature, flow velocity and pressure at final time for improved shape ( 𝜇 = 0.5 [Pa ∙ s] ).
(a) Finite element mesh. (b) Distribution of temperature. (c) Distributions of flow velocity and pressure.
Fig. 4 Finite element mesh and distributions of temperature, flow velocity and pressure at final time for improved shape ( 𝜇 = 1 [Pa ∙ s]).
Fig. 5 Variation of absolute value of objective functional and area for 𝜇 = 0.5 [Pa ∙ s] and 𝜇 = 1 [Pa ∙ s].
文 献
(1) 近藤力，他，“クライストロン電源における絶縁油の冷却効率の向上”，第4回加速器学会，TP42，和光，(2003)．
(2) 畔上秀幸，“領域最適化問題の一解法”，日本機械学会論文集A編，60 (1994)，No. 574，pp. 1479-1486．
(3) 山田崇恭，西脇眞二，伊賀淳郎，泉井一浩，吉村允孝，“レベルセット法に基づく熱拡散最大化問題のトポロジー 最適化”，日本機械学会論文集C編，75 (2009)，No. 759，pp. 2868-2876．
