不使用湍流模型，而是直接求解器中实现湍流计算
-
我正在修改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(); } -
@李东岳 在 不使用湍流模型，而是直接求解器中实现湍流计算 中说：
我测试的是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来调用伴随湍流模型，如何实现？ -
@李东岳 在 不使用湍流模型，而是直接求解器中实现湍流计算 中说：
如果pitzDaily没问题。再debug你的。
李老师，我发现问题了，在调用SA模型计算的时候，我忘记把turbulence->correct()；加进去了，并且算例那边使用的是层流进行计算，现在改成SA模型，计算结果是一样的了。
pitzDaily试过，没有问题。@李东岳 在 不使用湍流模型，而是直接求解器中实现湍流计算 中说：
我不太熟悉伴随湍流模型是什么 :-)
这边的伴随湍流模型只是一个叫法，可以这样理解，原始的湍流模型为湍流模型1，伴随湍流模型为湍流模型2（自定义的）。
在求解器中会计算控制方程（调用的是湍流模型1），计算完控制方程后会计算伴随方程（调用的是湍流模型2）。如下图所示：
由于constant/turbulenceProperties中只能调用一种湍流模型，我想在一个算例中调用两种不同湍流模型，要怎么实现呢？ -
@李东岳 在 不使用湍流模型，而是直接求解器中实现湍流计算 中说：
你这个是什么版本，storePrevIter()好久没遇到过了，还是说是extend版本？
李老师，使用的是OpenFOAM6.0版本。
我是在github的代码基础上进行移植的，它里面用到了storePrevIter()，但我还不清楚storePrevIter()有没有作用，删除了也可以正常计算。 -
@李东岳 在 不使用湍流模型，而是直接求解器中实现湍流计算 中说：
有一些代码不知道有啥作用，那就留着，留着总比删了更保险。
好像是的，对结果不影响的都会跟着作者留下来了。
@李东岳 在 不使用湍流模型，而是直接求解器中实现湍流计算 中说：
这个在extend那面很常见。
extend我只是听过，但没有用过。现在我常使用的版本一个是Org 6.0的，另外一个是ESI v2312的。
