CFD中文网

相同的设置在本地虚拟机上设置unrefineLevel是可以运行的，但是在服务器上会报错，很迷惑。

@zhangzh 您好请问您后面有解决动网格的问题嘛

各位老师好，我在使用自适应网格加密的时候当设置unrefineLevel阈值的时候会出现以下的报错（我的求解器耦合了cantera）

当我设置unrefineLevel阈值为0时，则不会出现上述问题。
想问一下大家有遇到过类似的问题吗

@李东岳

@Rachel0096 在 请教贴：拉格朗日耦合进密度基求解器，能量方程报错 中说：

AlParcels.Sh(e)

这个是自己写的铝颗粒的能量源项，类似于coalChemistryFoam中的coalParcels.Sh(he)和limestoneParcels.Sh(he)，表示的是颗粒的能量源项。

@李东岳

@Rachel0096 在 请教贴：拉格朗日耦合进密度基求解器，能量方程报错 中说：

• AlParcels.Sh(e)
• radiation->Sh(thermo, e)

李老师，是这个位置有问题，这两行注释了之后就求解器就可以运行了，我尝试过把e更改为rhoE也会报错

各位老师们好，目前我在OpenFOAM-v1912中尝试把一个压力基的拉格朗日两相反应求解器（基于coalChemstryFoam）耦合进detonationFoam（密度基）中。目前编译没有问题，可以顺利通过，但是在测试算例时运行会报错（使用的是NS_Sutherland模型）：

=== Start Sensible Enthalpy Transport ====


--> FOAM FATAL ERROR:
incompatible fields for operation
    [rhoE] - [h]

    From function void Foam::checkMethod(const Foam::fvMatrix<Type>&, const Foam::fvMatrix<Type>&, const char*) [with Type = double]
    in file /share/home/zhoufan/OpenFOAM/OpenFOAM-v1912/src/finiteVolume/lnInclude/fvMatrix.C at line 1337.

FOAM aborting

#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::error::abort() at ??:?
#2  void Foam::checkMethod<double>(Foam::fvMatrix<double> const&, Foam::fvMatrix<double> const&, char const*) at ??:?
#3  Foam::tmp<Foam::fvMatrix<double> > Foam::operator-<double>(Foam::tmp<Foam::fvMatrix<double> > const&, Foam::tmp<Foam::fvMatrix<double> > const&) at ??:?
#4  ? at ??:?
#5  __libc_start_main in /lib64/libc.so.6
#6  ? at ??:?
Aborted (core dumped)

这里显示动量方程和组分方程求解都没有问题，在求解能量方程时遇到了赋值不统一的问题。
这是原本detonationFoam中rhoEEqn的代码：

surfaceScalarField sigmaDotU
(
    "sigmaDotU",
    (
        fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
      + fvc::dotInterpolate(mesh.Sf(), tauMC)
    )
  & fvc::interpolate(U)    
);

solve
(
    fvm::ddt(rhoE)
  + fvc::div(rhoEPhi)
  - fvc::div(sigmaDotU)
  ==
    reaction->Qdot()
);

e = rhoE/rho - 0.5*magSqr(U);
e.correctBoundaryConditions();
thermo.correct();
rhoE.boundaryFieldRef() == rho.boundaryField()*(e.boundaryField()+0.5*magSqr(U.boundaryField()));


solve
(
    fvm::ddt(rho, e) - fvc::ddt(rho, e)
  + thermophysicalTransport->divq(e)
);
thermo.correct();
rhoE = rho*(e + 0.5*magSqr(U));

p.ref() = rho()/psi();
p.correctBoundaryConditions();
rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField();

Info<< "min/max(p) = "<< min(p).value() << ", " << max(p).value() << endl;
Info<< "min/max(T) = "<< min(T).value() << ", " << max(T).value() << endl;

这是我修改添加颗粒源项之后的代码：

clock_t t_begin_h = std::clock();
Info <<  "=== Start Sensible Enthalpy Transport ====" << endl;
auto dt = runTime.deltaT();

surfaceScalarField sigmaDotU
(
    "sigmaDotU",
    (
        fvc::interpolate(muEff)*mesh.magSf()*fvc::snGrad(U)
      + fvc::dotInterpolate(mesh.Sf(), tauMC)
    )
  & fvc::interpolate(U)    
);

