JP3926571B2 - Method and apparatus for measuring temperature conductivity etc. of solid - Google Patents

Description

又は前記固体が半無限体の場合の表面温度の近似式
【数１０】

を得る手段と、得た上述の近似式から推測される複数の箇所の表面温度が一致することから温度伝導率を求める手段とを備え、
前記固体の温度伝導率を計測することを特徴とする測定装置である。
また、液体あるいは気体内の測定点での質量濃度変化を得るための測定手段と、測定さ れた質量濃度変化から、近似式の係数を決定するための手段と、前記近似式をラプラス変換するための手段と、拡散方程式のラプラス変換式の解を、ラプラス変換された近似式の係数を用いて求める手段と、求めた解をｓ＝０の近傍で展開する手段と、展開された解のラプラスの逆変換を求めて、前記気体あるいは液体が有限体の場合の質量濃度の近似式
【数１１】

又は前記気体あるいは液体が半無限体の場合の質量濃度の近似式
【数１２】

を得る手段と、得た上述の近似式から推測される複数の箇所の質量濃度が一致することから拡散係数を求める手段とを備え、前記気体あるいは液体の拡散係数を計測することを特徴とする測定装置である。
その上、同様の測定方法や、この測定方法のプログラムおよびプログラムを格納した記録媒体も本発明である。
【０００４】
【発明の実施の形態】
本発明は、ラプラス変換法を用いて一次元逆問題を解く際に、測定された温度変化を熱伝導方程式の基本解の一部である時間の平方根の多項式で近似させ、温度変化の時間遅れを考慮することにより、計算を発散させず、予測精度を高くしようとするものである。
図１は、本発明の手法の概略を説明するフローチャートである。図１において、通常、微分方程式で表される熱伝導方程式をラプラス変換する（Ｓ１１０）。一方、物体内の適当な測定点での温度変化から近似式の各係数を決定する（Ｓ１８０）。この近似式としては、時間遅れを考慮した時間の平方根の多項式を持つ近似式を用いるとよい。近似式をラプラス変換して（Ｓ１９０）、ラプラス変換の係数を決定する。ラプラス変換をした熱伝導の解を近似式の係数を用いて与える（Ｓ１２０）。ラプラス変換形（積分形）で与えられた熱伝導解をｓ＝０近傍で展開する（Ｓ１４０）。展開した形式でラプラス逆変換を実行して、表面温度等を求める（Ｓ１６０）。
【０００５】
上述の表面温度等を求める手法を詳しく説明する。まず、以下の説明で使用する主な記号を以下に示す。
ｆ₁（τ），ｆ₂（τ）：測定点ξ₁，ξ₂での無次元温度
ｋ：温度上昇率（Ｔ／Ｔ₀＝ｋτあるいは２−ｋτ）
Ｌ：代表長さ
ｍｉｎ（θ）：測定器の最小感度あるいは有効数字の最小値
Ｎ：近似式の次数
Ｎ_sf：有効数字のけた
ｐ＝√ｓ
ｑ：熱流束
ｓ：ラプラス演算子（＝ｐ²）
ｔ：時間
ｘ：ｘ座標
Ｔ：温度
Ｔ₀：代表温度
ａ：温度伝導率
θ：無次元温度
λ：熱伝導率
φ：無次元熱流束
バーθ：無次元温度θのラプラス変換温度（像関数）
バーΦ：無次元熱流束φのラプラス変換熱流束（像関数）
ξ：無次元距離＝ｘ／Ｌ，ξ_l＜ξ₂
σ：標準偏差
τ：無次元時間（フーリエ数＝ａｔ／Ｌ²）
τ₁：最小予測可能時刻
τ_i ^*：無次元時間遅れ
Ｔｗ₁，Ｔｗ₂：測定点ｘ１，ｘ２からの推定表面温度
Ｔ_exact：厳密解温度
【数１３】

【０００６】
＜一次元熱伝導逆問題の定式化＞
物性値一定の一次元非定常熱伝導方程式は無次元形で次のように表される。
【数１４】

初期温度が一様（θ＝0としても一般性を失わないので、これ以降θ＝0とする）のとき、式（１）をラプラス変換する（Ｓ１１０）と
【数１５】

式（2）の一般解は容易に得られ、次式となる。
【数１６】

ここでＡ，Ｂは積分定数である。
【０００７】
（有限平板に対する解）
有限平板の逆問題を取扱う場合、少なくとも物体内の２点の温度が与えられないと式（3）は閉じた式とならない。いま、異なる２点（ξ＝ξ_l，ξ₂，ξ₂＞ξ₁＞０）での温度が次式で与えられていたとする。
【数１７】

未知定数Ａ，Ｂを式（4）から決定し、それを式（3）に代入するとバーθ（ξ，S）は
【数１８】

また、熱流束Φ［＝ｑ／（λＴ₀／Ｌ）＝−∂θ／∂ξ］は
【数１９】

となる。
逆問題では、一般にξ₁≦ξ≦ξ₂の範囲に含まれない点の温度や熱流束を予測することが目的となる。特に、表面（ξ＝０）での温度や熱流束を決定することが主目的となることが多い。そこで、ξ＝０を式（５），（６）に代入すると表面温度と表面熱流束はそれぞれ次式となる。
【数２０】

そして、式（７），（８）を逆ラプラス変換の定義（式（９）で与えられる）に従って実行すれば目的とする表面温度と表面熱流束が得られる。
【数２１】

【０００８】
＜測定温度の近似式＞
式（９）の積分を実行するためには、バーθ_w（ｓ）やバーΦ_w（ｓ）に含まれる既知関数、バーｆ_n（ｓ）、ｎ＝１，２の具体的な関数形を前もって与えておく必要がある。ところで、バーｆ_n（ｓ）、ｎ＝１，２は、ξ＝ξ₁，ξ₂での温度変化を示す関数ｆ_n（τ）をラプラス変換した関数であるから、ξ＝ξ_l，ξ₂での温度変化の測定から定まっている。そして、各点での温度変化をそれぞれ次の三つの関数形で近似されると仮定する（Ｓ１８０）。
【数２２】

ここで、係数ａ⁽ⁱ⁾ _k,nは、各測定点の温度の時間変化から例えば最小二乗近似によって決定される係数である。Ｎは近似式の次数を示す。
三つの関数形は、それぞれ次の特徴をもっている。式（１０ａ）は、温度変化を無作為に近似したものである。式（１０ｂ）は、初期条件ｆ_n（０）＝０を満足する近似関数である。また式（１０ｃ）は、測定点までの温度伝搬の時間遅れτ_n ^*を考慮した近似式である。この時間遅れτ_n ^*は、ｅｒｆｃ（ξ_n／２√τ_n ^*）＝ｍｉｎ（θ）で与えられる。
ところで、式（１０ａ），（１０ｂ），（１０ｃ）の３式を比較すると、式（１０ｂ）ではａ⁽²⁾ _0,n≡０また式（１０ｃ）ではτ_n ^*≠０となっていることを考慮すると、基本となる関数形はほぼ同一となっている。したがって、これ以降式（１０ｃ）を用いてのみ表記（これ以降、関数および係数の上付き添字を省略）することにする。
また、式（１０）の係数にガンマ関数を乗じた係数ｂ_k.i＝ａ_k.iΓ（ｋ／２＋１）ｉ＝１，２を用いると、式の表記上の煩雑さが避けられる。そして、式（１０）のラプラス変換は、次式となる（Ｓ１９０）。
【数２３】

