Source Term on Scalar Transport Equation - crashing simulation
 

Hello,

I have developed a scalar transport equation for a passive scalar + a source term in the equation. Unfortunately in the simulation I have some crashes.

I explain what I have done:

The equation is a tipical parabolic equation + a source term in the form of ST= C*(1-C)*a

where "C"= scalar variable, "a"= constant

creating a new file for the equation and attaching in the solver (editing simpleFoam), I obtain:
solve
(
fvm::ddt(C)
+fvm::div(phi, C)
-fvm::laplacian(nu, C)==
((1-C)*(C))*a*runTime.deltaT().value()/runTime.deltaT()
);


the "a" is declared on "createFields.H" as:

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

and on the case file on "constant/transportProperties" I add the following:

a a [ 0 0 0 0 0 0 0 ] 1;

--------------------
CRASHING PROBLEM: On case file the C b.c. are setted as:
1 in the inlet,
inletOutlet for outlet
zeroGradient in all walls.

When i change the value of "a" in the file and I increase up to 5 or more I obtain, after a few iterations, the following message:

-----------------------------------------------------

Time = 9

DILUPBiCG: Solving for Ux, Initial residual = 0.246096, Final residual = 0.0119745, No Iteration s 1
DILUPBiCG: Solving for Uy, Initial residual = 0.252583, Final residual = 0.0122985, No Iteration s 1
DILUPBiCG: Solving for Uz, Initial residual = 0.0510667, Final residual = 0.00160401, No Iterati ons 1
DICPCG: Solving for p, Initial residual = 0.716255, Final residual = 0.0714609, No Iterations 39
time step continuity errors : sum local = 194.423, global = -29.0713, cumulative = 39.8563
DILUPBiCG: Solving for C, Initial residual = 1, Final residual = 0.0981156, No Iterations 99
DILUPBiCG: Solving for epsilon, Initial residual = 0.00528578, Final residual = 0.000203886, No Iterations 1
DILUPBiCG: Solving for k, Initial residual = 0.0121429, Final residual = 0.000453808, No Iterati ons 1
ExecutionTime = 33.9 s ClockTime = 34 s

Time = 10

DILUPBiCG: Solving for Ux, Initial residual = 0.150043, Final residual = 0.00577888, No Iteratio ns 1
DILUPBiCG: Solving for Uy, Initial residual = 0.153792, Final residual = 0.00587662, No Iteratio ns 1
DILUPBiCG: Solving for Uz, Initial residual = 0.0578676, Final residual = 0.00175246, No Iterati ons 1
DICPCG: Solving for p, Initial residual = 0.625491, Final residual = 0.0608153, No Iterations 16 8
time step continuity errors : sum local = 115.09, global = 0.500666, cumulative = 40.3569
DILUPBiCG: Solving for C, Initial residual = 1, Final residual = 0.0858237, No Iterations 66
DILUPBiCG: Solving for epsilon, Initial residual = 0.00469556, Final residual = 0.000175706, No Iterations 1
DILUPBiCG: Solving for k, Initial residual = 0.0108158, Final residual = 0.000394342, No Iterati ons 1
ExecutionTime = 38.26 s ClockTime = 38 s

Time = 11

