OpenFOAM代码

非均一分布的速度进口

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
    #{
        //  调用mesh
        const fvMesh& mesh = this->patch().boundaryMesh().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");
        
        //  调用fvc
        volTensorField gradU = fvc::grad(U);
        
        // 计算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");
    #};
}

OpenFOAM中基本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())));                   // 形变率的双点积

fvMesh mesh1                                                // 创建网格1
(
    IOobject
    (
        "mesh1",
        runTime.timeName(),
        runTime,
        IOobject::MUST_READ
    )
);

fvMesh mesh2                                                // 创建网格2
(
    IOobject
    (
        "mesh2",
        runTime.timeName(),
        runTime,
        IOobject::MUST_READ
    )
);

IOdictionary physicalProperties                             // 声明字典文件
(
    IOobject
    (
        "physicalProperties",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
    )
);

bool training = 
    readBool(physicalProperties.lookup("training"));        // 去字典文件中读取关键词
    
scalar RR(readScalar(generalProperties.lookup("RR")));      // 去字典文件中读取关键词

const volScalarField& T                                     // 调取引用场
(
    mesh().lookupObject<volScalarField>("T")
);
    
volScalarField T                                            // 体心场
(
    IOobject
    (
        "T",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

volScalarField B                                            
(
    IOobject
    (
        "B",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar(dimless, 0.0),
    A.boundaryField().types()                                // 把A的边界条件赋值给B
);

volScalarField B
(
    IOobject
    (
        "B",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    mesh,
    dimensionedScalar(dimless, 0.0),
    zeroGradientFvPatchScalarField::typeName                // B的边界条件均为零梯度
);

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
);

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

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
        )                          
    );                            
} 

volScalarField y(wallDist::New(mesh).y());                 // 壁面距离，需要#include "wallDist.H"

U.primitiveFieldRef() = ...                                // 修改U的内部场的值
U.ref() = ...                                              // 修改U的内部场的值
U.boundaryFieldRef() = ...                                 // 修改U的边界场的值
bound(S, dimensionedScalar(S.dimensions(), 0));            // 对S场做有界处理

volScalarField& alphat =                                   // 强制修改const常量
    const_cast<volScalarField&>
    (
        mesh().lookupObject<volScalarField>("alphat")
    );  
    
List<label> markedCell;                                    // 声明一个数组并赋值
for (int i = 0; i < 10; i++)
{
    markedCell.append(i);
}

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

forAll(mesh.boundary(), patchi)                           // 判断壁面的边界类型
{                                                         // 需要#include "wallFvPatch.H"
    if (isA<wallFvPatch>(mesh.boundary()[patchi]))
    {

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

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

    }
}

禁用log输出

SolverPerformance<scalar>::debug = 0;

更改变量近壁网格体心值

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

矩阵操作

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;

显性、隐性离散

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)           \\显性梯度项离散

在算例层面进行后处理操作：

// 输出多相界面、界面高度要对libwaves.so进行包含
libs
(
    "libwaves.so"
);

// functions需要写入到算例的controlDict里
functions
{
    // 输出等值面
    #includeFunc isoSurface(isoField=alpha.water, isoValue=0.5, fields=())  
    
    // 输出波动速度场，注意，每个时间步文件夹下的UMean的数值要正确
    uPrime                                                
    {
        type            subtract;
        libs            ("libfieldFunctionObjects.so");
        fields          (U UMean);
        result          uPrime;
        executeControl  writeTime;
        writeControl    writeTime;
    }
    
    // 输出波动速度场自身的张量积
    uPrime2
    {
        type            multiply;
        libs            ("libfieldFunctionObjects.so");
        fields          (uPrime uPrime);
        result          uPrime2;
        executeControl  writeTime;
        writeControl    writeTime;
    }
    
    // 输出空气龄
    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;
    }
    
    // 输出热通量
    wallHeatFlux1
    {
        type        wallHeatFlux;
        libs        ("libfieldFunctionObjects.so");
        region      fluid;
        patches     (".*Wall");
    }
    
    // 输出形变率
    dymmy
    {
        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();
        #};
    }
    
    // 更改deltaT
    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)
                );
            }
        #};
    }
    
    // 输出缓存场
    cacheTemporaryObjects(kEpsilon:G);
    #includeFunc writeObjects(kEpsilon:G)
    
    // 输出失踪粒子
    particles
    {
        libs        ("liblagrangianFunctionObjects.so");
        type        particles;
    }
    
    // 输出场的最大值最小值
    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
    }
    
    // 输出流线
    #includeFunc streamlinesLine
    (    
        funcName=streamlines, start=(0 0.5 0), end=(9 0.5 0), nPoints=24, 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;
    }
    
    #includeFunc patchFlowRate(patch=outlet1)              // 计算边界流率
    
    #includeFunc faceZoneFlowRate(name=fz1)                // 计算某个faceZone的流率
    
    #includeFunc fieldAverage(U, p, prime2Mean = yes)      // 输出场的时间平均
    
    #includeFunc residuals                                 // 输出残差
    
    #includeFunc mag(UPrime2Mean)                          // 输出场的模
    
    #includeFunc components(U)                             // 输出速度分量
    
    #includeFunc MachNo                                    // 输出马赫数
  
    #includeFunc yPlus                                     // 输出y+
    
    #includeFunc Q                                         // 输出Q
    
    #includeFunc vorticity                                 // 输出涡量
    
    #includeFunc Lambda2                                   // 输出lambda
    
    #includeFunc reconstruct(phi)                          // 输出重组的速度场
    
    #includeFunc scalarTransport(T, alphal=1, alphat=1)    // 输出s标量传输
    
    #includeFunc patchAverage(patch=outlet, fields=(p U))  // 输出边界积分值
}

屏蔽迭代求解器输出

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];
            }
        }
    #};
}

limitT
{
    type            limitTemperature;
    active          yes;     
    selectionMode   all;   
    min             200;            
    max             20000;   
}

limitp
{
    type       limitPressure;

    minFactor  0.4;
    maxFactor  1.5;
}

拉格朗日算例输出相分数场：

cloudFunctions //需要放在cloudProperties
{ 
    alpha
    { 
        type voidFraction; 
    } 
}