論 文

熱対流場における熱源境界形状の最適化に関する数値実験- アダプティブリメッシングにおけるパラメータ設定が最適化後の形状に与える影響 -和田 海斗 1・倉橋 貴彦 2・河田 剛毅 3
1 長岡技術科学大学 工学部 工学課程 機械工学分野 (現:倉敷機械株式会社)（Mechanical Engineering, Bachelor’s Program, School of Engineering,
Nagaoka University of Technology）
2 長岡技術科学大学 技学研究院 機械系（Department of Mechanical Engineering, Institute of GIGAKU, Nagaoka University of Technology）
3 機械工学科(Department of Mechanical Engineering, Natural Institute of Technology, Nagaoka College)

Numerical experiments on optimization of heat source boundary shapein thermal convection field- Effect of parameter setting in adaptive remeshing on optimized shape Kaito WADA1, Takahiko KURAHASHI2 and Yoshitaka KAWADA3
Abstract
In this study, we perform shape optimization of a heat source in a thermal convection field considering
perimeter constraint condition, and investigate effect of parameter setting in adaptive remeshing on
optimized shape. The purpose of this study is to obtain an optimal shape of inner heat source so as tomaximize the radiation amount on outer boundary. The calculation for the optimization is carried out basedon the adjoint variable and the finite element methods, and the perimeter constraint condition is consideredas the constraint condition. We especially focus on the effect of parameter setting in adaptive remeshing on
optimized shape, and some numerical results are shown in this paper. The program code for the shape
optimization analysis is constructed by "FreeFem++".
Key Words : shape identification, finite element method, adjoint variable method, thermal convection
field, adaptive remeshing

１．はじめにの発熱量は増加している．機器の温度上昇は，動作効率の低下や誤作動につながることから冷却機構を

近年，電子機器の小型化に伴い，単位体積当たり用いて冷却することが必要となる．電子機器の冷却
21

方法には装置を用いた強制対流による冷却方法があ式，熱伝達方程式を考える．Navier-Stokes 方程式，る．しかし，電子機器の小型化において冷却装置を連続の式は非圧縮性とする．これらの支配方程式は設置することは望ましくない．また，装置にも電力式(3)～(5)のように表される．

を使用することで機器全体の消費電力が高くなる

1)．

そのため，装置を必要としない冷却方法として自然
対流による冷却方法があげられる．
自然対流を使用している電子機器の例としてクライストロン電源があげられる．近藤ら 2)が提案した

1
𝑢𝑢1̇ + 𝑢𝑢𝑗𝑗 𝑢𝑢𝑖𝑖,𝑗𝑗 + 𝑝𝑝,𝑖𝑖 − 𝜈𝜈𝑢𝑢𝑖𝑖,𝑗𝑗𝑗𝑗 − 𝛽𝛽𝛽𝛽(𝜃𝜃 − 𝜃𝜃0 )𝑒𝑒𝑖𝑖
𝜌𝜌= 0 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Ω
𝑢𝑢𝑖𝑖,𝑖𝑖 = 0 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Ω

クライストロン電源は絶縁油で満たした一つのタン
ク内に高電圧部品を収納しており，メンテナンス性𝜃𝜃̇ + 𝑢𝑢𝑖𝑖 𝜃𝜃,𝑖𝑖 − 𝛼𝛼𝜃𝜃,𝑖𝑖𝑖𝑖 = 0 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Ωやコンパクト性が求められるため，装置を必要としない自然対流での冷却を用いている．

温度，温度拡散率を表している．また，𝑒𝑒1 = 𝑒𝑒3 =

した形状最適化問題 4)ではリメッシュの際にメッシ

0, 𝑒𝑒2 = 1である．

ュのつぶれにより収束判定よりも前に解析が終了し
ている．

２．３

そこで，本論文では周長制約を考慮することで熱式(6)～(12)のように定義する．

化問題を行う．その中で，リメッシュにおいてパラメータ設定を変更することで最適形状にもたらす影
響について示す．数値シミュレーションはプログラ
ミング言語 FreeFem++ 5),6)を用いて行う．

