OpenFOAM代码汇总

下列代码需要在算例0文件夹中的某个场中进行植入

边界条件中定义随位置、时间变化的进口

//随位置变化
inlet
{
    type            codedFixedValue;
    name            dummy;
    value           uniform (0 0 0);
        
    redirectType    inletLaminarSquareProfile;
    code
    #{
         const vectorField& Cf = patch().Cf(); 

         const scalar H	= 0.0052; 
         const scalar Umax = 0.879542893; 

         forAll(Cf, faceI) // loop over all the patch faces
         {
            const scalar y = Cf[faceI].y() - 0.0049; 
            (*this)[faceI] = vector(Umax*(4*y/H-4*sqr(y/H)), 0, 0);
         }
     #};
}
 
//随时间变化
inlet
{
    type            uniformFixedValue;
    uniformValue    coded;

    name            pulse;

    codeInclude
    #{
        #include "mathematicalConstants.H"
    #};

    code
    #{
        return vector
        (
            0,
            0.5*(1 - cos(constant::mathematical::twoPi*x)),
            0
        );
    #};
}

边界条件中调用fvc、mesh、字典文件、时间步长

INLET
{
    type            codedFixedValue;
    value           uniform (10 0 0); 
    name            linearTBC1; 
    codeInclude
    #{
        #include "fvCFD.H"
    #};
    codeOptions
    #{
        -I$(LIB_SRC)/finiteVolume/lnInclude \
        -I$(LIB_SRC)/meshTools/lnInclude
    #};
    codeLibs
    #{
        -lmeshTools \
        -lfiniteVolume
    #};
    code
    #{
        const fvMesh& mesh = this->patch().boundaryMesh().mesh();\\调用mesh
        dictionary C = mesh.lookupObject<dictionary>("physicalProperties");
        scalar test(readScalar(C.lookup("test")));\\调用字典文件
        const fvPatch& inletPatch = this->patch();
        const scalar deltaT = mesh.time().deltaT().value();\\调用时间步长
        const volVectorField& U = mesh.lookupObject<volVectorField>("U");
        volTensorField gradU = fvc::grad(U);\\调用fvc
    #};
}

边界条件中进行计算、更改场的值

INLET
{
    type            codedFixedValue;
    value           uniform (10 0 0); 
    name            linearTBC1; 
    code
    #{
        const fvMesh& mesh = this->patch().boundaryMesh().mesh();
        // 计算inlet的面积
        label patchID = mesh.boundaryMesh().findPatchID("inlet");
        const polyPatch& myPatch = mesh.boundaryMesh()[patchID];
        scalar patchArea = 0.0;
        forAll(myPatch, faceI)
        {
            patchArea += mesh.magSf().boundaryField()[patchID][faceI];
        } 
        
        //更改U的值
        volVectorField& U = mesh.lookupObjectRef<volVectorField>("U");
    #};
}

下列代码需要在程序中植入

网格相关量、zone网格、常用CFD变量

const scalarField& V = mesh.V();                            // 网格体积

const surfaceVectorField& Sf = mesh.Sf();                   // 网格面矢量

const surfaceVectorField& Cf = mesh.Cf();                   // 网格面心位置

const surfaceScalarField& magSf = mesh.magSf();             // 网格面积的模

const surfaceVectorField& N = Sf/magSf;                     // 网格面法向量

const label& nCells = mesh.nCells();                        // 网格单元数量

const label& nInternalFaces = mesh.nInternalFaces();        // 网格内部面数量

const meshCellZones& cellZones = mesh.cellZones();          // 网格cellZone编号

wordList zoneNames = mesh.cellZones().names();              // 网格cellZone的名字

label zoneI =  mesh.cellZones().whichZone(celli);           // 判断celli在哪个cellZone里

label i = cellZones[0][i];                                  // 网格cellZone[0]的真实网格编号

const fvPatchList& patches = mesh.boundary();               // 边界网格，是一系列patch的集合

const polyBoundaryMesh& boundaryMesh = mesh.boundaryMesh(); // 边界网格，是一系列patch的集合

const label patchID = boundaryMesh.findPatchID("inlet");    // inlet边界的ID