DILUPBiCG: Solving for Ux, Initial residual = 0.188378, Final residual = 0.00755333, No Iteratio ns 1
DILUPBiCG: Solving for Uy, Initial residual = 0.194324, Final residual = 0.00780845, No Iteratio ns 1
DILUPBiCG: Solving for Uz, Initial residual = 0.0400986, Final residual = 0.00119253, No Iterati ons 1
DICPCG: Solving for p, Initial residual = 0.235323, Final residual = 0.0225155, No Iterations 95
time step continuity errors : sum local = 83.1461, global = -4.27269, cumulative = 36.0842
[6] #0 Foam::error::printStack(Foam::Ostream&)[4] #0 [7] Foam::error::printStack(Foam::Ostream& )#0 [5] Foam::error::printStack(Foam::Ostream&)#0 Foam::error::printStack(Foam::Ostream&)[1] #0 Foam::error::printStack(Foam::Ostream&)[0] [3] #0#0 Foam::error::printStack(Foam::Ostream&)F oam::error::printStack(Foam::Ostream&)[2] #0 Foam::error::printStack(Foam::Ostream&) in "/cvos/s hared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"
[1] #1 Foam::sigFpe::sigFpeHandler(int) in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux6 4GccDPOpt/libOpenFOAM.so"
[6] #1 Foam::sigFpe::sigFpeHandler(int) in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux6 4GccDPOpt/libOpenFOAM.so"
[5] #1 Foam::sigFpe::sigFpeHandler(int) in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux6 4GccDPOpt/libOpenFOAM.so"
[3] #1 Foam::sigFpe::sigFpeHandler(int) in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux6 4GccDPOpt/libOpenFOAM.so"
[0] #1 Foam::sigFpe::sigFpeHandler(int) in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux6 4GccDPOpt/libOpenFOAM.so"
[4] #1 Foam::sigFpe::sigFpeHandler(int) in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux6 4GccDPOpt/libOpenFOAM.so"
[2] #1 Foam::sigFpe::sigFpeHandler(int) in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux6 4GccDPOpt/libOpenFOAM.so"
[7] #1 Foam::sigFpe::sigFpeHandler(int) in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux6 4GccDPOpt/libOpenFOAM.so"
[1] #2 in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"
[5] #2 in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"
[3] #2 in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"
[0] #2 in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"
[6] #2 in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"
[7] #2 in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"
[4] #2 __restore_rt in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOA M.so"
[2] #2 __restore_rt__restore_rt__restore_rt__restore_rt__ restore_rt__restore_rt at sigaction.c:0
[1] #3 Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const __restore_rt at sigaction.c:0
[3] #3 Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at sigaction.c:0
[5] #3 Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at sigaction.c:0
[0] #3 Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at sigaction.c:0
[6] #3 Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const at sigaction.c:0
[7] #3 Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"
[3] #4 Foam::fvMatrix<double>::solve(Foam::dictionary const&) at sigaction.c:0
[4] #3 Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"
[5] #4 Foam::fvMatrix<double>::solve(Foam::dictionary const&) at sigaction.c:0
[2] #3 Foam::PBiCG::solve(Foam::Field<double>&, Foam::Field<double> const&, unsigned char) const in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"
[1] #4 Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/cvos/shared/apps/OpenFOAM/Ope nFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"
[0] #4 Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/cvos/shared/apps/OpenFOAM/Ope nFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"
[6] #4 Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/cvos/shared/apps/OpenFOAM/Ope nFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"
[7] #4 Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/cvos/shared/apps/OpenFOAM/Ope nFOAM-1.7.1/lib/linux64GccDPOpt/libfiniteVolume.so"
[3] #5 in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"
[4] #4 Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/cvos/shared/apps/OpenFOAM/Ope nFOAM-1.7.1/lib/linux64GccDPOpt/libfiniteVolume.so"
[5] #5 in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libfiniteVolume.so"
[1] #5 in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so"
[2] #4 Foam::fvMatrix<double>::solve(Foam::dictionary const&) in "/cvos/shared/apps/OpenFOAM/Ope nFOAM-1.7.1/lib/linux64GccDPOpt/libfiniteVolume.so"
[7] #5 in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libfiniteVolume.so"
[0] #5 in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libfiniteVolume.so"
[6] #5 in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libfiniteVolume.so"
[4] #5 in "/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libfiniteVolume.so"
[2] #5 mainmainmain in "/home/afancello/OpenFOAM/afancello-1.7.1/applications/bin/linux64GccDPOp t/SourceSimple"
[3] #6 __libc_start_mainmainmainmain in "/home/afancello/OpenFOAM/user-1.7.1/applications/b in/linux64GccDPOpt/SourceSimple"
[1] #6 __libc_start_mainmainmain in "/lib64/libc.so.6"
[3] #7 in "/home/user/OpenFOAM/user-1.7.1/applications/bin/linux64GccDPOpt/Source Simple"
[5] #6 __libc_start_main in "/home/user/OpenFOAM/user-1.7.1/applications/bin/linux64Gc cDPOpt/SourceSimple"
[7] #6 __libc_start_main in "/home/user/OpenFOAM/user-1.7.1/applications/bin/linux64Gc cDPOpt/SourceSimple"
[6] #6 __libc_start_main in "/home/user/OpenFOAM/user-1.7.1/applications/bin/linux64Gc cDPOpt/SourceSimple"
[0] #6 __libc_start_main in "/lib64/libc.so.6"
[1] #7 in "/home/user/OpenFOAM/user-1.7.1/applications/bin/linux64GccDPOpt/SourceSimple"
[4] #6 __libc_start_main in "/home/user/OpenFOAM/user-1.7.1/applications/bin/linux64Gc cDPOpt/SourceSimple"
[2] #6 __libc_start_main in "/lib64/libc.so.6"
[7] #7 Foam::regIOobject::writeObject(Foam::IOstream::str eamFormat, Foam::IOstream::versionNumbe r, Foam::IOstream::compressionType) const in "/lib64/libc.so.6"
[5] #7 in "/lib64/libc.so.6"
[0] #7 in "/lib64/libc.so.6"
[6] #7 Foam::regIOobject::writeObject(Foam::IOstream::str eamFormat, Foam::IOstream::versionNumbe r, Foam::IOstream::compressionType) const in "/lib64/libc.so.6"
[4] #7 in "/lib64/libc.so.6"
[2] #7 in "/home/user/OpenFOAM/user-1.7.1/applications/bin/linux64GccDPOpt/SourceSimple"
[node014:22802] *** Process received signal ***
[node014:22802] Signal: Floating point exception (8)
[node014:22802] Signal code: (-6)
[node014:22802] Failing at address: 0x20d00005912
[node014:22802] [ 0] /lib64/libc.so.6 [0x3fd6230280]
[node014:22802] [ 1] /lib64/libc.so.6(gsignal+0x35) [0x3fd6230215]
[node014:22802] [ 2] /lib64/libc.so.6 [0x3fd6230280]
[node014:22802] [ 3] /cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libOpenFOAM.so (_ZNK4Foam5PBiCG5solveERNS_5FieldIdEERKS2_h+0x111d ) [0x2aaaac4fc70d]
[node014:22802] [ 4] /cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/lib/linux64GccDPOpt/libfiniteVolum e.so(_ZN4Foam8fvMatrixIdE5solveERKNS_10dictionaryE +0x14b) [0x2aaaab9b313b]
[node014:22802] [ 5] SourceSimple [0x416d64]
[node014:22802] [ 6] /lib64/libc.so.6(__libc_start_main+0xf4) [0x3fd621d974]
[node014:22802] [ 7] SourceSimple(_ZNK4Foam11regIOobject11writeObjectEN S_8IOstream12streamFo rmatENS1_13versionNumberENS1_15compressionTypeE+0x d1) [0x4141f9]
[node014:22802] *** End of error message ***
Foam::regIOobject::writeObject(Foam::IOstream::str eamFormat, Foam::IOstream::versionNumber, Foam: :IOstream::compressionType) const---------------------------------------------------------------- ----------
mpirun noticed that process rank 3 with PID 22802 on node node014 exited on signal 8 (Floating po int exception).
--------------------------------------------------------------------------