２．放熱量最大化問題の定式化
目的汎関数に定義する．この問題では，ある任意の境界におけ(6)

𝑢𝑢𝑖𝑖 = 𝑢𝑢�𝑖𝑖 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γ1u and Γdesign(7)𝜃𝜃 = 𝜃𝜃� at 𝑡𝑡 = 𝑡𝑡0 in Ω

る放熱量が最大となることを考える．
(1)𝜃𝜃 = 𝜃𝜃� in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γ1θ

𝑏𝑏 = 𝑏𝑏� = −𝜃𝜃,𝑖𝑖 𝑛𝑛𝑖𝑖 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γ2θ

(2)

𝑞𝑞 = ℎ��𝜃𝜃 − 𝜃𝜃�𝑓𝑓 � = −𝛼𝛼𝜃𝜃,𝑖𝑖 𝑛𝑛𝑖𝑖
in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γ3θ

である．𝑄𝑄𝑖𝑖 は重み定数，𝑞𝑞は各計算ステップにおける放熱量， ℎ� は熱伝達係数，𝜃𝜃は温度，𝜃𝜃�𝑓𝑓 は外気温度を表している．𝑡𝑡0 は数値解析における初期時刻，
𝑡𝑡𝑓𝑓 は最終時刻，Γは数値計算境界を表している．また，重み定数𝑄𝑄𝑖𝑖 は放熱量を最大にする境界で 1 とし，他の境界では 0 として考える．
２．２

𝑢𝑢𝑖𝑖 = 𝑢𝑢�𝑖𝑖 at 𝑡𝑡 = 𝑡𝑡0 in Ω

1
𝑡𝑡𝑖𝑖 = 𝑡𝑡𝑖𝑖̂ = �− 𝑝𝑝𝛿𝛿𝑖𝑖𝑖𝑖 + 𝑢𝑢𝑖𝑖,𝑗𝑗 � 𝑛𝑛𝑗𝑗
𝜌𝜌
in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γ2u

放熱量最大化問題の目的汎関数𝐽𝐽を，式(1)のよう

𝑞𝑞 = ℎ��𝜃𝜃 − 𝜃𝜃�𝑓𝑓 �

境界条件および初期条件

非定常自然対流場における境界条件，初期条件を伝達部が一定の場合の放熱量増加のための形状最適

ここで，

(5)

密度，動粘性係数，体積膨張率，重力加速度，基準

形状最適化が行われてきた．特に，面積制約を考慮

1 𝑡𝑡𝑓𝑓
𝐽𝐽 = − � �𝑄𝑄𝑖𝑖 𝑞𝑞 2 𝑑𝑑Γdt
2 𝑡𝑡0 Γ(4)

ここで，𝑢𝑢𝑖𝑖 , 𝑝𝑝, 𝜌𝜌, 𝜈𝜈, 𝛽𝛽, 𝑔𝑔, 𝜃𝜃0 , 𝛼𝛼はそれぞれ流速，圧力，

先行研究として青木ら 3)による熱対流場における

２．１

(3)

(8)
(9)
(10)
(11)

(12)

� はそれぞれクロネッカーのデルタ，
ここで，𝛿𝛿𝑖𝑖𝑖𝑖 , 𝑛𝑛𝑗𝑗 , (∙)

境界における単位法線ベクトル，既知境界関数を表している．また，Γdesign , Γ1u , Γ2u , Γ1θ , Γ2θ , Γ3θ はそれぞれ 設 計 境 界 ， 流 れ 場 に お け る Dirichlet 境 界 ，Neumann 境界， 温 度場 に おける Dirichlet 境界 ，

支配方程式

Neumann 境界，Robin 境界を表している．

非定常自然対流場における支配方程式として，
Boussinesq 近似された Navier-Stokes 方程式，連続の
22

２．４

tf

随伴方程式および形状勾配関数の導出

任意の随伴変数𝑢𝑢𝑖𝑖∗ , 𝑝𝑝 ∗ , 𝜃𝜃 ∗ を用いることで，目的汎t0tf

関数𝐽𝐽は Lagrange 関数𝐽𝐽 に拡張され、式(13)で表され
∗

