imposing a temperature profile on a 0-d simulation

[gpsimms]

Hi there,

I see there have been a few posts on a similar topic, so I am sorry if this is redundant, but I am still having issues even after having read the other posts.  I have a plug flow reactor which I have measured a temperature profile for different flowrates and set temperatures.  Instead of using the energy equation, I plan to turn the energy equation off and impose the experimentally measured T profile.  I'm new to cantera, so the 0-d IdealGasConstPressureReactor is the method I was trying to use.

The temperature profile is discontinuous, and I am wondering if that is part of the problem.  I have attached my code, and some plots.

import sys

import numpy as np

import cantera as ct

import csv

path='JP10_1000K_0.5_SLPM_profile.txt'

prof=np.loadtxt(path, delimiter="\t")

reaction_mechanism='RMG_jp10_v0_19_SM.cti'

gas = ct.Solution(reaction_mechanism)

gas.TPX = prof[0,1], ct.one_atm, 'JP10(1):0.0025,Ar:0.9975'

step=1e-5

r = ct.IdealGasConstPressureReactor(gas, energy='off')

sim = ct.ReactorNet([r])

time = 0.0

states = ct.SolutionArray(gas, extra=['t'])

#xit=states

print('%10s %10s %10s %14s' % ('t [s]','T [K]','P [Pa]','u [J/kg]'))

tstep_0 = prof[0,0]

while time < tstep_0:

    sim.advance(time)

    time+=step

    states.append(r.thermo.state, t=time)

i=0

for temp in prof[1:4,1]:

    

    gas.TP = temp, None

    r.syncState()

    i=i+1

    while time < prof[i,0]:

        time+=step

        sim.advance(time)

        

        states.append(r.thermo.state, t=time)

        

    

The error I get is that ydot contains non-finite elements.  This code runs correctly if I keep the number of steps  very low, but 'blows up' quickly.  Here is a plot of temperature in the first few time steps:

I plan on linearizing/smoothing the temperature profile, but I don't think the discontinuities in temperature should be causing an issue. Am I wrong about that?  Here is a plot of C2H4 production in time:

It seems like if I take another step forward in time, many species blow up (in the negative direction). Is there something wrong with the way I am imposing energy=off, but also imposing a T profile?  Am I taking a completely wrong approach, do I need to instead use the PFR code?  

Thanks!

greg

[walterm...@me.com]

Hi Greg,

I didn't perfectly understand what you are trying to do but it seems that a class that I wrote with Cantera is suitable for your needs.
You can find it under my github https://github.com/waltermateriais/xphysics (you will need the whole package because I reuse some formats across files).
An example of PFR is presented on https://github.com/waltermateriais/xphysics/blob/master/tests/test_xcantera.py line 158.
The class allows for lambda function temperature profile, which can be set through the dictionary element rsetup['operation']['temperature'].
This lambda function (or any arbitrary function) shall be a function of the position only.

Please find below a snippet that could further illustrate how to use the class:

# ----------------------------------------------------------------------------

# class plug_flow_reactor

# ----------------------------------------------------------------------------

from xphysics.xcantera import plug_flow_reactor

# Wall temperature parameters (see function definition below)
pars = {
    '773'  : (0.04132785, 0.36586941, 1.92089872, 12.41516606),
    '873'  : (0.03457862, 0.39032227, 1.41582889, 9.79102679),
    '973'  : (0.02537489, 0.39703098, 0.99659743, 9.77523826),
    '1023' : (0.02528152, 0.40339555, 0.88494798, 10.55513796),
    '1073' : (0.02507178, 0.40847247, 0.81631547, 11.98899245),
    '1123' : (0.02497517, 0.40832661, 0.80065655, 11.97005813),
    '1173' : (0.02492942, 0.40810172, 0.78913918, 11.91548263),
    '1223' : (0.02596356, 0.40572591, 0.85168097, 11.01722351),
    '1273' : (0.02682903, 0.40342913, 0.91051192, 10.36909121)
}

# Conversion mbar->Pa
mbar = 100

# Reactor setup
rsetup = {
    'geometry'  : {
        'radius' : 0.014,
        'length' : 0.450,
        'slice'  : 0.001,
    },
    'operation' : {
        'pressure'    : None,
        'temperature' : None,
        'flowrate'    : None
    },
    'inletgas'  :{
        'N2'      : 1.0-0.36/0.98,
        'C2H2'    : 0.36,
        'CH3COCH3': 0.36/0.98*1.8e-02,
        'CH4'     : 0.36/0.98*2.0e-03
    },
    'hommech'   : 'norinaga2009.cti',
    'follow'    : ('C2H2', 'CH4', 'H2'),
    'overwrite' : False,
    'preheat'   : False,
    'energy'    : 'on',
    'transmod'  : 'Mix',
    'relerr'    : 1.0e-09,
    'abserr'    : 1.0e-20,
    'tstep'     : 5.0e-03,
    'epsx'      : 1.0e-09,
    'dxm'       : 0.005,
    'maxiter'   : 500
}


# Define the reactor temperature profile.
def Twall(x, Treg):
    a1, a2, m1, m2, k = *pars[str(int(Treg))], 1.0
    st_trm, nd_trm = np.exp(-(x/a1)**m1), np.exp(-(x/a2)**m2)
    return 300+(k*Treg-300)*(1.0-st_trm)-(k*Treg-400)*(1.0-nd_trm)


## P = 30 mbar
rsetup['operation']['flowrate'] = 222.0
rsetup['operation']['pressure'] = 30*mbar

for T in [1123, 1173, 1223]:
    try:
        rsetup['operation']['temperature'] = T
        rsetup['walltmp'] = lambda x : Twall(x, T)

pfr = plug_flow_reactor(rsetup)

        pfr.plot(['C2H2'])
    except Exception as err:
        print(err)
        pass

Sorry for the missing documentation, I just passed my thesis defense and I am working right now fixing some bugs and generalizing this class for my current needs.


Walter

[gpsimms]

Thanks Walter.

I had seen your other discussion on this topic before, but since I am not a python expert I thought my approach seemed 'simpler.'  But, since my approach was giving me issue I will give your package a try and let you know problems as I come across them. Congrats on passing your defense! No light at the end of the tunnel yet for me, but one can hope...

greg

[gpsimms]

Hey there Walter,

I tried sending you a direct email, but maybe you'd prefer I didn't do that. As I said, I'm basically completely new to python, and only just started learning it when our chemkinpro license expired.

So, off the bat, if I try to run xcantera.py, I get the message that there is no module named 'xphysics.'  When I look in the xphysics directory, I notice that __init__.py says that windows machines are not supported.  Is this my problem?

I am able to remotely log into a unix machine, but will need to install python and cantera to it.

Thanks,

greg