計測自動制御学会束北支部第232回研究集会（2006．11．24）
資料番号232−12

多孔質休内水流ベクトルの計測理論

TheoryfortheMeasurcmc11tOfVelocityVectorofWater
FlowinginPorous〜Iaterial
原道宏MichihiroHara
岩手大学
Iwa妃U11iverslty
キーワード：流速ベクトル（veloc圧yvector）．測温位置誤差（errorsin川eaSUringpolnts）．熟トレーサー法（heattracer method），多孔質体（porous materral），水流（water fovJ）

連絡先： 〒020−8550盛岡市卜円3−］．8・8 岩手大学農学部農林環境科学科リサイクル生物生産「二学喜諾座環境制御丁学研究豊 原
道象 TEL＆FAX：（q19）621−6125，E−nlail：些地＠通塾kユ！．往旦血

1．はじめに土壌のような多孔質体の中を流れる水など流体の速度ベクトルは土壌物理学、環境汚染、自然災害や農業などの分野において最も関心を持たれている事項の一つである、1それゆえ古くから多くの研究者が多様な方法でこの測定について研究している（Dirksen1972，1974；Byrneら1967，1968など）多孔質内の水流速計測への熱トレーサーの使用については、Renら（2000）は数値解法を計測に応用し、Wangら（2002）は解析解の利用を提案し、Ochsner ら（2005）は解析解により計測した。遠藤一原（2003a）は熱トレーサーに5極子センサ（QPHPセンサ）を導入し水流速ベクトル計測理論式を用いて計測精度を検討し、遠藤・原（2003b）はQP川⊃センサを用いて傾斜成層において計測し、他の計測値と比較した√

本論文の目的は、熱トレーサー法により多孔質内の水流速ベクトル計測する場合において，センサの設置位置に誤差があることを考慮に入れコ その対処方法をも含む計測理論を開発することである〔
2．材料および方法

無限揚における等温・均質・等方性の多孔質体を扱う√ 多孔質体の熱特性には，温度伝導率K［m2／s】，体積熱容量G［J／（Km3）】，熱伝導率g［J／（Ⅹms）］があり，多孔質体中の水流ベクトルuの旦ぷg成分は（ux，叫，U2）［1−1／s］であるが、それらは．いずれも，もちろん，未知である

熱トレけサー法を実施するため、点熱源または直線熱源（以下、線熱源と略称）を用いる。点熱源の位置を座標原点，また練熟源の位置をz軸とする座標軸上の，熱源に対して対称の位置に，温度センサを配置するし すなわ
ち，点熱源に対しては耶3座標軸卜の6箇所、線熱源に対してはギア2座標軸卜の4箇所に配置する。こうすることで‥点熱源＋6測温点により流速ベクトルの今成分（ux．恥払）が、また，練熱源＋4測温により2成分（ux，町）が，それぞれ計測されるものと期待される、なお− 卜記のように熱特性は未知と仮定しているので，必要な熱特性は同1

時に計測するものとする∴ すなわちJ計測理論は，必要な熱特性および流速ベクトルのすべてを求められるものでなければならない
測温点の位置は，いずれも，座標軸⊥の，原点からの距離がj・［m］とする。ただしプ 厳密には測点ごとにこの距離は異なり，厳密にrではなく，誤差が不可避と考えられる（：測温センサの設置位置誤差に対する対応は，本理論において必須■と考えられるので，それは本論において述べるものとする′
熱源（卜記のょうに，点熱源または線熱源）は期間血［s］の間，一定強度で発熱するものとする。熱源がこの期間に発生した熱量の総量を ￠（単位は，点熱源では［J］，線熱源では［J／m】）とする「時間軸の原点を発熱期間のL卜点にとる、．すなわち∫ 発熱期間は・ね／2＜t＜血／2である＿
温度測定は連続的に行われるものとしク 特に，温度卜昇が最高温度に達した時礼占n［slおよび最高卜昇温度n［degC］は計測結果から割り出されるヰ〕のとする．
3．結果