る．

=� �
t0
tf

𝑡𝑡𝑓𝑓
1
𝐽𝐽∗ = 𝐽𝐽 + � � 𝑢𝑢𝑖𝑖∗ �𝑢𝑢𝑙𝑙̇ + 𝑢𝑢𝑙𝑙 𝑢𝑢𝑖𝑖,𝑗𝑗 + 𝑝𝑝,𝑖𝑖 𝜈𝜈𝑢𝑢𝑖𝑖,𝑗𝑗𝑗𝑗
𝜌𝜌
𝑡𝑡0 Ω
𝑡𝑡𝑓𝑓

𝑡𝑡0
tf

Ω

+ � � 𝜃𝜃 ∗ �𝜃𝜃̇ + 𝑢𝑢𝑖𝑖 𝜃𝜃,𝑖𝑖 − 𝛼𝛼𝜃𝜃,𝑖𝑖𝑖𝑖 �dΩdt𝛤𝛤𝛤𝛤1𝑢𝑢+2𝑢𝑢+� �t0+ 𝛽𝛽𝛽𝛽(𝜃𝜃 − 𝜃𝜃0 )𝑒𝑒𝑖𝑖 � dΩdt

+ � � 𝑝𝑝∗ 𝑢𝑢𝑖𝑖,𝑖𝑖 𝑑𝑑Ω𝑑𝑑𝑑𝑑tf

� �𝑡𝑡𝑖𝑖∗ 𝛿𝛿𝑢𝑢𝑖𝑖 dΓdt + � �𝑏𝑏 ∗ 𝛿𝛿𝛿𝛿dΓdt

𝛤𝛤1𝑢𝑢+2𝑢𝑢
tf
t0
𝑏𝑏

れる．

(13)

�2 �𝜃𝜃 − 𝜃𝜃(𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡𝑡) �𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑 = 0+ � �𝑄𝑄ℎt0

Γin 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Ω

(15)

(𝐹𝐹 𝑏𝑏 + 𝐻𝐻𝑏𝑏 )dΓdt + 𝐺𝐺

(23)

(24)

𝐹𝐹 𝑑𝑑 = 𝑡𝑡𝑖𝑖∗ 𝑢𝑢𝑖𝑖,𝑗𝑗 𝛿𝛿𝑥𝑥𝑗𝑗

(25)

𝐻𝐻𝑏𝑏 = 𝑏𝑏 ∗ 𝛿𝛿𝛿𝛿(26)

𝐻𝐻𝑑𝑑 = 𝑏𝑏 ∗ 𝜃𝜃,𝑖𝑖 𝛿𝛿𝑥𝑥𝑖𝑖(27)

𝛤𝛤𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑∗�𝜈𝜈𝑢𝑢𝑗𝑗,𝑖𝑖𝑢𝑢𝑗𝑗,𝑖𝑖 + 𝛼𝛼𝜃𝜃,𝑖𝑖∗ �dΓ

(28)

勾配𝐺𝐺𝑖𝑖 は，線形弾性体の外力として考え，得られる(16)

反力を変位場𝑉𝑉𝑖𝑖 として用いる手法，力法を適用する．

線形弾性体の支配方程式を式(29)から式(31)，境界(17)

𝜃𝜃 ∗ = 0 at 𝑡𝑡 = 𝑡𝑡𝑓𝑓 in Ω(18)u∗i = 0 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γ1u and Γdesign(19)𝜃𝜃 ∗ = 0 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γ1θ

(20)

in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γ2u(21)

𝑏𝑏 ∗ = 𝑡𝑡�𝚤𝚤∗ = �𝜃𝜃 ∗ 𝑢𝑢𝑖𝑖 + 𝛼𝛼𝜃𝜃,𝑖𝑖∗ �𝑛𝑛𝑖𝑖
in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γ2θ𝛤𝛤𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑

𝐻𝐻𝑑𝑑 dΓdt

𝐹𝐹 𝑏𝑏 = 𝑡𝑡𝑖𝑖∗ 𝛿𝛿𝑢𝑢𝑖𝑖

𝐺𝐺 = �

