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