JP4349385B2 - 粒子追跡法を用いた流場計測方法 - Google Patents

Description

本発明は、流体中に混入させたトレーサ粒子の位置追跡を非接触で行い、トレーサ粒子の位置の時間変化を検出することにより流場を計測する粒子追跡法を用いた流場計測方法に関するものである。

従来から、流体中に微小粒子であるトレーサ粒子を混入させ、位置計測装置によって非接触でトレーサ粒子の位置を追跡し、トレーサ粒子の位置の時間変化を検出することにより流体の速度分布を求める技術が知られている。この種の技術において、個々のトレーサ粒子を追跡し、各トレーサ粒子の速度ベクトルを求める技術はＰＴＶ（Particle Tracking Velocimetry）と呼ばれている。ＰＴＶによって得られるトレーサ粒子の速度ベクトルからは速度分布が求められ、さらに速度分布からは、渦度、せん断ひずみ速度、流線などを求めることができる。このように、速度分布を求めたり、速度分布から流体の振る舞いに関する情報を抽出したりすることで流場に関する定量的な計測がなされる。つまり、流場が計測される。

ところで、位置計測装置としては主としてＴＶカメラが用いられており、トレーサ粒子の速度ベクトルを求める方法には、単独のＴＶカメラで撮像した時間差を持つ複数枚の２次元画像において対応するトレーサ粒子の位置関係を求める２次元ＰＴＶと、複数台のＴＶカメラで撮像した画像内で同じトレーサ粒子同士を対応付け、ＴＶカメラの視差を考慮してトレーサ粒子の３次元位置を求める３次元ＰＴＶとが知られている。

ただし、２次元ＰＴＶでは、トレーサ粒子の移動方向がＴＶカメラによる撮像平面に交差しているとトレーサ粒子の速度ベクトルを正確に抽出することができない。そこで、トレーサ粒子の速度ベクトルを求める必要があるときには３次元ＰＴＶが採用される（たとえば、特許文献１参照）。
特許第３５９９９３８号公報

ところで、３次元ＰＴＶではトレーサ粒子の３次元位置を追跡することができるから、速度ベクトルを正確に抽出することができるが、流場の特定箇所における速度分布を検出しようとすれば、多数個のトレーサ粒子について３次元位置を検出しなければならず、個々のトレーサ粒子の追跡が困難になる。一方、追跡が容易になるようにトレーサ粒子の個数を少なくすれば、位置計測装置の空間分解能に対して相対的にトレーサ粒子の個数が少なくなるから、位置計測装置の空間分解能に対して相対的に流場の特定箇所における速度分布の抽出精度が低下することになる。

この種の問題を解決するために、本発明者は、流場のすべての位置において幾何学的に同一の断面が得られる場合に、流場のどの断面においても速度ベクトルの分布が同一とみなせる点に着目した。言い換えると、幾何学的に同一とは、連続した流場において断面の形状および寸法が同一であるだけではなく、速度ベクトルの分布が同一であることを意味している。このような流場は、たとえば、中心軸が一直線であり断面積が一定である管路内に中心軸に沿って形成される流場、円筒形であって中心軸の周りに回転している流場などがある。

速度ベクトルの分布が同一である流場では、各断面で得られた速度ベクトルの分布は他の断面で得られた速度ベクトルの分布と等しいから、トレーサ粒子の追跡により各断面で得られた沿う速度ベクトルを１枚の断面の速度ベクトルとみなすことが可能である。つまり、各トレーサ粒子で求めた速度ベクトルを１枚の断面となる２次元平面に配置することにより、多数個のトレーサ粒子を用いて１枚の断面の速度ベクトルの分布を求めたことと等価になる。

上述のように、幾何学的に同一の断面が得られる流場においては、流場の幾何学的な性質を利用することにより、比較的少数のトレーサ粒子を用いながらも流場を精度よく計測することができる。ただし、流場の１枚の断面における速度ベクトルの分布を求めるだけでは、渦度、せん断ひずみ速度、流線などの諸量を求めることはできない。

本発明は上記事由に鑑みて為されたものであり、その目的は、位置計測装置で計測するトレーサ粒子の個数を比較的少なくしてトレーサ粒子の追跡を容易にしながらも、速度分布を抽出する際には位置計測装置の空間分解能に見合う程度に多数個のトレーサ粒子の位置を用いて速度ベクトルの分布を精度よく抽出し、しかも速度勾配テンソルを求めることにより流場に関する諸量を求めることを可能とする粒子追跡法を用いた流場計測方法を提供することにある。

請求項１の発明は、内側面が撹拌翼の回転軸を中心とする回転体形状に形成された槽であって撹拌翼の回転に伴って撹拌翼の回転軸を中心とする流動を生じる撹拌槽内の流体にトレーサ粒子を混入させ、非接触かつ所定の時間間隔で位置計測を行う位置計測装置によりトレーサ粒子の３次元位置を追跡し、トレーサ粒子の位置の時間変化により撹拌槽内の流場を計測する方法であって、撹拌翼の回転軸の延長方向を一つの座標軸に持つ２次元平面を基準平面として設定し、トレーサ粒子の３次元位置を追跡することにより撹拌槽内の全体で得られる速度ベクトルの原データを回転軸の回りに回転させて原データを基準平面の位置に配置した後、基準平面の上の着目位置の周囲に３次元の微小領域を設定するとともに微小領域内における基準平面の上の複数個の原データを用いて着目位置に関する微分演算を有限差分法で代用することにより着目位置の速度勾配テンソルの成分を求めることを特徴とする。

この方法によれば、流体に撹拌翼の回転軸を中心とする流動が生じていることを利用し、トレーサ粒子と回転軸の延長方向とを含む平面に対するトレーサ粒子の速度ベクトルの向きが、撹拌翼の回転方向における位置には依存しないとみなしている。言い換えると、トレーサ粒子の３次元位置の追跡で得られる速度ベクトルのような原データは、回転軸の回りに回転させても情報が失われることなく保存されることに着目している。したがって、回転軸を一つの座標軸に持つ２次元平面である基準平面を設定し、原データを回転軸の軸回りに回転させることにより設定した基準平面に原データを配置することにより、原データが持つ情報を保存しながらも、基準平面に多数の原データを集約して配置することができる。つまり、位置計測装置の空間分解能に対して流体に混入させるトレーサ粒子の個数を比較的少なくすることができるから、位置計測装置において計測したトレーサ粒子の位置の追跡が容易になる。しかも、１枚の基準平面に原データを集めて配置するから、流場を計測する際には基準平面に集められた多数個の原データを用いることができ、位置計測装置の空間分解能に見合う程度に多数個のトレーサ粒子を用いて速度分布を精度よく抽出することができる。

ただし、流場の速度勾配テンソルは、１枚の基準平面の上の速度ベクトルの成分のみでは求めることができない。そこで、基準平面の上の着目位置の周囲に３次元の微小領域を設定し、微小領域内の基準平面の上の複数個の原データを用いて着目位置に関する微分演算を有限差分法で代用している。微小領域内の任意の位置の速度ベクトルは、基準平面の上の速度ベクトルの原データを複数個用いた補間演算によって求めることができる。このことを利用して、着目位置の速度勾配テンソルに含まれる微分値を有限差分法で求めた差分値で代用し、基準平面上の任意の位置の速度勾配テンソルを求めることが可能になる。つまり、流場に関する諸量を求めることが可能になる。

