TOP > 国内特許検索 > 熱応答試験および揚水試験の解析方法および解析プログラム > 明細書

明細書 :熱応答試験および揚水試験の解析方法および解析プログラム

発行国 日本国特許庁(JP)
公報種別 公開特許公報(A)
公開番号 特開2016-017773 (P2016-017773A)
公開日 平成28年2月1日(2016.2.1)
発明の名称または考案の名称 熱応答試験および揚水試験の解析方法および解析プログラム
国際特許分類 G01N  25/18        (2006.01)
G01V   9/02        (2006.01)
FI G01N 25/18 L
G01V 9/02
請求項の数または発明の数 12
出願形態 OL
全頁数 49
出願番号 特願2014-139050 (P2014-139050)
出願日 平成26年7月4日(2014.7.4)
発明者または考案者 【氏名】藤縄 克之
出願人 【識別番号】504180239
【氏名又は名称】国立大学法人信州大学
個別代理人の代理人 【識別番号】100090170、【弁理士】、【氏名又は名称】横沢 志郎
【識別番号】100142619、【弁理士】、【氏名又は名称】河合 徹
【識別番号】100153316、【弁理士】、【氏名又は名称】河口 伸子
審査請求 未請求
テーマコード 2G040
2G105
Fターム 2G040AA01
2G040BA02
2G040BA03
2G040BA24
2G040CA02
2G040CB02
2G040CB16
2G040HA16
2G040HA18
2G040ZA05
2G105AA02
2G105LL07
要約 【課題】移流・分散現象を考慮した熱応答試験の解析方法および揚水試験の解析方法を提案すること。
【解決手段】熱応答試験の解析を行うにあたって、地盤の温度上昇式として、地下水が流動している流動場での移流・分散現象を考慮した熱移動方程式の解析解を用いる。そして、この解析解の近似式を用いて求められる地盤の温度上昇値の計算値と、熱応答試験で得られる地盤の温度上昇値の実測値との誤差が最小となるように、解析解に含まれる未知のパラメータ(等価体積熱容量(ρc)、縦熱分散率κ、横熱分散率κ、真流速)を同定する逆解析を行う。逆解析には、パウエルの共役傾斜法に基づく非線形最適化法を用いる。揚水試験の解析においても同様の解析方法を適用し、未知のパラメータを同定し、透水量係数T、貯留係数S、Bを算出する。
【選択図】図5.1
特許請求の範囲 【請求項1】
熱応答試験によって測定された調査対象の地盤における温度上昇量の経時変化に基づき前記地盤の熱交換特性を評価する熱応答試験の解析方法であって、
前記地盤の温度上昇式を、地下水が流動している流動場での移流・分散現象を考慮した熱移動方程式の解析解である式(A)で規定し、
JP2016017773A_000043t.gif パウエルの共役傾斜法に基づく非線形最適化法を用いて、前記式(A)で規定する解析解に含まれる未知のパラメータを所定の値に設定した場合に得られる温度上昇量の計算値と、前記熱応答試験によって測定された温度上昇量の実測値との誤差が最小となるように、前記パラメータを同定する逆解析を行い、
当該逆解析においては、同定対象の前記パラメータを、等価体積熱容量(ρc)、縦熱分散率κ、横熱分散率κ、真流速のうちの少なくとも1つとすることを特徴とする熱応答試験の解析方法。
【請求項2】
請求項1において、
前記逆解析に用いる目的関数として、前記計算値と前記実測値との2乗誤差の総和を1/2乗したユークリッドノルムに、係数Aおよび係数Aを乗じて重み付けした式(B)を使用し、
JP2016017773A_000044t.gif 前記係数Aの値を、次式(C1)に示すように、前記縦熱分散率κと前記横熱分散率κから算出した値aに応じて決定し、

=(κ-κ)/κ≧0.0の場合、A=1.0+a
=(κ-κ)/κ<0.0の場合、A=1.0 (C1)

前記係数Aの値を、次式(C2)に示すように、前記横熱分散率κと前記地盤の等価熱伝導率λから算出した値aに応じて決定することを特徴とする熱応答試験の解析方法。

=10(1/κ(ρc)-1/λ)≧0.0のとき、A=1.0+a
=10(1/κ(ρc)-1/λ)<0.0のとき、A=1.0 (C2)
【請求項3】
請求項1または2において、
式(A)で規定する前記解析解において、次式(D)のように変数を置き換えた場合に前記解析解に含まれるハンタッシュ-ヤコブの井戸関数W(u、B)を、次式(E)で規定する近似式で近似し、
JP2016017773A_000045t.gif 当該近似式に含まれる第1種変形ベッセル関数I(B)、第2種変形ベッセル関数K(B)、および、B=0とした場合に前記ハンタッシュ-ヤコブの井戸関数W(u、B)と一致する関数となるタイスの井戸関数W(u)について、それぞれのべき級数展開式を少なくとも第30項まで含む漸近解を使用して、前記地盤の温度上昇量を算出することを特徴とする熱応答試験の解析方法。
【請求項4】
請求項3において、
前記漸近解として、第1種変形ベッセル関数I(B)、第2種変形ベッセル関数K(B)のべき級数展開式を第50項まで含む漸近解を使用して、前記地盤の温度上昇量を算出することを特徴とする熱応答試験の解析方法。
【請求項5】
熱応答試験によって測定された調査対象の地盤における温度上昇量の経時変化に基づき前記地盤の熱交換特性を評価するために用いる熱応答試験の解析プログラムであって、
コンピュータを、地下水が流動している流動場での移流・分散現象を考慮した熱移動方程式の解析解である式(A)を用いて、地盤の温度上昇量を演算する演算手段として機能させ、
JP2016017773A_000046t.gif コンピュータを、パウエルの共役傾斜法に基づく非線形最適化法を用いて、前記解析解に含まれる未知のパラメータを所定の値に設定した場合に得られる温度上昇量の計算値と、前記熱応答試験によって測定された温度上昇量の実測値との誤差が最小となるように、前記パラメータの逆解析を行う逆解析手段として機能させ、
当該逆解析手段は、前記パラメータとして、等価体積熱容量(ρc)、縦熱分散率κ、横熱分散率κ、真流速のうちの少なくとも1つを同定する演算を行うことを特徴とする熱応答試験の解析プログラム。
【請求項6】
請求項5において、
前記逆解析手段は、前記逆解析に用いる目的関数として、前記計算値と前記実測値との2乗誤差の総和を1/2乗したユークリッドノルムに、係数Aおよび係数Aを乗じて重み付けした式(B)を使用し、
JP2016017773A_000047t.gif 前記係数Aの値を、次式(C1)に示すように、前記縦熱分散率κと前記横熱分散率κから算出した値aに応じて決定し、

=(κ-κ)/κ≧0.0の場合、A=1.0+a
=(κ-κ)/κ<0.0の場合、A=1.0 (C1)

前記係数Aの値を、次式(C3)、(C4)に示すように、前記横熱分散率κと前記地盤の等価熱伝導率λから算出した値aに応じて決定することを特徴とする熱応答試験の解析プログラム。

=10(1/κ(ρc)-1/λ)≧0.0の場合、A=1.0+a
=10(1/κ(ρc)-1/λ)<0.0の場合、A=1.0 (C2)
【請求項7】
請求項5または6において、
式(A)で規定する前記解析解において、次式(D)のようにパラメータを置き換えた場合に前記解析解に含まれるハンタッシュ-ヤコブの井戸関数W(u、B)を、式(E)で規定する近似式で近似し、
JP2016017773A_000048t.gif 前記演算手段は、当該近似式に含まれる第1種変形ベッセル関数I(B)、第2種変形ベッセル関数K(B)、および、B=0とした場合に前記ハンタッシュ-ヤコブの井戸関数W(u、B)と一致する関数となるタイスの井戸関数W(u)について、それぞれのべき級数展開式を少なくとも第30項まで含む漸近解を使用して、前記地盤の温度上昇量を演算することを特徴とする熱応答試験の解析プログラム。
【請求項8】
請求項7において、
前記演算手段は、前記漸近解として、第1種変形ベッセル関数I(B)、第2種変形ベッセル関数K(B)のべき級数展開式を第50項まで含む漸近解を使用して、前記地盤の温度上昇量を算出することを特徴とする熱応答試験の解析プログラム。
【請求項9】
揚水試験によって測定された調査対象の地盤における地下水位低下量の経時変化に基づき前記地盤の地下水流動特性を評価する揚水試験の解析方法であって、
前記地盤の地下水位低下量を、次式(F)で表されるハンタッシュ-ヤコブの式で規定し、
JP2016017773A_000049t.gif但し、
s:地下水位低下量
r:半径(m)
:揚水量(L/s)
:透水量係数(m/s)
:貯留係数

