[DustExplosion]

Hello All,

I am sure there is an answer to this on the forum but I could not get the correct keywords to find it. I am starting to setup a multiphase, reactive simulation with the coalChemistryFoam solver (OF 3.0.1). However; at first I just want to setup a 1D tube, closed at both ends, full of 100% CH4 (no particles). The goal is to have constant pressure, zero velocity, and the simulation just sit there and do nothing (static).

When I run the simulation all of the fields remain constant except the fluid velocity field, which has small perturbations. I am running the simulation 10 timesteps, and the perturbations come after the first one (see attached image).

The problem goes away when I comment out both

solve(UEqn == -fvc::grad(p)); in UEqn.H (in momentum predictor)

and

U = HbyA - rAU*fvc::grad(p); in pEqn.H (in momentum corrector)

I tried turning everything off including the particles, radiation, turbulence, and combustion, but the problem still persists. I think the mesh and boundaries are defined correctly (see both below).

At this point, I am stuck to as what the problem could be. The p-field in the output seems to be constant so I think grad(p) is zero, and the issue is with the HbyA term. However, I am not sure how to track it down further. Has anyone seen similar problems with the Lagrange or combustion solvers, or can point me to a tutorial demonstrating setup of a "static" simulation? I tried uploading the test case as an archive, but it wouldn't work (it is only 25 kb). What is the proper way to upload the case?

Thank you in advance for any help!

Chris

CH4
dimensions [0 0 0 0 0 0 0];
internalField uniform 1.00;
boundaryField
{
top
{
type zeroGradient;
}

bottom
{
type zeroGradient;
}

frontAndBack
{
type empty;
}
}

U
dimensions [0 1 -1 0 0 0 0];
internalField uniform (0 0 0);
boundaryField
{
top
{
type fixedValue;
value uniform (0 0 0);
}
bottom
{
type fixedValue;
value uniform (0 0 0);
}

frontAndBack
{
type empty;
}
}

p
dimensions [1 -1 -2 0 0 0 0];
internalField uniform 100000;
boundaryField
{
top
{
type zeroGradient;
}
bottom
{
type zeroGradient;
}

frontAndBack
{
type empty;
}
}

Mesh
convertToMeters 1.0;
vertices
(
(0 0 0)
(0.03 0 0)
(0.03 0.002 0)
(0 0.002 0)
(0 0 0.002)
(0.03 0 0.002)
(0.03 0.002 0.002)
(0 0.002 0.002)
);

blocks
(
hex (0 1 2 3 4 5 6 7) (300 1 1) simpleGrading (1 1 1)
);

edges
(
);

boundary
(
top
{
type wall;
faces
(
(4 7 3 0)
);
}

bottom
{
type wall;
faces
(
(1 2 6 5)
);
}

frontAndBack
{
type empty;
faces
(
(0 3 2 1)
(4 5 6 7)
(7 6 2 3)
(1 0 4 5)
);
}
);

[DustExplosion]

Hello,

I have been running an academic conference this week so I have not had a chance to get back to this.

I was wondering if anyone would be able to provide a hint to point me in the correct direction or an example of a "static" case they have setup?

Thanks again,

Chris

[DustExplosion]

Hello Again,

I realized that you can upload tarballs, so my test case is attached.

Thank you again if anyone can provide any help!

Chris

[KarenRei]

Originally Posted by DustExplosion

Hello Again,

I realized that you can upload tarballs, so my test case is attached.

Thank you again if anyone can provide any help!

Chris

Have you run checkMesh? Are you using a stable interpolation scheme in fvSchemes?

[DustExplosion]

CheckMesh mesh comes back with everything ok and the faces and volumes seem correct (although they are small as my domain is only 30 mm).

I think the schemes are all pretty standard. Do any of these look like they could cause problems?

ddtSchemes
{
default Euler;
}

gradSchemes
{
default Gauss linear;
}

divSchemes
{
default none;

div(phi,U) Gauss upwind;
div(phid,p) Gauss upwind;
div(phi,K) Gauss linear;
div(phi,h) Gauss upwind;
div(phi,k) Gauss upwind;
div(phi,epsilon) Gauss upwind;
div(U) Gauss linear;
div(((rho*nuEff)*dev2(T(grad(U))))) Gauss linear;
div(phi,Yi_h) Gauss upwind;
}

laplacianSchemes
{
default Gauss linear orthogonal;
}

interpolationSchemes
{
default linear;
}

snGradSchemes
{
default orthogonal;
}

[DustExplosion]

I started trying to set up the problem using different solvers and was able to reproduce something with a similar response with the incompressible solver icoFoam. If I take the cavity tutorial case and do the following I get the result attached:
1) Change the lid velocity to zero
2) Set the uniform pressure field to 100000
3) Set the pRefValue in fvSolution to 100000

I am not sure if this is related to my original problem, but there looks like there may be a precision issue when pRefValue is set to a large number and the flow is static.

Does anyone have any thoughts on this, or if it is related to my static simulation issue with coalChemistryFoam solver using pimple?

[DustExplosion]

Just a follow up on this.

The issue with the cavity case happened for pimpleFoam as well and appears to be related to pRefValue. I found a post here (http://www.cfd-online.com/Forums/ope...prefvalue.html) that describes the use of pRefValue and suggests that pRefCell cannot be located on the boundary (however, when I set is to the centre of the grid I was seeing the same issue).

All that being said, I do not think that this is the same issue I am seeing in coalChemistryFoam as it is compressible and does not require a reference pressure to be set. I made my domain open (inflow on one side, outflow on the other) with a constant velocity field and was seeing the same behaviour (fluctuating velocity solution).

I also noticed that the solvers for velocity, enthalpy, and pressure were reacting max iterations in this case. This occurred for both smoothSolver and GAMG. However, the preconditioning solver PBiCg was able to converge in one step.

I am still not sure why it is solving anything at all, as the other cases and incompressible solvers seem to take zero iterations when the flow is static. However, I am going to move ahead with the preconditioning solver until I can figure it out as it does not seem to diverge at least.