[henry]

You can do it that way or simply rewrite the model in terms of btilde.

[grtabor]

OK, good to know that.

I am experiencing problems with bounding the progress variable though. I'm trying to compute a backward-facing step; I've done a precomputation using turbFoam and have converted the data to run with my combustion code. My initial conditions for the ctilde variable are uniform 1.0 internal field, fixed value 0.0 at the inlet, zero gradient at the outlet and the walls (btilde is the inverse of this). However at the end of the first timestep the internal field ctilde has reached a value of 1.04 (pretty much throughout). This is clearly wrong, so...

1. Do my boundary conditions sound correct?

If so, I assume that either my source term is incorrect in the ctilde equation, or I need to do something fancy in terms of bounding the ctilde equation. Are those the only options, or am I likely to be doing something else wrong that cound have this effect?

(This is a rather difficult one to answer - I'm really just trying to pin down what might be going wrong with this solution).

Gavin

[henry]

I would guess your problem is to do with omegaDot. Have you tried to run without it? Did it remain bounded then? Given that you have implemented omegaDot as an explicit source term how do you propose to ensure ctilda remains bounded?

I solved this problem of Sigma-based models by reformulating and remodelling in terms of the flame wrinking Xi which leads to a source term which can be implemnted as an implicit transport term and hence ensures boundedness by suitable choice of scheme. All this is described in detail in the technical reports and papers I have written on the subject.

[grtabor]

Just tried running without omegaDot - it does indeed remain bounded, although there do seem to be problems beyond that. I will investigate this in more detail.

Gavin

[grtabor]

Further (slow) progress...

I've modified the omegaDot (source) term so that the first timestep remains bounded - ctilde now lies between 0 and 1 (chiefly around 1) and btilde between 1 and 0 (to within rounding error - see below). I've also included it as SuSp source rather than as purely explicit.

However, now on the second timestep I get the following error : this is outside the PISO loop.

Time = 0.50002
BICCG: Solving for Ux, Initial residual = 0.00778922, Final residual = 3.99268e-07, No Iterations 1
BICCG: Solving for Uy, Initial residual = 0.0476625, Final residual = 1.67977e-06, No Iterations 1
BICCG: Solving for ctilde, Initial residual = 0.215976, Final residual = 1.8554e-09, No Iterations 2
Got this far!
Got this far2!


--> FOAM FATAL ERROR : Maximum number of iterations exceeded

Function: specieThermo<thermo>::T(scalar f, scalar T0, scalar (specieThermo<thermo>::*F)(const scalar) const, scalar (specieThermo<thermo>::*dFdT)(const scalar) const) const
in file: /home/dm2/henry/OpenFOAM/OpenFOAM-1.0/src/thermophysicalModels/specie/lnInclude/ specieThermoI.H at line: 83.

FOAM aborting

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

AFAICT this traces through to a function in specieThermo which evaluates the temperature by inverting the equation of state using Newton-Raphson methods. Am I right?

It crops up because I call thermo->correct() directly after solving for ctilde and btilde. I assume it has something to do with the state of the btilde field.

btilde actually generates very small negative numbers (-O(1e-10)), presumably through rounding error. Is this likely to be what is causing the problems here? If so, what can I do about it? If not, what is the problem likely to be?

Gavin

[henry]

It looks like you are violating thermodynamic constraints. Are you saying ctilda goes from 0 to 1 in 1 timestep? With ctilde being chiefly around 1 after the first timestep you now have a fully-burnt system; that's an explosion not combustion and will blow the thermodynamics to pieces. Also I don't see the solution of the enthalpy equation.

[niklas]

Hi Gavin,

Yes, I would bet on that.
When we implemented the complex chemistry stuff I had similar problems when species fractions went negative ever so slightly.

Did you try the sledgehammer approach
btilde.max(0);
after the transport eq.

N

[grtabor]

Henry :

>It looks like you are violating thermodynamic >constraints. Are you saying ctilda goes from 0 to 1 >in 1 timestep? With ctilde being chiefly around 1 >after the first timestep you now have a fully-burnt >system; that's an explosion not combustion and will >blow the thermodynamics to pieces.

No, thats not what I mean. My initial condition for the ctilde field is uniform 1.0; i.e. the entire domain is filled with combusted gas. At t=0 my inlet starts introducing fresh gas (ie. ctilde=0.0). After one timestep this has not progressed very far, so the majority of the domain is still at ctilde=1.0, with just the region around the inlet dropping to ctilde=0.0. This seemed to be the most logical way of starting things up - please tell me if I am being stupid.

(I'm doing combustion behind a backward facing step, not (hopefully) detonation or deflagration or anything really galloping).

>Also I don't see >the solution of the enthalpy equation.

I'm not solving it. My impresion from Poinsot and Veynante was that for relatively 'gentle' combustion it was not necessary - this being equivalent to an assumption of adiabatic conditions. The temperature can be derived by working back from the progress variable ctilde, which is in effect a reduced temperature. (I assume I do need to use the thermo database still anyway).

Niklas : thanks, I'll try your suggestion in the first place.

Thanks

Gavin

[grtabor]

Niklas : just tried doing btilde.max(0). Now btilde is always 0 or greater - but I still get the same error on the second timestep.

Ho hum.

Gavin

[henry]

I don't think your ignition approach will work because you are creating an infinite gradient of ctilde at the inlet which corresponds to an infinte reaction rate over an infinitesimal region. If you want to pursue this approach I recommend you ramp the inlet value of ctilte from 1 down to 0 over a period of time which introduces a flame of a sensible/physical/numerically resolved width.

[grtabor]

That sounds feasible/sensible.

Your use of language suggests that there might be another way of doing this - is there? The only other option that I can think of would be to initialise the ctilde field with the value 0 - but then how does the combustion get started?

Gavin

[henry]

Yes the alternative is to initialise a flame within the domain, i.e. a profile going from 0 at the inlet to 1 somewhere downstream. Look at the ignition models proposed for the CFM model for SI engines, there are methods which involve the specification of Sigma and species distributions corresponding to reasonable (~5mm) sized kernel. However, these methods can also produce large initial perturbations due to the initial specification not being exactly consistent with the equations so I think the ramping method will be best/easiest if you choose an appropriate time-scale.

You could improve on the simple ramping by choosing some kind of sigmoid shape for the ctilde distribution being transported into the domain at the inlet. My model has an analytical solution for the 1D flame, you could try using that profile, actually as I remember it was you who first solved my model analytically after I cast it into an appropriate form :-)

[grtabor]

OK - a _really_ simple question this time.

I need to get the current time from the runTime database. How do I do that? I've found a function timeOutputValue() which claims to give the current time value, but it seems rather an odd name for what I want, and also it is returning a scalar rather than a dimensionedScalar (not really a problem, but I'm going to use startTime() which returns a dimensionedScalar, and I'd like to be consistent). Didn't there used to be a currTime() function which gave the current time?

Gavin

[henry]

runTime is a Time which is derived from TimeState which is derived from a dimensionedScalar which is the current time. If you want the time value without dimensions simply use runTime.value() otherwise even more simply use runTime.

[niklas]

Just a thought...

What is your initial temperature?
Seing how T0 and c=1 corresponds to a certain enthalpy H0, I just wonder what the temperature will
be for H0 and c=0.

N

[grtabor]

One further query.... how does one read a label (int) from a dictionary? I do seem to remember that there was a need for something like readLabel(), but I can't find the details.

Thanks

Gavin

[hjasak]

refCelli = readLabel(dict.lookup(refCellName));

[henry]

Yes there is both readLabel and readInt, for further details simply grep for them in say the OpenFOAM-1.2/src/OpenFOAM/lnInclude directory.

[mattos]

Hi Henry

Repeating my message for Tommaso, after long time working in many stupid matters I return to apply my attention and energy to use this wonderfull code: The FOAM!

Please, I need your help. I would like in a first step to use the combustion solvers available in the OpenFOAM pack. But I need the literature which they are based. Could your help me to say which are the references which each solvers are based? I hope also to construct another solver for combustion flows.

Many tanks in adavce for your help.

[rbw]

I would very much like to have these references also.

Having a repository of these papers would be extremely useful, it could be hosted on the OpenFOAM wiki. I think, however, that many if not all of the papers cannot be posted on a web site for copyright reasons. At the very least, a simple list of the more important references would already be very useful.

A lot of the work done on OpenFOAM seems to have been done during PhD thesis. Having all of them on the same spot would also be useful. I have Hrvoje's and Niklas' (both are stupendous works !)

It would be nice if someone would collect the publicly available references and put them up on the wiki. It would be a simple and useful non-programming way of contributing to OpenFOAM.
If no one else takes up the hint, I may volunteer myself in a couple of weeks time