パウエルの共役傾斜法に基づく非線形最適化法を用いて、前記ハンタッシュ-ヤコブの式に含まれる未知のパラメータを所定の値に設定した場合に得られる地下水位低下量の計算値と、前記揚水試験によって測定された地下水位低下量の実測値との誤差が最小となるように、前記パラメータを同定する逆解析を行い、
当該逆解析においては、同定対象の前記パラメータを、前記ハンタッシュ-ヤコブの式における未知量を用いて、次式(G1)、(G2)、(G3)で表されるα1、α2、およびBとし、
α1=r (G1)
α2=T (G2)
B=r/p (G3)
但し、
JP2016017773A_000050t.gif K:加圧層の透水係数
:層厚

前記逆解析によって同定された前記パラメータα1、α2、およびBから、前記地盤の透水量係数T、貯留係数S、およびBのうち、少なくとも透水量係数Tを算出することを特徴とする揚水試験の解析方法。
【請求項10】
請求項9において、
前記式(G3)のようにパラメータを置き換えた場合に前記ハンタッシュ-ヤコブの式に含まれるハンタッシュ-ヤコブの井戸関数W(u、B)を、次式(E)で規定する近似式で近似し、
JP2016017773A_000051t.gif 当該近似式に含まれる第1種変形ベッセル関数I(B)、第2種変形ベッセル関数K(B)、および、B=0とした場合に前記ハンタッシュ-ヤコブの井戸関数W(u、B)と一致する関数となるタイスの井戸関数W(u)について、それぞれのべき級数展開式を少なくとも第30項まで含む漸近解を使用して、前記地盤の地下水位低下量を算出することを特徴とする揚水試験の解析方法。
【請求項11】
揚水試験によって測定された調査対象の地盤における地下水位低下量の経時変化に基づき前記地盤の地下水流動特性を評価するために用いる揚水試験の解析プログラムであって、
コンピュータを、次式(F)で表されるハンタッシュ-ヤコブの式を用いて、地盤の地下水位低下量を演算する演算手段として機能させ、
JP2016017773A_000052t.gif但し、
s:地下水位低下量
r:半径(m)
:揚水量(L/s)
:透水量係数(m/s)
:貯留係数

コンピュータを、パウエルの共役傾斜法に基づく非線形最適化法を用いて、前記ハンタッシュ-ヤコブの式に含まれる未知のパラメータを所定の値にそれぞれ設定した場合に得られる地下水位低下量の計算値と、前記揚水試験によって測定された地下水位低下量の実測値との誤差が最小となるように、前記パラメータを同定する逆解析手段として機能させ、
当該逆解析手段は、同定対象の前記パラメータを、前記ハンタッシュ-ヤコブの式における未知量を用いて、次式(G1)、(G2)、(G3)で表されるα1、α2、およびBとし、
α1=r (G1)
α2=T (G2)
B=r/p (G3)
但し、
JP2016017773A_000053t.gif K:加圧層の透水係数
:層厚

さらに、コンピュータを、前記逆解析によって同定された前記パラメータα1、α2、およびBから、前記地盤の透水量係数T、貯留係数S、およびBのうち、少なくとも透水量係数Tを算出する算出手段として機能させることを特徴とする揚水試験の解析プログラム。
【請求項12】
請求項11において、
前記式(G3)のようにパラメータを置き換えた場合に前記ハンタッシュ-ヤコブの式に含まれるハンタッシュ-ヤコブの井戸関数W(u、B)を、次式(E)で規定する近似式で近似し、
JP2016017773A_000054t.gif 前記演算手段は、当該近似式に含まれる第1種変形ベッセル関数I(B)、第2種変形ベッセル関数K(B)、および、B=0とした場合に前記ハンタッシュ-ヤコブの井戸関数W(u、B)と一致する関数となるタイスの井戸関数W(u)について、それぞれのべき級数展開式を少なくとも第30項まで含む漸近解を使用して、前記地盤の地下水位低下量を演算することを特徴とする熱応答試験の解析プログラム。
発明の詳細な説明 【技術分野】
【0001】
本発明は熱応答試験および揚水試験の解析方法および解析プログラムに関する。さらに詳しくは、地下水流動による移流・分散の影響を考慮に入れた解析を行うことができ、パウエルの共役傾斜法に基づく非線形最適化手法を用いて未知のパラメータを同定することによって、短時間で精度良く地盤の熱交換特性および地下水流動特性を評価できる解析方法および解析プログラムに関する。
【背景技術】
【0002】
近年、原子力や火力に代わるクリーンなエネルギーとして、自然エネルギー(再生可能エネルギー)の利用が見直されてきている。その1つとして、地中熱の利用に関心が高まってきている。季節によって変動する外気温に対して、地中の温度は年間を通して約15℃で安定している。このような特性を持つ地中の熱を有効利用したものとして、地中熱ヒート・ポンプ(Geothermal Heat Pump、以下GeoHPと略す)システムがある。これは、地中熱を利用した冷暖房・給湯システムのことである。安定した熱の供給が行なえる地中熱の特性に着目し、夏場には冷房として地中に熱を放出し、冬場には暖房の熱源として利用する。
【0003】
しかしながら、日本ではGeoHPシステムの普及は、欧米に比べ、遅れているのが実情である。その理由としては、認知度の低さ、井戸の掘削等における初期費用の問題、および、日本では欧米とは違い、その土地の風土によって熱物性値が異なるため、経済的なシステムづくりが依然として確立されていない、などが挙げられる。そのため、GeoHPシステムの普及には初期費用を抑えるためのシステムの効率化が必要不可欠である。
【0004】
システムの効率化を図るための地盤調査試験として、熱応答試験が知られている。熱応答試験は主に地層の熱伝導率と熱交換井の熱抵抗を評価するための試験であり、地中熱交換井の熱交換挙動を予想する上で不可欠な情報である。熱応答試験を実施して、適正な井戸の本数・長さを決定することは、GeoHPシステムの初期費用の削減に極めて重要である。
【0005】
地盤調査のための熱応答試験は、主にクローズド型方式における垂直型地中熱交換井を対象としている。このための熱応答試験装置は、地上部の水タンク内で電気ヒーターによって熱媒体を加熱し、循環ポンプにより配管内を循環させる。加熱された熱媒体は地中部を通過する際に放熱するため、U字管の入口・出口の温度差を測定し、その地層の熱交換特性を評価する。このような熱応答試験装置は、特許文献1に開示されており、熱応答試験装置を用いた熱応答試験方法は非特許文献1に開示されている。
【0006】
熱応答試験の代表的な解析方法としては、地盤の温度上昇量を規定するKelvinの線源関数や円筒型熱源関数に基づく解析的手法、および有限差分法や有限要素法を用いた数値モデルによる解析方法がある。熱応答試験の解析的手法において最も一般的に用いられている手法としては、データ観測の簡便さや解析の容易さより、Kelvinの線源関数に基づいた作図法が挙げられる。作図法では循環時における熱媒体の温度データの経時変化を用いる方法と循環停止後の地中温度の回復データを用いる方法があるが、循環時のデータを用いる解析が一般的に行われている。
【0007】
本発明者らは、循環時データを用いた作図法で熱応答試験の解析を行う場合に、解析を行うデータ区間の選定に明確な基準が無いため、解析結果が人によって異なること、地中熱交換器や充填剤の影響を排除するため、試験期間の初期段階のデータは捨てるので試験
期間を長くとる必要があること、などの問題点に着目し、地盤の温度上昇量を規定するKelvinの線源関数の近似解における未知のパラメータを同定して、地盤の熱伝導率を算出する解析方法および解析プログラムを提案している(特許文献2参照)。特許文献2で提案した解析方法では、熱応答試験の初期段階において得られる実測値に対して、パウエルの共役傾斜法に基づく非線形最適化手法を用いて、未知のパラメータを所定の値に設定した場合に得られる温度上昇量の計算値と、熱応答試験によって測定された温度上昇量の実測値との誤差が最小となるように、パラメータの逆解析を行い、未知のパラメータを同定して、地盤の熱伝導率を算出する。
【0008】
また、本発明者らは、透水量係数や貯留係数といった地下水流動系パラメータを求めるための原位置試験として最も一般的に行われている揚水試験は、熱応答試験と同様に、熱伝導理論を基として解析解を導いていることに着目し、特許文献2において、熱応答試験の解析方法を揚水試験にも適用することを提案している。すなわち、地盤の地下水位低下量を規定するタイスの井戸関数の近似解(ヤコブの近似解)における未知のパラメータを同定して、地盤の透水量係数を算出するにあたって、揚水試験の初期段階において得られる実測値に対して、パウエルの共役傾斜法に基づく非線形最適化手法を用いて、未知のパラメータを所定の値に設定した場合に得られる地下水位低下量の計算値と、揚水試験によって測定された地下水位低下量の実測値との誤差が最小となるように、パラメータの逆解析を行い、未知のパラメータを同定して、地盤の透水量係数を算出している。
【先行技術文献】
【0009】

【特許文献1】特開2004-301750号公報
【特許文献2】特許第5334221号公報
【0010】