請求項２の発明では、請求項１の発明において、前記微小領域が前記基準平面を含む立方体の単位格子であって、有限差分法で用いる成分のうち基準平面に含まれない速度ベクトルを用いる成分には、前記回転軸の周りに基準平面を回転させたときの対応点の速度ベクトルの成分を適用することを特徴とする。

この方法によれば、微小領域を立方体の単位格子としているから、微小領域の形状が単純であって任意の位置の速度ベクトルを簡単な演算で求めることができる。また、基準平面で得られた速度ベクトルと基準平面との関係は、基準平面を回転軸の周りに回転させても保存されるから、速度勾配テンソルの成分を求める演算には、基準平面を回転軸の周りに回転させたときの対応点の速度ベクトルの成分を求めればよく、このことによって、撹拌槽内の任意位置の速度ベクトルを求めることが可能になる。その結果、撹拌槽内の所望の着目位置について速度勾配テンソルを求めることが可能になる。

請求項３の発明では、請求項２の発明において、前記単位格子の一面が前記基準平面であることを特徴とする。

この方法によれば、単位格子の一面が基準平面に含まれるから、単位格子の６面のうち基準平面に並行な一面における速度ベクトルの成分を基準平面の上の速度ベクトルから求めるだけで、着目位置に関する微分値（差分値）を求めることができ、比較的少ない演算量で速度勾配テンソルの成分を求めることができる。

請求項４の発明では、請求項２の発明において、前記着目位置が前記単位格子の中心であって、単位格子の６面のうちの４面の中心の位置の速度ベクトルの成分は前記基準平面の上で前記原データにより求め、単位格子の６面のうちの残りの２面の中心の位置の速度ベクトルの成分は前記回転軸の周りに基準平面を回転させたときの対応点の速度ベクトルの成分により求めることを特徴とする。

この方法によれば、着目位置を中心とする４点については基準平面内で求めることができるから、着目位置に関する微分値（差分値）の演算量を低減することができ、速度勾配テンソルの成分を処理負荷の少ない高速な演算で求めることが可能になる。

請求項５の発明では、請求項１ないし請求項４のいずれかの発明において、前記基準平面の上で求めた速度勾配テンソルの成分を回転軸の周りに回転させることにより撹拌槽内の所望の位置の流場を計測することを特徴とする。

この方法によれば、撹拌槽内であれば、どの位置であっても速度勾配テンソルの成分を求めることができるから、所望位置のせん断歪み速度や渦度のほか、流線も求めることが可能である。

本発明の方法によれば、回転軸を一つの座標軸に持つ２次元平面である基準平面を設定し、原データを回転軸の軸回りに回転させることにより設定した基準平面に原データを配置することにより、原データが持つ情報を保存しながらも、基準平面に多数の原データを集約して配置することができ、位置計測装置の空間分解能に対して流体に混入させるトレーサ粒子の個数を比較的少なくすることができるから、位置計測装置において計測したトレーサ粒子の位置の追跡が容易になるという利点がある。しかも、１枚の基準平面に原データを集めて配置するから、流場を計測する際には基準平面に集められた多数個の原データを用いることができ、位置計測装置の空間分解能に見合う程度に多数個のトレーサ粒子を用いて速度分布を精度よく抽出することができるという利点がある。さらに、基準平面の上の着目位置の周囲に３次元の微小領域を設定するとともに、微小領域内の基準平面の上の複数個の原データを用いて着目位置に関する微分演算を有限差分法で代用するから、基準平面上の任意の位置の速度勾配テンソルを求めることが可能になり、結果的に、流場に関する諸量を求めることが可能になるという利点がある。

以下に説明する実施形態では、図１に示すように、有底円筒状であって軸方向を上下方向に一致させて設置される撹拌槽１の中に、複数枚（図示例では６枚）の羽根２ａを有した撹拌翼２が配置される。撹拌翼２には撹拌槽１の軸方向に沿った中心線上に回転軸３が設けられる。回転軸３は上端部に結合されたモータ４により回転駆動力が与えられる。

撹拌槽１の内側面は回転軸３の軸方向に直交する断面が円形であるから、撹拌槽１の内側面は撹拌翼２の回転軸３を中心とする回転体形状に形成されていることになる。回転体形状としては、軸方向における直径が変化しない円筒状のほか、軸方向の位置に応じて直径が変化する形状を採用することも可能である。また、撹拌翼２は回転軸３を対称軸として対称性を有するように形成されている。具体的には、回転軸３の下端部から回転軸３の回転方向における等角度間隔で羽根２ａが放射状に突設される。したがって、撹拌槽１に流体７を注入し、モータ４により回転軸３を回転させると、撹拌槽１の中の流体７に撹拌翼２の回転軸３を中心とする流動が生じる。

本実施形態では３次元ＰＴＶの技術を用いて撹拌槽１の中の流体７に生じる流場を計測するるために、撹拌槽１の中の流体７には微粒子（粒子径が０．５〜１５０μｍ）であるトレーサ粒子が混入される。また、トレーサ粒子の３次元位置を計測するために位置計測装置５が設けられる。位置計測装置５は、２台のＴＶカメラ５ａ，５ｂを用いてトレーサ粒子を所定の時間間隔で撮像する。ＴＶカメラ５ａ，５ｂで撮像する時間間隔はトレーサ粒子の移動速度に応じて設定され（１秒間に数十枚から数千枚の画像が得られるように時間間隔が設定される）、撮像された画像は、動画像のように、異なる時刻に撮像した複数枚の静止画像を時系列に並べたものになる。また、両ＴＶカメラ５ａ，５ｂの視差を利用して多数のトレーサ粒子の３次元位置を非接触で計測する。すなわち、位置計測装置５ではトレーサ粒子について撹拌槽１の中での３次元位置の情報を持つ３次元画像を生成する。

ＴＶカメラ５ａ，５ｂは撹拌槽１の外部に配置されており、撹拌槽１の外部から撹拌槽１の内部のトレーサ粒子を撮像するから、撹拌槽１は透明材料により形成される。また、３次元位置は回転軸３に並行な方向をｚ軸とする直交座標系において求める。この場合、各ＴＶカメラ５ａ，５ｂを、光軸がｚ軸に直交し、かつ両ＴＶカメラ５ａ，５ｂの光軸が互いに直交するように配置すれば、各ＴＶカメラ５ａ，５ｂで撮像された画像の水平方向をそれぞれｘ方向，ｙ方向に対応付けることができる。

位置計測装置５には流場解析装置６が接続される。流場解析装置６は、コンピュータを用いた一種の画像処理装置であって、位置計測装置５で生成された３次元画像を対象として以下の処理を行う。３次元画像は所定時間ごとに得られているから、各３次元画像の中でトレーサ粒子の対応付けを行い、各トレーサ粒子を追跡する必要がある。トレーサ粒子の追跡には、時間軸方向において隣接する２枚ずつの３次元画像の間でトレーサ粒子を対応付ける。トレーサ粒子を追跡する技術については種々提案されているので、適宜の技術を採用すればよい。

流場解析装置６では、トレーサ粒子の追跡後に各トレーサ粒子ごとに速度ベクトルを原データとして求める。以下では、ベクトルとテンソルとの少なくとも一方を含む数式を表記する際に、ベクトルまたはテンソルを［Ａ］という形式で表記する。

