論 文

雪堤の変形抑制を目的とした
道路付帯構造物の形状最適化
広瀬 智史1・倉橋 貴彦2・芝崎 智貴3・杉原 幸信2・上村 靖司2・河田 剛毅4
1 長岡技術科学大学大学院 工学研究科 工学専攻 機械工学分野(現:株式会社メイテック)（Mechanical Engineering, Master’s Program, Graduate School of Engineering,
Nagaoka University of Technology）
2 長岡技術科学大学 技学研究院 機械系（Department of Mechanical Engineering, Institute of GIGAKU, Nagaoka University of Technology）
3 長岡技術科学大学大学院 工学研究科 工学専攻 機械工学分野(現:株式会社竹内製作所)（Mechanical Engineering, Master’s Program, Graduate School of Engineering,
Nagaoka University of Technology）
4 機械工学科（Department of Mechanical Engineering, National Institute of Technology, Nagaoka College）

Shape optimization of road incidental structures for deformation control of snowbank
Satoshi HIROSE1, Takahiko KURAHASHI2, Tomoki SIBASAKI3,
Yukinobu SUGIHARA2, Seiji KAMIMURA2 and Yoshitaka KAWADA4
Abstract
In this study, we present some numerical results of the shape optimization analysis of road incidental
structures for deformation control of snowbank. The performance function was defined by the energy due
to thermal deformation, and the shape optimization analysis was formulated based on the adjoint variable
and the finite element methods. The shape update was carried out by the gradient of the Lagrange function,i.e., the performance function extended by the adjoint variable and the governing equations, with respect tothe design variable and the traction method was applied to control the gradient oscillation. The shape of theroad incidental structure was updated by the gradient modified by the traction method so as to minimize
the performance function. The Python program was used to estimate the coefficient of thermal expansion,and the FreeFem program was used for the shape optimization analysis.
Key Words : Finite Element Method, FreeFem, Thermal elastic fields, Snow bank, Experiment for melting
phenomenon of snow bank

１．

はじめに

狭い間隔で等圧線が形成される西高東低型の気圧配置になる．冷え込んだシベリアから日本に流れ込む

日本の気候は冬期に着目すると，西に冷たいシベ

季節風は，日本海を流れる暖流の対馬海流によってリア高気圧，東に低気圧があり，日本近辺を南北に水蒸気が補給されるとともに，暖められ大気が不安7定となり積雲を発生させる．この積雲が，本州中央

𝑡𝑡𝑓𝑓の山々にぶつかることで，日本海側の多くの地域に

𝐽𝐽 = � �� �𝐶𝐶𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝛼𝛼𝑘𝑘𝑘𝑘 �𝜙𝜙 − 𝜙𝜙𝑓𝑓 �𝜀𝜀𝑖𝑖𝑖𝑖 (𝑢𝑢𝑖𝑖 ) + 𝑓𝑓𝑖𝑖 𝑢𝑢𝑖𝑖 �𝑑𝑑Ω
𝑡𝑡0

おいて雪を降らせる．このような気候のことを日本
海側気候あるいは日本海岸式気候と呼び，世界有数Ω+ � 𝑄𝑄𝑝𝑝 𝑃𝑃𝑖𝑖 𝑢𝑢𝑖𝑖 𝑑𝑑Γ� 𝑑𝑑𝑑𝑑
Γ1
𝜀𝜀𝑖𝑖𝑖𝑖 (𝑢𝑢𝑖𝑖 ) = (𝑢𝑢𝑖𝑖,𝑗𝑗 + 𝑢𝑢𝑗𝑗,𝑖𝑖 )2の豪雪地帯に分類される．その豪雪地帯では，市街地や高速道路において交通流の安定した確保のために，路面上の積雪の除去のために除雪トラックや除

(1)
(2)

ここで， 𝑄𝑄𝑝𝑝 は重み定数，𝜙𝜙は温度， 𝜙𝜙𝑓𝑓 は時刻𝑡𝑡 = 0
における領域内の温度， 𝐶𝐶𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 は剛性テンソル，𝛼𝛼𝑘𝑘𝑘𝑘

