Accessing diffusion velocities of a FreeFlame object

[Pablo Kandel]

Dear all,
I was wondering if there is a way to compute the diffusion velocities of each species in Cantera 3.0 in Python from a FreeFlame object. I didn't find any function doing that so far.
I would like to compute the diffusion velocities for the "mixture-averaged" transport model, the transport equation of the 1D flame problem are described here:
https://cantera.org/science/flames.html

Maybe it is possible to access the diffusive mass flux for each species along the 1D domain? I haven't succeed here neither ...

Here is a simplified piece of code computing a 1D freely-propagating flat flame using h2o2.yaml chemistry mechanism:

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

from pathlib import Path
import cantera as ct
import csv
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

# Simulation parameters
p = 101325  #pressure [Pa]
Tin = 300.0  #unburned gas temperature [K]
phi = 1
reactants = {'H2':0.01731720726, 'O2':0.2290635881, 'N2':0.7536192047} 
width = 0.04  # m
loglevel = 1  # amount of diagnostic output (0 to 8)

gas = ct.Solution('h2o2.yaml')
gas.TPY = Tin, p, reactants

# Set up flame object
f = ct.FreeFlame(gas, width=width)
f.set_refine_criteria(ratio=15.0, slope=0.01, curve=0.12)
f.show()

# Solve with mixture-averaged/unity-Lewis-number transport model
f.transport_model = 'mixture-averaged'
f.solve(loglevel=loglevel, auto=True)

output = Path() / "cantera.csv"
output.unlink(missing_ok=True)

# Solve with the energy equation enabled
f.save(output, name="Mix", description="solution with Mix transport")

reader = csv.reader(open("cantera.csv", "r"), delimiter=',')

f.show()
print(f"flamespeed = {f.velocity[0]:7f} m/s")

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

Thanks for your help!

Pablo

[Ray Speth]

Hi Pablo,

Unfortunately, this information isn’t accessible from the Python interface. There is internally an array at the C++ level, Flow1D::m_flux which stores this, but there’s currently no public accessor for it.

Regards,
Ray

​

[Pablo Kandel]

Hi Ray,
A simple "trick" to compute the diffusion velocity of H2 (here just with Fick's law and no correction velocity), first computing the mass fraction gradient and then multiplying by the diffusivity (red part of the code):



///////////////////////////////////////////////////////////////////////////////

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Tue Oct  8 15:02:31 2024

@author: pablo

"""
from pathlib import Path
import cantera as ct
import csv
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

phi = 1

# Simulation parameters
p = 101325

Tin = 300.0
reactants = {'H2':1, 'O2':0.5/phi, 'N2':0.5*3.76/phi} # premixed gas composition
width = 0.06  # m

loglevel = 1  # amount of diagnostic output (0 to 8)

# Solution object used to compute mixture properties, set to the state of the
# upstream fuel-air mixture

gas = ct.Solution('h2o2.yaml')

gas.TPX = Tin, p, reactants

# Set up flame object
f = ct.FreeFlame(gas, width=width)

f.set_refine_criteria(ratio=30.0, slope=0.0085, curve=0.12)
#f.set_refine_criteria(ratio=10.0, slope=0.05, curve=0.2, prune=0)

f.show()

# Solve with mixture-averaged/unity-Lewis-number transport model
f.transport_model = 'mixture-averaged'
f.solve(loglevel=loglevel, auto=True)

#Save results

output = Path() / "cantera.csv"
output.unlink(missing_ok=True)

f.save(output)

reader = csv.reader(open("cantera.csv", "r"), delimiter=',')

#Laminar flame speed

print(f"flamespeed = {f.velocity[0]:7f} m/s")

#Cantera data
file  = "cantera.csv"
flame = pd.read_csv(file)

#Species mass fraction gradient
size = len(flame['grid'])
grad_YH2 = np.zeros(size)
for n in range(size - 1):
grad_YH2[n] = (flame['Y_H2'][n + 1] - flame['Y_H2'][n]) / (flame['grid'][n + 1] - flame['grid'][n])

#Diffusion velocity, Fick's law, no correction velocity
VdH2 = -f.mix_diff_coeffs[gas.species_index('H2')]*grad_YH2/flame['Y_H2']

#Species
fig, ax = plt.subplots()
ax.plot(flame['grid']*1e2, VdH2, color="black", linestyle="-")
ax.set_xlim(4.1, 4.4)
#ax.set_ylim(0, 0.25)
ax.set_xlabel('x [cm]')
ax.set_ylabel('$Vd_{H2} [m/s]$')
ax.grid()
plt.title("H2 diffusion velocity")
plt.show()

///////////////////////////////////////////////////////////////////////////////