日本機械学会［No.207-1］北陸信越支部 第57期総会・講演会 講演論文集 [2020.3.8新潟県長岡市]
［No.207-1］日本機械学会北陸信越支部 第57期総会・講演会 講演論文集 [2020.3.8新潟県長岡市]
N034
随伴変数法および有限要素法に基づく熱対流場における形状最適化
青木 崇*1，倉橋 貴彦*2，片峯 英次*3
Shape optimization in thermal convection field based on the adjoint variable
and the finite element methods
Takashi AOKI
*1, Takahiko KURAHASHI, Eiji KATAMINE
*1
Nagaoka University of Technology, Department of Mechanical Engineering Kamitomioka 1603-1, Nagaoka-shi, Niigata, 160-0016 Japan
In this study, we present a shape identification analysis about the cooling efficiency of a heater based on the adjoint variable and the finite element methods. The purpose of this study is to obtain an identified shape so as to approach the target radiation amount of the arbitrary boundary in the target shape. Also, the shape gradient for this identification problem is derived by the first variation of the Lagrange function. Reshaping is accomplished using a traction method that is proposed as a solution to the domain optimization problems. The traction method is an algorithm for non-parametric shape optimization problems, and the smooth shape can be obtained by using this method. Based on these considerations, we have constructed a program code for shape identification analysis in thermal convection field. FreeFem++ was used to calculate the identified shape.
Key Words : Shape identification, Finite element method, Adjoint variable method, Thermal convection field, Natural
convection field.
1. 序 論
電子機器や熱交換器などにおいて，性能改善を目的とした形状設計は非常に重要である．例えば，電子機器に おいては発熱による電子機器本体の温度上昇が問題となっており，そのために冷却効率を高くする，つまり放熱 量を大きくすることが求められる．一般的に，電子機器においては温度上昇を抑えるために撹拌機などの冷却機 構が設置されている．しかし，機器のメンテナンス性やコンパクト性を考慮すると，そのような冷却機構の設置 はあまり望ましくないと言える．そこで，上記の性能を向上させるような例として，熱対流の中でも特に攪拌す る機構などを必要としない，自然対流による冷却を利用した機器の設計について取り上げる．
自然対流を利用した電子機器の例の
1
つとして，「クライストロン電源」と呼ばれるものがある．これについ て，近藤ら(1)の研究においては発熱する電源をヒーターと見做し，ヒーターの位置や形状について実験的に検討 を行っている．近年では，数値計算によるシミュレーションが発展してきており，前述した検討は数理設計の立 場から位置および形状を決定するような形状決定問題として捉えることができる．本研究では，ヒーターの形状 設計に応用するかたちで，自然対流場の形状設計問題について取り上げる．これまでに，自然対流場の形状設計 問題に関して，トポロジー最適化における流速最大化問題(2)や温度分布を規定する問題(3)が行われている．しか し，これらの問題は定常問題，あるいは断熱境界を設計境界として取り扱う問題であった．本論文では，ある任意の形状における境界の放熱量を目標値として設定し，ヒーターの位置および形状を変化 させることで，境界の放熱量との
2
乗誤差を最小化することを目的とした放熱量規定問題について検討する．ま*1 学生員，長岡技術科学大学（〒940-2188 新潟県長岡市上富岡町1603-1）
*2 正員，長岡技術科学大学
*3 正員，岐阜工業高等専門学校 E-mail: [email protected]
ず，非定常自然対流場の支配方程式を示し，随伴変数法および有限要素法を用いて，形状更新に必要な形状勾配 関数を導出する．次に，導出した形状勾配関数に基づいて力法(4)を適用し，ヒーターの位置および形状を決定す る．
2. 随伴変数法に基づく放熱量規定問題の定式化 2・1 目的汎関数
放熱量規定問題の目的汎関数 𝐽 を，式（1）のように定義する．この問題では，ある任意の境界における放熱量 を目標の放熱量に近づけることを考える．
𝐽 =
12
∫ ∫ 𝑄
𝑡𝑡𝑓 𝛤 𝑖(𝑞 − 𝑞
(𝑡𝑎𝑟𝑔𝑒𝑡))
20
𝑑𝛤𝑑𝑡
（1
）ただし，
𝑞 = ℎ̂(𝜃 − 𝜃̂
𝑓)
（2
）𝑞
(𝑡𝑎𝑟𝑔𝑒𝑡)= ℎ̂(𝜃
(𝑡𝑎𝑟𝑔𝑒𝑡)− 𝜃̂
𝑓)
（3）である．ここで，
𝑄
𝑖は重み定数， 𝑞 は各計算ステップにおける放熱量， 𝑞
(𝑡𝑎𝑟𝑔𝑒𝑡)は目標とする放熱量， ℎ̂ は熱伝達
係数，𝜃 は温度，𝜃(𝑡𝑎𝑟𝑔𝑒𝑡)は目標とする温度，𝜃̂𝑓は外気温度を表している．𝑡
0, 𝑡
𝑓, 𝛤 はそれぞれ数値解析における
初期時刻，最終時刻，数値計算境界を表している．重み定数 𝑄𝑖は，放熱量を規定する境界では 1
とし，他の境界 では0
として与える．2・2 支配方程式
非定常自然対流場における支配方程式として，
Boussinesq
近似されたNavier-Stokes
方程式，連続の式，熱伝達 方程式を考える．Navier-Stokes方程式，連続の式は非圧縮性とする．これらの支配方程式は式（4）～（6）のよ うに表される．𝑢
𝑖̇ + 𝑢
𝑗𝑢
𝑖,𝑗+
1𝜌
𝑝
,𝑖− 𝜈𝑢
𝑖,𝑗𝑗+ 𝛽𝑔(𝜃 − 𝜃
0)𝑒
𝑖= 0 in 𝑡 ∈ [𝑡
0, 𝑡
𝑓] on 𝛺
（4
）
𝑢
𝑖,𝑖= 0 in 𝑡 ∈ [𝑡
0, 𝑡
𝑓] on 𝛺
（5）𝜃̇ + 𝑢
𝑖𝜃
,𝑖− 𝛼𝜃
,𝑖𝑖= 0 in 𝑡 ∈ [𝑡
0, 𝑡
𝑓] on 𝛺
（6
） ここで，𝑢𝑖は流速，𝑝 は圧力である．𝜌, 𝜈, 𝛽, 𝑔, 𝜃
0, 𝛼 はそれぞれ密度，動粘性係数，体積膨張率，重力加速度，
基準温度，温度拡散率を表している．また，𝑒1
= 𝑒
3= 0, 𝑒
2= −1 である．
2・3 境界条件および初期条件
非定常自然対流場における境界条件，初期条件を式（7）～（14）のように定義する．
𝑢
𝑖= 𝑢 ̂
𝑖at 𝑡 = 𝑡
0in 𝛺
（7
）
𝑢
𝑖= 𝑢 ̂
𝑖in 𝑡 ∈ [𝑡
0, 𝑡
𝑓] on 𝛤
1𝑢and 𝛤
𝑑𝑒𝑠𝑖𝑔𝑛 （8）
𝑡
𝑖= 𝑡 ̂ = {−
𝑖 1𝜌
𝑝𝛿
𝑖𝑗+ 𝑢
𝑖,𝑗} 𝑛
𝑗in 𝑡 ∈ [𝑡
0, 𝑡
𝑓] on 𝛤
2𝑢 （9
）
𝜃 = 𝜃̂ at 𝑡 = 𝑡
0in 𝛺
（10）
𝜃 = 𝜃̂ in 𝑡 ∈ [𝑡
0, 𝑡
𝑓] on 𝛤
1𝜃 （11
）
𝑏 = 𝑏̂ = −𝛼𝜃
,𝑖𝑛
𝑖in 𝑡 ∈ [𝑡
0, 𝑡
𝑓] on 𝛤
2𝜃 （12
）𝑞 = ℎ̂(𝜃 − 𝜃̂
𝑓) = −𝛼𝜃
,𝑖𝑛
𝑖in 𝑡 ∈ [𝑡
0, 𝑡
𝑓] on 𝛤
3𝜃 （13）
𝑝 = 𝑝̂ at 𝑡 = 𝑡
0in 𝛺
（14）ここで，𝛿𝑖𝑗
, 𝑛
𝑗, (∙) ̂ はそれぞれクロネッカーのデルタ，境界における単位法線ベクトル，既知関数を表してい
る．また，𝛤1𝑢, 𝛤
2𝑢, 𝛤
1𝜃, 𝛤
2𝜃, 𝛤
3𝜃はそれぞれ流れ場における Dirichlet
境界，流れ場におけるNeumann
境界，温度場における
Dirichlet
境界，温度場におけるNeumann
境界，Robin
境界を表している．2・4 随伴方程式および形状勾配関数の導出
任意の随伴変数 𝑢𝑖∗
, 𝑝
∗, 𝜃
∗を用いることで，目的汎関数 𝐽 は Lagrange
関数 𝐽∗に拡張され，式（15）で表される．
𝐽
∗= 𝐽 + ∫ ∫ 𝑢
𝑖∗{𝑢
𝑖̇ + 𝑢
𝑗𝑢
𝑖,𝑗+
1𝜌
𝑝
,𝑖− 𝜈𝑢
𝑖,𝑗𝑗+ 𝛽𝑔(𝜃 − 𝜃
0)𝑒
𝑖}
𝛺 𝑡𝑓
𝑡0
𝑑𝛺𝑑𝑡 + ∫ ∫ 𝑝
𝑡𝑡𝑓 𝛺 ∗𝑢
𝑖,𝑖0
𝑑𝛺𝑑𝑡
+ ∫ ∫ 𝜃
𝑡𝑡𝑓 𝛺 ∗(𝜃̇ + 𝑢
𝑖𝜃
,𝑖− 𝛼𝜃
,𝑖𝑖)
0
𝑑𝛺𝑑𝑡
（15）𝑢
𝑖∗は随伴流速，𝑝
∗は随伴圧力，𝜃
∗は随伴温度である．
紙面の都合上，詳細は省略するが，Lagrange関数の停留状態（𝛿𝐽∗
= 0）を考慮すると，随伴方程式と随伴条
件は次のように決定される．−𝑢
𝑖̇ + 𝑢
∗ 𝑗,𝑖𝑢
𝑗∗− (𝑢
𝑗𝑢
𝑖∗)
,𝑗
− 𝑝
,𝑖∗− 𝜈𝑢
𝑖,𝑗𝑗∗+ 𝜃
∗𝜃
,𝑖= 0 in 𝑡 ∈ [𝑡
0, 𝑡
𝑓] on 𝛺
（16）−
1𝜌
𝑢
𝑖,𝑖∗= 0 in 𝑡 ∈ [𝑡
0, 𝑡
𝑓] on 𝛺
（17
）−𝜃
∗̇ − (𝜃
∗𝑢
𝑖)
,𝑖− 𝛼𝜃
,𝑖𝑖∗+ 𝛽𝑔𝑢
𝑖∗𝑒
𝑖+ ∫ ∫ 𝑄ℎ̂
𝑡𝑡𝑓 𝛤 2(𝜃 − 𝜃
(𝑡𝑎𝑟𝑔𝑒𝑡))
0
𝑑𝛤𝑑𝑡 = 0 in 𝑡 ∈ [𝑡
0, 𝑡
𝑓] on 𝛺
（18
）𝑢
𝑖∗= 0 at 𝑡 = 𝑡
𝑓in 𝛺
（19
）𝜃
∗= 0 at 𝑡 = 𝑡
𝑓in 𝛺
（20）𝑢
𝑖∗= 0 in 𝑡 ∈ [𝑡
0, 𝑡
𝑓] on 𝛤
1𝑢and 𝛤
𝑑𝑒𝑠𝑖𝑔𝑛 （21）𝜃
∗= 0 in 𝑡 ∈ [𝑡
0, 𝑡
𝑓] on 𝛤
1𝜃 （22）𝑡
𝑖∗= 𝑡 ̂ = {𝑢
𝑖∗ 𝑗𝑢
𝑖∗+ 𝑝
∗𝛿
𝑖𝑗+ 𝜈𝑢
𝑖,𝑗∗}𝑛
𝑗in 𝑡 ∈ [𝑡
0, 𝑡
𝑓] on 𝛤
2𝑢 （23
）𝑏
∗= 𝑏 ̂ = {𝜃
∗ ∗𝑢
𝑖+ 𝛼𝜃
,𝑖∗}𝑛
𝑖in 𝑡 ∈ [𝑡
0, 𝑡
𝑓] on 𝛤
2𝜃 （24
） さらに，設計境界の座標に対するLagrange
関数の勾配𝐺
𝑖は，式（25
）から導出される．
∫ ∫ 𝑡
𝛤 𝑖∗𝛿𝑢
𝑖 𝑡𝑓𝑡0
𝑑𝛤𝑑𝑡 + ∫ ∫ 𝑏
𝑡𝑡𝑓 𝛤 ∗𝛿𝜃
0
𝑑𝛤𝑑𝑡
= ∫ ∫
𝛤𝑡
𝑖∗𝛿𝑢
𝑖1𝑢+𝛤2𝑢 𝑡𝑓
𝑡0
𝑑𝛤𝑑𝑡 + ∫ ∫
𝛤𝑡
𝑖∗𝑢
𝑖,𝑗𝛿𝑥
𝑗𝑑𝑒𝑠𝑖𝑔𝑛 𝑡𝑓
𝑡0
𝑑𝛤𝑑𝑡 + ∫ ∫
𝛤𝑏
∗𝛿𝜃
1𝜃+𝛤2𝜃 𝑡𝑓
𝑡0
𝑑𝛤𝑑𝑡 + ∫ ∫
𝛤𝑏
∗𝜃
,𝑖𝛿𝑥
𝑖𝑑𝑒𝑠𝑖𝑔𝑛 𝑡𝑓
𝑡0
𝑑𝛤𝑑𝑡
= ∫ ∫
𝛤1𝑢+𝛤2𝑢𝑡
𝑖∗𝛿𝑢
𝑖 𝑡𝑓𝑡0
𝑑𝛤𝑑𝑡 + ∫ ∫
𝛤𝑏
∗𝛿𝜃
1𝜃+𝛤2𝜃 𝑡𝑓
𝑡0
𝑑𝛤𝑑𝑡 + ∫
𝛤(𝐺
𝑗𝛿𝑥
𝑗+ 𝐺
𝑖𝛿𝑥
𝑖)
𝑑𝑒𝑠𝑖𝑔𝑛
𝑑𝛤
（25）勾配 𝐺𝑖
は，線形弾性体の外力として考え，変位値を修正勾配 𝐺
𝑖∗として用いる手法，力法を適用することで変
更される．式（26）に示すように，修正された勾配 𝐺𝑖∗
を用いて，目標境界上の座標は更新される．ここで，𝑙
は形状更新ステップ回数，𝜂はステップ幅である．𝑥
𝑖(𝑙+1)= 𝑥
𝑖(𝑙)− 𝜂𝐺
𝑖∗(𝑙) （26）2・5 解析手順
本研究における形状最適化は，以下のステップを繰り返すことによって行われる．
Step 1
目標形状を与える．Step 2
状態方程式を用いて，時間 𝑡 = 𝑡0から 𝑡 = 𝑡
𝑓の方向へ，流速分布 𝑢
𝑖，圧力分布 𝑝 ，温度分布 𝜃 を解析す
る．このとき，目的汎関数の計算のために目標形状における温度分布
𝜃
を記録する．Step 3
初期形状を与える．Step 4
状態方程式を用いて，時間 𝑡 = 𝑡0から 𝑡 = 𝑡
𝑓の方向へ，流速分布 𝑢
𝑖，圧力分布 𝑝 ，温度分布 𝜃 を解析す
る．
Step 5
目的汎関数の収束判定を行う．式（27
）に示す収束判定条件式を満たした場合，解析を終了する．|𝐽
(𝑘)− 𝐽
(𝑘−1)| 𝐽 ⁄
(𝑘)< 𝜀
（27
）Step 6
上記のStep 4
で得られた流速分布 𝑢𝑖，圧力分布 𝑝 ，温度分布 𝜃 および随伴方程式を用いて，時間 𝑡 = 𝑡
𝑓から 𝑡 = 𝑡0
の方向へ，随伴流速分布 𝑢
𝑖∗，随伴圧力分布 𝑝
∗，随伴温度分布 𝜃
∗を解析する．
Step 7
これらの結果を用いて，形状勾配関数 𝐺𝑖を導出する．
Step 8
力法に基づいて形状更新を行い，Step 4に戻る．3. 解析例
3・1 解析モデル及び解析条件
放熱量規定問題の解析モデルとして，図
1
のようなヒーターを設定した．本論文の序論で述べたように，冷却 効率の向上を目的とするヒーターの形状の一例として，赤線に示す長方形を初期形状とし，赤色破線に示す三角 形を目標形状として設定した．初期形状および目標形状において，Robin
境界𝛤
3𝜃における目的汎関数を計算す る．また，温度既知境界 𝛤1𝜃を設計境界とし，すべての境界に対してノンスリップ境界条件 𝑢
𝑖= 0 を与えた．さ
らに，
Neumann
境界 𝛤2𝜃においては断熱境界として 𝑞 = 0 を与えた．図 1
中の領域の寸法は，𝐿 = 1 [m], 𝑙 = 0.05
[m]である．表 1
に，図1
の解析モデルにおける計算条件を示す．Table1 A Computational condition
Initial model Target model
Number of nodes 2769 2599
Number of elements 5284 4923
Time steps 𝑡
𝑓[s] 3000
Time increment Δ 𝑡 [s] 3
Convergence criterion 10
−4Kinematic viscosity 𝜈 [m
2/s] 1.16 × 10
−2Density 𝜌 [kg/m
3] 860
Thermal expansion coefficient 𝛽 [1/K] 7.5 × 10
−4Gravity 𝑔 [m/s
2] 9.8
Reference temperature 𝜃
0[K] 293.15 Temperature diffusivity 𝛼 [m
2/s] 8.48 × 10
−5Fig. 1 A computational model.
3・2 解析結果
図
2
に目標形状における有限要素メッシュ，温度分布，流速・圧力分布，図3
に同定形状における有限要素メ ッシュ，温度分布，流速・圧力分布，図4
に目的汎関数の収束履歴を示す．形状更新の際には，FreeFem++
上の関数「
adaptmesh
」を用いてリメッシュを行った．図
2, 3
を比較すると，得られた同定形状の温度分布と目標形状における温度分布は一致していることが確認で きる．また，図5
より，初期値を100 [%]とした目的汎関数が十分に収束しており，目標とした放熱量に同定して
いることが確認できる．同定形状において，三角形の下部が膨らむような形状になっているが，これは領域内で スムーズに対流が起こりやすくなるように変化したからであると考えられる．4. 結論
本研究では，ある目標形状において放熱量を規定し，形状を同定する問題について取り組んだ．また，得ら れた結果から，随伴変数法および有限要素法に基づいた定式化，作成したプログラムの妥当性を示した．今後の 課題として，
Navier-Stokes
方程式における浮力の項に対してBoussinesq
近似を適用しない，および温度による密 度変化を考慮した問題を予定している．謝 辞
本論文を執筆するにあたり，理化学研究所 放射光科学研究センターの近藤力様の助言を受けた．記して深く謝 意を表する．
文 献