𝑢𝑢𝑖𝑖∗ = 0 at 𝑡𝑡 = 𝑡𝑡𝑓𝑓 in Ω∗∗
𝑡𝑡𝑖𝑖∗ = 𝑡𝑡�𝚤𝚤 = �𝑢𝑢𝑗𝑗 𝑢𝑢𝑖𝑖∗ + 𝑝𝑝∗ 𝛿𝛿𝑖𝑖𝑖𝑖 + 𝜈𝜈𝑢𝑢𝑖𝑖,𝑗𝑗
�𝑛𝑛𝑗𝑗

𝑡𝑡0

(14)

∗
−𝜃𝜃 ∗̇ − (𝜃𝜃 ∗ 𝑢𝑢𝑖𝑖 ),𝑖𝑖 − 𝛼𝛼𝜃𝜃,𝑖𝑖𝑖𝑖+ 𝛽𝛽𝛽𝛽𝑢𝑢𝑖𝑖∗ 𝑒𝑒𝑖𝑖
tf

𝐻𝐻𝑏𝑏 dΓdt + � �

𝐹𝐹 𝑑𝑑 dΓdt

𝐹𝐹 𝑏𝑏 ，𝐹𝐹 𝑑𝑑 ，𝐻𝐻 ，𝐻𝐻 ，𝐺𝐺はそれぞれ式(24)～(28)で表さ随伴方程式および随伴条件は次のように決定される．

1 ∗− 𝑢𝑢𝑖𝑖,𝑖𝑖= 0 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Ωρ

𝑡𝑡0 𝛤𝛤𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑𝑑
𝑡𝑡𝑓𝑓

𝛤𝛤1𝑢𝑢+2𝑢𝑢
𝑑𝑑

Lagrange 関数の停留状態（𝛿𝛿𝐽𝐽 = 0）を考慮すると，= 0 in t ∈ �t 0 . t𝑓𝑓 � on Ω

𝐹𝐹 𝑏𝑏 dΓdt + � �

=� �

t0 Ω
𝑢𝑢𝑖𝑖∗ は随伴流速，𝑝𝑝∗ は随伴圧力，𝜃𝜃 ∗ 随伴温度である．
∗∗−𝑢𝑢𝑙𝑙̇ ∗ + 𝑢𝑢𝑗𝑗,𝑖𝑖 𝑢𝑢𝑗𝑗∗ − �𝑢𝑢𝑗𝑗 𝑢𝑢𝑖𝑖∗ �,𝑗𝑗 − 𝑝𝑝,𝑖𝑖∗ − 𝜈𝜈𝑢𝑢𝑖𝑖,𝑗𝑗𝑗𝑗+ 𝜃𝜃 ∗ 𝜃𝜃,𝑖𝑖

t0 𝛤𝛤
𝑡𝑡𝑓𝑓

条件を式(32)と式(33)に示す．
𝜎𝜎𝑖𝑖𝑖𝑖,𝑗𝑗 = 0 in Ω𝜀𝜀𝑖𝑖𝑖𝑖 =1
�𝑉𝑉 + 𝑉𝑉𝑗𝑗,𝑖𝑖 � in Ω
2 𝑖𝑖,𝑗𝑗(29)(30)

𝜎𝜎𝑖𝑖𝑖𝑖 = 𝐷𝐷𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝜀𝜀𝑘𝑘𝑘𝑘 in Ω

(31)

𝑉𝑉𝑖𝑖 = 0 in Γ ∖ Γdesign(32)

𝑇𝑇𝑖𝑖 = 𝐺𝐺𝑖𝑖 on Γdesign(33)

ここで，力法に基づき提案されている変形法 7)によ(22)

り周長を一定とするため線形弾性体の境界に一様な外力として定数λを加えた場合の解析を行う．この

設計境界の座標に対する Lagrange 関数の勾配𝐺𝐺𝑖𝑖 は，

解析で得られる変位場をΔ𝑉𝑉𝑖𝑖 として変位場𝑉𝑉𝑖𝑖 と合わ

式(23)から導出される．

せることで周長を一定にした形状の変形ができる．