いま、時刻ｔ_１と時刻ｔ_２（＞ｔ_１）とにおける３次元画像において、１個のトレーサ粒子について３次元位置がそれぞれ［Ａ］（ｔ_１）＝（ｘ_１，ｙ_１，ｚ_１）、［Ａ］（ｔ_２）＝（ｘ_２，ｙ_２，ｚ_２）であったとし、速度ベクトル［Ｕ］（ｔ）の時刻として始点位置の時刻ｔ_１を採用し、Δｔ＝ｔ_２−ｔ_１とおけば、速度ベクトル［Ｕ］（ｔ_１）は、次式で表される。
［Ｕ］（ｔ_１）＝｛［Ａ］（ｔ_２）−［Ａ］（ｔ_１）｝／Δｔ
＝（（ｘ_２−ｘ_１）／Δｔ，（ｙ_２−ｙ_１）／Δｔ，（ｚ_２−ｚ_１）／Δｔ）
ところで、本実施形態における流場は撹拌翼２の回転軸３を中心とする流動を生じるから、回転軸３を中心とする撹拌翼２の回転方向におけるどの位置であっても、回転軸３を含む断面に対する速度ベクトルは等価に扱うことができる。つまり、回転軸３を中心として回転方向のどの位置においても幾何学的には同一であると言える。そこで、回転方向の位置成分を除去しやすいように、トレーサ粒子の３次元位置を直交座標系から円筒座標系に変換する。この変換を容易にするために、直交座標系における原点の位置は回転軸３の上に設定しているものとする。

なお、とくに限定する趣旨ではないが、撹拌槽１と撹拌翼２との寸法関係は、以下のように設定した。すなわち、撹拌槽１の内径は２００ｍｍ、撹拌槽１内の流体７の高さは２００ｍｍ、撹拌翼２の外径は５０ｍｍ、撹拌翼２の回転軸３に沿う方向の幅は１０ｍｍ、撹拌槽１の内底面から撹拌翼２の下面までの高さは１００ｍｍとした。この条件では、回転軸３の回転方向のすべての位置において幾何学的に同一であるという条件を満たすことができる。つまり、撹拌槽１の内部は、回転軸を含むすべての断面において断面の形状および寸法が同一であり、また各断面における速度ベクトルの分布が同一になる。

時刻ｔにおけるトレーサ粒子の３次元位置を直交座標系で表したときに、［Ａ］（ｔ）＝（ｘ（ｔ），ｙ（ｔ），ｚ（ｔ））であり、円筒座標系で表したときに［Ａ］（ｔ）＝（ｒ（ｔ），θ（ｔ），ｚ（ｔ））であるとすれば、数１の関係が成立する。

上述したように、撹拌槽１の中の流場が回転方向のすべての位置において幾何学的に同一であって、速度ベクトルが回転軸３の回転方向の位置の影響を受けないとみなせる場合には、円筒座標系における角度成分θを無視することができる。つまり、回転軸３の回転方向の各位置の速度ベクトルについて、回転方向の各位置においてｚ軸を含むように設定した２次元平面に対する情報にのみ着目すればよいと言える。

そこで、図２、図３に示すように、ｚ軸を一つの座標軸に持つ基準の２次元平面である基準平面ＰＬを設定し、回転軸３の回転方向の各位置においてｚ軸を含む断面に対する各位置の速度ベクトルの関係を、基準平面ＰＬとの関係に置き換える。つまり、基準平面ＰＬの位置に、原データである速度ベクトルを配置するように速度ベクトル［Ｕ］（ｔ）の座標を変換する。座標変換にあたっては、速度ベクトルの始点を基準平面ＰＬの上に位置させる。時刻ｔ_１の速度ベクトル［Ｕ］（ｔ_１）であれば、始点である３次元位置［Ａ］（ｔ_１）を基準平面ＰＬの上に位置するように、着目するトレーサ粒子について座標変換を行う。

ここでは、基準平面ＰＬをｘｚ平面とし、回転軸３の回転方向の角度θをｘ軸の正の向きからｙ軸の正の向きに向かう回転方向で求めるものとする。この場合、時刻ｔ_１において、着目するトレーサ粒子の円筒座標系での位置が（ｒ（ｔ_１），θ（ｔ_１），ｚ（ｔ_１））であるとすれば、速度ベクトル［Ｕ］（ｔ_１）を角度θ（ｔ_１）だけ回転させることにより、基準平面ＰＬに位置させることができる。数１によれば、θ（ｔ_１）＝ｔａｎ^−１（ｙ（ｔ_１）／ｘ（ｔ_１））であるから、時刻ｔについて一般化し、角度がθ（ｔ）である位置の速度ベクトル［Ｕ］（ｔ）＝（ｕ，ｖ，ｗ）について、基準平面ＰＬ上の速度ベクトル［Ｕ］_ｂへの座標変換を行うとすれば、数２のようになる。なお、速度ベクトル［Ｕ］（ｔ）のうちｚ成分ｗについては座標変換の影響を受けない。

数２に示す座標変換を、すべてのトレーサ粒子について行えば、回転軸３の回転方向におけるトレーサ粒子の位置にかかわらず、着目するトレーサ粒子の速度ベクトル［Ｕ］（ｔ）を基準平面ＰＬに対する速度ベクトル［Ｕ］_ｂとして扱うことが可能になる。その結果、たとえば、１個のトレーサ粒子に着目するだけでも基準平面ＰＬの上では多数個の速度ベクトル［Ｕ］_ｂとして扱うことが可能になる。また、原データである速度ベクトルを基準平面ＰＬに集約することにより多数個のデータとして扱うことが可能になるから、ＴＶカメラ５ａ，５ｂの分解能に対するトレーサ粒子の個数を減らしてトレーサ粒子の追跡を容易にしながらも、流場の計測に用いるデータの密度を高めることによって、流場を精度よく計測することができるようになる。速度ベクトルからは、回転軸３を含む断面内での流体７の速度分布を求めることができ、また、速度分布の微分演算や積分演算によって、せん断ひずみ速度、流線などを求めることができる。これらの諸量を求める手順については後述する。

なお、上述した方法では、速度ベクトルを基準平面ＰＬに配置するから、回転軸３の回転方向における速度ベクトルの情報は失われるが、回転軸３に直交する断面内あるいは回転軸３に沿った断面内での速度分布の情報は保たれている。速度分布、せん断ひずみ速度、流線などを撹拌槽１の内部空間の各位置に対応付けるには、基準平面ＰＬの上での演算の後に、後述するように、撹拌槽１の各空間へのデータの再配置を行う。

また、上述の例では速度ベクトル［Ｕ］（ｔ）の始点位置を基準平面ＰＬに一致させる例を示したが、基準平面ＰＬに速度ベクトル［Ｕ］（ｔ）の終点を一致させてもよい。この場合、数２の角度θ（ｔ_１）として、ｔａｎ^−１（ｙ（ｔ_２）／ｘ（ｔ_２））を用いる。

上述した動作を図４に簡単にまとめる。まず、撹拌槽１内のトレーサ粒子について位置計測装置５を用いて３次元位置を計測する（Ｓ１）。得られた３次元画像において着目するトレーサ粒子を追跡し（Ｓ２）、さらに、回転軸３の回転方向におけるトレーサ粒子の位置と基準平面ＰＬとの間の角度差θ（ｔ）を求める（Ｓ３）。また、着目するトレーサ粒子について速度ベクトル［Ｕ］（ｔ）を求める（Ｓ４）。速度ベクトル［Ｕ］（ｔ）を基準平面ＰＬの上に配置するように座標変換を施す（Ｓ５）。基準平面ＰＬに配置された速度ベクトル［Ｕ］_ｂにより流場を計測し（Ｓ６）、速度分布やせん断ひずみ速度などを求める。基準平面ＰＬ上での流場計測の結果を撹拌槽１内の空間に展開することにより、撹拌槽１の中の流場を出力する（Ｓ７）。この結果は、コンピュータのモニタ装置のような出力装置に出力される。

