Reactor Networks: Gas Turbine Premixed Combustor- Division by zero error and Ignition problem

[Vaibhav Prakash]

Hello Cantera Users!

I am a new user to Cantera and I am using Cantera 2.3.0 and Python 3.4 in my Windows 10 System. I am working on my Masters thesis on building an emission model of industrial gas turbine combustors under Exhaust Gas Recirculation. I have looked into the examples of combustor.py and pfr.py examples provided by Cantera. I am trying to build the system with the premixer (upstream of combustor) as an Ideal Gas Reactor and the main combustor as a reactor Network (CSTR) and PFRs.

I am using Natural Gas as the fuel with GRI3.0 Mechanism. The premixed mixture's entry conditions into the combustor is at 700 K and 10.3 bar. As a start, I am trying to connect the premixer (with 2 inlets- Natural Gas and EGR oxidizer) with the main combustor (different volume). I am igniting the mixture as per combustor.py wiht a stream of H radicals. Initially, I have assumed the whole combustor volume as a single Ideal Gas reactor.

So the outline is basically:

'''
Fuel stream pipe----|
                     --------->  PREMIXER (Ideal Gas Reactor) ---------> Combustor (CSTR)
EGR oxidiser pipe---|

''' 

Problems and Errors for which I need help:

1.  After the combustor, an exhaust valve has been connected to calculate the residence time of the combustor.

#exhaust valve

v2=ct.Valve(combustor, exhaust2, K=100)

#reactor Network

sim=ct.ReactorNet([mixer2, combustor])

t_init2 = 0.0
Tprev = combustor.T
states = ct.SolutionArray(fuel2, extra=['t','tres'])

for n in range(100):
    tres2= combustor.mass/v2.mdot(t_init2)
    t_init2+=2*tres2
    sim.advance(t_init2)


print(combustor.thermo.report())

 But, the simulation loop shows  the following error:

Traceback (most recent call last):

  File "E:/Cantera/combustor_network.py", line 89, in <module>

    tres = combustor.mass/v2.mdot(tnow)

ZeroDivisionError: float division by zero

How to correct this error? I have correctly defined the mass flow rate as per the experimental work specifies. I have attached my code for your reference. 

2.  Also, in order to make the simulation run and to avoid the division by zero error, I filled in a value of 1.4 inside the loop to calculate the residence time of the reactor. But I am getting weird results and the Mixture is not getting ignited even if change the initial temperature of the igniter.

I get these results:

gri30:

       temperature             700  K

          pressure        1.03e+06  Pa

           density         4.95761  kg/m^3

  mean mol. weight         28.0135  amu

                          1 kg            1 kmol

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

          enthalpy      4.2636e+05        1.194e+07     J

   internal energy       2.186e+05        6.124e+06     J

           entropy          7050.3        1.975e+05     J/K

    Gibbs function     -4.5089e+06       -1.263e+08     J

 heat capacity c_p          1095.3        3.068e+04     J/K

 heat capacity c_v          798.55        2.237e+04     J/K

                           X                 Y          Chem. Pot. / RT

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

                N2              1                1         -21.7022

     [  +52 minor]    4.23751e-26      1.52475e-27

Process finished with exit code 0

This means that the combustor, initially filled with hydrogen is not ignited at all or is not properly connected by a mass flow controller from the premixer. How to solve this problem?

Kindly help me out with this problem. I have been trying hard to find solutions online, but so far there are no successful attempts. I need to tackle this problem to get going with my thesis! 

Kind Regards,

Vaibhav Prakash

[Nick Curtis]

Hi Vaibhav,

1.  I believe the issue here is that the transient mass flow through V2 is initially zero, hence the residence time is undefined here.

The easiest fix would be to take an initial step before computing the residence time.

If you're looking for the steady state behaviour, it may be simpler to just take a few time-steps and monitor (say) the delta T to look for convergence 

2.  This is because your mass flow controller between the mixer and the combustor has a mass flow rate of 0, check here.

You have to give it either a constant mass flow rate, or function, etc.

Nick

[Vaibhav Prakash]

Hello Nick and Hello all,

Thanks a lot for your reply. Your suggestion worked well and helped combust my mixture.

Now, I have a new problem, when I change the mass flow rate of the oxidizer, mixer volume and combustor volume.  My inlet pressure is at P=12.55 bar and inlet Temperature T=633 K. The internal time stepping sim.step() does not seem to work in this condition.

To re-iterate my problem statement:

 I am using Cantera 2.3.0 and Python 3.4 in my Windows 10 System. I am trying to build the system with the pre-mixer (upstream of combustor) as an Ideal Gas Reactor and the main combustor as a reactor Network (CSTR).

I am using Natural Gas as the fuel with GRI3.0 Mechanism. The premixed mixture's entry conditions into the combustor is now at 633 K and 12.55 bar. As a start, I am trying to connect the premixer (with 2 inlets- Natural Gas and EGR oxidizer) with the main combustor (different volume). I am igniting the mixture as per combustor.py wiht a stream of H radicals. Initially, I have assumed the whole combustor volume as a single Ideal Gas reactor.