solve
(
    fvm::ddt(rhoE)
  + fvc::div(rhoEPhi)
  - fvc::div(sigmaDotU)
  - fvm::laplacian(turbulence->mut()/Prt + turbulence->alpha(), e)
  ==
    combustion->Qdot()
  + AlParcels.Sh(e)
  + radiation->Sh(thermo, e)
);

e = rhoE/rho - 0.5*magSqr(U);
e.correctBoundaryConditions();
thermo.correct();
rhoE.boundaryFieldRef() == rho.boundaryField()*(e.boundaryField()+0.5*magSqr(U.boundaryField()));

solve
(
    fvm::ddt(rho, e) - fvc::ddt(rho, e)
);

thermo.correct();
radiation->correct(); //add by vv
rhoE = rho*(e + 0.5*magSqr(U));

p.ref() = rho()/psi();
p.correctBoundaryConditions();
rho.boundaryFieldRef() == psi.boundaryField()*p.boundaryField();

Info<< "min/max(p) = "<< min(p).value() << ", " << max(p).value() << endl;
Info<< "min/max(T) = "<< min(T).value() << ", " << max(T).value() << endl;

clock_t t_end_h = std::clock();
double elapsed_secs_h = double(t_end_h - t_begin_h) / CLOCKS_PER_SEC;

Info <<  "=== Sensible Enthalpy Transport [Done]: "
     << elapsed_secs_h
     << " s ====\n"
     << endl;

想请教一下大家有没有什么更改的建议。

@欧阳 上交大的xuleilei和中科院zhangjunqing也基于openfoam做了类似的工作

或许可以在createClouds.H中声明多个basicSprayCloud对象，并将相应的源项添加进去

感谢老师们！这里确实是多变量的计算，类似于把单个组分场拆分为很多不同的场来计算，但是输入相同的源项，所以会有两种不同的Y在同一个方程里。程序就是在方程这一步出现了错误，我试试各位老师的方法，多输出一些内容看看哪里出问题了。

@ice_flow 您好，我也遇到了类似的问题，可以向您请教吗？这是对openfoam2.3.x版本中coalChemistryFoam的改写。大致也是在一个方程中，有两个不同的变量，编译可以通过，但在计算式出现报错incompatible fields for operation [N2] == [N2]。个人感觉也是同样的问题，但不是很清楚要如何去更改。YiEqn部分代码如下：

    for(label i=0; i<Y.size(); i++)
    {
        {
            // volScalarField Yi   = Ysf[nsf][i];
            volScalarField Ymean = Y[i];

            volScalarField tempYi
            (
                IOobject
                (
                    Y[i].name(),
                    runTime.timeName(),
                    mesh,
                    IOobject::NO_READ,
                    IOobject::NO_WRITE
                ),
                0*Y[0],
                "zeroGradient"
            );

            tempYi = Ysf[nsf][i];

            WienerTermYi = pow(2*gamma,0.5)*( fvc::grad(tempYi) & dW_sf);
            fvScalarMatrix YiEqn
            (
                  fvm::ddt(rho, tempYi)
                + mvConvection->fvmDiv(phi, tempYi)
                - fvm::laplacian(turbulence->muEff()/Sct, tempYi)
                ==
                  coalParcels.SYi(i, Ymean)
                + combustion->R(Ymean)
                + fvOptions(rho, Ymean)
		        + (rho*WienerTermYi/deltaT)
            );
            YiEqn.relax();
            fvOptions.constrain(YiEqn);
            YiEqn.solve(mesh.solver("tempYi"));
            YiEqn.solve();
            fvOptions.correct(tempYi);

            tempYi.max(0.0);
	        Ysf[nsf][i] = tempYi;
            Yt += tempYi;
        }
    }

我在对coalChemistryFoam改写，想把它和随机场方程进行结合。简单来说是对于连续相用一系列随机场进行计算，颗粒相不用（单独计算）。但这就遇到了一个问题，我结合着输运方程来说明：

    for(label i=0; i<Y.size(); i++)
    {
        {
            volScalarField Yi   = Ysf[nsf][i];
            volScalarField Ymean = Y[i];

            WienerTermYi = pow(2*gamma,0.5)*( fvc::grad(Yi) & dW_sf);
            fvScalarMatrix YiEqn
            (
                  fvm::ddt(rho, Yi)
                + mvConvection->fvmDiv(phi, Yi)
                - fvm::laplacian(turbulence->muEff()/Sct, Yi)
                ==
                  coalParcels.SYi(i, Ymean)
                + combustion->R(Ymean)
                + fvOptions(rho, Ymean)
		        + (rho*WienerTermYi/deltaT)
            );
            YiEqn.relax();
            fvOptions.constrain(YiEqn);
            YiEqn.solve(mesh.solver("Yi"));
            YiEqn.solve();
            fvOptions.correct(Yi);

            Yi.max(0.0);
	        Ysf[nsf][i] = Yi;
            Yt += Yi;
        }
    }

