本文内容主要是不可压的雷诺平均方程(Reynolds-averaged Navier-Stokes equations,RANS)的推导。在推导雷诺方程之前先来总结一些关于时间平均的运算规律。

1. 关于时间平均的运算规律证明

  对于瞬时量 ϕ = Φ + ϕ ′ \phi=\Phi+\phi^\prime ϕ=Φ+ϕ ψ = Ψ + ψ ′ \psi=\Psi+\psi^\prime ψ=Ψ+ψ Φ \Phi Φ Ψ \Psi Ψ已经是与时间无关的量,所以根据时间平均的定义有
ϕ ′ ‾ = ψ ′ ‾ = 0 (1.1) \overline{\phi^\prime}=\overline{\psi^\prime}=0 \tag{1.1} ϕ=ψ=0(1.1)
Φ ‾ = 1 Δ t ∫ 0 Δ t Φ d t = Φ 1 Δ t ∫ 0 Δ t d t = Φ (1.2) \overline{\Phi}=\frac{1}{\Delta t} \int_0^{\Delta t} \Phi dt =\Phi \frac{1}{\Delta t} \int_0^{\Delta t} dt=\Phi \tag{1.2} Φ=Δt10ΔtΦdt=ΦΔt10Δtdt=Φ(1.2)
∂ ϕ ∂ s ‾ = 1 Δ t ∫ 0 Δ t ∂ ϕ ∂ s d t = ∂ ∂ s ( 1 Δ t ∫ 0 Δ t ϕ d t ) = ∂ ϕ ˉ ∂ s = ∂ Φ ∂ s (1.3) \begin{aligned} \overline{\frac{\partial \phi}{\partial s}} &= \frac{1}{\Delta t} \int_0^{\Delta t} \frac{\partial \phi}{\partial s} dt = \frac{\partial}{\partial s} \left(\frac{1}{\Delta t} \int_0^{\Delta t} \phi dt \right) = \frac{\partial \bar{\phi}}{\partial s} = \frac{\partial \Phi}{\partial s} \end{aligned} \tag{1.3} sϕ=Δt10Δtsϕdt=s(Δt10Δtϕdt)=sϕˉ=sΦ(1.3)
∫ ϕ d s ‾ = 1 Δ t ∫ 0 Δ t ( ∫ ϕ d s ) d t = ∫ ( 1 Δ t ∫ 0 Δ t ϕ d t ) d s = ∫ ϕ ˉ d s = ∫ Φ d s (1.4) \begin{aligned} \overline{\int \phi ds} &= \frac{1}{\Delta t}\int_0^{\Delta t}\left( \int \phi ds \right)dt =\int \left( \frac{1}{\Delta t}\int_0^{\Delta t} \phi dt\right)ds =\int \bar{\phi}ds = \int\Phi ds \end{aligned} \tag{1.4} ϕds=Δt10Δt(ϕds)dt=(Δt10Δtϕdt)ds=ϕˉds=Φds(1.4)
ϕ + ψ ‾ = 1 Δ t ∫ 0 Δ t ( ϕ + ψ ) d t = 1 Δ t ∫ 0 Δ t ϕ d t + 1 Δ t ∫ 0 Δ t ψ d t = Φ + Ψ (1.5) \begin{aligned} \overline{\phi + \psi} &= \frac{1}{\Delta t} \int_0^{\Delta t} (\phi + \psi)dt = \frac{1}{\Delta t} \int_0^{\Delta t} \phi dt + \frac{1}{\Delta t} \int_0^{\Delta t} \psi dt = \Phi + \Psi \end{aligned} \tag{1.5} ϕ+ψ=Δt10Δt(ϕ+ψ)dt=Δt10Δtϕdt+Δt10Δtψdt=Φ+Ψ(1.5)
ϕ Ψ ‾ = 1 Δ t ∫ 0 Δ t ( ϕ Ψ ) d t = Ψ 1 Δ t ∫ 0 Δ t ϕ d t = Φ Ψ (1.6) \overline{\phi\Psi} = \frac{1}{\Delta t } \int_0^{\Delta t} (\phi \Psi) dt = \Psi \frac{1}{\Delta t} \int_0^{\Delta t} \phi dt = \Phi\Psi \tag{1.6} ϕΨ=Δt10Δt(ϕΨ)dt=ΨΔt10Δtϕdt=ΦΨ(1.6)
ϕ ′ Ψ ‾ = 1 Δ t ∫ 0 Δ t ( ϕ ′ Ψ ) d t = Ψ 1 Δ t ∫ 0 Δ t ϕ ′ d t = Ψ ϕ ′ ˉ = 0 (1.7) \overline{\phi^\prime\Psi} = \frac{1}{\Delta t} \int_0^{\Delta t} (\phi^\prime \Psi) dt = \Psi \frac{1}{\Delta t} \int_0^{\Delta t} \phi^\prime dt = \Psi \bar{\phi^\prime} = 0 \tag{1.7} ϕΨ=Δt10Δt(ϕΨ)dt=ΨΔt10Δtϕdt=Ψϕˉ=0(1.7)
ϕ ψ ‾ = ( Φ + ϕ ′ ) ( Ψ + ψ ′ ) ‾ = Φ Ψ + Φ ψ ′ + ϕ ′ Ψ + ϕ ′ ψ ′ ‾ = Φ Ψ ‾ + Φ ψ ′ ‾ + ϕ ′ Ψ ‾ + ϕ ′ ψ ′ ‾ = Φ Ψ ‾ + 0 + 0 + ϕ ′ ψ ′ ‾ = Φ Ψ + ϕ ′ ψ ′ ‾ (1.8) \begin{aligned} \overline{\phi \psi} = \overline{(\Phi + \phi^\prime)(\Psi + \psi^\prime)} &= \overline{\Phi\Psi + \Phi\psi^\prime + \phi^\prime \Psi + \phi^\prime \psi^\prime} \\ &= \overline{\Phi\Psi} + \overline{\Phi\psi^\prime} + \overline{ \phi^\prime \Psi} + \overline{\phi^\prime \psi^\prime} \\ &= \overline{\Phi\Psi} + 0 + 0 + \overline{\phi^\prime \psi^\prime} \\ & = \Phi\Psi + \overline{\phi^\prime \psi^\prime} \end{aligned} \tag{1.8} ϕψ=(Φ+ϕ)(Ψ+ψ)=ΦΨ+Φψ+ϕΨ+ϕψ=ΦΨ+Φψ+ϕΨ+ϕψ=ΦΨ+0+0+ϕψ=ΦΨ+ϕψ(1.8)
  以上均为标量的运算,矢量的散度和梯度计算可以通过上述公式扩展而来。
  本文中均以粗体字母代表矢量,如速度矢量 u \bold u u和其他任意矢量 a \bold a a等,以 ∇ ⋅ u \nabla \cdot \bold u u表示矢量 u \bold u u的散度,以 ∇ ϕ \nabla \phi ϕ表示标量 ϕ \phi ϕ的梯度。根据雷诺分解,瞬时速度矢量 u = U + u ′ \bold u = \bold U + \bold u^\prime u=U+u和瞬时标量 ϕ = Φ + ϕ ′ \phi = \Phi + \phi^\prime ϕ=Φ+ϕ有如下运算规律:
