improved comments for control point picking, added more tests, updated example script

samples/python/onedim/diffusion_flame_continuation.py

@@ -6,17 +6,13 @@ Diffusion flame two-point control
 =================================
 
 This example uses the two-point flame control feature to march solutions
-down the unstable burning branch for a counterflow diffusion flame.
+down the stable and unstable burning branch for a counterflow diffusion flame.
 A hydrogen-oxygen diffusion flame at 1 bar is studied.
 
-The tutorial makes use of the scaling rules derived by Fiala and Sattelmayer
-(doi:10.1155/2014/484372). Please refer to this publication for a detailed
-explanation. Also, please don't forget to cite it if you make use of it.
-
 Requires: cantera >= 3.0, matplotlib >= 2.0
 
 .. tags:: Python, combustion, 1D flow, diffusion flame, strained flame, extinction,
-          saving output, plotting
+          saving output, plotting, two-point control
 """
 
 from pathlib import Path
@@ -78,137 +74,21 @@ file_name, entry = names("initial-solution")
 f.save(file_name, name=entry, description="Initial solution")
 
 
-# PART 2: COMPUTE EXTINCTION STRAIN
+# PART 2: CONTINUATION
 
-# Exponents for the initial solution variation with changes in strain rate
-# Taken from Fiala and Sattelmayer (2014)
-exp_d_a = - 1. / 2.
-exp_u_a = 1. / 2.
-exp_V_a = 1.
-exp_lam_a = 2.
-exp_mdot_a = 1. / 2.
-
-# Set normalized initial strain rate
-alpha = [1.]
-# Initial relative strain rate increase
-delta_alpha = 1.
-# Factor of refinement of the strain rate increase
-delta_alpha_factor = 50.
-# Limit of the refinement: Minimum normalized strain rate increase
-delta_alpha_min = .001
-# Limit of the Temperature decrease
-delta_T_min = 1  # K
-
-# Iteration indicator
-n = 0
-# Indicator of the latest flame still burning
-n_last_burning = 0
-# List of peak temperatures
-T_max = [np.max(f.T)]
-# List of maximum axial velocity gradients
-a_max = [np.max(np.abs(np.gradient(f.velocity) / np.gradient(f.grid)))]
-
-# Simulate counterflow flames at increasing strain rates until the flame is
-# extinguished. To achieve a fast simulation, an initial coarse strain rate
-# increase is set. This increase is reduced after an extinction event and
-# the simulation is again started based on the last burning solution.
-# The extinction point is considered to be reached if the abortion criteria
-# on strain rate increase and peak temperature decrease are fulfilled.
-while True:
-    n += 1
-    # Update relative strain rates
-    alpha.append(alpha[n_last_burning] + delta_alpha)
-    strain_factor = alpha[-1] / alpha[n_last_burning]
-    # Create an initial guess based on the previous solution
-    # Update grid
-    # Note that grid scaling changes the diffusion flame width
-    f.flame.grid *= strain_factor ** exp_d_a
-    normalized_grid = f.grid / (f.grid[-1] - f.grid[0])
-    # Update mass fluxes
-    f.fuel_inlet.mdot *= strain_factor ** exp_mdot_a
-    f.oxidizer_inlet.mdot *= strain_factor ** exp_mdot_a
-    # Update velocities
-    f.set_profile('velocity', normalized_grid,
-                  f.velocity * strain_factor ** exp_u_a)
-    f.set_profile('spread_rate', normalized_grid,
-                  f.spread_rate * strain_factor ** exp_V_a)
-    # Update pressure curvature
-    f.set_profile('lambda', normalized_grid, f.L * strain_factor ** exp_lam_a)
-    try:
-        f.solve(loglevel=0)
-    except ct.CanteraError as e:
-        print('Error: Did not converge at n =', n, e)
-
-    T_max.append(np.max(f.T))
-    a_max.append(np.max(np.abs(np.gradient(f.velocity) / np.gradient(f.grid))))
-    if not np.isclose(np.max(f.T), temperature_limit_extinction):
-        # Flame is still burning, so proceed to next strain rate
-        n_last_burning = n
-        file_name, entry = names(f"extinction/{n:04d}")
-        f.save(file_name, name=entry, description=f"Solution at alpha = {alpha[-1]}")
-
-        print('Flame burning at alpha = {:8.4F}. Proceeding to the next iteration, '
-              'with delta_alpha = {}'.format(alpha[-1], delta_alpha))
-    elif ((T_max[-2] - T_max[-1] < delta_T_min) and (delta_alpha < delta_alpha_min)):
-        # If the temperature difference is too small and the minimum relative
-        # strain rate increase is reached, save the last, non-burning, solution
-        # to the output file and break the loop
-        file_name, entry = names(f"extinction/{n:04d}")
-        f.save(file_name, name=entry, description=f"Flame extinguished at alpha={alpha[-1]}")
-
-        print('Flame extinguished at alpha = {0:8.4F}.'.format(alpha[-1]),
-              'Abortion criterion satisfied.')
-        break
-    else:
-        # Procedure if flame extinguished but abortion criterion is not satisfied
-        # Reduce relative strain rate increase
-        delta_alpha = delta_alpha / delta_alpha_factor
-
-        print('Flame extinguished at alpha = {0:8.4F}. Restoring alpha = {1:8.4F} and '
-              'trying delta_alpha = {2}'.format(
-                  alpha[-1], alpha[n_last_burning], delta_alpha))
-
-        # Restore last burning solution
-        file_name, entry = names(f"extinction/{n_last_burning:04d}")
-        f.restore(file_name, entry)
-
-
-# Print some parameters at the extinction point, after restoring the last burning
-# solution
-file_name, entry = names(f"extinction/{n_last_burning:04d}")
-f.restore(file_name, entry)
-
-print('----------------------------------------------------------------------')
-print('Parameters at the extinction point:')
-print('Pressure p={0} bar'.format(f.P / 1e5))
-print('Peak temperature T={0:4.0f} K'.format(np.max(f.T)))
-print('Mean axial strain rate a_mean={0:.2e} 1/s'.format(f.strain_rate('mean')))
-print('Maximum axial strain rate a_max={0:.2e} 1/s'.format(f.strain_rate('max')))
-print('Fuel inlet potential flow axial strain rate a_fuel={0:.2e} 1/s'.format(
-      f.strain_rate('potential_flow_fuel')))
-print('Oxidizer inlet potential flow axial strain rate a_ox={0:.2e} 1/s'.format(
-      f.strain_rate('potential_flow_oxidizer')))
-print('Axial strain rate at stoichiometric surface a_stoich={0:.2e} 1/s'.format(
-      f.strain_rate('stoichiometric', fuel='H2')))
-
-# Plot the maximum temperature over the maximum axial velocity gradient
-plt.figure()
-plt.semilogx(a_max, T_max)
-plt.xlabel(r'$a_{max}$ [1/s]')
-plt.ylabel(r'$T_{max}$ [K]')
-plt.savefig(output_path / "figure_T_max_a_max.png")
-plt.close()
-
-# Now restart and activate the two-point flame control feature
 f.two_point_control_enabled = True
 spacing = 0.95
-temperature_increment = 0.1 # Kelvin
+temperature_increment = 1.0 # Kelvin
 maximum_temperature = []
-for i in range(500):
+a_max = []
+for i in range(5000):
     control_temperature = np.min(f.T) + spacing*(np.max(f.T) - np.min(f.T))
+    print(f'Iteration {i}: Control temperature = {control_temperature} 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
 
@@ -219,8 +99,17 @@ for i in range(500):
         break
 
     maximum_temperature.append(np.max(f.T))
+    a_max.append(np.max(np.abs(np.gradient(f.velocity) / np.gradient(f.grid))))
 
 
+# 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")
+plt.close()
+
 # Plot maximum_temperature against number of iterations
 plt.figure()
 plt.plot(range(len(maximum_temperature)), maximum_temperature)
@@ -229,3 +118,5 @@ plt.ylabel('Maximum Temperature [K]')
 plt.savefig(output_path / "figure_max_temperature_iterations.png")
 plt.close()
 
+
+