ホーム > 特許ランキング > 国立大学法人北海道大学
特許7412786非接触型レオロジー物性計測装置、システム、プログラムおよび方法
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2024-01-04
(45)【発行日】2024-01-15
(54)【発明の名称】非接触型レオロジー物性計測装置、システム、プログラムおよび方法
(51)【国際特許分類】
G01N 11/04 20060101AFI20240105BHJP
G01N 29/02 20060101ALI20240105BHJP
【ＦＩ】
G01N11/04 Z
G01N29/02
【請求項の数】
12
(21)【出願番号】P 2021509496
(86)(22)【出願日】2020-03-25
(86)【国際出願番号】 JP2020013299
(87)【国際公開番号】W WO2020196612
(87)【国際公開日】2020-10-01
【審査請求日】2023-03-07
(31)【優先権主張番号】P 2019056493
(32)【優先日】2019-03-25
(33)【優先権主張国・地域又は機関】JP
(73)【特許権者】
【識別番号】504173471
【氏名又は名称】国立大学法人北海道大学
(74)【代理人】
【識別番号】100110766
【弁理士】
【氏名又は名称】佐川 慎悟
(74)【代理人】
【識別番号】100165515
【弁理士】
【氏名又は名称】太田 清子
(74)【代理人】
【識別番号】100169340
【弁理士】
【氏名又は名称】川野 陽輔
(74)【代理人】
【識別番号】100195682
【弁理士】
【氏名又は名称】江部 陽子
(74)【代理人】
【識別番号】100206623
【弁理士】
【氏名又は名称】大窪 智行
(72)【発明者】
【氏名】田坂 裕司
(72)【発明者】
【氏名】村井 祐一
(72)【発明者】
【氏名】芳田 泰基
【審査官】福田 裕司
(56)【参考文献】
【文献】米国特許出願公開第２００５／０１２６２６８（ＵＳ，Ａ１）
【文献】BIRKHOFER Beat et al.，IN-LINE RHEOMETRY BASED ON ULTRASONIC VELOCITY PROFILES: COMPARISON OF DATA PROCESSING METHODS，Applied Rheology，2012年，Volume.22,Issue.4，PP. 44701-1 - 44701-9，DOI:10.3933/ApplRheol-22-44701
(58)【調査した分野】(Int.Cl.，ＤＢ名)
Ｇ０１Ｎ １１／０４
Ｇ０１Ｎ ２９／０２
(57)【特許請求の範囲】
【請求項1】
管内を流れる流体のレオロジー物性を非接触で計測する非接触型レオロジー物性計測装置であって、超音波流速計測装置から超音波の照射方向に沿った複数の計測点における流速データを取得する流速取得部と、
この流速取得部によって取得された各計測点の流速データを周波数解析することにより前記流体における脈動の周波数から最大の振幅値を示すピーク周波数を検出するピーク周波数検出部と、
レオロジーモデルに基づく管内を流れる流体の運動方程式を時間に対してフーリエ変換して得られる下記式(１）に、前記ピーク周波数と複数の計測点における流速データとを代入して得られる複数の関数式から数値解を算出する数値解算出部と、
この数値解算出部により算出された数値解に基づきレオロジー物性値を決定するレオロジー物性決定部と
を有する、前記非接触型レオロジー物性計測装置。
［数１］
τはレオロジーモデルを用いて関係づけられるせん断応力、
Ａ，Ｂ，Ｃ，・・・はレオロジーモデルにおけるレオロジー物性、
αは圧力勾配、
ρは流体の密度、
ωは流体における脈動の角周波数、
ω０はピーク周波数、
ｒｎは計測点、
ｕは管内の時空間流速分布、
＾は関数のフーリエ変換、
添え字のｍは計測結果、をそれぞれ意味する。
【請求項2】
前記数値解算出部では、下記式（２）に示す費用関数の最小化により数値解を算出する、請求項１に記載の非接触型レオロジー物性計測装置。［数２］
【請求項3】
前記数値解算出部では、代入する流速データとして、連続して隣り合う３点以上の計測点における前記流速データの関数近似値を用いる、請求項１または請求項２に記載の非接触型レオロジー物性計測装置。【請求項4】
前記数値解算出部では、複数種のレオロジーモデルに基づく数値解を算出するとともに、当該数値解同士を比較して最も小さい数値解となる前記レオロジーモデルを前記流体のレオロジーモデルとして決定するレオロジーモデル決定部を有する、請求項１から請求項３のいずれかに記載の非接触型レオロジー物性計測装置。
【請求項5】
前記ピーク周波数検出部では、前記流体の粘性と前記流体における脈動の周波数の比により表される粘性層厚さが、下記式（３）の範囲を満たす周波数の中から最大の振幅値を示すピーク周波数を検出する、請求項１から請求項４のいずれかに記載の非接触型レオロジー物性計測装置。［数３］
（ν/ｋΔω）１/２は粘性層厚さ、
Δｒは計測点同士の距離、
νは流体の動粘性係数（動粘度）、
ｋΔωは検出された周波数、
Ｄは管の内径、をそれぞれ意味する。
【請求項6】
管外から管内に向けて超音波を照射するとともに前記管内から前記管外に向けて反射される超音波を受信して、前記超音波の照射方向に沿った複数の計測点における流速を時系列で計測する超音波流速計測装置と、この超音波流速計測装置から各計測点の流速データを取得して流体のレオロジー物性を決定する請求項１から請求項５のいずれかに記載の非接触型レオロジー物性計測装置と
を有する、非接触型レオロジー物性計測システム。
【請求項7】
管内を流れる流体のレオロジー物性を非接触で計測する非接触型レオロジー物性計測プログラムであって、超音波流速計測装置から超音波の照射方向に沿った複数の計測点における流速データを取得する流速取得部と、
この流速取得部によって取得された各計測点の流速データを周波数解析することにより前記流体における脈動の周波数から最大の振幅値を示すピーク周波数を検出するピーク周波数検出部と、
レオロジーモデルに基づく管内を流れる流体の運動方程式を時間に対してフーリエ変換して得られる下記式(１）に、前記ピーク周波数と複数の計測点における流速データとを代入して得られる複数の関数式から数値解を算出する数値解算出部と、
この数値解算出部により算出された数値解に基づきレオロジー物性値を決定するレオロジー物性決定部と
してコンピュータを機能させる、前記非接触型レオロジー物性計測プログラム。
［数１］
τはレオロジーモデルを用いて関係づけられるせん断応力、
Ａ，Ｂ，Ｃ，・・・はレオロジーモデルにおけるレオロジー物性、
αは圧力勾配、
ρは流体の密度、
ωは流体における脈動の角周波数、
ω０はピーク周波数、
ｒｎは計測点、
ｕは管内の時空間流速分布、
＾は関数のフーリエ変換、
添え字のｍは計測結果、をそれぞれ意味する。
【請求項8】
前記数値解算出部では、下記式（２）に示す費用関数の最小化により数値解を算出する、請求項７に記載の非接触型レオロジー物性計測プログラム。［数２］
【請求項9】
前記数値解算出部では、代入する流速データとして、連続して隣り合う３点以上の計測点における前記流速データの関数近似値を用いる、請求項７または請求項８に記載の非接触型レオロジー物性計測プログラム。【請求項10】
前記数値解算出部では、複数種のレオロジーモデルに基づく数値解を算出するとともに、当該数値解同士を比較して最も小さい数値解となる前記レオロジーモデルを前記流体のレオロジーモデルとして決定するレオロジーモデル決定部を有する、請求項７から請求項９のいずれかに記載の非接触型レオロジー物性計測プログラム。
【請求項11】
前記ピーク周波数検出部では、前記流体の粘性と前記流体における脈動の周波数の比により表される粘性層厚さが、下記式（３）の範囲を満たす周波数の中から最大の振幅値を示すピーク周波数を検出する、請求項７から請求項１０のいずれかに記載の非接触型レオロジー物性計測プログラム。［数３］
（ν/ｋΔω）１/２は粘性層厚さ、
Δｒは計測点同士の距離、
νは流体の動粘性係数（動粘度）、
ｋΔωは検出された周波数、
Ｄは管の内径、をそれぞれ意味する。
【請求項12】
管内を流れる流体のレオロジー物性を非接触で計測する非接触型レオロジー物性計測方法であって、超音波流速計測装置から超音波の照射方向に沿った複数の計測点における流速データを取得する流速取得ステップと、
この流速取得ステップによって取得された各計測点の流速データを周波数解析することにより前記流体における脈動の周波数から最大の振幅値を示すピーク周波数を検出するピーク周波数検出ステップと、
レオロジーモデルに基づく管内を流れる流体の運動方程式を時間に対してフーリエ変換して得られる下記式(１）に、前記ピーク周波数と複数の計測点における流速データとを代入して得られる複数の関数式から数値解を算出する数値解算出ステップと、
この数値解算出ステップにより算出された数値解に基づきレオロジー物性値を決定するレオロジー物性決定ステップと
を有する、前記非接触型レオロジー物性計測方法。
［数１］
τはレオロジーモデルを用いて関係づけられるせん断応力、
Ａ，Ｂ，Ｃ，・・・はレオロジーモデルにおけるレオロジー物性、
αは圧力勾配、
ρは流体の密度、
ωは流体における脈動の角周波数、
ω０はピーク周波数、
ｒｎは計測点、
ｕは管内の時空間流速分布、
＾は関数のフーリエ変換、
添え字のｍは計測結果、をそれぞれ意味する。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、管内を流れる流体のレオロジー物性を計測するための、非接触型レオロジー物性計測装置、システム、プログラムおよび方法に関するものである。
【背景技術】
【0002】
食品・材料・化学工学などの分野では、製品の品質管理、生産プロセスの最適化、工場プラントの保守点検等において、管内の流体のレオロジー物性を計測するニーズがある。
【0003】
例えば、店頭で販売されている牛乳は、牛から搾乳された生乳を殺菌処理するとともに、乳脂肪分などの成分が調整されている。生乳に含まれている乳脂肪分等の割合は、牛に与えた飼料、牛の体調、天候、および個々の牛によって異なるため一定に保つことは難しい。そこで、これまでは、乳脂肪分を計測等して常に管理するのではなく、必ず所定割合以上の乳脂肪分が含まれるように過剰に乳脂肪分を添加して調整する方法が採られている。
【0004】
しかし、過剰な乳脂肪分を添加したことによる費用コストが高くなることが問題となっている。そこで、管内を流れる牛乳のレオロジー物性を計測することで乳脂肪分の割合を見積もり、時々刻々と変化する成分をリアルタイムで調整したいというニーズが生まれている。
【0005】
そこで、これまでに管内の流れの流速分布を計測することのできる超音波流速計測装置と差圧計とにより、管内の流体のレオロジー物性を計測する手法が提案されている（非特許文献１）。この手法は、管内の流れを一方向定常流れとして仮定するとともにレオロジー物性を記述するモデル式（レオロジーモデル）に基づく流れの式を連立することにより得られる流速分布と、前記超音波流速計測装置により得られる流速分布とを比較することで、モデル式に含まれるレオロジー物性の定数を決定するものである。
【先行技術文献】
【非特許文献】
【0006】
【文献】Ouriev, B., and E. J. Windhab, "Rheological study of concentrated suspensions in pressure-driven shear flow using novel in-line ultrasound Doppler method", Exp. Fluids 32, 204-211 (2002).
【発明の概要】
【発明が解決しようとする課題】
【0007】
しかしながら、非特許文献１に記載された手法においては、管内流れが定常な流れでなければ計測原理上、精度を保証することができないという問題がある。つまり、前記手法は、管内の流れを一方向定常流れと仮定しているが、一般的な生産プロセスにおける管内の流れは、その流れを生み出すためのポンプや配管形状等によって発生する渦流等により常に脈動（変動）しているため、前記脈動による誤差が生じてしまう。
【0008】
また、レオロジーモデルに基づく流れの式には、流体を駆動する圧力勾配が含まれる。そのため、レオロジー物性を算出するには管内の圧力値を代入する必要がある。しかし、管内の圧力を計測するには、管内に圧力センサーを設けたり、管壁に孔を空けて圧力計に接続したりして、圧力センサーや圧力計を流体に接触させる必要がある。このため、例えば、衛生面が重要視される食品等では、常に圧力センサー等を清潔な状態に保たなければならない。流体によっては、圧力センサー等が直接触れること自体が不可のものがあり、適用対象が限られるという問題もある。また、圧力センサー等自体の計測精度がレオロジー物性の算出結果に与える影響も大きく、圧力センサー等の選定や当該圧力センサー等の設置箇所の選択も制限される。
【0009】
本発明は、以上のような問題点を解決するためになされたものであって、管内の流れにおける脈動を利用して非接触でレオロジー物性を計測することのできる、非接触型レオロジー物性計測装置、システム、プログラムおよび方法を提供することを目的としている。
【課題を解決するための手段】
【0010】
本発明に係る非接触型レオロジー物性計測装置および非接触型レオロジー物性計測プログラムは、管内の流れにおける脈動を利用することで、圧力計などによって管内の圧力を計測することなく、非接触でレオロジー物性を計測するという課題を解決するために、管内を流れる流体のレオロジー物性を非接触で計測する非接触型レオロジー物性計測装置であって、超音波流速計測装置から超音波の照射方向に沿った複数の計測点における流速データを取得する流速取得部と、この流速取得部によって取得された各計測点の流速データを周波数解析することにより前記流体における脈動の周波数から最大の振幅値を示すピーク周波数を検出するピーク周波数検出部と、レオロジーモデルに基づく管内を流れる流体の運動方程式を時間に対してフーリエ変換して得られる下記式(１）に、前記ピーク周波数と複数の計測点における流速データとを代入して得られる複数の関数式から数値解を算出する数値解算出部と、この数値解算出部により算出された数値解に基づきレオロジー物性値を決定するレオロジー物性決定部とを有する。
［数１］
ここで、
、
、
、
、τはレオロジーモデルを用いて関係づけられるせん断応力、Ａ，Ｂ，Ｃ，・・・はレオロジーモデルにおけるレオロジー物性、αは圧力勾配、ρは流体の密度、ωは流体における脈動の角周波数、ω０はピーク周波数、ｒｎは計測点、ｕは管内の時空間流速分布、＾は関数のフーリエ変換、添え字のｍは計測結果、をそれぞれ意味する。
［数１］
【0011】
また、本発明の一態様として、数値解算出部による数値解の算出処理に係る負荷を軽減してリアルタイムなレオロジー物性の計測を可能にするという課題を解決するために、前記数値解算出部では、下記式（２）に示す費用関数の最小化により数値解を算出するようにしてもよい。
［数２］
［数２］
【0012】
さらに、本発明の一態様として、超音波流速計測装置により得られる流速データの計測誤差により生じるレオロジー物性値の算出誤差を抑制するという課題を解決するために、前記数値解算出部では、代入する流速データとして、連続して隣り合う３点以上の計測点における前記流速データの関数近似値を用いるようにしてもよい。
【0013】
また、本発明の一態様として、管内を流れる流体のレオロジーモデルが不明の状態で計測をすることができるという課題を解決するために、前記数値解算出部では、複数種のレオロジーモデルに基づく数値解を算出するとともに、当該数値解同士を比較して最も小さい数値解となる前記レオロジーモデルを前記流体のレオロジーモデルとして決定するレオロジーモデル決定部を有するようにしてもよい。
【0014】
さらに、本発明の一態様として、計測精度を担保するために十分に厚い粘性層厚さとなるピーク周波数を検出するという課題を解決するために、前記ピーク周波数検出部では、前記流体の粘性と前記流体における脈動の周波数の比により表される粘性層厚さが、下記式（３）の範囲を満たす周波数の中から最大の振幅値を示すピーク周波数を検出するようにしてもよい。
［数３］
ここで、（ν/ｋΔω）１/２は粘性層厚さ、Δｒは計測点同士の距離、νは流体の粘性係数（動粘度）、ｋΔωは検出された周波数、Ｄは管の内径、をそれぞれ意味する。
［数３］
【0015】
本発明に係る非接触型レオロジー物性計測システムは、管内の流れにおける脈動を利用することで、圧力計などによって管内の圧力を計測することなく、非接触でレオロジー物性を計測するという課題を解決するために、管外から管内に向けて超音波を照射するとともに前記管内から前記管外に向けて反射される超音波を受信して、前記超音波の照射方向に沿った複数の計測点における流速を時系列で計測する超音波流速計測装置と、この超音波流速計測装置から各計測点の流速データを取得して流体のレオロジー物性を決定する前記非接触型レオロジー物性計測装置とを有する。
【0016】
本発明に係る非接触型レオロジー物性計測方法は、管内の流れにおける脈動を利用することで、圧力計などによって管内の圧力を計測することなく、非接触でレオロジー物性を計測するという課題を解決するために、管内を流れる流体のレオロジー物性を非接触で計測する非接触型レオロジー物性計測方法であって、超音波流速計測装置から超音波の照射方向に沿った複数の計測点における流速データを取得する流速取得ステップと、この流速取得ステップによって取得された各計測点の流速データを周波数解析することにより前記流体における脈動の周波数から最大の振幅値を示すピーク周波数を検出するピーク周波数検出ステップと、レオロジーモデルに基づく管内を流れる流体の運動方程式を時間に対してフーリエ変換して得られる上記式(１）に、前記ピーク周波数と複数の計測点における流速データとを代入して得られる複数の関数式から数値解を算出する数値解算出ステップと、この数値解算出ステップにより算出された数値解に基づきレオロジー物性値を決定するレオロジー物性決定ステップとを有する。
【発明の効果】
【0017】
本発明によれば、管内の流れにおける脈動を利用して非接触でレオロジー物性を計測することができる。
【図面の簡単な説明】
【0018】
【図1】本発明に係る非接触型レオロジー物性計測システムの一実施形態を示すブロック図である。
【図2】本実施形態の非接触型レオロジー物性計測システムによるレオロジー物性計測処理の流れを示すフロー図である。
【図3】実施例１において計算された脈動を伴う管内流れの流速データの計算結果を示すコンター図である。
【図4】本実施例１において式（１８）により定義されたノイズレベルa=5.0のノイズ速度の確率密度分布を示すグラフである。
【図5】本実施例１において式（１８）により計算されたノイズを加えた流速データを示すコンター図である。
【図6】本実施例１において式（７）により算出された実数値Reおよび虚数値Imを示すグラフである。
【図7】本実施例１においてランダムサーチ法による費用関数の計算結果を示す３次元グラフ図である。
【図8】本実施例１において計測誤差（ノイズレベル）に対する算出された粘性係数μおよび圧力勾配αの計測精度を示すグラフである。
【図9】本実施例１において周期の違いにおける算出された粘性係数μおよび圧力勾配αの計測精度を示すグラフである。
【図10】本実施例１において動粘性係数νに対する算出された粘性係数μおよび圧力勾配αの計測精度を示すグラフである。
【図11】本実施例１において周波数foに対する算出された粘性係数μおよび圧力勾配αの計測精度を示すグラフである。
【図12】本実施例１において振幅U1に対する算出された粘性係数μおよび圧力勾配αの計測精度を示すグラフである。
【図13】実施例２において用いられた実験装置を示す模式図である。
【図14】本実施例２の実験装置における超音波トランスデューサの設置状態を示す模式図である。
【図15】本実施例２の脈動を伴う管内流れにおける速度分布であって試験流体をニュートン流体とした計測結果を示すコンター図である。
【図16】本実施例２において試験流体をニュートン流体とした場合のランダムサーチ法による費用関数の計算結果を示す図である。
【図17】本実施例２において試験流体をニュートン流体とした場合の粘度および圧力振幅の計測結果を示す図である。
【図18】本実施例２の脈動を伴う管内流れにおける速度分布であって試験流体を非ニュートン流体とした計測結果のコンター図および式（１０）により算出された実数値Reおよび虚数値Imを示すグラフである。
【図19】本実施例２において試験流体を非ニュートン流体とした場合のランダムサーチ法による粘度および圧力（実部・虚部）の二乗和平方根の費用関数の計算結果を示す図である。
【図20】本実施例２において試験流体を非ニュートン流体とした場合のランダムサーチ法による費用関数の計算結果を示す３次元グラフ図である。
【図21】本実施例２において試験流体を非ニュートン流体とした場合に計測されたせん断応力とひずみ速度を示すグラフである。
【図22】本実施例２の試験流体を非ニュートン流体とした実験において計測された広範囲のひずみ速度に対するせん断応力を示すグラフである。
【発明を実施するための形態】
【0019】
以下、本発明に係る非接触型レオロジー物性計測装置、システム、プログラムおよび方法の一実施形態について図面を用いて説明する。
【0020】
本実施形態の非接触型レオロジー物性計測システム１は、図１に示すように、複数の計測点における流速を時系列で計測する超音波流速計測装置２と、この超音波流速計測装置２から各計測点の流速データを取得して流体のレオロジー物性を決定する非接触型レオロジー物性計測装置３とを有する。以下、各構成について詳細に説明する。
【0021】
超音波流速計測装置２は、管外から管内に向けて超音波を照射するとともに前記管内から前記管外に向けて反射される超音波を受信して、前記超音波の照射方向に沿った複数の計測点における流速を時系列で計測するものである。つまり、超音波流速計測装置２は、管内の時空間流速分布u(r,t)を非接触かつ非侵襲（流れを乱さず）に計測することができる。例えば、特開２００３－３４４１３１号公報に開示されている技術を用いることができる。
【0022】
非接触型レオロジー物性計測装置３は、超音波流速計測装置２から各計測点の流速データを取得して流体のレオロジー物性を決定するものである。本実施形態における非接触型レオロジー物性計測装置３は、コンピュータによって構成されており、図１に示すように、主として、各種の表示画面を表示するとともに各種のデータを入力する表示入力手段４と、各種のデータを記憶するとともに演算処理手段６が演算処理を行う際のワーキングエリアとして機能する記憶手段５と、記憶手段５にインストールされた非接触型レオロジー物性計測プログラム３ａを実行することにより、各種の演算処理を実行し後述する各構成部として機能する演算処理手段６とを有する。
【0023】
表示入力手段４は、入力機能と表示機能とを有するユーザインターフェースであり、本実施形態では、タッチパネル機能を備えたディスプレイによって構成されている。なお、表示入力手段４の構成は、タッチパネル式のディスプレイによるものに限定されるものではなく、表示機能のみを備えた表示手段、およびキーボードやマウスなどの入力機能のみを備えた入力手段をそれぞれ別個に有していてもよい。
【0024】
記憶手段５は、ハードディスク、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、フラッシュメモリ等で構成されており、非接触型レオロジー物性計測プログラム３１ａを記憶するプログラム記憶部５１と、レオロジーモデルに基づく管内を流れる流体の運動方程式を時間に対してフーリエ変換して得られる式を記憶するレオロジーモデル式記憶部５２と、各種のしきい値を記憶するしきい値記憶部５３とを有する。
【0025】
プログラム記憶部５１には、本実施形態の非接触型レオロジー物性計測装置３を制御するための非接触型レオロジー物性計測プログラム３ａがインストールされている。そして、演算処理手段６が、当該非接触型レオロジー物性計測プログラム３ａを読み出して実行することにより、コンピュータを非接触型レオロジー物性計測装置３における各構成部として機能させるようになっている。
【0026】
なお、非接触型レオロジー物性計測プログラム３ａの利用形態は、上記構成に限られるものではない。例えば、ＣＤ－ＲＯＭやＵＳＢメモリ等のように、コンピュータで読み取り可能な非一時的な記録媒体に非接触型レオロジー物性計測プログラム３ａを記憶させておき、当該記録媒体から直接読み出して実行してもよい。また、外部サーバ等からクラウドコンピューティング方式やＡＳＰ（Application Service Provider）方式等で利用してもよい。
【0027】
レオロジーモデル式記憶部５２は、レオロジーモデルに基づく管内を流れる流体の運動方程式を時間に対してフーリエ変換して得られる式を記憶するものである。本実施形態では、下記式（１）に基づきレオロジーモデルに導き出された式が記憶されている。
［数１］
ここで、
、
、
、
、
τはレオロジーモデルを用いて関係づけられるせん断応力、
Ａ，Ｂ，Ｃ，・・・はレオロジーモデルにおけるレオロジー物性、
αは圧力勾配、
ρは流体の密度、
ωは流体における脈動の角周波数、
ω０はピーク周波数、
ｒｎは計測点、
ｕは管内の時空間流速分布、
＾は関数のフーリエ変換、
添え字のｍは計測結果、をそれぞれ意味する。なお、レオロジー物性のＡ，Ｂ，Ｃ，・・・はレオロジーモデルに応じてその数が変わるものである。例えば、測定対象に適したレオロジーモデルがニュートン流体である場合のレオロジー物性の数は１であり、非ニュートン流体の場合のレオロジー物性の数は２以上である。
［数１］
τはレオロジーモデルを用いて関係づけられるせん断応力、
Ａ，Ｂ，Ｃ，・・・はレオロジーモデルにおけるレオロジー物性、
αは圧力勾配、
ρは流体の密度、
ωは流体における脈動の角周波数、
ω０はピーク周波数、
ｒｎは計測点、
ｕは管内の時空間流速分布、
＾は関数のフーリエ変換、
添え字のｍは計測結果、をそれぞれ意味する。なお、レオロジー物性のＡ，Ｂ，Ｃ，・・・はレオロジーモデルに応じてその数が変わるものである。例えば、測定対象に適したレオロジーモデルがニュートン流体である場合のレオロジー物性の数は１であり、非ニュートン流体の場合のレオロジー物性の数は２以上である。
【0028】
式（１）において圧力勾配αを表す項は時間によってフーリエ変換することにより流れにおける角周波数ωの関数となっている。このため、前記角周波数ωを特定すれば圧力勾配αは未知の定数として扱うことができる。ここで、前記角周波数ωは管内の時空間流速分布u(r,t)によって得ることができる。つまり、式（１）では、圧力勾配αの値を得るために、圧力センサーや圧力計などを用いて計測する必要がない。また、上述のとおり超音波流速計測装置２は、非接触かつ非侵襲で管内の時空間流速分布u(r,t)を計測することができる。つまり、本実施形態の非接触型レオロジー物性計測システム１は、計測原理的に非接触により管内を流れる流体のレオロジー物性を計測することができるシステムである。
【0029】
また、式（１）では、流体における脈動の角周波数ωが含まれる。一般的に、生産プロセスにおける管内の流れは、その流れを生み出すためのポンプや配管形状等によって発生する渦流等により常に脈動（変動）していることから、そのような管内を流れる流体のレオロジー物性の計測に使用することができる。よって、従来の定常流を仮定できる管内流れしか適用できないシステムに比べて、適用対象となる範囲は広い。
【0030】
なお、レオロジーモデルは、流れのせん断応力と流れの速度勾配とが比例する性質（線形）があるニュートン流体モデルと、前記せん断応力と前記速度勾配とが比例しない性質（非線形）がある非ニュートン流体モデルとに大別される。また、非ニュートン流体モデルには、塑性（ビンガム）流体モデル、準（擬）粘性流体モデル、準（擬）塑性流体モデル、ダイラタント流体モデル等が例示される。
【0031】
しきい値記憶部５３は、各種のしきい値を記憶するものである。本実施形態では、粘性層厚さを基準にピーク周波数ω０を検出するために用いられる下記式（３）およびそれに伴う各値が記憶されている。
［数３］
ここで、
（ν/ｋΔω）１/２は粘性層厚さ、
Δｒは計測点同士の距離、
νは流体の動粘性係数、
ｋΔωは検出された周波数、
Ｄは管の内径、をそれぞれ意味する。
［数３］
（ν/ｋΔω）１/２は粘性層厚さ、
Δｒは計測点同士の距離、
νは流体の動粘性係数、
ｋΔωは検出された周波数、
Ｄは管の内径、をそれぞれ意味する。
【0032】
つぎに、演算処理手段６について説明する。演算処理手段６は、ＣＰＵ（Central Processing Unit）等によって構成されており、記憶手段５にインストールされた非接触型レオロジー物性計測プログラム３ａを実行することにより、図１に示すように、超音波流速計測装置２から流速データを取得する流速取得部６１と、流体における脈動の周波数から最大の振幅値を示すピーク周波数ω０を検出するピーク周波数検出部６２と、ピーク周波数ω０と複数の計測点における流速データとを代入して得られる複数の関数式から数値解を算出する数値解算出部６３と、レオロジーモデルが不明の場合や変化する場合にレオロジーモデルを決定するレオロジーモデル決定部６４と、前記数値解算出部６３により算出された数値解に基づきレオロジー物性値を決定するレオロジー物性決定部６５として機能するようになっている。
【0033】
流速取得部６１は、超音波流速計測装置２から超音波の照射方向に沿った複数の計測点における流速データを取得するものである。流速データは、厳密には超音波流速計測装置２の測定限界に基づく時間間隔毎に計測される離散的なデータであるが、数ミリ秒～数十ミリ秒程度の流体における脈動の周波数よりも充分に短い時間間隔で計測できるため時間的に連続したデータとして扱うことができる。また、式（３）の計算に用いるため、本実施形態における流速取得部６１は、前記超音波流速計測装置２から測定限界に基づき決定される計測点同士の距離Δｒを取得し、しきい値記憶部５３に記憶させるようになっている。
【0034】
ピーク周波数検出部６２は、各計測点の流速データを周波数解析することにより前記流体における脈動の周波数から最大の振幅値を示すピーク周波数ω０を検出するものである。周波数解析は、スペクトル解析や波形分析ともいう。例えば、フーリエ変換やウェーブレット変換等が例示される。本実施形態では、フーリエ変換を用いてピーク周波数ω０を検出している。なお、ピーク周波数ω０は、優位周波数や卓越周波数と呼ばれることもある。
【0035】
また、本実施形態におけるピーク周波数検出部６２は、粘性層厚さを基準にピーク周波数を検出する機能を有する。具体的には、上記の式（３）を満たす周波数の中から最大の振幅値を示す周波数をピーク周波数ω０として検出する。ここで式（３）の最低値である３Δｒは、数値解算出部６３において数値解を算出する上で、選択されるレオロジーモデルによって定められる未知の定数の数が少なくとも３つ以上あり、それらの未知の定数を算出するために必要となる関数式の数が３本以上であること、また超音波流速計測装置２の計測誤差によるレオロジー物性の計測精度の低下を抑制するために流速データやレオロジー物性値の算出結果を空間差分や多項式近似するために必要なデータ点数（計測点の数）が３点以上であることに基づく。よって、粘性層厚さがこれよりも薄いと精度の高い数値解が得られない。また、最大値のＤ／２は管の半径を意味しており、軸対称となる管内流れにおいて粘性層厚さがこれ以上厚くなることがないことに基づく。
【0036】
数値解算出部６３は、上記の式（１）にピーク周波数ω０と複数の計測点における流速データとを代入して得られる複数の関数式から数値解を算出するものである。本実施形態における数値解算出部６３は、下記式（２）に示す費用関数の最小化により数値解を算出する。本実施形態では、費用関数を最小化するプログラムを作成し、その作成プログラムを実行することで数値解を算出している。費用関数を最小化するアルゴリズムは、特に限定されるものではないが、例えばランダムサーチ法や勾配降下法などが例示される。
［数２］
［数２］
【0037】
また、本実施形態における数値解算出部６３は、代入する流速データとして、連続して隣り合う３点以上の計測点における前記流速データの関数近似値を用いる。これは、超音波流速計測装置２により得られる流速データの計測誤差により生じるレオロジー物性値の算出誤差を抑制するためである。関数近似値の算出方法は、特に限定されるものではないが、ベッセル関数やチェビシェフ級数などが例示される。
【0038】
また、本実施形態における数値解算出部６３は、管内を流れる流体のレオロジーモデルが不明な場合や、温度または各種成分割合の変化によってレオロジーモデルが時々刻々と変化する場合には、複数種のレオロジーモデルに基づく数値解を算出する。具体的には、数値解算出部６３が、レオロジーモデル式記憶部５２から複数種のレオロジーモデルに基づく管内を流れる流体の運動方程式を時間に対してフーリエ変換して得られる式を用いて数値解を算出する。
【0039】
レオロジーモデル決定部６４は、レオロジーモデルが不明な場合などにおいて数値解算出部６３が複数種のレオロジーモデルに基づく数値解を算出した上で、レオロジーモデルを決定するものである。本実施形態におけるレオロジーモデル決定部６４は、数値解算出部６３が算出した複数の数値解同士を比較して最も小さい数値解となる前記レオロジーモデルを前記流体のレオロジーモデルとして決定する。
【0040】
なお、測定対象に適したレオロジーモデルが不明な場合のレオロジーモデルを決定する方法は、本実施形態のように数値解算出部６３により算出した複数の数値解同士を比較する方法に限定されるものではなく、例えば、数値解算出部６３により算出された数値解が所定のしきい値以下か否かを判別し、前記しきい値以下の場合には当該数値解を算出したレオロジーモデルを計測対象の流体のレオロジーモデルと決定し、前記しきい値より大きい場合には、レオロジーモデル式記憶部５２に記憶された他のレオロジーモデルに基づく式に変更した上で数値解を算出し、当該数値解が前記しきい値以下となってレオロジーモデルが決定されるまで当該判別処理を繰り返すようにしてもよい。
【0041】
レオロジー物性決定部６５は、数値解算出部６３により算出された数値解に基づきレオロジー物性値を決定するものである。本実施形態におけるレオロジー物性決定部６５は、算出された数値解から計測対象としている管内の流体の粘性係数μや動粘性係数νなどのレオロジー物性値を決定するようになっている。また、決定された動粘性係数νは、ピーク周波数検出部６２で算出される式（３）を算出する場合の動粘性係数νとして用いるため、しきい値記憶部５３に記憶される。ここで動粘性係数νが時々刻々と変化する場合は、前記しきい値記憶部５３に時系列で並べて記憶させるか、算出されるごとに書き換え（上書き）するようにしてもよい。
【0042】
なお、前記レオロジー物性値と、流体に含まれる特定成分の割合との間に相関などの関係性がある場合には、前記レオロジー物性値とともに、または前記レオロジー物性値に変えて当該成分割合を決定するようにしてもよい。
【0043】
つぎに、本実施形態の非接触型レオロジー物性計測装置３、システム１およびプログラム３ａにおける各構成の作用について、非接触型レオロジー物性計測方法とともに説明する。
【0044】
超音波流速計測装置２が、管内を流れる流体に対して管外から超音波の送受信を行い、受信した超音波を解析することにより、前記超音波の照射方向に沿った複数の計測点における流速を時系列に計測する。
【0045】
図２に示すように、非接触型レオロジー物性計測装置３では、流速取得部６１が、超音波流速計測装置２によって計測された流速データ、具体的には時空間流速分布u(ri,tj)、（i=1,2,・・・,N、j=1,2,・・・,M）を取得する（Ｓ１：流速取得ステップ）。本実施形態における流速取得部６１は、前記超音波流速計測装置２から流速データとともに、当該流速データを算出する時に定められた計測点同士の距離Δｒを取得し、しきい値記憶部５３に記憶させる。
【0046】
次に、ピーク周波数検出部６２が、流速取得部６１によって取得された各計測点の流速データを周波数解析し、流体における脈動の周波数から最大の振幅値を示すピーク周波数ω０を検出する（Ｓ２：ピーク周波数検出ステップ）。本実施形態では、下記時空間流速分布u(ri,tj)のフーリエ変換を行った式（４）と、各周波数における振幅値を示す下記の式（５）を用いて周波数解析（スペクトル解析）を行いピーク周波数ω０を検出する。
［数４］
［数５］
［数４］
【0047】
また、ピーク周波数検出部６２が、算出されたピーク周波数ω０が式（３）の範囲内の条件を満たすか否かを判別する（Ｓ３）。当該式（３）における計測点同士の距離Δｒ、流体の動粘性係数νおよび管の内径Ｄはしきい値記憶部６３から取得する。もし、ステップＳ２で算出されたピーク周波数ω０が式（３）の範囲内にない場合には（Ｓ３：ＮＯ）、ステップＳ２に戻って次に大きいピーク周波数ω０を検出し、式（３）の要件が満たすまでこの処理を繰り返す。そして、算出されたピーク周波数ω０が式（３）の範囲内である場合には（Ｓ３：ＹＥＳ）、その周波数をピーク周波数ω０に決定する。本実施形態におけるピーク周波数検出部６２では、このようにして式（３）の要件を満たす周波数の中から最大の振幅値をピーク周波数ω０として検出する。
【0048】
次に、数値解算出部６３が、式(１）に代入する流速データについて、超音波流速計測装置２における計測誤差の影響を抑制するために、流速データの関数近似値を算出する（Ｓ４）。本実施形態では、下記の式（６）を用いて５次関数近似値を計算しその値を用いる。
［数６］
［数６］
【0049】
そして、式(１）に基づきピーク周波数ω０と複数の計測点における流速データとを代入して得られる複数の関数式から数値解を算出する（Ｓ５：数値解算出ステップ）。本実施形態では、数値解を算出するために式（２）に示す費用関数を用いる。レオロジーモデルがニュートン流体の場合、式（２）は、下記の式（７）および式（８）に示すように書き換えることができる。
［数７］
［数８］
ここで、α（ｔｉ）は圧力勾配の時間変動を意味する。
［数７］
【0050】
また、レオロジーモデルが非ニュートン流体であり、レオロジーモデルを用いて関係づけられるせん断応力τを下記の式（９）で表した場合、式（２）は、下記の式（１０）ないし式（１２）に示すように書き換えることができる。
［数９］
［数１０］
［数１１］
［数１２］
ここで、Π1, Π2, Π3, ・・・は、非ニュートン流体におけるレオロジー物性を意味する。なお、式（１１）は上記式（８）と同じ式であり、α（ｔｉ）は圧力勾配の時間変動を意味する。
［数９］
【0051】
そして、この費用関数を、ランダムサーチ法や勾配降下法などを用いて最小化することで数値解を算出する。このとき、各式に代入されるのは前記ピーク周波数検出部６２により検出されたピーク周波数ω０や流速取得部６１により取得された流速データのみであり、管内の圧力データは不要である。また、費用関数を用いることで、数値解を算出するための計算負荷が抑制される。費用関数の最小化するためのアルゴリズムの種類や演算処理手段の処理能力に依存するが、１回の算出に対しおよそ数秒から数十秒間隔で数値解の算出が可能である。
【0052】
本実施形態では、ステップＳ５において、レオロジーモデルが不明な場合や温度や成分割合などに応じてせん断応力τが変化しレオロジーモデルの変更が必要な場合に対応するため、複数のレオロジーモデルを用いて数値解を算出する。具体的には、数値解算出部６３がレオロジーモデル式記憶部５２に記憶された複数種のレオロジーモデルに基づく管内を流れる流体の運動方程式を時間に対してフーリエ変換して得られる複数の式を読み出し、それぞれの式に対して数値解を算出する。
【0053】
そして、レオロジーモデル決定部６４では、複数の数値解同士を比較して最も小さい数値解となる前記レオロジーモデルを前記流体のレオロジーモデルとして決定する（Ｓ６：レオロジーモデル決定ステップ）。
【0054】
なお、ニュートン流体など、レオロジーモデルが予め判っている場合は、ステップＳ４における複数種のレオロジーモデルを対象にした数値解の算出やステップＳ５のレオロジーモデル決定ステップの処理を省略してもよい。また、不明であったレオロジーモデルがステップＳ６において決定された後、レオロジーモデルに変更がない（せん断応力τに殆ど変化がない）場合には、その後の処理において、ステップ５における他のレオロジーモデルに基づく数値解の算出ステップＳ５やレオロジーモデル決定ステップＳ６の処理を省略することができる。
【0055】
そして、レオロジー物性決定部６５が、数値解算出部６３によって得られた数値解に基づき流体のレオロジー物性値を決定する（Ｓ７：レオロジー物性決定ステップ）。もし、レオロジー物性値に基づく成分割合等が決定できる場合は、レオロジー物性値とともに成分割合も決定する。レオロジー物性値や成分割合などは、数値解を算出することで決定される。つまり、前記数値解と同様におよそ数秒から数十秒間隔でレオロジー物性値を決定することができる。よって、レオロジー物性値の変動や成分割合の変化などをほぼリアルタイムに計測することができる。
【0056】
以上のような本実施形態の非接触型レオロジー物性計測装置３、システム１、プログラム３ａおよび方法によれば、以下の効果を奏することができる。
１．管内を流れる流体の流速を非接触かつ非侵襲によって計測することのできる流速データを取得することによって、管内の圧力などを計測することなく、前記流体のレオロジー物性を計測することができる。
２．管内の流れにおける脈動を利用するため、適用範囲が広く、様々な管内を流れる流体に適用することができる。
３．費用関数などを用いることで計算にかかる処理負荷を軽くすることができるため、レオロジー物性値の変動や成分割合の変化などをリアルタイムに計測することができる。
４．流速データの関数近似値を用いることで、超音波流速計測装置２により計測される前記流速データの計測誤差に基づくレオロジー物性値等の計測誤差を抑制することができる。
５．ピーク周波数検出部６２では、粘性層厚さに基づきピーク周波数を検出することで、計測原理上の精度を確保することができる。
１．管内を流れる流体の流速を非接触かつ非侵襲によって計測することのできる流速データを取得することによって、管内の圧力などを計測することなく、前記流体のレオロジー物性を計測することができる。
２．管内の流れにおける脈動を利用するため、適用範囲が広く、様々な管内を流れる流体に適用することができる。
３．費用関数などを用いることで計算にかかる処理負荷を軽くすることができるため、レオロジー物性値の変動や成分割合の変化などをリアルタイムに計測することができる。
４．流速データの関数近似値を用いることで、超音波流速計測装置２により計測される前記流速データの計測誤差に基づくレオロジー物性値等の計測誤差を抑制することができる。
５．ピーク周波数検出部６２では、粘性層厚さに基づきピーク周波数を検出することで、計測原理上の精度を確保することができる。
【0057】
次に、本発明に係る非接触型レオロジー物性計測装置、システム、プログラムおよび方法の具体的な実施例について説明する。なお、本発明の技術的範囲は、以下の実施例によって示される特徴に限定されるものではない。
【実施例1】
【0058】
実施例１では、数値的に脈動を伴う管内流れを作り出し、作り出された流れを超音波流速計測装置で計測したと仮定して得られる流速データとして、本発明に係る非接触型レオロジー物性計測プログラムによるレオロジー物性値の算出を行った。また、管内流れにおける厳密解と、本発明に係る非接触型レオロジー物性計測プログラムで得られる数値解とを比較し、計測精度の検討を行った。
【0059】
＜脈動を伴う管内流れの数値計算について＞
ニュートン流体において一方向に流れる管内定常流れはハーゲン・ポアズイユ流れと呼ばれ、以下の支配方程式（１３）および式（１４）で表すことができる。
［数１３］
［数１４］
また、脈動する管内流れは、以下の支配方程式（１５）および式（１６）で表すことができる。
［数１５］
［数１６］
そして、脈動を伴う管内流れは、式（１２）～式（１６）により、以下の式（１７）として得られる。
［数１７］
ここで、
R は管の半径
U0 はハーゲン・ポアズイユ流れの最大値
U1は脈動の速度変動（振幅の大きさ）
ωは脈動の周波数（角周波数）
νは動粘性係数
ρは流体の密度
Ｊ０およびＪ１はベッセル関数、をそれぞれ意味する。
ニュートン流体において一方向に流れる管内定常流れはハーゲン・ポアズイユ流れと呼ばれ、以下の支配方程式（１３）および式（１４）で表すことができる。
［数１３］
［数１４］
また、脈動する管内流れは、以下の支配方程式（１５）および式（１６）で表すことができる。
［数１５］
［数１６］
そして、脈動を伴う管内流れは、式（１２）～式（１６）により、以下の式（１７）として得られる。
［数１７］
ここで、
R は管の半径
U0 はハーゲン・ポアズイユ流れの最大値
U1は脈動の速度変動（振幅の大きさ）
ωは脈動の周波数（角周波数）
νは動粘性係数
ρは流体の密度
Ｊ０およびＪ１はベッセル関数、をそれぞれ意味する。
【0060】
＜脈動を伴う管内の計算条件および計算結果＞
ここで、式（１７）により脈動を伴う管内流れを計算する際の計算条件としてU0 = 0.2 m/s、U1 = 0.1 m/s、r = 1000 kg/m3、R = 25.4 mm、Dt = 10 ms、Dr = 0.1 mmとした。
ここで、式（１７）により脈動を伴う管内流れを計算する際の計算条件としてU0 = 0.2 m/s、U1 = 0.1 m/s、r = 1000 kg/m3、R = 25.4 mm、Dt = 10 ms、Dr = 0.1 mmとした。
【0061】
また、異なる粘性（レオロジー物性）を有する流体の流れを解析対象とするため、動粘性係数νの条件は、n = 10 mm2/s、n = 50 mm2/s、n = 100 mm2/sなどの複数の条件とした。
【0062】
さらに、周波数ω、つまり脈動の速さの違いによる計算精度の影響について検討するためにfo （=ω／２π） =0.5 Hz、fo（=ω／２π）=1.0 Hzなどの複数の条件とした。
【0063】
各条件に基づき計算された脈動を伴う管内の計算結果を図３に示す。各図において、縦軸が管内の位置を表している。縦軸の中央が管の中心位置を表しており、上下の両端が管壁の位置を表している。具体的には、r/Rとして無次元化されており、管の中心位置である縦軸の中央がr/R=0、管壁に相当する値がr/R=1またはr/R=-1である。横軸は、時間の経過を表している。具体的には、tfo として無次元化されている。また、時間毎の管内の流速分布は色の濃淡で表している。色が黒いほど速度はゼロに近づき、色が白いほど速度は速くなる。具体的には、u/U0として無次元化されている。
【0064】
図３に示すように、各時間において中央位置の流速が早く、管壁に近づくにつれ速度が遅くなっている。これは、ハーゲン・ポアズイユ流れのように軸対象となる管内の流速分布を示している。また、例えば、中心位置は時間の経過とともに速度が速くなったり遅くなったりを繰りかえしている。よって、式（１７）により脈動を伴う管内流れを計算されていることが確認できる。
【0065】
＜超音波流速計測装置による計測シミュレーション＞
図３に示す計算結果は、脈動を伴う管内流れを計算した流速分布の厳密解である。一方、管内の流れを超音波流速計測装置を用いて非接触で計測する場合、得られる流速データには計測誤差（計測ノイズ）が含まれる。計測誤差の要因は様々であるが、例えば、流体に含まれており超音波を反射させる固形物の密度や流体内や管壁内における超音波の乱反射などによる影響が考えられる。
図３に示す計算結果は、脈動を伴う管内流れを計算した流速分布の厳密解である。一方、管内の流れを超音波流速計測装置を用いて非接触で計測する場合、得られる流速データには計測誤差（計測ノイズ）が含まれる。計測誤差の要因は様々であるが、例えば、流体に含まれており超音波を反射させる固形物の密度や流体内や管壁内における超音波の乱反射などによる影響が考えられる。
【0066】
そこで、本実施例１では、式（１７）に基づき計算された脈動を伴う管内流れの流速分布に対し計測精度を考慮した計測誤差（計測ノイズ）の含まれる流速分布を下記の式（１８）と定義した。
［数１８］
ここで、
n(ave, std)はノイズ関数、
aはノイズレベル、をそれぞれ意味する。
［数１８］
ここで、
n(ave, std)はノイズ関数、
aはノイズレベル、をそれぞれ意味する。
【0067】
例えば、ノイズレベルa = 5.0とすると、ノイズ関数は、図４に示すよう、真の値を中央値とした正規分布で表される。本実施例１では、この式（１８）を用いて、ノイズレベルaが1.0、5.0、10.0および15.0の場合について、流速データu(r,t)を代入して誤差を含む流速データu'(r,t)を計算した。図５は計算結果である。図５に示すように、ノイズレベルが大きくなるに従い、得られる流速分布の形状の滑らかさがなくなる。つまり、この流速分布は、計測誤差が含まれた流速データと見ることができる。本実施例１では、このように真の値として算出された流速データを計測誤差としてノイズを与えることで、超音波流速計測装置による計測誤差を含む流速データとした。
【0068】
＜本発明により得られる関数式の数値解について＞
次に、式（１７）および式（１８）によって算出された流速データを用いて式（１）により得られる関数式から数値解を算出した。レオロジーモデルは、ニュートン流体とした。そこで、ニュートン流体における脈動を伴う管内流れの費用関数を表す式（７）および式（８）を用いて、前記費用関数が最小化させるプログラムを作成し、当該プログラムにより数値解を算出した。費用関数を算出するアルゴリズムにはランダムサーチ法を用いた。
次に、式（１７）および式（１８）によって算出された流速データを用いて式（１）により得られる関数式から数値解を算出した。レオロジーモデルは、ニュートン流体とした。そこで、ニュートン流体における脈動を伴う管内流れの費用関数を表す式（７）および式（８）を用いて、前記費用関数が最小化させるプログラムを作成し、当該プログラムにより数値解を算出した。費用関数を算出するアルゴリズムにはランダムサーチ法を用いた。
【0069】
本実施例１では、脈動の周波数fo（角周波数ω）が数値解に与える影響について検討するため、周波数fo [Hz]を0.1、0.2、0.5、1.0、2.0、4.0のそれぞれの場合について計算した。
【0070】
図６に算出結果を示す。図の左側が式（７）における実数値Re、右側が虚数値Imの計算結果である。また、上側は動粘性係数n = 100 mm2/sの場合であり、下側は動粘性係数の値がそれよりも１／１０である動粘性係数n = 10 mm2/sの場合である。なお、n = 100 mm2/sは食用油程度の粘性であり、n = 10 mm2/sは水よりもやや高い粘性である。
【0071】
まずは、粘性の影響について検討する。図６の左側の実数値について、周波数fo =0.1[Hz]の場合に着目すると、上側の粘性が高い場合は、その大きさが管壁（r/R＝1.0）の近傍から管の中央（r/R＝0.0）に至るまでなだらかに実数値Reが増加している。これは、粘性が管の中央まで影響していることを示す。一方、下側の粘性が高い場合、管の中央（r/R＝0.0）からr/R＝0.25までの範囲はほぼ一定値である。つまり、管の中央では粘性の影響が弱いことを意味している。
【0072】
図６の右側の虚数値Imについても、周波数fo =0.1[Hz]の場合に着目すると、上側の粘性が高い場合は、実数値Reと同様、管壁（r/R＝1.0）の近傍から管の中央（r/R＝0.0）に至るまでなだらかに増加している。一方、下側の粘性が低い場合は、管の中央（r/R＝0.0）からr/R＝0.7までの範囲はほぼ一定値である。
【0073】
ここで、式（３）における粘性層厚さ（ν/ｋΔω）１/２を検討してみると、粘性の高い動粘性係数n = 100 mm2/sの場合、粘性層厚さは約12.6mmである。管の半径R = 25.4 mmであることから粘性層厚さは管壁（r/R＝1.0）からr/R＝0.5の範囲（管の半径の１／２）程度であった。一方、粘性の低い（さらさらしている）動粘性係数n = 10 mm2/sの場合、粘性層厚さは約4.0mmである。よって、粘性層厚さは管壁（r/R＝1.0）からr/R＝0.85程度であった。
【0074】
なお、粘性層厚さは約4.0mmの場合、本実施例１において式（３）における超音波流速計測装置における計測点同士の距離Δｒが0.1mmであることから、粘性層厚さはΔｒの４０倍程度あり、３Δｒより大きい。よって、粘性の影響の大きい管壁（r/R＝1.0）からr/R＝0.85の範囲内にある約４０点ある計測点のうち３点以上の流速データを用いれば式（１）に基づき十分に精度の良い数値解を算出することができる。
【0075】
次に、脈動の周波数fo（角周波数ω）の影響について検討する。図６の上側の実数値Reについて注目すると、比較的ゆっくりと振動している周波数fo（角周波数ω）が0.1や0.2の場合には、実数値Reは管壁（r/R＝1.0）の近傍から管の中央（r/R＝0.0）に至るまでなだらかに増加している。しかし、周波数fo（角周波数ω）が段々早くなるにつれて管の中央における実数値がほぼ一定値となる範囲が広がる。つまり、周波数fo（角周波数ω）が速くなると粘性の影響が管壁（r/R＝1.0）から離れるに従って伝わりにくくなっている。
【0076】
また、図６の上側の虚数値Imについても、同様に、比較的ゆっくりと振動している場合に数値解である実数値Reおよび虚数値Imは管壁（r/R＝1.0）の近傍から管の中央（r/R＝0.0）に至るまでなだらかに増加している。そして、周波数fo（角周波数ω）が段々早くなるにつれて管の中央における実数値がほぼ一定値となる範囲が広がる。
【0077】
さらに、粘性の低い図６の下側の実数値Reおよび虚数値Imについても、比較的ゆっくりと振動している場合には、管壁（r/R＝1.0）の近傍から管の中央（r/R＝0.0）に向けて比較的、解析値が変化しており、粘性の影響が現れているが、周波数fo（角周波数ω）が速くなると、管壁（r/R＝1.0）の近傍に粘性の影響が現れている。
【0078】
＜レオロジー物性の算出結果について＞
次に、本発明により得られる関数式の数値解に基づきレオロジー物性値を決定する。図７は、費用関数の最小化を行ったときの結果である。図７の高さ方向は算出された粘性係数μ、下の面は算出された実数値Re（αａ）および虚数値Im（αｂ）を表している。そして、各点はランダムサーチ法によりランダムに数値を入れたときの粘性係数μ、実数値および虚数値の組合せを表している。また、各点の色は費用関数の大きさを示しており、色が黒に近づくにつれて費用関数が小さくなっていく。
次に、本発明により得られる関数式の数値解に基づきレオロジー物性値を決定する。図７は、費用関数の最小化を行ったときの結果である。図７の高さ方向は算出された粘性係数μ、下の面は算出された実数値Re（αａ）および虚数値Im（αｂ）を表している。そして、各点はランダムサーチ法によりランダムに数値を入れたときの粘性係数μ、実数値および虚数値の組合せを表している。また、各点の色は費用関数の大きさを示しており、色が黒に近づくにつれて費用関数が小さくなっていく。
【0079】
図７に示すように、ランダムサーチ法により算出された費用関数の値はグラフ上の一点に集中していく。この集中されている点が最小化された費用関数となる。そして、費用関数の最小化された点の高さ方向の値が粘性係数μとして決定することができる。
【0080】
＜数値解の算出精度について＞
次に、数値解の算出結果から得られるレオロジー物性値と、真の値に相当する管内流れの計算に用いたレオロジー物性値との比較を行った。また、本発明では、レオロジー物性値とともに圧力勾配も算出することができる。そこで、数値解の算出結果から得られる圧力勾配と、真の値に相当する管内流れの計算に用いた圧力勾配との比較も行った。
次に、数値解の算出結果から得られるレオロジー物性値と、真の値に相当する管内流れの計算に用いたレオロジー物性値との比較を行った。また、本発明では、レオロジー物性値とともに圧力勾配も算出することができる。そこで、数値解の算出結果から得られる圧力勾配と、真の値に相当する管内流れの計算に用いた圧力勾配との比較も行った。
【0081】
図８は、横軸にノイズレベルを示している。そして、縦軸には、数値解の算出結果から得られる粘性係数μと、真の粘性係数μtrueとの比であるηと、数値解の算出結果から得られる圧力勾配αと、真の圧力勾配αtrueとの比であるλを表している。また、数値解の算出結果から得られる粘性係数μおよび圧力勾配αについては空間的（r方向）に連続する３点の平均値を用いている。縦軸は真の値との比であるため１に近いほど数値解の算出結果から得られるレオロジー物性値および真の値に近いことになる。
【0082】
図８に示すようにノイズレベルが５以下の場合には、数値解の算出結果から得られる結果は、ほぼ真の値と同じ値を示した。一方、ノイズレベルが１０では、η（粘性係数の比）が真の値に対してやや小さい値を示し、ノイズレベルが１５の場合にはさらに真の値からは離れる結果となった。また、λ（圧力勾配の比）については、ノイズレベル１５においてもほぼ真の値と同じ値を示し非常に高い精度であった。
【0083】
よって、本発明に係る非接触型レオロジー物性計測装置、システム、プログラムおよび方法において、超音波流速計測装置による計測精度を保つことで、非常に精度の高い結果が得られることが示された。
【0084】
次に、脈動の周期（半周期内）に対する精度の違いについて確認した。結果を図９に示す。横軸は脈動の周期であり、縦軸はηおよびλである。その結果、ηおよびλはいずれも１に近似しており、周期の違いによる精度には影響は確認できなかった。よって、本発明に係る非接触型レオロジー物性計測装置、システム、プログラムおよび方法では、脈動を伴う管内流れに対してどのタイミングで計測を行っても精度に大きな違いは生じないことがわかった。
【0085】
次に、動粘性係数νに対する精度の違いについて確認した。結果を図１０に示す。横軸は動粘性係数νであり、縦軸はηおよびλである。図１０に示すように、動粘性係数νが高い範囲では精度が落ちている。これは、圧力勾配を一定にしているため、速度変動が弱くなり、結果としてノイズが相対的に大きくなるためと考えられる。一方、動粘性係数νが100 mm2/s 以下においては、真の値とほぼ一致しており精度が高かった。
【0086】
次に、脈動の周波数fo（角周波数ω）に対する精度の違いについて確認した。結果を図１１に示す。横軸は周波数foであり、縦軸はηおよびλである。図１１に示すように、周波数foが低い範囲で精度が落ちている。これは、圧力勾配を一定にしているため、周波数foが低くなると振幅が小さくなり、結果としてノイズが相対的に大きくなるためと考えられる。一方、周波数foが1.0Hz以上においては、真の値とほぼ一致しており精度が高かった。
【0087】
次に、脈動の振幅U1に対する精度の違いについて確認した。結果を図１２に示す。横軸は振幅U1であり、縦軸はηおよびλである。図１２に示すように、振幅U1が低いと精度が落ちている。これは、脈動の周波数fo（角周波数ω）のときと同様に振幅が小さくなり、結果としてノイズが相対的に大きくなるためと考えられる。一方、振幅U1が0.1m/s以上においては、真の値とほぼ一致しており精度が高かった。
【0088】
以上より、本実施例１では、本発明に係る非接触型レオロジー物性計測装置、システム、プログラムおよび方法を用いることで、管内を流れる脈動を伴う流体のレオロジー物性を、非接触で精度よく計測することができることを数値実験により確かめることができた。
【実施例2】
【0089】
実施例２では、脈動を伴う管内流れを作り出す実験装置を作成し、本発明に係る本発明に係る非接触型レオロジー物性計測システムによってニュートン流体および非ニュートン流体のレオロジー物性値の計測を行った。また、カタログ値や市販されている回転式のレオメータによる計測値と、本発明に係る非接触型レオロジー物性計測システムで得られる計測値とを比較し、計測精度の検討を行った。
【0090】
＜脈動を伴う管内流れを作り出す実験装置について＞
本実施例２では、図１３に示すように脈動を伴う管内流を作り出す実験装置を作成した。計測対象となる管は、直線状のステンレス管と、交換可能に設置されるアクリル管とを有する。ステンレス管には、長さは３０００ｍｍ以上、外径５０．８ｍｍ（２インチ）、内径４０．８ｍｍのものを用いた。アクリル管は、ステンレス管の後流側に設置した。本実施例２では、内径約４８ｍｍのものと、内径約２２ｍｍのものを用いた。また、管内の流れを安定させるためステンレス管の入り口から約３０００ｍｍ以上離れた地点を計測位置とした。
本実施例２では、図１３に示すように脈動を伴う管内流を作り出す実験装置を作成した。計測対象となる管は、直線状のステンレス管と、交換可能に設置されるアクリル管とを有する。ステンレス管には、長さは３０００ｍｍ以上、外径５０．８ｍｍ（２インチ）、内径４０．８ｍｍのものを用いた。アクリル管は、ステンレス管の後流側に設置した。本実施例２では、内径約４８ｍｍのものと、内径約２２ｍｍのものを用いた。また、管内の流れを安定させるためステンレス管の入り口から約３０００ｍｍ以上離れた地点を計測位置とした。
【0091】
測定対象となる管の上流側には試験流体を貯留する貯留タンクを接続した。また、下流側には試験流体を流すためのロータリーポンプを接続した。本実施例２では、ロータリーポンプとして、主に食品や医薬品、化粧品などの製造工程で使用されるサニタリーポンプを用いた。また、前記ロータリーポンプはパソコンに接続されており、管内を流れる試験流体の振動数を制御できるようになっている。ロータリーポンプの排出口には貯留タンクに連結された管を設けた。よって、ロータリーポンプを作動させることで、試験流体は実験装置内で循環するようになっている。
【0092】
また、計測位置には、図１４に示すように、超音波流速計測装置の一構成である超音波トランスデューサと、この超音波トランスデューサを角度θ１で保持する保持治具とが設置されている。また、超音波トランスデューサと管と間には、超音波の伝達性を高めるための超音波ジェルを充填させた。
【0093】
超音波流速計測装置は、本発明に係る非接触型レオロジー物性計測装置として構成されるパソコンに接続されており、非接触型レオロジー物性計測システムを構築している。
【0094】
＜ニュートン流体として用いる試験流体について＞
ニュートン流体には、動粘度１０ｃＳｔのシリコーンオイル（粘度（粘性係数）：９．３５×１０－３［Ｐｓ／ｓ］）を用いた。本実施例２においては、シリコーンオイル内に超音波の反射を促すトレーサ粒子を含有させた。
ニュートン流体には、動粘度１０ｃＳｔのシリコーンオイル（粘度（粘性係数）：９．３５×１０－３［Ｐｓ／ｓ］）を用いた。本実施例２においては、シリコーンオイル内に超音波の反射を促すトレーサ粒子を含有させた。
【0095】
＜ニュートン流体の物性値の計測について＞
本実施例２では、０．１Ｈｚの振動流が発生するようにロータリーポンプを制御し、ニュートン流体（シリコーンオイル）のレオロジー物性値の計測を行った。まず、超音波流速計測装置が、前記超音波トランスデューサを介して超音波の送信および受信を行い、受信した超音波を解析することで脈動を伴う管内流の流速分布を計測した。なお、図１４に示すように、管と試験流体との屈折率の違いから超音波トランスデューサの設置角度θ１と超音波の伝播方向の角度θ２とは異なるが、本実施例２ではこの角度の違いを流速分布の計測結果に反映させている。
本実施例２では、０．１Ｈｚの振動流が発生するようにロータリーポンプを制御し、ニュートン流体（シリコーンオイル）のレオロジー物性値の計測を行った。まず、超音波流速計測装置が、前記超音波トランスデューサを介して超音波の送信および受信を行い、受信した超音波を解析することで脈動を伴う管内流の流速分布を計測した。なお、図１４に示すように、管と試験流体との屈折率の違いから超音波トランスデューサの設置角度θ１と超音波の伝播方向の角度θ２とは異なるが、本実施例２ではこの角度の違いを流速分布の計測結果に反映させている。
【0096】
図１５は、超音波流速計測装置により管内の流速分布を計測した結果である。縦軸が管内の位置を表しており、ξ＝０ｍｍが超音波トランスデューサを設置した管の内壁に相当している。流速分布はコンター（等高線）状の色の濃淡で表している。色が白いほど速度はゼロに近づき、色が黒いほど速度は速い。横軸は、時間であり約１０秒間のデータを示している。略中心地点（ξ＝約２８ｍｍ）の流速を時間に沿って見ていくと、この約１０秒間で約１周期の変動をしており、０．１Ｈｚの振動流が計測されていることがわかる。
【0097】
本実施例２では、実施例１と同様に、式（１）により得られる関数式から数値解を算出した。試験流体は、ニュートン流体であるためレオロジーモデルもニュートン流体とした。よって、ニュートン流体における脈動を伴う管内流れの費用関数を表す式（７）および式（８）を用いて、前記費用関数が最小化させるプログラムを作成し、当該プログラムにより数値解を算出した。費用関数を算出するアルゴリズムにはランダムサーチ法を用いた。
【0098】
図１６に費用関数の算出結果を示す。縦軸は粘度μを、横軸は圧力振幅（Pressure amplitude）を示す。ここで圧力振幅（Pressure amplitude）は、圧力の実部と虚部のノルム値である。費用関数の値はコンター（等高線）状の色の濃淡で表している。色が白いほど値が小さく、色が黒いほど値は大きい。図１６に示すように、ランダムサーチ法により算出された費用関数の値はグラフ上の一点に集中していく。この集中されている点が最小化された費用関数となる。そして、費用関数の最小化された点の高さ方向の値が粘性係数μとして決定することができる。
【0099】
図１７は、１０秒間隔で計測した粘度（Viscosity）および圧力振幅（Pressure amplitude）の約１０００秒間のデータを示すものである。また、図中の破線は、試験流体として用いたシリコーンオイルの粘度のカタログ値９．３５×１０－３［Ｐｓ／ｓ］である。図１７に示すように、計測した粘度（Viscosity）の値はカタログ値である９．３５×１０－３［Ｐｓ／ｓ］近辺の値に集約された。
【0100】
よって、本実施例２の非接触型レオロジー物性計測システムは、充分な精度でニュートン流体のレオロジー物性を計測できることが実証できた。
【0101】
＜非ニュートン流体の物性値の計測について＞
次に、非ニュートン流体の物性値についての計測を行った。
次に、非ニュートン流体の物性値についての計測を行った。
【0102】
＜非ニュートン流体として用いる試験流体について＞
非ニュートン流体には、カルボキシメチルセルロースの水溶液（以下、「ＣＭＣ水溶液」という。）を用いた。本実施例２ではＣＭＣ水溶液の濃度を０．５ｗｔ．％とした。カルボキシメチルセルロース（ＣＭＣ）は、一般に増粘剤として食品を含む加工製品に使用されるものである。ＣＭＣ水溶液は、せん断速度が大きくなると粘度が低下する擬塑性流体であって、本実施例２では粘度は約１００～４００［ｍＰａ・ｓ］のものを使用した。
非ニュートン流体には、カルボキシメチルセルロースの水溶液（以下、「ＣＭＣ水溶液」という。）を用いた。本実施例２ではＣＭＣ水溶液の濃度を０．５ｗｔ．％とした。カルボキシメチルセルロース（ＣＭＣ）は、一般に増粘剤として食品を含む加工製品に使用されるものである。ＣＭＣ水溶液は、せん断速度が大きくなると粘度が低下する擬塑性流体であって、本実施例２では粘度は約１００～４００［ｍＰａ・ｓ］のものを使用した。
【0103】
＜非ニュートン流体の物性値の計測について＞
図１８の左側は、超音波流速計測装置により管内の流速分布を計測した結果である。ここでは、図３と同様、縦軸がr/Rとして無次元化された管内の位置であり、r/R=0が管の中心位置を表している。横軸は、tf として無次元化された時間の経過を表しており、約２周期分のデータを示している。流速分布はコンター（等高線）状の色の濃淡で表しており、色が白いほど速度はゼロに近づき、色が黒いほど速度は速い。略中心地点（r/R=0）の流速を時間に沿って見ていくと、約１周期毎に変動をしており振動流が計測されていることがわかる。
図１８の左側は、超音波流速計測装置により管内の流速分布を計測した結果である。ここでは、図３と同様、縦軸がr/Rとして無次元化された管内の位置であり、r/R=0が管の中心位置を表している。横軸は、tf として無次元化された時間の経過を表しており、約２周期分のデータを示している。流速分布はコンター（等高線）状の色の濃淡で表しており、色が白いほど速度はゼロに近づき、色が黒いほど速度は速い。略中心地点（r/R=0）の流速を時間に沿って見ていくと、約１周期毎に変動をしており振動流が計測されていることがわかる。
【0104】
次に、式（９）ないし式（１２）を用いて非ニュートン流体における脈動を伴う管内流れの費用関数を表す前記費用関数が最小化させるプログラムを作成し、当該プログラムにより数値解を算出した。図１８の右側は、r/Rが０～０．６の範囲の流速分布を用いて式（１０）により算出された実数値Reと虚数値Imの計算結果である。実数値Reおよび虚数値Imは、中心位置に近づくにつれてなだらかに増加しており、０～０．６の全範囲において粘性の影響が現れているものと考えられる。
【0105】
図１９および図２０は、ランダムサーチ法により費用関数を算出した結果である。ニュートン流体を計測した結果の図１６と同様に、ランダムサーチ法により算出された費用関数の値はグラフ上の一点に集中していく。そして、費用関数の最小化された点の高さ方向の値が粘性係数μとして決定することができる。
【0106】
図２１は、ひずみ速度γに応じたせん断応力τの値を示すものであり、丸形のプロットは市販のレオメータによる計測値、三角形のプロットが本発明に係る非接触型レオロジー物性計測システムで計測された計測値である。本実施例２では、市販のレオメータとしてAnton Paar社製のMCR102を用いた。
【0107】
図２１に示すように、本実施例２で使用したＣＭＣ水溶液は、市販のレオメータで計測するとひずみ速度γが増加するにつれてせん断応力τも増加している。このときせん断応力τの増加は、ひずみ速度γが大きくなるにつれてやや鈍化するように増加している。せん断応力τとひずみ速度γとの関係がτ＝Ｋγｎのべき乗則で表されると仮定した場合、γがＡからＢの範囲においては、Ｋ＝０．５６９［Ｐａ・ｓ］、ｎ＝０．９７５１であった。
【0108】
これに対し、本発明に係る非接触型レオロジー物性計測システムで計測された計測値は、γがＡからＢの範囲においてＫ＝０．４６７［Ｐａ・ｓ］、ｎ＝０．９９８１であった。市販のレオメータの値と比較すると、ひずみ速度γに対するせん断応力τの増加も同様な値を示した。
【0109】
＜アクリル管による追加実験について＞
本実施例２では、ステンレス管の後流側に設置した交換式のアクリル管を用いて、流速を早めることでよりひずみ速度γの範囲を広範囲にした計測を行った。具体的には、上述するように内径約４８ｍｍのアクリル管と、内径約２２ｍｍのアクリル管を用いた。内径約４８ｍｍのアクリル管を用いた場合のおおよその最大流速は０．１２［ｍ／ｓ］であったのに対し、内径約２２ｍｍのアクリル管を用いた場合のおおよその最大流速は０．５［ｍ／ｓ］であった。そしてr/Rが約０～０．６の範囲の流速分布を用いてレオロジー物性値を算出した。
本実施例２では、ステンレス管の後流側に設置した交換式のアクリル管を用いて、流速を早めることでよりひずみ速度γの範囲を広範囲にした計測を行った。具体的には、上述するように内径約４８ｍｍのアクリル管と、内径約２２ｍｍのアクリル管を用いた。内径約４８ｍｍのアクリル管を用いた場合のおおよその最大流速は０．１２［ｍ／ｓ］であったのに対し、内径約２２ｍｍのアクリル管を用いた場合のおおよその最大流速は０．５［ｍ／ｓ］であった。そしてr/Rが約０～０．６の範囲の流速分布を用いてレオロジー物性値を算出した。
【0110】
図２２は、図２１と同様、ひずみ速度γに応じたせん断応力τの値を示すものであり、丸形のプロットは市販のレオメータによる計測値である。また、三角形のプロットは内径約４８ｍｍを用いた場合、およびバツ形のプロットは内径約２２ｍｍを用いた場合の本発明に係る非接触型レオロジー物性計測システムで計測された計測値である。図２２に示すように、市販のレオメータで計測できた全範囲に対応する計測結果が得られた。また、市販のレオメータの計測結果と比較して、全範囲で同様な値を示した。
【0111】
また、図２３に、ひずみ速度γに応じた粘度μ［Ｐａ・ｓ］についての計測結果を示す。丸形のプロットは市販のレオメータによる計測値である。また、三角形のプロットは内径約４８ｍｍを用いた場合、およびバツ形のプロットは内径約２２ｍｍを用いた場合の本発明に係る非接触型レオロジー物性計測システムで計測された計測値である。図２３に示すように、市販のレオメータで計測されたひずみ速度γが増加するにつれて粘度μは低下しており、擬塑性流体の特徴を示している。同様に、非接触型レオロジー物性計測システムで計測された計測値についても、ひずみ速度γに対して粘度μは低下傾向にある。粘度μとひずみ速度γとの関係をμ＝Ｋγｎ－１のべき乗則で表されると仮定した場合、図２４に示すように、Ｋ＝０．５５９［Ｐａ・ｓ］、ｎ＝０．８９となり擬塑性流体の特徴を示す値が計測されることが実証できた。
【0112】
以上より、本実施例２では、本発明に係る非接触型レオロジー物性計測装置、システム、プログラムおよび方法を用いることで、管内を流れる脈動を伴う流体のレオロジー物性を、非接触で精度よく計測することができることを実際の実験装置により確かめることができた。
【0113】
なお、本発明に係る非接触型レオロジー物性計測装置、システム、プログラムおよび方法は、前述した実施形態に限定されるものではなく、適宜変更することができる。
【0114】
例えば、ピーク周波数検出部６２が、式（３）の範囲内の条件を満たす周波数の範囲のみ周波数解析を行い、その中から最も大きい振幅を示すものをピーク周波数ω０と決定するように処理してもよい。
【符号の説明】
【0115】
１ 非接触型レオロジー物性計測システム
２ 超音波流速計測装置
３ 非接触型レオロジー物性計測装置
３ａ 非接触型レオロジー物性計測プログラム
４ 表示入力手段
５ 記憶手段
６ 演算処理手段
５１ プログラム記憶部
５２ レオロジーモデル式記憶部
５３ しきい値記憶部
６１ 流速取得部
６２ ピーク周波数検出部
６３ 数値解算出部
６４ レオロジーモデル決定部
６５ レオロジー物性決定部
２ 超音波流速計測装置
３ 非接触型レオロジー物性計測装置
３ａ 非接触型レオロジー物性計測プログラム
４ 表示入力手段
５ 記憶手段
６ 演算処理手段
５１ プログラム記憶部
５２ レオロジーモデル式記憶部
５３ しきい値記憶部
６１ 流速取得部
６２ ピーク周波数検出部
６３ 数値解算出部
６４ レオロジーモデル決定部
６５ レオロジー物性決定部
【図1】
【図2】
【図3】
【図4】
【図5】
【図6】
【図7】
【図8】
【図9】
【図10】
【図11】
【図12】
【図13】
【図14】
【図15】
【図16】
【図17】
【図18】
【図19】
【図20】
【図21】
【図22】
【図23】
【図24】