その場合の境界条件を式(34)と式(35)に示す．また，
感度𝐺𝐺𝑖𝑖 と定数λによって得られる変位場の模式図を図-1 から図-3 に示す．

23

Δ𝑉𝑉𝑖𝑖 = 0 in Γ ∖ Γdesign(34)

𝑇𝑇𝑖𝑖 = 𝜆𝜆 on Γdesign(35)

𝜎𝜎𝑖𝑖𝑖𝑖 は応力，𝜀𝜀𝑖𝑖𝑖𝑖 はひずみ，𝐷𝐷𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 は弾性係数行列，𝑇𝑇𝑖𝑖 は表面力を表す．無次元のヤング率𝐸𝐸 = 1，ポアソン比を 0 とし，弾性係数行列を式(36)に示す．
1 0 0[𝐷𝐷] = �0 1 0�
(36)
0 0 1
変位場𝑉𝑉𝑖𝑖 による周長の変化量𝑆𝑆𝑉𝑉 と変位場Δ𝑉𝑉𝑖𝑖 による

図-1

周長の変化量𝑆𝑆Δ𝑉𝑉 を式(37)と式(38)より求める．また，

感度𝐺𝐺𝑖𝑖 による形状更新

このとき定数λは未知の値であるため1と仮定したときの変位場Δ𝑉𝑉𝑖𝑖 (𝜆𝜆 = 1)とし，その模式図を図-4 に示す．

𝑆𝑆𝑉𝑉 = � 𝑑𝑑Γ − �Γ𝑉𝑉Γdesign

𝑆𝑆Δ𝑉𝑉(𝜆𝜆=1) = �ΓΔ𝑉𝑉(𝜆𝜆=1)𝑑𝑑Γ − �

Γdesign定数λは式(39)より求める．
𝑆𝑆𝑉𝑉 + 𝜆𝜆𝑆𝑆Δ𝑉𝑉(𝜆𝜆=1) = �𝑑𝑑Γ

Γ0design(37)𝑑𝑑Γ

𝑑𝑑Γ − �Γ𝑙𝑙design

(38)

𝑑𝑑Γ

(39)

図-2

定数λが正の値の場合の形状更新

図-3

定数λが負の値の場合の形状更新

𝑙𝑙0Γdesign，Γdesign
はそれぞれ形状更新𝑙𝑙回目の設計境界，

初期形状の設計境界を表す．ここで右辺は周長制約を考慮しているため0となる．したがって式(39)は二つの周長変化量の和が0となることを表す．二つの周長変化量の和が0となる定数λは式(40)より求める．
𝜆𝜆 = −∫Γ𝑙𝑙

design𝑑𝑑Γ − ∫Γ0design𝑆𝑆Δ𝑉𝑉(𝜆𝜆=1)𝑑𝑑Γ + 𝑆𝑆𝑉𝑉

(40)

変位場𝑉𝑉𝑖𝑖 とΔ𝑉𝑉𝑖𝑖 を用いて式(41)に示すように設計境界上の座標が更新される．ここで，𝑙𝑙は形状更新ステ
ップ数，𝜂𝜂はステップ幅である．
(𝑙𝑙+1)

𝑥𝑥𝑖𝑖(𝑙𝑙)= 𝑥𝑥𝑖𝑖 − 𝜂𝜂�𝑉𝑉𝑖𝑖𝑙𝑙 + Δ𝑉𝑉𝑖𝑖𝑙𝑙 � in Ω

(41)

図-4

24

定数λ = 1の場合の形状更新

２．５の動作を「adaptmesh」のコマンドにより自動的に

計算の流れ

行う．

本研究における形状最適化は，以下のステップを繰り返すことによって行われる．
Step1．初期形状を与える．
Step2．状態方程式を用いて，時間𝑡𝑡 = 𝑡𝑡0 から𝑡𝑡 = 𝑡𝑡𝑓𝑓の方向へ，流速分布𝑢𝑢𝑖𝑖 ，圧力分布𝑝𝑝，温度分
布𝜃𝜃を解析する．

