Fluid Governing Equation

流体の支配方程式

Last Update:2009年05月06日

流体の構成式

構成式に関する原理Ⅰ（局所作用の原理）より，点 $\b{X}$ における応力は $\b{X}$ の近傍の領域の物質店の集合の過去から現在時刻に至るまでの運動の履歴によって決定される．特に，ある物質が変形勾配 $\b{F}$ の履歴のみに依存して応力 $\b{T}$ の定まる単純体(simple material)とすると次の式が成り立つ．

${\b{T}}=\cal{F}_{s=0}^{\infty}({\b{F}}(t-s))$

単純物質の構成式においてもし，

$\cal{F}_{s=0}^{\infty}({\b{F}}(t-s))=\cal{F}_{s=0}^{\infty}({\b{F}}(t-s) \cdot {\b{P}})$

とするようななんらかの２階のテンソル $\b{P}$ が存在する場合，この物質は変換 $\b{P}$ に対して対称であると呼ぶ． $\b{P}$ は変形 $\b{F}$ に先立つ物体の変換を表すので，対称変換(symmetric transformation)と呼ばれる．液晶などのように流体に異方性があるような特別な場合を除いて，一般的に流体は密度が変化しなければ力学特性も同一のはずであるから， $\b{P}$ が

$\det{\b{P}}=1$

を満たす任意の変換であるとき，その物質を（単純）流体と定義することができる．なお上式を満たす変換をユニモジュラー変換(unimodular transformation)と呼ぶ．任意の回転は上式を満たすから流体は等方であることが分かる．上式の $\b{P}$ に

${\b{P}}=J\det{\b{F}^{-1}}(t)$

を代入すると流体の構成式は

${\b{T}}=\cal{F}_{s=0}^{\infty}\big({\b{F}}(t-s) \cdot J {\b{F}^{-1}}(t)\big)$

が成り立たつことが必要である．同時に上式が成り立っている場合,任意の $\det{\b{P}}\neq0$ となる変換 $\b{P}$ に対して $\tilde{\b{F}}(t-s)={\b{F}}(t-s)\cdot{\b{P}}$ を上式に代入すると

