HCCI Heat Transfer Modeling

[Shane Ryan Daly]

Hello, 

I have created an HCCI simulation that has successfully matched that of Chemkin perfectly (without heat transfer enabled).  Attached is a simple plot of in-cylinder pressure of PRF0 at 800 RPM, Ti=[350, 450, 550] K, Pi=1bar, phi =0.35 showing Chemkin (dotted line) overlaying perfectly with my Python/cantera model (file: PRF 0 350-550K...) .  Moving on from there I decided I would add in Woschni heat transfer, which is also a feature of Chemkin.  For brevity, I have attached the Chemkin theory manual PDF to show all the equations to do this (Pg. 150-152).  For simplicity, I greatly reduced the complexity of the correlation for the average gas cylinder velocity to be simply C11*Sp.  Therefore, Reynolds number can be calculated, from there the Nusselt number, and finally the heat transfer coefficient. For Re/Nu property values I used the air mechanism specified at the initial conditions of each simulation (each different temperature, as Chemkin does explained in the theory manual).  Therefore, the convective heat transfer coefficient is not a function of crank angle.  From here, a function for piston position vs. crank angle can be defined and used to calculate crank resolved heat transfer area.  Lastly, newtons law of cooling can be utilized, where the temperature of the cylinder wall has been set to a constant value (430K), and the in-cylinder gas temperature has been referenced from the reactor (r2.T, where r2 is my ideal gas reactor).  The heat transfer area has been divided by the piston cross sectional area to cancel out the fact that I have defined the piston cross sectional area as my wall area (to go along with the wall velocity to have the correct crank resolved reactor volume).  

All in all, I feel as if I have done everything correctly and matched that of Chemkin. Unfortunately, the results make no physical sense.  To help illustrate, I have attached a plot of the reactor temperature for an adiabatic case. Clearly, the temperature starts at 350K (below the wall temperature, so I will expect negative heat transfer), and from there it is going to increase, and slightly decrease during the expansion stroke (but always above the wall temperature, expecting positive heat transfer).  The heat loss--also attached (Python along with correct Chemkin results; again, Chemkin is the dotted line)--is seen to be highly negative after around 190 CAD, which makes no sense.  In addition to the results, I have attached a plot of the piston position versus time which is used for the area calculation --nothing wrong there, always a positive value that ranges from the stroke at BDC to the clearance volume length at TDC.  I think I can conclude that something wrong is going on with the r2.T temperature; how else could the heat flux be negative?  Even looking at the temperature profile for the non-adiabatic scenario the temperature is far above the wall temperature as expected (this is also attached).  I have ruled out that this is an issue with the forking parallel processing executor that I am using.  Simply running the function through the command prompt yields the same results. 

After running the simulation the performance parameters can be observed as follows: 

(with parallel process)

HCCI_solution_collected[Parallel process indice (use 0)][RPM indice (0 for 800 RPM) ][Temperature Indice (0 for 350K)][Phi indice (0 for 0.35)][Performance Parameter (15 for heat loss)]

here is a list of the performance parameters if you want to look at anything else other than heat transfer: 

output_variables=['Xf','Xo2','Xco','Xco2','Xn2' ,'Xh2o', 'pres','temp','volume','crank_rotation_angle','net_heat_release_rate_per_CA','mass,molecular_weight','mass_density','specific_heat','mixture_enthalpy','heat_loss']

If executing without parallel process executor: 

Solution=HCCI(RPM value)

Solution[Temperature Indice][Phi Indice][Performance Parameter Indice]

           

Lastly, I am finding a few odd bugs.  During the simulation time advance, I simply collect values such as temperature, pressure, and my heat loss function.  The first bug is that I cannot do any mathematical operations on q(time).  Meaning, that if I were to say multiply q(t) by 10^-6, it still will not show up in the results. I came about this discovery because I wanted to multiply q(t) by the piston cross sectional area.  I'm not recording the heat transfer after the cantera code multiplies by the piston cross sectional area; I'm recording my original q(t) function which is divided by that area to cancel the fact that cantera will multiply by the wall area.  I noticed that the heat loss was not being corrected by the area because it was orders of magnitude too large, so I tried multiplying by an incredibly small number (10^-6)--no change in the results!  The second issue is that I have to multiply all my other performance parameters (pressure, temperature, etc..) by one in order for everything to be collected properly.  Without this operation, only the values corresponding to my last iteration was saved for my entire study. In other words, all the performance parameters would correspond to 450K, and phi=0.45.  The other values for 350K and phi=.35 were essentially erased.  Funny enough, I stumbled upon this error because I was multiplying pressure by 10^-5 to convert from Pa to bar. All the values--except for pressure-- corresponded to the high temperature and phi case.  Pressure, however, was correct in that it represented each case I conducted.  So, I decided to multiply all my values by 1 and it fixed that issue...