在这个方程里，Yi指的是随机场的组分项（比如名为N2__01），Ymean指的是平均场的组分项（也就是组分场，名为N2）。颗粒的求解不考虑随机场，来自颗粒的源项被输入到每一个随机场中。但这里就会有一个问题，在计算时，会出现类似incompatible fields for operation [N2] + [N2__01]的报错，我的理解是因为这两个名称不同所导致的。这个问题要怎么去解决？我想的是，或许可以把这三项在方程外部导入到一个临时变量，再输入到方程里，但是不知道怎么用代码去实现。求各位大佬指点。

@李东岳 感谢李老师！原来没有特意写出来的都是同一个

我使用的是openfoam2.3.x，其reactingFoam目录是这样的：

其中没有rhoEqn.H，但是代码里有rhoEqn.H：

这样的话，我要怎么去找，它这里引用的rhoEqn.H在openfoam中的哪里？

@xuqiming 论坛里有个名为“【分享+搬运】自定义非均匀inlet U”的帖子，或许可以帮到你

是想用cantera算一维火焰做FGM表，用在openfoam计算里ToT

我参考cantera官网示例中的diffusion_flame_batch.py，设置一批应变率变化时的一维对冲火焰解。通过反复调用前一计算结果作为初始解，来逐渐增加或降低应变率，从而得到计算结果。当我按照代码设置应变率上升时（即和源代码中一样设置strain_factor = 1.25），边界结果无明显问题。但当我设置应变率下降时（将strain_factor改为0.8），在燃料侧的grid=0处，其结果中的组分组成随应变率下降而逐渐偏离边界条件。我的燃料设置为Y_O2:0.2，Y_NH3:0.1417，Y_N2:0.6583。取应变率下降循环第50次的结果而言，grid=0处的O2仅为0.15左右，NH3也偏差很大。这个可能是什么原因引起的？

示例中的源代码部分如下：

# PART 3: STRAIN RATE LOOP

# Compute counterflow diffusion flames at increasing strain rates at 1 bar
# The strain rate is assumed to increase by 25% in each step until the flame is
# extinguished
strain_factor = 1.25

# Exponents for the initial solution variation with changes in strain rate
# Taken from Fiala and Sattelmayer (2014)
exp_d_a = - 1. / 2.
exp_u_a = 1. / 2.
exp_V_a = 1.
exp_lam_a = 2.
exp_mdot_a = 1. / 2.

# Restore initial solution
file_name, entry = names("initial-solution")
f.restore(file_name, name=entry)

# Counter to identify the loop
n = 0
# Do the strain rate loop
while np.max(f.T) > temperature_limit_extinction:
    n += 1
    print('strain rate iteration', n)
    # Create an initial guess based on the previous solution
    # Update grid
    f.flame.grid *= strain_factor ** exp_d_a
    normalized_grid = f.grid / (f.grid[-1] - f.grid[0])
    # Update mass fluxes
    f.fuel_inlet.mdot *= strain_factor ** exp_mdot_a
    f.oxidizer_inlet.mdot *= strain_factor ** exp_mdot_a
    # Update velocities
    f.set_profile('velocity', normalized_grid,
                  f.velocity * strain_factor ** exp_u_a)
    f.set_profile('spread_rate', normalized_grid,
                  f.spread_rate * strain_factor ** exp_V_a)
    # Update pressure curvature
    f.set_profile('lambda', normalized_grid, f.L * strain_factor ** exp_lam_a)
    try:
        # Try solving the flame
        f.solve(loglevel=0)

已解决

当我使用原本的网格和0设置时，仅更换为新的机理（gri3.0，原机理是一个氢气的机理），此时会出现一些“attempt to use janafThermo<EquationOfState> out of temperature range 300 -> 3000; T = 298K”的warning，但是可以正常进行计算。