ところで、基準平面ＰＬの上に配置された原データ（速度ベクトル）を用いて流場に関する諸量を求めるには、原データを配置した基準平面ＰＬを諸量の演算が可能となるように撹拌槽１の適宜の位置に適用する。撹拌槽１の内部は幾何学的に同一であるから、回転軸３の周りで基準平面ＰＬをどの位置に回転させても、該位置の基準平面ＰＬの上での速度ベクトルの分布として用いることができる。

そこで、撹拌槽１内の流場に関する諸量を求めようとするときには、原データを配置した基準平面ＰＬを回転軸３の周りに回転させることによって、諸量を求めようとする所望位置に位置させる。このように基準平面ＰＬを所望位置に位置させることにより、ｘｚ平面以外の位置における諸量を求めることができる。つまり、基準となるｘｚ平面上に配置した原データに対して数２に示した座標変換の逆変換を行うのであって、所望位置の基準平面ＰＬと基準とする基準平面ＰＬとの角度差を求め、逆変換の際にこの角度差を適用することにより所望位置における速度ベクトルを求めることができる。

以下では、流場に関する諸量の一例としてせん断ひずみ速度を求める場合を例示する。せん断ひずみ速度は、速度勾配テンソル（トレーサ粒子の速度ベクトルの勾配＝grad［Ｕ］_ｂ）を、対称テンソルと非対称テンソルとの和の形式で表し、さらにこの対称テンソルを対角成分からなるテンソルと非対角成分からなるテンソルとの和で表したとき、非対角成分からなるテンソルと微小変位とのスカラー積に相当する。言い換えると、速度勾配テンソルに０ではない成分が含まれているときには、流場の近接した２点の速度が異なることを意味するのであって、当該２点を含む微小領域は時間の経過に伴って、形状および体積が変化しない変化（回転）、形状のみの変化（ずれ歪み）、形状および体積の変化（伸縮歪み）のいずれかを生じることになる。回転は非対称テンソルに対応し、ずれ歪みは対称テンソルの非対角成分からなるテンソルに対応し、伸縮歪みは対称テンソルの対角成分からなるテンソルに対応する。ずれ歪みに対応するテンソルを求めるとせん断ひずみ速度を求めることができ、回転に対応するテンソルを求めると渦度を求めることができる。

これらの関係を数式を用いて簡単に説明する。着目する流場については連続体として考える。流場の近接した２点の速度について考察するために、各点の位置をそれぞれ位置ベクトル［Ｘ］，［Ｘ＋ΔＸ］で表す。各点の速度ベクトルをそれぞれ［Ｕ］_ｂ，［Ｕ＋ΔＵ］_ｂとすれば、速度ベクトルの差［ｄＵ］_ｂは次式で表すことができる。
［ｄＵ］_ｂ＝［Ｕ（［Ｘ＋ΔＸ］）］_ｂ−［Ｕ（［Ｘ］）］_ｂ＝grad［Ｕ］_ｂ・［ｄＸ］
上述のように、grad［Ｕ］_ｂは速度勾配テンソルであって、テンソルの性質から対称成分と非対称成分の和として表すことができる。すなわち、grad［Ｕ］_ｂ＝［Ｌ］とおけば、次式が成立するから、次式の第１項と第２項とをそれぞれ［Ｄ］，［Ｗ］おけば、［Ｄ］は速度勾配テンソル［Ｌ］の対称テンソル（「伸長テンソル」または「変形速度テンソル」という）、［Ｗ］は速度勾配テンソル［Ｌ］の非対称テンソル（「渦度テンソル」または「スピンテンソル」という）になる。
［Ｌ］＝（１／２）（［Ｌ］＋［Ｌ］^Ｔ）＋（１／２）（［Ｌ］−［Ｌ］^Ｔ）
［Ｄ］＝（１／２）（［Ｌ］＋［Ｌ］^Ｔ）
［Ｗ］＝（１／２）（［Ｌ］−［Ｌ］^Ｔ）
具体的には、伸長テンソル［Ｄ］、渦度テンソル［Ｗ］は数３のように表される。

上述したように渦度テンソル［Ｗ］は変形を伴わない回転を意味しており、伸長テンソル［Ｄ］は変形を伴う伸縮歪みとずれ歪みとの和を意味している。伸長テンソル［Ｄ］は対角成分のテンソル［Ｄ１］と非対角成分のテンソル［Ｄ２］との和として表すことができる（すなわち、［Ｄ］＝［Ｄ１］＋［Ｄ２］）。対角成分からなるテンソル［Ｄ１］は単純伸長変形速度テンソル、非対角成分からなるテンソル［Ｄ２］はせん断変形速度テンソル（または、単純ずり変形速度テンソル）と呼ばれる。単純伸長変形速度テンソル［Ｄ１］とせん断変形速度テンソル［Ｄ２］とはそれぞれ数４のように表される。非対角成分によるせん断変形では、体積は変化しない。

ここまでの説明を簡単にまとめる。まず、［ｄＵ］_ｂ＝grad［Ｕ］_ｂ・［ｄＸ］＝［Ｌ］・［ｄＸ］であって、［Ｌ］＝［Ｄ］＋［Ｗ］＝［Ｄ１］＋［Ｄ２］＋［Ｗ］であるから、次式が得られる。
［ｄＵ］_ｂ＝（［Ｄ１］＋［Ｄ２］＋［Ｗ］）・［ｄＸ］
上式はスカラー積であるから、展開することができ、次式が得られる。
［ｄＵ］_ｂ＝［Ｄ１］・［ｄＸ］＋［Ｄ２］・［ｄＸ］＋［Ｗ］・［ｄＸ］
ここで、［Ｄ２］・［ｄＸ］＝［ｄＵ２］_ｂとおけば、［ｄＵ２］_ｂがせん断ひずみ速度になる。したがって、せん断ひずみ速度［ｄＵ２］_ｂを算出するには、数４に示したせん断変形速度テンソル［Ｄ２］の成分を演算する必要がある。ただし、トレーサ粒子により計測されるデータは基準平面ＰＬに配置した速度ベクトルの原データであるから、数４の演算に原データを用いるための工夫が必要になる。

いま、基準平面ＰＬ上では速度ベクトルの原データが十分に高い密度で得られていると仮定すると、図５のように、基準平面ＰＬの上で正方形の各頂点の位置に始点を持つ４個の速度ベクトルを選択することができると考えることができる（実際には４個の速度ベクトルの始点の位置が正方形の頂点に位置するとは限らないが、ここでは、理想化して説明する。また、始点が正方形の頂点に位置する４個の原データが得られないときには、後述する補間演算を行うことにより、始点が正方形の頂点に位置する４個の速度ベクトルを用いることができる）。以下では、これらの４個の速度ベクトルの原データを用い、正方形の中心の位置を着目位置として、着目位置におけるせん断変形速度テンソル［Ｄ２］の各成分を求める方法について説明する。