Any light that could be shed on these issues would be very much appreciated!  

My code can be seen below, where items highlighted are the lines of importance.  The Princeton chemical/therm mechanism that I am using is also attached.    

##############################################################################################

#Chemical Mechanism: Princeton n-heptane, iso-octane,ethanol, and Toluene reference fuel (10 atm) 

#Version: Python 2.7.

#OS: Windows 8

# -*- coding: utf-8 -*-

"""

Created on Wed Dec 10 12:00:03 2014

@author: Shane

#Contact: da...@onid.oregonstate.edu 

#0-D Homogeneous Charge Compression Ignition (HCCI) Engine Simulation

"""

from __future__ import division    #Utilize Python version 3 features

import numpy                #"Numpy" for Mathematical operations

import cantera              #import cantera library

from joblib import Parallel, delayed

num_cores=6 #cores used in parallel processing 

#-------------->Parametric Study Values<--------------

RPM=numpy.array([800.0, 1000.0, 1200.0, 1400.0, 1600.0, 1800.0])

Pi=100000 #Initial Pressure [Pa] =1 bar (at BDC, not inlet intake air)

mechanism=cantera.Solution('Princeton_10atm.cti') #mechanism import

#--------------------->HCCI Simulation Function Code<--------------------------------

def HCCI(RPM):

    phi=[0.35, 0.45]    

    Ti=[350.0,450.0]    

    length_phi=len(phi)    

    ##########################################################################################################        

    #------------->importing Princeton_10atm and air chemical/thermophysical mechanism<-------------------    

    air=cantera.Solution('air.xml') #air (the piston model is air separated by a wall(piston) to the reactor)  

    X_air=[0.0,0.21,0.0,0.0,0.0,0.0,0.79,0.0] #setting initial composition to O2 and N2 only (Like Chemkin)    

    io2 =mechanism.species_index('o2')  #Finding species indices 

    in2 = mechanism.species_index('n2') #Finding species indices 

    ico2 =mechanism.species_index('co2')#Finding species indices

    ico =mechanism.species_index('co')  #Finding species indices  

    ih2o = mechanism.species_index('h2o')#Finding species indices     

    #------------------->Engine and Heat Transfer Parameters<--------------------- 

    rc=9.5                              #Compression Ratio [Vmax/Vmin] (Given)     

    B=.0828                             #Cylinder Bore (Given)  

    Vd=616*10**(-6)                     #Displacement Volume of Engine [m^3] (Given)

    c=4.44                              #Connecting rod to crank radius ratio (Given) [m/m]

    s=4*Vd/(numpy.pi*B**2)              #Stroke [m]

    a=s/2                               #Crank radius [m]

    l=a*c                               #connecting rod length [m]

    Vmax=Vd*rc/(rc-1)                   #Maximum cylinder Volume [m^3]

            #Engine parameters for heat transfer model   

    Ap=(numpy.pi*B**2)/4                #Piston Area

    Vc=(Vmax/(rc-1))                    #Clearance Volume (For heat transfer model)

    Vc_length=Vc/Ap                     #length of clearence volume required given bore area (for heat transfer model)

    wdot=RPM*numpy.pi/30    

    #------------------->Setting resolution of simulation<-----------------------     

    timesteps=286

    t=numpy.linspace(0,60/RPM,timesteps) 

    theta=t*wdot*180/numpy.pi

    #slider crank rule for piston velocity. Cantera calls for piston(reactor wall) velocity, not Volume(t) or dV/dt (for the energy equation)

    def piston_speed(t):

        return (wdot*s/2*numpy.sin(t*wdot))*(1-(numpy.cos(t*wdot))/(numpy.sqrt(c**2-numpy.sin(t*wdot)**2)))      

    #Heat Transfer Model (Woschni) -->All values are taken from air at the specified inlet conditoins (Just as Chemkin Does)

    aa=.035                             #Nusselt Number Correlation Constant

    bb=.8                               #Nusselt Number Correlation Re power constant

    cc=0                                #Nusselt Number Correlation Pr power constant

    C11=2.28                            #from Chemkin

    Tw=430                              #Wall (Bore) Temperature             

    #Piston Position correlation (used for heat transfer Area calculation) 

    def position(t):

        return a*numpy.cos(wdot*t)+numpy.sqrt(l**2-a**2*numpy.sin(wdot*t)**2)

    #Heat Transfer Correlation 

    def q(t):

        return h*(r2.T-Tw)*(2*Ap+numpy.pi*B*(position(t)-l+a+Vc_length))/Ap  #piston and top area (not including clearance volume cylinder area) and the cylinder area affected from piston motion. Divide by Ap, cantera wants W/m^2, where canterea area is defined as Ap              

    #-------------->Fuels and composition Parameters (Solves for Stoich F/A micture)<-----------------------

    fuel_species = ['nc7h16']           #Naming the fuel species mixture  

    length_f=len(fuel_species)

    X_fuel=[1]        #Fuel composition (PRF0 Currently)

    ifuel=numpy.zeros(length_f)

    stoich_o2_fuel=numpy.zeros(length_f)

    for i in range(length_f):

        ifuel[i] =mechanism.species_index(fuel_species[i])

        stoich_o2_fuel[i]=X_fuel[i]*(mechanism.n_atoms(fuel_species[i],'C')+mechanism.n_atoms(fuel_species[i],'H')/4-mechanism.n_atoms(fuel_species[i],'O')/2)

    stoich_o2=numpy.sum(stoich_o2_fuel)    

    x = numpy.zeros((len(mechanism.X),length_phi))

    x[io2] = stoich_o2

    x[in2] = stoich_o2*(0.79/0.21)    

    #------------------->Air Reservoir and Pre-Allocating arrays for post-processing of solutions<-------------------

    env = cantera.Reservoir(air, name='Environment (Air at the specified heat transfer correlation temperature')

    pres = numpy.zeros((timesteps))

    Xf=numpy.zeros((timesteps))

    Xo2=numpy.zeros((timesteps))

    Xco=numpy.zeros((timesteps))

    Xco2=numpy.zeros((timesteps))

    Xn2=numpy.zeros((timesteps))

    Xh2o=numpy.zeros((timesteps))

    temp=numpy.zeros((timesteps))

    crank_rotation_angle = numpy.zeros((timesteps))

    #volumetric_heat_production_rate=

    net_heat_release_rate_per_CA = numpy.zeros((timesteps))

    mass = numpy.zeros((timesteps))

    #indicated_work= 

    molecular_weight = numpy.zeros((timesteps))

    mass_density = numpy.zeros((timesteps))

    specific_heat = numpy.zeros((timesteps))

    mixture_enthalpy = numpy.zeros((timesteps))

    volume = numpy.zeros((timesteps))

    heat_loss = numpy.zeros((timesteps))

    solution_RPM_T_phi=[]

    solution_RPM_T=[]

    #########################Initiation of Parametric Studies#################################

    for z in range(len(Ti)):            #Temperature Loop

        solution_RPM_T_phi=[] 

        #Heat Transfer Correlation Constants. Chemkin uses initial air properties for evaluation        

        air.TPX=Ti[z],Pi,X_air        

        mu=air.viscosity                #Viscosity of air

        k_g=air.thermal_conductivity    #Thermal Conductivity of air 

        rho=air.density                 #Density of Air 

        Pr=.7                           #Prandtl number (Doesnt matter at the moment, we set cc to zero...)

        Sm=2.0*a*wdot                   #Average Piston Speed

        w_bar=C11*Sm                    #Woschni average cylinder gas velocity

        Re=B*w_bar*rho/mu               #Reynolds Number

        Nu=aa*Re**bb*Pr**cc             #Nusselt Number

        h=Nu*k_g/(B)                    #Convective Heat Transfer Coefficient 

        for kk in range(len(phi)):      #Equivalence Ratio Loop

            for qq in range(length_f):  #The Fuel Composition Loop

                x[ifuel[qq],kk] = X_fuel[qq]*phi[kk]            

            mechanism.TPX=Ti[z],Pi,x[:,kk]

            r2 = cantera.IdealGasReactor(mechanism,  name='HCCI Reactor')

            r2.volume=Vmax #Specifying initial volume (BDC volume)    

            

            #####################HCCI CORE CODE#####################

                   

            cantera.Wall(env, r2,velocity=piston_speed,A=Ap,Q=q)

            sim = cantera.ReactorNet([r2]) 

            time=0

            n_steps=len(t)

            dt=t[2]-t[1] 

            sim.set_max_time_step(4*dt) #Sets maximum integration time-step             

            sim.rtol=1e-09              #Sets maximum integration relative tolerance 

            sim.atol=1e-20              #Sets maximum integration absolute tolerance         

            for n in range(n_steps):

                temp[n]=r2.thermo.T    

                pres[n]=r2.thermo.P

                volume[n] = r2.volume

                Xf[n]=r2.kinetics.X[ifuel[0]]

                Xo2[n]=r2.kinetics.X[io2]

                Xco[n]=r2.kinetics.X[ico]

                Xco2[n]=r2.kinetics.X[ico2]

                Xn2[n]=r2.kinetics.X[in2]

                Xh2o[n]=r2.kinetics.X[ih2o]

                crank_rotation_angle[n]=theta[n]

                net_heat_release_rate_per_CA[n]= r2.kinetics.enthalpy_mass

                mass[n]=r2.mass                

                molecular_weight[n]=r2.kinetics.mean_molecular_weight

                mass_density[n]=r2.kinetics.density

                specific_heat[n]= r2.kinetics.cp

                mixture_enthalpy[n]= r2.kinetics.enthalpy_mass

                heat_loss[n]=q(time)

                time += dt

                sim.advance(time)

                

            #####################END HCCI CODE#####################

    

            F= Xf*1,Xo2*1,Xco*1,Xco2*1,Xn2*1,Xh2o*1, pres*10**(-5),temp*1,volume*1,crank_rotation_angle*1,net_heat_release_rate_per_CA*1,mass,molecular_weight*1,mass_density*1,specific_heat*1,mixture_enthalpy*1,heat_loss*Ap

            solution_RPM_T_phi.append(F)

        solution_RPM_T.append(solution_RPM_T_phi)

    return  solution_RPM_T    