雪ドーザを用いた新雪除雪や，路肩に寄せた雪に対する車線幅確保のためのロータリー除雪車を用いたは熱膨張係数テンソル，𝑓𝑓𝑖𝑖 は𝑖𝑖方向の体積力，𝑢𝑢𝑖𝑖 は各拡幅除雪が行われ，中央分離帯および路肩に雪堤と計算ステップにおける𝑖𝑖方向の変位，𝑃𝑃𝑖𝑖 は𝑖𝑖方向の表
面力を表している．𝑡𝑡0 ，𝑡𝑡𝑓𝑓 ，Ω， Γはそれぞれ，数値呼ばれる，自然積雪と車によって踏み固められた圧雪によって雪の山が形成される 1)．この雪の山は雪

解析における初期時刻，最終時刻，数値計算領域，
数値計算境界を表している．重み係数 𝑄𝑄𝑝𝑝 は，表面堤と呼ばれ，雪堤は，堆積・成長あるいは日中の気温の上昇により崩壊することがある．その際，崩壊力を規定する境界では1とし，他の境界では0を与えした雪堤が車線を塞ぎ，自動車への物損被害や交通る．また，偏微分は， (∙),𝑖𝑖 = 𝜕𝜕(∙)/𝜕𝜕𝑥𝑥𝑖𝑖 を用いて表し

事故などの重大な事故につながる可能性がある．している．

かしながら，雪堤の崩壊メカニズムは未だ解明されておらず，雪堤処理のタイミングや処理方法は作業

２．２

員の経験則に依存しているのが実状である 2)．近年

熱流体構造連成場における支配方程式として，
Navier-Stokes方程式，連続の式，熱伝導方程式，熱では，数値計算によるシミュレーションが発達しており，雪堤を対象とした融解および崩壊の挙動を再支配方程式

膨張を考慮した弾性場の運動方程式を考える．これらの支配方程式は式(3) ~ (6)のように表される．

現する解析ソフトを開発することで定量的に処理を行うことができると考えられる．
著者らはこれまでに，雪堤を熱弾性体と仮定し，

𝜓𝜓1 �𝑣𝑣̇ 𝑖𝑖 + 𝑣𝑣𝑗𝑗 𝑣𝑣𝑖𝑖,𝑗𝑗 + 𝑝𝑝,𝑖𝑖 − 𝜈𝜈𝑣𝑣𝑖𝑖,𝑗𝑗𝑗𝑗 + 𝛽𝛽𝛽𝛽(𝜙𝜙 − 𝜙𝜙0 )� = 0(3)

𝜓𝜓1 𝑣𝑣𝑖𝑖,𝑖𝑖 = 0 in 𝑡𝑡 ∈ [𝑡𝑡0 , 𝑡𝑡𝑓𝑓 ] in Ω(4)𝜓𝜓1 �𝜙𝜙̇ + 𝑢𝑢𝑖𝑖 𝜙𝜙,𝑖𝑖 − 𝜅𝜅𝜙𝜙,𝑖𝑖𝑖𝑖 � + 𝜓𝜓2 �𝜙𝜙̇ − 𝑘𝑘𝜙𝜙,𝑖𝑖𝑖𝑖 − 𝑄𝑄� = 0

in 𝑡𝑡 ∈ [𝑡𝑡0 , 𝑡𝑡𝑓𝑓 ] in Ω

(5)𝜓𝜓2 𝐶𝐶𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 �𝜀𝜀𝑘𝑘𝑘𝑘,𝑗𝑗 (𝑢𝑢𝑖𝑖 ) − 𝛼𝛼𝑘𝑘𝑘𝑘 (𝜙𝜙 − 𝜙𝜙𝑓𝑓 ),𝑗𝑗 � + 𝜓𝜓2 𝑓𝑓𝑖𝑖 = 0(6)

in 𝑡𝑡 ∈ [𝑡𝑡0 , 𝑡𝑡𝑓𝑓 ] in Ω

