長岡工業高等専門学校研究紀要

論第 52 巻（2016）

文

凍結工法を用いた凍土形成時における凍結管温度の制御シミュレーション－ 随伴変数法および有限要素法による境界温度制御解析 －倉橋 貴彦1・石津 慎2・丸岡 宏太郎3・島村 淳4・土屋 勉5・河田 剛毅6
1 長岡技術科学大学大学院（Department of Mechanical Engineering, Nagaoka University of Technology）
2 長岡技術科学大学大学院（Department of Mechanical Engineering, Nagaoka University of Technology）
3 長岡技術科学大学大学院（Department of Mechanical Engineering, Nagaoka University of Technology）
4ケミカルグラウト株式会社（Chemical Grouting Co., Ltd.）
5ケミカルグラウト株式会社（Chemical Grouting Co., Ltd.）
6 機械工学科（Department of Mechanical Engineering, National Institute of Technology, Nagaoka College）

Numerical Temperature Determination of Frozen Ductin Frozen Soil Formation Based on Ground Freezing-Boundary Temperature Control Analysis
Based on Adjoint Variable and Finite Element MethodsTakahiko KURAHASHI1, Shin ISHIZU2, Kotaro MARUOKA3,Atsushi SHIMAMURA4, Tsutomu TSUCHIYA5 and Yoshitaka KAWADA6
Abstract
In this study, we focus on the Ground freezing that is one of the soil improvement method, and carry out
the boundary temperature control analysis on frozen duct in frozen soil formation based on the adjoint
variable and the finite element methods. The heat transfer equation is introduced to represent thetemperature diffusion in soil, and the finite element method is applied to simulate the temperature
distribution in the target domain. In addition, the adjoint variable method is employed to obtain theappropriate frozen duct temperature such that the temperature at the target point near the frozen duct isclose to the target temperature, i.e., the freezing temperature.
Key Words : ground freezing, heat transfer equation, finite element method, adjoint variable method

１．序論

工法はトンネル掘削時に適用され，強度の低い地盤周辺に凍結管を埋設し凍結管周辺の地盤を一時的に近年，環境にやさしい地盤改良工法として，凍土凍結し強度を増加させる目的で実施される．トンネ形成による不透水性の向上や地盤強度の増強を目的

ル掘削後は，自然に地下水へ再生できることもあり，とした工法（凍結工法 ）が注目されている．この環境負荷のない工法として知られている．しかし，

1)

- 1 -

倉橋 貴彦・石津 慎・丸岡 宏太郎・島村 淳・土屋 勉・河田 剛毅

地下水流速がある一定以上の地盤では遮水が困難であったとの報告もあり，現状では克服すべき課題が温度固定境界（凍結管表面を表す境界）を示す．また行列[M]は質量行列，[Hxx]，[Hyy]は拡散行列，{T}まだ残されている．そこで，本研究では，凍結管近は熱流束ベクトル，φ uは制御温度，bは熱流束，ハ傍において設定した目的点における温度を目標温度

ット記号は既知量を示す．にするような凍結管表面温度を随伴変数法 2)および有限要素法 により求めることを目的とする．

ここで，式(1)に示す評価関数の最小化問題において，式(2)，(3)を制約条件として付加するために，

3)

随伴変数λを導入し，ラグランジュ関数を誘導する．

２． 凍結管表面温度の計算アルゴリズム

J

式展開の始めとして，まず凍結管近傍において設

1 t f nel   obj. Te Qe    obj. e dt

t02 e 1

定した目的点における数値計算結果として得られる

tft0

(4)

 M   H xx   H yy  T T

温度φと目標温度φobj.の差の二乗値を時間積分した値を評価関数として定義する（式(1)）．

式(4)は汎関数であるため，停留条件を計算するために第一変分δJ*を計算すると式(5)に示す随伴方程J

1 t f nelT   obj. e Qe   obj. e dt
2 t0 e 1

式と，式(6)に示す随伴変数の終端条件と境界条件(1)が得られる．

   Q      S  in  t  t , t TTT M     H xx   H yy 

