CFD中文网

@李东岳 非常感谢李老师不厌其烦的指导！！tmp的问题已经解决了，不过又出现了量纲的问题，我自己再看看。

@李东岳 感谢李老师，前面整个correct函数是这样的

template<class BasicTurbulenceModel>
void kOmegaSSTIDDESTr<BasicTurbulenceModel>::correct()
{
    if (!this->turbulence_)
    {
        return;
    }

    // Local references
    const alphaField& alpha = this->alpha_;
    const rhoField& rho = this->rho_;
    const surfaceScalarField& alphaRhoPhi = this->alphaRhoPhi_;
    const volVectorField& U = this->U_;
    volScalarField& nut = this->nut_;
    fv::options& fvOptions(fv::options::New(this->mesh_));

    BasicTurbulenceModel::correct();

    volScalarField::Internal divU
    (
        fvc::div(fvc::absolute(this->phi(), U))()()
    );

    tmp<volTensorField> tgradU = fvc::grad(U);
    volScalarField S2(2*magSqr(symm(tgradU())));
    volScalarField::Internal GbyNu(dev(twoSymm(tgradU()())) && tgradU()());
    volScalarField::Internal G(this->GName(), nut()*GbyNu);
    tgradU.clear();

    // Update omega and G at the wall
    omega_.boundaryFieldRef().updateCoeffs();

    volScalarField CDkOmega
    (
        (2*alphaOmega2_)*(fvc::grad(k_) & fvc::grad(omega_))/omega_
    );

    volScalarField F1(this->F1(CDkOmega));
    volScalarField F23(this->F23());

    {
        volScalarField::Internal gamma(this->gamma(F1));
        volScalarField::Internal beta(this->beta(F1));

        // Turbulent frequency equation
        tmp<fvScalarMatrix> omegaEqn
        (
            fvm::ddt(alpha, rho, omega_)
          + fvm::div(alphaRhoPhi, omega_)
          - fvm::laplacian(alpha*rho*DomegaEff(F1), omega_)
         ==
            alpha()*rho()*gamma
           *min
            (
                GbyNu,
                (c1_/a1_)*betaStar_*omega_()
               *max(a1_*omega_(), b1_*F23()*sqrt(S2()))
            )
          - fvm::SuSp((2.0/3.0)*alpha()*rho()*gamma*divU, omega_)
          - fvm::Sp(alpha()*rho()*beta*omega_(), omega_)
          - fvm::SuSp
            (
                alpha()*rho()*(F1() - scalar(1))*CDkOmega()/omega_(),
                omega_
            )
          + Qsas(S2(), gamma, beta)
          + omegaSource()
          + fvOptions(alpha, rho, omega_)
        );

        omegaEqn.ref().relax();
        fvOptions.constrain(omegaEqn.ref());
        omegaEqn.ref().boundaryManipulate(omega_.boundaryFieldRef());
        solve(omegaEqn);
        fvOptions.correct(omega_);
        bound(omega_, this->omegaMin_);
    }
k方程。。。

@李东岳 李老师，前面有个epsilonByk函数，

template<class BasicTurbulenceModel>
tmp<volScalarField::Internal> kOmegaSSTDES<BasicTurbulenceModel>::epsilonByk
(
    const volScalarField& F1,
    const volTensorField& gradU
) const
{
    volScalarField CDES(this->CDES(F1));
    return sqrt(this->k_())/dTilda(mag(gradU), CDES)()();
}

但是我把tgradU换成gradU,编译报错：

kOmegaSSTIDDESTr.C: In member function ‘virtual void Foam::LESModels::kOmegaSSTIDDESTr<BasicTurbulenceModel>::correct()’:
kOmegaSSTIDDESTr.C:1248:59: error: ‘gradU’ was not declared in this scope
       - fvm::Sp(alpha()*rho()*betaStar_*omega_*dTilda(mag(gradU), CDES)(), k_)
                                                           ^~~~~
kOmegaSSTIDDESTr.C:1248:59: note: suggested alternative: ‘tgradU’
       - fvm::Sp(alpha()*rho()*betaStar_*omega_*dTilda(mag(gradU), CDES)(), k_)
                                                           ^~~~~
                                                           tgradU
/home/fxr/OpenFOAM/OpenFOAM-7/wmake/rules/General/transform:25: recipe for target 'Make/linux64GccDPInt64Opt/IDDESturbulentTransportModels.o' failed
make: *** [Make/linux64GccDPInt64Opt/IDDESturbulentTransportModels.o] Error 1

@李东岳 谢谢，按照您的方法试了，就是这行的问题，一去掉就可以算

- fvm::Sp(alpha()*rho()*betaStar_*omega_*dTilda(mag(tgradU), CDES)(), k_)

@李东岳 李老师我加了（）还是不行啊。。。纠结两天了。。。。
顺手贴下问题：主要是在of org版本上添加kOmegaSSTIDDES模型，但是org版本只有kOmegaSSTDES,所以我把openfoam +版本的kOmegaSSTIDDES移植过来，涉及到长度尺度dTilda()函数，目前跟露珠遇到的问题一样，编译通过但计算一步就崩，报错信息：

--> FOAM FATAL ERROR: 
tmp<N4Foam14GeometricFieldINS_6TensorIdEENS_12fvPatchFieldENS_7volMeshEEE> deallocated

    From function const T& Foam::tmp<T>::operator()() const [with T = Foam::GeometricField<Foam::Tensor<double>, Foam::fvPatchField, Foam::volMesh>]
    in file /home/fxr/OpenFOAM/OpenFOAM-7/src/OpenFOAM/lnInclude/tmpI.H at line 278.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::error::abort() at ??:?
#2  Foam::tmp<Foam::GeometricField<Foam::Tensor<double>, Foam::fvPatchField, Foam::volMesh> >::operator()() const at ??:?
#3  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::mag<Foam::Tensor<double>, Foam::fvPatchField, Foam::volMesh>(Foam::tmp<Foam::GeometricField<Foam::Tensor<double>, Foam::fvPatchField, Foam::volMesh> > const&) at ??:?
#4  Foam::LESModels::kOmegaSSTIDDESTr<Foam::IncompressibleTurbulenceModel<Foam::transportModel> >::correct() at ??:?
#5  ? in "/home/fxr/OpenFOAM/OpenFOAM-7/platforms/linux64GccDPInt64Opt/bin/pimpleFoam"
#6  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#7  ? in "/home/fxr/OpenFOAM/OpenFOAM-7/platforms/linux64GccDPInt64Opt/bin/pimpleFoam"
Aborted (core dumped)

dTilda代码：

template<class BasicTurbulenceModel>
tmp<volScalarField> kOmegaSSTIDDESTr<BasicTurbulenceModel>::dTilda
(
    const volScalarField& magGradU,
    const volScalarField& CDES
) const
{
    const volScalarField& k = this->k_;
    const volScalarField& omega = this->omega_;

    const volScalarField lRAS(sqrt(k)/(this->betaStar_*omega));
    const volScalarField lLES(CDES*this->delta());

    const volScalarField alpha(this->alpha());
    const volScalarField expTerm(exp(sqr(alpha)));

    tmp<volScalarField> fB = min(2*pow(expTerm, -9.0), scalar(1));
    tmp<volScalarField> fe1 =
        2*(pos0(alpha)*pow(expTerm, -11.09) + neg(alpha)*pow(expTerm, -9.0));
    tmp<volScalarField> fe2 = 1 - max(ft(magGradU), fl(magGradU));
    tmp<volScalarField> fe = max(fe1 - 1, scalar(0))*fe2;

    const volScalarField fdTilda(max(1 - fdt(magGradU), fB));

   
    return max
    (
        fdTilda*(1 + fe)*lRAS + (1 - fdTilda)*lLES,
        dimensionedScalar("SMALL", dimLength, SMALL)
    );
}

后面k方程中包含了dTilda的调用：

 // Turbulent kinetic energy equation
    volScalarField CDES(this->CDES(F1));
    tmp<fvScalarMatrix> kEqn
    (
        fvm::ddt(alpha, rho, k_)
      + fvm::div(alphaRhoPhi, k_)
      - fvm::laplacian(alpha*rho*DkEff(F1), k_)
     ==
        alpha()*rho()*Pk(G)
      - fvm::SuSp((2.0/3.0)*alpha()*rho()*divU, k_)
      - fvm::Sp(alpha()*rho()*betaStar_*omega_*dTilda(mag(tgradU), CDES)(), k_)
      + kSource()
      + fvOptions(alpha, rho, k_)
    );

各位老师有没有什么建议啊？？？谢谢了

@悬浮的猪 请问楼主解决了么？我修改湍流模型也遇到了一样的问题

@李东岳 谢谢李老师，计算之前controlDict已经事先挂载过了

libs ("test1_libIncompressibleTurbulenceModels.so");

请教各位老师，在DES模型的基础上做了修改，wmake后：

wmakeLnInclude: linking include files to ./lnInclude
Making dependency list for source file myTurbulentTransportModels.C
g++ -std=c++11 -m64 -DOPENFOAM=1806 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=64 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas -O3  -DNoRepository -ftemplate-depth-100 -I../turbulenceModels/lnInclude -I/home/fxr/OpenFOAM/OpenFOAM-v1806/src/transportModels -I/home/fxr/OpenFOAM/OpenFOAM-v1806/src/finiteVolume/lnInclude -I/home/fxr/OpenFOAM/OpenFOAM-v1806/src/meshTools/lnInclude -I/home/fxr/OpenFOAM/OpenFOAM-v1806/src/TurbulenceModels/turbulenceModels/lnInclude -I/home/fxr/OpenFOAM/OpenFOAM-v1806/src/TurbulenceModels/incompressible/lnInclude -IlnInclude -I. -I/home/fxr/OpenFOAM/OpenFOAM-v1806/src/OpenFOAM/lnInclude -I/home/fxr/OpenFOAM/OpenFOAM-v1806/src/OSspecific/POSIX/lnInclude   -fPIC -c turbulentTransportModels/myTurbulentTransportModels.C -o Make/linux64GccDPInt64Opt/turbulentTransportModels/myTurbulentTransportModels.o
g++ -std=c++11 -m64 -DOPENFOAM=1806 -Dlinux64 -DWM_ARCH_OPTION=64 -DWM_DP -DWM_LABEL_SIZE=64 -Wall -Wextra -Wold-style-cast -Wnon-virtual-dtor -Wno-unused-parameter -Wno-invalid-offsetof -Wno-attributes -Wno-unknown-pragmas -O3  -DNoRepository -ftemplate-depth-100 -I../turbulenceModels/lnInclude -I/home/fxr/OpenFOAM/OpenFOAM-v1806/src/transportModels -I/home/fxr/OpenFOAM/OpenFOAM-v1806/src/finiteVolume/lnInclude -I/home/fxr/OpenFOAM/OpenFOAM-v1806/src/meshTools/lnInclude -I/home/fxr/OpenFOAM/OpenFOAM-v1806/src/TurbulenceModels/turbulenceModels/lnInclude -I/home/fxr/OpenFOAM/OpenFOAM-v1806/src/TurbulenceModels/incompressible/lnInclude -IlnInclude -I. -I/home/fxr/OpenFOAM/OpenFOAM-v1806/src/OpenFOAM/lnInclude -I/home/fxr/OpenFOAM/OpenFOAM-v1806/src/OSspecific/POSIX/lnInclude   -fPIC -shared -Xlinker --add-needed -Xlinker --no-as-needed Make/linux64GccDPInt64Opt/turbulentTransportModels/myTurbulentTransportModels.o -L/home/fxr/OpenFOAM/OpenFOAM-v1806/platforms/linux64GccDPInt64Opt/lib \
    -lincompressibleTransportModels -lturbulenceModels -lfiniteVolume -lmeshTools -lincompressibleTurbulenceModels  -o /home/fxr/OpenFOAM/fxr-v1806/platforms/linux64GccDPInt64Opt/lib/test1_libIncompressibleTurbulenceModels.so

以上编译应该是没问题的，但是当我在算例中调用新的湍流模型后就报错：

Create time

--> FOAM Warning : 
    From function void* Foam::dlOpen(const Foam::fileName&, bool)
    in file POSIX.C at line 1516
    dlopen error : /home/fxr/OpenFOAM/fxr-v1806/platforms/linux64GccDPInt64Opt/lib/test1_libIncompressibleTurbulenceModels.so: undefined symbol: _ZNK4Foam9LESModels17kOmegaSSTIDDES_TrINS_29IncompressibleTurbulenceModelINS_14transportModelEEEE5GbyNuERKNS_16DimensionedFieldIdNS_7volMeshEEESA_SA_
--> FOAM Warning : 
    From function bool Foam::dlLibraryTable::open(const Foam::fileName&, bool)
    in file db/dynamicLibrary/dlLibraryTable/dlLibraryTable.C at line 109
    could not load "test1_libIncompressibleTurbulenceModels.so"
Create mesh for time = 0


PIMPLE: no residual control data found. Calculations will employ 5 corrector loops

Reading field p

Reading field U

Reading/calculating face flux field phi

Selecting incompressible transport model Newtonian
Selecting turbulence model type LES
[0] 
[0] 
[0] --> FOAM FATAL ERROR: 
[0] Unknown LESModel type banana

Valid LESModel types:

12
(
DeardorffDiffStress
Smagorinsky
SpalartAllmarasDDES
SpalartAllmarasDES
SpalartAllmarasIDDES
WALE
dynamicKEqn
dynamicLagrangian
kEqn
kOmegaSSTDDES
kOmegaSSTDES
kOmegaSSTIDDES
)

而且可以看到有效模型里面并没有出现自己修改的湍流模型，请问各位老师是什么原因？？？
make/files如下

turbulentTransportModels/myTurbulentTransportModels.C
LIB = $(FOAM_USER_LIBBIN)/test1_libIncompressibleTurbulenceModels

make/options如下

EXE_INC = \
-I../turbulenceModels/lnInclude \
-I$(LIB_SRC)/transportModels \
-I$(LIB_SRC)/finiteVolume/lnInclude \
-I$(LIB_SRC)/meshTools/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/turbulenceModels/lnInclude \
-I$(LIB_SRC)/TurbulenceModels/incompressible/lnInclude



LIB_LIBS = \
-lincompressibleTransportModels \
-lturbulenceModels \
-lfiniteVolume \
-lmeshTools \
-lincompressibleTurbulenceModels

用的openfoam1806，感谢

@星星星星晴 好的，感谢您的指点！

@星星星星晴 谢谢您的回复！
我打算在单相求解器上添加solidparticle实现粒子运动，还不涉及界面问题，所以选了pimpleFoam。
还有个问题想请教您，假如不考虑粒子的重力，那么该怎么调用粒子运动函数？
还是这样吗

particles.move(g);

根据chalmers教程在interFoam的基础上添加solidparticle可以编译通过，但是换成pimpleFoam之后就会报错：

pimpleFoam.C:99:24: error: 'g' was not declared in this scope
particles.move(g);

以下个人理解，这里particles是solidParticleCloud类的实例化对象，但是看了solidParticleCloud.H源码，g好像已经声明过了，不知道对不对，请大佬指点一下，谢谢！！！

public:
 
     // Constructors
 
         //- Construct given mesh
         solidParticleCloud
         (
             const fvMesh&,
             const word& cloudName = "defaultCloud",
             bool readFields = true
         );
 
 
     // Member Functions
 
         // Access
 
             bool hasWallImpactDistance() const;
 
             inline const fvMesh& mesh() const;
 
             inline scalar rhop() const;
             inline scalar e() const;
             inline scalar mu() const;
 
 
         // Edit
 
             //- Move the particles under the influence of the given
             //  gravitational acceleration
             void move(const dimensionedVector& g);

各位老师好！
报错信息：

solidParticleCloud.C的代码：

solidParticle::trackingData
        td(*this, rhoInterp, UInterp, nuInterp, g.value());

    Cloud<solidParticle>::move(*this, td, mesh_.time().deltaTValue());

    
    if(mesh_.time().value()> td.cloud().tInjStart_ &&                                   //
    mesh_.time().value()< td.cloud().tInjEnd_)
    {this->inject(td);}

    void Foam::solidParticleCloud::inject(solidParticle::trackingData &td)                  //.
    {                                                                                       //.
    label cellI = -1;                                                                       //.
    label tetFaceI = -1;
    label tetPtI = -1;
    mesh_.findCellFacePt(td.cloud().posB_, cellI, tetFaceI, tetPtI);
    if(cellI > 0){
    solidParticle* ptr1 = new solidParticle(mesh_, td.cloud().posB_,
    cellI, tetFaceI, tetPtI,td.cloud().dB_, td.cloud().UB_);
    Cloud<solidParticle>::addParticle(ptr1);
                 }
    }  

不应该是亚格子应力吗

@ibelief 多谢！

cfx计算用的SST?