u = ( u , v , w ) T \bold u = (u,v,w)^T u=(u,v,w)T

u ˉ = U = ( u ˉ , v ˉ , w ˉ , ) = ( U , V , W ) T \bar \bold u = \bold U =(\bar u,\bar v,\bar w,)= (U,V,W)^T uˉ=U=(uˉ,vˉ,wˉ,)=(U,V,W)T

u ′ = ( u ′ , v ′ , w ′ ) T \bold u^\prime = (u^\prime,v^\prime,w^\prime)^T u=(u,v,w)T
∇ ⋅ u ‾ = ∂ u ∂ x + ∂ v ∂ y + ∂ w ∂ z ‾ = ∂ u ˉ ∂ x + ∂ v ˉ ∂ y + ∂ w ˉ ∂ z = ∂ U ∂ x + ∂ V ∂ y + ∂ W ∂ z = ∇ ⋅ U = ∇ ⋅ u ˉ (1.9) \begin{aligned} \overline{\nabla \cdot \bold u} &= \overline{\frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} } = \frac{\partial \bar u}{\partial x} + \frac{\partial \bar v}{\partial y} + \frac{\partial \bar w}{\partial z}\\ &= \frac{\partial U}{\partial x} + \frac{\partial V}{\partial y} + \frac{\partial W}{\partial z} =\nabla \cdot \bold U = \nabla \cdot \bar \bold u \end{aligned} \tag{1.9} u=xu+yv+zw=xuˉ+yvˉ+zwˉ=xU+yV+zW=U=uˉ(1.9)
ϕ u ‾ = ( ϕ u ‾ , ϕ v ‾ , ϕ w ‾ ) T = [ Φ U + ϕ ′ u ′ ‾ Φ V + ϕ ′ v ′ ‾ Φ W + ϕ ′ w ′ ‾ ] = Φ U + ϕ ′ u ′ ‾ (1.10) \overline{\phi \bold u} = (\overline{\phi u},\overline{\phi v}, \overline{\phi w})^T= \left[ \begin{matrix} \Phi U + \overline{\phi^\prime u^\prime } \\ \\ \Phi V +\overline{\phi^\prime v^\prime } \\ \\ \Phi W +\overline{\phi^\prime w^\prime } \end{matrix} \right] = \Phi \bold U +\overline{\phi^\prime \bold u^\prime} \tag{1.10} ϕu=(ϕu,ϕv,ϕw)T=ΦU+ϕuΦV+ϕvΦW+ϕw=ΦU+ϕu(1.10)
∇ ⋅ ( ϕ u ) ‾ = ∇ ⋅ ϕ u ‾ = ∇ ⋅ ( Φ U ) + ∇ ⋅ ( ϕ ′ u ′ ‾ ) (1.11) \overline{\nabla \cdot (\phi \bold u)} = \nabla \cdot \overline{\phi \bold u} =\nabla \cdot (\Phi \bold U) + \nabla \cdot (\overline{\phi^\prime \bold u^\prime}) \tag{1.11} (ϕu)=ϕu=(ΦU)+(ϕu)(1.11)
∇ ϕ ‾ = ( ∂ ϕ ∂ x , ∂ ϕ ∂ y , ∂ ϕ ∂ z ) T ‾ = ( ∂ ϕ ˉ ∂ x , ∂ ϕ ˉ ∂ x , ∂ ϕ ˉ ∂ x ) T = ( ∂ Φ ∂ x , ∂ Φ ∂ x , ∂ Φ ∂ x ) T = ∇ Φ (1.12) \overline{\nabla \phi} =\overline{\left(\frac{\partial \phi}{\partial x},\frac{\partial \phi}{\partial y},\frac{\partial \phi}{\partial z}\right)^T} = \left( \frac{ \partial \bar \phi}{\partial x} ,\frac{ \partial \bar \phi}{\partial x},\frac{ \partial \bar \phi}{\partial x} \right)^T =\left( \frac{\partial \Phi}{\partial x} ,\frac{\partial \Phi}{\partial x},\frac{\partial \Phi}{\partial x}\right)^T = \nabla \Phi \tag{1.12} ϕ=(xϕ,yϕ,zϕ)T=(xϕˉ,xϕˉ,xϕˉ)T=(xΦ,xΦ,xΦ)T=Φ(1.12)
∇ ⋅ ∇ ϕ ‾ = ∇ ⋅ ∇ ϕ ‾ = ∇ ⋅ ∇ Φ (1.13) \overline{\nabla \cdot \nabla \phi} = \nabla \cdot \overline{\nabla \phi} = \nabla \cdot \nabla \Phi \tag{1.13} ϕ=ϕ=Φ(1.13)

