SprayFoam 只喷固体该如何设置？

chpjz0391

0_1458720018133_New Microsoft PowerPoint Presentation.jpg
模型如图所示，我需要模拟一个辆车的尾气颗粒在一个隧道里怎么分布，颗粒与颗粒的碰撞需要考虑，这辆车只喷颗粒物质一定速度喷（比如每秒20000个），然后隧道一段以一定的速度进风。之前@散漫守望2016 给我介绍可以看看coalchemistryFoam，但是我今天看了sprayFoam，我想着用sprayFoam应该也可以解决这个问题。今天就花了很久在看sprayFoam 但是发现了一个问题，我可以把所有的化学反应关闭，但是在设置sprayCloudProperties的时候遇见了问题。我是基于 tutorials/lagrangian/sprayFoam/aachenBomb 这个case直接改的。这个case是喷入C7H16，我这里如果想只喷入固体颗粒，需要在singlePhaseMixtureCoeffs 下面的phases 里面写成 solid就可以了，然后再在 thermophysicalProperties 里面也改。但是运行的时候出现这个的时候出现 segmentation fault。

请问什么问题。

chpjz0391

我之前把例子中喷入的液体换成水倒是没什么问题？但是就是纯固体不行

李东岳

@chpjz0391 说：

前把例子中喷入的液体换成水倒是没什么问题？但

Hi，

好图！非常直观

或许把 segmentation fault这部分的log贴出来更方便诊断 :)

散漫守望2016

之前写的一个pimpleCoupledKinematicParcelFoam求解器，给你共享一下，大致下面这样，这个solver的主要部分，可以求解你的那个问题。
主程序：
int main(int argc, char *argv[])
{
argList::addOption
(
"cloudName",
"name",
"specify alternative cloud name. default is 'kinematicCloud'"
);
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "readGravitationalAcceleration.H"
#include "readTransportProperties.H"
#include "createFields.H"
#include "initContinuityErrs.H"
#include "readTimeControls.H"
#include "CourantNo.H"
#include "setInitialDeltaT.H"

pimpleControl pimple(mesh);

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Info<< "\nStarting time loop\n" << endl;

while (runTime.run())
{

    #include "readTimeControls.H"
    #include "CourantNo.H"
    #include "setDeltaT.H"

    runTime++;

    Info<< "Time = " << runTime.timeName() << nl << endl;

    Info<< "Evolving " << kinematicCloud.name() << endl;
    kinematicCloud.evolve();
    kinematicCloud.info();

    // --- Pressure-velocity PIMPLE corrector loop
    while (pimple.loop())
    {
        #include "UEqn.H"

        // --- Pressure corrector loop
        while (pimple.correct())
        {
            #include "pEqn.H"
        }

        if (pimple.turbCorr())
        {
            turbulence->correct();
        }
    }
    
    runTime.write();

    Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
        << "  ClockTime = " << runTime.elapsedClockTime() << " s"
        << nl << endl;
}

Info<< "End\n" << endl;

return 0;

}