模擬雪堤の融解実験を撮影した画像とのフィッティ
ングにより，熱変形のファクターとなる熱膨張率の推定を行った 3)．また，その熱膨張率を用いて，雪堤周囲の流体からの熱の流入を考慮した熱流体構造連成問題の解析を行い撮影画像から求めた横幅の減
少量との比較し，融解過程を良好に再現した．
本論文では，熱変形に要したエネルギーであるひin 𝑡𝑡 ∈ [𝑡𝑡0 , 𝑡𝑡𝑓𝑓 ] in Ω

ずみエネルギーを評価関数に定義し，雪堤内部の構

ここで，𝜓𝜓1 ， 𝜓𝜓2 は重み定数を表しており𝜓𝜓1 = 1，

造物の形状を変化させることで，その総和を最小に𝜓𝜓2 = 0の場合は流れ場，𝜓𝜓1 = 0，𝜓𝜓2 = 1の場合は弾することを目的とした剛性最大化問題について検討する．まず，熱流体構造連成場における支配方程式性場となる．また，𝑣𝑣𝑖𝑖 ，𝑝𝑝，𝜈𝜈，𝛽𝛽，𝑔𝑔，𝜅𝜅，𝑘𝑘，𝑄𝑄はそれぞれ，各ステップにおける𝑖𝑖方向の流速，圧力，
動粘度，流体の熱膨張率，重力加速度，流体の熱拡を示し，随伴変数法および有限要素法を用いて，熱弾性場に着目した形状勾配関数を導出する．次に導出した形状勾配関数に基づいて力法 4)を適用し，内散率，弾性体の熱拡散率，発熱量を表している．時間微分は，(∙)̇ = 𝜕𝜕(∙)/𝜕𝜕𝜕𝜕を用いて表している．

部構造物の形状を決定する．

２．３

２．

剛性最大化問題の定式化

２．１

目的汎関数

境界条件および初期条件

自然対流場および熱弾性場における境界条件，初期条件を式(7) ~ (17)のように定義する．
𝑣𝑣𝑖𝑖 = 𝑣𝑣�𝑖𝑖 at 𝑡𝑡 = 𝑡𝑡0 in Ω

剛性最大化問題の目的汎関数 𝐽𝐽を，式(1)のように

定義する．

8

𝑣𝑣𝑖𝑖 = 𝑣𝑣�𝑖𝑖 in 𝑡𝑡 ∈ [𝑡𝑡0 , 𝑡𝑡𝑓𝑓 ] on Γ1v

(7)
(8)

𝑡𝑡𝑖𝑖 = 𝑡𝑡̂𝑖𝑖 = �𝑝𝑝𝛿𝛿𝑖𝑖𝑖𝑖 + 𝜈𝜈𝑣𝑣𝑖𝑖,𝑗𝑗 �𝑛𝑛𝑗𝑗 in 𝑡𝑡 ∈ [𝑡𝑡0 , 𝑡𝑡𝑓𝑓 ] on Γ2v(9)

𝑝𝑝 = 𝑝𝑝̂ at 𝑡𝑡 = 𝑡𝑡0 in Ω

(10)

𝜙𝜙 = 𝜙𝜙� in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γϕ(12)𝜙𝜙 = 𝜙𝜙� at 𝑡𝑡 = 𝑡𝑡0 in Ω(11)

𝑞𝑞 = 𝑞𝑞� = −𝑘𝑘𝜙𝜙,𝑖𝑖 𝑛𝑛𝑖𝑖 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γq(13)

𝑞𝑞 = ℎ��𝜙𝜙 − 𝜙𝜙�𝑜𝑜 � = −𝑘𝑘𝜙𝜙,𝑖𝑖 𝑛𝑛𝑖𝑖 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γh𝑓𝑓𝑖𝑖 = 𝑓𝑓̂𝑖𝑖 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � in Ω

𝑢𝑢𝑖𝑖 = 𝑢𝑢�𝑖𝑖 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γu

𝑃𝑃𝑖𝑖 = 𝐶𝐶𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 �𝜀𝜀𝑘𝑘𝑘𝑘 (𝑢𝑢𝑖𝑖 ) − 𝛼𝛼𝑘𝑘𝑘𝑘 (𝜙𝜙 − 𝜙𝜙𝑓𝑓 )�𝑛𝑛𝑗𝑗
in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on ΓP

