Changing reaction rate model

[Gary Sharpe]

Carlos, I am starting a new thread for this as the old
one is getting a bit longer and off the original topic.

Re: your email:
>Yeah, I have the changed the ReactiveEulerEquations.amr. I haven't made
>any modifications since you last looked at it, but it may be a good idea
>to cut the induction zone Q term out of the equation set now

Yes, as I first stab I would keep it simple. You will have to modify
the
specific ReactiveEulerEquations.amr I sent in a similar way you did
before.

>I think a good place to start would be to modify the ComputeZndProfile
>procedure.

Actually, this is probably the last piece of the jigsaw to consider -
all ComputeZNDprofile
is really needed for is to set up special initial conditions - it is
not worth
getting this right if we can't get the BCG/ReactiveEulerEquations
working
in the first place - other simple initial conditions can be used to
test this such as a top-hat
pressure profile (the standard 'direct ignition' startup).
I suggest starting by looking at the different files in the cloned bcg
directory
folder I sent and look at what would need to be changed in going from
the
3step to the 2step model.

[Carlos]

Gary,

I have modified the ReactiveEulerEquations.amr file and dealt away
with the Qi term, and now I simply have a Q term associated with the
second mass fraction.

My first question involves the SourceTerms.amr file. In this file I
modified P, T, DfDt and DyDt. My DfDt and DyDt read something like:
DfDt = -Ki*EXP(-E/T)
DyDt = -Ke*(1-y)**(ETA)

For our model of interest, would it be necessary to include a
switch? , for example:
If WW(5) < 1, then DfDt = as above
else DfDt = 0
If WW(5) < 1, then DyDt = 0
else DyDt = as above

My second question is about the Startup.amr file. I've noticed a
parameter, NU that's listed twice:
WRITE(AMR::LOGFILE(),*) 'NU = ',NU
WRITE(AMR::LOGFILE(),*) 'NU = ',NU
Is there any need for this to be listed twice or is it an extra line
that can be deleted?

Thanks,
Carlos

[Eric Marineau]

Thanks, I will look at the data next week
-Eric
---
Eric C. Marineau, Ph.D.
Postdoctoral Scholar in Aeronautics
T5 Hypervelocity Shock Tunnel
California Institute of Technology
MC 105-50
1200 East California Blvd
Pasadena CA 91125
(626) 395-8311 Fax: (626) 449-2677


-----Original Message-----
From: Carlos <cleu...@uottawa.ca>

Date: Wed, 2 Jul 2008 13:06:01
To: amrita-ebook<amrita...@googlegroups.com>
Subject: Re: Changing reaction rate model

[Gary Sharpe]

> My first question involves the SourceTerms.amr file. In this file I
> modified P, T, DfDt and DyDt. My DfDt and DyDt read something like:
> DfDt = -Ki*EXP(-E/T)
> DyDt = -Ke*(1-y)**(ETA)
>
> For our model of interest, would it be necessary to include a
> switch? , for example:
> If WW(5) < 1, then DfDt = as above
> else DfDt = 0
> If WW(5) < 1, then DyDt = 0
> else DyDt = as above

Looking at it, it should be something like

IF (f.GE.1.D0) THEN
DfDt=0.0
DyDt=...as above
ELSE
DfDt=...as above
DyDt=0.0
ENDIF

Note it appears to be doing a simple
two step integration so you will
have to do the same again, but
this time based on the half time
step value fs.

Note also keeping the activation
pressure Pact in is important to
prevent reaction in the initial quiescent
state.
Also, take care to get the expression
for T correct.

> My second question is about the Startup.amr file. I've noticed a
> parameter, NU that's listed twice:
> WRITE(AMR::LOGFILE(),*) 'NU = ',NU
> WRITE(AMR::LOGFILE(),*) 'NU = ',NU
> Is there any need for this to be listed twice or is it an extra line
> that can be deleted?

Yes, this looks like an duplication.
However, it isn't too important as
this just writes to a logfile.

Gary.

[Carlos]

That is for both progress variables going from 0 to 1, right? If I set
it up in a similar way as the 3cb model (f goes from 1 to 0, y goes
from 0 to 1), it should be something like