3．1測温センサの放置位置に誤差が無い場合の理論
ここでは3次元計測を中心に記す1・ただし，下記のpを3から2にすることにより線熱源の場合について対応できる。
瞬間点熱源￠【J］による温度卜昇は次の通りである。
V3【Ⅹ′y′ヱ′t］＝＿＿＿「 （Ⅹ−職t）2十（y一旬t）2十（ヱー鴨t）2

包（47TJくt）3／2cv1L4xt

期間一血／2く医ね／2の間J・定強ノ蜃￠／血［J／sjを発熱した後の温度卜昇は次の通りである二￣
T3［x′y′Z′t】＝∴∴Q3／to（X】tkて）2＋（y−uyて）2＋（ヱー職工）2（47rXT）3／2cv

4xて

1dT（t＞蔓）

温度上昇耶は次のように書き直される。

Q3／to

ワ潮y′‥］＝叩［器十崇血e r＝∨㌔＋㌔＋Z2訂du＝㌔＋u2t2（4刀∵X：て）3／2cv

4x・て］dて職2・u㌔＋uz2

凰酎こ対して対称の位置にある2点i旦．I；∂，‡笥署一カにおける温度上昇の差の和に対する比は、次のように双曲線正接となるが，その値は時間に関わらず一定である
ただし，双曲線正接tanh［x］は（ex−e￣う／（e・Y＋e−⇒であるし℡3［Ⅹ′y′Z′t］一丁3卜Ⅹ′−y′−Z′t］

T3［Ⅹ′y′Z′t］十T3卜Ⅹ，−y′−Z′t】

L2x 2x 2x′つ．＿」（LJx，Uy、Uユ）を求めるには次式がふさわしい ただし，tanll￣リカはtanh［元 の逆関数である
T3［x′y′Z′t］一丁3トX′−y′−Z′t］＋
2K 2x 2x彗空塑・竺ヱ＝由nh−1［ノー ＿＿′ヽ ＿＿′ヽ ＿＿LT3【x′y′Z′t］＋T3トⅩ′−y′−ヱ′t］

まず，無次元数り哀‥彷，月；）を次のようiこ定義する
粘≡・Hy≡′粘≡rlユ空

これらは，測温点における温度上昇7もから，次のように求められる＝
ワ3【r′0′0′t】一丁3卜r′0′0′t］‰≡票＝t訂叫T3［ご′0′0′セコ＋で3卜r′0′0′t】

T3［0′r′0′t］−℡3［0′−r′0．七】Hy≡芸＝血￣1［T3［0′r′0，t］＋T3［0′−r′0′t］
2

2x L

T3［0，0′r′t］−で3【0′0′−r′t］
T3［0′0′r′t］十で3［0′0′−r′t］

しかし「温度伝導率にが未知であるので】その求解を追究しなければならない、無次元数Hを次のように定義する。
ruH≡ −−−

2x

こ■の〃の値は，卜記の月立，．払，月妄から次のように計算される

H＝、／軋こ．Hy2．H＝ニ方、対称点における温度卜昇の相乗平均G［考ガろf］は次のように計算される「句【Ⅹ′y′Z′t】≡イ・で3【Ⅹ，y′Z′t］T3トⅩ′−y′一三′ヒ】

これは，伝熱理論においては次のようになる、t＋句／2 Q3／to

q3［Ⅹ′y′Z′七］＝叩卜

一旬／2（4刀∵六：T）3／2cv

血ご＝イ㌔＋㌔＋Z2弧du＝イ匝

職2＋uy2＋鴨2

ニこまでは点熱源を対象にしてきたが、次のようにpを設定すれば線熱源の場合に移行できる。

旬／2 Qp／亡0

Gp［r′t】＝0／2（4刀・入：て）針2cv

r r2＋u2t2

e呼卜二
L 4六：て二こ亡

forapointheat80tmeattheorigin，P＝3，r＝∨斎訴ardu＝イ両示′andforalineheatsourceonz−aXis，P＝2，r＝√百石乏ard u＝イ嘉二轟．
さて）にの求解を進めよう
対称点における温度の相乗平均の時間ピーク値は次式を満たす＝
dGp【r′tユ8七＝0 止七＝tm訂dち′m＝年［r′tm】

