【物质导数】一文搞懂流体仿真中的物质导数

发布于:2022-11-09 ⋅ 阅读:(9) ⋅ 点赞:(0) ⋅ 评论:(0)

【物质导数】一文搞懂流体仿真中的物质导数


物理意义

物质导数

物质导数(Material Derivative)又称全导数、实质导数、随体导数。例如在使用拉格朗日法某一确定流体微团进行分析时,我们希望观察视角随着该微团流动,并求得该确定微团中某一物理属性仅对于时间的变化率。简单来说,即 “跟随微团运动时,观察到的某一物理属性相对于时间的变化率”。我们可以将相对于时间的物质导数定义为:

D D t = ∂ ∂ t + ( V ⋅ ∇ ) \frac{\mathrm{D}}{\mathrm{D}t} = \frac{\partial}{\partial t} + (\mathbf{V}\cdot\nabla) DtD=t+(V)

此处 V \mathbf{V} V 为速度矩阵,挂载了各个流体微团在此时刻的速度信息;等式左侧为全导数算子,右侧两项分别为当地导数对流导数
例如对该流体微团的密度参量 ρ = ρ ( x , y , z , t ) \bm{\rho} = \rho(x, y, z, t) ρ=ρ(x,y,z,t) 进行物质导数分析:

D ρ D t = ∂ ρ ∂ t + ( V ⋅ ∇ ) ρ \frac{\mathrm{D}\bm{\rho}}{\mathrm{D}t} = \frac{\partial\bm{\rho}}{\partial t} + (\mathbf{V}\cdot\nabla)\bm{\rho} DtDρ=tρ+(V)ρ

等式左侧为该微团密度关于时间的变化率,右侧第一项为由当地时间变化引起的密度变化率,右侧第二项为由该微团位移引起的密度变化率。在下面章节会使用易懂且尽可能严谨的方式拆分推理各项,并向各位解释其物理意义。

注意从数学意义上讲, D D t ≠ d d t \frac{\mathrm{D}}{\mathrm{D}t} \neq \frac{\mathrm{d}}{\mathrm{d}t} DtD=dtd,在 Survey SPH Fluids in Computer Graphics 一文中存在数学符号书写错误。此外明显的, D D t ≠ ∂ ∂ t \frac{\mathrm{D}}{\mathrm{D}t} \neq \frac{\partial}{\partial t} DtD=t。如对 Nabla 算子 ∇ \nabla 的作用存在疑问,读者可移步至 一文读懂 Nabla 算子 一文补充前导数学知识。

当地导数

当地导数又称时变导数,即由当地时间变化引起的物理属性变化率。当观察者固定在某一世界空间(称之为当地)分析某一物理属性时,空间位置信息静态不变,该视角下物理属性的变化仅随着当地时间改变而改变。如将 ρ = ρ ( x , y , z , t ) \bm{\rho} = \rho(x, y, z, t) ρ=ρ(x,y,z,t) 中的 x , y , z x, y, z x,y,z 视为常参, t t t 视为变参,我们就可以对该物理属性求关于时间的偏导数,如 ∂ ρ ∂ t \frac{\partial\bm{\rho}}{\partial t} tρ ,来作为当地导数项。

对流导数

对流导数又称位变导数,即由此时刻该微团位移引起的物理属性变化率。当观察者凝滞在此时刻并跟随该微团移动分析该微团的某一物理属性时,忽略时间因素,该视角下物理属性的变化仅随着物质的位移而改变。如将 ρ = ρ ( x , y , z , t ) \bm{\rho} = \rho(x, y, z, t) ρ=ρ(x,y,z,t) 中的 x , y , z x, y, z x,y,z 视为变参, t t t 视为常参,我们就可以使用其速度点乘该微团的空间梯度,如 ( V ⋅ ∇ ) ρ (\mathbf{V}\cdot\nabla)\bm{\rho} (V)ρ ,来作为对流导数项。

