• 返回主页
  • scalarTransportFoam解析

    Reviewed by:
    求解器类型:稳态或瞬态标量传输求解器。

    scalarTransportFoam是OpenFOAM3个最基本的求解器之一,用于求解标量(passive scalar)传输问题。在这个求解器中,用户需要在constant/transportProperties当中设置常量DT(标量的扩散率)。由于标量的定义,其对流场不具有影响作用,因此我们没有必要计算流场。因此在0/U中提供一个不随时间变化的速度场即可。并且这个求解器不能用于求解流场随时间变化的问题。

    首先进入createFields.H:(大部分代码略,请参考laplacianFoam解析

    Info<< "Reading/calculating face flux field phi\n" << endl;
    
    surfaceScalarField phi
    (
        IOobject
        (
            "phi",
            runTime.timeName(),
            mesh,
            IOobject::READ_IF_PRESENT,
            IOobject::AUTO_WRITE
        ),
        linearInterpolate(U) & mesh.Sf()//公式(1)
    );
    

    phi是一个在OpenFOAM中非常常见的定义,即为通量(请参考:何为通量)。其在此处被定义为: \begin{equation}\nonumber \mathrm{phi}=u_f \cdot S_f \end{equation} 然后进入scalarTransportFoam.C

    #include "fvCFD.H"//在OpenFOAM的所有求解器中都可以看到这个头文件,它涉及到构建时间、组建矩阵、有限体积离散、组建网格、
    量纲设置等大量内容。
    建议在自定义的求解器中必备此项。初学可以忽略此头文件(下文中用“请忽略”来表示,意义为可以不去深究的头文件,
    来把精力放在重头)
    #include "fvOptions.H"//通过fvOption来修正源项,主要用于在运行时添加多孔介质、MRF多重参考系、隐形显性热源等,
    如果自定义求解器不需要源项,可以删去。请忽略
    #include "simpleControl.H"//定义SIMPLE循环,使用SIMPLE循环必须定义此文件。请忽略 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // int main(int argc, char *argv[]) { #include "setRootCase.H"//设置算例的根目录,必备头文件。请忽略 #include "createTime.H"//创建时间对象,大部分求解器都需要此文件。请忽略 #include "createMesh.H"//创建网格对象,必备头文件。请忽略 simpleControl simple(mesh);//对于采用SIMPLE算法的算例,必备此项,请忽略 #include "createFields.H"//包含上文分析过的createFields.H头文件 #include "createFvOptions.H"//创建源项,无需源项可删除。请忽略 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * // Info<< "\nCalculating scalar transport\n" << endl; #include "CourantNo.H" //计算库郎数,见下文分析 while (simple.loop()) //开始SIMPLE循环,采用SIMPLE算法必备语句 { Info<< "Time = " << runTime.timeName() << nl << endl; while (simple.correctNonOrthogonal()) { fvScalarMatrix TEqn ( fvm::ddt(T)//时间项,公式(1)左边第一项 + fvm::div(phi, T)//对流项,公式(1)左边第二项 - fvm::laplacian(DT, T)//扩散项,公式(1)左边第三项 == fvOptions(T)//源项,公式(1)右边 );//组建TEqn,公式(1) TEqn.relax();//对上述方程松弛 fvOptions.constrain(TEqn);//对方程系数进行源项限定 TEqn.solve();//对上述方程求解 fvOptions.correct(T);//对T进行源项修正 } runTime.write(); } Info<< "End\n" << endl; return 0; }

    现在我们分析scalarTransportFoam中植入的方程: \begin{equation} \frac{\partial T}{\partial t}+\nabla \cdot (uT)-\nabla \cdot(\nabla D_\mathrm{T} T) = S \end{equation} 参考上述代码,所有$T$均为隐形离散。

    最后我们分析Courant.H

    scalar CoNum = 0.0;//定义一个scalar(类似C++中的double),CoNum=0.0
    scalar meanCoNum = 0.0;//同上
    
    if (mesh.nInternalFaces())//如果是网格的内部面,不考虑边界面
    {
        scalarField sumPhi
        (
            fvc::surfaceSum(mag(phi))().internalField()//公式(2)
        );//scalarField在OpenFOAM中为一系列scalar组成的场,其区别于volScalarField主要为不具有边界条件等信息
    
        CoNum = 0.5*gMax(sumPhi/mesh.V().field())*runTime.deltaTValue();//公式(3)
    
        meanCoNum =
            0.5*(gSum(sumPhi)/gSum(mesh.V().field()))*runTime.deltaTValue();//公式(4)
    }
    
    Info<< "Courant Number mean: " << meanCoNum
        << " max: " << CoNum << endl;
    

    CFD计算中通常要求Courant数小于1,在某些多相流情况下小于0.5是最好的。然而某些特定的数值格式可以使用较大的Courant数。针对上述代码,我们有下述公式: \begin{equation} \mathrm{sumPhi}=\sum (|u_f \cdot S_f|) \end{equation} \begin{equation} \mathrm{CoNum}=0.5 \mathrm{max} \left( \frac{\mathrm{sumPhi}}{V} \right)\cdot \Delta t=0.5 \mathrm{max} \left( \frac{\sum (|u_f \cdot S_f|)}{V} \right)\cdot \Delta t \end{equation} 其中$\mathrm{CoNum}$即为Courant数。众所周知,文献中的Courant数定义为 \begin{equation} \mathrm{CoNum}= \frac{(|u_f|)}{\Delta x} \cdot \Delta t \end{equation} 仔细观察公式(3),(4)我们会发现OpenFOAM的Courant数被除了2,这是因为OpenFOAM对通量进行了加和来计算$\mathrm{sumPhi}$。比如,对于六面体网格的对立面,进入网格和出去网格的通量计算了2次,因此需要除2。

    如何引用:

    李东岳, 高翔宇. scalarTransportFoam解析.[OL]. http://dyfluid.com/scalarTransportFoam.html, 2016.3.10

    东岳流体
    info@dyfluid.com