すなわち次の通りである【Qp／b㌔．㌔t2，、tm＋句／2（4刀・Xt）針2cv

L 4x亡

＝0

七m−tO／2

したがって，次式が成り立つ（つ1＿．＿「（tm・号）P／2 」

工2・u2（tm・空）21

1

J「

4x（七m・号）（七m一号）p′2

4りtm一号）

これは，変形すれば〕次のようになる㌔・ぜ（七m一号）2 工2・㌔（七m・4x（七m一号）

4可tm・箸）…上‥：二二‡ 2‥

ここで，ノ刀＝（血／2）化nおよび腸1n［（1＋ノ刀）／（い′刀）］／2を定義し導入すると、次のようになる
3

㌔＋㌔七m2（1−m）2 ㌔＋㌔tm2（1＋m）2＝PM

4xtm（1＋m）

4x・tm（1−m）

これを占。（ピーク時刻）の方程式として整理すれば，次のようになる

r22xMptm

tm21㌔＋＝0l−m2m（1）

ここで，先に定義した関係：
2Hu＝ ＋ Kr

により，未知数をェ〜からKに変えれば，Kついての方程式を次のように得る二2 Mpr2r42鮮mtm 4H2（1−m2）tm2（2）無縁根（に＜0）を捨てれば，温度伝導率Kの解は次のようになる

ここで，M＝1n［（1＋m）／（1−m）］／2であるのでJこれは次のように書ける

Ic22ptm

2m

1
X ＝トm21n【壬慧］

叫壬莞】）2（苦）2

温度伝導率Kが求められたので流速uは次のように求められる
2H

1u＝ 一一一一 X ＝r2m

JII冊

1−が叫慧］（3b）

流速の各成分は，次のようになる

HY
Ⅰも射u′uy＝u′鴨＝uHヱ

（4）

速度ベクトルの成分（ux，Uy、Uz）がこのように求められた
体積熱容量αは次式により求められる。すなわち，積分の計算が必要である

Cv＝

Qp／to
（4刀∵K）P／2Gp′m‥∴二1r r2十l㌔て2げ／2 e坪卜L」

4xて（5）

なお、mの値が1より小さく十分0に近い場合には，2m／1n［（1＋m）ノh−m）】は（1−m2）1／3に非常に近いので次の近似式が成り育つr、この近似式の誤差は0＜nlく0．43で0．1％以下、0＜mく0．7でl％以下である．ノ2（1−m2）2／3

工2●）2p七m4−6。）

Jl≡（J（6b）さらに、Hが1に比べて小さい場合には次のように近似されるrl工2）2）2ptm（新「嘉十）

u＝ −−．K 之r折諒（妊諒＋（憲）2）Ptm

3・2 測温センサの設置位置に誤差がある場合
ここでは、線熱源の場合について述べる
まず，瞬間熱源の場合からはじめよう∴温度上昇は次の通りである（Ⅹ一職t）2＋（y−uyt）2圃V2［x′y′t］＝4xt

簡単化のため次のように定義する〔つp＝（Ⅹ，y）弧du≡イ冠■▲㍍；盲すると、温度上昇は次のように書ける

醐′t］＝甲［芸＋慧］「Ⅹ2＋㌔＋u2t2

勉
4刀・XtCvL

4x亡有限時間（一缶／2＜K血／2）における熱ぺ′レス発生後の温度上昇は次のようになる㌔＋㌔＋㌔て2 工

七＋句／2 Q2／to

T2［P′七］＝甲［芸・崇］上一

旬／2 4刀 六：tCv

4六：て］dり亡≧害）

測温点の座標は次の通り（1）p一＝（xり y一）＝（rり 0）and P2＝（Ⅹ2，y2）＝（，r2，0）（2）p。＝（Ⅹ。，y3）＝（r3，0）andP4＝（Ⅹ。，y4）＝（−r4，0）
測温点の原点からの距離の平均をrとし、平均との相対誤差をekと置く
rl＋r2＋ご3＋r4r≡4（1≦k≦4）′∑ek＝O′ek≡klすなわち、