∇ ⋅ ∇ \nabla \cdot \nabla 也可写作 ∇ 2 \nabla^2 2,即拉普拉斯算子。

2. 雷诺平均方程的推导过程

  有了上述准备,推导雷诺平均方程就容易多了,这里考虑三维笛卡尔坐标系下瞬态不可压的连续性方程和N-S方程,忽略体积力,方程如下:
∇ ⋅ u = 0 (2.1) \nabla \cdot \bold u = 0 \tag{2.1} u=0(2.1)

∂ u ∂ t + ∇ ⋅ ( u u ) = − 1 ρ ∂ p ∂ x + ν ∇ 2 u (2.2a) \frac{\partial u}{\partial t} + \nabla \cdot (u \bold u) = -\frac{1}{\rho} \frac{\partial p}{\partial x} + \nu \nabla^2 u \tag{2.2a} tu+(uu)=ρ1xp+ν2u(2.2a)

∂ v ∂ t + ∇ ⋅ ( v u ) = − 1 ρ ∂ p ∂ y + ν ∇ 2 v (2.2b) \frac{\partial v}{\partial t} + \nabla \cdot (v \bold u) = -\frac{1}{\rho} \frac{\partial p}{\partial y} + \nu \nabla^2 v \tag{2.2b} tv+(vu)=ρ1yp+ν2v(2.2b)