Hi there Alex,

By looking at your solver output, I would say the problem lies somewhere in your boundary conditions, time step continuity errors are really large. Try first to use simpleFoam without modifications (ie without the scalar transport equation), confirm every boundary condition, and see if it converges correctly.

Regards,
Jose

try to get rid of the relative solver tolerance in p.
since the mass conservation isn´t converging, your b.c. regarding u and p are wrong. how do you initialize u and p?

Hi
 

Hi there!

I tried to solve the problem using simpleFoam with an additional transport scalar equation without source term and there was a convergency. I setted the relaxation factor all to 0.5.

after I added the source term, the problem started when I increased the variable "a" which is a constant and I can give to it the value I want. My intention was to run some simulation increasing this value and see what is happening, but the results are the ones you saw

Boundary conditions are the following (keep in mind that this is a dump combustor, in which flow first passes through a small tube and after that there is a bigger chamber in which velocity decrease) and in the end there is the outlet.

below you will see the boundary conditions and also a slide of the RANS simulation without addition of the scalar.

-----------------------------------------------------------

object p;

dimensions [0 2 -2 0 0 0 0];

internalField uniform 0;

boundaryField
wall
{
type zeroGradient;
}
inlet
{
type zeroGradient;
}
wall-swirl-plate
{
type zeroGradient;
}
wall-outflow-dump
{
type zeroGradient;
}
wall-chamber
{
type zeroGradient;
}
wall-dump
{
type zeroGradient;
}
wall-inflow
{
type zeroGradient;
}
wall-outflow
{
type zeroGradient;
}
outlet
{
type fixedValue;
value uniform 0;
}

------------------------------------

object U;

dimensions [0 1 -1 0 0 0 0];

internalField uniform (0 0 133);