せん断変形速度テンソル［Ｄ２］の各成分には、∂ｕ／∂ｘ、∂ｕ／∂ｙ、∂ｕ／∂ｚ、∂ｖ／∂ｘ、∂ｖ／∂ｙ、∂ｖ／∂ｚ、∂ｗ／∂ｘ、∂ｗ／∂ｙ、∂ｗ／∂ｚの９種類の偏微分の要素が含まれているから、これらの各要素を４個の原データを用いて表せば、数４によりせん断変形速度テンソル［Ｄ２］を求めることができる。ただし、原データは離散値であるから、偏微分の値は差分で代用することになる。偏微分を差分で代用する演算方法には数種類が知られており、以下では演算が単純で負荷が少ない方法として、１次元の有限差分法のうちの中心差分法を用いる例を示すが、他の演算方法を採用することも可能である。

いま、３次元空間を立方体の単位格子で区分し、各単位格子の格子点の位置の座標を（ｉ，ｊ，ｋ）で表すものとする。ｉ，ｊ，ｋは原点位置からｘ，ｙ，ｚの各軸方向における単位格子の個数に相当する整数値であり、各軸方向の正の向きでは正の値になり、負の向きでは負の値になる。したがって、原点位置ではｉ，ｊ，ｋ＝０になる。着目する位置に関する偏微分の値を求めるにあたり、１次元有限差分法では各軸方向の一直線上に並ぶ値を用い、中心差分法では着目する位置の前後の値を用いて中心の値を求める。

たとえば、着目位置の座標を（ｉ，ｊ，ｋ）とし、この座標（ｉ，ｊ，ｋ）での関数ｆのｘに関する偏微分∂ｆ／∂ｘを差分で代用すれば次式で表すことができる。
（∂ｆ／∂ｘ）_ｉ＝｛ｆ（ｉ＋１，ｊ，ｋ）−ｆ（ｉ−１，ｊ，ｋ）｝／２Δｘ
ただし、Δｘはｘ軸方向における単位格子の格子定数であって、Δｘ＝（ｉ＋１）−ｉ＝ｉ−（ｉ−１）＝……＝１である。同様にして、∂ｆ／∂ｙ、∂ｆ／∂ｚは、それぞれ次式で求めることができる。
（∂ｆ／∂ｙ）_ｊ＝｛ｆ（ｉ，ｊ＋１，ｋ）−ｆ（ｉ，ｊ−１，ｋ）｝／２Δｙ
（∂ｆ／∂ｚ）_ｋ＝｛ｆ（ｉ，ｊ，ｋ＋１）−ｆ（ｉ，ｊ，ｋ−１）｝／２Δｚ

ここで、基準平面ＰＬの上の着目位置における速度勾配テンソルを求めるものとし、上述した単位格子を着目位置の周囲に微小領域として設定するものとする。微小領域としての単位格子は、着目位置の微分演算による微分値を有限差分法により求めた差分で代用するために設定した領域であって、図５に示すように、ｘｚ平面である基準平面ＰＬの上の着目位置を中心とする立方体として設定される。ここでは、単位格子を基準平面ＰＬで切り取った断面である正方形の頂点の位置の速度ベクトルが原データから求められるものとする。つまり、単位格子の６面のうちの２面は基準平面ＰＬに平行であって、この２面に平行で単位格子の中心を中心とする正方形の頂点における速度ベクトルを原データから求めるものとする（この速度ベクトルは原データであることが望ましいが、他の原データから補間演算により求めてもよい）
求めた４個の速度ベクトルを上式に適用すれば、以下の演算によりｘに関する偏微分の値とｚに関する偏微分の値とを求めることができる。

いま、時刻については考慮しないこととし（つまり、同時刻について考えることとし）、座標（ｉ，ｊ，ｋ）における速度ベクトルを［Ｕ］_{ｉ，ｊ，ｋ}と表記し、速度ベクトル［Ｕ］_{ｉ，ｊ，ｋ}の各成分を、ｕ_{ｉ，ｊ，ｋ}，ｖ_{ｉ，ｊ，ｋ}，ｗ_{ｉ，ｊ，ｋ}と表記するものとして、着目位置に関する演算を行うために選択した４個の原データが［Ｕ］_{ｉ，ｊ，ｋ}、［Ｕ］_{ｉ＋１，ｊ，ｋ}、［Ｕ］_{ｉ，ｊ，ｋ＋１}、［Ｕ］_{ｉ＋１，ｊ，ｋ＋１}であるとする。つまり、原データが単位格子の格子定数の間隔で求められるものとする。

この場合、着目位置である正方形の中心の座標は（ｉ＋０．５，ｊ，ｋ＋０．５）であり、単位格子の面心（一面の中心）の位置が着目位置になる。なお、基準平面ＰＬはｙ＝０であるからｊ＝０である。原データは単位格子の格子定数の間隔で求められると仮定しているから、着目位置の座標に対するｘ軸方向の前後の座標は、（ｉ，ｊ，ｋ＋０．５）と（ｉ＋１，ｊ，ｋ＋０．５）とになり、着目位置の座標に対するｚ軸方向の前後の座標は（ｉ＋０．５，ｊ，ｋ）と（ｉ＋０．５，ｊ，ｋ＋１）とになる。

したがって、着目位置の座標での速度ベクトル［Ｕ］_{ｉ＋０．５，ｊ，ｋ＋０．５}のｘ，ｚに関する偏微分（差分）の値は、それぞれ次式で与えられる。
∂［Ｕ］_{ｉ＋０．５，ｊ，ｋ＋０．５}／∂ｘ＝｛［Ｕ］_{ｉ＋１，ｊ，ｋ＋０．５}−［Ｕ］_{ｉ，ｊ，ｋ＋０．５}｝／２Δｘ
∂［Ｕ］_{ｉ＋０．５，ｊ，ｋ＋０．５}／∂ｚ＝｛［Ｕ］_{ｉ＋０．５，ｊ，ｋ＋１}−［Ｕ］_{ｉ＋０．５，ｊ，ｋ}｝／２Δｚ
基準平面ＰＬの上の速度ベクトルの原データの間の距離は等間隔とは限らないから、実際にはΔｘ、Δｚは原データの位置によって異なるが、後述するように補間演算を行うことによりΔｘ＝Δｚ＝１として演算することができる。

ところで、基準平面ＰＬの上の速度ベクトルの関係は以下のように考えることができる。
［Ｕ］_{ｉ，ｊ，ｋ＋０．５}＝（［Ｕ］_{ｉ，ｊ，ｋ＋１}＋［Ｕ］_{ｉ，ｊ，ｋ}）／２
［Ｕ］_{ｉ＋１，ｊ，ｋ＋０．５}＝（［Ｕ］_{ｉ＋１，ｊ，ｋ＋１}＋［Ｕ］_{ｉ＋１，ｊ，ｋ}）／２
［Ｕ］_{ｉ＋０．５，ｊ，ｋ}＝（［Ｕ］_{ｉ＋１，ｊ，ｋ}＋［Ｕ］_{ｉ，ｊ，ｋ}）／２
［Ｕ］_{ｉ＋０．５，ｊ，ｋ＋１}＝（［Ｕ］_{ｉ＋１，ｊ，ｋ＋１}＋［Ｕ］_{ｉ，ｊ，ｋ＋１}）／２
それゆえ、次式が成立する。
∂［Ｕ］_{ｉ＋０．５，ｊ，ｋ＋０．５}／∂ｘ＝｛（［Ｕ］_{ｉ＋１，ｊ，ｋ＋１}＋［Ｕ］_{ｉ＋１，ｊ，ｋ}）／２−（［Ｕ］_{ｉ，ｊ，ｋ＋１}＋［Ｕ］_{ｉ，ｊ，ｋ}）／２｝／２Δｘ
∂［Ｕ］_{ｉ＋０．５，ｊ，ｋ＋０．５}／∂ｚ＝｛（［Ｕ］_{ｉ＋１，ｊ，ｋ＋１}＋［Ｕ］_{ｉ，ｊ，ｋ＋１}）／２−（［Ｕ］_{ｉ＋１，ｊ，ｋ}＋［Ｕ］_{ｉ，ｊ，ｋ}）／２｝／２Δｚ