我在使用一个他人编写的of2.3.x上的求解器，使用其随带发来的算例可以正常运行。之后，主要对0、polymesh和机理进行了更换，尝试计算自己的算例，但是在开始计算时出现报错，看上去是在读取机理时出现问题，想请教各位老师，问题是出在哪里？

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Create time

Create mesh for time = 0


Reading g
Creating reaction model

Selecting chemistry type 
{
    chemistrySolver EulerImplicit;
    chemistryThermo psi;
}

Selecting thermodynamics package 
{
    type            hePsiThermo;
    mixture         reactingMixture;
    transport       sutherland;
    thermo          janaf;
    energy          sensibleEnthalpy;
    equationOfState perfectGas;
    specie          specie;
}

Selecting chemistryReader foamChemistryReader
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2   in "/lib64/libc.so.6"
#3  Foam::DimensionedField<double, Foam::volMesh>::operator/=(Foam::DimensionedField<double, Foam::volMesh> const&) at ??:?
#4  Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>::operator/=(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&) at ??:?
#5  Foam::multiComponentMixture<Foam::sutherlandTransport<Foam::species::thermo<Foam::janafThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy> > >::correctMassFractions() at ??:?
#6  Foam::multiComponentMixture<Foam::sutherlandTransport<Foam::species::thermo<Foam::janafThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy> > >::multiComponentMixture(Foam::dictionary const&, Foam::List<Foam::word> const&, Foam::HashPtrTable<Foam::sutherlandTransport<Foam::species::thermo<Foam::janafThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy> >, Foam::word, Foam::string::hash> const&, Foam::fvMesh const&) at ??:?
#7  Foam::reactingMixture<Foam::sutherlandTransport<Foam::species::thermo<Foam::janafThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy> > >::reactingMixture(Foam::dictionary const&, Foam::fvMesh const&) at ??:?
#8  Foam::heThermo<Foam::psiReactionThermo, Foam::SpecieMixture<Foam::reactingMixture<Foam::sutherlandTransport<Foam::species::thermo<Foam::janafThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy> > > > >::heThermo(Foam::fvMesh const&, Foam::word const&) at ??:?
#9  Foam::psiReactionThermo::addfvMeshConstructorToTable<Foam::hePsiThermo<Foam::psiReactionThermo, Foam::SpecieMixture<Foam::reactingMixture<Foam::sutherlandTransport<Foam::species::thermo<Foam::janafThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy> > > > > >::New(Foam::fvMesh const&, Foam::word const&) at ??:?
#10  Foam::autoPtr<Foam::psiReactionThermo> Foam::basicThermo::New<Foam::psiReactionThermo>(Foam::fvMesh const&, Foam::word const&) at ??:?
#11  Foam::psiReactionThermo::New(Foam::fvMesh const&, Foam::word const&) at ??:?
#12  Foam::psiChemistryModel::psiChemistryModel(Foam::fvMesh const&) at ??:?
#13  Foam::chemistryModel<Foam::psiChemistryModel, Foam::sutherlandTransport<Foam::species::thermo<Foam::janafThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy> > >::chemistryModel(Foam::fvMesh const&) at ??:?
#14  Foam::EulerImplicit<Foam::chemistryModel<Foam::psiChemistryModel, Foam::sutherlandTransport<Foam::species::thermo<Foam::janafThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy> > > >::EulerImplicit(Foam::fvMesh const&) at ??:?
#15  Foam::psiChemistryModel::addfvMeshConstructorToTable<Foam::EulerImplicit<Foam::chemistryModel<Foam::psiChemistryModel, Foam::sutherlandTransport<Foam::species::thermo<Foam::janafThermo<Foam::perfectGas<Foam::specie> >, Foam::sensibleEnthalpy> > > > >::New(Foam::fvMesh const&) at ??:?
#16  Foam::autoPtr<Foam::psiChemistryModel> Foam::basicChemistryModel::New<Foam::psiChemistryModel>(Foam::fvMesh const&) at ??:?
#17  Foam::psiChemistryModel::New(Foam::fvMesh const&) at ??:?
#18  
 at ??:?
#19  __libc_start_main in "/lib64/libc.so.6"
#20  
 at ??:?
Floating point exception (core dumped)

或许是我的表述不够清晰，我想模拟图一所示的构型，两股不同的喷雾分布从不同的喷雾入口注入，我看到之前DieselFoam(或许是sprayFoam的前身)可以通过指定不同的X来喷注两种液体，这应该在sprayFoam中如何实现呢？