boundaryField
{

wall
{
type fixedValue;
value uniform (0 0 0);
}
inlet
{

type cylindricalInletVelocity;
axialVelocity 133;
centre (0 0 0); //by default (0 0 0)
axis (0 0 1);
rpm 0;
radialVelocity 0;
tangVelo 77;
value uniform (0 77 133);
}
wall-swirl-plate
{
type fixedValue;
value uniform (0 0 0);
}
wall-outflow-dump
{
type fixedValue;
value uniform (0 0 0);
}
wall-chamber
{
type fixedValue;
value uniform (0 0 0);
}
wall-dump
{
type fixedValue;
value uniform (0 0 0);
}
wall-inflow
{
type fixedValue;
value uniform (0 0 0);
}
wall-outflow
{
type fixedValue;
value uniform (0 0 0);
}
outlet
{
type inletOutlet;
inletValue uniform (0 0 0);
value uniform (0 0 0);

Alex,

The results you shown indicate that the problem lies firstly in the p and U solution. It seems that C is a passive scalar, so the source term in C equation shouldnt affect p and U.

Regarding the value of 'a' in your source term, I think that it can affect convergence of C equation mostly because your source term is defined explicitly. Turn it (or at least part of it) implicit for better stability.

I just saw your figure, and it seems that the length of your computational domain is not sufficient to capture the flow characteristics of your system. Moreover, inletOutlet kills of your flow recirculation at the outlet, so you get a less realistic picture of your system.

Regards,
Jose

Hi Jose,

so you mean, instead of

ST = ((1-C)*(C))*a*runTime.deltaT().value()/runTime.deltaT()

with the use of the "Sp" you showed me last time?

regards,

alex

you could also try to ramp the inlet value "1" with the funkybc. starting with 0.

You can try that option, also reduce the under-relaxation coefficient for C (by the way, is your C equation implemented with under-relaxation in mind?).

Regards,
Jose

You mean the under relaxation factor that you modify on the system/fvSolution ?

All those values are setted as

p 0.5;
U 0.5;
k 0.5;
epsilon 0.5;
C 0.5;
R 0.5;
nuTilda 0.5;

Hi Jose,

I just added the Source Term with the other way you suggest me:

Ceqn.H

solve
(
fvm::ddt(C)
+fvm::div(phi, C)
-fvm::laplacian(nu, C)==
fvm::Sp(mySp, C)
);

and in CreateFields.H I added the following lines:


dimensionedScalar mySp
(
mySp=(C*(C-1))*a/runTime.deltaT().value()
);

but I cannot compile it properly.
do you think I should declare in a different way?

regards,

Alessio

Hi Alex,

As I said before, Sp(MySp,C) stands for MySp*C defined implicitly.

This translates into:
MySp*C = C*(C-1)*a*C, which is not what you want.

Try MySp = (C-1)*a/runTime.deltaT().value().

Regards,
Jose

Hi Jose,

I tried the following :

Ceqn.H
solve
(
fvm::ddt(C)
+fvm::div(phi, C)
-fvm::laplacian(nu, C)==
fvm::Sp(mySp, C)
);

adding in the createFields.H

dimensionedScalar mySp
(
mySp=(C-1)*a/runTime.deltaT().value().
);

and I receive the error:

.........................................
In file included from SourceSimple.C:43:
createFields.H: In function âint main(int, char**)â:
createFields.H:61: error: expected unqualified-id before â)â token
/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/src/finiteVolume/lnInclude/readSIMPLEControls.H:6: warning: unused variable âmomentumPredictorâ
/cvos/shared/apps/OpenFOAM/OpenFOAM-1.7.1/src/finiteVolume/lnInclude/readSIMPLEControls.H:9: warning: unused variable âtransonicâ
make: *** [Make/linux64GccDPOpt/SourceSimple.o] Error 1


regards,

alex

delete the "." in the end. this one is new comparing to your prev. posts.

createFields.H:60: error: no match for âoperator=â in âmySp = Foam::operator/(const Foam::tmp<Foam::GeometricField<Type , PatchField, GeoMesh> >&, const Foam::scalar&) [with Type = double, PatchField = Foam::fvPatchField, GeoMesh = Foam::v olMesh](((const Foam::scalar&)((const Foam::scalar*)Foam::TimeState::deltaT() const().Foam::dimensioned<Type>::value [w ith Type = double]())))â

jap according to http://foam.sourceforge.net/doc/Doxy...Scalar_8H.html there is no "=" operator .... but what do you need it for?
your scalar is named mySp and within the {} it is declared, so there is no need for another assigning the term within the bracelets to mySp
dimensionedScalar mySp
(
(C-1)*a/runTime.deltaT().value()
);

..or am i off the track??

I still get an error, I uploaded the files if you want to have a look.

regards,

alex

Alex,

You are declaring mySp as a dimensionedScalar, but it is determined from a volScalarField C (mySp=(1-C)*a). Therefore, mySp should also be declared as a volScalarField.

Regards,
Jose

Hi Jose, now it compiles without error,

I hope it is right (physically speaking)!

this is what I added on createFields.H :

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

Info << "Reading field mySp\n" << endl;
volScalarField mySp
(
(C-1)*a/runTime.deltaT().value()

);


many thanks !

regards, alex

Thats good news!

In terms of the physics, now its up to you ;-)

Regards,
Jose