たとえば、速度ベクトルの成分ｕに着目すると、着目位置の座標（ｉ＋０．５，ｊ，ｋ＋０．５）におけるｘに関する微分値∂ｕ_{ｉ＋０．５，ｊ，ｋ＋０．５}／∂ｘは、次式で表される。
∂ｕ_{ｉ＋０．５，ｊ，ｋ＋０．５}／∂ｘ＝｛（ｕ_{ｉ＋１，ｊ，ｋ＋１}＋ｕ_{ｉ＋１，ｊ，ｋ}）／２−（ｕ_{ｉ，ｊ，ｋ＋１}＋ｕ_{ｉ，ｊ，ｋ}）／２｝／２Δｘ
同様にして、着目位置の座標（ｉ＋０．５，ｊ，ｋ＋０．５）におけるｚに関する微分値∂ｕ_{ｉ＋０．５，ｊ，ｋ＋０．５}／∂ｚは、次式で表される。
∂ｕ_{ｉ＋０．５，ｊ，ｋ＋０．５}／∂ｚ＝｛（ｕ_{ｉ＋１，ｊ，ｋ＋１}＋ｕ_{ｉ，ｊ，ｋ＋１}）／２−（ｕ_{ｉ＋１，ｊ，ｋ}＋ｕ_{ｉ，ｊ，ｋ}）／２｝／２Δｚ
以上の説明から明らかなように、数４のうちｘに関する偏微分の値と、ｚに関する偏微分の値とは、基準平面ＰＬの上の速度ベクトルの原データにより求めることができる。ｙに関する偏微分の値についても、ｘ，ｚに関する偏微分の値と同様に、１次元有限差分法のうちの中心差分法で求めた差分で代用する。つまり、次式を用いる。
∂［Ｕ］_{ｉ＋０．５，ｊ，ｋ＋０．５}／∂ｙ＝｛［Ｕ］_{ｉ＋０．５，ｊ＋０．５，ｋ＋０．５}−［Ｕ］_{ｉ＋０．５，ｊ−０．５，ｋ＋０．５}｝／２Δｙ
ただし、Δｙ＝１とする。

ここにおいて、ｙに関する偏微分を差分で表すには、基準平面ＰＬに存在しない速度ベクトルのデータを用いる必要があるから、ｘ，ｚに関する偏微分のように基準平面ＰＬに配置した速度ベクトルの原データを用いるだけでは、ｙに関する偏微分の値を求めることはできない。そこで、ｚ軸を含むすべての断面において流場が幾何学的に同一であることを利用し、以下に説明するように、基準平面ＰＬに配置した速度ベクトルの原データを用いて、座標位置（ｉ＋０．５，ｊ＋０．５，ｋ＋０．５），（ｉ＋０．５，ｊ−０．５，ｋ＋０．５）のそれぞれの速度ベクトル［Ｕ］_{ｉ＋０．５，ｊ＋０．５，ｋ＋０．５}，［Ｕ］_{ｉ＋０．５，ｊ−０．５，ｋ＋０．５}として、基準平面ＰＬの上で座標位置（ｉ＋０．５，ｊ＋０．５，ｋ＋０．５），（ｉ＋０．５，ｊ−０．５，ｋ＋０．５）に対応する対応点における速度ベクトルで代用する。

いま、図６に示すように、ｘ軸上（つまり基準平面ＰＬ上）の点ｑをｚ軸（回転軸３）の周りにφだけ回転させたときに、点ｑが座標位置（ｉ＋０．５，ｊ＋０．５，ｋ＋０．５）に移動するものとする（つまり、点ｑが座標位置（ｉ＋０．５，ｊ＋０．５，ｋ＋０．５）の対応点であるものとする）。流場は幾何学的に同一なので、点ｑの位置における速度ベクトル［Ｕ］_ｑと基準平面ＰＬとの関係は、回転軸３と座標位置（ｉ＋０．５，ｊ＋０．５，ｋ＋０．５）とを含む平面と座標位置（ｉ＋０．５，ｊ＋０．５，ｋ＋０．５）における速度ベクトル［Ｕ］_{ｉ＋０．５，ｊ＋０．５，ｋ＋０．５}との関係と同じとみなせる。また、ｙに関する偏微分の値を求めることが目的であるから、ｘｙ平面について考慮すればよい。つまり、ｘｙ平面において、基準平面ＰＬ（ｘ軸）に対する速度ベクトル［Ｕ］_ｑの角度をψとすれば、回転軸３と座標位置（ｉ＋０．５，ｊ＋０．５，ｋ＋０．５）とを含む平面に対する速度ベクトル［Ｕ］_{ｉ＋０．５，ｊ＋０．５，ｋ＋０．５}の角度もψになる。

上述の性質から、速度ベクトル［Ｕ］_{ｉ＋０．５，ｊ＋０．５，ｋ＋０．５}，［Ｕ］_{ｉ＋０．５，ｊ−０．５，ｋ＋０．５}を求めるには、基準平面ＰＬの上でそれぞれ（ｉ＋０．５，ｊ＋０．５，ｋ＋０．５），（ｉ＋０．５，ｊ−０．５，ｋ＋０．５）に対応する点ｑを求め、点ｑの上の速度ベクトル［Ｕ］_ｑを求めればよいことがわかる。任意の位置の速度ベクトルは、近傍の速度ベクトルを用いた補間演算により求めることができると考えてよいから、点ｑが基準平面ＰＬ（ｘ軸）上において点（ｉ，ｊ，ｋ＋０．５）と点（ｉ＋１，ｊ，ｋ＋０．５）との間に位置していることを利用し、両点（ｉ，ｊ，ｋ＋０．５），（ｉ＋１，ｊ，ｋ＋０．５）における速度ベクトル［Ｕ］_{ｉ，ｊ，ｋ＋０．５}，［Ｕ］_{ｉ＋１，ｊ，ｋ＋０．５}を用いて、点ｑにおける速度ベクトル［Ｕ］_ｑを求めることができる。上述したように、速度ベクトル［Ｕ］_{ｉ，ｊ，ｋ＋０．５}，［Ｕ］_{ｉ＋１，ｊ，ｋ＋０．５}は、既知である４個の原データ［Ｕ］_{ｉ，ｊ，ｋ}、［Ｕ］_{ｉ＋１，ｊ，ｋ}、［Ｕ］_{ｉ，ｊ，ｋ＋１}、［Ｕ］_{ｉ＋１，ｊ，ｋ＋１}を用いて次式のように求めることができる。
［Ｕ］_{ｉ，ｊ，ｋ＋０．５}＝（［Ｕ］_{ｉ，ｊ，ｋ}＋［Ｕ］_{ｉ，ｊ，ｋ＋１}）／２
［Ｕ］_{ｉ＋１，ｊ，ｋ＋０．５}＝（［Ｕ］_{ｉ＋１，ｊ，ｋ}＋［Ｕ］_{ｉ＋１，ｊ，ｋ＋１}）／２