ところで、近似関数として時間の半値多項式を選定した理由は、一次元非定常熱伝導の解が、一般にＴ＝ｆ（ｘ／√（ａｔ））という関数形で与えられるためである。
【０００９】
＜近似式の係数の求め方＞
式（１０）の係数ａ_k.n（あるいはｂ_k.n）は、例えば各測定点での温度変化を最小二乗近似することによって決定されているが、使用される温度変化の有効けた数（実験では、温度の測定精度）に依存している。そこで、温度の有効けた数をＮ_sf＝２から４けたまで変化した場合（厳密解の温度を、θ（ξ_i，τ）＝Ｉｎｔ｛θ_exact（ξ_i，τ）ｘ１０^Nsf｝／１０^Nsfである有効けた数の温度とする）と、３けために外乱（分散σ＝１）がある場合（温度は、θ（ξ_n，τ）＝θ_exact（ξ_n，τ）＋０．００５ε（σ＝１），ｎ＝１，２）について検討した。ところで、熱電対での温度測定の精度は、たかだか３けた程度であることから、３けために外乱（分散σ＝１）がある場合について検討すれば、精度的には十分であろう。図２は、有限平板で表面温度がθ（０，τ）＝θ（１，τ）＝１に突然なったとき、外乱があるときの温度変化と近似曲線を示す。近似式は、いずれも温度変化を最小二乗近似することによって決定されている。
図２から、式（１０ａ），（１０ｃ）のいずれの近似式（式（１０ｂ）の近似は、式（１０ａ）とほぼ同様であったので省略）もその次数を５〜７まで増加するとよりよい近似を与えているが、次数Ｎをそれ以上増加しても近似式の精度の改善は期待できなくなる。したがって、近似式の次数は、たかだかＮ＝７もとれば十分であろう。また、近似式の予測精度は、いずれの式を用いても、次数ＮがＮ＝５〜７ではほとんど差違がみられなかった。したがって、近似式自身については、Ｎ＝５〜７程度であればどの近似式を用いても問題ないことになる。
有効けた数以下を切り捨てた場合と正規外乱を重ね合わせた場合の予測式の精度を比較すると正規外乱が重畳された場合のほうが、同一条件下ではつねに近似式の予測精度が悪くなることがわかった。この結果、得られた近似式を用いて決定される逆問題解の推定も正規外乱が重畳された場合のほうが当然悪くなる（結果は省略）。ここでは、予測精度が最も悪くなる正規外乱が重畳された場合の近似式を用いた逆問題解についてのみ検討することにする。
【００１０】
＜逆ラプラス変換とその解の性質＞
式（１１）を式（７）あるいは式（８）に代入し（Ｓ１２０）、さらにそれを式（９）に代入すると
【数２４】

式（１２），（１３）の逆変換（留数計算）を行うと表面温度と熱流束を厳密に求めることができる。
式（１２），（１３）の留数は、ｓの指数が整数の場合（式（１２）ではｋの値が偶数、一方、式（１３）ではｋの値が奇数）には容易に留数の和を計算できるが、ｓの指数が整数でない場合には、ｓ＝０に分枝をもつために初等関数で与えられず、級数展開式となる。ところで、積分を行うことによって求められた解は、ｓ＝０とそれ以外（ｓ≠０）から求められる解の集合に分割することができる。さらに、ｓ≠０の解集合は、時間とともに減衰する性質をもち、その基本形は、ｅｘｐ（−ｎ²τ）／ｎ^k（ｋ≧１）であるので、無次元時間（フーリエ数，τ）およびｎの値が大きくなると急速に零に収束する性質をもった関数の集合になっている。また、ｋの値が大きくなるほど、零への収束も早くなる。しかし、τが小さくなると収束性が劣化するのみならず、測定（あるいは推定）結果に含まれる誤差の影響を受けて発散する危険性さえ生ずる。
また、測定点ξ＝ξ₁，ξ₂は、任意に設定される点であるので、例えば、表面温度〔式（１２）］に注目すると、｛ξ₂／（ξ₂−ξ₁）｝や｛ξ₁／（ξ₂−ξ₁）｝の値が整数になるように選定すると式（１２）の留数はｓ＝０だけとなるが、一方、熱流束については、測定点での影響を最も受けることになる。したがって、推定しようとする量に応じて、測定点を選定することによって、ｓ＝０のみの解からθ_w（τ）やΦ_w（τ）を計算することが可能となる。
式（１２），（１３）の第１項の積分はｓ→ｉ∞で必ず発散するが、第２項は、２ξ₁＞ξ₂のとき収束する。したがって、表面温度と表面熱流束はある時刻（最小予測可能時刻）以上に対してのみ逆問題の解が存在することになる。この特徴は、逆問題固有のもので、数学的に証明されている。そこで、解の発散を避けるためのさまざまな工夫によって逆問題の解が得られることになる。本解析で提案された方法もその一つである。
【００１１】
＜表面温度と表面熱流束＞
式（１２），（１３）の解の性質から、無次元時間が大きくなるとｓ＝０のみの解でも十分近似できる可能性があることがわかった。そこで、式（１２），（１３）の被積分関数をまずｓ＝０の近傍で展開し（Ｓ１４０）、その後逆変換を行う（Ｓ１６０）と、表面温度と表面熱流束は陽に与えられ、次式となる。
【数２５】

【００１２】
＜半無限固体に対する解＞
半無限固体に対する一般解は、式（３）においてＡ＝０と置くことによって容易に求まる。
有限固体の場合、二つの未知数であったが、半無限固体の場合、一つの未知数のみとなる点が大きく異なっている。したがって、半無限固体の場合、未知数を決定するために必要な測定点は、１点のみで十分となる。今、ξ＝ξ₁の点での温度変化が式（４）のように与えられているとすると、未定係数Ｂの具体的な関数が定まる。それを用いて、表面温度変化および熱流束を有限固体の場合と同様に求めると最終的に次式のように与えられる。
【数２６】

式（１６），（１７）の複素積分は、ｓの指数が整数でない場合〔式（１７）では、ｋの値が偶数、式（１６）では、ｋの値が奇数〕ｓ＝０に分枝をもつため複雑な計算が必要となるが、ある時刻以上に対して解を陽に得ることができる。なお、ｓの指数が整数の場合の解は、当然ｓ＝０の近傍で展開して得られる解と一致している。
しかし、ここでも有限固体の場合と同様に被積分関数ｅ^{p ξ 1}をｓ＝０の近傍で展開し、その後逆変換を行うと、表面温度と表面熱流束は次式となる。
【数２７】

【００１３】
＜被積分関数の展開＞
式（１２），（１３）中の双曲関数をまずｓ＝０の近傍で展開すると、次式のような級数形で与えられる。
【数２８】

ここで、ｎ＝１，２はそれぞれの点ξ＝ξ₁，ξ₂の値を示す。また、式（２０），（２１）の係数ｃ_i,n，ｄ_i,nの一般項は、次式となる。
【数２９】

ここで、
【数３０】

式（１６）についても同様な展開を行うと、
【数３１】

となる。 式（２０），（２１）に示される級数を用いて、式（１２），（１３）の被積分関数をさらに変形すると次式となる。
【数３２】

Ｃ_-1,21，Ｃ_j,21については、添字の部分に関連した係数のみが入れ替えられる。
【数３３】

Ｄ_-1,21，Ｄ_j,21についても、Ｃ_-1,21，Ｃ_j,21と同様である。
【数３４】

【００１４】
＜温度伝導率等の物性値の測定＞
上記の温度測定を用いて、熱伝導率等の物性値を推定できる。以下にこれについて詳しく説明する。
（温度伝導率の推定）
さて、固体内の温度伝導率が一定の場合、ｘ＝ｘ₁と異なる点（ｘ＝ｘ₂）の温度変化を用いて、上述の半無限平板の場合から予測される表面温度Ｔ_w,2は、ｘ＝ｘ₁から予測された表面温度Ｔ_w,1と一致するはずである。
したがって、温度伝導率が未知であるときには、ｘ＝ｘ₁とｘ＝ｘ₂の異なる２点の温度変化から、式（１８）を用いて、以下のように、推定された表面温度が一致するように、温度伝導率の値を決定すればよい。
【数３５】

ここで、ｔ₁は、予測可能最小時刻、ｔ₂は、温度変化が固体の他の端面に到達するまでの時刻以下のある時刻である。
ところで、温度伝導率ａの値は、まだ未知数であるので、実際には適当な推定値を用いて、ｔ₁，ｔ₂を推定し、上述の式を満足するように繰り返し計算をすることになる。
具体的には、温度伝導率は、上述の式を用いて、以下の手順で推定される。
（１）温度伝導率ａを仮定し、ｔ₁，ｔ₂を推定する。
（２）各点の時間遅れｔ₁ ^*，ｔ₂ ^*を推定する。
（３）順問題の解から測定点ｘ₁，ｘ₂の温度変化を求める。
（４）各点の温度変化から上述の式の係数ｂ_{ｋ , １}，ｂ_{k, ２}を決定する。
（５）係数を計算する。
（６）ｄＦ（ａ）／ｄａ＝０となるａの値を決定する。
（７）得られたａの値が、仮定したａの値と一致するまで繰り返し計算する。
【００１５】
（熱伝導率の推定）
密度が既知である物質の熱伝導率と比熱は、温度伝導率が求まっていると、いずれか一方を求めればよい。しかし、これらの量は、いずれも熱量に関連する物理量であるので、温度の測定だけでは、決定することができない。このため、伝熱量（熱流束）、温度勾配あるいは、温度分布のいずれか１つの量を測定する必要がある。
ところで、上述した逆問題解析結果によれば、半無限固体に対する表面熱流束は、式（１９）に示されているように与えられることが分かっている。熱流束が測定されているとすると、式（１９）を熱伝導率が陽に与えられるように変形（有次元化）すると、次式のようになる。
【数３６】