Step3．目的汎関数の収束判定を行う．式(31)に示す
収束判定条件式を満たす場合，解析を終了
して，大きい場合は次のステップへ進む．
�

𝐽𝐽

(𝑙𝑙)

(𝑙𝑙−1)

− 𝐽𝐽
𝐽𝐽(1)

図-5

三角形メッシュが生成される例

図-6

三角形メッシュが作られない例

� < 𝜖𝜖

Step4．Step2 で得られた流速分布𝑢𝑢𝑖𝑖 ，圧力分布𝑝𝑝，

温度分布𝜃𝜃および随伴方程式を用いて時間
𝑡𝑡 = 𝑡𝑡𝑓𝑓 から𝑡𝑡 = 𝑡𝑡0 の方向へ，随伴流速分布𝑢𝑢𝑖𝑖∗ ，
随伴圧力分布𝑝𝑝∗ ，随伴温度分布𝜃𝜃 ∗ を解析する．

Step5．これらの結果を用いて勾配𝐺𝐺𝑖𝑖 を解析する．

Step6．勾配𝐺𝐺𝑖𝑖 に力法を適用し，修正された勾配𝐺𝐺𝑖𝑖∗

を計算し，形状更新を行い再度 Step2 に戻る．

４．数値解析例
３．アダプティブリメッシュ法
４．１

解析条件

アダプティブリメッシュ法は数値解析の結果をも放熱量最大化問題の解析モデルとして，図-6 のとに誤差の推定を行い，得られた結果から節点距離

ような円形を初期形状としたヒーターを設定した．

を決定して新たなメッシュを構成する手法である

8)．

この手法を用いることで値が急激に変化する領域のを計算する．また，温度既知境界𝛤𝛤1θ を設計境界
𝛤𝛤design として𝜃𝜃 = 298.15[K]を与え，すべての境界にみ有限要素サイズを小さくすることが可能になるた

対してノンスリップ境界条件𝑢𝑢𝑖𝑖 = 0[m/s]を与えた．

め領域全体の要素数を大きく増やすことなく精度のよい計算が可能となる．本解析では解析結果から節

また，Neumann 境界𝛤𝛤2𝜃𝜃 においては断熱境界として

𝑞𝑞 = 0[W/m2]を与えた．図-6 中の領域の寸法は，点位置を決定する．アダプティブリメッシュによる

𝐿𝐿 = 1[m],𝑟𝑟 = 0.2[m]である．表-1 に，図-6 の解析モ

節点の位置は式(42)を満たすように決定される．
sup|𝑢𝑢 − 𝑢𝑢ℎ | ≤ 𝑒𝑒

初期形状において Robin 境界𝛤𝛤3𝜃𝜃 における目的汎関数

デルにおける計算条件を示す．FreeFem++における(42)

アダプティブリメッシュは「adaptmesh」のコマンドにより自動的に行われる．その中でオプションを

𝑢𝑢 ， 𝑢𝑢ℎ ， 𝑒𝑒はそれぞれ厳密解，有限要素近似解，補設定することでリメッシュの制御を行うことができ

間誤差の上限を表す．式(42)は補間誤差が上限以下

る．三角形メッシュの辺の最小長さは「hmin」のオになることを示す．補間誤差の上限𝑒𝑒は初期設定値

プション，最大節点数の上限を「nbvx」のオプショ

では𝑒𝑒 = 0.01となっており，FreeFem++のオプションにより設定することができる．三角形メッシュのンにより変更することができる．決定した節点から最小長さは領域の大きさにより初期設定値が自動的

Delaunay 三角形分割 9)により新しい三角形メッシュに決定される．図-6 の解析モデルでは，初期設定の生成を行う．Delaunay 三角形分割は領域内の節点値は約0.016[m]となる．三角形メッシュの辺の最小を三点取り出して作られる外接円の中に他の節点が長さと最大節点数の上限を表-2 に示す組み合わせ

含まれていないとき，取り出した三点で三角形を作で解析を行った．

る手法である．三角形が作れる例と作ることができない例を図-5，図-6 に示す．FreeFem++ではこれら
25