rl＝r（l＋eL），r2＝r（1＋e2），r3＝r（1＋e3）andr4＝r（1＋e4）．
であるこのとき、各測温点における温度卜昇は次のとおりである1
恥［Pl′tい＝叩［lkr（1＋q）
2ズ℃［p2′t】＝＝e坪［

恥［p3′七】＝＝e呼【取【P4′t］＝＝e呼卜47rXてCv

ー1もご（1＋句）勉／to2x4刀1（てCvuyr（1十e3）

包／to2x47一光亡Cvuyr（1＋e4）
2x（1＋q））2＋㌦て2包／to包／to47rXT Cv

ここで、tが大きい場合には次のような近似が成り立つことを用いる
5叩卜］正二4大：て（ご（1＋句））2＋㌔て2e坪卜

4火：て］こ⊂（ご（1＋e3））2＋J㌧2巨l・e4））2＋ぜて2
4大：て］亡て叫′2叩卜慧］

t＋tO／2印［一旦浩三］て2すると、各測温点ペアにおける昇温差の昇温和に対する比は次のように近似される‖で2【Pl′t］一丁2［P2′七］（el一句）㌔だ叫（1・旦ヂ‥粘一T2［Pl′t］＋ワ2［P2′亡］

4Kt℡2［p3′t］一丁2［P4′t］（e3一句）㌔記叫（1＋誓叫y−T2［P3′t］＋T2【P4，t】

4xt

無次元数軋とHyは次のように解ける軋≡＝

T2［pl，t］一丁2［P2′t］】／（1・里芋））で2【pl′t］＋T2［p2′t］

且y≡芸＝碧T2［p3′七】−で2［P4′七］］／（1＋誓りT2［P3′t］＋で2［P4．七】

これらの計算を実施するにはekの値を知っている必要がある。相対誤差ekは各測温瓜こおける昇温のピーク州

現時刻tmk（1≦k≦4）と次のような近似関係がある、この関係を用いてekの値が求められる無次元数Hが小さい場合にはH＝0とおいても誤差は小さい。例えば，0≦H≦0．4であれば，左辺は誤差1−％以内でekに等しい．
tmk−

平凡十bセ＋叫＋b14
4

ekト芸−0叫＝三bq十血色十b喝十叫

fork＝1′ 2′ 3′ 443．3 数値例
物性定数および流速を次のように設定し，数値計算した（「に＝10￣6m2／S′Cv＝2＊106］／m3，uX＝3＊10．5m／S，uY＝4＊10→5m／S′U＝5＊10￣5m／S′r＝0・01m，Q2＝1000〕／m・
測温点の言呉差iこついては、つぎのように検討した。
（A）誤差の無い場合（B）誤差のある場合誤差は，−20％ と 20％の範囲とした（BトB3）．

発熱時間（ち）については、AとBlでは10s，B2では20s，B3では30sとした。、数値例の結果を図（Fig．1−5）と表（Tablel−8）に示すっ
4．おわりに測温位置に誤差のある場合を含め、一定時間一定強度の発熱により牛起される温度上昇を熱源からの対称位置における測定値から求める方法を示した√ 測温点の位置が厳密に対称にならないと考えられるので，それへの対処が必要である。本論文で示した対処方法は，数値例で例示されるように，有効なものであると判断される′、参考文献

l）Bristow，K．L．，］．R．Bilskie，G．］．KluitenbergandR．Horton二Comparisonoftechniquesforextractings
propertiesfromdual十PrObeheat−Pulsedata，SoilSci・160、1／7（1995）

2）Byrne，GF‥J．E．DrummondandC．W．Rose：Asensorforwaternuxinsoil・ Pointsoul．Ce 、instrument，WaterResoul−・
Res．3、1073／1078（1967）