【非特許文献1】「ボアホール型地中熱交換機に対する加熱法による熱応答試験の標準試験方法ver.2」、長野克則著、地下熱利用とヒートポンプシステム研究会編、財団法人ヒート・ポンプ・蓄熱センター
【発明の概要】
【発明が解決しようとする課題】
【0011】
従来から行われている作図法による熱応答試験の解析方法や、本発明者らが特許文献2で提案したパウエルの共役傾斜法に基づく非線形最適化手法を用いたパラメータ同定による熱応答試験の解析方法は、Kelvinの線源関数に基づいた熱応答試験の解析理論を前提としている。そして、この解析理論は、熱移動が熱伝導のみに支配されることを前提としている。このため、地下水流動による移流・分散の影響を考慮に入れた解析を行うことができない。従って、地下水流動がある場合には、熱応答試験の解析により求まる熱伝導率や体積熱容量が誤差を含んだものとなってしまう。また、地下水流速や縦・横分散長等のパラメータを求めることはできない。このため、地下水流動があるサイトでの熱交換井の熱交換挙動を精度良く予測することができない。
【0012】
同様に、揚水試験の解析においても、従来の解析方法では、隣接帯水層からの漏水がある場合には、解析から求まる透水量係数や貯留係数が誤差を含んだものとなってしまう。このため、隣接帯水層からの漏水があるサイトでの地下水位の挙動を精度良く予測することができない。
【0013】
本発明の課題は、このような点に鑑みて、地下水流動がある場合でも、熱伝導率や体積熱容量、地下水流速、縦・横分散長などのパラメータを精度良く同定することができる、熱応答試験の解析方法および解析プログラムを提案することにある。
【0014】
また、本発明の課題は、隣接帯水層からの漏水がある場合でも、透水量係数や貯留係数を精度良く求めることができる揚水試験の解析方法および解析プログラムを提案することにある。
【課題を解決するための手段】
【0015】
上記の課題を解決するために、本発明の熱応答試験の解析方法および解析プログラムでは、地盤の温度上昇式として、地下水が流動している流動場での移流・分散現象を考慮した熱移動方程式の解析解を用いる。そして、この解析解の近似式(高次漸近解)を用いて求められる地盤の温度上昇値の計算値と、熱応答試験で得られる地盤の温度上昇値の実測値との誤差が最小となるように、熱移動方程式の解析解に含まれる未知のパラメータを同定する逆解析を行う。逆解析には、パウエルの共役傾斜法に基づく非線形最適化法を用いる。
【0016】
また、本発明の揚水試験の解析方法および解析プログラムでは、揚水試験の解析において、地下水が隣接帯水層の影響を受けている流動場(隣接帯水層から漏水がある被圧帯水層)での地下水位低下量を、移流・分散現象を考慮した地下水流動方程式の解析解(ハンタッシュ-ヤコブの式)で規定し、これと上記の熱移動方程式の解析解との同一性に着目して、この解析解の近似式(高次漸近解)を用いて求められる地下水位低下量の計算値と、揚水試験で得られる地下水位低下量の実測値との誤差が最小となるように、地下水流動方程式の解析解に含まれる未知のパラメータを同定する逆解析を行う。逆解析には、パウエルの共役傾斜法に基づく非線形最適化法を用いる。
【0017】
(移流・分散現象を考慮した熱移動の基礎理論に基づく熱応答試験の解析方法)
以下に、地盤の温度上昇式として、地下水などの流体の流れに伴う熱移動を考慮した熱移動方程式の解析解を求めるにあたって、本発明者らが用いた解析理論を述べる。
【0018】
飽和・不飽和多孔体中における熱移動には熱伝導、熱分散、熱交換、熱対流といった現象を考慮する必要がある。図1は、空気や水などの流体で飽和した多孔体中を流体が流動している時の熱の流れを模式的に表したものである。図中のルートAは熱伝導による固相中の熱の流れ、ルートBは熱伝導による流体相中の熱の流れ、ルートCは固相と流体相の接触面を通して行われる熱交換、ルートDは流体の流れに伴う機構的分散による熱移動を表している。
【0019】
まず、空気や水などの流体によって飽和した多孔体中における熱移動において固相・液相・気相の温度が平衡状態にあるものとすると、飽和多孔体中における熱移動の基礎方程式は、移流項と分散項を含む式(1)の移流分散方程式で表される。
JP2016017773A_000003t.gif但し、 T:温度(℃)
v:間隙流速(m/s)
n:間隙率
(ρc):液相の体積熱容量(J/m/K)
(ρc):固相・液相・気相混合系の等価体積熱容量(J/m/K)
λ:固相・液相・気相の熱伝導率に機構的熱分散を加味した混合系の熱分散率
(W/m/K)
【0020】
ここで、X軸方向の地下水流動がある2次元流動場を考え、式(1)の両辺を(ρc)で割ると、次式(2)、(3)で表される飽和多孔体の熱移動方程式が得られる。
JP2016017773A_000004t.gif
【0021】
X軸方向に間隙流速vで地下水が流動している2次元流動場において、原点に与える熱負荷として、微少時間dτで熱量qが連続的に投入されるものとし、式(4)~(7)の初期条件・境界条件のもとで式(2)を解くと、2次元流動場の温度差分布として、次式(8)の解析解が得られる。
JP2016017773A_000005t.gif
【0022】
式(8)において、次式(9)のように置くと、式(8)は、次式(10)で表されるW(u、B)を含んでいる。W(u、B)はハンタッシュ-ヤコブ(Hantush-Jacob)の井戸関数と呼ばれ、隣接帯水層から漏水がある被圧帯水層中の透水量係数Tと貯留係数Sを求めるための揚水試験の解析解として知られている。
JP2016017773A_000006t.gif
【0023】
ここで、式(8)において、v=0とし、原点からの距離をrとすると、式(8)は次式(11)、式(12)のようになる。式(11)は、Kelvinの線源関数そのものである。Kelvinの線源関数は、従来から熱応答試験の解析に用いられてきた、移流・分散現象を考慮しない熱移動方程式の解析解である。
JP2016017773A_000007t.gif
【0024】
式(11)に含まれるW(u)は、タイス(Theis)の井戸関数である。式(10)で表されるW(u、B)において、B=0とすると、W(u、0)=W(u)となり、タイスの井戸関数と一致する。タイスの井戸関数は、従来から揚水試験の解析に用いられてきた、隣接帯水層の干渉を考慮しない地下水流動方程式の解析解であり、Kelvinの線源関数と式の形が同じである。このように、本発明で用いる、移流・分散現象を考慮した熱移動方程式の解析解は、隣接帯水層の干渉がない条件下では、従来用いられていた、地下水流れを考慮しない熱応答試験や揚水試験の解析解と一致する。
【0025】
また、地下水流動場において、原点に連続的に注入したトレーサの濃度差分布は、熱移動と同じ解析理論によって解析することが可能である。すなわち、上記の熱移動方程式の解析解のパラメータを、物質輸送式に用いるパラメータに置き換えることで、移流・分散現象を考慮した物質輸送式の解析解が得られる。
【0026】
(高次漸近解の適用)
式(10)で表されるハンタッシュ-ヤコブの井戸関数W(u、B)は、u≦1.0、B≦2.0の範囲において、次式(13)で近似される。
JP2016017773A_000008t.gif
【0027】
式(13)において、I(B)は第1種変形ベッセル(Bssel)関数、K(B)は第2種変形ベッセル(Bssel)関数である。また、上述したように、W(u)はタイスの井戸関数である。ここで、式(13)のI(B)、K(B)、W(u)は、それぞれ、次式(14)~(16)のべき級数に展開できる。
JP2016017773A_000009t.gifJP2016017773A_000010t.gifJP2016017773A_000011t.gif
【0028】
本発明者らは、式(14)~(16)のべき級数展開式の正確解と、べき級数展開式の第3項、第30項、第50項・・・などの項数までを用いて近似した近似解との一致度を検討した。図2.1にI(B)の正確解と第3項、第4項、第30項、第50項までを用いて近似した近似解を示す。また、図2.2にK(B)の正確解と第3項、第4項、第30項、第50項までを用いて近似した近似解を示す。I(B)では、第3項までの近似解はBが小さい場合は正確解と良く一致するが、Bが大きくなるにつれて正確解から乖離する。また、K(B)では、第3項、第4項までの近似解はBが小さい場合は正確解と良く一致するが、Bが大きくなるにつれて正確解から乖離する。I(B)、K(B)とも、第30項、第50項まで近似した場合はどの区間においても正確解と良く一致する。I(B)、K(B)は、ともに第50項まで近似した近似解を用いることが好ましい。
【0029】
図2.3にタイスの正確解と第2項、第6項、第30項、第100項までを用いて近似した近似解を示す。第2項までの近似解は、1/Uが大きい後半部では正確解と良く一致するが、1/Uが小さくなるにつれて正確解から乖離する。近似する項数を増加させると近似解は正確解に漸近し、第30項、第100項までで近似した場合はどの区間においても正確解と良く一致する。
【0030】
図2.4に、式(13)で示すW(u、B)の正確解(ハンタッシュ-ヤコブの標準曲線)と、I(B)、K(B)、W(u)のべき級数展開式を高次項(少なくとも、第30項)まで用いて近似した近似解(漸近解)との比較結果を示す。このように、漸近解の項数を増加させることで、正確解と非常に良く一致する計算値が得られる。本発明では、以下に説明する逆解析を行うにあたって、I(B)、K(B)、W(u)のべき級数展開式を、少なくとも第30項まで(I(B)、K(B)については好ましくは第50項まで)用いた解析プログラムを用いる。
【0031】
(パウエルの共益傾斜法)
本発明においては、熱応答試験あるいは揚水試験の解析を、最適化手法による逆解析を利用して行う。先に述べた地下水モデルによる熱応答試験の解析解や、物質輸送式の解析解、あるいは、後述する揚水試験の解析解は、透水量係数Tや熱伝導率λあるいは分散係数などを既知パラメータとして、地下水位、地下水温、物質の濃度などを計算する構造になっている。しかし、こうしたパラメータを直接精度よく測定することは難しく、むしろ解として得られるべき地下水位、地下水温、物質の濃度などの方が容易に、精度よく観測・測定が可能である。本発明では、地下水位や地下水温などの実測値に計算値が一致するように非線形最適化手法を用いて、各種パラメータを同定する逆解析を採用する。逆解析では、まず目的関数を設定し、それらを最小化することによって各パラメータを同定する。
【0032】
本発明では、最適化手法で用いる目的関数Fとして、熱応答試験や揚水試験で得られた実測値と、解析解から得られた計算値の2乗誤差の総和を1/2乗したユークリッドノル
ムに、係数Aおよび係数Aを乗じて重み付けした次式(17)を用いる。
JP2016017773A_000012t.gif
【0033】
また、最適化手法において、収束判定値として、ある探索点と、その1つ前の探索点に対して得られた目的関数Fの値の差分であるDIF=F-Fj-1を用いる。更に、データ1つあたりの誤差を示す評価尺度には、次式(18)で表すように、計算値と実測値の絶対値誤差の総和をデータ数Nで除した値f/Nを用いる。
JP2016017773A_000013t.gif
【0034】
本発明では、目的関数の未知のパラメータを最適値に収束させるための、非線形最適化手法として、パウエル(Powell)の共役傾斜法を採用する。パウエル(Powell)の共役傾斜法は、他の多くの最適化手法で必要とされる目的関数の各パラメータに関する微係数を必要とせず、目的関数のみを用いて、その最適パラメータを求めることができ、高い汎用性を持つ。
【0035】
図3.1に共役方向の概略図を示し、図3.2にパウエル(Powell)の共役傾斜法の基本アルゴリズムを示す。パウエル(Powell)の共役傾斜法では、2つの離れた点からvの方向に正値2次形式の最小点を見つけ出すとき、各々で求めた最小点x1,x2を結ぶベクトルは、vと共役であるという性質を利用している。この性質を利用して、n回直接探索を繰り返すことにより、最小点を見つけ出すことができる。
【0036】
(熱応答試験の解析における同定対象のパラメータ)
熱応答試験において温度T、加熱量q等は観測データとして既知である。式(8)の解析解において、等価体積熱容量(ρc)、縦熱分散率κ、横熱分散率κ、真流速が未知の状態である。本発明では、同定対象のパラメータをこれらの4つのパラメータとして、式(17)の目的関数Fを用いて、共役傾斜法を組み込んだ同定プログラムにより逆解析を行う。
【0037】
本発明の逆解析では、式(17)における係数Aの値は、次式(19)に示すように、同定対象のパラメータであるκ、κの設定値から算出した値aの正負に応じて決定される。aが正の値になるのは、κがκよりも大きい場合である。しかしながら、X軸方向(縦方向)に地下水の流れがあることを想定している場合には、κがκよりも小さいはずである。つまり、係数Aは、同定対象のパラメータのうち、κとκがありえない大小関係となるように設定された場合に、1.0よりも大きな値となって目的関数Fを増大させるように設定される。