in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 �

𝜙𝜙 ∗ = 0𝑢𝑢𝑖𝑖∗ = 0 at 𝑡𝑡 = 𝑡𝑡𝑓𝑓

𝜙𝜙 ∗ = 0𝑢𝑢𝑖𝑖∗ = 0

(14)
(15)

in Ω

on Γϕ

in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γu�∗ = �𝐶𝐶𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝛼𝛼𝑘𝑘𝑘𝑘 𝑢𝑢𝑖𝑖∗ − 𝑘𝑘𝜙𝜙,𝑖𝑖∗ �𝑛𝑛𝑖𝑖
𝑏𝑏 ∗ = 𝑏𝑏
in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γq and Γh

(16)

(21)
(22)
(23)
(24)
(25)
(26)

さらに，設計境界 Γdesign の座標に対するLagrange(17)

関数の勾配𝐺𝐺𝑖𝑖 は，式(27)から導出される．

熱伝達率，外気温度を表している．また，Γ1v ，Γ2v ，
Γϕ ，Γq ，Γh ，Γu ，ΓP はそれぞれ，流れ場におけるDirichlet 境界，流れ場における Neumann 境界，熱伝導 場に お け る Dirichlet 境界 ，熱 伝 導 場に お け る

𝑡𝑡𝑓𝑓

𝑡𝑡𝑓𝑓

� � 𝑡𝑡𝑖𝑖∗ 𝛿𝛿𝑢𝑢𝑖𝑖 𝑑𝑑Γ𝑑𝑑𝑑𝑑 + � � 𝑏𝑏 ∗ 𝛿𝛿𝛿𝛿 𝑑𝑑Γ𝑑𝑑𝑑𝑑
𝑡𝑡0

Γ
𝑡𝑡𝑓𝑓

=� �

Neumann 境界，熱伝導場における Robin 境界，熱弾

𝑡𝑡0

Γu +ΓP

𝑡𝑡0
𝑡𝑡𝑓𝑓

+� �
𝑡𝑡0

𝑡𝑡𝑖𝑖∗ 𝛿𝛿𝑢𝑢𝑖𝑖 𝑑𝑑Γ𝑑𝑑𝑑𝑑

Γϕ +Γq +Γh+�Γdesign

性 場に お け る Dirichlet 境界 ，熱 弾 性 場に お け るΓ𝑏𝑏 ∗ 𝛿𝛿𝛿𝛿 𝑑𝑑Γ𝑑𝑑𝑑𝑑(27)

𝐺𝐺𝑖𝑖 𝛿𝛿𝑥𝑥𝑖𝑖 𝑑𝑑Γ

𝑡𝑡𝑓𝑓

𝐺𝐺𝑖𝑖 = � �𝐶𝐶𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝛼𝛼𝑘𝑘𝑘𝑘 �𝜙𝜙 − 𝜙𝜙𝑓𝑓 � − 2𝐶𝐶𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝜀𝜀𝑘𝑘𝑘𝑘 (𝑢𝑢𝑖𝑖∗ )Neumann 境界を表している．

𝑡𝑡0

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

熱弾性体部分に着目し，任意の随伴変数𝜙𝜙 ∗ ，𝑢𝑢𝑖𝑖∗ を用いることで，目的汎関数𝐽𝐽はLagrange関数𝐽𝐽∗ に拡張+ 𝐶𝐶𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝛼𝛼𝑘𝑘𝑘𝑘 𝑢𝑢𝑖𝑖∗ − 𝑘𝑘𝜙𝜙,𝑖𝑖∗ � 𝑑𝑑𝑑𝑑𝑛𝑛𝑖𝑖(28)

滑らかな形状更新の勾配を得る力法を適用すること
で，勾配𝐺𝐺𝑖𝑖 を線形弾性体の外力として考え，形状更新の勾配𝐺𝐺𝑖𝑖∗ を得る．最急降下法に基づき，式(29)にされ，式(18)で表される．

示すように，修正された勾配𝐺𝐺𝑖𝑖∗ を用いて，設計境界

𝐽𝐽∗ = 𝐽𝐽