3）Byrne，GF，J．E．DrummondandC．W．Rose：Asensorforwaterfluxinsoil．2．．．LinesourcelへinstrResour．Res．4、607／611（1968）
64）Campbel，GS．1，C・CalissendorffandJ・H・Wiuia［TIS：ProbeformeasuringsoilspeCificheatusingaheatrpulsemethod，SoilSci．Soc．Am．］．55，291／293（1991）

5）Carslaw，H．S．andJ．C．Jaeger：Conductionofhea【insolids．211ded．．OxfordUniv．Press，London．266／26

6）Dirksen，C・：Aversatilesoilwaterfluxmeter，Proc・2ndSymp・OnFundamentalsofTransportPhel】OmenainPorousMedia，2，1AHR，ISSS，Canada，425／442，（1972）
7）Dirksen，C∴Fieldtestofsoilwaterfluxmeters．Trans．ASAE17，1038／1042（1974）

8）Melville，J・G，F．J．MoIzandO・Gtlven＝Laborat（）ryinvestigationandanalysisofaground−Waterflowmeter、GroundⅥrater23、486／495（1985）

9）Ochsner，TysonE・，RobertHorton，GerardJ．Kluitenberg，andQuanjiuWang：EvaluationoftheHeatPulseRatio
MethodforMeasuringSoilWaterFlux，SoilSci．Soc．Am．J・69：757−765（2005）．

10）Ren，T．，GJ．Kluitenberg and R Horton：Determinh〕g SOilwaterflux and pore water velocityby a heat pu】setechnique，SoilSci・Soc．Am．J．64，552／560（2000）

1し）Sisodia，S▲，0．J・Helweg：Modelingtheheatsenseilownleter．J．1rrig−Drain．Eng．124（3），148／151（1998）12）Wang，Q．，乃7SOn E．Ochsner andRobert Horton：Mathematicalarlalysis ofheatpulse signals forsoilwaterf
determination，WaterResourcesResearch、Vol．38，No．6，27−1／27r9．（20O2）

13）遠藤明・原道宏：有極子熱線パルス法を用いた2次元流れの場における砂質土壌の熱前線伝播速度ベクト
ルおよび熱特性の同時計測法およひその応用、計測自動制御学会産業論文集、ヽbl．2，No．1ユ，88／95（ユ003a）
14）遠藤明・原道宏：QPHPセンサを用し、た多重成層砂中の流体速度ベクトルと熱物性の同時測定、計測自動制御学会産業論文集、＼加．3，No．2，12／21（2003b）7

r＼＼＼＼こ こ二＝；≡40208060Fig．1T2［Pl′t］′T2［P2，t］and（T2［Pl，t］T2［P2，t］）1／2
er，lowerandmiddle′reSectivel

2040

Fig．3T2［Pl，t］′T2［P2′t］and（T2［Pl，t］T2［P2，t］）l／2
htandmiddle，「eS

ect：ive［0′0∩︶4︵エnU204C6080

Fig．5（T2［Pl′t］T2［P2／t］）1／2，（T2［P3，t］T2［P．′t］）1／2′and（T，［Pl，t］T2［P2，t］T2［P3′t］T，［P4，t］））1／4
Lowe「′∪e「′andmiddle，88060

Fig，2T2［P3，t］′T2［P．′t］and（T2［P3，t］T2［P．′t］）1J2,ectivel
er′Lowerandmiddle

20806040J】06080Fig，4T2［P3，t］′T2［P．′t］and（T2［P3，t］T2［P．，t］）1／2
htandmiddle，reS ectivel

＼408060Fig．7Tanh￣1［（T2［Pl′t］−T，［P2，t］）／（T2［Pl，t］＋T2［P，，t］）］and
Tanh 1［（T2［P3，叶T2［P4，t］）／（T2［P3，t］＋T2［P．，t］）］vsl／t
ectiveler′「eSLowerand uFig・6（T2［Pl′t］一丁2［P2，t］）／（T2［PIIt］＋T2［P2′t］）and（T2［P3，t］−T2［P4，t］）／（T2［P3，t］＋T2［P．′t］）vst
Lowerand uectivele「′「eS

0・2〔・ ＼ −−●、
￣ ＼ ＼、●−→−−叫−
0．05