volScalarField S(magSqr(symm(tgradU())));                   // 形变率的双点积
tmp<volTensorField> tgradU = fvc::grad(U);

声明字典文件、tmp体心场面心场、带单位的值、不带边界的体心场、壁面距离、场数组

IOdictionary physicalProperties
(
    IOobject
    (
        "physicalProperties",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
    )
);

bool training = readBool(physicalProperties.lookup("training"));

volScalarField T//体心场
(
    IOobject
    (
        "T",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

tmp<volScalarField> deltaLambdaT//tmp体心场
(
    volScalarField::New
    (
        "deltaLambdaT",
        mesh,
        dimensionedScalar(dimless, 1.0)
    )   
);

surfaceScalarField phi//面心场
(
    IOobject
    (
        "phi",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar(dimless, 0.0)
);

dimensionedScalar DT//带单位的值
(
    "DT",
    dimensionSet(0,0,0,0,0),
    1.0
);

// 一个scalarField
volScalarField::Internal TInternal("TInternal", T);//不带边界的体心场

#include "wallDist.H"//壁面距离需要包括此头文件
volScalarField y(wallDist::New(mesh).y());

PtrList<volScalarField> nutlist(BATCH);//场数组
forAll(nutlist, i)
{      
    nutlist.set
    ( 
        i, 
        new volScalarField
        (                
            IOobject    
            (          
                "nutlist" + name(i), 
                mesh.time().timeName(),
                mesh,                 
                IOobject::MUST_READ,
                IOobject::NO_WRITE   
            ),                      
            mesh
        )                          
    );                            
} 

修改内部场、修改边界值

U.primitiveFieldRef() = ...
U.ref() = ...

U.boundaryFieldRef() = ...

禁用log输出

SolverPerformance<scalar>::debug = 0;

创建2套网格

fvMesh mesh1
(
    IOobject
    (
        "mesh1",
        runTime.timeName(),
        runTime,
        IOobject::MUST_READ
    )
);

fvMesh mesh2
(
    IOobject
    (
        "mesh2",
        runTime.timeName(),
        runTime,
        IOobject::MUST_READ
    )
);

读取scalar的值、读取bool的值、寻找现存的某个场

scalar RR(readScalar(generalProperties.lookup("RR")));

bool RR(readBool(generalProperties.lookup("RR")));

const volScalarField& T(mesh().lookupObject<volScalarField>("T"));

更改变量近壁网格体心值

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

判断wall边界、fixedValue边界、processor边界

#include "wallFvPatch.H"

forAll(mesh.boundary(), patchi)
{
    if (isA<wallFvPatch>(mesh.boundary()[patchi]))
    {

    }
    
    if
    (
        isA<fixedValueFvPatchScalarField>
        (
            T.boundaryField()[patchi]
        ) 
    )
    {

    }
    
    if
    (
        isA<processorPolyPatch>
        (
            mesh.boundaryMesn()[patchi]
        ) 
    )
    {

    }
}

把A的边界条件类型赋给B、强制边界条件类型

volScalarField B
(
    IOobject
    (
        "B",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar(dimless, 0.0),
    A.boundaryField().types()
);

volScalarField B
(
    IOobject
    (
        "B",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar(dimless, 0.0),
    zeroGradientFvPatchScalarField::typeName
);

强制修改const变量

volScalarField& alphat =
    const_cast<volScalarField&>(mesh().lookupObject<volScalarField>("alphat"));  

声明一个动态调整的数组

List<label> markedCell;

for (int i = 0; i < 10; i++)
{
    markedCell.append(i);
}

矩阵操作

scalarSquareMatrix squareMatrix(3, Zero);
squareMatrix(0, 0) = 4;
squareMatrix(0, 1) = 12;
squareMatrix(0, 2) = -16;
squareMatrix(1, 0) = 12;
squareMatrix(1, 1) = 37;
squareMatrix(1, 2) = -43;
squareMatrix(2, 0) = -16;
squareMatrix(2, 1) = -43;
squareMatrix(2, 2) = 98;
const scalarField S(3, 1);

LUscalarMatrix L(squareMatrix);

scalarSquareMatrix inv(3); //矩阵求逆
L.inv(inv);

scalarField x(3, Zero);
scalarField Mx = squareMatrix*x; //矩阵乘以向量

scalarField x(L.solve(S)); //计算L \cdot x = S

监控代码计算时间

#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

IOField<scalar> utau
(
    IOobject
    (
        "utau",
         runTime.constant(),
        "../postProcessing",//输出到postProcessing文件夹中
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE 
    ),
    scalarField(totalFSize,0.0)
);

对某个变量做界限

bound(S, dimensionedScalar(S.dimensions(), 0));

显性、隐性离散

fvm::ddt(T)            \\隐性时间项离散
fvm::laplacian(T)      \\隐性拉普拉斯项离散
fvm::div(phi, T)       \\隐性对流项离散
fvm::Sp(coeff, T)      \\隐性源项离散
fvm::SuSp(coeff, T)    \\依据coeff的符号进行隐性或显性源项离散
fvc::ddt(T)            \\显性时间项离散
fvc::laplacian(T)      \\显性拉普拉斯项离散
fvc::div(phi, T)       \\显性对流项离散
fvc::grad(T)           \\显性梯度项离散

下面的代码需要放到算例文件的controlDict中执行：

输出时间平均、雷诺应力、空气龄、舒适度、缓存场、热通量、多相界面、界面高度

cacheTemporaryObjects//输出缓存场需要用
(
    kEpsilon:G
);

libs//输出多相界面、界面高度需要用
(
    "libwaves.so"
);

functions
{
    #includeFunc fieldAverage(U, p, prime2Mean = yes)\\场的时间平均
    
    #includeFunc mag(UPrime2Mean)
    #includeFunc multiply(half, mag(UPrime2Mean), result = k)
    
    age
    {
        libs            ("libfieldFunctionObjects.so");
        type            age;

        diffusion       on;

        writeControl    writeTime;
        executeControl  writeTime;
    }
    
    comfort
    {
        libs            ("libfieldFunctionObjects.so");
        type            comfort;

        clothing        0.5;
        metabolicRate   1.2;
        extWork         0;
        relHumidity     60;

        writeControl    writeTime;
        executeControl  writeTime;
    }
    

    #includeFunc writeObjects(kEpsilon:G)
    
    wallHeatFlux1
    {
        type        wallHeatFlux;
        libs        ("libfieldFunctionObjects.so");
        region      fluid;
        patches     (".*Wall");
    }
    
    #includeFunc isoSurface(isoField=alpha.water, isoValue=0.5, fields=())
    
    interfaceHeight1
    {
        type            interfaceHeight;
        libs            ("libfieldFunctionObjects.so");
        locations       ((300 0 0) (450 0 0) (600 0 0));
        alpha           alpha.water;
    }
}

输出速度分量、自定义变量、自定义时间步长、马赫数、yPlus、涡量、示踪颗粒

functions
{
    #includeFunc components(U)
    
    setDeltaT
    {
        type coded;
        libs            ("libutilityFunctionObjects.so");
        writeControl    writeTime;
        executeControl  writeTime;

        code
        #{
        #};

        codeExecute
        #{
            const volVectorField& U
            (
                mesh().lookupObject<volVectorField>("U")
            ); 
            volScalarField strainRate(sqrt(2.0)*mag(symm(fvc::grad(U))));
            strainRate.write();
        #};
    }
    
    setDeltaT
    {
        type coded;
        libs            ("libutilityFunctionObjects.so");
        codeExecute
        #{
            const Time& runTime = mesh().time();
            if (runTime.userTimeValue() >= -15.0)
            {
                const_cast<Time&>(runTime).setDeltaT
                (
                    runTime.userTimeToTime(0.025)
                );
            }
        #};
    }
    
    #includeFunc MachNo
    
    #includeFunc yPlus
    
    #includeFunc Q
    
    #includeFunc vorticity
    
    #includeFunc Lambda2
    
    particles
    {
        libs        ("liblagrangianFunctionObjects.so");
        type        particles;
    }
}
}

输出重组速度场、传输的标量、某个场的最大值、阻力系数、流线图、残差、监控流率

functions
{
    #includeFunc reconstruct(phi)\\重组速度场
     
    #includeFunc scalarTransport(T, alphal=1, alphat=1)
     
    minMaxp
    {
        type        fieldMinMax;
        functionObjectLibs ("libfieldFunctionObjects.so");
        fields
        (
             U
        );
        location        no;
        writeControl    timeStep;
        writeInterval   1;
    }
    
    forces
    {
        type            forceCoeffs;
        libs            ("libforces.so");
        writeControl    timeStep;
        writeInterval   1;
    
        patches         ("motorBike.*");
        rho             rhoInf;      // Indicates incompressible
        log             true;
        rhoInf          1;           // Redundant for incompressible
        liftDir         (0 0 1);
        dragDir         (1 0 0);
        CofR            (0.72 0 0);  // Axle midpoint on ground
        pitchAxis       (0 1 0);
        magUInf         20;
        lRef            1.42;        // Wheelbase length
        Aref            0.75;        // Estimated
    }
    
    streamlines1
    {
        type            streamlines;
    
        libs            ("libfieldFunctionObjects.so");
    
        // Output every
        writeControl    writeTime;
    
        // Write format
        setFormat       vtk;
    
        // Track forward (+U) or backward (-U) or both
        direction       forward;
    
        // Names of fields to sample. Should contain above velocity field!
        fields (p U);
    
        // Steps particles can travel before being removed
        lifeTime        10000;
    
        // Number of steps per cell (estimate). Set to 1 to disable subcycling.
        nSubCycle 5;
    
        // Cloud name to use
        cloudName       particleTracks;
    
        // Seeding method.
        seedSampleSet
        {
            type        lineUniform;
            start       (-1.001 1e-7 0.0011);
            end         (-1.001 1e-7 1.0011);
            nPoints     20;
        }
    }
    
    #includeFunc streamlinesLine(funcName=streamlines, start=(0 0.5 0), end=(9 0.5 0), nPoints=24, U)
    
    #includeFunc residuals
    
    #includeFunc patchFlowRate(patch=outlet1)
    #includeFunc faceZoneFlowRate(name=fz1)
    #includeFunc patchFlowRate(patch=outlet2)
    #includeFunc faceZoneFlowRate(name=fz2)
}

输出边界积分值、场的体积分、场的平均值、转换坐标系、屏蔽迭代求解器输出

functions
{
    #includeFunc patchAverage(patch=outlet, fields=(p U))
    #includeFunc patchAverage(patch=inlet, fields=(p U)) 
    
    volFieldValue1
    {
        type            volFieldValue;
        libs            ("libfieldFunctionObjects.so");
        writeControl    writeTime;
        log             yes;
        writeFields     no;
        regionType      all;
        name            outlet;
        operation       volIntegrate;//min, max, etc.
        //weightField     phi;
        fields
        (
            alpha
        );
    }
    
    volFieldValue1
    {
        type            volFieldValue;
        libs            ("libfieldFunctionObjects.so");
        writeControl    writeTime;
        log             yes;
        writeFields     no;
        regionType      all;
        name            outlet;
        operation       volAverage;//min, max, etc.
        fields
        (
            T
        );
    }
    
    cartesianToCylindrical
    {
        type        cylindrical;
        libs        ("libfieldFunctionObjects.so");

        origin      (0 0 0);
        axis        (1 0 0);

        field       U;

        writeControl    outputTime;
        writeInterval   1;
    }
}

屏蔽迭代求解器输出

DebugSwitches
{
    SolverPerformance   0;
}

fvModels自定义源项

//fvModels如下指定
USource
{
    type            coded;
    selectionMode   all;
    field           U;

    codeAddSup
    #{
        const volVectorField& U = eqn.psi();
        vectorField& USource = eqn.source();

        const fvMesh& mesh = U.mesh();
        const DimensionedField<scalar, volMesh>& V = mesh.V();
        const volVectorField& C = mesh.C();
        scalar A = 0.919;
        scalar Cd = 0.15;

        forAll(U, i)
        {
            vector position = C[i];
            scalar pz = position.z();
            if (pz < 10.0)
            {
                USource[i] -= -Cd*A*mag(U[i])*U[i]*V[i];
            }
        }
    #};
}