=(κ-κ)/κ≧0.0の場合、A=1.0+a
=(κ-κ)/κ<0.0の場合、A=1.0 (19)

【0038】
同様に、係数Aの値は、次式(20)に示すように、同定対称のパラメータであるκ、(ρc)と、λの設定値から算出した値aの正負に応じて決定される。aが正の値になるのは、κ=αv+λ/(ρc)と置くと、流速vが負の値になる場合であるが、流速vは正の値のはずである。つまり、係数Aは、同定対称のパラメータのうち、流速vがありえない範囲の値となるように設定された場合に、1.0よりも大きな値となり、目的関数Fを増大させるように設定される。

=10(1/κ(ρc)-1/λ)≧0.0のとき、A=1.0+a
=10(1/κ(ρc)-1/λ)<0.0のとき、A=1.0 (20)

【0039】
以上のように、係数Aおよび係数Aは、いずれも、パラメータの探索が間違った方向に向かった場合に、目的関数Fが増大する(ペナルティが加わる)ように設定されている。このような係数を設けることで、目的関数Fを最小値に収束させるまでの、最適化手法による探索回数を少なくすることが可能になる。
【0040】
(熱応答試験の解析方法の揚水試験への適用)
熱応答試験と揚水試験は、両者とも熱伝導理論を基として、解析解を導いている。そのため、本発明の熱応答試験の解析方法および解析プログラムは、揚水試験の解析にも適応することが可能である。揚水試験は現場においては、透水量係数Tや貯留係数Sといった地下水流動系パラメータを求めるための原位置試験として最も一般的に行われている。
【0041】
いま、隣接帯水層から漏水があり、無限の広がりを持つ被圧帯水層の非定常流の地下水流動方程式は、地下水位の降下量をs=H-hとし、隣接帯水層から漏水がある被圧帯水層へ完全貫入している井戸から一定量Qで地下水が揚水されているものとすると、次式(21)(22)で与えられる。また、初期条件・境界条件は次式(23)~(25)で与えられる。
JP2016017773A_000014t.gif但し、
h:ピエゾ水頭(m)
:隣接帯水層のピエゾ水頭(m)
r:半径(m)
:揚水している帯水層の透水量係数(m/s)
:層厚
K:加圧層の透水係数
:貯留係数
:揚水量(L/s)
【0042】
上記の初期条件、上記の初期条件・境界条件のもとで式(21)を解くと、次式(26)(27)の解析解が得られる。
JP2016017773A_000015t.gif
【0043】
式(26)は、ハンタッシュ-ヤコブの式と呼ばれている。r/λ=Bと置くと、式(26)は、次式(28)で表されるハンタッシュ-ヤコブの井戸関数W(u、B)を含む。
JP2016017773A_000016t.gif
【0044】
熱応答試験の解析で用いる、移流分散方程式から導出した式(8)の熱移動解析解においても、ハンタッシュ-ヤコブの井戸関数W(u、B)が用いられている。すなわち、式(26)の揚水試験に用いる解析解と、熱応答試験に用いる移流分散を考慮した解析解(式(8)は、式の形が良く似ている。そこで、熱応答試験の解析方法と同様に、揚水試験の解析でも、パウエル(Powell)の共役傾斜法を用いた逆解析を行う。
【0045】
(揚水試験の解析における同定対象のパラメータ)
次に、揚水試験においては、式(26)の解析解において、透水量係数T、貯留係数S、Bが未知の状態である。本発明では、求めるパラメータを次式(29)~(31)に示すα1、α2、Bとして、共役傾斜法を組み込んだ同定プログラムにより逆解析を行う。そして、同定されたパラメータα1、α2、Bから、式(29)~(31)より透水量係数T、貯留係数S、Bを算出する。
α1=r (29)
α2=T (30)
B=r/p (31)
但し、
JP2016017773A_000017t.gif
【0046】
以上のように、熱応答試験の解析解と揚水試験の解析解は、元を辿れば熱伝導理論によるものであり、パラメータは異なるが、式の形自体は同じである。これにより、本発明による熱応答試験の解析方法を揚水試験にも同様に適応できることが分かる。
【発明の効果】
【0047】
本発明によれば、地下水流動による移流・分散の影響を考慮に入れた解析解を用いて解析を行うことができる。また、パウエルの共役傾斜法を適用した逆解析により、熱応答試
験の測定データから未知のパラメータを同定でき、従来の解析では求めることのできなかったパラメータについても同定できる。具体的には、地下水流動があるサイトでの熱応答試験の解析により、等価体積熱容量(ρc)、縦熱分散率κ、横熱分散率κ、真流速などのパラメータを精度良く求めることができる。よって、地下水流動があるサイトでの熱交換井の熱交換挙動を従来よりも精度良く予測することが可能である。
【0048】
また、本発明によれば、揚水試験の解析においても同様な作用効果を得ることができる。すなわち、透水量係数T、貯留係数S、Bなどのパラメータを精度良く求めることができる。よって、地下水流動があるサイトでの地下水位の挙動を従来よりも精度良く予測することが可能である。
【図面の簡単な説明】
【0049】
【図1】飽和多孔体中を流体が流動している時の熱の流れを模式的に示す説明図である。
【図2.1】第1種変形ベッセル関数I(B)の正確解とその近似解(漸近解)である。
【図2.2】第2種変形ベッセル関数K(B)の正確解とその近似解(漸近解)である。
【図2.3】タイスの正確解W(u)とその近似解(漸近解)である。
【図2.4】W(u、B)の正確解と、その漸近解との比較結果である。
【図3.1】パウエルの共役傾斜法における共役方向を示す概念図である。
【図3.2】パウエルの共役傾斜法の基本アルゴリズムを示す説明図である。
【図4.1】温度回復データを用いた解析方法の検証結果(テストケース1)である。
【図4.2】温度回復データを用いた解析方法の検証結果(テストケース2)である。
【図4.3】熱応答試験結果(往き温度、還り温度、平均温度、加熱量)を示すグラフである。
【図4.4】熱伝導率λが異なる複数のテストケースの回復時の温度変化量の計算結果である。
【図4.5】r(ρc)が異なる複数のテストケースの回復時の温度変化量の計算結果である。
【図5.1】熱交換井直下でのテストケース1~3の温度変化量の計算結果である。
【図5.2】熱交換井直下での室内熱応答試験による計測値、および原位置熱応答試験による計測値を示すグラフである。
【図6.1】同定対象のパラメータの影響を検証する感度解析結果(座標1)である。
【図6.2】同定対象のパラメータの影響を検証する感度解析結果(座標2)である。
【図6.3】同定対象のパラメータの影響を検証する感度解析結果(座標3)である。
【図6.4】パラメータの探索過程(目的関数F)を示すグラフである。
【図6.5】パラメータの探索過程(収束判定値DIF)を示すグラフである。
【図6.6】パラメータの探索過程(等価体積熱容量(ρc))を示すグラフである。
【図6.7】パラメータの探索過程(縦熱分散率κ)を示すグラフである。
【図6.8】パラメータの探索過程(横熱分散率κ)を示すグラフである。
【図6.9】パラメータの探索過程(真流速)を示すグラフである。
【図6.10】逆解析より求めたパラメータを用いて温度変化量を再現計算した値と正確値との比較結果である。
【図6.11】室内熱応答試験装置を模式的に示す全体構成図である。
【図6.12】土槽の説明図である。
【図6.13】実験地盤の中心部に埋設された温度センサーの配置図(XY面上)である。
【図6.14】実験地盤の中心部に埋設された温度センサーの配置図(浸透流に垂直な断面上:YZ面上)である。
【図6.15】熱応答試験の計測結果(実験ケース1:U字管の出入口温度、平均温度、熱交換量)である。
【図6.16】熱応答試験の計測結果(実験ケース1:温度変化量の経時変化)である。
【図6.17】実験ケース1の同定結果を用いて演算した温度変化量の計算値と実測値との比較結果である。
【図6.18】熱応答試験の計測結果(実験ケース2:U字管の出入口温度、平均温度、熱交換量)である。
【図6.19】熱応答試験の計測結果(実験ケース2:温度変化量の経時変化)である。
【図6.20】実験ケース2の同定結果を用いて演算した温度変化量の計算値と実測値との比較結果である。
【図6.21】等価体積熱容量(ρc)の同定結果と簡易試験より得られた実測値である。
【図6.22】真流速の同定結果と簡易試験より得られた実測値である。
【図6.23】縦熱分散率κの同定結果と簡易試験より得られた実測値である。
【図6.24】横熱分散率κの同定結果と簡易試験より得られた実測値である。
【図6.25】縦分散長αの同定結果である。
【図6.26】縦分散長αの同定結果である。
【図6.27】逆解析より求めたパラメータを用いて水位降下量を再現計算した値と正確値との比較結果(テストケース1)である。
【図6.28】逆解析より求めたパラメータを用いて水位降下量を再現計算した値と正確値との比較結果(テストケース2)である。
【図6.29】揚水試験を実施したサイトの説明図である。
【図6.30】現場内の深井戸の掘削データをもとにした断面図である。
【図6.31】揚水井から23.5mの距離(第一帯水層)にある観測井の地下水位の変動の計測値である。
【図6.32】同定結果を用いて揚水試験を再現計算したときの観測井における水位降下量の計算値と実測値との比較結果である。
【図6.33】揚水井から24.5mの距離(第二帯水層)にある観測井の地下水位の変動の計測値である。
【図6.34】同定結果を用いて揚水試験を再現計算したときの観測井における水位降下量の計算値と実測値との比較結果である。
【発明を実施するための形態】
【0050】
[温度回復データを用いる解析方法]
まず、移流・分散現象を考慮した、本発明の熱応答試験の解析方法の妥当性検証に先立って、解析に用いるデータ区間の検討について述べる。本発明者らは、以前に開発した解析プログラム(移流・分散現象を考慮していない熱移動方程式の解析解(Kelvinの線源関数の解析解)に基づき、パウエルの共役傾斜法を用いた逆解析を行ってパラメータを同定するプログラム:特許第5334221号公報)では、熱応答試験において熱交換器に熱媒体を循環している間の計測値(循環時データ)を用いていた。ここでは、解析に用いるデータ区間を広げるため、循環停止後の温度回復データを用いることを検討する。

【0051】
(順解析による検証)
以下の解析方法は、移流・分散現象を考慮していない、Kelvinの線源関数の解析解を用いるものである。この場合、加熱停止後の温度変化量は、重ね合わせの原理により、次式(32)~(34)で表される。そこで、E(x)とE(x´)をべき級数展開式の第30項まで評価した高次漸近解を用いて、パウエルの共役傾斜法を適用したパラメータ同定による逆解析を行う解析プログラムを開発した。
JP2016017773A_000018t.gif

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

【0053】
図4.1と図4.2は、温度回復データを用いた解析方法の検証結果であり、テストケース1、2について、逆解析より求めたパラメータを用いて温度変化量を再現計算した値と正確値を示す。また、表1.1と表1.2はパラメータの同定結果を示す一覧表である。このように、同定されたパラメータを用いた計算結果と正確値は非常に良く一致することが確認できた。また、f/N値(データ1つ当たりの誤差)もほぼ0に等しく、精度の高い結果を得ることが出来た。
【表1.1】
JP2016017773A_000019t.gif
【表1.2】
JP2016017773A_000020t.gif

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

【0055】
逆解析による同定結果を表1.3に示す。「全データ」は加熱時間と同程度の温度回復データを用いた場合、「前半部分」は温度回復初期のデータを用いた場合、「循環時法」は循環時データを用いた場合である。温度回復データの全データと温度回復初期データのどちらを用いた場合も、熱伝導率λは循環時法で得られた同定結果2.51(W/m/K)とほとんど等しいが、r(ρc)は、温度回復初期データを用いた方が近い値となった。なお、f/N値はそれぞれ0.184、0.133で、温度回復初期データを用いた結果の方が実測値との適合性が高いことが確認された。
【表1.3】
JP2016017773A_000021t.gif

【0056】
また、本発明者らは、温度回復初期データを用いた場合に、解析解のパラメータの値(熱伝導率λ、r(ρc))が温度変化量に与える影響を検討するため、パラメータの値が異なる複数のテストケースで計算を行った。表1.4は熱伝導率λを変化させたケース、表1.5はr(ρc)を変化させたケースである。
【表1.4】
JP2016017773A_000022t.gif
【表1.5】
JP2016017773A_000023t.gif

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

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

【0059】
[移流・分散現象を考慮した熱応答試験の解析方法]
続いて、本発明の実施の形態に係る移流・分散現象を考慮した熱応答試験の解析方法の妥当性検証について述べる。本発明者らは、以下の1~5の順で妥当性検証を行った。
1.順解析による、地下水流動の影響の検証
2.物質移動との比較
3.パラメータの感度解析
4.テストケースによるパラメータ同定精度の検証
5.室内熱応答試験データによるパラメータ同定精度の検証

【0060】
1.順解析による、地下水流動の影響の検証
本発明者らは、本発明の移流・分散現象を考慮した熱移動方程式の解析解に基づいて温度上昇量を演算するプログラムをFORTRANコードで開発した。そして、表2に示す3つのテストケース(テストケース1~3)のパラメータを用いて、温度上昇量を演算する順解析を行った。テストケースのパラメータは、いずれも、飽和砂を用いた室内実験値および文献値を用いた。テストケース1~3は、流速vの値だけがケースによって異なり、他のパラメータは全て同じである。また、縦分散長αと横分散長αの比が10.0となるようにパラメータを設定した。
【表2】
JP2016017773A_000024t.gif

【0061】
表3は、テストケース1~3における等価熱拡散率κと機構的分散α×vを示している。テストケース1から3の順に、流速vの値が大きくなるにつれて、機構的分散α×vの方が等価熱拡散率κよりも大きくなっている。
【表3】
JP2016017773A_000025t.gif

【0062】
演算結果により、以下の事項が確認された。
(1)テストケース1(流速vが最小)
等価熱拡散率κの方が機構的分散α×vよりも大きいことを反映し、同心円状に熱が拡がる結果となった。熱の拡がり範囲は、計算時間1000hrでも、最大で5.0mとなった。
(2)テストケース2(流速vがテストケース1よりも1桁大きい)
機構的分散α×vの方が等価熱拡散率κよりも大きくなっており、熱分散が支配的になったことを反映して、楕円状に熱が拡がる結果となった。また、移流効果も伴って、計算時間1000hrでは10.0mまで熱が拡がった。
(3)テストケース3(流速vがテストケース2よりもさらに1桁大きい)
テストケース2よりも更に移流・分散効果が卓越し、より広域に熱が拡がっており、計算時間1000hrでは30.0mまで熱が拡がった。

【0063】
図5.1は、熱交換井直下(座標:x=1.0m、y=0.0m:熱注入位置から1m下流の位置)におけるテストケース1~3の温度変化量(温度上昇量)の計算結果を示すグラフである。この図に示すように、流速が小さいとき、熱伝導が支配的であるため、温度の立ち上がりは緩やかであるが、時間の経過に伴って温度は上昇し続ける。流速が大きくなるにつれて、加熱開始直後に急激に温度が変化し、時間の経過とともに温度変化が緩慢になる。これは、図2.4に示した、ハンタッシュ-ヤコブの標準曲線からもわかるように、流速が小さくなるとタイスの曲線に漸近し、流速が大きくなると早く定常状態に達するようになるためである。

【0064】
ここで、図5.2に、図5.1と同じ座標(x=1.0m、y=0.0m)における室内熱応答試験による計測値、および原位置熱応答試験による計測値を示す。本発明者らは、後述するように、室内熱応答実験装置を用いて熱応答試験を実施しており、その結果、図5.1で説明した順解析の結果と同様の傾向がみられることが確認された。また、現場サイトの2地点で実施した熱応答試験の結果から、地下水流動が少ないとみられるサイトBと、地下水流動が多いとみられるサイトAでの計測値は、順解析と同様の傾向となっていることが確認された。

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

【0066】
2:物質移動との比較
上述したように、地下水流動場において、原点に連続的に注入したトレーサの濃度差分布は、熱移動と同じ解析理論によって解析できることから、熱移動方程式の解析解に基づいて温度上昇量を演算するプログラムは、物質輸送式の解析解についても適用できる。そ
こで、熱移動と物質移動との違いを比較検証するために、熱移動のテストケース1、2と同じ計算条件(表2参照)を用いて、濃度差分布を演算する順解析を行った。表4に、物質移動の演算に用いる2つのテストケースのパラメータを示す。また、表5に、各ケースの有効分子拡散係数Defと機構的分散α×vを示す。表4では、熱移動と物質移動の演算結果を比較しやすくするために、q/(ρc)=M/nとなるように、M(質量)の値を設定した。一般的に,溶質の有効分子拡散係数 は10(cm/s))のオーダーで非常に小さいため、物質移動では移流・分散現象が卓越することが多い。
【表4】
JP2016017773A_000026t.gif
【表5】
JP2016017773A_000027t.gif

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