５． 結論
本論文では，放熱量の最大化を目的として周長一定の制約のもと，熱対流場において発熱源の形状最適化を行った．その中で，FreeFem++のコマンドである「adaptmesh」によるアダプティブリメッシュのパラメータ設定を変更することによって最適形状に及ぼす影響について示した．また，今回の結果から設計境界が重なる可能性が確認されたため，初期形状を変えて解析を行うことを今後の課題とする．

図-7

謝辞：本論文を執筆するにあたり, 長岡技術科学

数値解析モデル大学機械創造工学専攻 2020 年度修了の青木崇氏と
表-1

節点数
要素数
時間ステップ𝑡𝑡𝑓𝑓 [s]
時間刻み幅Δ𝑡𝑡 [s]密度𝜌𝜌 [kg/m3]
体積膨張率𝛽𝛽 [1/K]
重力加速度𝑔𝑔 [m/s2]
基準温度𝜃𝜃0 [K]
温度拡散率α [m2/s]
収束判定値𝜖𝜖 [-]
表-2
Case1
Case2
Case3

４．２

解析結果

岐阜工業高等専門学校の片峯英次教授の助言を受け

数値解析条件

初期形状
2524
4786
3000
3
860
7.5 × 10−4
9.8
293.158.48 × 10−5
1.0 × 10−5

た．記して深く謝意を表する．

各解析条件

最小長さ [m]8.0 × 10−3
1.6 × 10−3
1.6 × 10−3

最大節点数[-]
9000
9000
4000

図-8

初期形状の有限要素メッシュ

図-9

初期形状の温度-流速分布

初期形状での有限要素メッシュを図-8，温度・流速分布を図-9 に，解析結果を図-10 から図-12 に示す．また，本解析時間を図-13 に，最終形状での収束判定値を表-3 に示す．表-3 よりリメッシュの際の三角形メッシュのサイズを小さくすることで収束判定値に近づけることができた．また，図-13 に示す解析時間の履歴から三角形メッシュを小さくする
ことで起こる解析時間の増加は最大節点数を抑制することで改善することができた．今回得られた最終形状である図-10 から図-12 の形状は上部分が細く，下部分が太い形状であったが，図-12 では最終形状の上部分の境界が重なる傾向が見られた．三角形メッシュのサイズを変えることで境界付近の解析精度が変わるため，最終形状に差異が生じたと考える．

26

図-13
表-3
図-10

各最終形状での収束判定値

Case1 の温度-流速分布
Case1
Case2
Case3

参考文献
1)

各解析時間

収束判定値[-]
1.09 × 10−3
0.927 × 10−3
0.847 × 10−3

吉川実,“電子機器における省エネルギー冷却技術の研究開発”,NEC 技術,Vol.63,No4,pp.55-58,2010.

2)

近藤力,他“クライストロン電源における絶縁油の冷却効率の向上”,第 4 回加速器学会,TP42,2008,和光.

3)

青木崇,他“随伴変数法および有限要素法に基づく熱対流場における熱源境界形状の最適化”,長岡工業高等専門学校研究紀要,Vol.56,pp.17-23,2020.

4)

青木崇,他“熱対流場において粘性効果が形状最適化の結果に与える影響に関する考察”,日本機械学会,計図-11

力スクウェア研究報告書,Report,No.3-06,2020.

Case2の温度-流速分布

5)

大 塚 厚 二,他“有 限 要 素 法 で 学 ぶ 現 象 と 数 理FreeFem++数理思考プログラミング-”,共立出版株式会社,2014.

6)

Hecht F,“New development in FreeFem++”,Journal of
Numerical Mathematics,20(3-4) PP.251-265,2012.

7)

能井宏弥,他“力法の原理を用いたシェル構造の形状
最 適 化 ”, 日 本 建 築 学 会 構 造 系 論 文集,Vol.72,No.616,pp.121-126,2007.

8)

手塚明,他“アダプティブ有限要素法”,丸善株式会
社,2003.

9)

谷口健男,“FEM のための要素自動分割 デローニー三角分割法の利用”森北出版社株式会社,1992.
(2022. 8. 8 受付)

図-12

Case3 の温度-流速分布27