したがって、原点（回転軸３）から点ｑまでの距離をＱとすれば、点ｑにおける速度ベクトル［Ｕ］_ｑは次式で表すことができる。
［Ｕ］_ｑ＝［｛（Ｑ−（ｉ＋０．５）＋（ｉ＋０．５）−ｉ）｝／｛（ｉ＋１）−ｉ｝］［Ｕ］_{ｉ＋１，ｊ，ｋ＋０．５}＋［｛（ｉ＋１）−（ｉ＋０．５）−（Ｑ−（ｉ＋０．５）｝／｛（ｉ＋１）−ｉ｝］［Ｕ］_{ｉ，ｊ，ｋ＋０．５}
＝｛（Ｑ−ｉ）［Ｕ］_{ｉ＋１，ｊ，ｋ＋０．５}＋（（ｉ＋１）−Ｑ）［Ｕ］_{ｉ，ｊ，ｋ＋０．５}｝
＝｛（Ｑ−ｉ）（［Ｕ］_{ｉ＋１，ｊ，ｋ}＋［Ｕ］_{ｉ＋１，ｊ，ｋ＋１}）＋（（ｉ＋１）−Ｑ）（［Ｕ］_{ｉ，ｊ，ｋ}＋［Ｕ］_{ｉ，ｊ，ｋ＋１}）｝／２
ここに、距離Ｑは次式で表される。
Ｑ＝［（ｉ＋０．５）^２＋｛（（ｊ＋０．５）−（ｊ−０．５））／２｝^２］^０．５
さらに、図５に示した関係で理想化すれば、（ｊ＋０．５）−（ｊ−０．５）＝１であるから、これらの関係を適用すれば距離Ｑは次式で表すことができる。
Ｑ＝［｛（ｉ＋０．５｝^２＋（０．５）^２］^０．５

以上の結果を用い、点ｑにおける速度ベクトル［Ｕ］_ｑの絶対値をｒとすると、図６に示す関係から、着目位置の座標での速度ベクトル［Ｕ］_{ｉ＋０．５，ｊ，ｋ＋０．５}の成分ｕ_{ｉ＋０．５，ｊ，ｋ＋０．５}のｙに関する偏微分（差分）の値は次式で与えられる。
∂ｕ_{ｉ＋０．５，ｊ，ｋ＋０．５}／∂ｙ＝（ｕ_{ｉ＋０．５，ｊ＋０．５，ｋ＋０．５}−ｕ_{ｉ＋０．５，ｊ−０．５，ｋ＋０．５}）／２Δｙ
＝ｒ・ｃｏｓ（ψ＋φ）−ｒ・ｃｏｓ（ψ−φ）／２Δｙ
＝−２ｒ・ｓｉｎψｓｉｎφ／２Δｙ
ここで、速度ベクトル［Ｕ］_ｑの各成分をｕ_ｑ，ｖ_ｑ，ｗ_ｑと表記すれば、
ｓｉｎψ＝ｖ_ｑ／ｒ
ｓｉｎφ＝｛（ｊ＋０．５）−（ｊ−０．５）｝／２／Ｑ＝１／２Ｑ
になる。Δｙ＝０．５であるから、上式を次式のように変形することができる。
∂ｕ_{ｉ＋０．５，ｊ，ｋ＋０．５}／∂ｙ＝−２ｒ（ｖ_ｑ／ｒ）｛１／２Ｑ｝
＝−ｖ_ｑ／Ｑ
ここで、［Ｕ］_ｑ＝｛（Ｑ−ｉ）（［Ｕ］_{ｉ＋１，ｊ，ｋ}＋［Ｕ］_{ｉ＋１，ｊ，ｋ＋１}）＋（（ｉ＋１）−Ｑ）（［Ｕ］_{ｉ，ｊ，ｋ}＋［Ｕ］_{ｉ，ｊ，ｋ＋１}）｝／２であるから、ｖ_ｑは次式のようになる。
ｖ_ｑ＝｛（Ｑ−ｉ）（ｖ_{ｉ＋１，ｊ，ｋ}＋ｖ_{ｉ＋１，ｊ，ｋ＋１}）＋（（ｉ＋１）−Ｑ）（ｖ_{ｉ，ｊ，ｋ}＋ｖ_{ｉ，ｊ，ｋ＋１}）｝／２
それゆえ、次式が成立する。
∂ｕ_{ｉ＋０．５，ｊ，ｋ＋０．５}／∂ｙ＝｛（ｉ−Ｑ）（ｖ_{ｉ＋１，ｊ，ｋ}＋ｖ_{ｉ＋１，ｊ，ｋ＋１}）＋（Ｑ−（ｉ＋１））（ｖ_{ｉ，ｊ，ｋ}＋ｖ_{ｉ，ｊ，ｋ＋１}）｝／２Ｑ

同様にして速度ベクトル［Ｕ］_{ｉ＋０．５，ｊ，ｋ＋０．５}の成分ｖ_{ｉ＋０．５，ｊ，ｋ＋０．５}のｙに関する偏微分（差分）の値は次式で与えられる。
∂ｖ_{ｉ＋０．５，ｊ，ｋ＋０．５}／∂ｙ＝｛ｖ_{ｉ＋０．５，ｊ＋０．５，ｋ＋０．５}−ｖ_{ｉ＋０．５，ｊ−０．５，ｋ＋０．５}｝／２Δｙ
＝ｒ・ｓｉｎ（ψ＋φ）−ｒ・ｓｉｎ（ψ−φ）／２Δｙ
＝２ｒ・ｃｏｓψｓｉｎφ／２Δｙ
＝２ｒ（ｕ_ｑ／ｒ）｛１／２Ｑ｝
＝ｕ_ｑ／Ｑ
＝−｛（ｉ−Ｑ）（ｕ_{ｉ＋１，ｊ，ｋ}＋ｕ_{ｉ＋１，ｊ，ｋ＋１}）＋（Ｑ−（ｉ＋１））（ｕ_{ｉ，ｊ，ｋ}＋ｕ_{ｉ，ｊ，ｋ＋１}）｝／２Ｑ
速度ベクトル［Ｕ］_{ｉ＋０．５，ｊ，ｋ＋０．５}の成分ｗ_{ｉ＋０．５，ｊ，ｋ＋０．５}はｚ軸（回転軸３）の周りに回転しても変化しないから、ｙに関する偏微分の値は０になる。つまり、次式が得られる。
∂ｗ_{ｉ＋０．５，ｊ，ｋ＋０．５}／∂ｙ＝０

以上まとめると、速度ベクトル［Ｕ］_{ｉ＋０．５，ｊ，ｋ＋０．５}の成分各ｕ_{ｉ＋０．５，ｊ，ｋ＋０．５}，ｖ_{ｉ＋０．５，ｊ，ｋ＋０．５}，ｗ_{ｉ＋０．５，ｊ，ｋ＋０．５}のｙに関する偏微分の値は、次式のようになる。
∂ｕ_{ｉ＋０．５，ｊ，ｋ＋０．５}／∂ｙ＝｛（ｉ−Ｑ）（ｖ_{ｉ＋１，ｊ，ｋ}＋ｖ_{ｉ＋１，ｊ，ｋ＋１}）＋（Ｑ−ｉ−１）（ｖ_{ｉ，ｊ，ｋ}＋ｖ_{ｉ，ｊ，ｋ＋１}）｝／２Ｑ
∂ｖ_{ｉ＋０．５，ｊ，ｋ＋０．５}／∂ｙ＝−｛（ｉ−Ｑ）（ｕ_{ｉ＋１，ｊ，ｋ}＋ｕ_{ｉ＋１，ｊ，ｋ＋１}）＋（Ｑ−ｉ−１）（ｕ_{ｉ，ｊ，ｋ}＋ｕ_{ｉ，ｊ，ｋ＋１}）｝／２Ｑ
∂ｗ_{ｉ＋０．５，ｊ，ｋ＋０．５}／∂ｙ＝０