IF (f.LE.0.D0) THEN
DfDt = 0
DyDt = Ke*(1-y)**(ETA)
ELSE
DfDt = -Ki*EXP(-E/T)
DyDt = 0
END IF

ie. LE 0 instead of GE 1

[Carlos]

I've gone back and ran the old pulsate file using the 3cb model, but
the simulation seems to get stuck when it tries to march in time, even
after changing the node spacing and letting it run for a whole day. It
seems to be able to build the znd code and the solver code fine. I
simply took one part of the code and modified it to look like:

plugin amr_sol
utilize 2step::Machinery
ReactiveEulerEquations {
space = 1D
model = 3cb
}

I was wondering whether I didn't implement the model properly, or
whether the 3cb model is really that much more expensive than the 1a
model

[James Quirk]

Carlos,

On Fri, 4 Jul 2008, Carlos wrote:

>
> I've gone back and ran the old pulsate file using the 3cb model, but
> the simulation seems to get stuck when it tries to march in time, even
> after changing the node spacing and letting it run for a whole day. It
> seems to be able to build the znd code and the solver code fine. I
> simply took one part of the code and modified it to look like:
>
> plugin amr_sol
> utilize 2step::Machinery
> ReactiveEulerEquations {
> space = 1D
> model = 3cb
> }
>
> I was wondering whether I didn't implement the model properly, or
> whether the 3cb model is really that much more expensive than the 1a
> model
>

There is clearly a problem with your implementation.
Post a mailit and I'll take a look at it.

James

[cleu...@uottawa.ca]

When I try to make the mailit, I get the same problem. The process freezes
and the mailit can't get created. I've attached my script.

In the directory where the script is contained, I have the sub-directory
/3step/Machinery, where the cloned solver/znd/equation set files are
stored, as well as a AmritaUtilize.amr

Thanks,
Carlos

[James Quirk]

Carlos,

On Sun, 6 Jul 2008, cleu...@uottawa.ca wrote:

> When I try to make the mailit, I get the same problem. The process freezes
> and the mailit can't get created. I've attached my script.

You could have made a mailit of your directory tree
without actually running the script. That way I'd have
a complete record of what you're doing.

>
> In the directory where the script is contained, I have the sub-directory
> /3step/Machinery, where the cloned solver/znd/equation set files are
> stored, as well as a AmritaUtilize.amr

AMRITA is a system that helps those who help themselves.
Here you've sent me your driver script but the details
I need to look at are in 3step/Machinery directory
and the AmritaUtilize.amr procedure, both of which
are missing.

Have another go. But this time, copy all the relevant
file to a directory xyz, preserving whatever directory
structure you're using then run:

amrita -mailit xyz

You could simply run amrita -mailit on the directory
tree where you're working, but that would pack binary
files as well as text files and so, depending on
what you've done, may result in a file that is
too big to mail.

James

[Gary Sharpe]

Carlos,
I am a bit surprised you are ready to run the script already.
Have you modified all the files already then? I was hoping you were
going to
post each one in turn here, so we could discuss (and if necessary,
iterate) them - so that there would be a record for the next person.

Have you already re-done the ComputeZNDwave file too?
If so have you plotted out what your initiation actually gives
using the along command for P, RHO, U, Zi, Ze
- i.e. after setting up the ZND wave on the grid but
BEFORE trying to march - this may illuminate whether
something is wrong?

Gary

[cleu...@uottawa.ca]

James,
I've attached the mailit for my script
Thanks in advance

Gary,
I haven't finished modifying the files yet. I was just testing the 3 step
chain branching files that you sent me so that I know how to implement it
into the 1-D galloping detonation script. I'll post what I've done and the
list of modifications I've made here

ReactiveEulerEquations.amr - I put in all of the parameters that the model
uses here. I also got rid of the Q in the induction zone that I had when
you last saw it

Startup.amr/model.h - Modified so that it grabs the problem specific
parameters for the 2 step model instead of the 3 step model

PartialDerivatives.amr - I basically modified it so that the pressure reads:
P = (GMz-1)*(D*i-D*Q*y)
and did the appropriate partial derivatives based on that

SourceTerms.amr - I modified the P and T equations and put in the
equations DfDt and DyDt as described in previous posts

I believe the ComputeZndProfile.amr is the only file that needs to be
modified that I haven't yet touched.

