1D flame negative left boundary during 2-point control

[wandadars]

In the following script, during iteration 15, the solve() fails. During the failure, the script loops back to the top and tries again. Upon printing the data for the flame at the start of the loop, it appears that the previously successful flame state was not preserved.

On iteration 18 there doesn't appear to be any issues, but the final solution seems to converge to a negative velocity at the right boundary. On the next iteration it recovers. I'm using the python interface on the main branch of Cantera.

I'm trying to figure out if I've done something obviously wrong here, but I just wanted to report it for now.

from pathlib import Path
import cantera as ct
import numpy as np
import matplotlib.pyplot as plt
import imageio
from pathlib import Path


p = 5.4e6 # pressure [Pascals]

fuel_inlet_temp= 800.0  # fuel inlet temperature [Kelvin]
oxidizer_inlet_temp= 711.0  # oxidizer inlet temperature [Kelvin]
fuel_mdot= 0.5 # kg/m^2/s
width= 50.0e-3 # Distance between inlets [m]

oxidizer_composition= 'O2:1.0'  # oxidizer composition (right)
fuel_composition= 'H2:1.0'  # fuel composition (left)

# Reaction mechanism
mechanism= 'h2o2.yaml'

# Define output locations
output_path = Path() / "testing_output"
output_path.mkdir(parents=True, exist_ok=True)

gas = ct.Solution(mechanism)

gas.TPY = fuel_inlet_temp, p, oxidizer_composition
density_f = gas.density

gas.TPY = oxidizer_inlet_temp, p, oxidizer_composition
density_o = gas.density

#Unity Lewis number testing
gas.transport_model = 'unity-Lewis-number'

f = ct.CounterflowDiffusionFlame(gas, width=width)

f.set_refine_criteria(ratio=15, slope= 0.15, curve= 0.08, prune= 0.04)
#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*4
f.oxidizer_inlet.X = oxidizer_composition
f.oxidizer_inlet.T = oxidizer_inlet_temp

# Generate initial condition
f.solve(auto=True, loglevel=4)
f.save('backup.yaml', name="solution", overwrite=True)

print('mdot info:')
print('Fuel mdot: ' + str(f.fuel_inlet.mdot))
print('Oxidizer mdot: ' + str(f.oxidizer_inlet.mdot))
print('BC State:')
print('Left U: ' + str(f.velocity[0]))
print('Right U: ' + str(f.velocity[-1]))


print('Starting two-point control')
f.two_point_control_enabled = True

spacing = 0.75
temperature_increment = 10 # Initial temperature increment
maximum_temperature = []
a_max = []
filenames = []
for i in range(80):
    print('Iteration: ' + str(i))
    print('mdot info:')
    print('Fuel mdot: ' + str(f.fuel_inlet.mdot))
    print('Oxidizer mdot: ' + str(f.oxidizer_inlet.mdot))

    print('BC State:')
    print('Left U: ' + str(f.velocity[0]))
    print('Right U: ' + str(f.velocity[-1]))

    control_temperature = np.min(f.T) + spacing*(np.max(f.T) - np.min(f.T))
    print(f'Iteration {i}: Control temperature = {control_temperature} K')
    print(f'Increment size: {temperature_increment} K')
    f.set_left_control_point(control_temperature)
    f.set_right_control_point(control_temperature)
    # This decrement is what drives the two-point control. If failure
    # occurs, try decreasing the decrement.
    f.left_control_point_temperature -= temperature_increment
    f.right_control_point_temperature -= temperature_increment

    try:
        f.solve(loglevel=1)
        #f.solve(loglevel=0,refine_grid=False)
    except ct.CanteraError as e:
        print('Solver did not converge.')
        print('Error:')
        print(e)

        if temperature_increment < 1:
            print('Smallest increment reached, stopping')
            break

        f.left_control_point_temperature += temperature_increment
        f.right_control_point_temperature += temperature_increment
        temperature_increment *= 0.5
        continue

    # Increase the increment if the solution worked out well
    if temperature_increment < 50:
        temperature_increment *= 1.1

    maximum_temperature.append(np.max(f.T))
    a_max.append(np.max(np.abs(np.gradient(f.velocity) / np.gradient(f.grid))))

    print(f'Plotting data for iteration {i}')
    # Plot the heat release rate as well as the locations of the control points using the
    # f.left_control_point_coordinate and f.right_control_point_coordinate attributes. Plot
    # The locations as vertical lines. Generate a gif of the heat release rate as a function
    # of position with the markers and save it as heat_release_rate.gif
    fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(16, 8))  # 16x8 is a fairly large and wide size
    # Set the title as having the iteration number in it
    fig.suptitle(f"Iteration {i}")

    # Plot the U
    ax1.plot(f.grid, f.velocity,marker='o', label='Velocity')
    ax1.axvline(f.left_control_point_coordinate, color='r', linestyle='--')
    ax1.axvline(f.right_control_point_coordinate, color='r', linestyle='--')
    ax1.set_xlabel('Position [m]')
    ax1.set_ylabel('Axial Velocity [m/s]')
    ax1.legend()

    # Plot the V
    ax2.plot(f.grid, f.spread_rate,marker='o', label='Spread Rate')
    ax2.axvline(f.left_control_point_coordinate, color='r', linestyle='--')
    ax2.axvline(f.right_control_point_coordinate, color='r', linestyle='--')
    ax2.set_xlabel('Position [m]')
    ax2.set_ylabel('Spread Rate [1/s]')
    ax2.legend()

    # Plot the temperature
    ax3.plot(f.grid, f.T,marker='o', label='Temperature')
    ax3.axvline(f.left_control_point_coordinate, color='r', linestyle='--')
    ax3.axvline(f.right_control_point_coordinate, color='r', linestyle='--')
    ax3.set_xlabel('Position [m]')
    ax3.set_ylabel('Temperature [K]')
    ax3.legend()

    # Constrain the x axis to be between 0.02 and 0.03
    #ax1.set_xlim(0.02, 0.03)
    #ax2.set_xlim(0.02, 0.03)
    #ax3.set_xlim(0.02, 0.03)


    # Save each frame as a separate image file
    frame_filename = output_path / f"flame_{i}.png"
    plt.savefig(frame_filename)
    plt.close(fig)  # Close the figure to avoid memory leaks

    # Append the filename to the list
    filenames.append(frame_filename)

# Generate the gif(duration is in hundreths of a second per frame)
with imageio.get_writer(output_path / "heat_release_rate.gif", mode='I', duration=500, loop=0) as writer:
    for filename in filenames:
        image = imageio.imread(filename)
        writer.append_data(image)

# Plot the maximum temperature versus the maximum axial velocity gradient
plt.figure()
plt.plot(a_max, maximum_temperature)
plt.xlabel('Maximum Axial Velocity Gradient [1/s]')
plt.ylabel('Maximum Temperature [K]')
plt.savefig(output_path / "figure_max_temperature_vs_max_velocity_gradient.png")


# Plot maximum_temperature against number of iterations
plt.figure()
plt.plot(range(len(maximum_temperature)), maximum_temperature)
plt.xlabel('Number of Continuation Steps')
plt.ylabel('Maximum Temperature [K]')
plt.savefig(output_path / "figure_max_temperature_iterations.png")

[speth]

I had seen a couple instances of this behavior as well. One way to fix this is to require the spread_rate to be positive by setting

f.flame.set_bounds(spread_rate=(-1e-5, 1e20))

I had added this to the diffusion_flame_continuation.py example, but forgot to push that change before #1779 was merged.