#---------------------->Parallel Process job Executor for HCCI Simulation<--------------------

#Send num_cores jobs at once, collect, then send another batch of num_cores studies

if __name__ == '__main__': #No idea what this does, but Windows machines need it to "protect" the main loop (whatever that is...) for parallel processing                      

    HCCI_solution_collected=[]   

    for j in range(int(len(RPM)/num_cores)): #Send num_cores RPM's Loop    

        HCCI_solution=Parallel(n_jobs=num_cores, verbose=100)(delayed(HCCI)(RPM[o+num_cores*j]) for o in range(num_cores)) #parallel process executor

        HCCI_solution_collected.append(HCCI_solution)

##############################################################################################################################

[Shane Ryan Daly]

I forgot to mention I am using Cantera-Python-2.1.2-x64-py2.7 and Cantera-2.1.2-x64

[Ray Speth]

Hi Shane,

I have answers for some of your questions, but there are a few others I'm not clear on.

1) Your heat transfer function has the wrong sign.

When you add the Wall object, you have it with "env" on the left and "r2" on the right, which means a positive heat transfer is from env to r2, but your heat transfer function is defined as q = h*(T(r2) - T(env)). You can see that there's a heat transfer into the cylinder from the fact that T is increasing as you approach the end of the expansion stroke.