[Gary Sharpe]

> Gary,
> I haven't finished modifying the files yet. I was just testing the 3 step
> chain branching files that you sent me so that I know how to implement it
> into the 1-D galloping detonation script.

Do you mean you are just trying to run
with the original, unmodified files I sent?
If so I have uploaded a driver script for the
3step pulsating case. See if that works.

[cleu...@uottawa.ca]

Gary,

I've had a look at the ComputeZnd file

Do you think we can get rid of the chunk of code that re-scales the value
for RS? I believe we can scale Ki without all those calculations

Is the QR term simply the local chemical heat energy term? So for us, for
our second progress variable YR, going from 0 to 1, it should read
something like
QR = - YR * Q ?

[Gary Sharpe]

Have you sorted out the other problems you were having then.
Did the 3step files work with the driver script I uploaded?

Yes, for the 2step model, Ki can be set algebrically to get the
induction lenght
scale of the ZND to be unity, so no need to integrate. Ke is then
a parameter the user will supply anyway.
For the Q term, won't you need a switch on Zi again, as the heat
release only occurs once induction step is over?
Hmmm, though YR will be 0 up to that point anyway won't it
so you are probably correct.
Hence what the code should do is write a constant state
between the shock and x=-1, corresponding to the induction zone,
and then, starting from x=-1 integrate the equations you have
for the main heat release stage.
Note that the 3step model is scaled such that P_shock=rho_shock=1
- is this the way you want it scaled for the 2step case?

Gary.

[cleu...@uottawa.ca]

The 3 step script you uploaded worked. Thanks

I'm happy with how the problem is scaled currently. I believe the speed,
U, right now is scaled with the post-shock sound speed

I've modified the ComputeZndProfile script and uploaded what I've changed
so far. I've deleted what was previously lines 144-199, with the scaling
of Ki. I've deleted RI = ..., RB = ..., RC = ... and I've changed QR as
discussed above.

I've also changed FX1, YX1, QX1 to:
IF (FR.LE.0.D0) THEN
FX1 = 0
YX1 = (Ke*(1-YR)**ETA)/UR
ELSE
FX1 = -(Ki*EXP(-E/TR))/UR
YX1 = 0
END IF
QX1 = -Q*YX1

and likewise with FX2, YX2 and QX2.

However, when I try to run the script, it writes the data for one point,
but stops the computation there.

I suspect it has something to do with these lines:

IF((FR.GE.Ftol).OR.(YR.GE.Ytol)) GO TO 300
IF(X.LT.Xc) GO TO 300

and that the program doesn't want to loop, but looking at the data
inputted, it seems that the program should loop (FR > FTol so it should
loop).

[cleu...@uottawa.ca]

I've had a stab at running the znd script (attached below). It has several
issues when I try to run it.

When the mass fraction of the fuel in the induction zone approaches 0, all
of the profiles start to blow up (P goes from 1 to -315000). When I make a
some computations by hand, I don't run into the same problem, so this may
be a matter of weak programming on my part

For the fuel mass fraction, I have something like:
FR = MAX(0.0 D0,FRL+(FX1+FX2)*DX/2.)

However when I look at the znd profile produced, the mass fraction dips to
negative numbers. Any idea what may be causing this?

Carlos

[Gary Sharpe]

> I'm happy with how the problem is scaled currently. I believe the speed,
> U, right now is scaled with the post-shock sound speed

Yes.

> I've modified the ComputeZndProfile script and uploaded what I've changed
> so far. I've deleted what was previously lines 144-199, with the scaling
> of Ki. I've deleted RI = ..., RB = ..., RC = ... and I've changed QR as
> discussed above.
>
> I've also changed FX1, YX1, QX1 to:
> IF (FR.LE.0.D0) THEN
> FX1 = 0
> YX1 = (Ke*(1-YR)**ETA)/UR
> ELSE
> FX1 = -(Ki*EXP(-E/TR))/UR
> YX1 = 0
> END IF
> QX1 = -Q*YX1
>
> and likewise with FX2, YX2 and QX2.

As I hinted at in previous post, for the ZND model, you don't need to
do this.
For ZND you have initial state (ahead of shock), shock state, constant
induction zone (apart from F which varies linearly - I assume you are
scaling length so this induction zone is length 1?), main reaction
stage
(F=0). So you can calculate and write out these stages separately.