ここに，中カッコはベクトル，大カッコは行列，T(5)T

は転置を示しており，t0，tfは解析を始める時間およobj,

び終端時間を示す．eは各要素の番号を示しており，[Q]は重み対角行列を示す．重み対角行列は，対角0 t f   0 in   0 on d t  t 0 , t f
  0 on cont. t  t 0 , t f*s  J  0 on  t  t , tn0f

値において目的点において1，それ以外の点では0が設定される．ここで，計算の目的は式(1)に示す評
価関数の値が最小となるように凍結管表面の温度を求めることになる．評価関数の値が最小となるということは，目的点において数値計算結果として得ら

f(6)



れた温度φと目標温度φobj.の値が近づくことを意味している．ここに，数値計算結果として得られた温ここに，sは随伴変数に対する流束を示しており，度φは有限要素法による熱伝導解析から得られる温ラグランジュ関数の温度に対する勾配として表され度とすると，式(1)に対する制約条件として式(2)，(3)を導入する必要がある．

る．また，随伴変数に対しては終端時間tfにおける条件が付加されるため，終端時間から逆時間方向に随伴方程式（式(5)）を解くことになる．

M   H xx   H yy   T  in  t  t 0 , t f 

以上に示す式(2)を計算することで，式(5)の計算を行うことができ，式(5)の計算結果から凍結管表(2)面を表す境界Γ cont. 上におけるラグランジュ関数の温度に対する勾配sの値を求めることにより，凍結

 t 0   ˆt 0  in ˆ   on d t  t 0 , t f
  u on cont. t  t 0 , t fb  ,i ni  bˆ on n t  t 0 , t f




 

管表面温度を逐次更新するという流れで計算が行われる．本論文においては，Sakawa-Shindo法4)を適用(3)し，式(7)のように凍結管表面温度の更新を行う．

    W   J  on ( l 1)
u(l )
u

1

*  u (l )con.

t  t0 , t f

ここに，Ωは全計算領域，Γ dは温度固定境界，Γ n(7)は熱流束を与える境界，Γ cont. は制御温度を与える

本論文においては，式(2)および式(5)の時間方向の- 2 -凍結工法を用いた凍土形成時における凍結管温度の制御シミュレーション

離散化には陽的オイラー法を適用する．以下に，計が得られ，左側から入力している正弦波の逆の形

算アルゴリズムを示す．

（正弦波にマイナスを乗じた波形）が得られていることがわかる．計算領域中央の点（目的点）における温度の経時変化を図-4に示す．最終イタレーショ[凍結管表面温度を求めるための計算アルゴリズム]

ン時の目的点における温度（図-4の実線）は零度に1.

反復計算における初期制御温度φuの入力なっており，目的点における温度の制御が適切に行

2.

有限要素法により離散化された熱伝導方程式えていることがわかる．また，図-5に0イタレーシ（式(2)）を計算する．

ョン時（右側の境界温度を零度とした場合）と，最3.

評価関数J(0)を計算する．

終イタレーション時（図-3に示す実線の波形を入力

4.

有限要素法により離散化された随伴方程式（式(5)）を計算する．

した場合）の計算領域における温度分布を示す．温凍結管表面を表す境界Γcont.上におけるラグランジ ュ 関 数 の温 度 に対 す る勾 配 sの 値 を 用 いて制御ありの場合には計算領域中央の点において温度

Sakawa-Shindo法による温度の更新を行う．（式

検証をもとに，凍土形成時における凍結管温度の制(7)）

御シミュレーションへ本解析手法を適用する．

5.

6.

度分布が正しく得られていることを確認でき，また，が常に零度になっていることを確認できる．以上の反復回数(l+1)と(l)における評価関数Jの差の絶対値を計算し，各問題において設定する許容判

表-1 計算条件

定値εより小さい場合は計算を終了する．そう

総節点数
総要素数
時間増分量Δt，sec.
熱拡散率，m2/ sec.
時間ステップ
Sakawa-Shindo法における1イタレーション目の重み関数の対角値
1イタレーション目の制御境界における温度，deg.
収束判定定数εでない場合は，ステップ7へ進む．
7.