So the outline is basically:

'''
Fuel stream pipe----|
                     --------->  PREMIXER (Ideal Gas Reactor) ---------> Combustor (CSTR)
EGR oxidiser pipe---|

''' 

When I try to run the program, I get the following error,

Traceback (most recent call last):

  File "E:/Cantera/bin/NG_air combustion.py", line 113, in <module>

    tnow = sim.step()

  File "interfaces\cython\cantera\reactor.pyx", line 937, in cantera._cantera.ReactorNet.step (interfaces\cython\cantera\_cantera.cpp:68465)

cantera._cantera.CanteraError: 

***********************************************************************

CanteraError thrown by CVodesIntegrator::step:

CVodes error encountered. Error code: -3

At t = 0.578136 and h = 2.32606e-09, the error test failed repeatedly or with |h| = hmin.

Components with largest weighted error estimates:

56: 175.074

58: 51.2686

65: -0.0183529

78: 0.0173774

63: 0.00922252

71: -0.00739355

76: -0.0064666

64: -0.00483891

60: 0.00416556

83: -0.00177641

***********************************************************************

Process finished with exit code 1

Irrespective of what equivalence ratio I put it in, the program does not calculate the internal time step sim.step(). Can you please help me and suggest some changes on where I am going wrong? What is the reason for the time stepping to fail when the network is created right? The pre-mixer reactor seems to work fine when run independently.

Kind Regards,

Vaibhav 

[Nick Curtis]

Hi Vaibhav,

This is a recoverable CVODE error, usually due to very fast (i.e. stiff) components in the solution.
The reason is that your mass flow is nearly equal to the mass in the reactor (right around when simulation breaks).

Putting the following print statement in the for loop:

print(tnow, combustor.mass, mdot_fuel)

Results in the following just before the break

0.5894000656800016 0.0014579644914642497 0.0014108698319224636

See how the combustor mass and the mdot are almost equal?  This forces some very small timescales in the solution.
Note that you can also get that info from the CVODE error:

print(sim.component_name(56))

-> IdealGasReactor_2: mass

I was able to get your simulation to run by fiddling with tolerances, but there is some very non-physical behavior (the temperature shoots up to ~4600K at points), hence you probably need to recheck your flow rates, combustor size, etc.

Nick

[Ray Speth]

Vaibhav,

The main problem, which seems to be quite a common one, is the value you have specified for the exhaust valve coefficient (K=100). This value is in (kg/s)/Pa -- i.e. a 1 Pa pressure difference generates a 1 kg/s mass flow out of the reactor. This introduces a lot of unnecessary stiffness to your system, and reducing this by a couple of orders of magnitude will help a lot.

The other issue is the mass flow rate for your "igniter" which is orders of magnitude larger than your fuel and air flow rates, which means that the dynamics of this igniter and its mass flow rate as a function of time dominate the computation. A better (and actually far simpler) solution is to (1) just use one reactor and set the composition of the incoming fuel/air mixture using the 'set_equivalence_ratio' method of the Solution class and (2) Set the initial condition for the reactor volume to correspond to the fuel/air mixture at equilibrium at constant enthalpy and pressure, and remove the simulated "igniter".

Regards,

Ray

On Friday, March 17, 2017 at 9:19:02 AM UTC-4, Vaibhav Prakash wrote:

[Vaibhav Prakash]

Dear Nick and Ray,

Thanks a lot for your suggestions. Your tips were really useful.

Ray,

 

To reduce the stiffness in the problem, I have scaled down the amplitude of the Gaussian Pulse by 100. This in turn gives me realistic temperatures. Hope this is also a right approach.

New Question:

I am now trying to integrate my PSR network (Mixer + Combustor + Recirculation) with a PFR.  The program works fine. But I am not getting the desired Temperature profile in the PFR i.e, a reducing temperature with increase in position. I just get a plot that reaches the equilibrium temperature and stays constant.

The program is split into 4 functions. main() , combustor_network_air() -- [PSR network] , pfr() --- PFR , T_ad_air [Adiabatic Flame Temperature]

1. Unlike the definition in pfr.py, I do not have data of the inflow velocity. In that example, the mass flow rate has been calculated using the inflow velocity. On the other hand, in my case, I know the incoming mass flow rate (mdot_fuel+mdot_oxidizer). Why is the desired reducing temperature trend not occurring? I have attached the image for the trend I get. Is it because of the mass flow rate? How to fix this problem? Where is it going on?

Kind Regards,

Vaibhav

[arief...@gmail.com]

.

Dear Vaibhav,

To cope with the stiffness problem you need two things:

(1) A Rosenbrock method with a proper Milne-device for adaptive timestep adjustment.

(2) Singular Value Decomposition to cope with the linear algebra of the ODE-solver.

Arief

From: Vaibhav Prakash <vaibha...@gmail.com>
Date: Mon, Apr 3, 2017 at 8:13 PM
Subject: [cantera-users] Re: Reactor Networks: Gas Turbine Premixed Combustor- Division by zero error and Ignition problem
To: Cantera Users' Group <canter...@googlegroups.com>

