wandadars
commented
1 year ago For reference, this is the control file that I've been using for testing the overall execution of the two-point flame control method. It's still not very stable, which I believe it should at least be somewhat stable based on the lack of any caveats present in the 1994 paper from Nishioka, Law, and Takeno (the paper that serves as the basis for this implementation).
Also attached is the mechanism file I am using to debug the implementation. I had to rename the YAML file to have a .txt suffix to make Github happy. burke.txt
import os
import time
import matplotlib.pyplot as plt
import numpy as np
import cantera as ct
def plot_solution(f, solution_num, control_points=False):
variables = [('temperature', f.T, 'Temperature (K)'), ('Uo', f.Uo, 'Oxidizer Velocity (m/s)'),
('L', f.L, 'Lambda (m/s)'), ('u', f.velocity, 'Axial Velocity (m/s)'),
('V', f.spread_rate, 'Spread Rate (1/s)')
]
for folder, data, ylabel in variables:
fig, ax = plt.subplots()
ax.plot(f.grid, data, 'o-')
if control_points and folder == 'temperature':
plt.axvline(x=f.zLeft, color='r', linestyle='--', linewidth=2, label='Left Point')
plt.axvline(x=f.zRight, color='g', linestyle='-.', linewidth=2, label='Right Point')
ax.legend()
ax.set_ylabel(ylabel)
plt.savefig(f'output/{folder}/flamelet_{solution_num}.png')
plt.close()
# Combined plot
fig, axs = plt.subplots(3, 2, figsize=(10,10))
for i, (folder, data, ylabel) in enumerate(variables + [('temperature_zoom', f.T, 'Temperature (K) Zoomed')]):
ax = axs.flat[i]
ax.plot(f.grid, data, 'o-')
ax.set_xlabel('Distance, m')
ax.set_ylabel(ylabel)
if folder == 'temperature_zoom':
high_temp_indices = np.where(data > 310)[0]
if len(high_temp_indices) > 0:
zoom_range_min = f.grid[high_temp_indices[0]]
zoom_range_max = f.grid[high_temp_indices[-1]]
ax.set_xlim([zoom_range_min, zoom_range_max])
plt.axvline(x=f.zLeft, color='r', linestyle='--', linewidth=2, label='Left Point')
plt.axvline(x=f.zRight, color='g', linestyle='-.', linewidth=2, label='Right Point')
ax.legend()
# Add a title to the figure
fig.suptitle(f'Flamelet Number: {solution_num}', fontsize=16)
plt.subplots_adjust(wspace=0.4, hspace=0.6)
plt.savefig(f'output/combined/flamelet_{solution_num}.png')
plt.close()
def plot_mixture_fraction_solution(f, solution_num, control_points=False):
variables = [('temperature_z', f.T, 'Temperature (K)'), ('Uo_z', f.Uo, 'Oxidizer Velocity (m/s)'),
('L_z', f.L, 'Lambda (m/s)'), ('u_z', f.velocity, 'Axial Velocity (m/s)'),
('V_z', f.spread_rate, 'Spread Rate (1/s)')
]
# Get the practical range for the mixture fraction for plotting only valuable data
epsilon = 1e-3
grid = f.mixture_fraction('Bilger')
# Find indices where grid crosses epsilon threshold
greater_than_epsilon_indices = np.where(grid > epsilon)[0]
less_than_epsilon_indices = np.where(grid < 1 - epsilon)[0]
if greater_than_epsilon_indices.size > 0 and less_than_epsilon_indices.size > 0:
first_index = greater_than_epsilon_indices[0]
last_index = less_than_epsilon_indices[-1]
subset_indices = range(first_index, last_index + 1)
else:
print("No indices found for the specified epsilon threshold")
return
for folder, data, ylabel in variables:
fig, ax = plt.subplots()
ax.plot(grid[subset_indices], data[subset_indices], 'o-')
if control_points and folder == 'temperature':
plt.axvline(x=f.zLeft, color='r', linestyle='--', linewidth=2, label='Left Point')
plt.axvline(x=f.zRight, color='g', linestyle='-.', linewidth=2, label='Right Point')
ax.legend()
ax.set_ylabel(ylabel)
ax.set_xlabel('Mixture Fraction(Bilger), Z')
plt.savefig(f'output/{folder}/flamelet_{solution_num}.png')
plt.close()
# Combined plot
fig, axs = plt.subplots(3, 2, figsize=(10,10))
for i, (folder, data, ylabel) in enumerate(variables + [('temperature_control', f.T, 'Temperature (K)')]):
ax = axs.flat[i]
ax.plot(grid[subset_indices], data[subset_indices], 'o-')
ax.set_xlabel('Mixture Fraction(Bilger), Z')
ax.set_ylabel(ylabel)
if folder == 'temperature_control':
plt.axvline(x=f.zLeft, color='r', linestyle='--', linewidth=2, label='Left Point')
plt.axvline(x=f.zRight, color='g', linestyle='-.', linewidth=2, label='Right Point')
ax.legend()
# Add a title to the figure
fig.suptitle(f'Flamelet Number: {solution_num}', fontsize=16)
plt.subplots_adjust(wspace=0.4, hspace=0.6)
plt.savefig(f'output/combined_z/flamelet_{solution_num}.png')
plt.close()
def create_output_folders():
folders = ['temperature', 'Uo', 'L', 'u', 'V', 'combined','temperature_z',
'Uo_z', 'L_z', 'u_z', 'V_z', 'combined_z', 'flamelets']
for folder in folders:
path = f'output/{folder}'
if not os.path.isdir(path):
os.makedirs(path, exist_ok=True) # exist_ok=True will prevent any error if the directory exists.
#Callback testing
def plot_variable(grid, data, title, xlabel, ylabel, filename):
plt.figure()
plt.plot(grid, data, 'o-', lw=2, label=ylabel)
plt.title(title)
plt.xlabel(xlabel)
plt.ylabel(ylabel)
plt.legend()
plt.savefig(filename)
plt.close()
def steady_callback(t, f):
variables = [
(f.T, 'Temperature (K)', 'temperature_plot_at_'),
(f.velocity, 'Axial Velocity (m/s)', 'axial_velocity_plot_at_'),
(f.L, 'Lambda', 'lambda_plot_at_'),
(f.Uo, 'Oxidizer Velocity Uo (m/s)', 'oxidizer_velocity_plot_at_')
]
for data, ylabel, prefix in variables:
plot_variable(f.grid, data, f'State at t={t}', 'Grid Points', ylabel, f'{prefix}{t}.png')
#input("Press Enter to continue...") # Pause for user input before continuing
return 0
def unsteady_callback(t, f):
variables = [
(f.T, 'Temperature (K)', 'temperature_plot_unsteady'),
(f.velocity, 'Axial Velocity (m/s)', 'axial_velocity_plot_unsteady'),
(f.L, 'Lambda', 'lambda_plot_unsteady'),
(f.Uo, 'Oxidizer Velocity Uo (m/s)', 'oxidizer_velocity_plot_unsteady')
]
for data, ylabel, prefix in variables:
plot_variable(f.grid, data, f'State at t={t}', 'Grid Points', ylabel, f'{prefix}.png')
#input("Press Enter to continue...") # Pause for user input before continuing
return 0
def plot_process(process_data):
#Plot process
# 1/a is often used as the parameter to plot flame solution data against in literature
inverse_strain = [1.0/a for a in flamelet_process_plot['a_mean']]
fig, ax = plt.subplots()
ax.plot(inverse_strain, flamelet_process_plot['T_max'], 'o-')
plt.savefig('T_max_vs_strainRate.png')
fig, ax = plt.subplots()
ax.plot(inverse_strain, flamelet_process_plot['Q'], 'o-')
plt.savefig('Q_vs_strainRate.png')
fig, ax = plt.subplots()
ax.plot(np.diff(flamelet_process_plot['T_max']), 'o-')
plt.savefig('Temperature_change.png')
#-----------------------
def setup_flame(configs):
p = configs['p']
fuel_inlet_temp = configs['fuel_inlet_temp']
oxidizer_inlet_temp = configs['oxidizer_inlet_temp']
fuel_mdot = configs['fuel_mdot']
width = configs['width']
oxidizer_composition = configs['oxidizer_composition']
fuel_composition = configs['fuel_composition']
mechanism = configs['mechanism']
gas = ct.Solution(mechanism)
# Set the composition using mass fractions
gas.Y = oxidizer_composition
# Retrieve and print the mole fractions for the oxidizer
oxidizer_composition = gas.X
print("Oxidizer composition in mole fractions:", oxidizer_composition)
# Repeat for the fuel composition
gas.Y = fuel_composition
# Retrieve and print the mole fractions for the fuel
fuel_composition = gas.X
print("Fuel composition in mole fractions:", fuel_composition)
gas.TPX = fuel_inlet_temp, p, fuel_composition
density_f = gas.density
gas.TPX = oxidizer_inlet_temp, p, oxidizer_composition
density_o = gas.density
f = ct.CounterflowDiffusionFlame(gas, width=width)
# Set the state of the two inlets
f.fuel_inlet.mdot = fuel_mdot
f.fuel_inlet.X = fuel_composition
f.fuel_inlet.T = fuel_inlet_temp
#Create a guestimate for what the oxidizer mdot would be
f.oxidizer_inlet.mdot = (fuel_mdot / density_f) * density_o * 0.5
f.oxidizer_inlet.X = oxidizer_composition
f.oxidizer_inlet.T = oxidizer_inlet_temp
#Solver settings
f.set_refine_criteria(ratio=10, slope=0.6, curve=0.6, prune=0.05)
# Adjust solver tolerances
f.flame.set_transient_tolerances(default=(1e-8, 1e-12))
f.flame.set_steady_tolerances(default=(1e-8, 1e-12))
# Attach callback functions to the solver
if callable(configs['steady_callback']):
steady_callback = configs['steady_callback']
f.set_steady_callback(lambda t: steady_callback(t, f))
#if callable(configs['unsteady_callback']):
# unsteady_callback = configs['unsteady_callback']
# f.set_time_step_callback(lambda t: unsteady_callback(t, f))
return f
def generate_flamelets(configs):
print('Solving first solution with no two-point control active')
loglevel = configs['loglevel']
f = setup_flame(configs)
# Solve the initial problem with the standard counterflow solver settings
f.solve(loglevel)
print('Current Flamelet State:')
print('Flame: Tmin: ' + str(np.min(f.T)) + ', Tmax: ' + str(np.max(f.T)))
# Dict for storing useful data during execution
flamelet_process = {'T_max': [max(f.T)]}
flamelet_process['Q'] = [np.trapz(f.heat_release_rate, f.grid)]
flamelet_process['a_mean'] = [f.strain_rate('mean')]
#Point control
f.two_point_control_enabled = True
spacing = configs['spacing']
maximum_temperature_increment = configs['maximum_temperature_increment']
Tmax_threshold = 25
temperature_increment_growth_rate = configs['temperature_increment_growth_rate']
failed_previous_attempt = False
temperature_increment = configs['temperature_increment'] #Kelvin
flamelet_solution_num = 1
iterations = 0
create_output_folders()
while ((np.max(f.T)-np.min(f.T))/np.max(f.T)) > 0.005 and iterations < 1000:
control_temperature = np.min(f.T) + spacing*(np.max(f.T) - np.min(f.T))
print('New Control Temperature: ' + str(control_temperature))
f.set_left_control_point(control_temperature) #Sets a location based on the temperature of the current flame
f.set_right_control_point(control_temperature)
print(f'Current Flamelet {flamelet_solution_num} State:')
print('Flame: Tmin: ' + str(np.min(f.T)) + ', Tmax: ' + str(np.max(f.T)))
print('tLeft= ' + str(f.tLeft) + ' , tRight= ' + str(f.tRight))
print('ZLeft= ' + str(f.zLeft) + ' , ZRight= ' + str(f.zRight))
print('Solving flowfield without any temperature decrement to adjust for shifting temperature control points')
f.solve(loglevel)
if iterations == 0:
iterations +=1
continue # We want to run this first part onces using two-point algorithm with no temperature decrements before beginning to march
f.save(f'output/flamelets/flamelet_{flamelet_solution_num}.yaml',overwrite=True)
f.save('previous_solution.yaml',overwrite=True) #Backup in case of failure
# Lower the control point temperatures
f.tLeft -= temperature_increment
f.tRight -= temperature_increment
print(f'Attempting to solve flame with: TLeft = {f.tLeft}, TRight= {f.tRight}, Tincrement= {temperature_increment}')
iterations +=1
try:
f.solve(loglevel)
except Exception as err:
print(f'\nerror when using control points: TLeft = {f.tLeft}, TRight= {f.tRight}')
print(err)
# Reverse failed temperature control point decrement
f.tLeft += temperature_increment
f.tRight += temperature_increment
print(f'Resetting solution to previous state: TLeft = {f.tLeft}, TRight= {f.tRight}')
f.restore('previous_solution.yaml') #Backup in case of failure
control_temperature = np.min(f.T) + spacing*(np.max(f.T) - np.min(f.T))
print('Control Temperature: ' + str(control_temperature))
f.set_left_control_point(control_temperature)
f.set_right_control_point(control_temperature)
print(f'Attempting to solve flame with: TLeft = {f.tLeft}, TRight= {f.tRight}, Tincrement= {temperature_increment}')
f.solve(loglevel)
print('Successfully restored state & shifted control points' )
print('tLeft= ' + str(f.tLeft) + ' , tRight= ' + str(f.tRight))
print('ZLeft= ' + str(f.zLeft) + ' , ZRight= ' + str(f.zRight))
#Update increment
print('Reducing temperature increment to: ' + str(0.5*temperature_increment))
temperature_increment *= 0.5
if(temperature_increment < 5.0e-2):
print('Smallest temperature increment reached, exiting.')
break
failed_previous_attempt = True
continue
'''
if abs(np.max(f.T) - previous_Tmax) > Tmax_threshold: #Do not allow a large maximum temperature jump between solutions
print('Maximum temperature change violated, ajusting step size')
print('Reducing temperature increment to: ' + str(0.5*temperature_increment))
#Update increment
temperature_increment *= 0.5
if(temperature_increment < 5.0e-1):
spacing *= 1.01
print('No more solutions reached at temperature increment: ' + str(temperature_increment) + ', adjusting spacing to: ' + str(spacing))
temperature_increment = original_temperature_increment # Reset the increment
failed_previous_attempt = True
f.restore('previous_solution.yaml')
continue
'''
if failed_previous_attempt == False:
if temperature_increment < maximum_temperature_increment:
temperature_increment *= temperature_increment_growth_rate
failed_previous_attempt = False
flamelet_process['T_max'].append(max(f.T))
flamelet_process['Q'].append(np.trapz(f.heat_release_rate, f.grid))
flamelet_process['a_mean'].append(f.strain_rate('mean'))
print(f"Flame Information: {flamelet_process['T_max'][-1]:.4g} K (a = {flamelet_process['a_mean'][-1]:.4g}, Q = {flamelet_process['Q'][-1]:.4g})")
plot_solution(f, flamelet_solution_num, control_points=True)
plot_mixture_fraction_solution(f, flamelet_solution_num, control_points=True)
flamelet_solution_num += 1
if iterations == 1000:
print('Exited due to maximum iterations being exceeded')
if __name__ == '__main__':
configs = {}
#'''
# Input parameters
configs['p'] = 8.0e5 # pressure [Pascals]
configs['fuel_inlet_temp'] = 1200.0 # fuel inlet temperature [Kelvin]
configs['oxidizer_inlet_temp'] = 1200.0 # oxidizer inlet temperature [Kelvin]
configs['fuel_mdot'] = 0.2 # kg/m^2/s
configs['width'] = 20e-3 # Distance between inlets [m]
configs['loglevel'] = 1 # amount of diagnostic output (0 to 9)
#Mass Fraction composition
configs['oxidizer_composition'] = 'O2:0.8, N2:0.2' # oxidizer composition (right)
configs['fuel_composition'] = 'H2:0.9, N2: 0.1' # fuel composition (left)
# Reaction mechanism
configs['mechanism'] = 'mechanism/burke/burke.yaml'
# Two-Point Controls
configs['spacing'] = 0.95
configs['temperature_increment'] = 10
configs['maximum_temperature_increment'] = 25
configs['temperature_increment_growth_rate'] = 1.03 #Give this as a factor where 1.02 would be a 2% increase per iteration
#Callbacks if they need to be used
configs['steady_callback'] = steady_callback
configs['unsteady_callback'] = unsteady_callback
#'''
'''
# Ray Speths' Test
# Input parameters
configs['p'] = ct.one_atm # pressure [Pascals]
configs['fuel_inlet_temp'] = 300.0 # fuel inlet temperature [Kelvin]
configs['oxidizer_inlet_temp'] = 300.0 # oxidizer inlet temperature [Kelvin]
configs['fuel_mdot'] = 0.2 # kg/m^2/s
configs['width'] = 50e-3 # Distance between inlets [m]
configs['loglevel'] = 1 # amount of diagnostic output (0 to 9)
#Mass Fraction composition
configs['oxidizer_composition'] = 'O2:0.21, AR:0.78' # oxidizer composition (right)
configs['fuel_composition'] = 'H2:1, AR:1' # fuel composition (left)
# Reaction mechanism
configs['mechanism'] = 'h2o2.yaml'
# Two-Point Controls
configs['spacing'] = 0.95
configs['temperature_increment'] = 10
configs['maximum_temperature_increment'] = 25
#Callbacks if they need to be used
configs['steady_callback'] = steady_callback
configs['unsteady_callback'] = unsteady_callback
'''
'''
# Input parameters
configs['p'] = 5.4e6 # pressure [Pascals]
configs['fuel_inlet_temp'] = 800.0 # fuel inlet temperature [Kelvin]
configs['oxidizer_inlet_temp'] = 711.0 # oxidizer inlet temperature [Kelvin]
configs['fuel_mdot'] = 0.13 # kg/m^2/s
configs['width'] = 60e-3 # Distance between inlets [m]
configs['oxidizer_composition'] = 'O2:0.9449, H2O:0.0551' # oxidizer composition (right)
configs['fuel_composition'] = 'H2:0.4018, H2O: 0.5982' # fuel composition (left)
configs['loglevel'] = 1 # amount of diagnostic output (0 to 5)
configs['mechanism'] = 'mechanism/H2.yaml'
configs['steady_callback'] = steady_callback
'''
generate_flamelets(configs)
ischoegl
speth
codecov[bot]
This pull request supersedes #1510 and is based on the newer coding in the StFlow class. This feature aims to provide a capability for users to generate 1D flame solutions that include solutions along the unstable burning branch of the typical s-curve that is used for flamelet progress variable methods.
Checklist
scons build&scons test) and unit tests address code coverage