October 15, 2025 Last updated: October 15, 2025
流体力学の基礎方程式を導出する場合、弾性力学と対応するような形で理解しておくと見通しがよい。どちらも同じ連続体の運動量の式（運動方程式）からスタートして、
ラグランジュ（物質）微分を用いるか、オイラー（空間）微分を用いるか
物質の物理的特性（応力と変位・速度の関係）を表す構成則をどう仮定するか
によって、異なる物理考察が生まれてくる。このことを意識しつつナビエ・ストークス方程式（Navier-Stokes Equations）を導出するところまでを見ていこう。
今回の記事では、レイノルズの輸送定理(1)と連続の式(2)は既知として利用するので、詳細については以前の記事（レイノルズの輸送定理と連続の式 ）を参照してほしい。
D f D t = ∫ V ( t ) ( ∂ g ∂ t + ∇ ⋅ ( g v ) ) d V where f = ∫ V ( t ) g ( x 1 , x 2 , x 3 , t ) d V \begin{align}
&\frac{Df}{Dt} = \int_{V(t)} \left( \frac{\partial g}{\partial t} + \nabla \cdot (g \bm{v}) \right) dV \\
&\text{where} \quad f = \int_{V(t)} g(x_1, x_2, x_3, t) dV \notag
\end{align} D t D f = ∫ V ( t ) ( ∂ t ∂ g + ∇ ⋅ ( g v ) ) d V where f = ∫ V ( t ) g ( x 1 , x 2 , x 3 , t ) d V ∂ ρ ∂ t + ∇ ⋅ ( ρ v ) = 0 \begin{equation}
\frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \bm{v}) = 0
\end{equation} ∂ t ∂ ρ + ∇ ⋅ ( ρ v ) = 0 コーシーの第一運動法則
まず、物質に作用する力として物体力ρ f \rho \boldsymbol{f} ρ f t \boldsymbol{t} t
D D t ( ∫ ρ v d V ) = ∫ ρ f d V + ∫ t d s \begin{equation}
% \label{eq:Hisada1992_4.13}
\frac{D}{Dt} \left( \int \rho \boldsymbol{v} dV \right) = \int \rho \boldsymbol{f} dV + \int \boldsymbol{t} ds
\end{equation} D t D ( ∫ ρ v d V ) = ∫ ρ f d V + ∫ t d s 運動量のラグランジュ微分を、Reynoldsの輸送定理(1)および質量保存の式(2)を用いて以下のように書き直す。
D D t ( ∫ ρ v i d V ) = ∫ ( ρ ∂ v i ∂ t + v i ∂ ρ ∂ t + ∇ ⋅ ( ρ v i v ) ) d V = ∫ ( ρ ∂ v i ∂ t + v i ∂ ρ ∂ t + v i ∇ ⋅ ( ρ v ) + ρ v ⋅ ∇ v i ) d V = ∫ ( ρ ∂ v i ∂ t + ρ v ⋅ ∇ v i ) d V \begin{align}
&\frac{D}{Dt} \left( \int \rho v_i dV \right) = \int \left( \rho \frac{\partial v_i}{\partial t} + v_i \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho v_i \boldsymbol{v}) \right) dV \notag \\
&= \int \left( \rho \frac{\partial v_i}{\partial t} + v_i \frac{\partial \rho}{\partial t} + v_i \nabla \cdot (\rho \bm{v}) + \rho \bm{v} \cdot \nabla v_i \right) dV \notag \\
&= \int \left( \rho \frac{\partial v_i}{\partial t} + \rho \bm{v} \cdot \nabla v_i \right) dV
\end{align} D t D ( ∫ ρ v i d V ) = ∫ ( ρ ∂ t ∂ v i + v i ∂ t ∂ ρ + ∇ ⋅ ( ρ v i v ) ) d V = ∫ ( ρ ∂ t ∂ v i + v i ∂ t ∂ ρ + v i ∇ ⋅ ( ρ v ) + ρ v ⋅ ∇ v i ) d V = ∫ ( ρ ∂ t ∂ v i + ρ v ⋅ ∇ v i ) d V 被積分関数の中身がオイラー微分の形になっていることに注意すると、これをラグランジュ微分に戻すことができる。
D D t ( ∫ ρ v i d V ) = ∫ ρ D v i D t d V \begin{align}
\frac{D}{Dt} \left( \int \rho v_i dV \right) = \int \rho \frac{Dv_i}{Dt} dV
\end{align} D t D ( ∫ ρ v i d V ) = ∫ ρ D t D v i d V D v D t = a \frac{D\bm{v}}{Dt} = \bm{a} D t D v = a t \boldsymbol{t} t T \boldsymbol{T} T t = T ⋅ n \boldsymbol{t} = \boldsymbol{T} \cdot \boldsymbol{n} t = T ⋅ n
∫ ρ ( a − f ) d V = ∫ d i v T d V \begin{equation}
\int \rho (\boldsymbol{a} - \boldsymbol{f}) dV = \int \mathrm{div}~ \boldsymbol{T} dV
\end{equation} ∫ ρ ( a − f ) d V = ∫ div T d V これが物体の任意の一部分について成り立つので、Cauchy’s First Law of Motionが(7)のように得られる。
この式は、ラグランジュ微分を用いた弾性力学のための運動方程式であり、注目する物質の変形を追いかけていく形で表されている。
ρ a = d i v T + ρ f \begin{equation}
% \label{Hisada1992_4.28}
\rho \boldsymbol{a} = \mathrm{div}~ \boldsymbol{T} + \rho \boldsymbol{f}
\end{equation} ρ a = div T + ρ f 一方で、先ほど出てきたオイラー微分の形を用いると、ある領域に注目した、流体力学における運動方程式が得られる。
ρ ∂ v ∂ t + ρ v ⋅ ∇ v = d i v T + ρ f \begin{equation}
% \label{Hisada1992_4.29}
\rho \frac{\partial \boldsymbol{v}}{\partial t} + \rho \boldsymbol{v} \cdot \nabla \boldsymbol{v} = \mathrm{div}~ \boldsymbol{T} + \rho \boldsymbol{f}
\end{equation} ρ ∂ t ∂ v + ρ v ⋅ ∇ v = div T + ρ f 構成則
次に重要なのが物質の物理特性を表す構成則をどう仮定するかだ。これによって応力がどう発生するか、つまり具体的なT T T T \boldsymbol{T} T
T = T i j e i ⊗ e j \begin{equation}
\boldsymbol{T} = T_{ij} \boldsymbol{e}_i \otimes \boldsymbol{e}_j
\end{equation} T = T ij e i ⊗ e j 弾性力学では、物質の変形や変形履歴に応じて物質内部に発生する応力を応力テンソルとして表し、これは物理現象をもとに仮定される。最も単純なものの例として、線形等方弾性体の構成式のHookeの法則（要するに力に比例して伸びたり縮んだりする）がある。
流体を扱う場合、流体の速度勾配に応じた応力が発生すると考える。速度勾配テンソル（velocity gradient tensor）L \boldsymbol{L} L
d v = L ⋅ d x \begin{equation}
d\boldsymbol{v} = \boldsymbol{L} \cdot d\boldsymbol{x}
\end{equation} d v = L ⋅ d x これより、具体的には以下のように表される。
L = v ⊗ ∇ = v i e i ⊗ ∂ ∂ x j e j = ∂ v i ∂ x j e i ⊗ e j = L i j e i ⊗ e j \begin{equation}
\boldsymbol{L} = \boldsymbol{v} \otimes \nabla = v_i \boldsymbol{e}_i \otimes \frac{\partial}{\partial x_j} \boldsymbol{e}_j = \frac{\partial v_i}{\partial x_j} \boldsymbol{e}_i \otimes \boldsymbol{e}_j = L_{ij} \boldsymbol{e}_i \otimes \boldsymbol{e}_j
\end{equation} L = v ⊗ ∇ = v i e i ⊗ ∂ x j ∂ e j = ∂ x j ∂ v i e i ⊗ e j = L ij e i ⊗ e j これを対称成分の変形速度テンソル（deformation rate tensor or strain rate tensor）D \boldsymbol{D} D W \boldsymbol{W} W L = D + W \boldsymbol{L} = \boldsymbol{D} + \boldsymbol{W} L = D + W
D = 1 2 ( ∂ v i ∂ x j + ∂ v j ∂ x i ) e i ⊗ e j , W = 1 2 ( ∂ v i ∂ x j − ∂ v j ∂ x i ) e i ⊗ e j \begin{equation}
\boldsymbol{D} = \frac{1}{2} \left( \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i} \right) \boldsymbol{e}_i \otimes \boldsymbol{e}_j, \quad
\boldsymbol{W} = \frac{1}{2} \left( \frac{\partial v_i}{\partial x_j} - \frac{\partial v_j}{\partial x_i} \right) \boldsymbol{e}_i \otimes \boldsymbol{e}_j
\end{equation} D = 2 1 ( ∂ x j ∂ v i + ∂ x i ∂ v j ) e i ⊗ e j , W = 2 1 ( ∂ x j ∂ v i − ∂ x i ∂ v j ) e i ⊗ e j スピンテンソルはある軸周りでの回転を表すもので、このような渦流れは応力を生じないと考える。これよりStokes流体の構成式として以下のように仮定する。
T = p I + f ( D ) \begin{equation}
% \label{Hisada1992_5.119}
\boldsymbol{T} = p \boldsymbol{I} + \boldsymbol{f} (\boldsymbol{D})
\end{equation} T = p I + f ( D ) p p p f ( D ) \boldsymbol{f} (\boldsymbol{D}) f ( D ) D \boldsymbol{D} D f ( O ) = 0 \boldsymbol{f} (\boldsymbol{O})= 0 f ( O ) = 0 f ( D ) \boldsymbol{f} (\boldsymbol{D}) f ( D ) D \boldsymbol{D} D
T = { − p + ( κ − 2 3 μ ) t r D } I + 2 μ D \begin{equation}
% \label{Hisada1992_5.122}
\boldsymbol{T} = \left\{ -p + \left( \kappa - \frac{2}{3} \mu \right) \mathrm{tr}~\boldsymbol{D} \right\} \boldsymbol{I} + 2 \mu \boldsymbol{D}
\end{equation} T = { − p + ( κ − 3 2 μ ) tr D } I + 2 μ D ここでκ \kappa κ μ \mu μ p p p κ ( t r D ) \kappa (\mathrm{tr}~\boldsymbol{D}) κ ( tr D ) − 2 3 μ ( t r D ) I - \frac{2}{3} \mu (\mathrm{tr}~\boldsymbol{D}) \boldsymbol{I} − 3 2 μ ( tr D ) I 2 μ D 2\mu \boldsymbol{D} 2 μ D
d i v v = ∂ v i ∂ x i = t r L = t r D = 0 \begin{equation}
\mathrm{div}\ \boldsymbol{v} = \frac{\partial v_i}{\partial x_i}= \mathrm{tr}~\boldsymbol{L} = \mathrm{tr}~\boldsymbol{D} = 0
\end{equation} div v = ∂ x i ∂ v i = tr L = tr D = 0 T = − p I + 2 μ D \begin{equation}
% \label{Hisada1992_5.124}
\boldsymbol{T} = -p \boldsymbol{I} + 2 \mu \boldsymbol{D}
\end{equation} T = − p I + 2 μ D ナビエ・ストークス方程式
オイラー微分を使った運動方程式(8)に、圧縮性を考慮したNewton流体の構成式(14)適用すると、Navier-Stokes方程式が得られる。
ρ ∂ v ∂ t + ρ v ⋅ ∇ v = − ∇ p + μ Δ v + ( κ + 1 3 μ ) ∇ ( ∇ ⋅ v ) + ρ f \begin{equation}
\rho \frac{\partial \boldsymbol{v}}{\partial t} + \rho \boldsymbol{v} \cdot \nabla \boldsymbol{v} = - \nabla p + \mu \Delta \boldsymbol{v} + \left( \kappa + \frac{1}{3} \mu \right) \nabla (\nabla \cdot \boldsymbol{v}) + \rho \boldsymbol{f}
\end{equation} ρ ∂ t ∂ v + ρ v ⋅ ∇ v = − ∇ p + μ Δ v + ( κ + 3 1 μ ) ∇ ( ∇ ⋅ v ) + ρ f もちろん以下のように書いてもよい。
ρ [ ∂ v ∂ t + ( v ⋅ g r a d v ) ] = − g r a d p + μ Δ v + ( κ + 1 3 μ ) g r a d d i v v + ρ f \begin{equation}
% \label{eq:LandauLifshitzVol6_15.6}
\rho \left[ \frac{\partial \boldsymbol{v}}{\partial t} + (\boldsymbol{v} \cdot \mathrm{grad}\ \boldsymbol{v}) \right]
= - \mathrm{grad}\ p + \mu \Delta \boldsymbol{v} + \left(\kappa + \frac{1}{3} \mu \right) \mathrm{grad}\ \mathrm{div}\ \boldsymbol{v} + \rho \boldsymbol{f}
\end{equation} ρ [ ∂ t ∂ v + ( v ⋅ grad v ) ] = − grad p + μ Δ v + ( κ + 3 1 μ ) grad div v + ρ f ただし、応力テンソルのdivergenceは次のように整理した。
∇ ⋅ T = e l ∂ ∂ x l ⋅ [ { − p + ( κ − 2 3 μ ) t r D } I + 2 μ D ] = − e l ∂ p ∂ x l ⋅ δ m n ( e m ⊗ e n ) + e l ∂ ∂ x l ⋅ δ m n ( κ − 2 3 μ ) D k k ( e m ⊗ e n ) + e l ∂ ∂ x l ⋅ 2 μ D i j ( e i ⊗ e j ) = − δ m n ∂ p ∂ x l ( e l ⋅ e m ) e n + δ m n ( κ − 2 3 μ ) ∂ D k k ∂ x l ( e l ⋅ e m ) e n + 2 μ ∂ D i j ∂ x l ( e l ⋅ e i ) e j = − ∂ p ∂ x n e n + ( κ − 2 3 μ ) ∂ D k k ∂ x n e n + 2 μ ∂ D i j ∂ x i e j = − ∂ p ∂ x n e n + ( κ − 2 3 μ ) ∂ ∂ x n ( ∂ v k ∂ x k ) e n + μ ∂ ∂ x i ( ∂ v i ∂ x j + ∂ v j ∂ x i ) e j = − ∂ p ∂ x n e n + ( κ + 1 3 μ ) ∂ ∂ x n ( ∂ v k ∂ x k ) e n + μ ∂ 2 v j ∂ x i 2 e j = − ∇ p + μ Δ v + ( κ + 1 3 μ ) ∇ ( ∇ ⋅ v ) \begin{align*}
&\nabla \cdot \boldsymbol{T} = \boldsymbol{e}_l \frac{\partial}{\partial x_l} \cdot \left[ \left\{ -p + \left( \kappa - \frac{2}{3} \mu \right) \mathrm{tr}~\boldsymbol{D} \right\} \boldsymbol{I} + 2 \mu \boldsymbol{D} \right] \\
&= - \boldsymbol{e}_l \frac{\partial p}{\partial x_l} \cdot \delta_{mn} (\boldsymbol{e}_m \otimes \boldsymbol{e}_n) + \boldsymbol{e}_l \frac{\partial}{\partial x_l} \cdot \delta_{mn} \left( \kappa - \frac{2}{3} \mu \right) D_{kk} (\boldsymbol{e}_m \otimes \boldsymbol{e}_n) \\
&\hspace{12pt}+ \boldsymbol{e}_l \frac{\partial}{\partial x_l} \cdot 2\mu D_{ij} (\boldsymbol{e}_i \otimes \boldsymbol{e}_j) \\
&= - \delta_{mn} \frac{\partial p}{\partial x_l} (\boldsymbol{e}_l \cdot \boldsymbol{e}_m) \boldsymbol{e}_n + \delta_{mn} \left( \kappa - \frac{2}{3} \mu \right) \frac{\partial D_{kk}}{\partial x_l} (\boldsymbol{e}_l \cdot \boldsymbol{e}_m) \boldsymbol{e}_n + 2 \mu \frac{\partial D_{ij}}{\partial x_l} (\boldsymbol{e}_l \cdot \boldsymbol{e}_i) \boldsymbol{e}_j \\
&= - \frac{\partial p}{\partial x_n} \boldsymbol{e}_n + \left( \kappa - \frac{2}{3} \mu \right) \frac{\partial D_{kk}}{\partial x_n} \boldsymbol{e}_n + 2 \mu \frac{\partial D_{ij}}{\partial x_i} \boldsymbol{e}_j \\
&= - \frac{\partial p}{\partial x_n} \boldsymbol{e}_n + \left( \kappa - \frac{2}{3} \mu \right) \frac{\partial}{\partial x_n} \left( \frac{\partial v_k}{\partial x_k} \right) \boldsymbol{e}_n + \mu \frac{\partial}{\partial x_i} \left( \frac{\partial v_i}{\partial x_j} + \frac{\partial v_j}{\partial x_i} \right) \boldsymbol{e}_j \\
&= - \frac{\partial p}{\partial x_n} \boldsymbol{e}_n + \left( \kappa + \frac{1}{3} \mu \right) \frac{\partial}{\partial x_n} \left( \frac{\partial v_k}{\partial x_k} \right) \boldsymbol{e}_n + \mu \frac{\partial^2 v_j}{\partial {x_i}^2} \boldsymbol{e}_j \\
&= - \nabla p + \mu \Delta \boldsymbol{v} + \left( \kappa + \frac{1}{3} \mu \right) \nabla (\nabla \cdot \boldsymbol{v})
\end{align*} ∇ ⋅ T = e l ∂ x l ∂ ⋅ [ { − p + ( κ − 3 2 μ ) tr D } I + 2 μ D ] = − e l ∂ x l ∂ p ⋅ δ mn ( e m ⊗ e n ) + e l ∂ x l ∂ ⋅ δ mn ( κ − 3 2 μ ) D kk ( e m ⊗ e n ) + e l ∂ x l ∂ ⋅ 2 μ D ij ( e i ⊗ e j ) = − δ mn ∂ x l ∂ p ( e l ⋅ e m ) e n + δ mn ( κ − 3 2 μ ) ∂ x l ∂ D kk ( e l ⋅ e m ) e n + 2 μ ∂ x l ∂ D ij ( e l ⋅ e i ) e j = − ∂ x n ∂ p e n + ( κ − 3 2 μ ) ∂ x n ∂ D kk e n + 2 μ ∂ x i ∂ D ij e j = − ∂ x n ∂ p e n + ( κ − 3 2 μ ) ∂ x n ∂ ( ∂ x k ∂ v k ) e n + μ ∂ x i ∂ ( ∂ x j ∂ v i + ∂ x i ∂ v j ) e j = − ∂ x n ∂ p e n + ( κ + 3 1 μ ) ∂ x n ∂ ( ∂ x k ∂ v k ) e n + μ ∂ x i 2 ∂ 2 v j e j = − ∇ p + μ Δ v + ( κ + 3 1 μ ) ∇ ( ∇ ⋅ v ) 非圧縮を仮定した場合∇ ⋅ v = 0 \nabla \cdot \boldsymbol{v}=0 ∇ ⋅ v = 0
ρ [ ∂ v ∂ t + ( v ⋅ ∇ v ) ] = − ∇ p + μ Δ v + ρ f \begin{equation}
% \label{eq:LandauLifshitzVol6_15.7}
\rho \left[ \frac{\partial \boldsymbol{v}}{\partial t} + (\boldsymbol{v} \cdot \nabla \boldsymbol{v}) \right] = - \nabla p + \mu \Delta \boldsymbol{v} + \rho \boldsymbol{f}
\end{equation} ρ [ ∂ t ∂ v + ( v ⋅ ∇ v ) ] = − ∇ p + μ Δ v + ρ f このとき、粘性係数はμ \mu μ ν = μ / ρ \nu = \mu / \rho ν = μ / ρ
References
久田 俊明, 野口 裕久, “非線形有限要素法の基礎と応用”, 丸善, 1996
久田 俊明, “非線形有限要素法のためのテンソル解析の基礎”, 丸善 , 1999
Karan S. Surana, “Classical Continuum Mechanics, Second Edition”, CRC Press, 2022