【0068】
3.パラメータの感度解析
本発明者らは、熱移動方程式の解析解を用いた逆解析において、同定対象のパラメータに設定した等価体積熱容量(ρc)、縦熱分散率κ、横熱分散率κ、真流速が解析解に与える影響を検証する感度解析を行った。感度解析には、表6.1に示すパラメータの値を用いた。表6.1に示す3点の座標(x=2.0m、y=1.0m)、(x=2.0m、y=0.0m)、(x=0.0m、y=1.0m)において、同定対象の各パラメータを、表6.1に示す値(正確値)、正確値の10倍、正確値の1/10、に設定した場合の温度上昇量を演算した。
【表6.1】
JP2016017773A_000029t.gif

【0069】
図6.1、図6.2、図6.3に各座標における感度解析結果を示す。図6.1に示すように、座標1(x=2.0m、y=1.0m)では、等価体積熱容量(ρc)が最も大きく解析解へ影響を与えており、パラメータを正確値(表6.1の値)とした場合の演算結果との乖離が最も大きい。次いで縦熱分散率κ、横熱分散率κ、真流速の順に、影響(正確値を用いた演算結果との乖離)が大きい。一方、図6.2に示すように、座標2(x=2.0m、y=0.0m)では、等価体積熱容量(ρc)が最も大きく解析解へ影響を与えており、次いで横熱分散率κ、縦熱分散率κ、真流速の順に影響が大きい。また、座標3(x=0.0m、y=1.0m)は流動方向と直交する点であるが、等価体積熱容量(ρc)、縦熱分散率κ、横熱分散率κ、真流速の順に影響が大きい。

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

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