以上説明したように、基準平面ＰＬに集約された速度ベクトルの原データのうち着目位置の近傍に位置する４個の原データを用い、着目位置を中心とした基準平面ＰＬ内の正方形を設定してｘ，ｚに関する偏微分の値を求める。また、着目位置からｙ方向に所定距離だけ離れた位置と回転軸３との距離Ｑを求め、基準平面ＰＬの上で回転軸３から距離Ｑに位置する速度ベクトル［Ｕ］_ｑと基準平面ＰＬとの位置関係を用いることによってｙに関する偏微分の値を求めることができる。つまり、基準平面ＰＬの上で着目位置を中心とした４個の速度ベクトルを用いることにより、基準平面ＰＬの上の任意の位置のせん断ひずみ速度を求めることができる。したがって、基準平面ＰＬの上でせん断変形速度テンソル［Ｄ２］を求めることができる。これにより、せん断変形速度テンソルを求めることができるのである。また、同様にして流場を記述する他のテンソルも求めることができる。つまり、せん断ひずみ速度のほか、渦度も求めることができる。さらに、回転軸３を含むすべての断面において流場は幾何学的に同一であるから、基準平面ＰＬにおいて求めた流場は、撹拌槽１の中の任意の位置に適用することが可能であるから、流線も求めることができる。

なお、上述の演算は着目画素を中心とする正方形の頂点の位置の速度ベクトルを用いているが、基準平面ＰＬの上で求められる速度ベクトルが正方形の頂点位置に位置するとは限らない。したがって、実際には上述した演算で用いたように、速度ベクトルについて適宜の補間を行うことにより、正方形の頂点位置に相当する速度ベクトルを求め、この速度ベクトルを用いればよい。

また、上述の例では１次元の有限差分法のうち中心差分を用いたが、１次元の有限差分法のうちの前進差分、リチャードソン法、最小二乗法などを用いることができ、また、２次元の有限差分法である中心差分、最小二乗法などを用いることができる。これらの演算についてｘに関する偏微分の値を差分で代用する場合について概念のみ示しておく。これらの差分を用いる場合も上述した手順に準じる。
前進差分法（１次元）：（∂ｆ／∂ｘ）_{ｉ＋０．５}＝｛ｆ（ｉ＋１，ｊ，ｋ）−ｆ（ｉ，ｊ，ｋ）｝／Δｘ
リチャードソン法（１次元）：（∂ｆ／∂ｘ）_ｉ＝｛ｆ（ｉ−２，ｊ，ｋ）−８ｆ（ｉ−１，ｊ，ｋ）＋８ｆ（ｉ＋１，ｊ，ｋ）−ｆ（ｉ＋２，ｊ，ｋ）｝／１２Δｘ
最小二乗法（１次元）：（∂ｆ／∂ｘ）_ｉ＝｛２ｆ（ｉ＋２，ｊ，ｋ）＋ｆ（ｉ＋１，ｊ，ｋ）−ｆ（ｉ−１，ｊ，ｋ）−ｆ（ｉ−２，ｊ，ｋ）｝／１０Δｘ
中心差分法（２次元）：（∂ｆ／∂ｘ）_ｉ，ｊ＝｛ｆ（ｉ＋１，ｊ＋１，ｋ）−ｆ（ｉ−１，ｊ＋１，ｋ）＋２（ｆ（ｉ＋１，ｊ，ｋ）−ｆ（ｉ−１，ｊ，ｋ））＋ｆ（ｉ＋１，ｊ−＋１，ｋ）−ｆ（ｉ−１，ｊ−１，ｋ）｝／８Δｘ
最小二乗法（２次元）：（∂ｆ／∂ｘ）_ｉ，ｊ＝｛２（ｆ（ｉ＋１，ｊ＋１，ｋ）＋ｆ（ｉ＋１，ｊ−１，ｋ）−ｆ（ｉ−１，ｊ＋１，ｋ）−ｆ（ｉ−１，ｊ−１，ｋ））＋ｆ（ｉ＋１，ｊ，ｋ）−ｆ（ｉ−１，ｊ，ｋ）｝／１０Δｘ

本発明に用いる装置の概略構成図である。 同上における直交座標系と円柱座標系との関係を示す斜視図である。 同上における直交座標系と円柱座標系との関係を示す平面図である。 同上の動作説明図である。 同上における演算の概念を示す図である。 同上においてｙに関する偏微分の値を求める演算の概念を示す図である。

符号の説明

１ 撹拌槽
２ 撹拌翼
３ 回転軸
４ モータ
５ 位置計測装置
６ 流場解析装置
７ 流体
ＰＬ 基準平面

Claims (5)

1. 内側面が撹拌翼の回転軸を中心とする回転体形状に形成された槽であって撹拌翼の回転に伴って撹拌翼の回転軸を中心とする流動を生じる撹拌槽内の流体にトレーサ粒子を混入させ、非接触かつ所定の時間間隔で位置計測を行う位置計測装置によりトレーサ粒子の３次元位置を追跡し、トレーサ粒子の位置の時間変化により撹拌槽内の流場を計測する方法であって、撹拌翼の回転軸の延長方向を一つの座標軸に持つ２次元平面を基準平面として設定し、トレーサ粒子の３次元位置を追跡することにより撹拌槽内の全体で得られる速度ベクトルの原データを回転軸の回りに回転させて原データを基準平面の位置に配置した後、基準平面の上の着目位置の周囲に３次元の微小領域を設定するとともに微小領域内における基準平面の上の複数個の原データを用いて着目位置に関する微分演算を有限差分法で代用することにより着目位置の速度勾配テンソルの成分を求めることを特徴とする粒子追跡法を用いた流場計測方法。
2. 前記微小領域は前記基準平面を含む立方体の単位格子であって、有限差分法で用いる成分のうち基準平面に含まれない速度ベクトルを用いる成分には、前記回転軸の周りに基準平面を回転させたときの対応点の速度ベクトルの成分を適用することを特徴とする請求項１記載の粒子追跡法を用いた流場計測方法。
3. 前記単位格子の一面が前記基準平面であることを特徴とする請求項２記載の粒子追跡法を用いた流場計測方法。
4. 前記着目位置が前記単位格子の中心であって、単位格子の６面のうちの４面の中心の位置の速度ベクトルの成分は前記基準平面の上で前記原データにより求め、単位格子の６面のうちの残りの２面の中心の位置の速度ベクトルの成分は前記回転軸の周りに基準平面を回転させたときの対応点の速度ベクトルの成分により求めることを特徴とする請求項２記載の粒子追跡法を用いた流場計測方法。
5. 前記基準平面の上で求めた速度勾配テンソルの成分を回転軸の周りに回転させることにより撹拌槽内の所望の位置の流場を計測することを特徴とする請求項１ないし請求項４のいずれか１項に記載の粒子追跡法を用いた流場計測方法。

Images