上の座標は更新される．ここで，𝑙𝑙は形状更新ステ

𝑡𝑡𝑓𝑓

+� �

𝑡𝑡0 Ω
𝑡𝑡𝑓𝑓

̇ 𝚤𝚤𝚤𝚤 − 𝑄𝑄�𝑑𝑑Ω𝑑𝑑𝑑𝑑
𝜙𝜙 ∗ �𝜙𝜙 − 𝑘𝑘𝜙𝜙+ � � 𝑢𝑢𝑖𝑖∗ �𝐶𝐶𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 �𝜀𝜀𝑘𝑘𝑘𝑘,𝑗𝑗 (𝑢𝑢𝑖𝑖 ) − 𝛼𝛼𝑘𝑘𝑘𝑘 �𝜙𝜙 − 𝜙𝜙𝑓𝑓 �,𝑗𝑗 �
𝑡𝑡0

in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 �in Ω𝑡𝑡𝑖𝑖∗ = 𝑡𝑡�𝐼𝐼∗ = �𝐶𝐶𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝛼𝛼𝑘𝑘𝑘𝑘 �𝜙𝜙 − 𝜙𝜙𝑓𝑓 � − 2𝐶𝐶𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝜀𝜀𝑘𝑘𝑘𝑘 (𝑢𝑢𝑖𝑖∗ )�𝑛𝑛𝑗𝑗
in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on ΓP� ，𝑛𝑛𝑖𝑖 ，𝑞𝑞，ℎ，𝜙𝜙𝑜𝑜 はそれぞれ，既知関数，
ここで，(∙)
境界における𝑖𝑖方向の単位法線ベクトル，熱流束，

２．４

at 𝑡𝑡 = 𝑡𝑡𝑓𝑓

in Ω

Ω

+ 𝑓𝑓𝑖𝑖 � 𝑑𝑑Ω𝑑𝑑𝑑𝑑

ップ回数，𝜂𝜂はステップ幅である．
(𝑙𝑙+1)

𝑥𝑥𝑖𝑖

(18)

(𝑙𝑙)

∗(𝑙𝑙)

= 𝑥𝑥𝑖𝑖 − 𝜂𝜂𝐺𝐺𝑖𝑖(29)

ここで，𝜙𝜙 ∗ は随伴温度，𝑢𝑢𝑖𝑖∗ は随伴変位を表す．紙

２．５

計算の流れ

一変分𝛿𝛿𝐽𝐽 は，目的汎関数の停留条件を満たすため方程式の計算）を実施し，得られた状態量を用いてに0となる必要があり，停留条件より随伴方程式お目的汎関数を計算することで，目的汎関数が収束す

よび随伴条件は次のように決定される．

るまで繰り返し計算を実施する．具体的な計算の流本研究では，構造の形状更新ごとに順解析（状態面の都合上，詳細は省略するが，Lagrange関数の第∗れを以下に示す．また，式(30)においてεは収束判定
𝜙𝜙 ∗ + 𝑘𝑘𝜙𝜙𝑖𝑖𝑖𝑖∗ − 𝐶𝐶𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝛼𝛼𝑘𝑘𝑘𝑘 𝜀𝜀𝑘𝑘𝑘𝑘 (𝑢𝑢𝑖𝑖 ) − 𝐶𝐶𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝛼𝛼𝑘𝑘𝑘𝑘 𝜀𝜀𝑖𝑖𝑖𝑖 (𝑢𝑢𝑖𝑖∗ ) = 0
in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � in Ω
𝐶𝐶𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝜀𝜀𝑘𝑘𝑘𝑘,𝑗𝑗 (𝑢𝑢𝑖𝑖∗ ) − 𝐶𝐶𝑖𝑖𝑖𝑖𝑖𝑖𝑖𝑖 𝛼𝛼𝑘𝑘𝑘𝑘 �𝜙𝜙 − 𝜙𝜙𝑓𝑓 �,𝑗𝑗 + 𝑓𝑓𝑖𝑖 = 0(19)定数を示す．

(20)
9

Step 1.
Step 2.
Step 3.