例えば、表面熱流束が、時間と独立に一定である場合には、次式から簡単に熱伝導率を求めることができる。
【数３７】

一方、熱流束が与えられていない場合には、基準となる有限固体内の温度変化から熱流束を推定し、それを上式に代入すればよい。
【００１６】
【実施例】
図２および図３に本発明の手法の評価の例を示す。
＜具体的な例＞
表１に示される五つの具体的な例について、本発明の境界値問題（初期条件は、いずれもθ＝０）を検討した。
【表１】
表１ 境界条件

【００１７】
＜逆問題の計算結果＞
図２は、具体例（Ｃａｓｅ１）に対して、二つの近似式（１０ａ）と（１０ｂ）からそれぞれ推定された表面温度（逆問題解）と真の表面温度（厳密解）との比較を示す。なお、図２中のＮの値は近似式の次数、また、ξ₁、ξ₂の値は、測定位置を示す。また、各位置での測定精度は、いずれも３けために正規外乱が重畳されたものとなっている。
図２から、いずれの推定値も次数Ｎが大きくなると予測精度が改善されて、厳密解に漸近していることがわかる。また、各近似式間の予測精度を比較すると、時間遅れを考慮した近似式（１０ｃ）から推定された値は、近似式（１０ａ）を用いた場合と比較して、いずれの次数においても非常によい近似を与えていることがわかる。なお、近似式（１０ｂ）を用いて計算した結果は、省略されているが、その結果は、式（１０ａ）の場合とほぼ同様な結果であることを記しておく。
【００１８】
＜推定結果の評価と最適近似法＞
（最小予測可能時刻）
逆問題の解は、τ＝０では、存在しない（発散する）ことが数学的に証明されている。したがって、予測可能な最小時刻τ₁も、解析法を評価するための重要な因子となる。そこで、予測値が、厳密解の９９％（誤差１％以内）に初めてなった時刻を予測可能な最小時刻τ₁とした。
表２は、最小時刻を各近似式およびその次数に対してまとめたものである。なお、表２中には、これから検討する境界条件についても示されている。
【表２】
表２ 最小予測可能時間

表２から、最小時刻τ₁についても時間遅れを考慮した近似式（１０ｃ）を用いた場合が最も小さいことおよび最小時刻の値は、次数による影響を受けにくいことがわかる。
【００１９】
（計測精度）
逆問題の解を次式で与えられる標準偏差で評価することにする。
【数３８】

ここで、τ₂は、測定終了時刻である。なお、Ｃａｓｅ１，２では、近似式の終端効果で精度が劣化するので全測定時間の９０％までの時刻とした。
【００２０】
＜最適な近似法＞
測定点での温度変化の近似式として、三つの異なった方法があるが、表１のＣａｓｅ１（表面温度がステップ上昇）の場合、予想値の精度および最小予測可能時刻τ₁の値から、近似式（１０ｃ）を用いた解析結果が最もよい推定値を与え得ることおよびその近似式の次数も、Ｎ＝５，６程度で十分であることも明らかになった。
（測定結果に及ぼす測定点の影響）
測定位置が、表面に近くなればなるほど予測精度が高くなる。特に、表面に近い測定点については、ξ₁＝０．１以下であることが望まれる。一方、最小予測可能時刻については、ξ₁＝０．２以下であれば、十分であることが確認されている。
【００２１】
＜表面熱流束の推定結果＞
図３は、ここで推奨された近似式（１０ｃ）を用いて、表面熱流束を推定した結果と厳密解との比較を示す。ここでも、表面温度の値は、３けために正規外乱が重畳されたものが使用されている。
図３から、次数Ｎ＝５，７での推定値は、厳密解と非常によく一致していることがわかる。この結果、ここで提案された方法［具体的には式（１５）］は、熱流束の予測にも十分対応できることがわかる。
【００２２】
＜他の具体的問題の計算結果と評価＞
時間の遅れを考慮した近似式（１０ｃ）を用いて、表１に示される残り四つの例に対して同様な検討を行った結果の予測精度をもとめたものを表３に示す。
【表３】
表3 予測精度

表３に示される推定値の標準偏差から、Ｃａｓｅ５以外に対しては、本逆問題解析で得られた結果は、測定精度から判断して十分な精度で推定されていると考えることができるだろう。しかし、Ｃａｓｅ５については、予測精度が他の場合と比較して、かなり劣っている。これは、Ｃａｓｅ５では表面温度の変化がτ＝１で折れ（１階微分係数が不連続）ているために、測定点の温度変化がその点（τ＝１）近くで鈍ってしまうためである。その結果、τ＝１近くに推定値が、大きくずれるという結果となる。なお、測定値の近似式の精度については、いずれの場合に対してもほぼ同等な標準偏差内にあった。したがって、温度変化の１階微分係数が不連続となるような温度変化に対しては、逆問題解の予測精度が著しく劣化することになる。
また、表２，３から、近似式の次数を大きくしても、その精度や最小予測可能時刻の改善があまり期待されなくなり、適当な次数（表１に示される境界条件では、Ｎ＝５，６程度）で十分予測可能となる。
【００２３】
＜初期温度分布の影響＞
固体内に初期定常温度分布が存在する場合、式（１）及びその解を与える式（３）は、温度分布に関連した補正項が追加されることになる。その補正項を追加して得られた結果は、初期温度θ＝０に対する結果と本質的に同じであり、同程度の予測精度が得られる。従って、本手法は、初期定常温度分布が固体内に存在する場合にも、補正項を追加することによってそのまま適用可能となる。
＜他の座標系＞
上述に示された計算は、直角座標系を用いて説明されているが、例えば、円筒座標系や球座標系に対しても、等しく適用可能である。また、他の座標系でも、精度等についても全く同程度となる。
＜２次元逆問題への拡張＞
ここで示された処理を、２次元逆問題に適用することも可能である。２次元問題では、空間座標（例えば、ｘ座標）に付いての温度分布をフーリエ展開するという計算が新たに必要となる。この結果、２次元逆問題の解は、温度分布の近似の程度に影響されるようになるため、その精度は、１次元における解の結果と比較して悪くなる。
【００２４】
＜実際の温度測定のための装置構成例＞
図４は、本発明の温度測定等に用いる機器の構成例を示す図である。
図４において、断熱材２７０に覆われた試験材料２６０内に、温度測定のために、少なくとも２つの熱電対２８０を配置している。熱電対２８０は、基準温度を与えるアイスボックス２４０および増幅器２３０を介して、Ａ／Ｄ変換器２２０でデジタル化されて、コンピュータ２１０に入力され、温度を測定する。試験材料２６０にはヒータ２５０から熱を加えている。コンピュータ・システム２１０はラプラス変換やラプラス逆変換等の本発明の処理も行っている。
【００２５】
＜温度伝導率の推定＞
温度伝導率は、加熱開始時刻（ｔ＝０）と、前に説明した手法で、測定点での温度変化を得ることで、求めることができる。数値実験の結果を図５に示す。図５は、実際の各物質の温度伝導率（ａ_true）と本手法で求めた温度伝導率（定熱流束で熱を与えた場合の温度変化を与える）（ａ_pre）との比較したものである。この図で示されているように、固体の温度伝導率を広い範囲にわたって、２〜３％以内の精度で十分に予測可能であることが分かる。
加熱開始時刻を定めることは、数値実験の場合とは異なり、実際の計測器の測定感度と関連して容易ではないが、以下の手順で決定することができる。
（１）ｘ₁の測定点での時間遅れｔ^*を
【数３９】