∂ w ∂ t + ∇ ⋅ ( w u ) = − 1 ρ ∂ p ∂ z + ν ∇ 2 w (2.2c) \frac{\partial w}{\partial t} + \nabla \cdot (w \bold u) = -\frac{1}{\rho} \frac{\partial p}{\partial z} + \nu \nabla^2 w \tag{2.2c} tw+(wu)=ρ1zp+ν2w(2.2c)
  雷诺平均方程的思想就是先将方程中的流动变量分解成平均量和脉动量之和,然后对方程两边同时取时间平均。根据雷诺分解,方程中的各瞬时量分解如下:
u = U + u ′ u = U + u ′ v = V + v ′ w = W + w ′ p = P + p ′ (2.3) \bold u = \bold U + \bold u^\prime \quad u=U+u^\prime \quad v=V+v^\prime \quad w=W+w^\prime \quad p=P + p^\prime \tag{2.3} u=U+uu=U+uv=V+vw=W+wp=P+p(2.3)
( 2.3 ) (2.3) (2.3)各式代入到式 ( 2.1 ) (2.1) (2.1)~ ( 2.2 c ) (2.2c) (2.2c)中,然后对方程取时间平均,再根据第1节中的运算规律对方程进行整理。

先来推导方程 ( 2.1 ) (2.1) (2.1),根据式 ( 1.9 ) (1.9) (1.9)
∇ ⋅ u ‾ = ∇ ⋅ u ‾ = ∇ ⋅ U \overline{\nabla \cdot \bold u} = \nabla \cdot \overline \bold u = \nabla \cdot \bold U u=u=U
  于是得到平均流的连续性方程:
∇ ⋅ U = 0 (2.4) \nabla \cdot \bold U = 0 \tag{2.4} U=0(2.4)
  其中 U \bold U U是平均流速度矢量,写成分量形式为
∂ U ∂ x + ∂ V ∂ y + ∂ W ∂ z = 0 \frac{\partial U}{\partial x} + \frac{\partial V}{\partial y} +\frac{\partial W}{\partial z} = 0 xU+yV+zW=0

推导方程 ( 2.2 a ) (2.2a) (2.2a)
  方程左边第一项,根据式 ( 1.3 ) (1.3) (1.3)
∂ u ∂ t ‾ = ∂ U ∂ t \overline{\frac{\partial u}{\partial t}} = \frac{\partial U}{\partial t} tu=tU
  左边第二项,根据式 ( 1.11 ) (1.11) (1.11)
∇ ⋅ ( u u ) ‾ = ∇ ⋅ ( U U ) + ∇ ⋅ ( u ′ u ′ ‾ ) \overline{\nabla \cdot (u\bold u)} = \nabla \cdot (U\bold U) + \nabla \cdot (\overline{u^\prime \bold u^\prime}) (uu)=(UU)+(uu)
  方程右边第一项,
