• 返回主页
  • sonicLiquidFoam解析


    2016.11.23,Update:统一HbyA符号;

    求解器类型:液体可压缩层流求解器。

    在某些情况下如密闭容器内注入液体、高压容器内放出液体等工况需要考虑液体的可压缩性。虽然在这些工况下密度的变化是很小的,但是却会引起压力波的传递。sonicLiquidFoam为一个可压缩液体层流的求解器。在sonicLiquidFoam中,液体的密度和压力的关系用正压法则来表示。正压法则表明液体的密度仅仅和压力有关而和温度无关。在液体中正压法则是适用的。然而在气体中则不然,气体中的密度和温度均有很大的关联性。正压法则在CFD中的具体体现为并不需要求解温度方程。因此sonicLiquidFoam为一个仅仅适用于液体的可压缩求解器。

    首先有可压缩流体的NS方程: \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} 具体的有关方程($\ref{NS1}$)以及方程($\ref{NS2}$)中的操作符含义以及其是如何离散的请参考icoFoam解析。下面对其进行简述。

    方程($\ref{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} 将方程($\ref{Upreal}$)代入到方程($\ref{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} 对方程($\ref{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} 依据正压法则(在这里假定$p_0$的值相对$p$足够小): \begin{equation} \rho=\rho_0+\psi (p-p_0) \approx \rho_0+\psi p \label{bar} \end{equation} 方程($\ref{NS4}$)可以继续化简为: \begin{equation} \frac{\partial \psi p}{\partial t}+\nabla \cdot \left( \left(\rho_0+\psi p \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 \left(\frac{\rho_0}{\psi}\psi \mathbf{HbyA}_\mathrm{P} \right)+\nabla \cdot \left( \psi \mathbf{HbyA}_\mathrm{P} \cdot p \right)-\nabla \cdot \left(\rho \frac{1}{{{A_\mathrm{P}}}}\nabla p \right)=0 \label{NS6} \end{equation} 方程($\ref{NS6}$)即为求解的压力方程。其中的$\psi \mathbf{HbyA}_\mathrm{P}$即为代码中的phid

    下面进入代码分析,首先进入createFields
    volScalarField rho
    (
        IOobject
        (
            "rho",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        rhoO + psi*p//对应方程($\ref{bar}$)
    );
    

    下面进入主程序:

    int main(int argc, char *argv[])
    {
        #include "postProcess.H"
    
        #include "setRootCase.H"
        #include "createTime.H"
        #include "createMesh.H"
        #include "createControl.H"
        #include "createFields.H"
        #include "initContinuityErrs.H"
    
        // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    
        Info<< "\nStarting time loop\n" << endl;
    
        while (runTime.loop())
        {
            Info<< "Time = " << runTime.timeName() << nl << endl;
    
            #include "compressibleCourantNo.H"
    
            solve(fvm::ddt(rho) + fvc::div(phi));
    
            // --- Pressure-velocity PIMPLE corrector loop
            while (pimple.loop())
            {
                fvVectorMatrix UEqn
                (
                    fvm::ddt(rho, U)
                  + fvm::div(phi, U)
                  - fvm::laplacian(mu, U)
                );//对应方程($\ref{Up}$)
    
                solve(UEqn == -fvc::grad(p));
    
                // --- Pressure corrector loop
                while (pimple.correct())
                {
                    volScalarField rAU("rAU", 1.0/UEqn.A());
                    surfaceScalarField rhorAUf
                    (
                        "rhorAUf",
                        fvc::interpolate(rho*rAU)
                    );
    
                    U = rAU*UEqn.H();
    
                    surfaceScalarField phid
                    (
                        "phid",
                        psi
                       *(
                           fvc::flux(U)
                         + rhorAUf*fvc::ddtCorr(rho, U, phi)/fvc::interpolate(rho)
                        )
                    );
    
                    phi = (rhoO/psi)*phid;
    
                    fvScalarMatrix pEqn
                    (
                        fvm::ddt(psi, p)
                      + fvc::div(phi)
                      + fvm::div(phid, p)
                      - fvm::laplacian(rhorAUf, p)
                    );//对应方程($\ref{NS6}$)
    
                    pEqn.solve();
    
                    phi += pEqn.flux();
    
                    solve(fvm::ddt(rho) + fvc::div(phi));
                    #include "compressibleContinuityErrs.H"
    
                    U -= rAU*fvc::grad(p);
                    U.correctBoundaryConditions();
                }
            }
    
            rho = rhoO + psi*p;
    
            runTime.write();
    
            Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
                << "  ClockTime = " << runTime.elapsedClockTime() << " s"
                << nl << endl;
        }
    
        Info<< "End\n" << endl;
    
        return 0;
    }
    

    需要注意的几下几点:
  • 代码中的phi += pEqn.flux();实际形式为phi += pEqn.flux() + fvc::interpolate(p)*phid;
  • 方程($\ref{NS6}$)中的第二项需要显性离散
  • 第四项直接采用当前的密度,这样的操作可以视为一种线性化,否则此项会变为关于压力的非线性项。因此需要在压力循环内的压力方程求解后对密度再次求解进行更新
  • sonicLiquidFoam中,由于压力波的传递,为了很好的捕获这种压力波,需要把时间步长设置的足够小。压力波在水中的传递以声速进行。水中的声速和$\psi$有关。在运行算例的时候,确保要保证依据声速计算的库朗数足够小。例如当$\psi=4e-7$的情况下,库朗数约为5e-7。在这样的时间步下,依据速度计算的库朗数会非常小。


  • 如何引用:
    李东岳. sonicLiquidFoam解析[OL]. http://dyfluid.com/sonicLiquidFoam.html, 2016.09.22

    东岳流体dyfluid.com
    勘误、讨论、补充内容请前往CFD中国