から求める。
（２）測定された温度変化とｍｉｎ｛（Ｔ−Ｔ_ｉ）／（Ｔ_L−Ｔ_ｉ）｝＝０．０１との交点がｔ＝ｔ^*となる。
（３）最初に記録された測定温度と時刻（測定上の零点）の点と（２）で定まった点を直線で結び、Ｔ＝Ｔ_ｉとなる点を加熱開始時間とする。
図４と同様の装置構成で、測定点の温度変化を計測することで、上述の手法および加熱開始時間を用いて、色々な物質について温度伝導率を推定することができる。
【００２６】
＜熱伝導率の推定＞
図６に示した装置構成は、熱伝導率を推定するためのものである。図４に示した温度伝導率等を推定するための装置構成と比較すると、既知の基準素材２６５を用いており、その温度変化を計測していることが異なっている。これは、物性が既知である基準素材２６５を用いることで、通過する熱量を求めることができ、この熱量を用いて、前に述べたように熱伝導率を求めることができる。
【００２７】
＜気体や液体の質量濃度と物質拡散流束への応用＞
上述の測定方法を、気体や液体の質量濃度と物質拡散流束に対して応用することが可能である。これを以下に説明する。
熱伝導と拡散に対する支配方程式は、全く同一であるから、上述の熱伝導で提案された推定方法や結果は、そのまま拡散問題に適用でき、拡散係数を推定することができる。
このため、上述の熱伝導の逆問題を解く方法は、温度を質量濃度に、熱流束を物質拡散流束に置き換えて、液体あるいは気体内の物質濃度変化からこれらの表面状態を予測する物質の拡散の問題を解く場合にも適用できる。例えば、アンモニアの蒸気がアンモニア水に吸収される過程における溶液表面でのアンモニア吸収量や質量濃度を直接測定することは困難であるが、これは熱伝導における表面温度直接測定の困難さと同じである。
気体や液体の質量濃度と物質拡散流束に対して適用する場合、上述の説明において、以下に示すように、置き換えることにより、質量濃度と物質拡散流束の場合に対して適用することができる。
温度 −＞ 質量濃度
熱流束 −＞ 物質拡散流束
温度伝導率 −＞ 拡散係数（率）
固体 −＞ 液体あるいは気体
熱伝導 −＞ 拡散
そして、上述の式の記号の意味を以下のように読み替える。
ｆ₁（τ），ｆ₂（τ）：測定点ξ₁，ξ₂での無次元質量濃度
ｋ：質量濃度上昇率（Ｔ／Ｔ₀＝ｋτあるいは２−ｋτ）
Ｌ：代表長さ
ｍｉｎ（θ）：測定器の最小感度あるいは有効数字の最小値
Ｎ：近似式の次数
Ｎ_sf：有効数字のけた
ｐ＝√ｓ
ｑ：物質拡散流束
ｓ：ラプラス演算子（＝ｐ²）
ｔ：時間
ｘ：ｘ座標
Ｔ：質量濃度
Ｔ₀：代表質量濃度
ａ：拡散係数
θ：無次元質量濃度
φ：無次元物質拡散流束
バーθ：無次元質量濃度θのラプラス変換質量濃度（像関数）
バーΦ：無次元物質拡散流束φのラプラス変換物質拡散流束（像関数）
ξ：無次元距離＝ｘ／Ｌ，ξ_l＜ξ₂
σ：標準偏差
τ：無次元時間（フーリエ数＝ａｔ／Ｌ²）
τ₁：最小予測可能時刻
τ_i ^*：無次元時間遅れ
Ｔｗ₁，Ｔｗ₂：測定点ｘ１，ｘ２からの推定表面質量濃度
Ｔ_exact：厳密解質量濃度
このように、上述の式を読み替えることにより、熱伝導と拡散に対する支配方程式は全く同一であるから、上述の熱伝導で提案された推定方法や結果は、そのまま拡散問題に適用でき、拡散係数を推定することができる。
【００２８】
さて、上述の予測法を適用する場合、液体や気体内部の濃度変化を測定する必要がある。液体の場合の濃度変化を測定する機器構成例を、図７を用いて以下に説明する。なお、気体の場合も同様な機器構成で測定することが可能である。
図７（ａ）は、容器３１０中に格納している液体３２０内部の複数測定点における濃度測定をガスクロマトグラフィー等の濃度測定器３４０に接続したセンサ３４２によって測定していることを示している。図７（ｂ）は、ガラス窓３１４を有する容器３１２中に格納している液体３２０にレーザ光を照射して、液体３２０の濃度測定を光干渉計３５０によって測定していることを示している。濃度測定の結果はコンピュータ・システム２１０に送られ、ラプラス変換やラプラス逆変換を行って、必要とする結果を得ることができる。
【００２９】
本発明の手法は、スタンド・アローンのコンピュータ・システムばかりではなく、複数のシステムから構成される例えばクライアント・サーバ・システム等で実行してもよい。
本発明に関するプログラムを格納した記憶媒体から、プログラムをシステムで読み出して実行することにより、本発明の手法を実装することができる。この記録媒体には、フロッピー・ディスク、ＣＤ−ＲＯＭ、磁気テープ、ＲＯＭカセット等がある。このプログラムを電子通信（例えばインターネット）等を利用して伝送することにより、プログラムを伝送されたコンピュータ・システムで実行して、本発明を実施することができる。
【００３０】
【発明の効果】
本発明では、一次元非定常熱伝導の逆問題解をラプラス変換法によって解析し、表面温度と表面熱流束を陽に与えており、逆問題の解は、表面温度を高い予測精度で推定可能であり、測定精度に対して高い安定性を有している。
この発明は、従来の方法にくらべて、装置が簡単で、精度が良く、非定常状態の測定が可能であるなど、多くの特長をもっている。適用分野も、璧温モニタ、熱設計（ＣＡＥ）ツール、温度シミュレータ、熱物性測定器などが考えられ、加熱や冷却が行われる分野で広く利用できる可能性がある。
また、このラプラス変換法を質量濃度の測定に応用することもでき、同様の高い測定精度を得ることができる。これは、濃度モニタ、流体設計（ＣＡＥ）ツール、濃度シミュレータ、物性測定器等が考えられ、液体設計や解析の分野で広く利用できる。
【図面の簡単な説明】
【図１】 本発明の手法を説明するためのフローチャートである。
【図２】 本発明で求めた表面温度の評価結果を示す図である。
【図３】 表面熱流束の評価結果を示す図である。
【図４】 温度変化を測定する機器構成例を示す図である。
【図５】 本発明で求めた温度伝導率の評価を示す図である。
【図６】 熱伝導率を測定する機器構成例を示す図である。
【図７】 液体の拡散係数を測定する機器構成例を示す図である。[0001]
BACKGROUND OF THE INVENTION
  The present invention relates to a technique for predicting physical quantities on these surfaces from changes in internal physical quantities that can be measured when it is difficult to directly measure physical quantities on the surface such as surface temperature and surface heat flux, For example, it is used in the field where heating, cooling or absorption is performed.
[0002]
[Prior art]
  Research on inverse heat conduction problems makes it difficult to measure surface temperature and surface heat flux directly and unsteadily, such as in severe environmental conditions such as spacecraft re-entry and fusion reactor water jetting accidents. In some cases, this is done in connection with the prediction of these surface states from temperature changes in the solid.
  As an analytical method of the one-dimensional inverse problem, for example, there is a method using a Laplace transform method (see Masahiro Shoji, “Research on inverse problems of unsteady heat conduction, theory, 44-381 (1978), 1633-1643”). In this method, there is a problem that the calculation tends to diverge because the solution is given in the differential form of the measured temperature, and the inverse transformation method using the Laplace transform is an extension of a simple system to a two- to three-dimensional problem. However, there is a method extended to the two-dimensional problem (see Imber et al “Prediction of Transient Temperature Distributions with Embedded thermocouples” AIAA J., 12-8 (1978) 1089-1093). There is a problem that prediction accuracy is low due to approximation by a polynomial of time.
  Furthermore, there are conventional methods for measuring thermal properties such as thermal conductivity and thermal conductivity of solids, such as pulse heating, step heating, periodic heating, and continuous heating, but the accuracy of the measurement depends on the boundary conditions. However, there are problems such as difficulty in setting the heating method and the complexity of the measuring apparatus. Further, the arbitrary heating method has a problem that the measurement accuracy is inferior because the accurate measurement time is not known.
[Problems to be solved by the invention]
  The object of the present invention is to solve the inverse problem using the Laplace transform method, without diverging the calculation, to increase the prediction accuracy, and from the change in the internal measurement points, the surface temperature and the surface heat flux (unit area, unit The amount of heat passing per hour) is obtained. It is also an object of the present invention to measure thermophysical properties such as temperature conductivity using the surface temperature obtained in this way.