初期形状を与える．
計算条件を入力する．
状態方程式と初期条件および境界条件から有限要素解析を行う．
目的汎関数を計算する．
目的汎関数の収束判定を行う．式(30)に示す収束判定条件式を満たした場合，解析を終了する．

Step 4.
Step 5.

�𝐽𝐽(𝑙𝑙) − 𝐽𝐽(𝑙𝑙−1) �< 𝜖𝜖
𝐽𝐽(𝑙𝑙)

Step 6.

(30)

随伴方程式と終端条件および境界条件から有限要素解析を行う．
形状勾配関数𝐺𝐺𝑖𝑖 を計算する．
形状勾配関数𝐺𝐺𝑖𝑖 に対し，修正された形
状勾配関数𝐺𝐺𝑖𝑖∗ を計算する．
力法に基づいて形状更新を行い，Step 3に戻る．

Step 7.
Step 8.
Step 9.

図-1

３．

数値解析例

３．１

解析条件

表-1

-1のような雪堤を仮定し，雪堤の内部構造物の形状を力法によって修正された勾配𝐺𝐺𝑖𝑖∗ に基づいて更新す

る．

on Γv

𝑣𝑣1 = 𝑣𝑣2 = 0 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γsℎ = 0.001 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γh1 and Γh2

熱膨張率(流体)in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 �

(37)

外気温１
外気温２
時間ステップ
時間刻み幅

𝑢𝑢1 = 𝑢𝑢2 = 0 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γ \ Γs1 and Γs(38)

on Γs1

基準温度(34)

on Γv , Γs1 and Γdesign

𝑢𝑢2 = 0

熱拡散率

(33)

(35)

in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 �

プラントル数(32)(36)

𝑞𝑞 = 0粘性係数(31)𝑘𝑘 = 2.22 in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γs

Element
62544
67616
58744

た．図-2に測定結果を示し，推定した熱膨張率を式(39)，式(40)に示す．

数を加えて示す．また，境界条件は次のように与え

in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 �

Node
31613
34149
29713

定方法を基に各層に対して14回計測し，推定を行っし，各Caseのメッシュ作成時の節点数ならびに要素

𝑣𝑣2 = 0

b [m]
0.040
0.040
0.012

す．また，解析に使用する熱膨張率は文献(3)の推するため図-1中のa，b は，表-1に示すように設定

in 𝑡𝑡 ∈ �𝑡𝑡0 , 𝑡𝑡𝑓𝑓 � on Γh1 and Γh2

a [m]
0.020
0.010
0.040

表-2に，図-1の解析モデルにおける計算条件を示ることで最適化を行う．また，最適化の傾向を確認

𝑣𝑣1 = 0解析モデルの形状に関するパラメータ

Case 1Case 2
Case 3

雪堤内部に構造物がある解析モデルの例として図数値解析モデル

熱流束
ポアソン比
ヤング率
密度(上層)
密度(下層)
重力加速度

10

表-2
𝜈𝜈

数値解析条件𝑃𝑃𝑃𝑃𝜅𝜅𝛼𝛼𝐴𝐴𝐴𝐴𝐴𝐴

𝜙𝜙𝑏𝑏𝑏𝑏𝑏𝑏𝑏𝑏
𝜙𝜙1

𝜙𝜙2
𝑡𝑡𝑓𝑓

∆𝑡𝑡
𝑞𝑞

𝜈𝜈𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝𝑝
𝐸𝐸𝜌𝜌𝑢𝑢𝑢𝑢𝑢𝑢𝑢𝑢𝑢𝑢
𝜌𝜌𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙
𝑔𝑔

1.333 × 10−50.71
𝜈𝜈/𝑃𝑃𝑃𝑃3.59 × 10−3
273
274
278
20000
200
41.8680.15

1.64 × 10−1 × 𝑔𝑔 × 𝜌𝜌 5

0.5
0.3
9.80665[Pa ⋅ s][−][m2 /s]
[1/K]
[K]
[K]
[K]
[s]
[s]

[W/m2 ][−]

[Pa]

