# premixed_counterflow_twin_flame.py (Source)

```# coding: utf-8

"""
Simulate two counter-flow jets of reactants shooting into each other. This
simulation differs from the similar premixed_counterflow_flame.py example as the
latter simulates a jet of reactants shooting into products.
"""

import cantera as ct
import numpy as np
import sys

# Differentiation function for data that has variable grid spacing Used here to
# compute normal strain-rate
def derivative(x, y):
dydx = np.zeros(y.shape, y.dtype.type)

dx = np.diff(x)
dy = np.diff(y)
dydx[0:-1] = dy/dx

dydx[-1] = (y[-1] - y[-2])/(x[-1] - x[-2])

return dydx

def computeStrainRates(oppFlame):
# Compute the derivative of axial velocity to obtain normal strain rate
strainRates = derivative(oppFlame.grid, oppFlame.u)

# Obtain the location of the max. strain rate upstream of the pre-heat zone.
# This is the characteristic strain rate
maxStrLocation = abs(strainRates).argmax()
minVelocityPoint = oppFlame.u[:maxStrLocation].argmin()

# Characteristic Strain Rate = K
strainRatePoint = abs(strainRates[:minVelocityPoint]).argmax()
K = abs(strainRates[strainRatePoint])

return strainRates, strainRatePoint, K

def computeConsumptionSpeed(oppFlame):

Tb = max(oppFlame.T)
Tu = min(oppFlame.T)
rho_u = max(oppFlame.density)

integrand = oppFlame.heat_release_rate/oppFlame.cp

I = np.trapz(integrand, oppFlame.grid)
Sc = I/(Tb - Tu)/rho_u

return Sc

# This function is called to run the solver
def solveOpposedFlame(oppFlame, massFlux=0.12, loglevel=1,
ratio=2, slope=0.3, curve=0.3, prune=0.05):
"""
Execute this function to run the Oppposed Flow Simulation This function
takes a CounterFlowTwinPremixedFlame object as the first argument
"""

oppFlame.reactants.mdot = massFlux
oppFlame.set_refine_criteria(ratio=ratio, slope=slope, curve=curve, prune=prune)

oppFlame.show_solution()
oppFlame.solve(loglevel, auto=True)

# Compute the strain rate, just before the flame. This is not necessarily
# the maximum We use the max. strain rate just upstream of the pre-heat zone
# as this is the strain rate that computations comprare against, like when
# plotting Su vs. K
strainRates, strainRatePoint, K = computeStrainRates(oppFlame)

return np.max(oppFlame.T), K, strainRatePoint

# Select the reaction mechanism
gas = ct.Solution('gri30.cti')

# Create a CH4/Air premixed mixture with equivalence ratio=0.75, and at room
# temperature and pressure.
gas.set_equivalence_ratio(0.75, 'CH4', {'O2':1.0, 'N2':3.76})
gas.TP = 300, ct.one_atm

# Set the velocity
axial_velocity = 2.0 # in m/s

# Domain half-width of 2.5 cm, meaning the whole domain is 5 cm wide
width = 0.025

# Done with initial conditions
# Compute the mass flux, as this is what the Flame object requires
massFlux = gas.density * axial_velocity # units kg/m2/s

# Create the flame object
oppFlame = ct.CounterflowTwinPremixedFlame(gas, width=width)

# Uncomment the following line to use a Multi-component formulation. Default is
# mixture-averaged
#oppFlame.transport_model = 'Multi'

# Now run the solver. The solver returns the peak temperature, strain rate and
# the point which we ascribe to the characteristic strain rate.

(T, K, strainRatePoint) = solveOpposedFlame(oppFlame, massFlux, loglevel=1)

# You can plot/see all state space variables by calling oppFlame.foo where foo
# is T, Y[i], etc. The spatial variable (distance in meters) is in oppFlame.grid
# Thus to plot temperature vs distance, use oppFlame.grid and oppFlame.T

Sc = computeConsumptionSpeed(oppFlame)

print("Peak temperature: {0:.1f} K".format(T))
print("Strain Rate: {0:.1f} 1/s".format(K))
print("Consumption Speed: {0:.2f} cm/s".format(Sc*100))
oppFlame.write_csv("premixed_twin_flame.csv", quiet=False)

# Generate plots to see results, if user desires
if '--plot' in sys.argv:

import matplotlib.pyplot as plt

plt.figure(figsize=(8,6), facecolor='white')

# Axial Velocity Plot
plt.subplot(1,2,1)
plt.plot(oppFlame.grid, oppFlame.u, 'r', lw=2)
plt.xlim(oppFlame.grid, oppFlame.grid[-1])
plt.xlabel('Distance (m)')
plt.ylabel('Axial Velocity (m/s)')

# Identify the point where the strain rate is calculated
plt.plot(oppFlame.grid[strainRatePoint], oppFlame.u[strainRatePoint],'gs')
plt.annotate('Strain-Rate point',
xy=(oppFlame.grid[strainRatePoint], oppFlame.u[strainRatePoint]),
xytext=(0.001, 0.1),
arrowprops={'arrowstyle':'->'})

# Temperature Plot
plt.subplot(1,2,2)
plt.plot(oppFlame.grid, oppFlame.T, 'b', lw=2)
plt.xlim(oppFlame.grid, oppFlame.grid[-1])
plt.xlabel('Distance (m)')
plt.ylabel('Temperature (K)')

plt.tight_layout()
plt.show()

else:
print('************')
print('Plotting option not enabled. Re-run script with --plot to see key plots.')
print('************')
```