[0003]
[Means for Solving the Problems]
    In order to achieve the above object, the present invention provides a measuring means for obtaining a temperature change at a measurement point in a solid, a means for determining a coefficient of an approximate expression from the measured temperature change, and the approximation. Means for Laplace transform of the equation, means for obtaining a solution of the Laplace transform equation of the heat conduction equation using the coefficients of the approximate expression subjected to Laplace transform, and means for expanding the obtained solution in the vicinity of s = 0 ,
Find inverse Laplace transform of expanded solutionThe approximate expression of the surface temperature when the solid is a finite field
[Equation 9]

Or an approximate expression of the surface temperature when the solid is a semi-infinite body
[Expression 10]

GetAnd meansMeans for obtaining the temperature conductivity from the fact that the surface temperatures of a plurality of locations inferred from the obtained approximate expression are the same, andWith
  The solidTemperature conductivitymeasurementYouThis is a measuring device.
  Also, change in mass concentration at measurement points in liquid or gasMeasuring means to obtain and measured Change in mass concentrationTo determine the coefficient of the approximate expressionMeans forA means for performing Laplace transform of the approximate expression and a solution of the Laplace transform expression of the diffusion equation are obtained using coefficients of the approximate expression subjected to Laplace transform.meansAnd expand the obtained solution in the vicinity of s = 0meansAnd the inverse Laplace transform of the expanded solutionThe approximate expression of mass concentration when the gas or liquid is a finite body
## EQU11 ##

Or approximate expression of mass concentration when the gas or liquid is a semi-infinite body
[Expression 12]

GetAnd meansMeans for obtaining a diffusion coefficient from the fact that the mass concentrations at a plurality of locations inferred from the obtained approximate expression are the same.WithThe diffusion coefficient of the gas or liquidmeasurementYouMeasuring deviceIt is.
  In addition, the same measurement method,This measurement method program and a recording medium storing the program are also the present invention.
[0004]
DETAILED DESCRIPTION OF THE INVENTION
  When solving the one-dimensional inverse problem using the Laplace transform method, the present invention approximates the measured temperature change with a polynomial of the square root of time, which is a part of the basic solution of the heat conduction equation, and the time delay of the temperature change. By taking into account, the calculation accuracy is not diverged and the prediction accuracy is increased.
  FIG. 1 is a flowchart for explaining the outline of the method of the present invention. In FIG. 1, a Laplace transform is usually performed on a heat conduction equation represented by a differential equation (S110). On the other hand, each coefficient of the approximate expression is determined from the temperature change at an appropriate measurement point in the object (S180). As this approximate expression, an approximate expression having a polynomial of a square root of time considering time delay may be used. Laplace transform is performed on the approximate expression (S190), and a Laplace transform coefficient is determined. A solution of heat conduction that has undergone Laplace transformation is given using the coefficient of the approximate expression (S120). The heat conduction solution given in the Laplace transform form (integration form) is developed near s = 0 (S140). Laplace inverse transformation is executed in the developed form to obtain the surface temperature and the like (S160).
[0005]
  The method for obtaining the surface temperature and the like will be described in detail. First, main symbols used in the following description are shown below.
f₁(Τ), f₂(Τ): Measurement point ξ₁, Ξ₂Dimensionless temperature at
                k: Temperature rise rate (T / T₀= Kτ or 2-kτ)
                L: Representative length
        min (θ): Minimum sensitivity of the measuring instrument or minimum value of significant figures
                N: degree of approximate expression
              N_{science fiction}: Significant digits
                p = √s
                q: Heat flux
                s: Laplace operator (= p²)
                t: time
                x: x coordinate
                T: Temperature
              T₀: Representative temperature
                a: Temperature conductivity
                θ: dimensionless temperature
                λ: thermal conductivity
                φ: Dimensionless heat flux
            Bar θ: Laplace conversion temperature (image function) of dimensionless temperature θ
            Bar Φ: Laplace transform heat flux (image function) of dimensionless heat flux φ
                ξ: dimensionless distance = x / L, ξ_l<Ξ₂
                σ: Standard deviation
                τ: dimensionless time (Fourier number = at / L²)
              τ₁: Minimum predictable time
              τ_i ^*: Dimensionless time delay
      Tw₁, Tw₂: Estimated surface temperature from measurement points x1 and x2
          T_exact: Exact solution temperature
【number13]

[0006]
<Formulation of inverse one-dimensional heat conduction problem>
  A one-dimensional unsteady heat conduction equation with a constant physical property value is expressed in a dimensionless form as follows.
【number14]

  When the initial temperature is uniform (since generality is not lost even when θ = 0), the equation (1) is converted to Laplace (S110).
【number15]

The general solution of equation (2) is easily obtained and becomes
【number16]

  Here, A and B are integral constants.
[0007]
(Solution for finite plate)
  When dealing with the inverse problem of a finite plate, Equation (3) cannot be a closed equation unless at least two temperatures in the object are given. Now, two different points (ξ = ξ_l, Ξ₂, Ξ₂> Ξ₁Suppose the temperature at> 0) is given by:
【number17]

  When unknown constants A and B are determined from equation (4) and substituted into equation (3), bar θ (ξ, S) is
【number18]

Further, the heat flux Φ [= q / (λT₀/ L) =-∂θ / ∂ξ] is
【number19]

It becomes.
  For inverse problems, in general ξ₁≦ ξ ≦ ξ₂The purpose is to predict the temperature and heat flux of points not included in the range. In particular, the main purpose is often to determine the temperature and heat flux at the surface (ξ = 0). Therefore, when ξ = 0 is substituted into the equations (5) and (6), the surface temperature and the surface heat flux are respectively expressed by the following equations.
【number20]

Then, if the equations (7) and (8) are executed according to the inverse Laplace transform definition (given by the equation (9)), the target surface temperature and surface heat flux can be obtained.
【number21]

[0008]
<Approximation formula of measurement temperature>
  To perform the integration of equation (9), the bar θ_w(S) and bar Φ_wThe known function included in (s), the bar f_nIt is necessary to give a specific function form of (s), n = 1, 2 in advance. By the way, bar f_n(S), n = 1 and 2, ξ = ξ₁, Ξ₂Function f indicating temperature change at_nSince (τ) is a Laplace transformed function, ξ = ξ_l, Ξ₂It is determined from the measurement of the temperature change at. Then, it is assumed that the temperature change at each point is approximated by the following three function forms (S180).
【number22]

  Where coefficient a⁽ⁱ⁾ _{k, n}Is a coefficient determined by, for example, least square approximation from the time change of the temperature at each measurement point. N indicates the order of the approximate expression.
  Each of the three functional forms has the following characteristics. Equation (10a) is a random approximation of the temperature change. Equation (10b) is the initial condition f_nIt is an approximate function that satisfies (0) = 0. Equation (10c) is a time delay τ of temperature propagation to the measurement point._n ^*Is an approximate expression that takes into account This time delay τ_n ^*Is erfc (ξ_n/ 2√τ_n ^*) = Min (θ).
  By the way, when the three formulas (10a), (10b), and (10c) are compared, in the formula (10b), a⁽²⁾ _{0, n}≡0 or τ in equation (10c)_n ^*Considering that ≠ 0, the basic function form is almost the same. Therefore, hereinafter, only the expression (10c) is used (hereinafter, superscripts of functions and coefficients are omitted).
  Further, a coefficient b obtained by multiplying the coefficient of the expression (10) by the gamma function_ki= A_kiIf Γ (k / 2 + 1) i = 1, 2 is used, the complexity of notation of the formula can be avoided. And the Laplace transform of Formula (10) becomes following Formula (S190).
【number23]

  By the way, the reason why the half-value polynomial of time is selected as the approximate function is that the solution of the one-dimensional unsteady heat conduction is generally given in a functional form of T = f (x / √ (at)).