其中:
( V ⋅ ∇ ) ρ = ( ( V x , V y , V z ) ⋅ ( ∂ ∂ x , ∂ ∂ y , ∂ ∂ z ) ) ρ = ( V x ∂ ∂ x + V y ∂ ∂ y + V z ∂ ∂ z ) ρ = V x ∂ ρ ∂ x + V y ∂ ρ ∂ y + V z ∂ ρ ∂ z (\mathbf{V}\cdot\nabla)\bm{\rho} = \big((\mathbf{V_x}, \mathbf{V_y}, \mathbf{V_z})\cdot(\frac{\partial}{\partial x}, \frac{\partial}{\partial y}, \frac{\partial}{\partial z})\big)\bm{\rho} = (\mathbf{V_x}\frac{\partial}{\partial x} + \mathbf{V_y}\frac{\partial}{\partial y} + \mathbf{V_z}\frac{\partial}{\partial z})\bm{\rho} = \mathbf{V_x}\frac{\partial\bm{\rho}}{\partial x} + \mathbf{V_y}\frac{\partial\bm{\rho}}{\partial y} + \mathbf{V_z}\frac{\partial\bm{\rho}}{\partial z} (V)ρ=((Vx,Vy,Vz)(x,y,z))ρ=(Vxx+Vyy+Vzz)ρ=Vxxρ+Vyyρ+Vzzρ
或者对于某一微团:
( v ⋅ ∇ ) ρ = u ∂ ρ ∂ x + v ∂ ρ ∂ y + w ∂ ρ ∂ z v = u i ^ + v j ^ + w k ^ (\mathbf{v}\cdot\nabla)\rho = u\frac{\partial\rho}{\partial x} + v\frac{\partial\rho}{\partial y} + w\frac{\partial\rho}{\partial z} \\ \mathbf{v} = u\mathbf{\hat{i}} + v\mathbf{\hat{j}} + w\mathbf{\hat{k}} (v)ρ=uxρ+vyρ+wzρv=ui^+vj^+wk^


实例类比

以宏观世界举例子,假设在一辆汽车上,一个温度计从北京被运往哈尔滨,中午出发晚上抵达。现在我们以拉格朗日法,分析流体微团“汽车”的物理属性“温度计读数”物质导数可以简单理解为 汽车内“温度计读数”的变化:其中由于中午到晚上时间温差带来的读数变化为当地导数;由于北京哈尔滨地域差异带来的读数变化为对流导数

当然,这不是微观数学条件下的例子。若逐渐缩短北京到哈尔滨的距离,使得运动时间趋向于无穷小,那么中午到晚上温差可完全视为北京当地导数


简要推导