【0072】
図6.4~6.10はパラメータの探索過程を示すグラフであり、図6.4は目的関数F、図6.5は収束判定値DIFである。また、図6.6は等価体積熱容量(ρc)、図6.7は縦熱分散率κ、図6.8は横熱分散率κ、図6.9は真流速である。また、図6.10は逆解析より求めたパラメータを用いて温度変化量を再現計算した値と正確値(図6.1参照)との比較である。

【0073】
正確値に対し、適当な初期推定値を与えて逆解析を行った。目的関数Fが最小値を目指すよう探索を繰り返し、これ以上値が変化せず、DIF値が限りなく0に近い値を示したところで探索を終了した。図6.6~6.10から、探索を繰り返すごとに、パラメータの値が正確値に徐々に近づき、同定されたパラメータを用いた計算結果と正確値(図6.10)は、非常に良く一致することが確認できた。表6.3は各パラメータの同定結果を
示す一覧表である。表6.3に示すように、f/N値(データ1つ当たりの誤差)もほぼ0に等しく、精度の高い結果を得ることが出来た。以上の検証より、本発明の逆解析方法の妥当性を確認することが出来た。
【表6.3】
JP2016017773A_000031t.gif

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

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

【0076】
図6.11に示すように、水位設定部12は、土槽11の上流側に配置された上流側水位調節部12Aと、土槽11の下流側に配置された下流側水位調節部12Bを備える。上流側水位調節部12Aは給水容器14Aおよび上流側水位調整容器14Bを備え、下流側水位調節部12Bは下流側水位調整容器15を備える。上流・下流の水位調節部12A、12Bの高さを調整することで、土槽11内に所定の水位が維持される。

【0077】
図6.12は土槽11の説明図である。土槽11は、厚さ20mmの透明アクリル板製の箱型の容器であり、その内法寸法は、高さ480mm、幅960mm、奥行き460mmである。図6.12に示すように、土槽11の内部には幅700mmの試料充填部16が設けられ、試料充填部16の幅方向Xの両側に貯水部17A、17Bが設けられている。試料充填部16と貯水部17A、17Bは、アクリル板にパンチ穴を設け、さらに金網で被覆した透水性スリット18により分離されている。土槽11はその幅方向Xに傾斜するように配置され、土槽11内では幅方向Xに浸透流が流れる。試料充填部16の上流側に設けられた貯水部17Aの幅は175mmであり、下流側に設けられた貯水部17Bの幅は75mmである。

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

【0079】
熱応答実験装置20は、図6.11に示すように、熱媒体(水)を加熱し、循環を行う循環系と、熱媒体温度・消費電力・循環流量などを計測する計測系を備える。循環系は、実験地盤10に設置された熱交換器としてのU字管21を含む循環流路22、循環流路22で熱媒体を循環させる循環ポンプ23と、循環流路22を通る熱媒体を加熱する電気ヒーター24、循環流路22に熱媒体を供給する給水容器25等を備える。計測系は、U字管21の入口・出口の温度を計測する温度センサー26a、26b、実験地盤10に埋設される複数の温度センサー27、流量計28、データロガー29、パーソナルコンピューター30、室温計31等を備える。

【0080】
図6.13および図6.14は、実験地盤10の中心部に埋設された温度センサー27の配置図である。図6.12に示すように、本明細書では、X方向を土槽11の幅方向とする。また、Y方向(図6.12の紙面に垂直な方向:図示省略)を土槽11の奥行き方向とし、Z方向を土槽11の高さ方向とする。X方向、Y方向、Z方向は互いに直交する方向である。図6.13は温度センサー27のXY面上の配置図であり、図6.14はU字管21を通り浸透流に垂直な断面(YZ面)上の温度センサー27の配置図である。