[0009]
<How to find approximate formula coefficients>
  Coefficient a in equation (10)_kn(Or b_kn) Is determined by, for example, least square approximation of the temperature change at each measurement point, but depends on the effective number of temperature changes used (temperature measurement accuracy in the experiment). Therefore, the effective number of temperatures is N_{science fiction}= 2 When changing from 2 to 4 digits (The exact solution temperature is expressed as θ (ξ_i, Τ) = Int {θ_exact(Ξ_i, Τ) × 10^Nsf} / 10^NsfAnd a case where there is a disturbance (dispersion σ = 1) due to 3 digits (the temperature is θ (ξ_n, Τ) = θ_exact(Ξ_n, Τ) + 0.005ε (σ = 1), n = 1, 2). By the way, since the accuracy of temperature measurement with a thermocouple is only about 3 digits, it will be sufficient in terms of accuracy to consider the case where there is a disturbance (dispersion σ = 1) due to 3 digits. FIG. 2 shows an approximate curve and temperature change when there is a disturbance when the surface temperature suddenly becomes θ (0, τ) = θ (1, τ) = 1 on a finite flat plate. All approximate expressions are determined by least square approximation of the temperature change.
  From FIG. 2, any approximate expression of the expressions (10a) and (10c) (the approximation of the expression (10b) was omitted because it was almost the same as the expression (10a)), and the order was increased to 5-7. Although a good approximation is given, the accuracy of the approximate expression cannot be improved even if the order N is increased further. Therefore, it is sufficient that the order of the approximate expression is N = 7 at most. In addition, the prediction accuracy of the approximate expression was almost the same when the order N was N = 5 to 7 regardless of which expression was used. Therefore, as for the approximate expression itself, any approximate expression can be used as long as N = about 5 to 7.
  Comparing the accuracy of prediction formulas when the number of significant digits is rounded down and when normal disturbances are superimposed, it is clear that the prediction accuracy of approximate equations always worsens under the same conditions when normal disturbances are superimposed. It was. As a result, the estimation of the inverse problem solution determined using the obtained approximate expression is naturally worse when the normal disturbance is superimposed (the result is omitted). Here, only the inverse problem solution using the approximate expression when the normal disturbance with the worst prediction accuracy is superimposed will be considered.
[0010]
<Inverse Laplace transform and its properties>
  Substituting equation (11) into equation (7) or equation (8) (S120) and further substituting it into equation (9)
【number24]

  When the inverse transformation (residue calculation) of Equations (12) and (13) is performed, the surface temperature and the heat flux can be determined strictly.
  The residues in equations (12) and (13) are easily retained when the exponent of s is an integer (the value of k is even in equation (12), whereas the value of k is odd in equation (13)). Although the sum of numbers can be calculated, if the exponent of s is not an integer, it has a branch at s = 0, so it is not given by an elementary function and becomes a series expansion formula. By the way, the solution obtained by performing the integration can be divided into a set of solutions obtained from s = 0 and others (s ≠ 0). Furthermore, a solution set of s ≠ 0 has a property of decaying with time, and its basic form is exp (−n²τ) / n^kSince (k ≧ 1), the dimensionless time (Fourier number, τ) and the value of n become a set of functions having the property of rapidly converging to zero. Also, the larger the value of k, the faster the convergence to zero. However, when τ becomes small, not only the convergence is deteriorated, but also there is a risk of divergence due to the influence of an error included in the measurement (or estimation) result.
  Also, measurement point ξ = ξ₁, Ξ₂Is an arbitrarily set point. For example, when attention is paid to the surface temperature [formula (12)], {ξ₂/ (Ξ₂−ξ₁)} Or {ξ₁/ (Ξ₂−ξ₁When the value of)} is selected to be an integer, the residue in equation (12) is only s = 0, whereas the heat flux is most affected by the measurement point. Therefore, by selecting the measurement point according to the quantity to be estimated, the solution from s = 0 only_w(Τ) and Φ_w(Τ) can be calculated.
  The integral of the first term in equations (12) and (13) always diverges from s → i∞, but the second term is 2ξ₁> Ξ₂Converge when Therefore, the solution of the inverse problem exists only for a certain time (minimum predictable time) for the surface temperature and the surface heat flux. This feature is inherent to the inverse problem and has been mathematically proven. Therefore, the solution to the inverse problem can be obtained by various devices for avoiding the divergence of the solution. The method proposed in this analysis is one of them.
[0011]
<Surface temperature and surface heat flux>
  From the properties of the solutions of equations (12) and (13), it has been found that if the dimensionless time is increased, a solution with only s = 0 may be sufficiently approximated. Therefore, when the integrands of equations (12) and (13) are first developed in the vicinity of s = 0 (S140) and then subjected to inverse transformation (S160), the surface temperature and the surface heat flux are given positively, The following formula.
【number25]

[0012]
<Solutions for semi-infinite solids>
  The general solution for a semi-infinite solid can be easily obtained by setting A = 0 in equation (3).
  In the case of a finite solid, there are two unknowns, but in the case of a semi-infinite solid, there is a significant difference in that there is only one unknown. Therefore, in the case of a semi-infinite solid, only one measurement point is sufficient to determine the unknown. Now ξ = ξ₁Assuming that the temperature change at the point is given by the equation (4), a specific function of the undetermined coefficient B is determined. Using this, the surface temperature change and the heat flux are obtained in the same manner as in the case of a finite solid, and finally given by the following equation.
【number26]

The complex integrals of equations (16) and (17) are divided into s = 0 when the exponent of s is not an integer [the value of k is an even number in equation (17) and the value of k is an odd number in equation (16)). Since it has branches, a complicated calculation is required, but a solution can be obtained explicitly for a certain time or more. Note that the solution when the exponent of s is an integer is of course consistent with the solution obtained by developing in the vicinity of s = 0.
  However, as in the case of the finite solid, the integrand e^{p ξ 1}Is developed in the vicinity of s = 0, and then the inverse transformation is performed, the surface temperature and the surface heat flux are as follows.
【number27]

[0013]
<Development of integrand>
  When the hyperbolic functions in the equations (12) and (13) are first expanded in the vicinity of s = 0, they are given in the series form as shown in the following equation.
【number28]

Here, n = 1 and 2 are the respective points ξ = ξ₁, Ξ₂Indicates the value of. Further, the coefficient c in the equations (20) and (21)_{i, n}, D_{i, n}The general term is given by
【number29]

here,
【number30]

If the same expansion is performed for equation (16),
【number31]

It becomes. Using the series shown in Equations (20) and (21), the integrand of Equations (12) and (13) is further transformed to the following equation.
【number32]

C_-1,21, C_{j, 21}For, only the coefficients associated with the subscript part are replaced.
【number33]

D_-1,21, D_{j, 21}Also about C_-1,21, C_{j, 21}It is the same.
【number34]

[0014]
<Measurement of physical properties such as temperature conductivity>
  Physical property values such as thermal conductivity can be estimated using the above temperature measurement. This will be described in detail below.
(Estimation of temperature conductivity)
  Now, if the temperature conductivity in the solid is constant, x = x₁(X = x₂), The surface temperature T predicted from the case of the semi-infinite flat plate described above._{w, 2}X = x₁Predicted surface temperature T_{w, 1}Should match.
  Therefore, when the temperature conductivity is unknown, x = x₁And x = x₂From the temperature change at two different points, the value of temperature conductivity may be determined using Equation (18) so that the estimated surface temperatures coincide as follows.
【number35]

Where t₁Is the minimum predictable time, t₂Is a time less than or equal to the time until the temperature change reaches the other end face of the solid.
  By the way, since the value of the temperature conductivity a is still unknown, in practice, an appropriate estimated value is used to₁, T₂Is repeatedly calculated so as to satisfy the above formula.
  Specifically, the temperature conductivity is estimated by the following procedure using the above formula.
(1) Assuming temperature conductivity a, t₁, T₂Is estimated.
(2) Time delay t at each point₁ ^*, T₂ ^*Is estimated.
(3) Measurement point x from the solution of the forward problem₁, X₂Obtain the temperature change.
(4) The coefficient b in the above formula from the temperature change at each point_{k , 1}, B_{k, 2}To decide.
(5) Calculate the coefficient.
(6) The value of a that satisfies dF (a) / da = 0 is determined.
(7) Repeat the calculation until the obtained value of a matches the assumed value of a.
[0015]
(Estimation of thermal conductivity)
  Any one of the thermal conductivity and specific heat of a substance having a known density may be obtained when the temperature conductivity is obtained. However, since these quantities are both physical quantities related to the amount of heat, they cannot be determined only by measuring the temperature. For this reason, it is necessary to measure any one of heat transfer amount (heat flux), temperature gradient, or temperature distribution.
  By the way, according to the inverse problem analysis result described above, it is known that the surface heat flux for the semi-infinite solid is given as shown in the equation (19). Assuming that the heat flux is measured, the following equation is obtained when equation (19) is deformed (made dimensional) so that the thermal conductivity is positively given.
