return 0;
wmake

sonicLiquidFoam解析

如果打不开图像,请右键在新标签页打开图像后刷新几次 如果打不开图像,请右键在新标签页打开图像后刷新几次


引言

在某些情况下,如密闭容器内注入液体、高压容器内泄压放出液体等工况需要考虑液体的可压缩性。虽然在这些工况下密度的变化是很小的,但是却会引起压力波的传递。如果液体流动方向在某一时刻发生骤变,也会引起压力波动以声速进行传递,即水锤。水锤是大量工程问题产生的主要原因,如噪音、震动损耗等,严重情况下,这种极高的压力传递会将管道击穿。一些工业器件的使用可以减少水锤的产生,如急泄阀、膨胀箱等。

sonicLiquidFoam为一个可压缩液体层流的求解器,可用于模拟水锤(Wang et al., 2017)、气动活塞(Suponitsky et al., 2017)等现象。在sonicLiquidFoam中,液体被认为是正压流体(Barotropic fluid)。其意味着液体的密度仅仅和压力有关而和温度无关(Massey and Ward-Smith, 1998)。在液体中,正压法则是适用的。然而在气体中则不然,气体中的密度和温度均有很大的关联性。正压法则在CFD中的具体体现为并不需要求解温度方程。因此sonicLiquidFoam为一个仅仅适用于液体的可压缩求解器。

如果打不开图像,请右键在新标签页打开图像后刷新几次 如果打不开图像,请右键在新标签页打开图像后刷新几次

图1. 左图:Suponitsky等模拟的气动活塞(Suponitsky et al., 2017)。右图:水锤的产生。当流体自由流动的时候,不会有压力波的传递。在阀门突然截止的时候,突然产生压力波的传递,即水锤。

方程离散以及求解

首先有可压缩流体的N-S方程:

\begin{equation} \frac{\partial \rho}{\partial t}+\nabla \cdot \left(\rho \mathbf{U} \right)=0 \label{NS1} \end{equation}
\begin{equation} \frac{\partial \rho \mathbf{U} }{\partial t}+\nabla \cdot \left(\rho \mathbf{U} \mathbf{U} \right)-\nabla \cdot \left(\mu \nabla \mathbf{U} \right)=-\nabla p \label{NS2} \end{equation}

有关方程\eqref{NS1}及方程\eqref{NS2}中的操作符含义以及其是如何离散的请参考icoFoam解析。下面对其进行简述。方程\eqref{NS1}在进行半离散化之后有:

\begin{equation} \mathbf{U}_\mathrm{P} = \frac{1}{{{A_\mathrm{P}}}}\left( { - \sum {{A_\mathrm{N}}\mathbf{U}_\mathrm{N}} + E_\mathrm{P}} \right)=\mathbf{HbyA}_\mathrm{P} \label{Up} \end{equation}

附加压力梯度之后有:

\begin{equation} \mathbf{U}_\mathrm{P} = \mathbf{HbyA}_\mathrm{P} -\frac{1}{{{A_\mathrm{P}}}}\nabla p \label{Upreal} \end{equation}

将方程\eqref{Upreal}代入到方程\eqref{NS1}有:

\begin{equation} \frac{\partial \rho}{\partial t}+\nabla \cdot \left(\rho \left( \mathbf{HbyA}_\mathrm{P} -\frac{1}{{{A_\mathrm{P}}}}\nabla p \right)\right)=0 \label{NS3} \end{equation}

对方程\eqref{NS3}进行分解:

\begin{equation} \frac{\partial \rho}{\partial t}+\nabla \cdot \left(\rho \mathbf{HbyA}_\mathrm{P} \right)-\nabla \cdot \left(\rho \frac{1}{{{A_\mathrm{P}}}}\nabla p \right)=0 \label{NS4} \end{equation}

依据正压法则有:

\begin{equation} \rho=A+\psi (p-B) \label{bar} \end{equation}

其中$A$为参考密度,$B$为参考压力。方程\eqref{NS4}可以继续化简为:

\begin{equation} \frac{\partial \psi p}{\partial t}+\nabla \cdot \left( \left(A+\psi (p-B) \right) \mathbf{HbyA}_\mathrm{P} \right)-\nabla \cdot \left(\rho \frac{1}{{{A_\mathrm{P}}}}\nabla p \right)=0 \label{NS5} \end{equation}
\begin{equation} \frac{\partial \psi p}{\partial t}+\nabla\cdot ((A-\psi B)\mathbf{HbyA}_\mathrm{P})+\nabla\cdot(\psi\mathbf{HbyA}_\mathrm{P}p)-\nabla \cdot \left(\rho \frac{1}{{{A_\mathrm{P}}}}\nabla p \right)=0 \label{NS6} \end{equation}

方程\eqref{NS6}即为求解的压力方程。其中的$\psi \mathbf{HbyA}_\mathrm{P}$即为代码中的phid。

代码分析

下面进入代码分析,首先有createFields:

volScalarField rho
(
    IOobject
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    rhoO + psi*p//对应方程\eqref{bar}
);

下面进入速度方程:

fvVectorMatrix UEqn
(
    fvm::ddt(rho, U)
  + fvm::div(phi, U)
  - fvm::laplacian(mu, U)
);//对应方程\eqref{NS2}

下面是压力方程:

fvScalarMatrix pEqn//对应方程\eqref{NS6}
(
    fvm::ddt(psi, p)
  + fvc::div(phi)
  + fvm::div(phid, p, "div(phid,p)")
  - fvm::laplacian(rho/A, p)
);

pEqn.solve();

需要注意,方程\eqref{NS6}第四项直接采用当前的密度,这样的操作可以视为一种线性化,否则此项会变为关于压力的非线性项。因此需要在压力循环内的压力方程求解后对密度再次求解进行更新。同时,由于压力波的传递,为了很好的捕获这种压力波,需要把时间步长设置的足够小。压力波在水中的传递以声速进行。水中的声速和$\psi$有关。在运行算例的时候,确保要保证依据声速计算的库朗数足够小。例如当$\psi=4e-7$的情况下,库朗数约为5e-7。在这样的时间步下,依据速度计算的库朗数会非常小。

参考文献

Massey, B. S., Ward-Smith, J., 1998. Mechanics of fluids. Crc Press.

Suponitsky, V., Plant, D., Avital, E. J., & Munjiza, A., 2017. Pressure Wave in Liquid Generated by Pneumatic Pistons and Its Interaction with a Free Surface. International Journal of Applied Mechanics, 9, 3710-3737.

Wang, C., Nilsson, H., Yang, J., Petit, O., 2017. 1D–3D coupling for hydraulic system transient simulations. Computer Physics Communications, 210, 1-9.

更新历史
2018.06.07对文章改版

东岳流体®版权所有
勘误、讨论、补充内容请前往CFD中文网

');