CFD中文网

@Afr1yne 好的好的非常感谢 我已经找到了

各位好
我在sprayCloudProperties中添加了以下语句，想要调用CloudFunctionObjects功能，但我在计算结果中没有看到相关变量的输出，这是什么原因呢？

cloudFunctions
{
    VoidFraction1
    {
        type voidFraction;
    }
}

求助各位大佬们，我计算了一个冷态的流场，中心氮气携带颗粒，周围为旋转流入的氧气。在计算一段时间后，我想用setFields进行点火。一开始尝试给到1800K，但很快发散了。目前尝试先加到600K。
但是在测试中，会出现奇怪的升温（关闭了combustion、chemistry）。之后会一直升温，直到温度超出物理限制后报错。测试了很多方面也不知道是为什么。我观测到的是，这个升温的区域和速度场有一些近似，但不知道有没有什么关系。求助各位老师同学们指点一下迷津，实在是想不明白了ToT。

冷态流场速度图：

全场升到600K后，计算一步的温度图（入口和壁面温度边界都给的300K；出口给的zeroGradient但出口在很下游，图中并未截取到）：

相同的设置在本地虚拟机上设置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)

已解决