【number36]

  For example, when the surface heat flux is constant independently of time, the thermal conductivity can be easily obtained from the following equation.
【number37]

  On the other hand, when the heat flux is not given, the heat flux is estimated from the temperature change in the finite solid serving as a reference, and is substituted into the above equation.
[0016]
【Example】
  2 and 3 show examples of evaluation of the method of the present invention.
<Specific example>
  Regarding the five specific examples shown in Table 1, the boundary value problem of the present invention (initial conditions are all θ = 0) was examined.
[Table 1]
Table 1 Boundary conditions

[0017]
<Calculation result of inverse problem>
  FIG. 2 shows a comparison between the surface temperature (inverse problem solution) estimated from two approximate equations (10a) and (10b) and the true surface temperature (exact solution) for the specific example (Case 1). . Note that the value of N in FIG. 2 is the order of the approximate expression, and ξ₁, Ξ₂The value of indicates the measurement position. In addition, since the measurement accuracy at each position is 3 digits, normal disturbance is superimposed.
  As can be seen from FIG. 2, the prediction accuracy is improved as the order N increases, and asymptotically approaches the exact solution. In addition, when the prediction accuracy between the approximate expressions is compared, the value estimated from the approximate expression (10c) considering the time delay is much higher in any order than in the case of using the approximate expression (10a). It can be seen that this gives a good approximation. It should be noted that although the result calculated using the approximate expression (10b) is omitted, the result is almost the same as the case of the expression (10a).
[0018]
<Evaluation of estimation results and optimal approximation method>
(Minimum predictable time)
  It has been mathematically proven that the solution to the inverse problem does not exist (diverges) at τ = 0. Therefore, the smallest predictable time τ₁Is also an important factor for evaluating analytical methods. Therefore, the minimum time τ that can predict the time when the predicted value is 99% of the exact solution (within 1% error)₁It was.
  Table 2 summarizes the minimum time for each approximate expression and its order. Table 2 also shows boundary conditions to be studied.
[Table 2]
Table 2 Minimum predictable time

  From Table 2, minimum time τ₁It can also be seen that the approximate expression (10c) considering the time delay is the smallest and the value of the minimum time is less affected by the order.
[0019]
(Measurement accuracy)
  The solution to the inverse problem is evaluated with the standard deviation given by the following equation.
【number38]

Where τ₂Is the measurement end time. In Cases 1 and 2, since the accuracy deteriorates due to the termination effect of the approximate expression, the time is set to 90% of the total measurement time.
[0020]
<Optimal approximation method>
  There are three different methods for approximating the temperature change at the measurement point. In the case 1 of Table 1 (surface temperature is stepped up), the accuracy of the predicted value and the minimum predictable time τ₁From the above values, it was found that the analysis result using the approximate expression (10c) can give the best estimated value, and that the order of the approximate expression is about N = 5,6.
(Influence of measurement points on measurement results)
  The closer the measurement position is to the surface, the higher the prediction accuracy. Especially for measuring points close to the surface, ξ₁= 0.1 or less is desired. On the other hand, for the minimum predictable time,₁= 0.2 or less has been confirmed to be sufficient.
[0021]
<Estimated results of surface heat flux>
  FIG. 3 shows a comparison between the result of estimating the surface heat flux and the exact solution using the approximate expression (10c) recommended here. Here too, the surface temperature value is 3 and a normal disturbance is superimposed on it.
  It can be seen from FIG. 3 that the estimated values at the orders N = 5 and 7 agree very well with the exact solution. As a result, it can be seen that the proposed method [specifically, equation (15)] can sufficiently cope with the prediction of heat flux.
[0022]
<Calculation results and evaluation of other specific problems>
  Table 3 shows the prediction accuracy obtained as a result of performing the same examination on the remaining four examples shown in Table 1 using the approximate expression (10c) considering the time delay.
[Table 3]
Table 3 Prediction accuracy

  From the standard deviation of the estimated values shown in Table 3, for cases other than Case 5, it can be considered that the results obtained by this inverse problem analysis are estimated with sufficient accuracy as judged from the measurement accuracy. Let's go. However, Case 5 is considerably inferior in prediction accuracy compared to other cases. This is because in Case 5, the change in the surface temperature is broken at τ = 1 (the first-order differential coefficient is discontinuous), and the temperature change at the measurement point becomes dull near that point (τ = 1). . As a result, the estimated value is greatly shifted near τ = 1. The accuracy of the approximate expression of the measured value was within a standard deviation that was almost the same in any case. Therefore, the prediction accuracy of the inverse problem solution significantly deteriorates with respect to a temperature change in which the first-order differential coefficient of the temperature change is discontinuous.
  Further, from Tables 2 and 3, even if the order of the approximate expression is increased, improvement in accuracy and minimum predictable time is not expected so much, and an appropriate order (N = 5 under the boundary conditions shown in Table 1). 6), it becomes sufficiently predictable.
[0023]
<Influence of initial temperature distribution>
  When the initial steady temperature distribution exists in the solid, the correction term related to the temperature distribution is added to the expression (1) and the expression (3) that gives the solution. The result obtained by adding the correction term is essentially the same as the result for the initial temperature θ = 0, and the same degree of prediction accuracy is obtained. Therefore, this method can be applied as it is by adding a correction term even when the initial steady-state temperature distribution exists in the solid.
<Other coordinate systems>
  The calculations shown above have been described using a rectangular coordinate system, but are equally applicable to, for example, a cylindrical coordinate system or a spherical coordinate system. In other coordinate systems, the accuracy and the like are almost the same.
<Extension to 2D inverse problem>
  It is also possible to apply the processing shown here to a two-dimensional inverse problem. In the two-dimensional problem, a new calculation for Fourier expansion of the temperature distribution with respect to the space coordinates (for example, the x coordinate) is required. As a result, since the solution of the two-dimensional inverse problem is affected by the degree of approximation of the temperature distribution, the accuracy is worse than the result of the one-dimensional solution.
[0024]
<Example of device configuration for actual temperature measurement>
  FIG. 4 is a diagram illustrating a configuration example of a device used for temperature measurement or the like according to the present invention.
  In FIG. 4, at least two thermocouples 280 are arranged in the test material 260 covered with the heat insulating material 270 for temperature measurement. The thermocouple 280 is digitized by the A / D converter 220 via the ice box 240 and the amplifier 230 that provide a reference temperature, and is input to the computer 210 to measure the temperature. Heat is applied to the test material 260 from the heater 250. The computer system 210 also performs processing of the present invention such as Laplace transform and Laplace inverse transform.
[0025]
<Estimation of temperature conductivity>
  The temperature conductivity can be obtained by obtaining the temperature change at the measurement point by the heating start time (t = 0) and the method described above. The result of the numerical experiment is shown in FIG. FIG. 5 shows the actual temperature conductivity (a_true) And the temperature conductivity obtained by this method (giving temperature change when heat is applied with constant heat flux) (a_pre). As shown in this figure, it can be seen that the temperature conductivity of the solid can be sufficiently predicted over a wide range with an accuracy within 2-3%.
  Unlike the numerical experiment, it is not easy to determine the heating start time in relation to the measurement sensitivity of an actual measuring instrument, but it can be determined by the following procedure.
(1) x₁Time delay t at the measurement point^*The
【number39]

Ask from.
(2) Measured temperature change and min {(T−T_i) / (T_L-T_i)} = 0.01 is the intersection of t = t^*It becomes.
(3) Connect the first recorded measured temperature and time (measurement zero point) and the point determined in (2) with a straight line, and T = T_iThis point is the heating start time.
  By measuring the temperature change at the measurement point with the same apparatus configuration as in FIG. 4, it is possible to estimate the temperature conductivity for various substances using the above-described method and the heating start time.
[0026]
<Estimation of thermal conductivity>
  The apparatus configuration shown in FIG. 6 is for estimating the thermal conductivity. Compared to the apparatus configuration for estimating the temperature conductivity and the like shown in FIG. 4, a known reference material 265 is used, and the temperature change is measured. This is because the amount of heat passing through can be obtained by using the reference material 265 whose physical properties are known, and the heat conductivity can be obtained as described above by using this amount of heat.