.

[Vaibhav Prakash]

Hi Arief,

Thank you for the swift response. I have managed to eliminate the stiffness problem by playing around with the amplitude of the pulse. Do you have an idea of what is going wrong in my PFR code that I have attached in my previous post?

Kind Regards,

Vaibhav

[Ray Speth]

Vaibhav,

Why do you expect the temperature to decrease along the length of the PFR? For an adiabatic PFR, the mixture should tend toward the equilibrium composition and temperature along its length.

Regards,

Ray

[arief...@gmail.com]

.

Hello Vaibhav,

I checked your code and noticed the following two code-fragments.

(A) An implementation of a Milne-device for solution-adaptive stepsize adjustment.

    # time stepping
    tnow = 0.0
    tfinal = 6.0
    Tprev = combustor.T
    tprev = tnow
    states = ct.SolutionArray(combustor.thermo)
    k = 0
    while tnow < tfinal:
        tnow = sim.step()
        Tnow = combustor.T
        if abs(Tnow - Tprev) > 1.0 or tnow - tprev > 2e-2:
            tprev = tnow
            Tprev = Tnow

(B) Application of an integrator

    # iterate through the PFR cells
    for n in range(n_reactors):
        # Set the state of the reservoir to match that of the previous reactor
        gas.TDY = r.thermo.TDY
        upstream.syncState()
        # integrate the reactor forward in time until steady state is reached
        sim2.reinitialize()
        sim2.advance_to_steady_state()

First we need to clarify which method is applied for the time-integration. And then figure out what might be a proper Milne-device for that method.
 

The PFR-problem you're working on is similar to Test Problem C in the file at http://www.hysafe.org/science/eAcademy/docs/FlowTurbulenceAndCombustion_v91_p281to317.pdf

Only the 2nd order and 7th order Gear method (as implemented in the SUNDIALS package) and the Rosenbrock method (as described in the paper) were able to solve this test problem. All other methods I investigated failed on this problem. And this is still a simple problem compared to the one you're working on.

Perhaps you could first clarify which time-integration method you're applying. And maybe we should resort to a Milne-device as per equation (81) of the paper.

Arief

From: Vaibhav Prakash <vaibha...@gmail.com>
Date: Tue, Apr 4, 2017 at 8:20 AM
Subject: Re: [cantera-users] Re: Reactor Networks: Gas Turbine Premixed Combustor- Division by zero error and Ignition problem
To: Cantera Users' Group <canter...@googlegroups.com>
Cc: vaibha...@gmail.com

Hi Arief,

Thank you for the swift response. I have managed to eliminate the stiffness problem by playing around with the amplitude of the pulse. Do you have an idea of what is going wrong in my PFR code that I have attached in my previous post?

Kind Regards,

Vaibhav

On Monday, 3 April 2017 23:56:12 UTC+2, Arief Dahoe wrote:

.

Dear Vaibhav,

To cope with the stiffness problem you need two things:

(1) A Rosenbrock method with a proper Milne-device for adaptive timestep adjustment.

(2) Singular Value Decomposition to cope with the linear algebra of the ODE-solver.

Arief

From: Vaibhav Prakash <vaibha...@gmail.com>
Date: Mon, Apr 3, 2017 at 8:13 PM
Subject: [cantera-users] Re: Reactor Networks: Gas Turbine Premixed Combustor- Division by zero error and Ignition problem
To: Cantera Users' Group <canter...@googlegroups.com>

.

Dear Nick and Ray,

Thanks a lot for your suggestions. Your tips were really useful.

Ray,

 

To reduce the stiffness in the problem, I have scaled down the amplitude of the Gaussian Pulse by 100. This in turn gives me realistic temperatures. Hope this is also a right approach.

New Question:

I am now trying to integrate my PSR network (Mixer + Combustor + Recirculation) with a PFR.  The program works fine. But I am not getting the desired Temperature profile in the PFR i.e, a reducing temperature with increase in position. I just get a plot that reaches the equilibrium temperature and stays constant.

The program is split into 4 functions. main() , combustor_network_air() -- [PSR network] , pfr() --- PFR , T_ad_air [Adiabatic Flame Temperature]

1. Unlike the definition in pfr.py, I do not have data of the inflow velocity. In that example, the mass flow rate has been calculated using the inflow velocity. On the other hand, in my case, I know the incoming mass flow rate (mdot_fuel+mdot_oxidizer). Why is the desired reducing temperature trend not occurring? I have attached the image for the trend I get. Is it because of the mass flow rate? How to fix this problem? Where is it going on?

Kind Regards,

Vaibhav

--
You received this message because you are subscribed to the Google Groups "Cantera Users' Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cantera-users+unsubscribe@googlegroups.com.
To post to this group, send email to canter...@googlegroups.com.
Visit this group at https://groups.google.com/group/cantera-users.
For more options, visit https://groups.google.com/d/optout.