OpenFOAM小代码


  • 管理员

    这些是我感觉对初学者有用的一些代码,收藏在此,持续更新。有问题回帖提问,我会对我贡献的代码负责回答。

    计算Patch的面积

        // 计算某个patch的面积 http://www.cfd-china.com/topic/3498
        {
            label patchID = mesh.boundaryMesh().findPatchID("inlet");
    
            // Create a polyPatch for looping
            const polyPatch& myPatch = mesh.boundaryMesh()[patchID];
            
            // Initialize patchArea
            scalar patchArea = 0.0;
            
            // Loop trhough all faces on the polyPatch, adding their magnitude surface
            // area vectors
            forAll(myPatch, faceI)
            {
                patchArea += mesh.magSf().boundaryField()[patchID][faceI];
            } 
        }
    

    获得Patch的名字

    Info << "patch name is " << U.boundaryField()[patch].patch().name() << nl;
    

    网格相关量

    const scalarField& V = mesh.V();                            // 网格体积
    const surfaceVectorField& Sf = mesh.Sf();                   // 网格面矢量
    const surfaceScalarField& magSf = mesh.magSf();             // 网格面积的模
    const surfaceVectorField& N = Sf/magSf;                     // 网格面法向
    

    更改壁面一层网格值

            label patchID = mesh.boundaryMesh().findPatchID("wall");
            const scalarField TfaceCell = 
                T.boundaryField()[patchID].patchInternalField();
    
            k.boundaryFieldRef()[patchID] == sqrt(TfaceCell);
    

    给定patch调用网格

    如果p是一个fvPatch

    const fvMesh& mesh = p.boundaryMesh().mesh();
    

    类内调用时间步长

    假定mesh_是类内可以调用的成员变量。

    const scalar deltaT = mesh_.time().deltaT().value();
    

    如果mesh_不可获取,U_可获取:

    const scalar deltaT = U_.mesh().time().deltaT().value();
    

    声明一个tmp变量

    tmp<volScalarField> deltaLambdaT   
    (
        volScalarField::New
        (
            "deltaLambdaT",
            mesh_,
            dimensionedScalar(dimless, 1.0)
        )   
    );   
    

    简写边界场

    volScalarField::Boundary& meanUb = meanU_.boundaryFieldRef();
    
    forAll(meanUb, patchi)
    {
        fvPatchScalarField test = meanUb[patchi];
    }
    

    判断是否是固定值边界、processor边界

        forAll(mesh_.boundary(), P)                                              
        {                                                                        
            if
            (
                isA<fixedValueFvPatchScalarField>
                (
                    T.boundaryField()[P]
                ) 
             ||
                isA<processorPolyPatch>(mesh_.boundaryMesh()[P])
            )
            {
            }
        }
    
    

    codedFixedValue

    wall
        {
            type            codedFixedValue;
            value           uniform 300; //default value
            redirectType    linearTBC; //name of new BC type
            code
            #{
            const vectorField& Cf = patch().Cf(); // get face center coordinate;
            //assume the wall is a vertical wall       
            scalar ymax = max(Cf&vector(0,1,0)); // `&` is dot product
            scalar ymin = min(Cf&vector(0,1,0));
            // Info<<"ymax="<<ymax<<",ymin="<<ymin<<nl;
            
            vectorField& vf = *this; // get the boundary field (List of vectors defined at face centers) to fill.
            //temporal coordinate, type: scalar
            scalar t =this->db().time().value(); // get time
            // temporal term
            // scalar tt = 1+0.1*sin(5*t);
            scalar tt = 1.0;
            forAll(Cf,faceI)
            {
            	//spatial coordinates, type: scalar
            	//scalar x = Cf[faceI].x();
            	scalar y = Cf[faceI].y(); 
            	//scalar z = Cf[faceI].z();
            
            	vf[faceI] = ((y-ymin)/(ymax-ymin)*500+300)*tt; 
            }
            #};
            
            //I do not know why I need to add those things
            //codeInclude
            //#{
            //    #include "fvCFD.H"
            //#};
            
            //codeOptions
            //#{
            //    -I$(LIB_SRC)/finiteVolume/lnInclude
            //#};
        }
    

    努塞尔数

    下列代码可以添加到controlDict中执行

    
    functions
    {
        coded
        {
            functionObjectLibs ( "libutilityFunctionObjects.so" );
            enabled         true;
            type            coded;
            redirectType    printMinU;
            writeControl   timeStep;
            writeInterval  1;
    
            codeOptions
            #{
                -I$(LIB_SRC)/meshTools/lnInclude
            #};
    
            codeExecute
            #{
                const volScalarField& T
                (
                    mesh().lookupObject<volScalarField>("T")
                );
    
                const fvPatchList& patches = mesh().boundary();
    
                forAll(patches, patchi)
                {
                    const fvPatch& currPatch = patches[patchi];
    
                    if (mesh().time().value() == 30001)
                    {
                        if (currPatch.name() == "downwall")
                        {
                            fvPatchScalarField nus = T.boundaryField()[patchi];
    
                            scalar L = 0.014231;
                            scalar Tref = 309.1;
                            scalar Timp = 347.6;
                            
                            nus = T.boundaryField()[patchi].snGrad()*L/(Timp - Tref);
    
                            forAll(T.boundaryField()[patchi], facei)
                            {
                                Pout << mesh().C().boundaryField()[patchi][facei].x()/0.014231
                                     << " " << nus[facei] << nl;
                            }
                        }
                    }
                } 
           #};
        }
    }
    

    获取某个函数的计算时间

    #include <chrono>
    
    auto start = std::chrono::steady_clock::now();                                                   
    // Functions here
    auto end = std::chrono::steady_clock::now();                                                     
    auto diff = end - start;
    Info<< "Calculate nodes and weights, using " 
        << std::chrono::duration <double, std::milli> (diff).count()                                 
        << " ms" << endl;
    

    在某个文件夹下输出scalarField

    by @Samuel-Tu

        IOField<scalar> utau
        (
            IOobject
            (
                "utau",
                 runTime.constant(),
                "../postProcessing",
                mesh,
                IOobject::NO_READ,
                IOobject::AUTO_WRITE 
            ),
            scalarField(totalFSize,0.0)
        );
    

    BlockMesh-simpleGrading

    处理非均一网格

    blocks
    (
        hex (0 1 2 3 4 5 6 7) (100 300 100)
        simpleGrading
        (
            1                  // x-direction expansion ratio
            (
                (0.2 0.3 4)    // 20% y-dir, 30% cells, expansion = 4
                (0.6 0.4 1)    // 60% y-dir, 40% cells, expansion = 1
                (0.2 0.3 0.25) // 20% y-dir, 30% cells, expansion = 0.25 (1/4)
            )
            3                  // z-direction expansion ratio
        )
    );
    

    输出yplus

    functions
    {
           yplus
          {
             type                   yPlus;	 
             functionObjectLibs      ("libutilityFunctionObjects.so");
             outputControl           outputTime;
             enabled                 true;
           }
    }
    

    by @Sloan

    在流场中添加颗粒失踪

    controlDict添加下述内容,同时需要在constant文件夹添加kinematicProperties字典文件

    functions
    {
        tracks
        {
            libs        ("liblagrangianFunctionObjects.so");
            type        icoUncoupledKinematicCloud;
            U           U.water;//速度场名字
        }
    }
    
    

    在流场中添加欧拉标量传输

    controlDict添加下述内容,同时需要在0文件夹下添加s标量场,以及fvScheme、fvSolution都要做适配

    functions
    {
        s
        {
            type            scalarTransport;
            libs            ("libsolverFunctionObjects.so");
            resetOnStartUp  no;
            field           s;
            phi             phi.water;
        }
    }
    
    

    对网格进行赋值

    forAll(T, C)//遍历网格,T为场
    {
        T[C] = min(nu[C], 100);//对T的第C个网格,赋值第C个网格的nu与100的最小值
    }
    

    声明多个场PtrList

    PtrList<surfaceScalarField> abPosCoeff(4);                                                           
                                                                                                         
    forAll(abPosCoeff, i)
    {                                                                                                    
        abPosCoeff.set
        (                                                                                                
            i,                                                                                           
            new surfaceScalarField                                                                       
            (                                                                                            
                IOobject                                                                                 
                (                                                                                        
                    "abPosCoeff" + name(i),                                                              
                    mesh.time().timeName(),                                                              
                    mesh,                                                                                
                    IOobject::NO_READ,                                                                   
                    IOobject::NO_WRITE                                                                   
                ),                                                                                       
                mesh,                                                                                    
                dimensionedScalar(inv(dimLength), 0.0)                                                   
            )                                                                                            
        );                                                                                               
    } 
    

    OpenFOAM的二维数组

        PtrList<PtrList<scalar>> te(4);
                                
        forAll(te, i)
        {            
            te.set(i, new PtrList<scalar>(4));
                                  
            forAll(te[i], j)   
            {             
                te[i].set(j, new scalar (0.0));
            }                   
        } 
    

    输出最大值最小值

    functions
    {
        minMaxp
        {
            type        fieldMinMax;
            functionObjectLibs ("libfieldFunctionObjects.so");
            fields
            (
                 U
            );
            location        no;
            writeControl    timeStep;
            writeInterval   1;
        }
    }
    

  • 讲师

    补充一个CodeStream实现充分发展的湍流入口,请东岳老师不要建议
    $\begin{array}{l}
    \delta_{T}>\delta_{w} \rightarrow u=u_{\infty}\left(\frac{d_{w}}{\delta_{T}}\right)^{1 / 7} \
    \delta_{T}<\delta_{w} \rightarrow u=u_{\infty}
    \end{array}
    \delta_{T}=0.37\left(v / u_{\infty}\right)^{0.2} L_{T}^{0.8}
    $
    dw为到壁面的最短距离

    OpenFoam 7:

    gasinlet
    {
        type            fixedValue;
        value           #codeStream
        {
            codeInclude
            #{
                #include "fvCFD.H"
                #include "DynamicList.H"
            #};
    
            codeOptions
            #{
                -I$(LIB_SRC)/finiteVolume/lnInclude \
                -I$(LIB_SRC)/meshTools/lnInclude
            #};
    
    	    //libs needed to visualize BC in paraview
    	    codeLibs
    	    #{
        	    -lmeshTools \
        	    -lfiniteVolume
    	    #};
    
            code
            #{
                const IOdictionary& d = static_cast<const IOdictionary&>
    			(
                    dict.parent().parent()
                );
                const fvMesh& mesh = refCast<const fvMesh>(d.db());
                const label id = mesh.boundary().findPatchID("gasinlet");
                const fvPatch& patch = mesh.boundary()[id];
    
                vectorField U(patch.size(), vector(0, 0, 0));
    
                const scalar U_inf=93.5;
                const scalar LT=0.2;
                const scalar visco=2.3183558929634845e-06;
                const scalar deltaT=0.37*pow(visco/U_inf,0.2)*pow(LT,0.8);
                scalar dw=0.0;
                forAll(U, i)
                {
                    const scalar x = patch.Cf()[i][0];
                    const scalar y = patch.Cf()[i][1];
                    const scalar z = patch.Cf()[i][2];
                    DynamicList<scalar>  distance(4);
                    distance.append(mag(y));
                    distance.append(mag(y-30e-6));
                    distance.append(mag(z-23e-6));
                    distance.append(mag(z+23e-6));
                    dw=min(distance);
                    distance.clearStorage();
                    scalar Ux=0.0;
                    if(deltaT>dw)
                    {
                        Ux=U_inf*pow(dw/deltaT,1.0/7.0);
                    }
                    else
                    {
                        Ux=U_inf;
                    }
    	            U[i] = vector(Ux, 0, 0);
                }
    
                writeEntry(os, "", U);
            #};
        };
    }


  • @yfclark 谢谢贡献代码。有一个疑问还请教,这里没有脉动量,只有平均速度的径向分布。在LES或者DNS中,脉动量比较重要,所以关于脉动量的话有什么看法么?还是说这个入口可以作为一个入口,然后在一个长管道内比较容易自然发展成充分发展的湍流?代码里的U_inf可以理解为平均速度么?


  • 副教授

    wordHashSet unMatchedkeys(dict.toc());
    forAll(names, nameI)
    {
        const word& name = names[nameI];
        const entry* ePtr = dict.findEntry(name, keyType::REGEX);
        if(ePtr)
           unMatchedKeys.erase(ePtr->keyWord());
    }
    

  • 副教授

    class base
    {
    private:
        bool ownedByRegistry_;
    public:
        template<class Type>
        static void store(Type* p)
        {
    	//p->base::ownedByRegistry_ =true;
    	p->ownedByRegistry_ = true;
    	p->printInfo();
        }
    };
    
    class derived
    :public base
    {
    public:
        derived() {}
        void printInfo()
        {
    	Info << "derived class" << endl;
        }
    };
    
    int main(int argc, char *argv[])
    {
    ...
        derived d;
        base::store(&d);
    ...
    }
    


  • @东岳OpenFOAM小代码 中说:

    东岳老师,请教几个小问题:

    1. 在编译过程中报错:'mesh_' was not declared,所以mesh_是不可以直接使用的,但是我不太理解什么是类内可以调用的成员变量,感觉速度场U,温度场T等都是类内可以调用的成员变量?
    2. 如果mesh_不可获取,U_可获取。上面说mesh_是类内可以调用的成员变量,然而U_又是什么?
    3. 声明一个tmp变量,mesh_表示什么意思呢?在声明U、T、p等场的时候,mesh_的位置常常写的是mesh

    期待回复,祝好!


  • 管理员

    @D-Benjamin

    我不太理解什么是类内可以调用的成员变量

    我说的这个东西,就是类的成员变量。如果存在成员变量mesh, 它在OpenFOAM类内都被定义为mesh_,在类外那肯定是mesh,后面有没有_只是写法风格,OpenFOAM类里面的变量都带_,类外面的变量都没有_

    OpenFOAM这面类继承会涉及到好几层。你可以一层一层的找有没有mesh_的声明,你可以直接用mesh_看看是否报错,没报错就是声明了,报错了就是没生命。



  • @Zhong-combustion 你好,我目前也在用LES做加热管流的湍流研究,关于湍流入口的设置看了一些资料,如有兴趣,可以加个qq:1191364394,共同学习探讨。



  • 东岳老师您好,请问如果是两相流中相分数(比如alpha.particle)随时间变化,它在codeFixedValue中应该怎么定义呢?


  • 管理员

    @Zhy2022 codedFixedValue上面不是有嘛,你举一反三一下咩



  • @东岳 东岳老师,没能找出报错原因:135: ,只能再来麻烦一下您:

    INLET1
        {
            type            codedFixedValue;
            value           $internalField;
            name            stepInjection;
            redictType      variedValue;
            code
            #{
                 scalar alphapar = this->patch().fluid().phase();
    
                 scalar t = this->db().time().value();
           
                 if(t >= 0 && t<=5)
                 {
                      alphapar = 0;
                 }
                 
                 else
    	     {
    	          alphapar = 0.15;
    	     }
    
    	     operator==(alphapar)
    
             #};
            
        }
    

    56e49c26-63f1-479c-bd8e-d479760394d9-image.png



  • @东岳 东岳老师,问题已解决!无法编译与文件路径有关
    ee3e6b1d-91b1-4633-879c-41502c56cf20-image.png



  • fvOption小工具

    根据GitHub上一个python脚本改写了一个生成fvOption的小工具 https://github.com/Veenxz/fvSchemes
    用Foam没有多久,各位大佬可以多提提意见。

    # version 2.0
    # Based on https://github.com/Carlopasquinucci/fvSchemes
    # Generation by Veenxz
    # relaesed under license GPL GNU 3.0
    
    steady = True
    pseudo_transient = False
    precision = 2    # First order 1 or Second order 2
    unbounded = False
    LUST = False
    secondorder = True
    
    maxOrtho = 80
    maxSkew = 20
    
    # Header and Footer
    h = [
        "/*--------------------------------*- C++ -*----------------------------------*|",
        "| =========                 |                                                 |",
        "| \\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |",
        "|  \\\    /   O peration     | Website:  https://openfoam.org                  |",
        "|   \\\  /    A nd           | Version:  7                                     |",
        "|    \\\/     M anipulation  |                                                 |",
        "\*---------------------------------------------------------------------------*/",
        "FoamFile", "{", "    version     2.0;", "    format      ascii;",
        "    class       dictionary;", '    location    "system";', "    object      fvSchemes;",
        "}",
        "// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //",
        ""
    ]
    footer = '\n// ************************************************************************* //'
    
    nonOrthogonalCorrectors = 1
    if (maxOrtho) > 80:
        print(
            'Warning: mesh is not so nice. Use 3 nonOrthogonalCorrectors in fvSolution file'
        )
        nonOrthogonalCorrectors = 3
    if (maxSkew) > 8:
        print('Warning: mesh is not so nice')
    
    if (maxOrtho) > 80:
    
        gradSchemes = (
            '{\n    default          cellLimited Gauss linear 0.5;\n'
               '    grad(U)          faceLimited Gauss linear 1.0;\n}'
        )
        divSchemes = (
            '{\n    div(phi,U)       Gauss linearUpwind grad(U);\n'
               '    div(phi,omega)   Gauss upwind;\n'
               '    div(phi,k)       Gauss upwind;\n'
               '    div(phi,e)       Gauss upwind;\n'
               '    div((nuEff*dev(T(grad(U))))) Gauss linear;\n}'
        )
        laplacianSchemes = (
            '{\n    default          Gauss linear limited 0.333;\n}'
        )
        snGradSchemes = (
            '{\n    default          Gauss linear limited 0.333;\n}'
        )
    
        blending = 0.2
        nonOrthogonalCorrectors = 3
    
    if (maxOrtho) > 70:
    
        gradSchemes = (
            '{\n    default          cellLimited Gauss linear 0.5;\n'
               '    grad(U)          cellLimited Gauss linear 1.0;\n}'
        )
        divSchemes = (
            '{\n    div(phi,U)       Gauss linearUpwind grad(U);\n'
               '    div(phi,omega)   Gauss linearUpwind grad(omega);\n'
               '    div(phi,k)       Gauss linearUpwind grad(k);\n'
               '    div(phi,e)       Gauss linearUpwind grad(e);\n'
               '    div((nuEff*dev(T(grad(U))))) Gauss linear;\n}'
        )
        laplacianSchemes = (
            '{\n    default          Gauss linear limited 0.5;\n}'
        )
        snGradSchemes = (
            '{\n    default          Gauss linear limited 0.5;\n}'
        )
    
        blending = 0.5
        nonOrthogonalCorrectors = 3
    
    if (maxOrtho) > 60:
    
        gradSchemes = (
            '{\n    default          cellMDLimited Gauss linear 0.5;\n'
               '    grad(U)          cellMDLimited Gauss linear 0.5;\n}'
        )
        divSchemes = (
            '{\n    div(phi,U)       Gauss linearUpwind grad(U);\n'
               '    div(phi,omega)   Gauss linearUpwind grad(omega);\n'
               '    div(phi,k)       Gauss linearUpwind grad(k);\n'
               '    div(phi,e)       Gauss linearUpwind grad(e);\n'
               '    div((nuEff*dev(T(grad(U))))) Gauss linear;\n}'
        )
        laplacianSchemes = (
            '{\n    default          Gauss linear limited 0.777;\n} '
        )
        snGradSchemes = (
            '{\n    default          Gauss linear limited 0.777;\n} '
        )
    
        blending = 0.7
        nonOrthogonalCorrectors = 2
    
    if (maxOrtho) > 0:
    
        gradSchemes = (
            '{\n    default          cellMDLimited Gauss linear 0;\n'
               '    grad(U)          cellMDLimited Gauss linear 0.333;\n}'
        )
        divSchemes = (
            '{\n    div(phi,U)       Gauss linearUpwind grad(U);\n'
               '    div(phi,omega)   Gauss linearUpwind grad(omega);\n'
               '    div(phi,k)       Gauss linearUpwind grad(k);\n'
               '    div(phi,e)       Gauss linearUpwind grad(e);\n'
               '    div((nuEff*dev(T(grad(U))))) Gauss linear;\n}'
        )
        laplacianSchemes = (
            '{\n    default          Gauss linear limited 0.95;\n}'
        )
        snGradSchemes = (
            '{\n    default          Gauss linear limited 0.95;\n}'
        )
    
        blending = 0.8
        nonOrthogonalCorrectors = 1
    
    if (LUST):
        gradSchemes = (
            '{\n    default          Gauss linear;\n'
               '    grad(U)          cellMDLimited leastSquares 1;\n}'
        )
        divSchemes = (
            '{\n    div(phi,U)       Gauss LUST grad(U);\n'
               '    div(phi,omega)   Gauss LUST grad(omega);\n'
               '    div(phi,k)       Gauss LUST grad(k);\n'
               '    div(phi,e)       Gauss LUST grad(e);\n'
               '    div((nuEff*dev(T(grad(U))))) Gauss linear;\n}'
        )
        laplacianSchemes = (
            '{\n    default          Gauss linear corrected;\n}'
        )
        snGradSchemes = (
            '{\n    default          Gauss linear corrected;\n}'
        )
        
        blending = 0.9
        nonOrthogonalCorrectors = 1
    
    if (steady):
        ddtSchemes = (
            '{\n    default          steadyState;\n}'
        )
        if (pseudo_transient):
                    ddtSchemes = (
                        '{\n    default          localEuler;\n}'
                    )
    else:
        ddtSchemes = (
            '{\n    default           CrankNicolson ' + str(blending) + ' ;\n}'
        )
        if precision == 1:
            ddtSchemes = (
                '{\n    default           Euler;\n}'
            )
        if (unbounded):
            ddtSchemes = (
                '{\n    default           backward;\n}'
            )
    
    
    wallDist = (
        "{\n    method           meshWave;\n}"
    )
    
    #open fvSchemes and write inside
    
    f = open("fvSchemes", "w")
    
    for i in h:
    	f.write(i + "\n")
    
    f.write("ddtSchemes" + "\n")
    f.write(ddtSchemes + "\n")
    f.write("\n")
    f.write("gradSchemes" + "\n")
    f.write(gradSchemes + "\n")
    f.write("\n")
    f.write("divSchemes" + "\n")
    f.write(divSchemes + "\n")
    f.write("\n")
    f.write("laplacianSchemes" + "\n")
    f.write(laplacianSchemes + "\n")
    f.write("\n")
    f.write("snGradSchemes" + "\n")
    f.write(snGradSchemes + "\n")
    f.write("\n")
    f.write("wallDist" + "\n")
    f.write(wallDist + "\n")
    f.write(footer)
    
    f.close()
    
    print('File fvSchemes created')
    
    

  • 管理员

    @veen 厉害厉害 :146: :146:



  • @veen 这个代码还有些格式没法用……我回头测试了再放到github大家看仓库里边的代码就好了。


Log in to reply
 


CFD中文网 | 东岳流体学术 | 东岳流体商业 | 吉ICP备20003622号-1