[0027]
<Application to mass concentration of gas and liquid and mass diffusion flux>
  The above-described measurement method can be applied to the mass concentration of gas or liquid and the mass diffusion flux. This will be described below.
  Since the governing equations for heat conduction and diffusion are exactly the same, the above estimation method and result proposed for heat conduction can be directly applied to the diffusion problem, and the diffusion coefficient can be estimated.
  For this reason, the method of solving the above inverse problem of heat conduction replaces the temperature with mass concentration and the heat flux with the substance diffusion flux, and predicts the surface state of the substance from the change in the substance concentration in the liquid or gas. It can also be applied to solving diffusion problems. For example, it is difficult to directly measure the amount of ammonia absorbed and the mass concentration at the solution surface in the process where ammonia vapor is absorbed by aqueous ammonia, which is the same as the difficulty of direct surface temperature measurement in heat conduction. .
  When applied to the mass concentration and mass diffusion flux of gas or liquid, in the above description, it can be applied to the case of mass concentration and mass diffusion flux by replacing as shown below. .
Temperature-> Mass concentration
Heat flux-> Mass diffusion flux
Temperature conductivity-> Diffusion coefficient (rate)
Solid-> Liquid or gas
Heat conduction-> diffusion
And the meaning of the symbol of the above-mentioned formula is read as follows.
f₁(Τ), f₂(Τ): Measurement point ξ₁, Ξ₂Dimensionless mass concentration at
                k: Mass concentration increase rate (T / T₀= Kτ or 2-kτ)
                L: Representative length
        min (θ): Minimum sensitivity of the measuring instrument or minimum value of significant figures
                N: degree of approximate expression
              N_{science fiction}: Significant digits
                p = √s
                q: Mass diffusion flux
                s: Laplace operator (= p²)
                t: time
                x: x coordinate
                T: Mass concentration
              T₀: Representative mass concentration
                a: Diffusion coefficient
                θ: dimensionless mass concentration
                φ: Dimensionless mass diffusion flux
            Bar θ: Laplace transform mass concentration (image function) of dimensionless mass concentration θ
            Bar Φ: Laplace transform mass diffusion flux (image function) of dimensionless mass diffusion flux φ
                ξ: dimensionless distance = x / L, ξ_l<Ξ₂
                σ: Standard deviation
                τ: dimensionless time (Fourier number = at / L²)
              τ₁: Minimum predictable time
              τ_i ^*: Dimensionless time delay
      Tw₁, Tw₂: Estimated surface mass concentration from measurement points x1 and x2
          T_exact: Exact mass concentration
  Thus, by replacing the above equation, the governing equations for heat conduction and diffusion are exactly the same, so the estimation methods and results proposed for the above heat conduction can be applied to the diffusion problem as they are, and the diffusion coefficient is Can be estimated.
[0028]
  Now, when applying the above-mentioned prediction method, it is necessary to measure the concentration change inside the liquid or gas. An example of a device configuration for measuring a concentration change in the case of a liquid will be described below with reference to FIG. In addition, in the case of gas, it is possible to measure with the same equipment configuration.
  FIG. 7A shows that the concentration measurement at a plurality of measurement points inside the liquid 320 stored in the container 310 is measured by a sensor 342 connected to a concentration measuring device 340 such as gas chromatography. . FIG. 7B shows that the liquid 320 stored in the container 312 having the glass window 314 is irradiated with laser light, and the concentration measurement of the liquid 320 is measured by the optical interferometer 350. . The result of the concentration measurement is sent to the computer system 210, and Laplace transform or Laplace inverse transform can be performed to obtain the required result.
[0029]
  The method of the present invention may be executed not only by a stand-alone computer system but also by a client server system or the like including a plurality of systems.
  The method of the present invention can be implemented by reading the program from the storage medium storing the program related to the present invention and executing it. Examples of the recording medium include a floppy disk, a CD-ROM, a magnetic tape, and a ROM cassette. By transmitting the program using electronic communication (for example, the Internet) or the like, the program can be executed by the transmitted computer system to implement the present invention.
[0030]
【The invention's effect】
  In the present invention, the inverse problem solution of one-dimensional unsteady heat conduction is analyzed by the Laplace transform method, and the surface temperature and the surface heat flux are given explicitly, and the solution of the inverse problem can estimate the surface temperature with high prediction accuracy. And has high stability with respect to measurement accuracy.
  Compared with the conventional method, the present invention has many features such as a simple apparatus, high accuracy, and measurement in an unsteady state. Possible application fields include a wall temperature monitor, a thermal design (CAE) tool, a temperature simulator, a thermophysical property measuring instrument, and the like, and may be widely used in fields where heating and cooling are performed.
  Moreover, this Laplace conversion method can also be applied to the measurement of mass concentration, and the same high measurement accuracy can be obtained. This may be a concentration monitor, a fluid design (CAE) tool, a concentration simulator, a physical property measuring instrument, etc., and can be widely used in the field of liquid design and analysis.
[Brief description of the drawings]
FIG. 1 is a flowchart for explaining a method of the present invention.
FIG. 2 is a diagram showing the evaluation results of the surface temperature obtained in the present invention.
FIG. 3 is a diagram showing evaluation results of surface heat flux.
FIG. 4 is a diagram showing an example of a device configuration for measuring a temperature change.
FIG. 5 is a diagram showing evaluation of temperature conductivity obtained in the present invention.
FIG. 6 is a diagram showing an example of a device configuration for measuring thermal conductivity.
FIG. 7 is a diagram showing an example of a device configuration for measuring a liquid diffusion coefficient.

Claims (6)

Measuring means for obtaining temperature changes at a plurality of measuring points in the solid;
Means for determining a coefficient of the approximation from a plurality of measured temperature changes;
Means for Laplace transforming the approximate expression;
Means for obtaining a solution of a Laplace transform equation of the heat conduction equation using a coefficient of an approximate equation obtained by Laplace transform;
Means for expanding the obtained solution in the vicinity of s = 0;
Approximate expression of surface temperature when the solid is a finite field by finding the inverse Laplace transform of the developed solution

Or an approximate expression of the surface temperature when the solid is a semi-infinite body

Means for obtaining
A plurality of surface temperatures inferred from the obtained approximate expression, and a means for obtaining a temperature conductivity from a match,
A measuring apparatus for measuring the temperature conductivity of the solid. A measuring means for obtaining a change in mass concentration at a plurality of measuring points in a liquid or gas;
Means for determining a coefficient of the approximate expression from a plurality of measured mass concentration changes;
Means for Laplace transforming the approximate expression;
Means for obtaining a solution of a Laplace transform equation of a diffusion equation using a coefficient of an approximate equation obtained by Laplace transform;
Means for expanding the obtained solution in the vicinity of s = 0;
Calculating the inverse Laplace transform of the developed solution, approximating mass concentration when the gas or liquid is finite

Or approximate expression of mass concentration when the gas or liquid is a semi-infinite body

Means for obtaining
Means for obtaining a diffusion coefficient from the fact that a plurality of mass concentrations estimated from the obtained approximate expression match,
A measuring apparatus for measuring a diffusion coefficient of the gas or liquid. Determining a coefficient of an approximate expression from a plurality of measured temperature changes;
Laplace transforming the approximate expression;
Obtaining a solution of the Laplace transform equation of the heat conduction equation using a coefficient of the approximate equation transformed by Laplace,
Expanding the obtained solution in the vicinity of s = 0;
Approximate expression of surface temperature when the solid is a finite field by finding the inverse Laplace transform of the developed solution

Or an approximate expression of the surface temperature when the solid is a semi-infinite body

And getting the steps
Obtaining a temperature conductivity from a plurality of surface temperatures inferred from the obtained approximate expression, and
A method for measuring the temperature conductivity of the solid. Determining a coefficient of an approximate expression from a plurality of measured mass concentration changes;
Laplace transforming the approximate expression;
Obtaining a solution of a Laplace transform equation of the diffusion equation using a coefficient of an approximate equation obtained by Laplace transform;
Expanding the obtained solution in the vicinity of s = 0;
Calculating the inverse Laplace transform of the developed solution, approximating mass concentration when the gas or liquid is finite

Or approximate expression of mass concentration when the gas or liquid is a semi-infinite body

And getting the steps
Obtaining a diffusion coefficient because a plurality of mass concentrations inferred from the obtained approximate expression match,
A measurement method characterized by measuring a diffusion coefficient of the gas or liquid.   The recording medium which stored the program which performs the measuring method of Claim 3 or 4 with a computer system.   A program that causes a computer system to execute the measurement method according to claim 3 or 4.

Images