Initial state point at x=0
shock state point at x=0 (P=rho=F=1, Y=0)
end of induction zone point at x=-1 (P=rho=1, F=Y=0)
then several points in
main reaction stage (all varying apart from F=0).

> However, when I try to run the script, it writes the data for one point,
> but stops the computation there.
>
> I suspect it has something to do with these lines:
>
> IF((FR.GE.Ftol).OR.(YR.GE.Ytol)) GO TO 300
> IF(X.LT.Xc) GO TO 300
>
> and that the program doesn't want to loop, but looking at the data
> inputted, it seems that the program should loop (FR > FTol so it should
> loop).

end of induction zone occurs at FR=0, so it will never get there
unless
Ftol was negative (this is not like 3step case where formally F only
asymptotes to zero).

[Gary Sharpe]

> I've had a stab at running the znd script (attached below). It has several
> issues when I try to run it.

See previous post. In ZND wave FR is simple to deal with. It is unity
ahead of shock,
varies linearly from 1 to 0 in induction zone, and is zero in main
reaction zone.
Hence there is no need (and you shouldn't be) integrating it
numerically.
Also, for ETA=1/2, the main reaction zone I believe has an analytical
solution
for Y (and hence P, rho, U) as a function of x from the conservation
equations.
It might be an idea as a first test to use this analytical solution,
i.e. just
get the ComputeZND code to write out the analytical structure as
function
of x, without any numerical integration involved.
We can go from there depending on whether or not that works....

Cheers
Gary

> ComputeZndProfile.amr
> 17KDownload

[cleu...@uottawa.ca]

I have fixed up my Znd profile. I now have a working ComputeZndProfile
procedure attached below. It turned out that I simply had to delete the
existing code that bcg created so that the updated and proper code could
be generated. I also had to make sure the second mass fraction did not
exceed 1, as that causes some problems too

I am now running a simulation to see if the model actually runs. Depending
on the results, I may have to go back and modify the files in the bcg
directory

Carlos

[Gary Sharpe]

Carlos,
You can force it to recompile using
amrita -sys redo ....
It looks OK now. You should change the model name from 3cb
to something like 2ie to avoid confusion.
Can you send the znd.RHO, etc, files that the procedure
should produce?
Gary

> ComputeZndProfile.amr
> 17KDownload

[cleu...@uottawa.ca]

I've attached the znd profiles for the parameters d = 1, ETA = 0.5, Ke =
0.2, Q = 4, E = 3.205

[Gary Sharpe]

Carlos, well they certainly look OK and reasonable.
Did you manage to get the script to run with these
initial conditions?

> znd.P
> 135KDownload
>
> znd.RHO
> 183KDownload
>
> znd.T
> 129KDownload
>
> znd.U
> 183KDownload
>
> znd.Ze
> 114KDownload
>
> znd.Zi
> 108KDownload

[cleu...@uottawa.ca]

Gary,

The script runs just fine, but there is definitely something wrong with
the simulation. The pressure does not behave as it ought to. I've attached
a .ps file of the P history. I may have mistyped an equation or forgot to
change something in the solver files. I'll need to peruse through the
modified files in the bcg directory with a closer eye

[cleu...@uottawa.ca]

I can't seem to pinpoint the bug in the code. I've made some fixes to my
SourceTerms.amr file. The equations for DfDt, DyDt, P and T I'm pretty
confident are correct now.

[Gary Sharpe]

Can you send a mailit with everything in,
and I will take a look

[cleu...@uottawa.ca]

I've attached the mailit.

Thanks,
Carlos

[Gary Sharpe]

James, can you assist on this one at this stage...?
Carlos, I'm am not sure what solver is getting built in your script
..in build solver you still have
scheme = roe'3cb.
However, when I change this to roe'2ie, it reports
***********************************************
BCG can't produce a 1d'cartesian'roe'2ie scheme
***********************************************
However, if I look at the source files created in directory
"code" with scheme = roe'3cb,
it does seem to have built the solver with your modified
2ie bits. So I am unclear what is happening here.
James could you advise?

Carlos, some observations:
(i)the ZND setup looks to be working fine.
(ii) you need to change your refinement criteria,
if you look at the profiles, you will see that is refining
everywhere behind the shock not just the reaction zone,
so that your calculations are very expensive
(iii) for initial testing you don't need to run with
the very high resolution .. I would set Npts=4
and lmax=1 initially. I would also set E to be small
or zero such that the ZND wave is definetly stable.
(iv) when it runs, the Zi and Ze profiles remain
looking reasonable, while the flow variables go
'haywire'. This suggests the reaction source terms
are correct, but that the coupling through Q isn't
quite right, but until the solver issue above is resolved
its unclear why....

Gary

> 2step.mailit
> 56KDownload

[James Quirk]

Gary,

On Thu, 24 Jul 2008, Gary Sharpe wrote:

>
> James, can you assist on this one at this stage...?
> Carlos, I'm am not sure what solver is getting built in your script
> ..in build solver you still have
> scheme = roe'3cb.
> However, when I change this to roe'2ie, it reports
> ***********************************************
> BCG can't produce a 1d'cartesian'roe'2ie scheme
> ***********************************************
> However, if I look at the source files created in directory
> "code" with scheme = roe'3cb,
> it does seem to have built the solver with your modified
> 2ie bits. So I am unclear what is happening here.
> James could you advise?

What you're seeing is a limitation of the current implementation
of BCG. As it stands, BCG allows you to overload its innards,
but it does not allow you to add to its innards, at least
not in a convenient fashion. Thus here Carlos is more or less
forced to retain the "3cb" handle when requesting a solver,
even though the code-generating routines now produce a 2ie code.

As it happens, BCG has now risen to the top of my todo list
in terms of revamping AMRITA's sub-systems and I hope to provide
a much cleaner means of users adding reaction models, thus
removing many of the hoops that Carlos is currently jumping through.

James

>
> Carlos, some observations:
> (i)the ZND setup looks to be working fine.
> (ii) you need to change your refinement criteria,
> if you look at the profiles, you will see that is refining
> everywhere behind the shock not just the reaction zone,
> so that your calculations are very expensive
> (iii) for initial testing you don't need to run with
> the very high resolution .. I would set Npts=4
> and lmax=1 initially. I would also set E to be small
> or zero such that the ZND wave is definetly stable.
> (iv) when it runs, the Zi and Ze profiles remain
> looking reasonable, while the flow variables go
> 'haywire'. This suggests the reaction source terms
> are correct, but that the coupling through Q isn't
> quite right, but until the solver issue above is resolved
> its unclear why....

There is no solver issue other than a mismatch in nomenclature.

[Gary Sharpe]

> What you're seeing is a limitation of the current implementation
> of BCG. As it stands, BCG allows you to overload its innards,
> but it does not allow you to add to its innards, at least
> not in a convenient fashion. Thus here Carlos is more or less
> forced to retain the "3cb" handle when requesting a solver,
> even though the code-generating routines now produce a 2ie code.

OK, I guess what I'm not clear about is why it still finds and builds
the 2ie routines even though you have specified 3cb handle .. i.e.
how does it know to look in the right place with the wrong handle?

> As it happens, BCG has now risen to the top of my todo list
> in terms of revamping AMRITA's sub-systems and I hope to provide
> a much cleaner means of users adding reaction models, thus
> removing many of the hoops that Carlos is currently jumping through.

OK thats good, in the meantime Carlos and should perservere so we have
a working model.

[Gary Sharpe]

> There is no solver issue other than a mismatch in nomenclature.

OK, Carlos, we will have to hunt down the problem.
I have modified your script, tidying up, and removing some extraneous
bits that aren't needed at the moment. I have also fixed the
refinenemt
criteria issue, and setup it up for debugging - i.e. it now generates
profiles
after every time step, with a moderate resolution,
so you can see what is happening and going wrong.
I will upload this as 2step_ZND_1D.

Now, I don't think it is an issue with the initial conditions,
I have checked that mass and momentum are conserved throughout
your initial profile - but you may want to do the same for total
energy as a check (see below).
[By the way, you may want to change
DXtarget =1.D0/float(Npts)
to
DXtarget = 0.5D0/float(Npts)
in ComputeZNDProfile - otherwise it interpolates a point
inside the shock instead of starting with the shock as
a discontinuity - and the interpolation is such that
the mas flux, etc is wrong at this point. I am not
sure what this will do, beyond changing the startup
error due to initial shock smearing, but in some sense
it is a 'non-standard' starting condition].

Looking at the evolution, one notes that intially
RHO, U, Ze and Zi all remain pretty much unchanged
for the first few time steps (in shock rest frame) and
the issue seems to mainly be with P. For P (apart
from the startup error causes by shock smearing),
the induction zone remains OK initially, but you can
see that behind it in the reaction zone, P drops
pretty much everywhere straightaway.
I suspect therefore the issue is with your
defintion of internal energy somehow. I would
check that the EOS dependent terms, and partial
derivatives in the solver which gets produced
(code/2step.src) are correct and what you intended.
If I get time to day (doubtfull) I will go through the maths
for this case, and check if I agree with what you
currently have.
Cheers
Gary

[James Quirk]

Gary,

>
> > What you're seeing is a limitation of the current implementation
> > of BCG. As it stands, BCG allows you to overload its innards,
> > but it does not allow you to add to its innards, at least
> > not in a convenient fashion. Thus here Carlos is more or less
> > forced to retain the "3cb" handle when requesting a solver,
> > even though the code-generating routines now produce a 2ie code.
>
> OK, I guess what I'm not clear about is why it still finds and builds
> the 2ie routines even though you have specified 3cb handle .. i.e.
> how does it know to look in the right place with the wrong handle?

BCG is just a collection of code-generating routines which happen
to be arranged in a tree. When you run:

plugin amr_sol
ReactiveEulerEquations {
model = 3cb
}
BasicCodeGenerator {
scheme = roe'3cb
code = abc
}

BCG switches to amr_sol's code-generating tree then jumps
to the sub-tree used for the Reactive Euler equations. It then
uses the model identifier, or handle, to jump to the sub-sub-tree
of generating routines. At that point it looks for a file,
bcg.schemes, which lists the available schemes. In this instance
the file looks like:

def NamingConvention
convention for reactive-euler-code is space'grid'solver'model
where space = {1d:one-dimensional|2d:two-dimensional}
where grid = {c:cartesian|b:body-fitted:curvilinear}
where solver = {roe|contrib}
where model = {1a:1step-Arrhenius:OneStepArrhenius:Arrhenius}
where model .= {sd:sqrt-Depletion:SquareRootDepletion}
where model .= {3cb:3step-ChainBranching:ThreeStepChainBranching:ChainBranching}
where model .= {axs}
exclude names 1d'b'*'*
exclude names *'*'roe'axs
end def

which explains why you can't run:

utilize 2steplib
plugin amr_sol
ReactiveEulerEquations {
model = 2ie
}
BasicCodeGenerator {
scheme = roe'2ie
code = abc
}

The utilize command does rejig AMRITA's search path such
that the call to ReactiveEulerEquations invokes Carlos's
version rather that the one in:

$AMRITA/stdlib/equations

However it does not inform BCG that it has a private
code-generating tree and so BCG still sees the original
bcg.schemes file and so fails to recognize 2ie as a
known model.

Now when you run:

utilize 2steplib
plugin amr_sol
ReactiveEulerEquations {
model = 2ie
}
BasicCodeGenerator {
scheme = roe'3cb
code = abc
}

BCG sees the known model '3cb' and starts generating a 3cb code.
But since the utilize command includes an:

autopath <stdlib> ++$amrita::thisdir/bcg/2ie

when BCG invokes a routine such as PartialDerivatives,
it is tricked into finding Carlos's version for 2ie rather
than the original 3cb routine. Hence your confusion.
The code produced is the one you want, but the nomenclature
is messed up.

>
> > As it happens, BCG has now risen to the top of my todo list
> > in terms of revamping AMRITA's sub-systems and I hope to provide
> > a much cleaner means of users adding reaction models, thus
> > removing many of the hoops that Carlos is currently jumping through.
>
> OK thats good, in the meantime Carlos and should perservere so we have
> a working model.

Yes.

BCG was cobbled together back in 1995 and has provided
a reasonable service, for a wide range of people, ever since.
But for some time it has been clear that the implementation
needed to be revamped to keep pace with peoples' needs.
Therefore if you and Carlos keep plugging away I will
use whatever you produce as a test in ensuring that
the revamped BCG does indeed provide a more convenient
means of being overloaded with different reaction models
and the like.

James
>
>
> >

[Gary Sharpe]

Carlos, I believe you have the sign in front of the Q wrong in your
pressure routines in the solver.
You have
P=(G-1)*D*(i-Q*Z2)
giving
i=P/(G-1)/D+Q*Z2
for the internal energy
shoudln't it be
i=P/(G-1)/D-Q*Z2?
i.e. internal energy gets converted into heat
and hence i decreases as QZ2 increases?
(see CTM paper by Short and Sharpe)

When I change the signs appropriately in all the bcg EOS routines
(pressure and partial derivatives), and run, it does at least
tend to a suitable looking steady state ZND profile. However, this
ZND profile is quite different to the original one, so I suspect
there is a similar error in the ComputeZNDProfile procedure
...in particular the (CJ) state just behind the long time
reaction zone structure is not that predicted by your ZND
solution.

[Carlos]

> Carlos, I believe you have the sign in front of the Q wrong in your
> pressure routines in the solver.
> You have
> P=(G-1)*D*(i-Q*Z2)
> giving
> i=P/(G-1)/D+Q*Z2
> for the internal energy
> shoudln't it be
> i=P/(G-1)/D-Q*Z2?
> i.e. internal energy gets converted into heat
> and hence i decreases as QZ2 increases?
> (see CTM paper by Short and Sharpe)

Ah yes, that's right. I made a mistake there. I believe I made the
same mistake in the equation set

> ...in particular the (CJ) state just behind the long time
> reaction zone structure is not that predicted by your ZND
> solution.

What do you mean by this? When I ran it, the state behind the reaction
zone was the same as the state that the ComputeZnd procedure produced
(P = 0.57005644 for both cases)

Carlos

[Gary Sharpe]

> Ah yes, that's right. I made a mistake there. I believe I made the
> same mistake in the equation set

> What do you mean by this? When I ran it, the state behind the reaction
> zone was the same as the state that the ComputeZnd procedure produced
> (P = 0.57005644 for both cases)

Carlos, OK I missed the equation set mistake. If you make the sign
correction
throughout PartialDerivatives, ReactiveEulerEquations AND SourceTerms
then your code more less works now, in that it propagates your
original
ZND profile.

Looking closely at the late time solution, there is a slight
difference
between the numerical structure and the initial ZND profile. Looking
closely, Ze actually starts to rise slightly above 0 before
Zi reaches zero (I attached a plot - ZeZi58.jpg)
It is probably something to do with the way you have done the source
terms
...we should look at this in more detail next week...

[Gary Sharpe]

Carlos, also, there is some differnece in the reaction zone
structure from the original ZND profile, near the end of the reaction
zone.
See figure P58 which compares the numerical structure with the ZND
profile... the numerical solution appears to have Ze tending to zero
asympotically as if ETA=1 rather than the slope discontinuity for
ETA=0.5 in the
ZND profile... again something we need to look at next week
Gary

[cleu...@uottawa.ca]

Gary,
I can't seem to see the figures that you uploaded. Maybe there was an
error in uploading, as it shows the file sizes as 0 bytes

> Looking closely at the late time solution, there is a slight
> difference
> between the numerical structure and the initial ZND profile. Looking
> closely, Ze actually starts to rise slightly above 0 before
> Zi reaches zero (I attached a plot - ZeZi58.jpg)
> It is probably something to do with the way you have done the source
> terms

I notice that this happens when Zi is very close to 0. See attached ZiZe.
I think that it may be due to the 2 step integration in the SourceTerms
file. It's possible that the first step of the integration predicts that
Zi would reach 0, so that in the second step of the integration, Ze starts
prematurely, before Zi reaches 0 completely.

I haven't noticed this. In the ZiZe graph that I've attached, Ze actually
hits 1. What were the variables that you used?

[Gary Sharpe]

> I notice that this happens when Zi is very close to 0. See attached ZiZe.
> I think that it may be due to the 2 step integration in the SourceTerms
> file. It's possible that the first step of the integration predicts that
> Zi would reach 0, so that in the second step of the integration, Ze starts
> prematurely, before Zi reaches 0 completely.

I have played around with doing this differently and making sure
the above doesn't happen. I can alleviate the problem slightly, but
not cure it.
I think that numerical diffusion is playing a role here - in the model
you essentially have a discontinuity in the gradient of Ze here (and
in
the flow variables). In the numerics these large changes in the
gradient
would I assume result in large artificial diffusion locally.

> I haven't noticed this. In the ZiZe graph that I've attached, Ze actually
> hits 1. What were the variables that you used?

It was in the pressure and other flow variables. Again for ETA=0.5,
there is a discontinuity in the flow variables at this point, which
the
numerics will find it hard to handle properly. In particular you need
to make sure that this point is refined on the finest grid - I found
if
you set Zetol to something like 0.999 it alleviates the problem
somewhat.

I have done a resolution study to see if resolution improves these
issues.
However, I am still finding noticeable differences between the late
time
solution and the ZND solution - see pressure.jpg I uploaded. In
particular
the induction zone is somewhat longer than unity, and the solution
doesn't appear to be converging to the 'right' answer.

[Gary Sharpe]

> I have done a resolution study to see if resolution improves these
> issues.
> However, I am still finding noticeable differences between the late
> time
> solution and the ZND solution - see pressure.jpg I uploaded. In
> particular
> the induction zone is somewhat longer than unity, and the solution
> doesn't appear to be converging to the 'right' answer.

OK, the error is in the ComputeZNDProfile, not in SourceTerms
- if you look at what it writes out for dimensionless U
for the ZND file, it is actually (Ds-UR)*SQRT(G),
hence UR that appears in your routine is not your
dimensionless fluid speed, but is a factor of sqrt(gamma)
out (UR is essentially the fluid speed scaled with the
post-shock sound speed, while our non-dimensionizlation
for speeds is with sqrt(Pvn/RHOvn) for consistency).
Hence your Ki is out by a factor of SQRT(G) - which
also means SourceTerms is using this incorrect Ki,
but then your reaction rate terms (and hence length scales
which get calculated) in your ZND profile are also out by SQRT(G).
Fixing all these, then the comparison between ZND and
simulation profiles agree well, so I think we are in business.
I will upload a mailit with all my changes and edits in
so that we are working with the same set of files...

[Gary Sharpe]

Carlos,
I have uploaded new 2step.mailit
which contains all the changes - you should
probably use these so we are using consistent files.
I have changed the source terms to either do
an induction step or an heat release step
during one chemical time step - I think
this is the most consistent way to do it.
I have tried other possibilities (SourceTerms2
and your SourceTerms_Orig) but it doesn't make
too much difference.

By careful of Zetol - I have set this to 0.9999
for ETA=0.5 to try and refine whole reaction zone,
but if you use this for ETA=1 it will probably refine
very far back from shock since Y then only tends
asymptotically to 1.

I think the next step is to see if you can recover the
stabilitiy results in Short and Sharpe CTM paper,
produced by a different code, with your new solver.
You will have to be careful about going between the
different non-dimensionilzations used (E in your solver
and "epsilon" in the paper)

Cheers
Gary

[Carlos]

Thanks a lot for all the help. I'll take a look at the new solver

[Carlos]

Gary,
I've played around with the parameters (both the physical parameters
and the refinement parameters) and it seems that no matter how I
fiddle around with it, I can't get anything reasonable. I've uploaded
an image of the shock pressure history. This particular simulation
used E = 10, ETA = 0.5, Q = 1.7, and Ke = 0.2, but it seems that just
about any set of parameters produce this shape of a profile. Do you
think there's still a mistake we haven't caught in the bcg files, or
maybe there's a glitch in the pressure recording routine?

On Jul 29, 10:01 am, Gary Sharpe <men...@leeds.ac.uk> wrote:

[Gary Sharpe]

Carlos, I think we sorted this one out as the issue with
not recompiling code.
Any progress on reproducing previous results?

[cleu...@uottawa.ca]

Ah yes, I posted that before we sorted it out.

I was able to reproduce E=10 almost identically to what was in the paper.
Right now, I'm running tests to try to find the neutral stability boundary
for different activation energies. There are small differences between the
results, but I suspect it may be a rounding issue when converting from one
normalization to the other

Carlos