JP6230025B2 - 熱応答試験の解析方法および解析プログラム - Google Patents

Description

本発明は熱応答試験の解析方法および解析プログラムに関する。さらに詳しくは、地下水流動による移流・分散の影響を考慮に入れた解析を行うことができ、パウエルの共役傾
斜法に基づく非線形最適化手法を用いて未知のパラメータを同定することによって、短時間で精度良く地盤の熱交換特性を評価できる解析方法および解析プログラムに関する。

近年、原子力や火力に代わるクリーンなエネルギーとして、自然エネルギー(再生可能エネルギー）の利用が見直されてきている。その１つとして、地中熱の利用に関心が高まってきている。季節によって変動する外気温に対して、地中の温度は年間を通して約１５℃で安定している。このような特性を持つ地中の熱を有効利用したものとして、地中熱ヒート・ポンプ（Ｇｅｏｔｈｅｒｍａｌ Ｈｅａｔ Ｐｕｍｐ、以下ＧｅｏＨＰと略す）システムがある。これは、地中熱を利用した冷暖房・給湯システムのことである。安定した熱の供給が行なえる地中熱の特性に着目し、夏場には冷房として地中に熱を放出し、冬場には暖房の熱源として利用する。

しかしながら、日本ではＧｅｏＨＰシステムの普及は、欧米に比べ、遅れているのが実情である。その理由としては、認知度の低さ、井戸の掘削等における初期費用の問題、および、日本では欧米とは違い、その土地の風土によって熱物性値が異なるため、経済的なシステムづくりが依然として確立されていない、などが挙げられる。そのため、ＧｅｏＨＰシステムの普及には初期費用を抑えるためのシステムの効率化が必要不可欠である。

システムの効率化を図るための地盤調査試験として、熱応答試験が知られている。熱応答試験は主に地層の熱伝導率と熱交換井の熱抵抗を評価するための試験であり、地中熱交換井の熱交換挙動を予想する上で不可欠な情報である。熱応答試験を実施して、適正な井戸の本数・長さを決定することは、ＧｅｏＨＰシステムの初期費用の削減に極めて重要である。

地盤調査のための熱応答試験は、主にクローズド型方式における垂直型地中熱交換井を対象としている。このための熱応答試験装置は、地上部の水タンク内で電気ヒーターによって熱媒体を加熱し、循環ポンプにより配管内を循環させる。加熱された熱媒体は地中部を通過する際に放熱するため、Ｕ字管の入口・出口の温度差を測定し、その地層の熱交換特性を評価する。このような熱応答試験装置は、特許文献１に開示されており、熱応答試験装置を用いた熱応答試験方法は非特許文献１に開示されている。

熱応答試験の代表的な解析方法としては、地盤の温度上昇量を規定するＫｅｌｖｉｎの線源関数や円筒型熱源関数に基づく解析的手法、および有限差分法や有限要素法を用いた数値モデルによる解析方法がある。熱応答試験の解析的手法において最も一般的に用いられている手法としては、データ観測の簡便さや解析の容易さより、Ｋｅｌｖｉｎの線源関数に基づいた作図法が挙げられる。作図法では循環時における熱媒体の温度データの経時変化を用いる方法と循環停止後の地中温度の回復データを用いる方法があるが、循環時のデータを用いる解析が一般的に行われている。

本発明者らは、循環時データを用いた作図法で熱応答試験の解析を行う場合に、解析を行うデータ区間の選定に明確な基準が無いため、解析結果が人によって異なること、地中熱交換器や充填剤の影響を排除するため、試験期間の初期段階のデータは捨てるので試験
期間を長くとる必要があること、などの問題点に着目し、地盤の温度上昇量を規定するＫｅｌｖｉｎの線源関数の近似解における未知のパラメータを同定して、地盤の熱伝導率を算出する解析方法および解析プログラムを提案している（特許文献２参照）。特許文献２で提案した解析方法では、熱応答試験の初期段階において得られる実測値に対して、パウエルの共役傾斜法に基づく非線形最適化手法を用いて、未知のパラメータを所定の値に設定した場合に得られる温度上昇量の計算値と、熱応答試験によって測定された温度上昇量の実測値との誤差が最小となるように、パラメータの逆解析を行い、未知のパラメータを同定して、地盤の熱伝導率を算出する。

また、本発明者らは、透水量係数や貯留係数といった地下水流動系パラメータを求めるための原位置試験として最も一般的に行われている揚水試験は、熱応答試験と同様に、熱伝導理論を基として解析解を導いていることに着目し、特許文献２において、熱応答試験の解析方法を揚水試験にも適用することを提案している。すなわち、地盤の地下水位低下量を規定するタイスの井戸関数の近似解（ヤコブの近似解）における未知のパラメータを同定して、地盤の透水量係数を算出するにあたって、揚水試験の初期段階において得られる実測値に対して、パウエルの共役傾斜法に基づく非線形最適化手法を用いて、未知のパラメータを所定の値に設定した場合に得られる地下水位低下量の計算値と、揚水試験によって測定された地下水位低下量の実測値との誤差が最小となるように、パラメータの逆解析を行い、未知のパラメータを同定して、地盤の透水量係数を算出している。

特開２００４−３０１７５０号公報 特許第５３３４２２１号公報

「ボアホール型地中熱交換機に対する加熱法による熱応答試験の標準試験方法ｖｅｒ．２」、長野克則著、地下熱利用とヒートポンプシステム研究会編、財団法人ヒート・ポンプ・蓄熱センター

従来から行われている作図法による熱応答試験の解析方法や、本発明者らが特許文献２で提案したパウエルの共役傾斜法に基づく非線形最適化手法を用いたパラメータ同定による熱応答試験の解析方法は、Ｋｅｌｖｉｎの線源関数に基づいた熱応答試験の解析理論を前提としている。そして、この解析理論は、熱移動が熱伝導のみに支配されることを前提としている。このため、地下水流動による移流・分散の影響を考慮に入れた解析を行うことができない。従って、地下水流動がある場合には、熱応答試験の解析により求まる熱伝導率や体積熱容量が誤差を含んだものとなってしまう。また、地下水流速や縦・横分散長等のパラメータを求めることはできない。このため、地下水流動があるサイトでの熱交換井の熱交換挙動を精度良く予測することができない。

本発明の課題は、このような点に鑑みて、地下水流動がある場合でも、熱伝導率や体積熱容量、地下水流速、縦・横分散長などのパラメータを精度良く同定することができる、熱応答試験の解析方法および解析プログラムを提案することにある。

上記の課題を解決するために、本発明の熱応答試験の解析方法および解析プログラムでは、地盤の温度上昇式として、地下水が流動している流動場での移流・分散現象を考慮した熱移動方程式の解析解を用いる。そして、この解析解の近似式（高次漸近解）を用いて求められる地盤の温度上昇値の計算値と、熱応答試験で得られる地盤の温度上昇値の実測値との誤差が最小となるように、熱移動方程式の解析解に含まれる未知のパラメータを同定する逆解析を行う。逆解析には、パウエルの共役傾斜法に基づく非線形最適化法を用いる。

（移流・分散現象を考慮した熱移動の基礎理論に基づく熱応答試験の解析方法）
以下に、地盤の温度上昇式として、地下水などの流体の流れに伴う熱移動を考慮した熱移動方程式の解析解を求めるにあたって、本発明者らが用いた解析理論を述べる。

飽和・不飽和多孔体中における熱移動には熱伝導、熱分散、熱交換、熱対流といった現象を考慮する必要がある。図１は、空気や水などの流体で飽和した多孔体中を流体が流動している時の熱の流れを模式的に表したものである。図中のルートＡは熱伝導による固相中の熱の流れ、ルートＢは熱伝導による流体相中の熱の流れ、ルートＣは固相と流体相の接触面を通して行われる熱交換、ルートＤは流体の流れに伴う機構的分散による熱移動を表している。

まず、空気や水などの流体によって飽和した多孔体中における熱移動において固相・液相・気相の温度が平衡状態にあるものとすると、飽和多孔体中における熱移動の基礎方程式は、移流項と分散項を含む式（１）の移流分散方程式で表される。
但し、 Ｔ：温度（℃）
ｖ：間隙流速（ｍ／ｓ）
ｎ：間隙率
（ρｃ）_Ｗ：液相の体積熱容量（Ｊ／ｍ^３／Ｋ）
（ρｃ）_ｅ：固相・液相・気相混合系の等価体積熱容量（Ｊ／ｍ^３／Ｋ）
λ_ｈ：固相・液相・気相の熱伝導率に機構的熱分散を加味した混合系の熱分散率
（Ｗ／ｍ／Ｋ）

ここで、Ｘ軸方向の地下水流動がある２次元流動場を考え、式（１）の両辺を（ρｃ）_ｅで割ると、次式（２）、（３）で表される飽和多孔体の熱移動方程式が得られる。

Ｘ軸方向に間隙流速ｖで地下水が流動している２次元流動場において、原点に与える熱負荷として、微少時間ｄτで熱量ｑ_Ｈが連続的に投入されるものとし、式（４）〜（７）の初期条件・境界条件のもとで式（２）を解くと、２次元流動場の温度差分布として、次式（８）の解析解が得られる。

式（８）において、次式（９）のように置くと、式（８）は、次式（１０）で表されるＷ（ｕ、Ｂ）を含んでいる。Ｗ（ｕ、Ｂ）はハンタッシュ−ヤコブ（Ｈａｎｔｕｓｈ−Ｊａｃｏｂ）の井戸関数と呼ばれ、隣接帯水層から漏水がある被圧帯水層中の透水量係数Ｔ_ｍと貯留係数Ｓ_ｃを求めるための揚水試験の解析解として知られている。

ここで、式（８）において、ｖ＝０とし、原点からの距離をｒとすると、式（８）は次式（１１）、式（１２）のようになる。式（１１）は、Ｋｅｌｖｉｎの線源関数そのものである。Ｋｅｌｖｉｎの線源関数は、従来から熱応答試験の解析に用いられてきた、移流・分散現象を考慮しない熱移動方程式の解析解である。

式（１１）に含まれるＷ（ｕ）は、タイス（Ｔｈｅｉｓ）の井戸関数である。式（１０）で表されるＷ（ｕ、Ｂ）において、Ｂ＝０とすると、Ｗ（ｕ、０）＝Ｗ（ｕ）となり、タイスの井戸関数と一致する。タイスの井戸関数は、従来から揚水試験の解析に用いられてきた、隣接帯水層の干渉を考慮しない地下水流動方程式の解析解であり、Ｋｅｌｖｉｎの線源関数と式の形が同じである。このように、本発明で用いる、移流・分散現象を考慮した熱移動方程式の解析解は、隣接帯水層の干渉がない条件下では、従来用いられていた、地下水流れを考慮しない熱応答試験や揚水試験の解析解と一致する。

また、地下水流動場において、原点に連続的に注入したトレーサの濃度差分布は、熱移動と同じ解析理論によって解析することが可能である。すなわち、上記の熱移動方程式の解析解のパラメータを、物質輸送式に用いるパラメータに置き換えることで、移流・分散現象を考慮した物質輸送式の解析解が得られる。

（高次漸近解の適用）
式（１０）で表されるハンタッシュ−ヤコブの井戸関数Ｗ（ｕ、Ｂ）は、ｕ≦１．０、Ｂ≦２．０の範囲において、次式（１３）で近似される。

式（１３）において、Ｉ_０（Ｂ）は第１種変形ベッセル（Ｂｓｓｅｌ）関数、Ｋ_０（Ｂ）は第２種変形ベッセル（Ｂｓｓｅｌ）関数である。また、上述したように、Ｗ（ｕ）はタイスの井戸関数である。ここで、式（１３）のＩ_０（Ｂ）、Ｋ_０（Ｂ）、Ｗ（ｕ）は、それぞれ、次式（１４）〜（１６）のべき級数に展開できる。

本発明者らは、式（１４）〜（１６）のべき級数展開式の正確解と、べき級数展開式の第３項、第３０項、第５０項・・・などの項数までを用いて近似した近似解との一致度を検討した。図２．１にＩ_０（Ｂ）の正確解と第３項、第４項、第３０項、第５０項までを用いて近似した近似解を示す。また、図２．２にＫ_０（Ｂ）の正確解と第３項、第４項、第３０項、第５０項までを用いて近似した近似解を示す。Ｉ_０（Ｂ）では、第３項までの近似解はＢが小さい場合は正確解と良く一致するが、Ｂが大きくなるにつれて正確解から乖離する。また、Ｋ_０（Ｂ）では、第３項、第４項までの近似解はＢが小さい場合は正確解と良く一致するが、Ｂが大きくなるにつれて正確解から乖離する。Ｉ_０（Ｂ）、Ｋ_０（Ｂ）とも、第３０項、第５０項まで近似した場合はどの区間においても正確解と良く一致する。Ｉ_０（Ｂ）、Ｋ_０（Ｂ）は、ともに第５０項まで近似した近似解を用いることが好ましい。

図２．３にタイスの正確解と第２項、第６項、第３０項、第１００項までを用いて近似した近似解を示す。第２項までの近似解は、１／Ｕが大きい後半部では正確解と良く一致するが、１／Ｕが小さくなるにつれて正確解から乖離する。近似する項数を増加させると近似解は正確解に漸近し、第３０項、第１００項までで近似した場合はどの区間においても正確解と良く一致する。

図２．４に、式（１３）で示すＷ（ｕ、Ｂ）の正確解（ハンタッシュ−ヤコブの標準曲線）と、Ｉ_０（Ｂ）、Ｋ_０（Ｂ）、Ｗ（ｕ）のべき級数展開式を高次項（少なくとも、第３０項）まで用いて近似した近似解（漸近解）との比較結果を示す。このように、漸近解の項数を増加させることで、正確解と非常に良く一致する計算値が得られる。本発明では、以下に説明する逆解析を行うにあたって、Ｉ_０（Ｂ）、Ｋ_０（Ｂ）、Ｗ（ｕ）のべき級数展開式を、少なくとも第３０項まで（Ｉ_０（Ｂ）、Ｋ_０（Ｂ）については好ましくは第５０項まで）用いた解析プログラムを用いる。

（パウエルの共益傾斜法）
本発明においては、熱応答試験あるいは揚水試験の解析を、最適化手法による逆解析を利用して行う。先に述べた地下水モデルによる熱応答試験の解析解や、物質輸送式の解析解、あるいは、後述する揚水試験の解析解は、透水量係数Ｔ_ｍや熱伝導率λあるいは分散係数などを既知パラメータとして、地下水位、地下水温、物質の濃度などを計算する構造になっている。しかし、こうしたパラメータを直接精度よく測定することは難しく、むしろ解として得られるべき地下水位、地下水温、物質の濃度などの方が容易に、精度よく観測・測定が可能である。本発明では、地下水位や地下水温などの実測値に計算値が一致するように非線形最適化手法を用いて、各種パラメータを同定する逆解析を採用する。逆解析では、まず目的関数を設定し、それらを最小化することによって各パラメータを同定する。

本発明では、最適化手法で用いる目的関数Ｆとして、熱応答試験や揚水試験で得られた実測値と、解析解から得られた計算値の２乗誤差の総和を１／２乗したユークリッドノル
ムに、係数Ａ_１および係数Ａ_２を乗じて重み付けした次式（１７）を用いる。

また、最適化手法において、収束判定値として、ある探索点と、その１つ前の探索点に対して得られた目的関数Ｆの値の差分であるＤＩＦ＝Ｆ_ｊ−Ｆ_ｊ−１を用いる。更に、データ１つあたりの誤差を示す評価尺度には、次式（１８）で表すように、計算値と実測値の絶対値誤差の総和をデータ数Ｎで除した値ｆ／Ｎを用いる。

本発明では、目的関数の未知のパラメータを最適値に収束させるための、非線形最適化手法として、パウエル（Ｐｏｗｅｌｌ）の共役傾斜法を採用する。パウエル（Ｐｏｗｅｌｌ）の共役傾斜法は、他の多くの最適化手法で必要とされる目的関数の各パラメータに関する微係数を必要とせず、目的関数のみを用いて、その最適パラメータを求めることができ、高い汎用性を持つ。

図３．１に共役方向の概略図を示し、図３．２にパウエル（Ｐｏｗｅｌｌ）の共役傾斜法の基本アルゴリズムを示す。パウエル（Ｐｏｗｅｌｌ）の共役傾斜法では、２つの離れた点からｖの方向に正値２次形式の最小点を見つけ出すとき、各々で求めた最小点ｘ１，ｘ２を結ぶベクトルは、ｖと共役であるという性質を利用している。この性質を利用して、ｎ回直接探索を繰り返すことにより、最小点を見つけ出すことができる。

（熱応答試験の解析における同定対象のパラメータ）
熱応答試験において温度Ｔ、加熱量ｑ_Ｈ等は観測データとして既知である。式（８）の解析解において、等価体積熱容量（ρｃ）_ｅ、縦熱分散率κ_Ｌ、横熱分散率κ_Ｔ、真流速が未知の状態である。本発明では、同定対象のパラメータをこれらの４つのパラメータとして、式（１７）の目的関数Ｆを用いて、共役傾斜法を組み込んだ同定プログラムにより逆解析を行う。

本発明の逆解析では、式（１７）における係数Ａ_１の値は、次式（１９）に示すように、同定対象のパラメータであるκ_Ｔ、κ_Ｌの設定値から算出した値ａ_１の正負に応じて決定される。ａ_１が正の値になるのは、κ_Ｔがκ_Ｌよりも大きい場合である。しかしながら、Ｘ軸方向（縦方向）に地下水の流れがあることを想定している場合には、κ_Ｔがκ_Ｌよりも小さいはずである。つまり、係数Ａ_１は、同定対象のパラメータのうち、κ_Ｔとκ_Ｌがありえない大小関係となるように設定された場合に、１．０よりも大きな値となって目的関数Ｆを増大させるように設定される。

ａ_１＝（κ_Ｔ−κ_Ｌ）／κ_Ｌ≧０．０の場合、Ａ_１＝１．０＋ａ_１
ａ_１＝（κ_Ｔ−κ_Ｌ）／κ_Ｌ＜０．０の場合、Ａ_１＝１．０ （１９）

同様に、係数Ａ_２の値は、次式（２０）に示すように、同定対称のパラメータであるκ_Ｔ、（ρｃ）_ｅと、λ_ｅの設定値から算出した値ａ_２の正負に応じて決定される。ａ_２が正の値になるのは、κ_Ｔ＝α_Ｔｖ＋λ_ｅ／（ρｃ）_ｅと置くと、流速ｖが負の値になる場合であるが、流速ｖは正の値のはずである。つまり、係数Ａ_２は、同定対称のパラメータのうち、流速ｖがありえない範囲の値となるように設定された場合に、１．０よりも大きな値となり、目的関数Ｆを増大させるように設定される。

ａ_２＝１０（１／κ_Ｔ（ρｃ）_ｅ−１／λ_ｅ）≧０．０のとき、Ａ_２＝１．０＋ａ_２
ａ_２＝１０（１／κ_Ｔ（ρｃ）_ｅ−１／λ_ｅ）＜０．０のとき、Ａ_２＝１．０ （２０）

以上のように、係数Ａ_１および係数Ａ_２は、いずれも、パラメータの探索が間違った方向に向かった場合に、目的関数Ｆが増大する（ペナルティが加わる）ように設定されている。このような係数を設けることで、目的関数Ｆを最小値に収束させるまでの、最適化手法による探索回数を少なくすることが可能になる。

（熱応答試験の解析方法の揚水試験への適用）
熱応答試験と揚水試験は、両者とも熱伝導理論を基として、解析解を導いている。そのため、本発明の熱応答試験の解析方法および解析プログラムは、揚水試験の解析（本発明の参考例）にも適応することが可能である。揚水試験は現場においては、透水量係数Ｔ_ｍや貯留係数Ｓ_ｃといった地下水流動系パラメータを求めるための原位置試験として最も一般的に行われている。

いま、隣接帯水層から漏水があり、無限の広がりを持つ被圧帯水層の非定常流の地下水流動方程式は、地下水位の降下量をｓ＝Ｈ_０−ｈとし、隣接帯水層から漏水がある被圧帯水層へ完全貫入している井戸から一定量Ｑ_ｗで地下水が揚水されているものとすると、次式（２１）（２２）で与えられる。また、初期条件・境界条件は次式（２３）〜（２５）で与えられる。
但し、
ｈ：ピエゾ水頭（ｍ）
Ｈ_０：隣接帯水層のピエゾ水頭（ｍ）
ｒ：半径（ｍ）
Ｔ_ｍ：揚水している帯水層の透水量係数（ｍ^２／ｓ）
ｍ_ｃ：層厚
Ｋ：加圧層の透水係数
Ｓ_ｃ：貯留係数
Ｑ_ｗ：揚水量（Ｌ／ｓ）

上記の初期条件、上記の初期条件・境界条件のもとで式（２１）を解くと、次式（２６）（２７）の解析解が得られる。

式（２６）は、ハンタッシュ−ヤコブの式と呼ばれている。ｒ／λ＝Ｂと置くと、式（２６）は、次式（２８）で表されるハンタッシュ−ヤコブの井戸関数Ｗ（ｕ、Ｂ）を含む。

熱応答試験の解析で用いる、移流分散方程式から導出した式（８）の熱移動解析解においても、ハンタッシュ−ヤコブの井戸関数Ｗ（ｕ、Ｂ）が用いられている。すなわち、式（２６）の揚水試験に用いる解析解と、熱応答試験に用いる移流分散を考慮した解析解（式（８）は、式の形が良く似ている。そこで、熱応答試験の解析方法と同様に、揚水試験の解析でも、パウエル（Ｐｏｗｅｌｌ）の共役傾斜法を用いた逆解析を行う。

（揚水試験の解析における同定対象のパラメータ）
次に、揚水試験においては、式（２６）の解析解において、透水量係数Ｔ_ｍ、貯留係数Ｓ_ｃ、Ｂが未知の状態である。本発明では、求めるパラメータを次式（２９）〜（３１）に示すα１、α２、Ｂとして、共役傾斜法を組み込んだ同定プログラムにより逆解析を行う。そして、同定されたパラメータα１、α２、Ｂから、式（２９）〜（３１）より透水量係数Ｔ_ｍ、貯留係数Ｓ_ｃ、Ｂを算出する。
α１＝ｒ^２Ｓ_ｃ （２９）
α２＝Ｔ_ｍ （３０）
Ｂ＝ｒ／ｐ （３１）
但し、

以上のように、熱応答試験の解析解と揚水試験の解析解は、元を辿れば熱伝導理論によるものであり、パラメータは異なるが、式の形自体は同じである。これにより、本発明による熱応答試験の解析方法を揚水試験にも同様に適応できることが分かる。

本発明によれば、地下水流動による移流・分散の影響を考慮に入れた解析解を用いて解析を行うことができる。また、パウエルの共役傾斜法を適用した逆解析により、熱応答試
験の測定データから未知のパラメータを同定でき、従来の解析では求めることのできなかったパラメータについても同定できる。具体的には、地下水流動があるサイトでの熱応答試験の解析により、等価体積熱容量（ρｃ）_ｅ、縦熱分散率κ_Ｌ、横熱分散率κ_Ｔ、真流速などのパラメータを精度良く求めることができる。よって、地下水流動があるサイトでの熱交換井の熱交換挙動を従来よりも精度良く予測することが可能である。

飽和多孔体中を流体が流動している時の熱の流れを模式的に示す説明図である。 第１種変形ベッセル関数Ｉ_０（Ｂ）の正確解とその近似解（漸近解）である。 第２種変形ベッセル関数Ｋ_０（Ｂ）の正確解とその近似解（漸近解）である。 タイスの正確解Ｗ（ｕ）とその近似解（漸近解）である。 Ｗ（ｕ、Ｂ）の正確解と、その漸近解との比較結果である。 パウエルの共役傾斜法における共役方向を示す概念図である。 パウエルの共役傾斜法の基本アルゴリズムを示す説明図である。 温度回復データを用いた解析方法の検証結果（テストケース１）である。 温度回復データを用いた解析方法の検証結果（テストケース２）である。 熱応答試験結果（往き温度、還り温度、平均温度、加熱量）を示すグラフである。 熱伝導率λが異なる複数のテストケースの回復時の温度変化量の計算結果である。 ｒ^２（ρｃ）が異なる複数のテストケースの回復時の温度変化量の計算結果である。 熱交換井直下でのテストケース１〜３の温度変化量の計算結果である。 熱交換井直下での室内熱応答試験による計測値、および原位置熱応答試験による計測値を示すグラフである。 同定対象のパラメータの影響を検証する感度解析結果（座標１）である。 同定対象のパラメータの影響を検証する感度解析結果（座標２）である。 同定対象のパラメータの影響を検証する感度解析結果（座標３）である。 パラメータの探索過程（目的関数Ｆ）を示すグラフである。 パラメータの探索過程（収束判定値ＤＩＦ）を示すグラフである。 パラメータの探索過程（等価体積熱容量（ρｃ）_ｅ）を示すグラフである。 パラメータの探索過程（縦熱分散率κ_Ｌ）を示すグラフである。 パラメータの探索過程（横熱分散率κ_Ｔ）を示すグラフである。 パラメータの探索過程（真流速）を示すグラフである。 逆解析より求めたパラメータを用いて温度変化量を再現計算した値と正確値との比較結果である。 室内熱応答試験装置を模式的に示す全体構成図である。 土槽の説明図である。 実験地盤の中心部に埋設された温度センサーの配置図（ＸＹ面上）である。 実験地盤の中心部に埋設された温度センサーの配置図（浸透流に垂直な断面上：ＹＺ面上）である。 熱応答試験の計測結果（実験ケース１：Ｕ字管の出入口温度、平均温度、熱交換量）である。 熱応答試験の計測結果（実験ケース１：温度変化量の経時変化）である。 実験ケース１の同定結果を用いて演算した温度変化量の計算値と実測値との比較結果である。 熱応答試験の計測結果（実験ケース２：Ｕ字管の出入口温度、平均温度、熱交換量）である。 熱応答試験の計測結果（実験ケース２：温度変化量の経時変化）である。 実験ケース２の同定結果を用いて演算した温度変化量の計算値と実測値との比較結果である。 等価体積熱容量（ρｃ）_ｅの同定結果と簡易試験より得られた実測値である。 真流速の同定結果と簡易試験より得られた実測値である。 縦熱分散率κ_Ｌの同定結果と簡易試験より得られた実測値である。 横熱分散率κ_Ｔの同定結果と簡易試験より得られた実測値である。 縦分散長α_Ｌの同定結果である。 縦分散長α_Ｔの同定結果である。 逆解析より求めたパラメータを用いて水位降下量を再現計算した値と正確値との比較結果（テストケース１）である。 逆解析より求めたパラメータを用いて水位降下量を再現計算した値と正確値との比較結果（テストケース２）である。 揚水試験を実施したサイトの説明図である。 現場内の深井戸の掘削データをもとにした断面図である。 揚水井から２３．５ｍの距離（第一帯水層）にある観測井の地下水位の変動の計測値である。 同定結果を用いて揚水試験を再現計算したときの観測井における水位降下量の計算値と実測値との比較結果である。 揚水井から２４．５ｍの距離（第二帯水層）にある観測井の地下水位の変動の計測値である。 同定結果を用いて揚水試験を再現計算したときの観測井における水位降下量の計算値と実測値との比較結果である。

［温度回復データを用いる解析方法］
まず、移流・分散現象を考慮した、本発明の熱応答試験の解析方法の妥当性検証に先立って、解析に用いるデータ区間の検討について述べる。本発明者らは、以前に開発した解析プログラム（移流・分散現象を考慮していない熱移動方程式の解析解（Ｋｅｌｖｉｎの線源関数の解析解）に基づき、パウエルの共役傾斜法を用いた逆解析を行ってパラメータを同定するプログラム：特許第５３３４２２１号公報）では、熱応答試験において熱交換器に熱媒体を循環している間の計測値（循環時データ）を用いていた。ここでは、解析に用いるデータ区間を広げるため、循環停止後の温度回復データを用いることを検討する。

（順解析による検証）
以下の解析方法は、移流・分散現象を考慮していない、Ｋｅｌｖｉｎの線源関数の解析解を用いるものである。この場合、加熱停止後の温度変化量は、重ね合わせの原理により、次式（３２）〜（３４）で表される。そこで、Ｅ（ｘ）とＥ（ｘ´）をべき級数展開式の第３０項まで評価した高次漸近解を用いて、パウエルの共役傾斜法を適用したパラメータ同定による逆解析を行う解析プログラムを開発した。

式（３２）の高次漸近解を用いた解析プログラムの妥当性検証を以下の手順で行った。（１）架空の試験条件を用いて算出した温度変化量を検証データとする。
（２）同定対象のパラメータ（熱伝導率λ、ｒ^２（ρｃ））の任意の初期推定値を出発点として、これらのパラメータを逆解析によって同定する。
（３）同定したパラメータを用いて温度変化量の再現計算を行い、（１）で算出した検証データと比較することにより、その精度を確認する。

図４．１と図４．２は、温度回復データを用いた解析方法の検証結果であり、テストケース１、２について、逆解析より求めたパラメータを用いて温度変化量を再現計算した値と正確値を示す。また、表１．１と表１．２はパラメータの同定結果を示す一覧表である。このように、同定されたパラメータを用いた計算結果と正確値は非常に良く一致することが確認できた。また、ｆ／Ｎ値（データ１つ当たりの誤差）もほぼ０に等しく、精度の高い結果を得ることが出来た。

（原位置熱応答試験による検証）
本発明者らは、現場サイトで熱応答試験を行い、計測結果に対して逆解析を行い、熱伝導率λ、ｒ^２（ρｃ）を同定した。実施条件、測定機器などの詳細は省略する。７２時間かけて加熱を行い、その後に得られた地盤温度回復データを解析に用いた。図４．３は、熱応答試験結果（往き温度、還り温度、平均温度、加熱量）を示すグラフである。加熱量、流量ともに安定して推移している。加熱前の地層温度は約１４．１℃で、加熱開始から７２時間後には往き温度約２７．０℃、還り温度約２４．０℃となり、加熱停止後は緩やかに加熱開始前の温度へ回復している。

逆解析による同定結果を表１．３に示す。「全データ」は加熱時間と同程度の温度回復データを用いた場合、「前半部分」は温度回復初期のデータを用いた場合、「循環時法」は循環時データを用いた場合である。温度回復データの全データと温度回復初期データのどちらを用いた場合も、熱伝導率λは循環時法で得られた同定結果２．５１（Ｗ／ｍ／Ｋ）とほとんど等しいが、ｒ^２（ρｃ）は、温度回復初期データを用いた方が近い値となった。なお、ｆ／Ｎ値はそれぞれ０．１８４、０．１３３で、温度回復初期データを用いた結果の方が実測値との適合性が高いことが確認された。

また、本発明者らは、温度回復初期データを用いた場合に、解析解のパラメータの値（熱伝導率λ、ｒ^２（ρｃ））が温度変化量に与える影響を検討するため、パラメータの値が異なる複数のテストケースで計算を行った。表１．４は熱伝導率λを変化させたケース、表１．５はｒ^２（ρｃ）を変化させたケースである。

図４．４は表１．４のテストケース（熱伝導率λが異なる複数のケース）に対する回復時の温度変化量の計算結果である。熱伝導率λの値によってグラフの形状は大きく異なり，熱伝導率λが解析解へ与える影響が非常に強いことがわかる。一方、図４．５は表１．５のテストケース（ｒ^２（ρｃ）が異なる複数のケース）に対する回復時の温度変化量の計算結果である。ｒ^２（ρｃ）の値を変化させてもあまりグラフの形状は変わらず、ｒ^２（ρｃ）が解析解へ与える影響が小さいことが分かる。計算結果より、加熱時間と同程度の温度回復データを用いた逆解析では、ｒ^２（ρｃ）が解析解へ与える影響が小さくなり、循環時データを用いた場合の結果と大きく乖離したものと考えられる。また、温度回復初期データを用いた場合では、ｒ^２（ρｃ）の影響が考慮されるため、循環時データを用いた同定結果に近い値が得られたと考えられる。

以上のように、本発明者らは、熱応答試験の加熱停止後の温度回復データを用いる場合の逆解析方法を開発し、現場サイトで実施された熱応答試験データの解析を行った。その結果、非常に高精度でパラメータ同定を行えることが確認できた。また、熱応答試験の加熱停止後初期の温度回復データを用いることで、適正なｒ^２（ρｃ）を評価できることが確認できた。これは、移流・分散現象を考慮していない熱移動方程式を用いた解析の場合であるが、以下に説明する本発明の解析方法についても、同様に、重ね合わせの原理等を
用いて、加熱停止後の温度変化量を示す解析解を導き、温度回復データを用いる場合の逆解析方法を開発することで、パラメータ同定の精度を向上させることのできる可能性がある。

［移流・分散現象を考慮した熱応答試験の解析方法］
続いて、本発明の実施の形態に係る移流・分散現象を考慮した熱応答試験の解析方法の妥当性検証について述べる。本発明者らは、以下の１〜５の順で妥当性検証を行った。
１．順解析による、地下水流動の影響の検証
２．物質移動との比較
３．パラメータの感度解析
４．テストケースによるパラメータ同定精度の検証
５．室内熱応答試験データによるパラメータ同定精度の検証

１．順解析による、地下水流動の影響の検証
本発明者らは、本発明の移流・分散現象を考慮した熱移動方程式の解析解に基づいて温度上昇量を演算するプログラムをＦＯＲＴＲＡＮコードで開発した。そして、表２に示す３つのテストケース（テストケース１〜３）のパラメータを用いて、温度上昇量を演算する順解析を行った。テストケースのパラメータは、いずれも、飽和砂を用いた室内実験値および文献値を用いた。テストケース１〜３は、流速ｖの値だけがケースによって異なり、他のパラメータは全て同じである。また、縦分散長α_Ｌと横分散長α_Ｔの比が１０．０となるようにパラメータを設定した。

表３は、テストケース１〜３における等価熱拡散率κ_ｅと機構的分散α_Ｌ×ｖを示している。テストケース１から３の順に、流速ｖの値が大きくなるにつれて、機構的分散α_Ｌ×ｖの方が等価熱拡散率κ_ｅよりも大きくなっている。

演算結果により、以下の事項が確認された。
（１）テストケース１（流速ｖが最小）
等価熱拡散率κ_ｅの方が機構的分散α_Ｌ×ｖよりも大きいことを反映し、同心円状に熱が拡がる結果となった。熱の拡がり範囲は、計算時間１０００ｈｒでも、最大で５．０ｍとなった。
（２）テストケース２（流速ｖがテストケース１よりも１桁大きい）
機構的分散α_Ｌ×ｖの方が等価熱拡散率κ_ｅよりも大きくなっており、熱分散が支配的になったことを反映して、楕円状に熱が拡がる結果となった。また、移流効果も伴って、計算時間１０００ｈｒでは１０．０ｍまで熱が拡がった。
（３）テストケース３（流速ｖがテストケース２よりもさらに１桁大きい）
テストケース２よりも更に移流・分散効果が卓越し、より広域に熱が拡がっており、計算時間１０００ｈｒでは３０．０ｍまで熱が拡がった。

図５．１は、熱交換井直下（座標：ｘ＝１．０ｍ、ｙ＝０．０ｍ：熱注入位置から１ｍ下流の位置）におけるテストケース１〜３の温度変化量（温度上昇量）の計算結果を示すグラフである。この図に示すように、流速が小さいとき、熱伝導が支配的であるため、温度の立ち上がりは緩やかであるが、時間の経過に伴って温度は上昇し続ける。流速が大きくなるにつれて、加熱開始直後に急激に温度が変化し、時間の経過とともに温度変化が緩慢になる。これは、図２．４に示した、ハンタッシュ−ヤコブの標準曲線からもわかるように、流速が小さくなるとタイスの曲線に漸近し、流速が大きくなると早く定常状態に達するようになるためである。

ここで、図５．２に、図５．１と同じ座標（ｘ＝１．０ｍ、ｙ＝０．０ｍ）における室内熱応答試験による計測値、および原位置熱応答試験による計測値を示す。本発明者らは、後述するように、室内熱応答実験装置を用いて熱応答試験を実施しており、その結果、図５．１で説明した順解析の結果と同様の傾向がみられることが確認された。また、現場サイトの２地点で実施した熱応答試験の結果から、地下水流動が少ないとみられるサイトＢと、地下水流動が多いとみられるサイトＡでの計測値は、順解析と同様の傾向となっていることが確認された。

以上より、本発明の解析方法は、地下水流速の熱移動への影響を適切に表現できることが確認された。

２：物質移動との比較
上述したように、地下水流動場において、原点に連続的に注入したトレーサの濃度差分布は、熱移動と同じ解析理論によって解析できることから、熱移動方程式の解析解に基づいて温度上昇量を演算するプログラムは、物質輸送式の解析解についても適用できる。そ
こで、熱移動と物質移動との違いを比較検証するために、熱移動のテストケース１、２と同じ計算条件（表２参照）を用いて、濃度差分布を演算する順解析を行った。表４に、物質移動の演算に用いる２つのテストケースのパラメータを示す。また、表５に、各ケースの有効分子拡散係数Ｄ_ｅｆと機構的分散α_Ｌ×ｖを示す。表４では、熱移動と物質移動の演算結果を比較しやすくするために、ｑ_Ｈ／（ρｃ）_ｅ＝Ｍ／ｎとなるように、Ｍ（質量）の値を設定した。一般的に，溶質の有効分子拡散係数 は１０^５（ｃｍ^２／ｓ））のオーダーで非常に小さいため、物質移動では移流・分散現象が卓越することが多い。

物質移動の演算結果を熱移動の演算結果と比較すると、流速が小さい場合は、熱移動に比べて物質移動の方が注入点付近における変化が大きく、あまり広がらないことが確認された。これは、溶質の分子拡散係数が水の熱拡散率よりも２桁程度小さいため、熱移動の方が熱伝導によってより広域に移動するからである。逆に、流速が大きくなると、物質の方が熱よりも広域に移動している。これは、物質移動では流速ｖで移動するのに対して、熱移動では次式（３５）で示す速度（真流速）で熱が移動するため、物質移動の方がより流速の影響を顕著に受け、移流現象が卓越するためである。

３．パラメータの感度解析
本発明者らは、熱移動方程式の解析解を用いた逆解析において、同定対象のパラメータに設定した等価体積熱容量（ρｃ）_ｅ、縦熱分散率κ_Ｌ、横熱分散率κ_Ｔ、真流速が解析解に与える影響を検証する感度解析を行った。感度解析には、表６．１に示すパラメータの値を用いた。表６．１に示す３点の座標（ｘ＝２．０ｍ、ｙ＝１．０ｍ）、（ｘ＝２．０ｍ、ｙ＝０．０ｍ）、（ｘ＝０．０ｍ、ｙ＝１．０ｍ）において、同定対象の各パラメータを、表６．１に示す値（正確値）、正確値の１０倍、正確値の１／１０、に設定した場合の温度上昇量を演算した。

図６．１、図６．２、図６．３に各座標における感度解析結果を示す。図６．１に示すように、座標１（ｘ＝２．０ｍ、ｙ＝１．０ｍ）では、等価体積熱容量（ρｃ）_ｅが最も大きく解析解へ影響を与えており、パラメータを正確値（表６．１の値）とした場合の演算結果との乖離が最も大きい。次いで縦熱分散率κ_Ｌ、横熱分散率κ_Ｔ、真流速の順に、影響（正確値を用いた演算結果との乖離）が大きい。一方、図６．２に示すように、座標２（ｘ＝２．０ｍ、ｙ＝０．０ｍ）では、等価体積熱容量（ρｃ）_ｅが最も大きく解析解へ影響を与えており、次いで横熱分散率κ_Ｔ、縦熱分散率κ_Ｌ、真流速の順に影響が大きい。また、座標３（ｘ＝０．０ｍ、ｙ＝１．０ｍ）は流動方向と直交する点であるが、等価体積熱容量（ρｃ）_ｅ、縦熱分散率κ_Ｌ、横熱分散率κ_Ｔ、真流速の順に影響が大きい。

以上より、座標によって、各パラメータが解析解へ与える影響の傾向が異なることが確認できた。また、どの座標においても、等価体積熱容量（ρｃ）_ｅが最も解析解へ大きく影響することが確認できた。

４．テストケースによるパラメータ同定精度の検証
本発明者らは、本発明の移流・分散現象を考慮した熱移動方程式の解析解の高次漸近解を用いてパラメータの逆解析を行う解析プログラムを開発した。そして、これを用いた場
合のパラメータの同定精度の検証を、以下の手順で行った。
（１）表６．２に示す試験条件（テストケース）を与え、開発した解析プログラムを用いて温度変化量を算出し、これを検証データとする。
（２）同定対象の４つのパラメータ（等価体積熱容量（ρｃ）_ｅ、縦熱分散率κ_Ｌ、横熱分散率κ_Ｔ、真流速）の任意の初期推定値を出発点として、これらのパラメータを逆解析によって同定する。
（３）同定したパラメータを用いて温度変化量の再現計算を行い、（１）で算出した検証データと比較することにより、その精度を確認する。

図６．４〜６．１０はパラメータの探索過程を示すグラフであり、図６．４は目的関数Ｆ、図６．５は収束判定値ＤＩＦである。また、図６．６は等価体積熱容量（ρｃ）_ｅ、図６．７は縦熱分散率κ_Ｌ、図６．８は横熱分散率κ_Ｔ、図６．９は真流速である。また、図６．１０は逆解析より求めたパラメータを用いて温度変化量を再現計算した値と正確値（図６．１参照）との比較である。

正確値に対し、適当な初期推定値を与えて逆解析を行った。目的関数Ｆが最小値を目指すよう探索を繰り返し、これ以上値が変化せず、ＤＩＦ値が限りなく０に近い値を示したところで探索を終了した。図６．６〜６．１０から、探索を繰り返すごとに、パラメータの値が正確値に徐々に近づき、同定されたパラメータを用いた計算結果と正確値（図６.１０）は、非常に良く一致することが確認できた。表６．３は各パラメータの同定結果を
示す一覧表である。表６．３に示すように、ｆ／Ｎ値（データ１つ当たりの誤差）もほぼ０に等しく、精度の高い結果を得ることが出来た。以上の検証より、本発明の逆解析方法の妥当性を確認することが出来た。

５．室内熱応答試験データによるパラメータ同定精度の検証
本発明者らは、室内熱応答実験装置を用いて熱応答試験を実施し、得られたデータに本発明による逆解析方法を適用し、パウエルの共役傾斜法を用いたパラメータの同定プログラム（本発明を適用した熱応答試験の解析プログラム）を用いて解析を行った。以下に、室内熱応答実験装置の構成および試験の実施手順を説明する。

図６．１１は、室内熱応答試験装置を模式的に示す全体構成図である。室内熱応答試験装置１は、実験地盤１０と熱応答実験装置２０から構成される。実験地盤１０は、土槽を用いて、原位置熱応答試験を実地に行う際の地盤環境を模擬的に表現するためのものである。図６．１１に示すように、実験地盤１０は、土槽１１と、水位設定部１２と、浸透量調節部１３等を備える。浸透量調節部１３は、ジャッキ、敷石、滑り止めマット、土槽補強材等から構成される。

図６．１１に示すように、水位設定部１２は、土槽１１の上流側に配置された上流側水位調節部１２Ａと、土槽１１の下流側に配置された下流側水位調節部１２Ｂを備える。上流側水位調節部１２Ａは給水容器１４Ａおよび上流側水位調整容器１４Ｂを備え、下流側水位調節部１２Ｂは下流側水位調整容器１５を備える。上流・下流の水位調節部１２Ａ、１２Ｂの高さを調整することで、土槽１１内に所定の水位が維持される。

図６．１２は土槽１１の説明図である。土槽１１は、厚さ２０ｍｍの透明アクリル板製の箱型の容器であり、その内法寸法は、高さ４８０ｍｍ、幅９６０ｍｍ、奥行き４６０ｍｍである。図６．１２に示すように、土槽１１の内部には幅７００ｍｍの試料充填部１６が設けられ、試料充填部１６の幅方向Ｘの両側に貯水部１７Ａ、１７Ｂが設けられている。試料充填部１６と貯水部１７Ａ、１７Ｂは、アクリル板にパンチ穴を設け、さらに金網で被覆した透水性スリット１８により分離されている。土槽１１はその幅方向Ｘに傾斜するように配置され、土槽１１内では幅方向Ｘに浸透流が流れる。試料充填部１６の上流側に設けられた貯水部１７Ａの幅は１７５ｍｍであり、下流側に設けられた貯水部１７Ｂの幅は７５ｍｍである。

図６．１２において、実験時に想定される浸透流経路を模式的に示す。図６．１１に示すように、土槽１１内には、熱応答実験装置２０のＵ字管２１が配置されている。実験時の浸透量の調整を行うには、図６．１２に示すように土槽１１の上流側をジャッキアップさせて実験地盤１０を傾斜させて、上流側水位と下流側水位が同一深となるように調整す
る。これにより、Ｕ字管２１に直交する浸透流が発生する。

熱応答実験装置２０は、図６．１１に示すように、熱媒体（水）を加熱し、循環を行う循環系と、熱媒体温度・消費電力・循環流量などを計測する計測系を備える。循環系は、実験地盤１０に設置された熱交換器としてのＵ字管２１を含む循環流路２２、循環流路２２で熱媒体を循環させる循環ポンプ２３と、循環流路２２を通る熱媒体を加熱する電気ヒーター２４、循環流路２２に熱媒体を供給する給水容器２５等を備える。計測系は、Ｕ字管２１の入口・出口の温度を計測する温度センサー２６ａ、２６ｂ、実験地盤１０に埋設される複数の温度センサー２７、流量計２８、データロガー２９、パーソナルコンピューター３０、室温計３１等を備える。

図６．１３および図６．１４は、実験地盤１０の中心部に埋設された温度センサー２７の配置図である。図６．１２に示すように、本明細書では、Ｘ方向を土槽１１の幅方向とする。また、Ｙ方向（図６．１２の紙面に垂直な方向：図示省略）を土槽１１の奥行き方向とし、Ｚ方向を土槽１１の高さ方向とする。Ｘ方向、Ｙ方向、Ｚ方向は互いに直交する方向である。図６．１３は温度センサー２７のＸＹ面上の配置図であり、図６．１４はＵ字管２１を通り浸透流に垂直な断面（ＹＺ面）上の温度センサー２７の配置図である。

この構成の室内熱応答試験装置１では、電気ヒーター２４によって熱媒体を加熱し、循環ポンプ２３によって循環流路２２に循環させる。加熱された熱媒体は、試験対象の実験地盤１０に形成されているボアホール孔を通過する際に放熱するため、Ｕ字管２１の入口・出口の温度差を温度センサー２６ａ、２６ｂでそれぞれ測定し、実験地盤１０の熱交換特性を評価する。

パーソナルコンピューター３０には、キーボード、マウスなどの入力装置と、液晶表示器などのモニターと、プリンターなどの出力装置等が接続されている。パーソナルコンピューター３０には、本発明による熱応答試験の解析プログラム（以下の説明では、これを「同定プログラム」と呼ぶ場合もある。）がインストールされている。解析プログラムを実行することにより、データロガー２９から取り込んだ計測データに基づき、試験対象の地盤の熱交換特性の評価（パラメータの同定処理や、演算処理等）を行う。

本発明者らは、室内熱応答試験装置１を２４℃設定の恒温室内に設置して熱応答試験を行い、本発明の逆解析方法を適用してパラメータの同定を行った。表６．４に実験条件（実験ケース）を示す。本試験では、実験地盤１０の傾斜角度が異なる２つの実験ケース１，２を設定した。実験地盤１０に傾斜をつけることで地下水流れ場を発生させ、各実験ケースの浸透流量、浸透流が通過する断面積（２２０８ｃｍ^２）、間隙率（０．４６３）からダルシー流速、真流速を算出した。２つの実験ケースの単位長さ当たりの加熱量ｑ_Ｈはそれぞれ、１２２．９２５（Ｗ／ｍ）、１２２．６７５（Ｗ／ｍ）とした。逆解析時には，加熱量は常に一定であるという解析解の境界条件を満たすために、加熱開始より３０分以降のデータを使用した。

図６．１５、６．１６は実験ケース１に対する熱応答試験の計測結果である。図６．１
５はＵ字管２１の出入口温度、平均温度、熱交換量である。加熱開始前の地盤温度は約１７．４℃で安定しており、加熱開始から１２０分後には約３３．４９℃まで上昇した。熱交換量は、加熱開始から３０分まではあまり安定しないが、それ以降ではほとんど一定で推移している。また、図６．１６は２つの座標（座標１：ｘ＝０．０３ｍ、ｙ＝０．０６ｍ）、（座標２：ｘ＝０．０３ｍ、ｙ＝０．１２５ｍ）での温度変化量の経時変化を示す。なお，加熱時間は１２０分で、各座標における深さ方向の温度測定結果の平均値をとっている。このグラフより、熱源から近い座標１（ｘ＝０．０３ｍ、ｙ＝０．０６ｍ）の方が遠い座標２（ｘ＝０．０３ｍ、ｙ＝０．１２５ｍ）よりも温度変化量が大きくなることが分かる。

実験ケース１の場合の、逆解析によるパラメータの同定結果を表６．５に示す。また、図６．１７は、同定結果を用いて演算した温度変化量の計算値と実測値を示すグラフである。図６．１７より、どの座標においても計算値と実測値は非常に良く一致していることが確認できる．各座標における同定結果を比較すると、縦熱分散率κ_Ｌはほぼ等しいものの、横熱分散率κ_Ｔと真流速は、座標１（ｘ＝０．０３ｍ、ｙ＝０．０６ｍ）の方が座標２（ｘ＝０．０３ｍ、ｙ＝０．１２５ｍ）よりも小さく、等価体積熱容量（ρｃ）_ｅは座標１（ｘ＝０．０３ｍ、ｙ＝０．０６ｍ）の方が座標２（ｘ＝０．０３ｍ、ｙ＝０．１２５ｍ）よりも大きい。なお、逆解析に用いたデータは８８個で、データ１つ当たりの誤差を示すｆ／Ｎ値はそれぞれ６．３０×１０^−３（Ｋ）、４．２０×１０^−３（Ｋ）で非常に小さい。

次に、図６．１８、６．１９は実験ケース２に対する熱応答試験の計測結果であり、図６．１８はＵ字管２１の出入口温度、平均温度、熱交換量である。実験ケース２では、上記の２つの座標に加えて、座標３（ｘ＝０．０５ｍ、ｙ＝０．０５ｍ、ｚ＝０．１５ｍ）についても測定した。加熱開始前の地盤温度は約１７．１８℃で安定しており、加熱開始から１２０分後には約３３．１４℃まで上昇した。熱交換量は、実験ケース１と同様に、加熱開始から３０分以降ではほとんど一定で推移している。また、図６．１９は３つの座標１（ｘ＝０．０３ｍ、ｙ＝０．０６ｍ）、座標２（ｘ＝０．０３ｍ、ｙ＝０．１２５ｍ）、座標３（ｘ＝０．０５ｍ、ｙ＝０．０５ｍ）での温度変化量の経時変化を示す。なお、加熱時間は１２０分で、各座標における深さ方向の温度測定結果の平均値をとっている。この結果を実験ケース１の結果と比較すると、流速の大きいケース２の方が温度変化量が小さくなっている。これは、前章で述べたように、流速が大きくなると早く定常状態に達するようになるためである。

実験ケース２の場合の、逆解析によるパラメータの同定結果を表６．６に示す。また、
図６．２０は、同定結果を用いて演算した温度変化量の計算値と実測値を示すグラフである。図６．２０より、どの座標においても計算値と実測値は非常に良く一致していることが確認できる。各座標における同定結果を比較すると、縦熱分散率κ_Ｌはほとんど等しいものの、横熱分散率κ_Ｔは座標の関係もあり値がばらついているが、真流速は、どの座標でもほとんど等しい。また、等価体積熱容量（ρｃ）_ｅの同定値は、３つの座標で、それぞれ、１００５．０４（Ｗ・ｈｒ／ｍ^３／Ｋ）、９２４．７３１（Ｗ・ｈｒ／ｍ^３／Ｋ）、９９７．６７５（Ｗ・ｈｒ／ｍ^３／Ｋ）で、座標２（ｘ＝０．０３ｍ、ｙ＝０．１２５ｍ）における値が一番小さくなっている。なお、逆解析に用いたデータは８８個で、データ一つ当たりの誤差を示すｆ／Ｎ値はそれぞれ６．８０×１０^−３（Ｋ）、６．７０×１０^−３（Ｋ）２．３１×１０^−３（Ｋ）で非常に小さい。

ここで、逆解析によって得られた等価体積熱容量（ρｃ）_ｅの確からしさを検証するために，簡易試験を行って、実験的に等価体積熱容量（ρｃ）_ｅを求めた。簡易試験の手順は、概ね以下のようになる。断熱容器（魔法瓶）を用意し、これを、室内熱応答試験を行った恒温室内温度と同程度の水温、密度、比熱の脱気水で満たしておく。そして、室内熱応答試験装置１と同じく間隙率ｎ＝０．４６３となるように試料（硅砂４号）を充填した小型容器を加温し、これを断熱容器（魔法瓶）内に入れて密閉し、試料および水の初期温度、平衡温度等を計測する。計測結果から、次式（３６）に示す熱量保存則に基づいて等価体積熱容量（ρｃ）_ｅを算出する。

本発明者らは、上記の簡易試験を３回行った。表６．７に、試験結果と、得られた等価体積熱容量（ρｃ）_ｅの値を示す。３回の試験結果の平均値は、８６２．７９（Ｗ・ｈｒ／ｍ^３／Ｋ）となった。

図６．２１に、室内熱応答試験の実験ケース１、２における等価体積熱容量（ρｃ）_ｅの同定結果と簡易試験より得られた等価体積熱容量（ρｃ）_ｅの実測値を示す。簡易試験より得られた実測値よりも本発明の逆解析より得られた同定値の方がやや大きく評価されており、座標２（ｘ＝０．０３ｍ、ｙ＝０．１２５ｍ）における同定結果が最も実測値と近くなった。本来、等価体積熱容量（ρｃ）_ｅは地下水流動や座標等の影響を受けない固定の値を示すが、本逆解析方法より得られた同定結果はややばらついている。

また、本発明者らは、実験ケース１、２の座標２（ｘ＝０．０３ｍ、ｙ＝０．１２５ｍ）のデータに基づき、従来の解析方法（移流・分散現象を考慮しない解析方法：Ｋｅｌｖｉｎの線源関数の解析解を用いる方法）で逆解析を行い、パラメータを同定した。その結果、等価体積熱容量（ρｃ）_ｅとして、以下の値が得られた。
実験ケース１：１１２８．７２３（Ｗ・ｈｒ／ｍ^３／Ｋ）
実験ケース２：１２０９．８５（Ｗ・ｈｒ／ｍ^３／Ｋ）
また、ｆ／Ｎ値は、それぞれ以下のようになった。
実験ケース１：１．３×１０^−２
実験ケース２：１．２２×１０^−２

図６．２１に示す、等価体積熱容量（ρｃ）_ｅの実測値と本発明の解析方法による同定値を、従来の解析方法による同定値とを比較すると、本発明の解析方法では、従来の解析方法よりも実測値に近い同定値が得られ、従来の解析方法よりも格段に精度良くパラメータを同定できることが確認できた。また、ｆ／Ｎ値についても、本発明の解析方法では、従来の解析方法よりも小さくなることが確認できた。

図６．２２に、実験ケース１、２における真流速の同定結果と室内熱応答試験の結果より得られた真流速の実測値を示す。この図から、実測値よりも本発明の逆解析より得られた同定値の方がやや大きく評価されており、ケース毎で比較すると、実験ケース２の方が実験ケース１よりも実測値に近い値を示している。また、真流速においても等価体積熱容量（ρｃ）_ｅと同様に、座標２（ｘ＝０．０３ｍ、ｙ＝０．１２５ｍ）における同定結果が最も実測値と近くなった。

図６．２３、６．２４に、実験ケース１、２における縦熱分散率κ_Ｌ、横熱分散率κ_Ｔの同定結果を示す。どのケースにおいても縦熱分散率κ_Ｌの方が横熱分散率κ_Ｔより大きくなっている。したがって、一般的に縦熱分散率κ_Ｌに含まれる縦分散長α_Ｌの方が、横熱分散率κ_Ｔに含まれる横分散長α_Ｔよりも大きくなるという特徴が正確に評価されたといえる。縦熱分散率κ_Ｌは、同じ座標であれば地下水流速の大きい実験ケース２の方が実
験ケース１よりも大きくなる。図６．２３をみると、座標２（ｘ＝０．０３ｍ、ｙ＝０．１２５ｍ）では実験ケース２の方が実験ケース１よりも大きくなっているが、座標１（ｘ＝０．０３ｍ、ｙ＝０．０６ｍ）では実験ケース２の方が実験ケース１よりも小さくなっている。この原因は、微妙な実験誤差や同定誤差が原因であると考えられる。また、横熱分散係数も同様に、より流速の大きい実験ケース２の方が実験ケース１よりも大きくなるが、座標２（ｘ＝０．０３ｍ、ｙ＝０．１２５ｍ）ではその傾向を示しているが、座標１（ｘ＝０．０３ｍ、ｙ＝０．０６ｍ）では逆の傾向が見られる。この原因についても、微妙な実験誤差や同定誤差が原因であると考えられる。

図６．２５、６．２６に、各実験ケースにおける縦分散長α_Ｌ、横分散長α_Ｔの算出結果を示す。縦分散長α_Ｌ、横分散長α_Ｔは、パラメータの同定結果を用いて、本発明の解析理論で用いた式（３）より算出した。算出用の等価体積熱容量（ρｃ）_ｅ、真流速は同定値を用い、熱伝導率λは硅砂４号完全飽和の傾斜０％における同定結果（２．４６６Ｗ／ｍ／Ｋ）を用いた。一般的に、縦分散長α_Ｌと横分散長α_Ｔの関係はα_Ｔ≦α_Ｌとなることが知られており、図６．２４、６．２５より、各実験ケースにおける縦分散長α_Ｌ、横分散長α_Ｔを比較すると、この関係が満足されていることがわかる。

［揚水試験の解析方法への適用の妥当性検証］
次に、本発明の移流・分散現象を考慮した熱応答試験の解析方法を、揚水試験の解析方法に適用することの妥当性検証について述べる。本発明者らは、以下の６〜７の順で妥当性検証を行った。
６．テストケースによるパラメータ同定精度の検証
７．揚水試験データによるパラメータ同定精度の検証

６．テストケースによるパラメータ同定精度の検証
本発明者らは、本発明の移流・分散現象を考慮した地下水流動方程式の解析解の高次漸近解を用いてパラメータの逆解析を行う解析プログラムを開発した。そして、これを用いた場合のパラメータの同定精度の検証を、以下の手順で行った。
（１）表６．８に示す試験条件（テストケース）を与え、開発した解析プログラムを用いて水位降下量算出し、これを検証データとする。
（２）同定対象の３つのパラメータ（透水量係数Ｔ_ｍ、貯留係数Ｓ_ｃ、Ｂ）の任意の初期推定値を出発点として、これらのパラメータを逆解析によって同定する。
（３）同定したパラメータを用いて水位降下量の再現計算を行い、（１）で算出した検証データと比較することにより、その精度を確認する。

表６．８に示す２つの試験条件のうち、テストケース１に対する検証結果を表６．９および図６．２７に示す。表６．９は逆解析による同定結果であり、図６．２７は逆解析より求めたパラメータを用いて水位降下量を再現計算した値と正確値との比較である。テストケース１では、正確値に対して初期推定値を２倍の大きさにして解析を行った。図６.２７より、計算値と正確値は非常に良く一致している。また。表６.９に示すようにｆ／Ｎ値（データ１つ当たりの誤差）もほぼ０に等しく、精度の高い結果を得ることが出来た。

テストケース２に対する検証結果を表６．１０および図６．２８に示す。表６．１０は逆解析による同定結果であり、図６．２８は逆解析より求めたパラメータを用いて水位降下量を再現計算した値と正確値との比較である。テストケース２では、設定値に対して初期推定値を１／２倍の大きさにして解析を行った。図６.２８より、計算値と設定値は非常に良く一致しており、表６.１０に示すようにｆ／Ｎ値（データ１つ当たりの誤差）もほぼ０に等しく、テストケース１と同様に精度の高い結果を得ることが出来た。以上より、本逆解析方法の妥当性を確認することが出来た。

７．揚水試験データによるパラメータ同定精度の検証
本発明者らは、現場サイトで実施した揚水試験のデータに本発明による逆解析方法を適用し、パウエルの共役傾斜法を用いたパラメータの同定プログラム（本発明の参考例に係る揚水試験の解析プログラム）を用いて解析を行った。以下に、揚水試験を実施した現場の地質条件や試験条件等を説明する。

（地質条件）
図６．２９に揚水試験を実施したサイトの位置を示す。揚水試験を実施したサイトは信州大学工学部キャンパス内にあり、犀川から約７６１ｍ離れた扇状地扇端部に位置する。現場内で掘削された複数の井戸から、サイトの帯水層構造についてはおおむね明らかになっている。図６.３０は現場内の深井戸の掘削データをもとにした断面図である。調査地付近の地質は主に砂礫層から構成されており、ところどころに粘土層・粘土混り砂層砂礫層を挟んでいる。このうち砂礫層が帯水層となり、粘土層・粘土混り砂層砂礫層が不（難）透水層となっている。図６.３０に示す電気検層結果より、同じ帯水層内にも異なる比抵抗値を示す部分があり、地盤の透水性に不均一性が見られる。

（試験条件）
上記のサイトにおいて、第１帯水層と第２帯水層の井戸能力を把握するため、層別揚水試験を行った。表６.１１に各帯水層における揚水時間、揚水量、帯水層厚、揚水井と観測井との距離について示す。揚水は観測井の水位変化が見られなくなるまで行い、それぞれ１８０ｍｉｎと１７５ｍｉｎとなった。揚水量はどちらも４１７（Ｌ／ｍｉｎ）である。また、帯水層厚は第１帯水層が３２．７ｍ、第２帯水層が２１．２ｍである。

揚水試験結果とパラメータ同定結果は、以下のようになった。
（１）第一帯水層
揚水井から２３．５ｍの距離にある観測井の地下水位の変動を図６.３１に示す。地下水位は試験開始１ｍｉｎ後から減少を続け、１２０ｍｉｎで０．２１５の水位降下が確認されている。１２０ｍｉｎ以降は水位の変化がほとんど見られなくなったため、１８０ｍｉｎで試験を終了している。

表６.１２に逆解析による同定結果を、図６.３２に同定結果を用いて揚水試験を再現計算したときの観測井における水位降下量の計算値と実測値を示す。逆解析より得られた透水量係数Ｔ_ｍは１．０２（ｍ^２／ｍｉｎ）、貯留係数Ｓ_ｃは３．５９×１０^−４、Ｂは０．０３４で、この値を用いた計算値と実測値は非常によく一致している。透水量係数Ｔ_ｍを帯水層厚３２．７ｍで除して、透水係数Ｋは０．０３１２（ｍ／ｍｉｎ）となった。この値は一般的な透水係数Ｋの値である。また、ｆ／Ｎ値は１．６２×１０^−３（ｍ）と小さく、精度の高い結果が得られた。なお、逆解析時に使用したデータは１８０個である。

（２）第二帯水層
揚水井から２４．５ｍの距離にある観測井の地下水位の変動を図６．３３に示す。地下水位は試験開始１ｍｉｎ後から減少を続け、１４０ｍｉｎで０．５８７ｍの水位降下が確認されている。１４０ｍｉｎ以降は水位の変化がほとんど見られなくなったため、１７５ｍｉｎで試験を終了している。

表６．１３に逆解析による同定結果を、図６．３４に同定結果を用いて揚水試験を再現計算したときの観測井における水位降下量の計算値と実測値を示す。逆解析より得られた透水量係数Ｔ_ｍは０．４５９（ｍ^２／ｍｉｎ）、貯留係数Ｓ_ｃは３．１８×１０^−５、Ｂは４．８４×１０^−５で、この値を用いた計算値と実測値は非常によく一致している。透水量係数Ｔ_ｍを帯水層厚２１．２ｍで除して、透水係数Ｋは０．０２１６（ｍ／ｍｉｎ）となった。この値は一般的な透水係数Ｋの値である。また、ｆ／Ｎ値は５．３２×１０^−３（ｍ）と小さく、第一帯水層と同様に精度の高い結果が得られた。なお、逆解析時に使用したデータは１７５個である。以上のように、隣接帯水層からの漏水がある場合でも、精度よく解析を行えることが確認できた。

［検証により得られた知見］
以上のように、本発明者らは、地下水流動場における熱応答試験の解析方法として、熱伝導に加えて移流・分散現象を考慮した解析解を導出し、パウエル（Ｐｏｗｅｌｌ）の共役傾斜法を用いた逆解析方法を開発してパラメータ同定を行った。また、この解析方法を揚水試験の解析方法にも適用し、パラメータ同定を行った。その結果と過程により得られた知見は以下の通りである。

（１）開発した解析解を用いて地下水流速の異なる３ケースで順解析を実施したところ、流速が小さく熱伝導現象が卓越するケースでは熱が円形に広がるのに対して、流速が大きくなるにつれて熱分散が卓越し、温度分布が楕円形になった。すなわち、本発明の解析方法は、地下水流速の熱移動への影響を適切に表現できる。
（２）同一条件下で熱移動と物質移動を比較すると、流速の小さいケースでは、熱伝導によって熱の方がより広範囲に移動するのに対して、流速が大きくなると地下水流動の影響を強く受ける物質の方が熱よりもより広域に移動する。
（３）開発した解析プログラムを室内熱応答試験データの解析に適用した結果、同定された等価体積熱容量（ρｃ）_ｅ、真流速は実測値と近い値が得られた。また、従来のＫｅｌｖｉｎの線源関数の解析解に基づく解析方法で同定した場合よりも、格段に精度良くパラメータを同定できることが確認できた。従って、地下水流動があるサイトでの熱交換井の熱交換挙動を従来よりも精度良く予測することが可能になると期待できる。
（４）漏水性を伴う帯水層への揚水試験解析解（ハンタッシュ−ヤコブの式）を用いた逆解析方法を開発し、現場サイトで実施された揚水試験へ適用した。その結果、非常に精度の高いパラメータ同定を行えることが確認できた。従って、地下水流動があるサイトでの地下水位の挙動を従来よりも精度良く予測することが可能になると期待できる。

１…室内熱応答試験装置
１０…実験地盤
１１…土槽
１２…水位設定部
１２Ａ…上流側水位調節部
１２Ｂ…下流側水位調節部
１３…浸透量調節部
１４Ａ…給水容器
１４Ｂ…上流側水位調整容器
１５…下流側水位調整容器
１６…試料充填部
１７Ａ、１７Ｂ…貯水部
１８…透水性スリット
２０…熱応答実験装置
２１…Ｕ字管
２２…循環流路
２３…循環ポンプ
２４…電気ヒーター
２５…給水容器
２６ａ、２６ｂ…温度センサー
２７…温度センサー
２８…流量計
２９…データロガー
３０…パーソナルコンピューター
３１…室温計

Claims (6)

1. 熱応答試験によって測定された調査対象の地盤における温度上昇量の経時変化に基づき前記地盤の熱交換特性を評価する熱応答試験の解析方法であって、
   前記地盤の温度上昇式を、地下水が流動している流動場での移流・分散現象を考慮した熱移動方程式の解析解である式（Ａ）で規定し、
   
   パウエルの共役傾斜法に基づく非線形最適化法を用いて、前記式（Ａ）で規定する解析解に含まれる未知のパラメータを所定の値に設定した場合に得られる温度上昇量の計算値と、前記熱応答試験によって測定された温度上昇量の実測値との誤差が最小となるように、前記パラメータを同定する逆解析を行い、
   当該逆解析においては、同定対象の前記パラメータを、等価体積熱容量（ρｃ）_ｅ、縦熱分散率κ_Ｌ、横熱分散率κ_Ｔ、真流速のうちの少なくとも１つとし、
   前記逆解析に用いる目的関数として、前記計算値と前記実測値との２乗誤差の総和を１／２乗したユークリッドノルムに、係数Ａ _１ および係数Ａ _２ を乗じて重み付けした式（Ｂ）を使用し、
   前記係数Ａ _１ の値を、次式（Ｃ１）に示すように、前記縦熱分散率κ _Ｌ と前記横熱分散率κ _Ｔ から算出した値ａ _１ に応じて決定し、
   
   ａ _１ ＝（κ _Ｔ −κ _Ｌ ）／κ _Ｌ ≧０．０の場合、Ａ _１ ＝１．０＋ａ _１
   ａ _１ ＝（κ _Ｔ −κ _Ｌ ）／κ _Ｌ ＜０．０の場合、Ａ _１ ＝１．０ （Ｃ１）
   
   前記係数Ａ _２ の値を、次式（Ｃ２）に示すように、前記横熱分散率κ _Ｔ と前記地盤の等価熱伝導率λ _ｅ から算出した値ａ _２ に応じて決定することを特徴とする熱応答試験の解析方法。
   
   ａ _２ ＝１０（１／κ _Ｔ （ρｃ） _ｅ −１／λ _ｅ ）≧０．０のとき、Ａ _２ ＝１．０＋ａ _２
   ａ _２ ＝１０（１／κ _Ｔ （ρｃ） _ｅ −１／λ _ｅ ）＜０．０のとき、Ａ _２ ＝１．０ （Ｃ２）
2. 請求項１において、
   式（Ａ）で規定する前記解析解において、次式（Ｄ）のように変数を置き換えた場合に前記解析解に含まれるハンタッシュ−ヤコブの井戸関数Ｗ（ｕ、Ｂ）を、次式（Ｅ）で規定する近似式で近似し、
   当該近似式に含まれる第１種変形ベッセル関数Ｉ _０ （Ｂ）、第２種変形ベッセル関数Ｋ _０ （Ｂ）、および、Ｂ＝０とした場合に前記ハンタッシュ−ヤコブの井戸関数Ｗ（ｕ、Ｂ）と一致する関数となるタイスの井戸関数Ｗ（ｕ）について、それぞれのべき級数展開式を少なくとも第３０項まで含む漸近解を使用して、前記地盤の温度上昇量を算出することを特徴とする熱応答試験の解析方法。
3. 請求項２において、
   前記漸近解として、第１種変形ベッセル関数Ｉ _０ （Ｂ）、第２種変形ベッセル関数Ｋ _０ （Ｂ）のべき級数展開式を第５０項まで含む漸近解を使用して、前記地盤の温度上昇量を算出することを特徴とする熱応答試験の解析方法。
4. 熱応答試験によって測定された調査対象の地盤における温度上昇量の経時変化に基づき前記地盤の熱交換特性を評価するために用いる熱応答試験の解析プログラムであって、
   コンピュータを、地下水が流動している流動場での移流・分散現象を考慮した熱移動方程式の解析解である式（Ａ）を用いて、地盤の温度上昇量を演算する演算手段として機能させ、
   コンピュータを、パウエルの共役傾斜法に基づく非線形最適化法を用いて、前記解析解に含まれる未知のパラメータを所定の値に設定した場合に得られる温度上昇量の計算値と、前記熱応答試験によって測定された温度上昇量の実測値との誤差が最小となるように、前記パラメータの逆解析を行う逆解析手段として機能させ、
   当該逆解析手段は、前記パラメータとして、等価体積熱容量（ρｃ） _ｅ 、縦熱分散率κ _Ｌ 、横熱分散率κ _Ｔ 、真流速のうちの少なくとも１つを同定する演算を行い、
   前記逆解析手段は、前記逆解析に用いる目的関数として、前記計算値と前記実測値との２乗誤差の総和を１／２乗したユークリッドノルムに、係数Ａ _１ および係数Ａ _２ を乗じて重み付けした式（Ｂ）を使用し、
   前記係数Ａ _１ の値を、次式（Ｃ１）に示すように、前記縦熱分散率κ _Ｌ と前記横熱分散率κ _Ｔ から算出した値ａ _１ に応じて決定し、
   
   ａ _１ ＝（κ _Ｔ −κ _Ｌ ）／κ _Ｌ ≧０．０の場合、Ａ _１ ＝１．０＋ａ _１
   ａ _１ ＝（κ _Ｔ −κ _Ｌ ）／κ _Ｌ ＜０．０の場合、Ａ _１ ＝１．０ （Ｃ１）
   
   前記係数Ａ _２ の値を、次式（Ｃ３）、（Ｃ４）に示すように、前記横熱分散率κ _Ｔ と前記地盤の等価熱伝導率λ _ｅ から算出した値ａ _２ に応じて決定することを特徴とする熱応答試験の解析プログラム。
   
   ａ _２ ＝１０（１／κ _Ｔ （ρｃ） _ｅ −１／λ _ｅ ）≧０．０の場合、Ａ _２ ＝１．０＋ａ _２
   ａ _２ ＝１０（１／κ _Ｔ （ρｃ） _ｅ −１／λ _ｅ ）＜０．０の場合、Ａ _２ ＝１．０ （Ｃ２）
5. 請求項４において、
   式（Ａ）で規定する前記解析解において、次式（Ｄ）のようにパラメータを置き換えた場合に前記解析解に含まれるハンタッシュ−ヤコブの井戸関数Ｗ（ｕ、Ｂ）を、式（Ｅ）で規定する近似式で近似し、
   前記演算手段は、当該近似式に含まれる第１種変形ベッセル関数Ｉ _０ （Ｂ）、第２種変形ベッセル関数Ｋ _０ （Ｂ）、および、Ｂ＝０とした場合に前記ハンタッシュ−ヤコブの井戸関数Ｗ（ｕ、Ｂ）と一致する関数となるタイスの井戸関数Ｗ（ｕ）について、それぞれのべき級数展開式を少なくとも第３０項まで含む漸近解を使用して、前記地盤の温度上昇量を演算することを特徴とする熱応答試験の解析プログラム。
6. 請求項５において、
   前記演算手段は、前記漸近解として、第１種変形ベッセル関数Ｉ _０ （Ｂ）、第２種変形ベッセル関数Ｋ _０ （Ｂ）のべき級数展開式を第５０項まで含む漸近解を使用して、前記地盤の温度上昇量を算出することを特徴とする熱応答試験の解析プログラム。