0．1

0，15

0．2

−−−1トー、

−一寸トー

0．25

Fig．8．P10tSOfeachek（1≦k≦4）vs．m2．
Dotsarefo「Bl′B2andB3．
Linesare「eg「ess70nOFthedots．

9

Tablel Setvalues forthe measurements．

Table2 Measuredva山esforeach setvarueofto．
tm2

t叩．

tm3

tm4 ≡Ave［tmk］／sel［t。］

e2［t。］

25．28125．2810．00000

0．00000

35．68820．805e3［t。］

e．［t。］≧sum［ek］

25．28125．281i25．281

0，00000

37．03223．00016．84530．247i25．896−0．17476

0．08400；0．00000

19．45フ31．818i27．827

Table3 Calculation forC、′．
C、

／degC

／（106J／m3）

／degC
e［Cv］

0・57250 董e［∂］2．00000．5フ2500．000％

Bl

0．556750．00000
1．98870．55991−0．56卵昌一0．00270

1．9797
B20，54738

0．55299−1．015ウら−0．00485

1，9738B3／s／degC

25．281 0．57250

25．811 0．55675

27．475

B2120

A

G2．．¶

0．00000 ≡ 0．00000

Bl 100．18906 【0．09830

tm

0．533640．54071−1．3089ム

ー0．00625

100．54738

TabFe4 CaEcu［ationofH／一HyandH・

0．96066

Table5a CalcuEatedvaEuesof FC，u，ux、uyVVitheq・（3）・m〟uxUuyC＝（to／2 ／（10￣6m2／s）

A 0．19778

Bl 0．19372

B2 0．36397／（105m／s）

1．00000

00

3．00000．0009昌

0％0．000％

0．9786495

2．99834．00000．000％
3．8116−2．136予も−3．00鋸も ー0．0559も ー4．711ウノら

0．98573363．1755−1．4279もB3 0．49990

3．59098％ 5．850％ 一10．23ワも44．7551％

一4．898％ 14．27％＋＋−17．62サム

3．4282

3．2952

Table5b Calculatedvaluesof／（，uTuン，uyWitheq・（7）
リ〝

uyリズ

lm＝（t。／2）／tm毒／（10−6m2／s）1
A 0．19778

0．99973／（10L5m／s）

4，99862．9992

3．9989−0．027サム−0．0279昌 一0．027ウ竜 一0．027！昌

0．9783983

2．99763．8■106

日1 0．19372
5％ −0．082％ −4．7369も−2．1619′もヨ ‡4・7904 3・1734 3・5885
B2 0．36397
2％ 5．780％ …10．29％−1．493％B3 0．49990

0．99363H55

3．4212

3．2886−0．637％ −5．091％ 14．04うも 一17．79％

11

Table6 Theindexeobtainedwitheq．（8）andtheh2qextrapolationmethod．

TabLe7 TheHr−HyandHusingthe（mprovedevalueswiththeh2qextrapo］ationinTable6・

仙1 宰 0−17498 0．95257

1．04743 ，−・…・叫・

0．16005

0．95257

Table8a ThevaLuesofu，uゝ．uyVV；theq・（3）usingimprovedHvaEuesinTable7

Bl

64 〜 4，8524 2．9925 3．8198

0．193726％‖巨2■953％ −0・25用 一4・506％

73 蕾 4．8003 3．1510 3．6214
B20．3639了
7％ i−3・993％ 5・0329もー9・469も山 64 ぎ 4．7561 3．3804 3．3458
B30．4u−0．4369も⊆一4．877％12．6鍋 −16．36％

Table8b ThevaIuesofu，ux，uyWitheq・（7）usingimprovedHvaluesinTable7
疋muニ（㌔／2）／t仰 ≡／（106m2／s）

Bl0．19372【／（10￣5m／s）−2．161％0．3639739右B3uy【0．9783907B2uゝ0．49990 0．99363≒ 4．7971 3．1489 3．6190−4．057％ 4．962チ昌一9．529后

4．7465 3．3735 3．3390−5．06酌も 12．459竜 一16．539占12