ステップ6において更新された凍結管表面温度の値を使用し，有限要素法により離散化された熱伝導方程式（式(2)）を計算する．

8.

評価関数J(l+1)を計算する．

9.

重みパラメータW(l)を更新する．J(l+1)<J(l)の場合
はWl+1)= W(l)×0.9としてステップ4へ戻る．また，
J> J の場合はW

(l+1)

(l)

= W ×2.0としてステップl+1)(l)5へ戻る．

33
40
1.0
0.002
2000
20.0

0.0
10-6

３．数値計算例; 目的点

3.1 境界温度制御解析の検証例題
まず，図-1に示す矩形領域の温度制御の問題に対n計算領域してプログラムの妥当性について検証を行う．左側から周期は2000sec.の正弦波による温度波形を与え，


矩形領域の中心線上の節点に目的点を設定し，目的点における温度を零度にする右側境界からの制御温cont.

d度を求める．計算条件を表-1に示す． なお，前章の計算アルゴリズムのステップ6にて述べた収束判定については，本検証では各反復回数における右側
境界温度の差の最大ノルムと収束判定定数との比較

nにより行う．
数値実験による検討を行った結果，図-2に示すよL2
L

うに評価関数の値（0イタレーション時における評価関数の値により正規化した値）は短調に減少し，収束していることがわかる．この反復計算により，
図-3に示す制御温度（実線：最終イタレーション時）- 3 -L2

図-1

計算モデル図u ?

倉橋 貴彦・石津 慎・丸岡 宏太郎・島村 淳・土屋 勉・河田 剛毅

図-2

図-3

評価関数の収束履歴(a-1) 0イタレーション時：Time=400sec.

（図-5-1）(a-2) 0イタレーション時：Time=800sec.

（図-5-2）(a-3) 0イタレーション時：Time=1200sec.

（図-5-3）境界 Γcont.における制御温度（実線：最終イタレーション時）と境界Γd における流入温度（点線）の比較

図-4

目的点(x=0.5m，y=0.1m)における温度の経時変化（実線：最終イタレーション時，点線：0 イタレーション時）- 4 -凍結工法を用いた凍土形成時における凍結管温度の制御シミュレーション(a-4) 0イタレーション時：Time=1600sec.

（図-5-4）(b-2) 最終イタレーション時：Time=800sec.

（図-5-7）(a-5) 0イタレーション時：Time=2000sec.

（図-5-5）(b-3) 最終イタレーション時：Time=1200sec.

（図-5-8）(b-1) 最終イタレーション時：Time=400sec.

（図-5-6）(b-4) 最終イタレーション時：Time=1600sec.

（図-5-9）- 5 -倉橋 貴彦・石津 慎・丸岡 宏太郎・島村 淳・土屋 勉・河田 剛毅となった．これも地表面で与えている温度の影響と考えられ，この結果より地表面近くは地中内の温度を制御することが困難であることを想像できる．また，図-11に0イタレーション時および最終イタレーション時の計算領域における温度分布を示す． 図の温度表示は，温度分布を分かりやすく表示するた
めに－20度を下限として表示しており，凍結管表面においては，－20度以下の箇所もあることに注意する．結果より0イタレーション時に比べ最終イタレーション時は凍結管表面において温度が低くなっており，また最終イタレーション時の結果より凍結管表面の上部近傍（図-7における点A近傍）は，凍結管表面のその他の場所に比べて温度が低くなっていることが確認できる．また，図-10および図-11（図-11-8）において，解析時間終端(T=86400秒の時)の(b-5) 最終イタレーション時：Time=2000sec. （図-5-10）

最終イタレーション時に，凍結管表面において温度図-5が極端に低く得られているが，これは随伴変数の終0 イタレーション時と最終イタレーション時端条件（式(6)）に起因するものであると考えられ，の温度分布図

この点に対する検討は今後の課題とする．
3.2 凍結管表面温度制御解析への適用
次に，図-6に示す矩形領域の地盤内に凍結管を埋表-2 計算条件