CreateFields.H文件
Info<< "Reading normalised dynamic presure field p_rgh\n" << endl;
volScalarField p_rgh
(
IOobject
(
"p_rgh",
runTime.timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
);

Info<< "Reading field U\n" << endl;
volVectorField U
(
    IOobject
    (
        "U",
        runTime.timeName(),
        mesh,
        IOobject::MUST_READ,
        IOobject::AUTO_WRITE
    ),
    mesh
);

Info<< "Creating phi [m3/s] " << endl;
#include "createPhi.H"

Info<< "Creating turbulence model\n" << endl;
singlePhaseTransportModel laminarTransport(U, phi);
autoPtr<incompressible::turbulenceModel> turbulence
(
    incompressible::turbulenceModel::New(U, phi, laminarTransport)
);

//Info<< "Creating source list\n" << endl;
fv::IOoptionList fvOptions(mesh);


Info<< "Calculating field g.h\n" << endl;
volScalarField gh("gh", g & mesh.C());
surfaceScalarField ghf("ghf", g & mesh.Cf());

Info<< "Creating absolute pressure field p [Pa]\n" << endl;
volScalarField p
(
    IOobject
    (
        "p",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    rhoVal*(p_rgh + gh)
);

Info<< "Reading field U\n" << endl;

label pRefCell = 0;
scalar pRefValue = 0.0;
setRefCell
(
    p,
    p_rgh,
    mesh.solutionDict().subDict("PIMPLE"),
    pRefCell,
    pRefValue
);

if (p_rgh.needReference())
{
    p += dimensionedScalar
    (
        "p",
        p.dimensions(),
        pRefValue - getRefCellValue(p, pRefCell)
    );
}

Info<< "Creating mass flux density surface field rhoPhi [kg/s] " << endl;
surfaceScalarField rhoPhi
(
    IOobject
    (
        "rhoPhi",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    rhoVal*linearInterpolate(U) & mesh.Sf()
);

// dummy-fields to satisfy the requirements of the kinematicCloud
Info<< "Creating rho field " << endl;
volScalarField rho
(
    IOobject
    (
        "rho",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    rhoVal,
    "zeroGradient"
);

Info<< "Creating mu field " << endl;
volScalarField mu
(
    IOobject
    (
        "mu",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::NO_WRITE
    ),
    mesh,
    rhoVal*nu,
    "zeroGradient"
);

Info<< "Constructing kinematicCloud" << endl;
basicKinematicCollidingCloud kinematicCloud
(
    "kinematicCloud",
    rho,
    U,
    mu,
    g
);

**
readTransportProperties文件**：
Info<< "\nReading transportProperties\n" << endl;

IOdictionary transportProperties
(
    IOobject
    (
        "transportProperties",
        runTime.constant(),
        mesh,
        IOobject::MUST_READ_IF_MODIFIED,
        IOobject::NO_WRITE
    )
);

dimensionedScalar nu
(
    transportProperties.lookup("nu")
);

dimensionedScalar rhoVal
(
    transportProperties.lookup("rho")
);
Info<< "rhoVal is " << rhoVal << endl;

dimensionedScalar irhoVal("irhoVal",(1.0/rhoVal));

**
UEqn.H文件：**
// Solve the momentum equation
fvVectorMatrix UEqn
(
fvm::ddt(U)

fvm::div(phi, U)
turbulence->divDevReff(U)
==
irhoVal*kinematicCloud.SU(U)
fvOptions(U)
);

UEqn.relax();

fvOptions.constrain(UEqn);

if (pimple.momentumPredictor())
{
solve
(
UEqn
==
fvc::reconstruct
(
-fvc::snGrad(p_rgh)*mesh.magSf()
)
);
}

pEqn.H文件：
{
volScalarField rAU("rAU", 1.0/UEqn.A());
surfaceScalarField rAUf("(1|A(U))", fvc::interpolate(rAU));

U = rAU*UEqn.H();

phi = (fvc::interpolate(U) & mesh.Sf());
//+ fvc::ddtCorr(U, phi);

adjustPhi(phi, U, p_rgh);

while (pimple.correctNonOrthogonal())
{
fvScalarMatrix p_rghEqn
(
fvm::laplacian(rAU, p_rgh) == fvc::div(phi)
);

  p_rghEqn.setReference(pRefCell, getRefCellValue(p_rgh, pRefCell));

  p_rghEqn.solve(mesh.solver(p_rgh.select(pimple.finalInnerIter())));

  if (pimple.finalNonOrthogonalIter())
  {
      // Calculate the conservative fluxes
      phi -= p_rghEqn.flux();

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

      // Correct the momentum source with the pressure gradient flux
      // calculated from the relaxed pressure
      U -= rAU*fvc::reconstruct((p_rghEqn.flux())/rAUf);
      U.correctBoundaryConditions();
  }

}

#include "continuityErrs.H"

p = rhoVal*(p_rgh + gh);

if (p_rgh.needReference())
{
p += dimensionedScalar
(
"p",
p.dimensions(),
pRefValue - getRefCellValue(p, pRefCell)
);
p_rgh = irhoVal*p - gh;
}
}

chpjz0391

这几天放假没在学校，才看到回复。谢谢分享。我研究一下。万分感谢