2) Verifying the function of q(t):

Even with your original heat loss function, I can't reproduce the q(t) that you show in "Heat Loss.PNG". What conditions is that for? You can verify that r2.T is returning the correct value by capturing it inside your q(t) function (create empty lists for t, T, and Q outside your loop, then append to them inside the q(t) function).

3) why multiplying by one affects things:

In Python, when you write:

x = np.zeros(5)
y = np.zeros(5)
z = x,y

the tuple "z" contains references to the original "x" and "y" objects, not copies of them. By multiplying an array by a factor, you force the creation of a copy. What I would suggest is moving the initialization of these arrays into your "for" loop so that you get a new array for each simulation.

4) What do you mean by "cannot do any mathematical operations on q(time)."?

What exactly are you trying to do? Please show your code.

Regards,

Ray

[Shane Ryan Daly]

Thanks for the response Ray,

1) Even when I flip the sign in my heat transfer function, I still have the same heat transfer profile in that it is still indicating positive heat transfer.  Although, the answer to this mystery is answered in my followup to the suggestion you made in 2. 

2) First, the conditions that I had represented in my heat transfer function was Ti=350K and phi=0.35.  Anyway, I created a temperature array and then appended it with the referenced r2.T within the q(t) function as you suggested.  The results are shown in "r2.T collected temperature.PNG".  It looks as the the temperature is the correct profile, but it is inverted!  No wonder your profile is not matching mine, assuming you're not having the same issue. What could be causing this? 

3) Thank you for addressing how a tuple works, I understand now. 

4) By doing mathematical operations on q(t), I meant something like the following within the simulation"for" loop: 

heat_loss[n]=q(time)*( any number, perhaps piston cross sectional area) 

No matter the multiple it's not represented in the solution "heat_loss".

Thanks for your previous help and suggestions! I really appreciate the time you took to read my elaborate scenario. Next time I'll condense my problem down.  Again, thanks.  

-Shane     

[Ray Speth]

Shane,

The graph you attached looks like q(t) to me, not T(t), given that it contains a zero crossing and the non-monotonic behavior during the expansion stroke. The signs are correct assuming that this is the case with 350 as the initial temperature. How do you figure you're still getting (net) positive heat transfer?

Writing "heat_loss[n]=q(time)*( any number, perhaps piston cross sectional area)" will certainly modify the value of heat_loss[n], but it of course will not affect the results of the simulation. I'm not sure what's leading you to the conclusion that the multiplier is having no effect. Maybe you have multiple statements that set heat_loss and you're overwriting the value?

Regards,

Ray