総節点数
総要素数
時間増分量Δt，sec.
熱拡散率，m2/ sec.
時間ステップ

設したモデルについて検討を行う．x=0.2m上の中心線に対して対象モデルのため，半分の領域のモデル(1/2モデル)に対して検討を行う．図-7に計算領域および初期・境界条件の定義，また制御点を設定する
境界，目的点を設定する境界の代表点において結果を出力する箇所（点A～F）の位置を示す．計算条

Sakawa-Shindo法における1イタレーション目の重み関数の対角値
1イタレーション目の制御境界における温度，deg.
収束判定定数ε件を表-2に示す．
結果として，図-8に示すように評価関数の値（0イタレーション時における評価関数の値により正規化した値）は短調に減少し収束した．評価関数収束時における図-7の代表点A，B，Cにおける制御温度

2846
5470
1.0102×10-886400(24時間)
20000.0

0.010-6の経時変化を図-9-1に示す．また0～200秒における0.4m

制御温度の経時変化を図-9-2に示す．地表温度を200.2m

度と温度固定しているため，点Aにおける温度が一番低く得られていることがわかる．また，点B，点Cにおいては，点Aの結果を踏まえると，地表面かFrozen duct

ら一番遠い方（点C）が制御温度は一番高くなる（零度に近くなる）と想定されるが，解析結果では，0.35m

点Bの方が点Cに比べて制御温度が高くなる結果となった．図-10-1に0イタレーション・最終イタレー
ション時における点D，E，Fにおける温度の経時変化を示す．また，0～200秒における温度の経時変化を図-10-2に示す．結果より最終イタレーション時
は，点B，Cでは目標とする温度である零度に近づけられていることを確認できるが，点Aの結果をその他のB，C点に比べて零度より若干差がある結果- 6 -y0.025m0.01mx

図-6計算モデル図

凍結工法を用いた凍土形成時における凍結管温度の制御シミュレーション

Initial condition in computational domain: 20deg.

Boundary on fixed temperature points (20deg.)
Point A ((x,y)=(0.195m,0.345m))
Point D ((x,y)=(0.180m,0.345m))
Boundary on target points (x=0.180m)

Boundary on control points (x=0.195m)
Point B ((x,y)=(0.195m,0.185m))
Point E ((x,y)=(0.180m,0.185m))

Point C ((x,y)=(0.195m,0.025m))
Point F ((x,y)=(0.180m,0.025m))

Adiabatic boundary(Heat flux = 0W/m2)

図-7(b)0～200秒における制御温度の変化（図-9-2）
0.02m (2.0cm)

図-9計算領域および初期・境界条件の定義

図-8

評価関数の収束履歴

点 A・B・C における制御温度の経時変化(a)全時間における温度変化（図-10-1）(a)全時間における制御温度の変化（図-9-1）(b)0～200秒における温度変化（図-10-2）
図-10

- 7 -

点 D・E・F における温度の経時変化

倉橋 貴彦・石津 慎・丸岡 宏太郎・島村 淳・土屋 勉・河田 剛毅(a-1) 0イタレーション時：Time=0.5h.

（図-11-1）(a-4) 0イタレーション時：Time=24h.

(a-2) 0イタレーション時：Time=1h.

（図-11-2）(b-1) 最終イタレーション時：Time=0.5h.

（図-11-5）(a-3) 0イタレーション時：Time=12h.

（図-11-3）(b-2) 最終イタレーション時：Time=1h.

（図-11-6）- 8 -（図-11-4）

凍結工法を用いた凍土形成時における凍結管温度の制御シミュレーション[矩形領域のモデル]
・随伴変数法に基づく制御計算より，目的点における温度を目標温度とするような制御温度が適切に計算できることがわかった．
[凍結工法のモデル]
・制御温度を地表近くと，凍結管の中間と末端において比較したところ，地表に近いほど温度が一番低く得られており，次に低いのは末端で，一番高い温度は凍結管の中間部で得られていることがわかった．
・地表面付近においては，最終イタレーション時における制御温度を入力しても，目標としている温度である零度近傍に制御はできなかったが，地表面から遠い目的点では，目標値の零度近傍に温度を制御(b-3) 最終イタレーション時：Time=12h.