− 1 ρ ∂ p ∂ x ‾ = − 1 ρ ∂ p ˉ ∂ x = − 1 p ∂ P ∂ x \overline{-\frac{1}{\rho} \frac{\partial p}{\partial x}} = -\frac{1}{\rho} \frac{\partial \bar p}{\partial x} = -\frac{1}{p} \frac{\partial P}{\partial x} ρ1xp=ρ1xpˉ=p1xP
  右边第二项,根据式 ( 1.13 ) (1.13) (1.13)
ν ∇ 2 u ‾ = ν ∇ 2 U \overline{\nu \nabla^2 u} = \nu \nabla^2 U ν2u=ν2U
  综上,对方程 ( 2.2 a ) (2.2a) (2.2a)取时间平均得到 x x x分量的时均动量方程:

∂ U ∂ t + ∇ ⋅ ( U U ) + ∇ ⋅ ( u ′ u ′ ‾ ) = − 1 ρ ∂ P ∂ x + ν ∇ 2 U ( 1 ) ( 2 ) ( 3 ) ( 4 ) ( 5 ) (2.5a) \begin{aligned} \frac{\partial U}{\partial t} + \nabla \cdot (U \bold U) + \nabla \cdot (\overline{u^\prime \bold u^\prime}) = -\frac{1}{\rho} \frac{\partial P}{\partial x} + \nu \nabla^2U \tag{2.5a} \\ (1) \quad \qquad (2) \qquad \qquad (3) \qquad \qquad (4) \qquad \quad(5) \end{aligned} tU+(UU)+(uu)=ρ1xP+ν2U(1)(2)(3)(4)(5)(2.5a)

同样的推导过程可以得到 y y y z z z分量的时均动量方程:
∂ V ∂ t + ∇ ⋅ ( V U ) + ∇ ⋅ ( v ′ u ′ ‾ ) = − 1 ρ ∂ P ∂ y + ν ∇ 2 V (2.5b) \frac{\partial V}{\partial t} + \nabla \cdot (V \bold U) + \nabla \cdot (\overline{v^\prime \bold u^\prime}) = -\frac{1}{\rho} \frac{\partial P}{\partial y} + \nu \nabla^2V \tag{2.5b} tV+(VU)+(vu)=ρ1yP+ν2V(2.5b)

