CFD中文网

你试试这个安装，我用的openfoam v2312，安装没有问题。

sudo apt-get install mercurial
hg clone http://hg.code.sf.net/p/openfoam-extend/swak4Foam
cd swak4Foam
hg update develop
// 编译swak4foam，需要较长时间
./AllwmakeAll

下面是我需要计算的方程：

fvScalarMatrix TbEqn
(
    - fvm::div(-phi, Tb) 
    + fvm::laplacian(DT, Tb)
);
TbEqn -= Q;
TbEqn.solve();

其中，DT为热传导系数，Q为数值1的标量场，phi通过SIMPLE求解。
设置出口边界条件为：

outlet
{
    type                 groovyBC;
    variables            "Tk=DT;h=U&normal();Tinf=0;f=1/(1+Tk/(h*mag(delta())));";
    valueExpression      "Tinf";
    gradientExpression   "0";
    fractionExpression   "f";
    value                uniform 0;
}

设置的divSchemes为：

div(-phi,Tb)     Gauss linearUpwind grad(Tb);

使用的solvers为：

"(Tb)"
{
     solver                  PBiCGStab;
     preconditioner          DILU;
     tolerance               1e-6;
     relTol                  0;
 }

上述方程是用于拓扑优化求解的，但是优化过程中，经常出现浮点数错误“Signal: Floating point exception (8)”，如何有效解决这个问题呢？
尝试过在h*mag(delta()中添加一个小数，但是也同样出现浮点数报错。

大佬们，想请教下为什么OpenFOAM自定义并行计算时，初始化（这边初始化是指下述的代码）有时候很快2-3s，有时候需要很久1-2min（测试的时候没有其他程序运行）。
下述是我自定义的并行计算初始化代码：

scalarField MyField(5*N,0);
List<scalarField> AllMyField_List(Pstream::nProcs());
AllMyField_List[Pstream::myProcNo()] = MyField;
Pstream::gatherList(AllMyField_List);
Pstream::scatterList(AllMyField_List);
scalarField AllMyField(ListListOps::combine<scalarField>(AllMyField_List, accessOp<scalarField>()));

其中，N为所有网格单元的数量，N=86400。程序是有个变量，我将其分配到各个线程进行计算，在总线程中进行收集。
此外，程序中不止这一个变量，有五六个需要同样操作的变量。之前测试过，如果N很小的时候，初始化过程很快，但是当N增大时，初始化时间就完全不一样了。想请教下如何解决这个bug。

可以在求解器层面实现你想要的功能。
@zhangK 在 paraview提取部分结构 中说：

相分布稳定之后导出第一相的结构

在计算稳定后直接输出第一项结果，比如

firstphasefield.write();

@zhangK 在 paraview提取部分结构 中说：

进一步计算scalar transport

将上述的firstphasefield文件放置在0文件夹下，在第二个程序计算时读取该文件用于新的方程计算：

volScalarField firstphasefield
(
    IOobject
    (
        "firstphasefield",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    mesh
);
firstphasefield.read();

@李东岳 李老师 动图中iter是指迭代步数。目前设置的是计算速度场计算200步完成后，再继续伴随速度场，这边为了减少gif的大小，我只生成了前50步的结果。
在10-15步的时候，此时计算的Ua是正常的，但是随着迭代步增加，计算就不收敛了。一直没搞明白为什么。

这是其他密度场下计算得到的Ua场：

现在的问题是，有些密度场下可以计算，有些就异常了。（密度场是通过0时刻读入的，作为一个初始值）

volScalarField xh
(
    IOobject
    (
        "xh",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    mesh,
    scalar(voluse),
    zeroGradientFvPatchScalarField::typeName
);
xh.read();

我尝试用其他的密度场进行测试，也发现同样的问题。
密度场如下：

速度场如下：

但是伴随速度场Ua计算还是异常，下图是Ua迭代计算前50步的计算情况，似乎计算后越来越不收敛了：

想请教下各位老师，有没有解决方法？

下图是我基于不同的密度场(xh)进行速度场(U)和伴随速度场(Ua)的计算结果，xh有些小差异，但是伴随速度场计算出来完全不一样，右图计算出错了。

速度场计算是基于SA模型进行，在UEqn方程中增加了源项：

tmp<fvVectorMatrix> tUEqn
(
    fvm::div(phi, U)
  + turbulence->divDevReff(U)
  + fvm::Sp(alpha, U)
 ==
    fvOptions(U)
);
fvVectorMatrix& UEqn = tUEqn.ref();
UEqn.relax();

其中alpha是由密度场xh插值进行的，alphaMax=2.5e5，qu=0.01。

alpha=alphaMax*qu*(1-xh)/(qu+xh);

对于伴随速度场的计算，和UEqn类似，如下所示：

volVectorField adjointTransposeConvectiona((fvc::grad(U) & Ua));
tmp<fvVectorMatrix> tUaEqn
  (
      fvm::div(-phi, Ua)
    + adjointTransposeConvectiona
    + turbulence->divDevReff(Ua)
    + fvm::Sp(alpha, Ua)
   ==
      fvOptions(Ua)
  );
  fvVectorMatrix& UaEqn = tUaEqn.ref();
  UaEqn.relax();

请大家帮忙看看，会有哪些可能性导致代码计算出现异常。

@1064168551 在 openfoam自定义边界条件编译时报错，恳请论坛内各位老师解惑 中说：

const incompressible::RASModel& rasModel =
db().lookupObjectFoam::incompressible::RASModel("RASProperties");

scalarField nueff = rasModel.nuEff()().boundaryField()[patch().index()];

你这边伴随边界条件中并没有用到nueff的值，可以直接把这几行代码注释掉。

@李东岳 在 不使用湍流模型，而是直接求解器中实现湍流计算 中说：

有一些代码不知道有啥作用，那就留着，留着总比删了更保险。

好像是的，对结果不影响的都会跟着作者留下来了。

@李东岳 在 不使用湍流模型，而是直接求解器中实现湍流计算 中说：

这个在extend那面很常见。

extend我只是听过，但没有用过。现在我常使用的版本一个是Org 6.0的，另外一个是ESI v2312的。

@李东岳 在 不使用湍流模型，而是直接求解器中实现湍流计算 中说：

你这个是什么版本，storePrevIter()好久没遇到过了，还是说是extend版本？

李老师，使用的是OpenFOAM6.0版本。
我是在github的代码基础上进行移植的，它里面用到了storePrevIter()，但我还不清楚storePrevIter()有没有作用，删除了也可以正常计算。

@李东岳 是的是的 目前打算通过在求解器中植入的方式实现

@李东岳 在 不使用湍流模型，而是直接求解器中实现湍流计算 中说：

如果pitzDaily没问题。再debug你的。

李老师，我发现问题了，在调用SA模型计算的时候，我忘记把turbulence->correct()；加进去了，并且算例那边使用的是层流进行计算，现在改成SA模型，计算结果是一样的了。
pitzDaily试过，没有问题。

@李东岳 在 不使用湍流模型，而是直接求解器中实现湍流计算 中说：

我不太熟悉伴随湍流模型是什么 :-)

这边的伴随湍流模型只是一个叫法，可以这样理解，原始的湍流模型为湍流模型1，伴随湍流模型为湍流模型2（自定义的）。
在求解器中会计算控制方程（调用的是湍流模型1），计算完控制方程后会计算伴随方程（调用的是湍流模型2）。如下图所示：

由于constant/turbulenceProperties中只能调用一种湍流模型，我想在一个算例中调用两种不同湍流模型，要怎么实现呢？

@李东岳 在 不使用湍流模型，而是直接求解器中实现湍流计算 中说：

我测试的是pitzDaily。你测试的是什么算例。

李老师，我这边测试的是自己的算例，不是官方自带的。
下图是我的模型，由于模型的对称性，算例miny边界设置的是对称边界条件，冷却液从左侧流入，从右侧流出。这边设置目的是获取直管道的粘性耗散值作为拓扑优化的约束函数值，因此模型也不是最原始的直管模型，而是通过gamma场来区分固体域(gamma=0)和流体域(gamma=1)，如下图所示：

在求解器层面，通过在UEqn添加源项来实现固体域和流体域的区分：

fvm::Sp(alpha, U)

对应的nuTildaEqn源项为：

fvm::Sp(alpha, nuTilda)

其中，alpha为自定义的参数，用来惩罚速度：

volScalarField alpha(alphaMax*qu*(1-gamma)/(qu+gamma));

上述方法在层流模型中已经实现拓扑优化了，现在想将其应用到湍流中实现湍流传热的拓扑优化。因此我打算在求解器层面实现湍流求解，以便于后续将推导的伴随方程通过同样的方式移植进去。
但目前在求解器层面计算的结果和直接调用湍流模型的计算结果不一样。
求解器层面计算结果如下图所示：

调用湍流模型计算结果如下图所示：

@李东岳 在 不使用湍流模型，而是直接求解器中实现湍流计算 中说：

storePrevIter()这种操作好久都没看过了。如果是老版本，你要不要试试把对流项跟扩散项两行兑换一下。

我尝试了上述的修改，换成新版本代码或者使用旧代码对换对流项和扩散项，计算出来的结果都不正确，不知道是不是有啥bug。

另外想请教新的问题，如果是调用湍流模型的话，我想要写一个伴随湍流模型，但是要怎么在一个算例中即能实现湍流模型调用，又能实现伴随湍流模型计算呢？
即湍流模型调用是通过constant/turbulenceProperties进行选择和调用，我想写一个constant/adjointturbulenceProperties来调用伴随湍流模型，如何实现？

@xpqiu 在 不使用湍流模型，而是直接求解器中实现湍流计算 中说：

先求解 UEqn，pEqn，然后再用修正的速度来构建和求解nuTildaEqn

我尝试把求解nuTildaEqn的代码移到后面，计算出来的结果也不正确。
想请教下，您这边提到的先求解UEqn和pEqn，是它们完全收敛后再求解nuTildaEqn嘛？

我正在修改SA湍流的模型计算方式（通过修改turbulence->divDevReff(U)项，直接在求解器中表示，以便后续我进行其他算法开发），但计算结果和直接使用SA模型计算结果不一样，请教下大家如何解决？

{       
    const volScalarField chi = SpalartAllmaras::chi(nuTilda, nu);
    const volScalarField fv1 = SpalartAllmaras::fv1(chi, cv1);
    const volScalarField Stilda = SpalartAllmaras::Stilda(chi, fv1, nuTilda, U, d, kappa, cs);

    tmp<fvScalarMatrix> nuTildaEqn
    (
        fvm::div(phi, nuTilda)
      - fvm::laplacian(SpalartAllmaras::nuTildaEff(nuTilda, sigma, nu), nuTilda)
      - cb2/sigma*magSqr(fvc::grad(nuTilda))
     ==
        cb1*Stilda*nuTilda
      - fvm::Sp(cw1*SpalartAllmaras::fw(Stilda, nuTilda, d, kappa, cw2, cw3)*nuTilda/sqr(d), nuTilda)
      + fvOptions(nuTilda)
    );

    nuTildaEqn.ref().relax(0.7); 
    fvOptions.constrain(nuTildaEqn.ref());  
    solve(nuTildaEqn);
    fvOptions.correct(nuTilda);
    bound(nuTilda, dimensionedScalar("0", nuTilda.dimensions(), 0.0));
    nuTilda.correctBoundaryConditions();
    nut = nuTilda*fv1;
    nut.correctBoundaryConditions();

    MRF.correctBoundaryVelocity(U);
    tmp<fvVectorMatrix> tUEqn
    (
        fvm::div(phi, U)
        - fvm::laplacian((nu+nut), U)
	- fvc::div((nu+nut)*dev2(T(fvc::grad(U))))
        ==
        fvOptions(U)
    );
    fvVectorMatrix& UEqn = tUEqn.ref();

    UEqn.relax();
    fvOptions.constrain(UEqn);
    if(simple.momentumPredictor())
    {
        solve(UEqn == -fvc::grad(p));
        fvOptions.correct(U);
    }

    //****************************************//
    volScalarField rAU(1.0/UEqn.A());
    volVectorField HbyA(constrainHbyA(rAU*UEqn.H(), U, p));
    surfaceScalarField phiHbyA("phiHbyA", fvc::flux(HbyA));
    MRF.makeRelative(phiHbyA);
    adjustPhi(phiHbyA, U, p);

    tmp<volScalarField> rAtU(rAU);

    if (simple.consistent())
    {
        rAtU = 1.0/(1.0/rAU - UEqn.H1());
        phiHbyA +=
            fvc::interpolate(rAtU() - rAU)*fvc::snGrad(p)*mesh.magSf();
        HbyA -= (rAU - rAtU())*fvc::grad(p);
    }

    tUEqn.clear();

    // Update the pressure BCs to ensure flux consistency
    constrainPressure(p, U, phiHbyA, rAtU(), MRF);

    // Non-orthogonal pressure corrector loop
    while (simple.correctNonOrthogonal())
    {
        fvScalarMatrix pEqn
        (
            fvm::laplacian(rAtU(), p) == fvc::div(phiHbyA)
        );

        pEqn.setReference(pRefCell, pRefValue);
        pEqn.solve();

        if(simple.finalNonOrthogonalIter())
        {
            phi = phiHbyA - pEqn.flux();
        }
    }

    // Explicitly relax pressure for momentum corrector
    p.relax();

    // Momentum corrector
    U = HbyA - rAtU()*fvc::grad(p);
    U.correctBoundaryConditions();
    fvOptions.correct(U);

    U.storePrevIter();
    p.storePrevIter();
    phi.storePrevIter();
}

有研究湍流传热拓扑优化的嘛，我看OpenFOAM v2312版本中adjointoptimisation中增加了拓扑优化功能，但我尝试修改库文件，并加入温度场和伴随温度场，并没有得到想要的效果，想请教大家有没有这块的学习资料，以便我能正确算出优化构型。

@Do1975 在 paraview图例设置 中说：

@吴建民 您这边的问题是不是可以理解为，把colorbar按最大值和最小值分割成等间距的刻度。之前我通过pvpython来出图，可以实现这种效果，但是是连续型的，离散型没有试过。

这是效果图：

对应的处理代码：

tRange = appfoam.CellData[field].GetRange()
tLUT.RescaleTransferFunction(tRange[0], tRange[1])
tPWF.RescaleTransferFunction(tRange[0], tRange[1])
            
# Set custom labels for the color bar
numberOfLabels = 5
labelPositions = [tRange[0] + j * (tRange[1] - tRange[0]) / (numberOfLabels - 1) for j in range(numberOfLabels)]
tLUTColorBar.CustomLabels = labelPositions

@吴建民 您这边的问题是不是可以理解为，把colorbar按最大值和最小值分割成等间距的刻度。之前我通过pvpython来出图，可以实现这种效果，但是是连续型的，离散型没有试过。

@李东岳 好勒好勒，谢谢李老师回答。感觉又接触到很多没学过的知识了

@李东岳 在 OpenFOAM散热器自然散热问题 中说：

1. 在一个时间步里面压力循环执行一次从理论上来说也是不正常的。
2. 瞬态求解器不能添加松弛因子。

李老师，有个不解的地方，为什么执行一次不正常，像simplefoam之类的求解器，每个时间步内压力循环执行不也是只有一次吗？
另外程序是基于buoyantBoussinesqSimpleFoam求解器进行修改的，应该属于稳态求解。