• 返回主页
  • pimpleDyMFoam解析

    Reviewed by: 2016.09.07,Update:修正体积守恒方程以及连续性方程;规范方程书写
    求解器类型:瞬态大时间步长不可压缩单相流体求解器。支持动网格。

    pimpleDyMFoam求解器的大体思想就是在时间步内采取piso算法进行时间步进,同时也可以对速度矩阵系数、非最终迭代步的压力场进行低松弛求解以提高稳定性。相对于pimpleFoampimpleDyMFoam添加了动网格支持。另外一个和pimpleFoam的区别为所有基于动网格的控制方程和非动网格的控制方程有所不同,本文首先对动网格控制方程简单介绍,然后对求解器进行分析。

    考虑一个控制质量(control mass,其可以看做可变形的控制体(control volume))内的流体,我们有质量守恒$^2$: \begin{equation} \frac{\mathrm{d} m}{\mathrm{d} t}=0 \end{equation} 其中$m$表示流体的质量。上述方程表明,控制质量内的流体的质量(不管控制体怎样变形)对时间的变化等于0。其也可以表示为$^2$: \begin{equation} m=\int_{\Omega_\mathrm{CM}} \rho \mathrm{d} \Omega \equiv \mathrm{Constant} \end{equation} 其中$\Omega_\mathrm{CM}$表示控制质量,把公式(2)代入到公式(1)并应用Reynold's transport theorem,我们有$^2$: \begin{equation} \frac{\mathrm{d}}{dt}\int_{\Omega_\mathrm{CM}} \rho \mathrm{d} \Omega = \frac{\mathrm{d}}{dt}\int_{\Omega_\mathrm{CV}} \rho \mathrm{d} \Omega +\frac{\mathrm{d}}{\mathrm{d}t}\int_{\Omega_\mathrm{CV}} \rho (\mathbf{U}-\mathbf{U}_b)\mathrm{d} \Omega=0 \end{equation} 其中$\Omega_\mathrm{CV}$表示控制体,$\mathbf{U}_b$表示控制体面上的速度。其也可以表示为$^2$: \begin{equation} \frac{\partial \rho}{\partial t} + \nabla \cdot \left(\rho \left(\mathbf{U} - \mathbf{U}_b \right) \right)=0 \end{equation} ● 对于静态网格,控制体不移动,$u_b=0$,有常见的可压缩流体连续性方程$^2$: \begin{equation} \frac{\partial \rho}{\partial t} + \nabla \cdot \left(\rho \mathbf{U} \right)=0 \end{equation} ● 对于动网格,$u_b \neq 0$,公式(4)即为动网格的可压缩流体连续性方程。进一步的如果考虑不可压缩流体,提取出密度,进一步演化为(16.09.06删除)对于不可压缩流体,由于网格运动导致的体积变化不可以简单的把密度导数项约去。同时,网格运动产生的质量通量和密度非稳态项相互抵消$^2$,这样我们有动网格下的不可压缩流体连续性方程: \begin{equation} \nabla \cdot \left(\mathbf{U} \right)=0 \end{equation} 类似的,我们有可压缩流体动网格的动量方程以及不可压缩流体的动量方程$^2$: \begin{equation} \frac{\partial \rho \mathbf{U}}{\partial t}+\nabla \cdot \left(\rho \mathbf{U}\left(\mathbf{U}-\mathbf{U}_b \right) \right)=-\nabla p+\nabla \cdot(\rho \nabla \mathbf{U}) \end{equation} \begin{equation} \frac{\partial \mathbf{U}}{\partial t}+\nabla \cdot \left(\mathbf{U}\left(\mathbf{U}-\mathbf{U}_b \right) \right)=-\frac{\nabla p}{\rho}+\nabla \cdot( \nabla \mathbf{U}) \label{mom} \end{equation} 观察不可压缩流体的控制方程(6)和($\ref{mom}$),我们发现动网格的控制方程和普通的控制方程区别仅在于对流项,因此我们求解公式(6)和(8)就可以了。但是同时还需要一个不可压缩流体的附加方程,即不可压缩体积守恒方程(space conservative low)$^{1,2}$: \begin{equation} \frac{\partial 1}{\partial t}-\nabla \cdot \mathbf{U}_b=0 \end{equation} 公式(6),(8)和(9)即为不可压缩流体涉及到网格变形的控制方程。需要注意的是$\mathbf{U}_b$只在运动的网格处不为0,在静止的网格$\mathbf{U}_b=0$。因此可以看出,动网格求解器的速度方程($\ref{mom}$)中的phi实际为相对速度组建。在组建压力泊松方程的时候,调用的为绝对速度。

    pimpleDyMFoamcreateFields.H文件分析请参考scalarTransportFoam解析。下面先介绍求解器流程,打开pimpleDyMFoam.C

    #include "fvCFD.H"
    #include "dynamicFvMesh.H"//动网格求解器需要包含此文件
    #include "singlePhaseTransportModel.H"
    #include "turbulentTransportModel.H"
    #include "pimpleControl.H"//pimple控制需要包含此文件
    #include "CorrectPhi.H"//动网格求解器需要包含此文件
    #include "fvOptions.H"
    
    // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    
    int main(int argc, char *argv[])
    {
        #include "setRootCase.H"
        #include "createTime.H"
        #include "createDynamicFvMesh.H"//动网格求解器需要包含此文件
        #include "initContinuityErrs.H"
    
        pimpleControl pimple(mesh);
    
        #include "createFields.H"
        #include "createUf.H"//从网格体心速度插值计算面心速度
        #include "createMRF.H"
        #include "createFvOptions.H"
        #include "createControls.H"
        #include "CourantNo.H"
        #include "setInitialDeltaT.H"
    
        turbulence->validate();
    
        // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
    
        Info<< "\nStarting time loop\n" << endl;
    
        while (runTime.run())
        {
            #include "readControls.H"
            #include "CourantNo.H"
    
            #include "setDeltaT.H"
    
            runTime++;
    
            Info<< "Time = " << runTime.timeName() << nl << endl;
    
            mesh.update();//这里计算网格运动,例如网格滑移旋转、扭曲变形、6DOF运动、震荡等,具体网格如何变形的,还需参考具体的网格运动具体分析。
    首先,mesh.update()调用solidBodyMotionFvMesh.update(),然后调用fvMesh中的movePoints()函数,进而调用运动求解器中的transformation()函数。
    同时fvMesh还会创建网格面通量meshPhi(mesh motion flux)。此meshphi会在makeAbsolute()以及makeRelative()函数中用到。
    // Calculate absolute flux from the mapped surface velocity phi = mesh.Sf() & Uf;//绝对通量计算,可参考《何为通量》。另外,如果是从下一个时间步计算而来,最后写入的phi为相对通量,因此在此处需要转换为绝对通量。 if (mesh.changing() && correctPhi) { #include "correctPhi.H"//公式(9)。 }

    在动网格算法中,需要进行通量修正以保证连续性。举一个简单的例子:考虑一个气球为一个控制体(网格),气球的边界即为控制体的边界,其足够的柔软因此气球边界的速度等于内部流体流动的速度。气球内部的水的流动会导致气球的变形,但是并没有水从气球内流出。因此,在这个气球内,虽然速度在流动,但是并面通量为0。在每个时间步的第一次求解压力方程之前都要进行修正。现在我们看如何对通量进行修正,截取CorrectPhi.H中的必要片段:

        while (pimple.correctNonOrthogonal())
        {
            // Solve for pcorr such that the divergence of the corrected flux
            // matches the divU provided (from previous iteration, time-step...)
            fvScalarMatrix pcorrEqn
            (
                fvm::laplacian(rAUf, pcorr) == fvc::div(phi) - divU//公式(9),然而实际上公式(9)对应于fvm::laplacian(rAUf, pcorr) == fvc::div(mesh.phi()) - divU,
    进一步讨论:点击此处
    ); pcorrEqn.setReference(0, 0); pcorrEqn.solve(); if (pimple.finalNonOrthogonalIter()) { phi -= pcorrEqn.flux();//参考icoFoam解析 } }

    上述代码即为求解公式(9),在求解之后phi即为守恒的$\mathrm{U}_f \cdot \mathrm{S}_f$场。在有了$\mathrm{U}$之后,继续看pimpleDyMFoam的代码:

            fvc::makeRelative(phi, U);//转换为相对通量。 
    
            if (mesh.changing() && checkMeshCourantNo)
            {
                #include "meshCourantNo.H"
            }
    
            // --- Pressure-velocity PIMPLE corrector loop
            while (pimple.loop())
            {
                #include "UEqn.H"//公式(8)。在这里采用相对通量计算绝对速度 
                ...............
    

    第一行,在fvc::makeRelative(phi, U)执行之后,phi则转化为相对通量$\mathbf{U}_r = \mathbf{U}-\mathbf{U}_b$,因此下文中的phi则表示$\mathbf{U}_r$。其中$\mathbf{U}$被传入meshPhi(U)函数来计算网格面通量。然后,速度方程通过相对通量phi来构建关于$\mathbf{U}$的矩阵系数。随后求得绝对速度$\mathbf{U}$。进入压力方程:

    .............
    while (pimple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
        (
            fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
        );//公式(6)
    
        pEqn.setReference(pRefCell, pRefValue);
    
        pEqn.solve(mesh.solver(p.select(pimple.finalInnerIter())));
    
        if (pimple.finalNonOrthogonalIter())
        {
            phi = phiHbyA - pEqn.flux();//phi在此处计算后为绝对通量,前文之所以要把phi转换为相对通量,是因为需要通过动量方程构建离散的速度矩阵方程组。
        }
    }
    
    #include "continuityErrs.H"
    
    // Explicitly relax pressure for momentum corrector
    p.relax();
    
    U = HbyA - rAtU*fvc::grad(p);//更新速度
    U.correctBoundaryConditions();
    fvOptions.correct(U);
    
    {
        Uf = fvc::interpolate(U);
        surfaceVectorField n(mesh.Sf()/mesh.magSf());
        Uf += n*(phi/mesh.magSf() - (n & Uf));//0? Leave your comments below
    }
    
    // Make the fluxes relative to the mesh motion
    fvc::makeRelative(phi, U);//再次将phi变为相对通量,因为湍流方程中对流项中的显性离散项也为相对速度
    

    [1]. Demirdžić I, Perić M. Finite volume method for prediction of fluid flow in arbitrarily shaped domains with moving boundaries[J]. International Journal for Numerical Methods in Fluids, 1990, 10(7): 771-790.
    [2]. Ferziger J.H., Perić M. Computational methods for fluid dynamics[M]. Springer Science & Business Media, 2002.


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

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