物质导数
我们选取如图所示的流体微团,分析其密度 ρ \rho ρ 从时刻 t = t 1 t = t_1 t=t1 t = t 2 t = t_2 t=t2 的物质导数。如第一节,我们可以将该微团的速度表示为矢量 v = u i ^ + v j ^ + w k ^ \mathbf{v} = u\mathbf{\hat{i}} + v\mathbf{\hat{j}} + w\mathbf{\hat{k}} v=ui^+vj^+wk^。不难得出, ρ , u , v , w \rho, u, v, w ρ,u,v,w 均是关于三维空间与一维时间的四参函数: ρ = ρ ( x , y , z , t ) ; u = u ( x , y , z , t ) ; ⋯   ; w = w ( x , y , z , t ) \rho = \rho(x, y, z, t); u = u(x, y, z, t); \cdots; w = w(x, y, z, t) ρ=ρ(x,y,z,t);u=u(x,y,z,t);;w=w(x,y,z,t)

  • t = t 1 t = t_1 t=t1 时, ρ 1 = ρ ( x 1 , y 1 , z 1 , t 1 ) \rho_1 = \rho(x_1, y_1, z_1, t_1) ρ1=ρ(x1,y1,z1,t1)
  • t = t 2 t = t_2 t=t2 时, ρ 2 = ρ ( x 2 , y 2 , z 2 , t 2 ) \rho_2 = \rho(x_2, y_2, z_2, t_2) ρ2=ρ(x2,y2,z2,t2)
  • ρ 2 \rho_2 ρ2 ρ 1 \rho_1 ρ1 处进行泰勒展开:
    • ρ 2 = ρ 1 + ( ∂ ρ ∂ x ) 1 ( x 2 − x 1 ) + ( ∂ ρ ∂ y ) 1 ( y 2 − y 1 ) + ( ∂ ρ ∂ z ) 1 ( z 2 − z 1 ) + ( ∂ ρ ∂ t ) 1 ( t 2 − t 1 ) + o ( ⋅ ) \rho_2 = \rho_1 + (\frac{\partial\rho}{\partial x})_1(x_2-x_1) + (\frac{\partial\rho}{\partial y})_1(y_2-y_1) + (\frac{\partial\rho}{\partial z})_1(z_2-z_1) + (\frac{\partial\rho}{\partial t})_1(t_2-t_1) + o(\cdot) ρ2=ρ1+(xρ)1(x2x1)+(yρ)1(y2y1)+(zρ)1(z2z1)+(tρ)1(t2t1)+o()
    • o ( ⋅ ) o(\cdot) o() 为与 x , y , z , t x, y, z, t x,y,z,t 相关的高阶无穷小。
  • 移项并忽略高阶无穷小,等式两边同时除以 t 2 − t 1 t_2 - t_1 t2t1
    • ρ 2 − ρ 1 = ( ∂ ρ ∂ x ) 1 ( x 2 − x 1 ) + ( ∂ ρ ∂ y ) 1 ( y 2 − y 1 ) + ( ∂ ρ ∂ z ) 1 ( z 2 − z 1 ) + ( ∂ ρ ∂ t ) 1 ( t 2 − t 1 ) \rho_2 - \rho_1 = (\frac{\partial\rho}{\partial x})_1(x_2-x_1) + (\frac{\partial\rho}{\partial y})_1(y_2-y_1) + (\frac{\partial\rho}{\partial z})_1(z_2-z_1) + (\frac{\partial\rho}{\partial t})_1(t_2-t_1) ρ2ρ1=(xρ)1(x2x1)+(yρ)1(y2y1)+(zρ)1(z2z1)+(tρ)1(t2t1)
    • ρ 2 − ρ 1 t 2 − t 1 = ( ∂ ρ ∂ x ) 1 x 2 − x 1 t 2 − t 1 + ( ∂ ρ ∂ y ) 1 y 2 − y 1 t 2 − t 1 + ( ∂ ρ ∂ z ) 1 z 2 − z 1 t 2 − t 1 + ∂ ρ ∂ t \frac{\rho_2 - \rho_1}{t_2 - t_1} = (\frac{\partial\rho}{\partial x})_1\frac{x_2-x_1}{t_2 - t_1} + (\frac{\partial\rho}{\partial y})_1\frac{y_2-y_1}{t_2 - t_1} + (\frac{\partial\rho}{\partial z})_1\frac{z_2-z_1}{t_2 - t_1} + \frac{\partial\rho}{\partial t} t2t1ρ2ρ1=(xρ)1t2t1x2x1+(yρ)1t2t1y2y1+(zρ)1t2t1z2z1+tρ
  • t 2 t_2 t2 右趋近于 t 1 t_1 t1 时:
    • D ρ D t = lim ⁡ t 2 → t 1 + ρ 2 − ρ 1 t 2 − t 1 = ( ∂ ρ ∂ x ) 1 u + ( ∂ ρ ∂ y ) 1 v + ( ∂ ρ ∂ z ) 1 w + ∂ ρ ∂ t = ( v ⋅ ∇ ) ρ + ∂ ρ ∂ t \frac{\mathrm{D}\rho}{\mathrm{D}t} = \lim\limits_{t_2 \to t_1^+} \frac{\rho_2 - \rho_1}{t_2 - t_1} = (\frac{\partial\rho}{\partial x})_1u + (\frac{\partial\rho}{\partial y})_1v + (\frac{\partial\rho}{\partial z})_1w + \frac{\partial\rho}{\partial t} = (\mathbf{v}\cdot\nabla)\rho + \frac{\partial \rho}{\partial t} DtDρ=t2t1+limt2t1ρ2ρ1=(xρ)1u+(yρ)1v+(zρ)1w+tρ=(v)ρ+tρ
    • D ρ D t = u ∂ ρ ∂ x + v ∂ ρ ∂ y + w ∂ ρ ∂ z + ∂ ρ ∂ t = ( v ⋅ ∇ ) ρ + ∂ ρ ∂ t \frac{\mathrm{D}\rho}{\mathrm{D}t} = u\frac{\partial\rho}{\partial x} + v\frac{\partial\rho}{\partial y} + w\frac{\partial\rho}{\partial z} + \frac{\partial\rho}{\partial t} = (\mathbf{v}\cdot\nabla)\rho + \frac{\partial \rho}{\partial t} DtDρ=uxρ+vyρ+wzρ+tρ=(v)ρ+tρ

可推得物质导数为:
D ρ D t = ∂ ρ ∂ t + ( V ⋅ ∇ ) ρ \frac{\mathrm{D}\bm{\rho}}{\mathrm{D}t} = \frac{\partial\bm{\rho}}{\partial t} + (\mathbf{V}\cdot\nabla)\bm{\rho} DtDρ=tρ+(V)ρ