$\tilde{\b{T}}=\cal{F}_{s=0}^{\infty}\big(\tilde{\b{F}}(t-s) \cdot J \tilde{\b{F}^{-1}}(t)\big)\\ = \cal{F}_{s=0}^{\infty}\big({\b{F}}(t-s)\cdot{\b{P}} \cdot J ({\b{F}}(t)\cdot{\b{P}})^{-1}\\ = \cal{F}_{s=0}^{\infty}\big({\b{F}}(t-s)\cdot{\b{P}} \cdot J {\b{P}^{-1}} \cdot {\b{F}^{-1}}(t)\big)\\ = \cal{F}_{s=0}^{\infty}\big({\b{F}}(t-s) \cdot J {\b{F}^{-1}}(t)\big)\\ = {\b{T}}$

となり，変換 $\b{P}$ に依存しない形で表せるので,流体の構成式として十分満足する．以上から(単純)流体の構成式は次のように書けることが必要十分条件となる．

${\b{T}}=\cal{F}_{s=0}^{\infty}\big({\b{F}}(t-s) \cdot J {\b{F}^{-1}}(t)\big)$

次に上式右辺の( )内が

${\b{F}}(t-s) \cdot J {\b{F}^{-1}}(t) = {\b{F}}(t) \cdot J {\b{F}^{-1}}(t)\\ + \frac{\partial{\b{F}}(t)}{\partial t} \cdot J {\b{F}^{-1}}(t) ds\\ + \frac{1}{2!}\frac{\partial^2{\b{F}}(t)}{\partial^2t} \cdot J{\b{F}^{-1}}(t)ds^2\\ + \cdots\\ = J + J {\b{L}} ds + J \frac{1}{2!} {\b{L}_{(2)}} ds^2 + \cdots$

のようにTaylor展開できるものとする．但し， $\b{L}$ は速度勾配テンソルであり, $\b{L}_{(n)}$ はn階の速度勾配テンソルである．もし上式右辺２項までとって残余項が十分小さくなるなら，構成式は次のように書くことができる．

${\b{T}} = {\b{f}}\big(J,J{\b{L}}\big)\\ = {\b{f}}\big(J,{\b{L}}\big)$

また基準配置からの体積変化率，つまりヤコビアンと密度は上式によって対応付けられるので上式からヤコビアンを除き，代わりに現在の時刻での密度の関数とすることができる．つまり，

${\b{T}}={\b{f}}\big(\rho,{\b{L}}\big)$

となる

ここで，流体の速度が一様，つまり速度勾配が０である場合の応力 $\tilde{\b{T}}$ を考える．速度勾配が０である場合は $\b{L}=0$ であるから，上式より

$\tilde{\b{T}}={\b{f}}(\rho)$

となる．また，物質客観性の原理から $\b{Q}$ を任意の対称テンソルとすると，

$\tilde{\b{T}}^* = {\b{Q}} \tilde{\b{T}} {\b{Q}^{T}} = {\b{f}}(\rho^*) = {\b{f}}(\rho) = \tilde{\b{T}}$

任意の対称テンソル $\b{Q}$ による回転に対して不変なテンソルは単位テンソル $\b{I}$ の定数倍に限られる．よってこの場合の応力 $\tilde{\b{T}}$ は

$\tilde{\b{T}}=-p(\rho){\b{I}}$

とかける．このように流体の速度が一様であるときの応力 $\tilde{\b{T}}$ を静水圧応力(hydrostatic stress)，を静水圧(hydrostatic pressure)と呼ぶ．また上式のように圧力が密度だけから決定されるような流体をバロトロピー流体(Barotoropic fluid)または順圧流体と呼ぶ．また，圧力が密度以外に温度などに依存する流体を斜圧流体(baroclinic fluid)と呼ぶ．また，非圧縮性流体のように圧力が非圧縮性を満たすための拘束力として働く場合，この圧力を不定静水圧という．

応力を決定する関数 $\b{f}$ からこの静水圧応力 $\tilde{\b{T}}$ を除いた関数を新たに $\b{f}$ と定義すると，

${\b{T}}=-p(\rho){\b{I}}+{\b{f}}(\rho,{\b{L}})$

となる．但し， $\b{f}$ は流体の速度が一様，つまり速度勾配が０であるときの応力を応力 $\b{T}$ から除いたものなので，速度勾配が0である場合に $\b{f}=0$ を満たさなければならない．つまり,

${\b{f}}(\rho,0)=0$

である．また常に $\b{f}=0$ が成り立つ場合，つまり

${\b{T}}=-p(\rho){\b{I}}$

である場合，この構成式は理想流体を表す．流体の基本的物理特性は並行にあるとき剪断応力を保持できないということであるが，このような構成式で表される流体はそれが運動中にこの性質を保てるという意味で理想的なのである．理想流体において剪断応力が存在しないということは，単純剪断運動において隣り合う層が抵抗なくお互いに滑りあうことができることを意味している．粘性流体においてはこの種の相対運動は摩擦抵抗をうけ，剪断力を発生させるので理想流体は非粘性流体(inviscid fluid)といっても差し支えない.

非理想流体の場合は $\b{f}\ne0$ であり， $\b{f}$ は理想流体からのずれの応力を表すので $\b{f}$ は過剰応力と呼ばれる．一般的に過剰応力は粘性による応力を表している．

速度勾配テンソル $\b{L}$ の対称成分を $\b{D}$ ，反対称成分を $\b{W}$ とすると，

${\b{L}} = {\b{D}} + {\b{W}}$

となる． $\b{D}$ は変形速度テンソル(deformation rate tensor) $\b{W}$ はスピンテンソル(spin tensor)と呼ばれる． $\b{W}$ は剛体回転を表しており，応力に関係しないつまり，構成式は変形速度テンソル $\b{D}$ のみで表すことができると考えられる．実際のところ，物質客観性の原理より，関数 $\b{f}$ は任意の直行テンソル $\b{Q}$ とその時間微分で任意の反対称テンソルとなる $\dot{\b{Q}}$ に対して

${\b{T}^*}={\b{Q}}{\b{T}}{\b{Q}^T}\\ = -p(\rho){\b{I}}+{\b{Q}}{\b{f}}\big(\rho,{\b{L}}\big){\b{Q}^T}\\ = -p(\rho){\b{I}}+{\b{f}}\big(\rho,{\b{L}^*}\big)\\ = -p(\rho){\b{I}}+{\b{f}}\big(\rho,{\b{Q}}{\b{L}}{\b{Q}^T}+\dot{\b{Q}}{\b{Q}^T})$

を満たさなければならない． $\b{Q}$ ， $\dot{\b{Q}}$ は任意であったことから，ここで $\b{Q}={\b{I}$ ， $\dot{\b{Q}}$ を $\b{L}$ の反対称成分にマイナスしたもの $\dot{\b{Q}}=-({\b{L}})_a=-{\b{W}}=-\frac{1}{2}({\b{L}}-{\b{L}^T})$ とすると

${\b{f}}\big(\rho,{\b{L}}\big)={\b{f}}\big(\rho,{\b{L}}-({\b{L}})_a\big)={\b{f}}\big(\rho,({\b{L}})_s\big)={\b{f}}\big(\rho,{\b{D}}\big)$

であるから，

${\b{T}}=-p({\rho}){\b{I}}+{\b{f}}\big(\rho,{\b{D}}\big)$

となり，過剰応力 $\b{f}$ は歪速度テンソル $\b{D}$ の関数となる

速度勾配テンソル $\b{L}$ の反対称成分 $\b{W}$ が出てきたが，一般に反対称テンソル $\b{B}$ は軸性ベクトル(axial vector)

${\b{\omega}}=-\frac{1}{2}{\epsilon_{ijk}}{B_{ij}}{\b{e}_k}$

を持ち任意のベクトル $\b{a}$ に対して $\b{B}\cdot{\b{a}}=\omega\times{\b{a}$ となる．スピンテンソル $\b{W}$ に対して軸性ベクトル $\b{\omega}$ は

${\b{\omega}} = -\frac{1}{2}{\epsilon_{ijk}}{W_{jk}}{\b{e}_i} = \frac{1}{2}{\epsilon_{ijk}}{\frac{\partial\b{v}_k}{\partial\b{x}_j}}{\b{e}_i} = \frac{1}{2}{\nabla_x \times {\b{v}}}$

で表され，流体が剛体回転する角速度ベクトルを表している．上式より角速度ベクトル $\b{\omega}$ は $\b{\Omega} = \nabla_x \times {\b{v}}$ で定義される渦度ベクトル(vorticity)の半分になる．

非圧縮性の式

流体に非圧縮性の条件を仮定すれば，物質点 $\b{X}$ の質量密度は時間に依存せず，一定であるから以下の条件を満たすことになる．

$\left.\frac{\partial\rho}{\partial t}\right|_{\b{X}} = 0$

したがってこれを連続の式に適応すれば，Lagrange表示，Eular表示に対してそれぞれ非圧縮性の連続の式は以下のようになる．

$\nabla_x \cdot {\b{v}} = 0$

$\left.\frac{\partial\rho}{\partial t}\right|_{\b{x}}+{\b{v}}\cdot(\nabla_x\rho)=0$

さらに単一流体，もしくは空間的に流体の密度が一定の場合

$\rho = Const$

であるから，非圧縮性の条件は常に満たされ，連続の式はLagrange表示形式のみ意味を持つ．すなわち，連続の式のEular表示は常に満たされる．以上より密度一定流体の非圧縮性の式は以下のようになる．

$\nabla_x \cdot {\b{v}} = 0$

非圧縮性の式
$\nabla_x \cdot {\b{v}} = 0$

流体の圧力

物質の変形に制約条件が存在する場合がある．非圧縮性の条件などに代表されるこのような制約は内部束縛(internal constraint)と呼ばれる．物体に内部束縛が存在する場合には構成式に関する原理Ⅰ（応力決定の原理）は次のように変形される必要がある．

${\b{T}}=\cal{F}_{s=0}^{\infty}({\b{F}}(t-s))+\cal{T}$

この応力 $\cal{T}$ は内部仕事に無関係であり，物質点の運動，すなわち $\b{F}$ の履歴からは決定することができない応力なので非決定応力(interminate stress)と呼ばれる．非決定応力は仕事をしないので

$\cal{T} : {\b{D}} = \tr( \cal{T} \cdot {\b{D}}) = 0$

が成り立つ．

代表的な内部拘束はスカラー値テンソル関数

$\gamma({\b{F}}) = 0$

のように表される．たとえばその代表例としては非圧縮の拘束条件

$\gamma({\b{F}}) = \det{\b{F}} - 1 = J-1 = 0$

があげられる．

内部仕事も物質客観性の原理を満たさなければならないことから，

$\gamma^*({\b{F}}) = \gamma({\b{F}}) = \gamma({\b{Q}} \cdot {\b{F}})$

を満たさなければならない．ここで $\b{Q}={\b{R}^T}$ とおくと

$\gamma({\b{F}}) = \gamma({\b{R}^T} \cdot {\b{R}} \cdot {\b{U}}) = \gamma({\b{U}})$

となり $\gamma$ は右ストレッチテンソル $\b{U}$ のみの関数となることが必要である．一方で右ストレッチテンソルは基準枠の回転に対して不変つまり， $\b{U}^*=\b{U}$ であるから

$\gamma^*({\b{U}})=\gamma({\b{U}}^*)$

となり十分客観性を満たしていることがわかる．よって拘束条件 $\gamma$ は右ストレッチテンソル $\b{U}$ のみの関数でなければならない．

ここで非圧縮性の拘束条件にたいしてどのような非決定応力が発生するかを調べる． $\gamma$ は０で定数であったので $\lambda$ を時間微分したものも０であるはずである．よって

$\dot{\gamma} = (\det{\b{F}}-1)^{\dot{}} = \dot{J} = J\tr{\b{D}} = \tr{\b{D}} = 0$

ここで拘束条件から J=1 であることを用いた．上式から

$\tr{\b{D}}={\b{I}}:{\b{D}}=\tr({\b{I}}\cdot{\b{D}})=0$

となる．よって上式から非圧縮性の制約条件が課された場合に単位テンソルの定数倍で表される応力は仕事をしないということがわかる．よって非決定応力 $\cal{T}$ は未知数を用いて

$\cal{T}=-p{\b{I}}$

のように表すことができる．ここで導入したは明らかに不定静水圧である．

Newton流体の構成式

過剰応力が密度に依存しないような流体をStorks流体(Storkesian fluid)と呼び，その構成式は

${\b{T}} = -p{\b{I}} +{\b{f}}({\b{D}})$

により与えられる．は圧力, $\b{I}$ は単位テンソルである．テンソル $\b{D}$ が客観性のあるテンソルであることから，物質客観性の原理により関数 $\b{f}$ は

${\b{Q}} \cdot {\b{f}}({\b{D}}) \cdot {\b{Q}^T} = {\b{f}}( {\b{Q}} \cdot {\b{D}} \cdot {\b{Q}^T} )$

を満たさなければならない．つまり $\b{f}$ は等方テンソル関数となるから，表示定理により， $\phi_i$ を $\b{D}$ の主不変量のスカラー関数として

${\b{f}}({\b{D}}) = \phi_0{\b{I}} + \phi_1 {\b{D}} + \phi_2 {\b{D}}^2$

と表すことができる．上式を過剰応力とする流体はReiner-Rivlin流体とも呼ばれる．上式第３項は外力条件が急激に変化する場合に影響を持つことが知られているが，一般的に強い非線形性を示すような流体においてもほぼ無視できるほど小さいことが分かっている．

過剰応力 $\b{f}(\b{D})$ が $\b{D}$ に関して線形同次であるような流体をNewton流体という．上式において $\phi_1=\lambda(tr{\b{D}})$ , $\phi_2=2\mu$ , $\phi_3=0$ とおくと

${\b{f}}({\b{D}})=\lambda(tr{\b{D}}){\b{I}} + 2\mu{\b{D}}$

となるのでNewton流体の構成式は以下のように表すことができる．

${\b{T}} = -p {\b{I}} + \lambda(tr{\b{D}}){\b{I}} + 2\mu{\b{D}}\\= -p {\b{I}} + \lambda(\nabla_x \cdot {\b{v}}){\b{I}} + 2\mu{\b{D}}$

ここで $\mu$ はせん断粘性率(shear viscosity), $\lambda$ は第二粘性率(second viscosity)と呼ばれる．第二粘性率 $\lambda$ は線形弾性体におけるLameの第一定数のようなものである。

第二粘性率については以下のリンクが詳しい

International College of FEM : FEM : Fluid Dynamics : Second Coefficient of Viscosity: http://www.fem.gr.jp/fem/fluid/2ndviscosity/2ndviscosity.html

Newton流体の構成式
${\b{T}} = -p {\b{I}} + \lambda(\nabla_x \cdot {\b{v}}){\b{I}} + 2\mu{\b{D}}$

上式を $\b{D}$ の偏差成分 $\tilde{\b{D}$

$\tilde{\b{D}} = {\b{D}} - \frac{1}{3}(tr{\b{D}}){\b{I}}\\ = {\b{D}} - \frac{1}{3}(\nabla_x \cdot {\b{v}}){\b{I}}$

についてみると

${\b{T}} = -p{\b{I}} + (\lambda+\frac{2}{3}\mu)(tr{\b{D}}){\b{I}}+2\mu\tilde{\b{D}}$

となる．ここで現れる $\lambda+\frac{2}{3}\mu$ を体積粘性率(bulk viscosity) $\kappa$ と呼ぶ．体積粘性率は流体の体積変化に伴なう粘性応力の目安である．上式は $\kappa$ を用いて

${\b{T}} = -p{\b{I}} + \kappa(\nabla_x \cdot {\b{v}}){\b{I}}+2\mu\tilde{\b{D}}$

と現わすことができる．多くの流体では $\kappa$ は非常に小さいことが知られている.そこで， $\kappa=0$ とするStokesの近似を用いると圧縮性流体の構成式は次のようになる

${\b{T}} = -p{\b{I}} + 2\mu\tilde{\b{D}}\\ = -p{\b{I}}+2\mu( {\b{D}} - \frac{1}{3}(\nabla_x \cdot {\b{v}}){\b{I}} )$

これは体積粘性係数を無視した圧縮性Newton流体の構成式である．さらに流体に非圧縮性を仮定すると

$\nabla_x \cdot {\b{v}} = 0$

だから，非圧縮性のNewton流体の構成式は

${\b{T}} = -p{\b{I}} + 2\mu{\b{D}}$

となる．

非圧縮性Newton流体の構成式
${\b{T}} = -p{\b{I}} + 2\mu{\b{D}}$

Navier-Storks方程式

Newton流体の構成式は以下のとおり

${\b{T}}=-p{\b{I}}+\lambda(tr{\b{D}}){\b{I}}+ 2 \mu {\b{D}}$

$\mu$ は粘性係数である．これをオイラー座標系で記述されたCauchyの第一運動法則

$\rho\left. \frac{\partial \b{v}}{\partial t} \right|_{\b{x}}+\rho({\b{v}} \otimes \nabla ) \cdot {\b{v}}=\rho {\b{g}}+ \nabla_x \cdot {\b{T}}$

に代入して

$\rho\left. \frac{\partial \b{v}}{\partial t} \right|_{\b{x}}+\rho({\b{v}} \otimes \nabla_x ) \cdot {\b{v}}=\rho {\b{g}}+ \nabla_x \cdot (-p{\b{I}}+\lambda(tr{\b{D}}){\b{I}}+2\mu{\b{D}})$

ここで

$\nabla_x \cdot (-p{\b{I}}) = -\frac{\partial}{\partial x_k}{\b{e}}_k\cdot( p \delta_{ji}{\b{e}}_j \otimes {\b{e}}_i ) = -\frac{\partial p}{\partial x_k}({\b{e}}_k \cdot {\b{e}}_j )\delta_{ji}{\b{e}}_i \\ = -\frac{\partial p}{\partial x_k}\delta_{kj}\delta_{ij}{\b{e}}_i = -\frac{\partial p}{\partial x_i}{\b{e}}_i = -\nabla_x p$

を代入すると

$\rho\left.\frac{\partial {\b{v}}}{\partial t} \right|_x + \rho({\b{v}} \otimes \nabla_x ) \cdot {\b{v}} = \rho {\b{g}}+-\nabla_x p+\nabla \cdot (\lambda(tr{\b{D}}){\b{I}}+2\mu{\b{D}})$

これを応力発散形式のNavier-Storks方程式と呼ぶ．

ここで粘性係数 $\mu$ ,第二粘性係数 $\lambda$ が流体中で一定とした場合について，この応力発散形式のNavier-Storks方程式を変形する．

$\nabla_x \cdot ( \nabla \otimes {\b{v}} ) = \frac{\partial}{\partial x_k}{\b{e}}_k \cdot ( \frac{\partial v_i}{\partial x_l} {\b{e}}_l \otimes {\b{e}}_i ) = \frac{\partial}{\partial x_k}\frac{\partial v_i}{\partial x_l}( {\b{e}}_l \cdot {\b{e}}_i ){\b{e}}_i \\ = \frac{\partial}{\partial x_k}\frac{\partial v_i}{\partial x_l}( \delta_{kl} ){\b{e}}_i = \frac{\partial}{\partial x_k}\frac{\partial v_i}{\partial x_k}{\b{e}}_i \\ = \nabla_x^2 {\b{v}}$

$\nabla \cdot ( {\b{v}} \otimes \nabla ) = \frac{\partial}{\partial x_k}{\b{e}}_k \cdot ( \frac{\partial v_l}{\partial x_i} {\b{e}}_l \otimes {\b{e}}_i )\\ = \frac{\partial}{\partial x_k}\frac{\partial v_k}{\partial x_i}{\b{e}}_i = \frac{\partial}{\partial x_i} \cdot \frac{\partial v_k}{\partial x_k}{\b{e}}_i \\= \frac{\partial}{\partial x_i}{\b{e}}_i( \frac{\partial}{\partial x_k}{\b{e}}_k \cdot {\b{v}}_l ) \\ = \nabla( \nabla \cdot {\b{v}} )$

よって， $\b{D}$ の発散は以下のように変形される．

$\nabla_x \cdot {\b{D}} = \nabla_x \cdot\frac{1}{2}({\b{v}} \otimes \nabla_x+\nabla_x \otimes {\b{v}}) \\ = \frac{1}{2}( \nabla^2 {\b{v}}+\nabla ( \nabla_x \cdot {\b{v}} ) )$

さらに，

$\nabla_x \cdot (tr{\b{D}}){\b{I}} = \nabla_x ( \nabla_x \cdot {\b{v}} )$

これを粘性係数，第二粘性係数が流体内で一定とした応力発散形式のNavier-Stokes方程式に代入するとNavier-Stokes方程式は

$\rho\left.\frac{\partial {\b{v}}}{\partial t} \right|_x + \rho({\b{v}} \otimes \nabla_x ) \cdot {\b{v}} = \rho {\b{g}}+ -\nabla_x p+\mu \nabla^2 {\b{v}}+(\mu+\lambda) \nabla_x ( \nabla_x \cdot {\b{v}} )$

$\kappa=\lambda+\frac{2}{3}\mu$ を代入して

$\rho\left.\frac{\partial {\b{v}}}{\partial t} \right|_x + \rho({\b{v}} \otimes \nabla_x ) \cdot {\b{v}} = \rho {\b{g}}+ -\nabla_x p+\mu \nabla^2 {\b{v}}+(\kappa+\frac{1}{3}\mu) \nabla_x ( \nabla_x \cdot {\b{v}} )$

となる．ここでStokesの仮定 $\kappa=0$ が成り立つ場合，Navier-Stokes方程式は次のようにかける．

$\rho\left.\frac{\partial {\b{v}}}{\partial t} \right|_x + \rho({\b{v}} \otimes \nabla_x ) \cdot {\b{v}} = \rho {\b{g}} + -\nabla_x p + \mu \nabla^2 {\b{v}}+\frac{1}{3}\mu \nabla_x ( \nabla_x \cdot {\b{v}} )$

さらに流体に非圧縮性が成り立つ場合 $\nabla_x \cdot {\bf{v}}=0$ が成り立つので

$\rho\left.\frac{\partial {\b{v}}}{\partial t} \right|_x + \rho({\b{v}} \otimes \nabla_x ) \cdot {\b{v}} = \rho {\b{g}}+ -\nabla_x p+\mu \nabla^2 {\b{v}}$

これはラプラシアン形式の非圧縮性Navier-Storks方程式と呼ばれる．

非圧縮性Navier-Storkes方程式
$\rho\left.\frac{\partial {\b{v}}}{\partial t} \right\|_x + \rho({\b{v}} \otimes \nabla_x ) \cdot {\b{v}} = \rho {\b{g}}+ -\nabla_x p+\mu \nabla^2 {\b{v}}$

参考文献

Books

非線形有限要素法のためのテンソル解析の基礎	久田俊明　著
「連続体力学―簡明な理論と例題」	Ｐ．チャドウィック　著

Last Update:2009年05月06日

Made by Nobuyuki UMETANI 梅谷　信行

n.umetani@gmail.com