[g/cm3 ][g/cm3 ][m/s 2 ]定式化によって導出した勾配に対して力法を適用した，滑らかな勾配を用いて，雪堤と接する構造物の形状を変更し，目的汎関数の総和が最小となる剛性最大化問題について検討した．初期形状の異なる3
種類の雪堤モデルに対し，形状最適化を行い，それぞれの雪堤内部の構造物の形状が高密度の雪層がある上側に伸びるまたは，横に広がることによって雪堤全体の変形が抑制され，目的汎関数の減少が確認
できた．
図-2

各層における横幅の減少量
謝辞：本研究を行うにあたり，長岡技術科学大学

αupper = −2.658 × 10−10 𝑡𝑡 2 − 3.617 × 10−6 𝑡𝑡
α𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙𝑙 = −2.252 × 10−10 𝑡𝑡 2 − 7.645 × 10−6 𝑡𝑡

３．２

技術開発センタープロジェクト「高速道路雪堤の形(39)

成と不安定化のメカニズム解明と管理手法の検討」
における共同研究費（担当：長岡技術科学大学，長(40)

岡工業高等専門学校，㈱ネクスコ・エンジニアリン
グ新潟，㈱ネクスコ・メンテナンス新潟，町田建設㈱）の援助を受けた．ここに謝意を表す．

解析結果

図-3，図-4，図-5にそれぞれCase 1の初期形状における有限要素メッシュ，最終形状における有限要
素メッシュ，初期形状の目的汎関数と面積を1とし，形状更新の回数を横軸に取った時の収束履歴を示す．
また，図-6，図-7，図-8はCase 2，図-9，図-10，
図-11はCase 3の結果をCase1と同様に示す．
Case 1の形状では，初期形状に比べ根元部が細く，先端部が横に広がる形状が得られた．また，収束履
歴からわかるように，面積を一定に保ちながら，
10[%]の目的汎関数の改善が見られた．Case 2の形
状はCase 1と比べ細いモデルとなっているが，Case1と同様に根元部が細く，先端が横に広がるような形

図-3

初期形状における有限要素メッシュ [Case 1]

図-4

最終形状における有限要素メッシュ [Case 1]状が得られた．収束履歴では面積を一定に保ちながら，わずかに目的汎関数の改善が見られた．Case 3はこれまでの最適化の傾向とは異なり，内部構造物
が上部に伸びた形状が得られた．また，収束履歴を見ると面積を一定に保ちながら最終的に2.5[%]の目的汎関数の改善が見られた．これらの結果から，内部構造物の形状が上側にシフトすることにより高密度である上層の沈み込みが抑制され，雪堤全体の変形が小さくなったと考えられる．

４．

結論

本論文では，雪堤の崩壊の抑制のため，雪堤領域
内の構造物に着目し，熱変形に要するひずみエネル
ギーを目的汎関数に定義し，随伴変数法および有限要素法に基づいて形状最適化問題の定式化を行った．
11

図-5

目的汎関数と面積の収束履歴 [Case 1]
図-9

図-6

初期形状における有限要素メッシュ [Case 2]

図-7

最終形状における有限要素メッシュ [Case 2]

初期形状における有限要素メッシュ [Case 3]

図-10 最終形状における有限要素メッシュ [Case 3]

図-11

目的汎関数と面積の収束履歴 [Case 3]

参考文献
1)

芝崎智貴

他, “道路雪堤の崩壊現象の観察と分類”,

寒地技術論文・報告集, Vol.37, pp.169-174, 2021.
2)

西信衛, 武山浩二, “雪堤切りアタッチメントによる
中央分離帯雪堤処理の効率化”, (株)ネクスコ・メン
テナンス新潟, ゆき(97), 雪センター, pp.86-89, 2014.

3)

広瀬智史

他, “雪堤の融解現象に対する熱変形解析の適用（融解過程の撮影画像と解析結果のフィッテ
ィングによる検討）”, 日本機械学会第 34 回計算力

図-8

目的汎関数と面積の収束履歴 [Case 2]

学講演会講演論文集, No.21-36, 2021.
4)

畔上秀幸, “領域最適化問題の一解法”, 日本機械学会論文集 A 編, Vol.60, No.574, pp.1479-1486, 1994.
(2022. 8. 8 受付)

12