SprayFoam 只喷固体该如何设置?



  • 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。

    请问什么问题。



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


  • 网格教授 OpenFOAM教授 管理员

    @chpjz0391 说:

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

    Hi,

    好图!非常直观

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


  • 离散相副教授

    之前写的一个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;
    }
    }



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