【0081】
この構成の室内熱応答試験装置1では、電気ヒーター24によって熱媒体を加熱し、循環ポンプ23によって循環流路22に循環させる。加熱された熱媒体は、試験対象の実験地盤10に形成されているボアホール孔を通過する際に放熱するため、U字管21の入口・出口の温度差を温度センサー26a、26bでそれぞれ測定し、実験地盤10の熱交換特性を評価する。

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

【0083】
本発明者らは、室内熱応答試験装置1を24℃設定の恒温室内に設置して熱応答試験を行い、本発明の逆解析方法を適用してパラメータの同定を行った。表6.4に実験条件(実験ケース)を示す。本試験では、実験地盤10の傾斜角度が異なる2つの実験ケース1,2を設定した。実験地盤10に傾斜をつけることで地下水流れ場を発生させ、各実験ケースの浸透流量、浸透流が通過する断面積(2208cm)、間隙率(0.463)からダルシー流速、真流速を算出した。2つの実験ケースの単位長さ当たりの加熱量qはそれぞれ、122.925(W/m)、122.675(W/m)とした。逆解析時には,加熱量は常に一定であるという解析解の境界条件を満たすために、加熱開始より30分以降のデータを使用した。
【表6.4】
JP2016017773A_000032t.gif

【0084】
図6.15、6.16は実験ケース1に対する熱応答試験の計測結果である。図6.1
5はU字管21の出入口温度、平均温度、熱交換量である。加熱開始前の地盤温度は約17.4℃で安定しており、加熱開始から120分後には約33.49℃まで上昇した。熱交換量は、加熱開始から30分まではあまり安定しないが、それ以降ではほとんど一定で推移している。また、図6.16は2つの座標(座標1:x=0.03m、y=0.06m)、(座標2:x=0.03m、y=0.125m)での温度変化量の経時変化を示す。なお,加熱時間は120分で、各座標における深さ方向の温度測定結果の平均値をとっている。このグラフより、熱源から近い座標1(x=0.03m、y=0.06m)の方が遠い座標2(x=0.03m、y=0.125m)よりも温度変化量が大きくなることが分かる。

【0085】
実験ケース1の場合の、逆解析によるパラメータの同定結果を表6.5に示す。また、図6.17は、同定結果を用いて演算した温度変化量の計算値と実測値を示すグラフである。図6.17より、どの座標においても計算値と実測値は非常に良く一致していることが確認できる.各座標における同定結果を比較すると、縦熱分散率κはほぼ等しいものの、横熱分散率κと真流速は、座標1(x=0.03m、y=0.06m)の方が座標2(x=0.03m、y=0.125m)よりも小さく、等価体積熱容量(ρc)は座標1(x=0.03m、y=0.06m)の方が座標2(x=0.03m、y=0.125m)よりも大きい。なお、逆解析に用いたデータは88個で、データ1つ当たりの誤差を示すf/N値はそれぞれ6.30×10-3(K)、4.20×10-3(K)で非常に小さい。
【表6.5】
JP2016017773A_000033t.gif

【0086】
次に、図6.18、6.19は実験ケース2に対する熱応答試験の計測結果であり、図6.18はU字管21の出入口温度、平均温度、熱交換量である。実験ケース2では、上記の2つの座標に加えて、座標3(x=0.05m、y=0.05m、z=0.15m)についても測定した。加熱開始前の地盤温度は約17.18℃で安定しており、加熱開始から120分後には約33.14℃まで上昇した。熱交換量は、実験ケース1と同様に、加熱開始から30分以降ではほとんど一定で推移している。また、図6.19は3つの座標1(x=0.03m、y=0.06m)、座標2(x=0.03m、y=0.125m)、座標3(x=0.05m、y=0.05m)での温度変化量の経時変化を示す。なお、加熱時間は120分で、各座標における深さ方向の温度測定結果の平均値をとっている。この結果を実験ケース1の結果と比較すると、流速の大きいケース2の方が温度変化量が小さくなっている。これは、前章で述べたように、流速が大きくなると早く定常状態に達するようになるためである。

【0087】
実験ケース2の場合の、逆解析によるパラメータの同定結果を表6.6に示す。また、
図6.20は、同定結果を用いて演算した温度変化量の計算値と実測値を示すグラフである。図6.20より、どの座標においても計算値と実測値は非常に良く一致していることが確認できる。各座標における同定結果を比較すると、縦熱分散率κはほとんど等しいものの、横熱分散率κは座標の関係もあり値がばらついているが、真流速は、どの座標でもほとんど等しい。また、等価体積熱容量(ρc)の同定値は、3つの座標で、それぞれ、1005.04(W・hr/m/K)、924.731(W・hr/m/K)、997.675(W・hr/m/K)で、座標2(x=0.03m、y=0.125m)における値が一番小さくなっている。なお、逆解析に用いたデータは88個で、データ一つ当たりの誤差を示すf/N値はそれぞれ6.80×10-3(K)、6.70×10-3(K)2.31×10-3(K)で非常に小さい。
【表6.6】
JP2016017773A_000034t.gif

【0088】
ここで、逆解析によって得られた等価体積熱容量(ρc)の確からしさを検証するために,簡易試験を行って、実験的に等価体積熱容量(ρc)を求めた。簡易試験の手順は、概ね以下のようになる。断熱容器(魔法瓶)を用意し、これを、室内熱応答試験を行った恒温室内温度と同程度の水温、密度、比熱の脱気水で満たしておく。そして、室内熱応答試験装置1と同じく間隙率n=0.463となるように試料(硅砂4号)を充填した小型容器を加温し、これを断熱容器(魔法瓶)内に入れて密閉し、試料および水の初期温度、平衡温度等を計測する。計測結果から、次式(36)に示す熱量保存則に基づいて等価体積熱容量(ρc)を算出する。
JP2016017773A_000035t.gif

【0089】
本発明者らは、上記の簡易試験を3回行った。表6.7に、試験結果と、得られた等価体積熱容量(ρc)の値を示す。3回の試験結果の平均値は、862.79(W・hr/m/K)となった。
【表6.7】
JP2016017773A_000036t.gif

【0090】
図6.21に、室内熱応答試験の実験ケース1、2における等価体積熱容量(ρc)の同定結果と簡易試験より得られた等価体積熱容量(ρc)の実測値を示す。簡易試験より得られた実測値よりも本発明の逆解析より得られた同定値の方がやや大きく評価されており、座標2(x=0.03m、y=0.125m)における同定結果が最も実測値と近くなった。本来、等価体積熱容量(ρc)は地下水流動や座標等の影響を受けない固定の値を示すが、本逆解析方法より得られた同定結果はややばらついている。

【0091】
また、本発明者らは、実験ケース1、2の座標2(x=0.03m、y=0.125m)のデータに基づき、従来の解析方法(移流・分散現象を考慮しない解析方法:Kelvinの線源関数の解析解を用いる方法)で逆解析を行い、パラメータを同定した。その結果、等価体積熱容量(ρc)として、以下の値が得られた。
実験ケース1:1128.723(W・hr/m/K)
実験ケース2:1209.85(W・hr/m/K)
また、f/N値は、それぞれ以下のようになった。
実験ケース1:1.3×10-2
実験ケース2:1.22×10-2

【0092】
図6.21に示す、等価体積熱容量(ρc)の実測値と本発明の解析方法による同定値を、従来の解析方法による同定値とを比較すると、本発明の解析方法では、従来の解析方法よりも実測値に近い同定値が得られ、従来の解析方法よりも格段に精度良くパラメータを同定できることが確認できた。また、f/N値についても、本発明の解析方法では、従来の解析方法よりも小さくなることが確認できた。

【0093】
図6.22に、実験ケース1、2における真流速の同定結果と室内熱応答試験の結果より得られた真流速の実測値を示す。この図から、実測値よりも本発明の逆解析より得られた同定値の方がやや大きく評価されており、ケース毎で比較すると、実験ケース2の方が実験ケース1よりも実測値に近い値を示している。また、真流速においても等価体積熱容量(ρc)と同様に、座標2(x=0.03m、y=0.125m)における同定結果が最も実測値と近くなった。

【0094】
図6.23、6.24に、実験ケース1、2における縦熱分散率κ、横熱分散率κの同定結果を示す。どのケースにおいても縦熱分散率κの方が横熱分散率κより大きくなっている。したがって、一般的に縦熱分散率κに含まれる縦分散長αの方が、横熱分散率κに含まれる横分散長αよりも大きくなるという特徴が正確に評価されたといえる。縦熱分散率κは、同じ座標であれば地下水流速の大きい実験ケース2の方が実
験ケース1よりも大きくなる。図6.23をみると、座標2(x=0.03m、y=0.125m)では実験ケース2の方が実験ケース1よりも大きくなっているが、座標1(x=0.03m、y=0.06m)では実験ケース2の方が実験ケース1よりも小さくなっている。この原因は、微妙な実験誤差や同定誤差が原因であると考えられる。また、横熱分散係数も同様に、より流速の大きい実験ケース2の方が実験ケース1よりも大きくなるが、座標2(x=0.03m、y=0.125m)ではその傾向を示しているが、座標1(x=0.03m、y=0.06m)では逆の傾向が見られる。この原因についても、微妙な実験誤差や同定誤差が原因であると考えられる。