∂ W ∂ t + ∇ ⋅ ( W U ) + ∇ ⋅ ( w ′ u ′ ‾ ) = − 1 ρ ∂ P ∂ z + ν ∇ 2 W (2.5c) \frac{\partial W}{\partial t} + \nabla \cdot (W \bold U) + \nabla \cdot (\overline{w^\prime \bold u^\prime}) = -\frac{1}{\rho} \frac{\partial P}{\partial z} + \nu \nabla^2W \tag{2.5c} tW+(WU)+(wu)=ρ1zP+ν2W(2.5c)
以方程 ( 2.5 a ) (2.5a) (2.5a)为例,可以看到第(2)、(3)项都是来自原方程的对流项,第(3)项是湍流脉动量的乘积,是时间平均后多出来的项,其代表的物理含义是湍流涡的对流所引起动量输运,这就类似于流体在粘性剪应力的作用下又附加了湍流剪应力。所以习惯上,这一项会移到方程的右边和粘性项放在一起,以反应其物理属性。为了展示湍流应力的结构,把第(3)项中的脉动速度矢量写出分量形式,并移到等号右边,则有
∂ U ∂ t + ∇ ⋅ ( U U ) = − 1 ρ ∂ P ∂ x + ν ∇ 2 U + 1 ρ [ ∂ ( − ρ u ′ 2 ‾ ) ∂ x + ∂ ( − ρ u ′ v ′ ‾ ) ∂ y + ∂ ( − ρ u ′ w ′ ‾ ) ∂ z ] (2.6a) \begin{aligned} \frac{\partial U}{\partial t} + \nabla \cdot (U \bold U) =& -\frac{1}{\rho} \frac{\partial P}{\partial x}+\nu \nabla^2U \\ &+\frac{1}{\rho}\left[ \frac{\partial (-\rho \overline{u^{\prime 2}})}{\partial x} + \frac{\partial (-\rho \overline{u^{\prime }v^\prime})}{\partial y} +\frac{\partial (-\rho \overline{u^{\prime }w^\prime})}{\partial z} \right] \tag{2.6a} \end{aligned} tU+(UU)=ρ1xP+ν2U+ρ1[x(ρu2)+y(ρuv)+z(ρuw)](2.6a)
∂ V ∂ t + ∇ ⋅ ( V U ) = − 1 ρ ∂ P ∂ y + ν ∇ 2 V + 1 ρ [ ∂ ( − ρ u ′ v ′ ‾ ) ∂ x + ∂ ( − ρ v ′ 2 ‾ ) ∂ y + ∂ ( − ρ v ′ w ′ ‾ ) ∂ z ] (2.6b) \begin{aligned} \frac{\partial V}{\partial t} + \nabla \cdot (V \bold U) =& -\frac{1}{\rho} \frac{\partial P}{\partial y}+\nu \nabla^2V \\ & +\frac{1}{\rho}\left[ \frac{\partial (-\rho \overline{u^{\prime} v^\prime})}{\partial x} + \frac{\partial (-\rho \overline{v^{\prime2} })}{\partial y} +\frac{\partial (-\rho \overline{v^{\prime }w^\prime})}{\partial z} \right] \tag{2.6b} \end{aligned} tV+(VU)=ρ1yP+ν2V+ρ1[x(ρuv)+y(ρv2)+z(ρvw)](2.6b)
∂ W ∂ t + ∇ ⋅ ( W U ) = − 1 ρ ∂ P ∂ z + ν ∇ 2 W + 1 ρ [ ∂ ( − ρ u ′ w ′ ‾ ) ∂ x + ∂ ( − ρ v ′ w ′ ‾ ) ∂ y + ∂ ( − ρ w ′ 2 ‾ ) ∂ z ] (2.6c) \begin{aligned} \frac{\partial W}{\partial t} + \nabla \cdot (W \bold U) =& -\frac{1}{\rho} \frac{\partial P}{\partial z}+\nu \nabla^2W \\ & +\frac{1}{\rho}\left[ \frac{\partial (-\rho \overline{u^{\prime} w^\prime})}{\partial x} + \frac{\partial (-\rho \overline{v^{\prime }w^\prime})}{\partial y} +\frac{\partial (-\rho \overline{w^{\prime2 }})}{\partial z} \right] \tag{2.6c} \end{aligned} tW+(WU)=ρ1zP+ν2W+ρ1[x(ρuw)+y(ρvw)+z(ρw2)](2.6c)
湍流应力一共有6个独立分量,包括3个法向应力分量:
τ x x = − ρ u ′ 2 ‾ τ y y = − ρ v ′ 2 ‾ τ z z = − ρ w ′ 2 ‾ \tau_{xx}=-\rho \overline{u^{\prime 2}} \quad \tau_{yy}=-\rho \overline{v^{\prime 2}} \quad \tau_{zz}=-\rho \overline{w^{\prime 2}} τxx=ρu2τyy=ρv2τzz=ρw2
和3个切应力分量:
τ x y = τ y x = − ρ u ′ v ′ ‾ τ x z = τ z x = − ρ u ′ w ′ ‾ τ z y = τ y z = − ρ w ′ v ′ ‾ \tau_{xy}=\tau_{yx}=-\rho \overline{u^\prime v^\prime} \quad \tau_{xz}=\tau_{zx}=-\rho \overline{u^\prime w^\prime} \quad \tau_{zy}=\tau_{yz}=-\rho \overline{w^\prime v^\prime} \quad τxy=τyx=ρuvτxz=τzx=ρuwτzy=τyz=ρwv
这些湍流应力称为雷诺应力,方程式 ( 2.4 ) (2.4) (2.4) ( 2.6 a ) (2.6a) (2.6a)~ ( 2.6 c ) (2.6c) (2.6c)称为雷诺平均N-S方程(RANS)。

参考资料:
  1. Versteeg H K , Malalasekera W . An introduction to computational fluid dynamics : the finite volume method = 计算流体动力学导论[M]. 世界图书出版公司, 2010.
Logo

旨在为数千万中国开发者提供一个无缝且高效的云端环境,以支持学习、使用和贡献开源项目。

更多推荐