できていることがわかった．

（図-11-7）

付録 凍土形成時における温度履歴の測定
砂の粒径の違いに対する凍土の形成過程を比較するために，マイナス19℃の低温室において，水を含ませた桐生砂と川砂に対する凍土の形成実験を行った．本実験において用いた各砂の粒径を表-3に示す．
実験では，幅400mm×奥行250mm×深さ300mmの水槽を2つ用意し，周囲と底面を断熱材で覆い，各水槽に100mmの深さまで砂を敷き詰め，水を砂表面まで充填した後に，凍土形成の実験を行った．
各砂に対して，砂表面から20mmの所の温度を熱電対により計測した．温度履歴の結果を図-12に示す．図-12は最初の温度をゼロ度としたものであり，(b-4) 最終イタレーション時：Time=24h.

基準温度からの変化に対する履歴を示している．結（図-11-8）

果より川砂の方が桐生砂に比べて温度の下がるスピ
図-11ードが速いことがわかる．また，150分後の桐生砂，0 イタレーション時と最終イタレーション時

川砂の表面の様子を図-13，図-14に示す．結果より，の温度分布図川砂（粒径が細かい）の方が桐生砂（粒径が粗い）に比べて凍り方が速い様子も確認できる．これらの４．結論

結果より，粒径の細かい方が凍土形成には有利であると考えられる．

本論文では，凍結工法を用いた凍土形成時における凍結管温度の制御シミュレーションに着目し，随表-3 実験に用いた砂の粒径伴変数法および有限要素法に基づき凍結管表面境界における制御温度について検討を行った．まず初めに，矩形領域における温度の制御問題に関する解析について検証を行い，凍結工法における凍結管周りの温度の制御問題について検討を行った．本研究により得られた知見を以下に示す．
- 9 -

桐生砂
川砂

粒径20mm以下
2mm以下

倉橋 貴彦・石津 慎・丸岡 宏太郎・島村 淳・土屋 勉・河田 剛毅

1

開始温度をゼロ度とした時の温度[℃]

0-1 0

2000

4000

6000

8000

謝辞：本論文を執筆するにあたり「公益財団法

10000

人

-2-3-4-5

内田エネルギー科学振興財団」および「次世代川砂（粒径2mm以下）ものづくり技術の基盤となる超高信頼性材料創成事

桐生砂（粒径20mm以下)業（長岡技術科学大学）」の援助を受けた．本論文-6の解析結果は，九州大学情報基盤研究開発センター-7の高性能アプリケーションサーバSR16000を使用し-8-9図-12経過時間[s]計算を行ったものであり，また，凍土の形成実験は表面からの深さ 20mm に対する長岡工業高等専門学校の低温室を借り行ったもので各砂の温度変化の履歴ある．計算機および低温室の使用について，関係者の方へ謝意を表す．
参考文献
1)

凍土分科会，凍土の知識－人工凍土壁の技術－，日本雪氷学会誌

2)

雪氷，Vol. 76, No. 2, pp.179-192, 2014.

T.Kurahashi and M.Kawahara, Examinations for Terminal
Condition of Lagrange Multiplier for Heat Transfer ControlProblems, International Journal for Numerical Methods inEngineering, Vol.73, pp.982-1009, 2008.

図-13

桐生砂

150min 後の砂表面の様子

3)

竹内 則雄, 樫山 和男, 寺田 賢二郎，”計算力学 有限要素法の基礎”，森北出版, pp.175, 2003.

4)

Y.Sakawa and Y.Shindo, On Global Convergence of An
Algorithm for Optimal Control, IEEE Transactions on
Automatic Control, Vol. ac-25, No.6, pp.1149-1153, 1980.
(2016. 9. 16 受付)

図-14

川砂

150min 後の砂表面の様子- 10 -