【0095】
図6.25、6.26に、各実験ケースにおける縦分散長α、横分散長αの算出結果を示す。縦分散長α、横分散長αは、パラメータの同定結果を用いて、本発明の解析理論で用いた式(3)より算出した。算出用の等価体積熱容量(ρc)、真流速は同定値を用い、熱伝導率λは硅砂4号完全飽和の傾斜0%における同定結果(2.466W/m/K)を用いた。一般的に、縦分散長αと横分散長αの関係はα≦αとなることが知られており、図6.24、6.25より、各実験ケースにおける縦分散長α、横分散長αを比較すると、この関係が満足されていることがわかる。

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

【0097】
6.テストケースによるパラメータ同定精度の検証
本発明者らは、本発明の移流・分散現象を考慮した地下水流動方程式の解析解の高次漸近解を用いてパラメータの逆解析を行う解析プログラムを開発した。そして、これを用いた場合のパラメータの同定精度の検証を、以下の手順で行った。
(1)表6.8に示す試験条件(テストケース)を与え、開発した解析プログラムを用いて水位降下量算出し、これを検証データとする。
(2)同定対象の3つのパラメータ(透水量係数T、貯留係数S、B)の任意の初期推定値を出発点として、これらのパラメータを逆解析によって同定する。
(3)同定したパラメータを用いて水位降下量の再現計算を行い、(1)で算出した検証データと比較することにより、その精度を確認する。
【表6.8】
JP2016017773A_000037t.gif

【0098】
表6.8に示す2つの試験条件のうち、テストケース1に対する検証結果を表6.9および図6.27に示す。表6.9は逆解析による同定結果であり、図6.27は逆解析より求めたパラメータを用いて水位降下量を再現計算した値と正確値との比較である。テストケース1では、正確値に対して初期推定値を2倍の大きさにして解析を行った。図6.27より、計算値と正確値は非常に良く一致している。また。表6.9に示すようにf/N値(データ1つ当たりの誤差)もほぼ0に等しく、精度の高い結果を得ることが出来た。
【表6.9】
JP2016017773A_000038t.gif

【0099】
テストケース2に対する検証結果を表6.10および図6.28に示す。表6.10は逆解析による同定結果であり、図6.28は逆解析より求めたパラメータを用いて水位降下量を再現計算した値と正確値との比較である。テストケース2では、設定値に対して初期推定値を1/2倍の大きさにして解析を行った。図6.28より、計算値と設定値は非常に良く一致しており、表6.10に示すようにf/N値(データ1つ当たりの誤差)もほぼ0に等しく、テストケース1と同様に精度の高い結果を得ることが出来た。以上より、本逆解析方法の妥当性を確認することが出来た。
【表6.10】
JP2016017773A_000039t.gif

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

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

【0102】
(試験条件)
上記のサイトにおいて、第1帯水層と第2帯水層の井戸能力を把握するため、層別揚水試験を行った。表6.11に各帯水層における揚水時間、揚水量、帯水層厚、揚水井と観測井との距離について示す。揚水は観測井の水位変化が見られなくなるまで行い、それぞれ180minと175minとなった。揚水量はどちらも417(L/min)である。また、帯水層厚は第1帯水層が32.7m、第2帯水層が21.2mである。
【表6.11】
JP2016017773A_000040t.gif

【0103】
揚水試験結果とパラメータ同定結果は、以下のようになった。
(1)第一帯水層
揚水井から23.5mの距離にある観測井の地下水位の変動を図6.31に示す。地下水位は試験開始1min後から減少を続け、120minで0.215の水位降下が確認されている。120min以降は水位の変化がほとんど見られなくなったため、180minで試験を終了している。

【0104】
表6.12に逆解析による同定結果を、図6.32に同定結果を用いて揚水試験を再現計算したときの観測井における水位降下量の計算値と実測値を示す。逆解析より得られた透水量係数Tは1.02(m/min)、貯留係数Sは3.59×10-4、Bは0.034で、この値を用いた計算値と実測値は非常によく一致している。透水量係数Tを帯水層厚32.7mで除して、透水係数Kは0.0312(m/min)となった。この値は一般的な透水係数Kの値である。また、f/N値は1.62×10-3(m)と小さく、精度の高い結果が得られた。なお、逆解析時に使用したデータは180個である。
【表6.12】
JP2016017773A_000041t.gif

【0105】
(2)第二帯水層
揚水井から24.5mの距離にある観測井の地下水位の変動を図6.33に示す。地下水位は試験開始1min後から減少を続け、140minで0.587mの水位降下が確認されている。140min以降は水位の変化がほとんど見られなくなったため、175minで試験を終了している。

【0106】
表6.13に逆解析による同定結果を、図6.34に同定結果を用いて揚水試験を再現計算したときの観測井における水位降下量の計算値と実測値を示す。逆解析より得られた透水量係数Tは0.459(m/min)、貯留係数Sは3.18×10-5、Bは4.84×10-5で、この値を用いた計算値と実測値は非常によく一致している。透水量係数Tを帯水層厚21.2mで除して、透水係数Kは0.0216(m/min)となった。この値は一般的な透水係数Kの値である。また、f/N値は5.32×10-3(m)と小さく、第一帯水層と同様に精度の高い結果が得られた。なお、逆解析時に使用したデータは175個である。以上のように、隣接帯水層からの漏水がある場合でも、精度よく解析を行えることが確認できた。
【表6.13】
JP2016017773A_000042t.gif

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

【0108】
(1)開発した解析解を用いて地下水流速の異なる3ケースで順解析を実施したところ、流速が小さく熱伝導現象が卓越するケースでは熱が円形に広がるのに対して、流速が大きくなるにつれて熱分散が卓越し、温度分布が楕円形になった。すなわち、本発明の解析方法は、地下水流速の熱移動への影響を適切に表現できる。
(2)同一条件下で熱移動と物質移動を比較すると、流速の小さいケースでは、熱伝導によって熱の方がより広範囲に移動するのに対して、流速が大きくなると地下水流動の影響を強く受ける物質の方が熱よりもより広域に移動する。
(3)開発した解析プログラムを室内熱応答試験データの解析に適用した結果、同定された等価体積熱容量(ρc)、真流速は実測値と近い値が得られた。また、従来のKelvinの線源関数の解析解に基づく解析方法で同定した場合よりも、格段に精度良くパラメータを同定できることが確認できた。従って、地下水流動があるサイトでの熱交換井の熱交換挙動を従来よりも精度良く予測することが可能になると期待できる。
(4)漏水性を伴う帯水層への揚水試験解析解(ハンタッシュ-ヤコブの式)を用いた逆解析方法を開発し、現場サイトで実施された揚水試験へ適用した。その結果、非常に精度の高いパラメータ同定を行えることが確認できた。従って、地下水流動があるサイトでの地下水位の挙動を従来よりも精度良く予測することが可能になると期待できる。
【符号の説明】
【0109】
1…室内熱応答試験装置
10…実験地盤
11…土槽
12…水位設定部
12A…上流側水位調節部
12B…下流側水位調節部
13…浸透量調節部
14A…給水容器
14B…上流側水位調整容器
15…下流側水位調整容器
16…試料充填部
17A、17B…貯水部
18…透水性スリット
20…熱応答実験装置
21…U字管
22…循環流路
23…循環ポンプ
24…電気ヒーター
25…給水容器
26a、26b…温度センサー
27…温度センサー
28…流量計
29…データロガー
30…パーソナルコンピューター
31…室温計
図面
【図1】
0
【図2.1】
1
【図2.2】
2
【図2.3】
3
【図2.4】
4
【図3.1】
5
【図3.2】
6
【図4.1】
7
【図4.2】
8
【図4.3】
9
【図4.4】
10
【図4.5】
11
【図5.1】
12
【図6.1】
13
【図6.2】
14
【図6.3】
15
【図6.4】
16
【図6.5】
17
【図6.6】
18
【図6.7】
19
【図6.8】
20
【図6.9】
21
【図6.10】
22
【図6.15】
23
【図6.16】
24
【図6.17】
25
【図6.18】
26
【図6.19】
27
【図6.20】
28
【図6.21】
29
【図6.22】
30
【図6.23】
31
【図6.24】
32
【図6.25】
33
【図6.26】
34
【図6.27】
35
【図6.28】
36
【図6.31】
37
【図6.32】
38
【図6.33】
39
【図6.34】
40
【図5.2】
41
【図6.11】
42
【図6.12】
43
【図6.13】
44
【図6.14】
45
【図6.29】
46
【図6.30】
47