From 9fd067dbaae76ddedc738fd787ce2f96b74be802 Mon Sep 17 00:00:00 2001 From: Christopher Neal Date: Wed, 19 Jun 2024 15:52:17 -0400 Subject: [PATCH] formatting changes & documentation/samples tweaks --- .../reference/onedim/governing-equations.md | 66 +++++----- include/cantera/oneD/StFlow.h | 119 +++--------------- interfaces/cython/cantera/_onedim.pxd | 2 +- .../onedim/diffusion_flame_continuation.py | 48 ++----- src/oneD/Sim1D.cpp | 53 ++++---- src/oneD/StFlow.cpp | 115 ++++++++++++++++- test/python/test_onedim.py | 6 +- 7 files changed, 211 insertions(+), 198 deletions(-) diff --git a/doc/sphinx/reference/onedim/governing-equations.md b/doc/sphinx/reference/onedim/governing-equations.md index cb9833b95..645d58d8f 100644 --- a/doc/sphinx/reference/onedim/governing-equations.md +++ b/doc/sphinx/reference/onedim/governing-equations.md @@ -224,15 +224,17 @@ number. ## Counterflow Two-Point Flame Control -A two-point temperature control feature is available for counterflow diffusion flames. This -feature allows users to set a control points on both sides of a flame and incrementally lower -the flame temperature. This allows for the simulation of the stable burning branch as well -as the unstable burning branch of the standard flamelet "S-curve". The implementation is based -on the method discussed in {cite:t}`nishioka1996` and {cite:t}`huo2014`. The diagram below shows -the general concept of the two-point flame control method, with control points located on either -side of the peak flame temperature. An initial flame solution is used as a starting point, and -the temperatures at the control points are lowered to produce a new flame solution that satisfies -the governing equations and passes through the new temperatures at the control points. +A two-point temperature control feature is available for counterflow diffusion flames. +This feature allows users to set a control point on each side of a flame and +incrementally lower the flame temperature. This allows for the simulation of the stable +burning branch as well as the unstable burning branch of the standard flamelet +"S-curve". The implementation is based on the method discussed in +{cite:t}`nishioka1996` and {cite:t}`huo2014`. The diagram below shows the general +concept of the two-point flame control method, with control points located on either +side of the peak flame temperature. An initial flame solution is used as a starting +point, and the temperatures at the control points are lowered to produce a new flame +solution that satisfies the governing equations and passes through the new temperatures +at the control points. ```{image} /_static/images/two_point_control_diagram.svg :width: 50% @@ -240,12 +242,13 @@ the governing equations and passes through the new temperatures at the control p :align: center ``` -For the two-point control method, one governing equation was modified ($\Lambda$), and a new -governing equation for the axial oxidizer velocity was added ($U_o$). The fuel and oxidizer velocity boundary -conditions are modified when the two-point control is active. These equations allow -for the temperature reduction to be performed in a numerically consistent manner (preventing -any issues of over-defining the system of governing equations). The two equations that are activated -when two-point control is turned on are: +For the two-point control method, one governing equation is modified ($\Lambda$), and +a new governing equation for the axial oxidizer velocity is added ($U_o$). The fuel +and oxidizer velocity boundary conditions are modified when the two-point control is +active. These equations allow for the temperature reduction to be performed in a +numerically consistent manner (preventing any issues of over-defining the system of +governing equations). The two equations that are activated when two-point control is +turned on are: $$ \frac{d\Lambda}{dz} = 0 @@ -257,40 +260,41 @@ $$ \frac{dU_o}{dz} = 0 $$ -These equations are zero everywhere, except at their respective control points. At the left control point -the residual for the $\Lambda$ equation is: +At the left control point the residual for the $\Lambda$ equation is: $$ -residual = T(z=Z_L) - T_{L, control} +T(z=Z_L) - T_{L, \t{control}} $$ At the left control point the residual for the $U_o$ equation is: $$ -residual = T(z=Z_R) - T_{R, control} +T(z=Z_R) - T_{R, \t{control}} $$ -Where T(z=Z_L) is the temperature of the flowfield at the left control point, T(z=Z_R) is the temperature of -the flowfield at the right control point, T_{L, control} is the left control point desired temperature, and -T_{R, control} is the right control point desired temperature. +Where $T(z=Z_L)$ is the temperature of the flowfield at the left control point, +$T(z=Z_R)$ is the temperature of the flowfield at the right control point, +$T_{L, \t{control}}$ is the left control point desired temperature, and +$T_{R, \t{control}}$ is the right control point desired temperature. -The values of $\Lambda$ and $U_o$ are influenced by the left and right control points, respectively. -A residual error is induced because of the difference between the flow's temperature at that point and -the desired control point temperature. -In order to drive this error to zero, the solver adjusts the flow rates at the boundaries, which changes +The values of $\Lambda$ and $U_o$ are influenced by the left and right control points, +respectively. A residual error is induced because of the difference between the flow's +temperature at that point and the desired control point temperature. In order to drive +this error to zero, the solver adjusts the flow rates at the boundaries, which changes the temperature distribution, which in turn affects the values of $\Lambda$ and $U_o$. -At the right boundary, the boundary condition for the continuity equation is imposed by using -the solution from the oxidizer velocity equation. At the left boundary, the boundary condition for -the continuity equation is imposed by using the value of the axial velocity at the left boundary. +At the right boundary, the boundary condition for the continuity equation is imposed +by using the solution from the oxidizer velocity equation. At the left boundary, the +boundary condition for the continuity equation is imposed by using the value of the +axial velocity at the left boundary. $$ -mdot(z=0) = rho(z=0)*U(z=0) +\dot{m}_\t{in,left} = \rho(z_{\t{in,left}}) U(z_{\t{in,left}}) $$ $$ -mdot(z=L) = rho(z=L)*U_0(z=L) +\dot{m}_\t{in,right} = \rho(z_{\t{in,right}}) U_0(z_{\t{in,right}}) $$ diff --git a/include/cantera/oneD/StFlow.h b/include/cantera/oneD/StFlow.h index 57b01bfcb..f8bf85911 100644 --- a/include/cantera/oneD/StFlow.h +++ b/include/cantera/oneD/StFlow.h @@ -248,10 +248,14 @@ public: void setBoundaryEmissivities(double e_left, double e_right); //! Return emissivity at left boundary - double leftEmissivity() const { return m_epsilon_left; } + double leftEmissivity() const { + return m_epsilon_left; + } //! Return emissivity at right boundary - double rightEmissivity() const { return m_epsilon_right; } + double rightEmissivity() const { + return m_epsilon_right; + } void fixTemperature(size_t j=npos); @@ -278,124 +282,31 @@ public: //! @{ //! Returns the temperature at the left control point - double leftControlPointTemperature() const { - if (m_twoPointControl) { - if (m_zLeft != Undef){ - return m_tLeft; - } else { - throw CanteraError("StFlow::leftControlPointTemperature", - "Invalid operation: left control point location is not set"); - } - } else { - throw CanteraError("StFlow::leftControlPointTemperature", - "Invalid operation: two-point control is not enabled."); - } - } + double leftControlPointTemperature() const; //! Returns the z-coordinate of the left control point - double leftControlPointCoordinate() const { - if (m_twoPointControl) { - if (m_zLeft != Undef) - return m_zLeft; - else { - throw CanteraError("StFlow::leftControlPointCoordinate", - "Invalid operation: left control point location is not set"); - } - } else { - throw CanteraError("StFlow::leftControlPointCoordinate", - "Invalid operation: two-point control is not enabled."); - } - } + double leftControlPointCoordinate() const; //! Sets the temperature of the left control point - void setLeftControlPointTemperature(double temperature) { - if (m_twoPointControl) { - if (m_zLeft != Undef){ - m_tLeft = temperature; - } else { - throw CanteraError("StFlow::setLeftControlPointTemperature", - "Invalid operation: left control point location is not set"); - } - } else { - throw CanteraError("StFlow::setLeftControlPointTemperature", - "Invalid operation: two-point control is not enabled."); - } - } + void setLeftControlPointTemperature(double temperature); //! Sets the coordinate of the left control point - void setLeftControlPointCoordinate(double z_left) { - if (m_twoPointControl) { - m_zLeft = z_left; - } else { - throw CanteraError("StFlow::setLeftControlPointCoordinate", - "Invalid operation: two-point control is not enabled."); - } - } + void setLeftControlPointCoordinate(double z_left); //! Returns the temperature at the right control point - double rightControlPointTemperature() const { - if (m_twoPointControl) { - if (m_zRight != Undef) { - return m_tRight; - } else { - throw CanteraError("StFlow::rightControlPointTemperature", - "Invalid operation: right control point location is not set"); - } - } else { - throw CanteraError("StFlow::rightControlPointTemperature", - "Invalid operation: two-point control is not enabled."); - } - } + double rightControlPointTemperature() const; //! Returns the z-coordinate of the right control point - double rightControlPointCoordinate() const { - if (m_twoPointControl) { - if (m_zRight != Undef){ - return m_zRight; - } else { - throw CanteraError("StFlow::rightControlPointCoordinate", - "Invalid operation: right control point location is not set"); - } - } else { - throw CanteraError("StFlow::rightControlPointCoordinate", - "Invalid operation: two-point control is not enabled."); - } - } + double rightControlPointCoordinate() const; //! Sets the temperature of the right control point - void setRightControlPointTemperature(double temperature) { - if (m_twoPointControl) { - if (m_zRight != Undef){ - m_tRight = temperature; - } else { - throw CanteraError("StFlow::setRightControlPointTemperature", - "Invalid operation: right control point location is not set"); - } - } else { - throw CanteraError("StFlow::setRightControlPointTemperature", - "Invalid operation: two-point control is not enabled."); - } - } + void setRightControlPointTemperature(double temperature); //! Sets the coordinate of the right control point - void setRightControlPointCoordinate(double z_right) { - if (m_twoPointControl) { - m_zRight = z_right; - } else { - throw CanteraError("StFlow::setRightControlPointCoordinate", - "Invalid operation: two-point control is not enabled."); - } - } + void setRightControlPointCoordinate(double z_right); //! Sets the status of the two-point control - void enableTwoPointControl(bool twoPointControl) { - if (isStrained()){ - m_twoPointControl = twoPointControl; - } else { - throw CanteraError("StFlow::enableTwoPointControl", - "Invalid operation: two-point control can only be used with axisymmetric flames."); - } - } + void enableTwoPointControl(bool twoPointControl); //! Returns the status of the two-point control bool twoPointControlEnabled() const { diff --git a/interfaces/cython/cantera/_onedim.pxd b/interfaces/cython/cantera/_onedim.pxd index 10284fb8e..cf2800bf3 100644 --- a/interfaces/cython/cantera/_onedim.pxd +++ b/interfaces/cython/cantera/_onedim.pxd @@ -92,7 +92,7 @@ cdef extern from "cantera/oneD/StFlow.h": ThermoBasis fluxGradientBasis() void setFreeFlow() void setAxisymmetricFlow() - void enableTwoPointControl(cbool) + void enableTwoPointControl(cbool) except +translate_exception cbool twoPointControlEnabled() double leftControlPointTemperature() except +translate_exception double leftControlPointCoordinate() except +translate_exception diff --git a/samples/python/onedim/diffusion_flame_continuation.py b/samples/python/onedim/diffusion_flame_continuation.py index fa47effdd..6aec5e830 100644 --- a/samples/python/onedim/diffusion_flame_continuation.py +++ b/samples/python/onedim/diffusion_flame_continuation.py @@ -1,18 +1,15 @@ -# This file is part of Cantera. See License.txt in the top-level directory or -# at https://cantera.org/license.txt for license and copyright information. - """ -Diffusion flame two-point control -================================= +Diffusion flame unstable branch +=============================== This example uses the two-point flame control feature to march solutions down the stable and unstable burning branch for a counterflow diffusion flame. A hydrogen-oxygen diffusion flame at 1 bar is studied. -Requires: cantera >= 3.0, matplotlib >= 2.0 +Requires: cantera >= 3.1, matplotlib >= 2.0 .. tags:: Python, combustion, 1D flow, diffusion flame, strained flame, extinction, - saving output, plotting, two-point control + saving output, plotting """ from pathlib import Path @@ -21,8 +18,9 @@ import matplotlib.pyplot as plt import cantera as ct - -# PART 1: INITIALIZATION +# %% +# Flame Initialization +# -------------------- # Set up an initial hydrogen-oxygen counterflow flame at 1 bar and low strain # rate (maximum axial velocity gradient = 2414 1/s) @@ -44,10 +42,6 @@ f.oxidizer_inlet.T = 500 # K # Set refinement parameters f.set_refine_criteria(ratio=3.0, slope=0.1, curve=0.2, prune=0.03) -# Define a limit for the maximum temperature below which the flame is -# considered as extinguished and the computation is aborted -temperature_limit_extinction = max(f.oxidizer_inlet.T, f.fuel_inlet.T) - # Initialize and solve print('Creating the initial solution') f.solve(loglevel=0, auto=True) @@ -56,29 +50,13 @@ f.solve(loglevel=0, auto=True) output_path = Path() / "diffusion_flame_continuation_data" output_path.mkdir(parents=True, exist_ok=True) -hdf_output = "native" in ct.hdf_support() -if hdf_output: - file_name = output_path / "flame_data.h5" - file_name.unlink(missing_ok=True) - -def names(test): - if hdf_output: - # use internal container structure for HDF - file_name = output_path / "flame_data.h5" - return file_name, test - # use separate files for YAML - file_name = output_path / f"{test}.yaml".replace("-", "_").replace("/", "_") - return file_name, "solution" - -file_name, entry = names("initial-solution") -f.save(file_name, name=entry, description="Initial solution") - - -# PART 2: CONTINUATION +# %% +# Flame Continuation +# ------------------ f.two_point_control_enabled = True spacing = 0.95 -temperature_increment = 1.0 # Kelvin +temperature_increment = 10.0 # Kelvin maximum_temperature = [] a_max = [] for i in range(5000): @@ -108,7 +86,7 @@ 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() @@ -116,7 +94,7 @@ 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") -plt.close() + diff --git a/src/oneD/Sim1D.cpp b/src/oneD/Sim1D.cpp index c1963354b..db3b085e4 100644 --- a/src/oneD/Sim1D.cpp +++ b/src/oneD/Sim1D.cpp @@ -742,19 +742,24 @@ void Sim1D::setLeftControlPoint(double temperature) if (!d_axis.twoPointControlEnabled()) { continue; } - two_point_domain_found = true; // At least one domain with two-point control enabled was found + two_point_domain_found = true; + double current_val, next_val; for (size_t m = 0; m < np-1; m++) { - if ((value(n,c_offset_T,m) - temperature) * (value(n,c_offset_T,m+1) - temperature) < 0.0) { - // Pick the coordinate of the point with the temperature closest to the desired temperature - int closest_index = 0; - if (std::abs(value(n,c_offset_T,m) - temperature) < std::abs(value(n,c_offset_T,m+1) - temperature)) { - closest_index = m; + current_val = value(n,c_offset_T,m); + next_val = value(n,c_offset_T,m+1); + if ((current_val - temperature) * (next_val - temperature) < 0.0) { + // Pick the coordinate of the point with the temperature closest + // to the desired temperature + int index = 0; + if (std::abs(current_val - temperature) < + std::abs(next_val - temperature)) { + index = m; } else { - closest_index = m+1; + index = m+1; } - d_axis.setLeftControlPointCoordinate(d_axis.grid(closest_index)); - d_axis.setLeftControlPointTemperature(value(n,c_offset_T,closest_index)); + d_axis.setLeftControlPointCoordinate(d_axis.grid(index)); + d_axis.setLeftControlPointTemperature(value(n,c_offset_T,index)); return; } } @@ -765,7 +770,8 @@ void Sim1D::setLeftControlPoint(double temperature) "No domain with two-point control enabled was found."); } else { throw CanteraError("Sim1D::setLeftControlPoint", - "No control point with temperature {} was able to be found in the flame's temperature range.", temperature); + "No control point with temperature {} was able to be found in the" + "flame's temperature range.", temperature); } } @@ -788,20 +794,24 @@ void Sim1D::setRightControlPoint(double temperature) if (!d_axis.twoPointControlEnabled()) { continue; } - two_point_domain_found = true; // At least one domain with two-point control enabled was found - + two_point_domain_found = true; + double current_val, next_val; for (size_t m = np-1; m > 0; m--) { - if ((value(n,c_offset_T,m) - temperature) * (value(n,c_offset_T,m-1) - temperature) < 0.0) { - // Pick the coordinate of the point with the temperature closest to the desired temperature - int closest_index = 0; - if (std::abs(value(n,c_offset_T,m) - temperature) < std::abs(value(n,c_offset_T,m-1) - temperature)) { - closest_index = m; + current_val = value(n,c_offset_T,m); + next_val = value(n,c_offset_T,m-1); + if ((current_val - temperature) * (next_val - temperature) < 0.0) { + // Pick the coordinate of the point with the temperature closest + // to the desired temperature + int index = 0; + if (std::abs(current_val - temperature) < + std::abs(next_val - temperature)) { + index = m; } else { - closest_index = m-1; + index = m-1; } - d_axis.setRightControlPointCoordinate(d_axis.grid(closest_index)); - d_axis.setRightControlPointTemperature(value(n,c_offset_T,closest_index)); + d_axis.setRightControlPointCoordinate(d_axis.grid(index)); + d_axis.setRightControlPointTemperature(value(n,c_offset_T,index)); return; } } @@ -812,7 +822,8 @@ void Sim1D::setRightControlPoint(double temperature) "No domain with two-point control enabled was found."); } else { throw CanteraError("Sim1D::setRightControlPoint", - "No control point with temperature {} was able to be found in the flame's temperature range.", temperature); + "No control point with temperature {} was able to be found in the" + "flame's temperature range.", temperature); } } diff --git a/src/oneD/StFlow.cpp b/src/oneD/StFlow.cpp index 28b211f5f..8ef4cd343 100644 --- a/src/oneD/StFlow.cpp +++ b/src/oneD/StFlow.cpp @@ -1044,7 +1044,7 @@ void StFlow::setMeta(const AnyMap& state) m_zRight = cm["right-location"].asDouble(); m_tLeft = cm["left-temperature"].asDouble(); m_tRight = cm["right-temperature"].asDouble(); - } else{ + } else { warn_user("StFlow::setMeta", "Unknown continuation method '{}'.", cm["type"].asString()); } @@ -1155,4 +1155,117 @@ void StFlow::grad_hk(const double* x, size_t j) } } +// Two-point control functions +double StFlow::leftControlPointTemperature() const { + if (m_twoPointControl) { + if (m_zLeft != Undef){ + return m_tLeft; + } else { + throw CanteraError("StFlow::leftControlPointTemperature", + "Invalid operation: left control point location is not set"); + } + } else { + throw CanteraError("StFlow::leftControlPointTemperature", + "Invalid operation: two-point control is not enabled."); + } +} + +double StFlow::leftControlPointCoordinate() const { + if (m_twoPointControl) { + if (m_zLeft != Undef) + return m_zLeft; + else { + throw CanteraError("StFlow::leftControlPointCoordinate", + "Invalid operation: left control point location is not set"); + } + } else { + throw CanteraError("StFlow::leftControlPointCoordinate", + "Invalid operation: two-point control is not enabled."); + } +} + +void StFlow::setLeftControlPointTemperature(double temperature) { + if (m_twoPointControl) { + if (m_zLeft != Undef){ + m_tLeft = temperature; + } else { + throw CanteraError("StFlow::setLeftControlPointTemperature", + "Invalid operation: left control point location is not set"); + } + } else { + throw CanteraError("StFlow::setLeftControlPointTemperature", + "Invalid operation: two-point control is not enabled."); + } +} + +void StFlow::setLeftControlPointCoordinate(double z_left) { + if (m_twoPointControl) { + m_zLeft = z_left; + } else { + throw CanteraError("StFlow::setLeftControlPointCoordinate", + "Invalid operation: two-point control is not enabled."); + } +} + +double StFlow::rightControlPointTemperature() const { + if (m_twoPointControl) { + if (m_zRight != Undef) { + return m_tRight; + } else { + throw CanteraError("StFlow::rightControlPointTemperature", + "Invalid operation: right control point location is not set"); + } + } else { + throw CanteraError("StFlow::rightControlPointTemperature", + "Invalid operation: two-point control is not enabled."); + } +} + +double StFlow::rightControlPointCoordinate() const { + if (m_twoPointControl) { + if (m_zRight != Undef){ + return m_zRight; + } else { + throw CanteraError("StFlow::rightControlPointCoordinate", + "Invalid operation: right control point location is not set"); + } + } else { + throw CanteraError("StFlow::rightControlPointCoordinate", + "Invalid operation: two-point control is not enabled."); + } +} + +void StFlow::setRightControlPointTemperature(double temperature) { + if (m_twoPointControl) { + if (m_zRight != Undef){ + m_tRight = temperature; + } else { + throw CanteraError("StFlow::setRightControlPointTemperature", + "Invalid operation: right control point location is not set"); + } + } else { + throw CanteraError("StFlow::setRightControlPointTemperature", + "Invalid operation: two-point control is not enabled."); + } +} + +void StFlow::setRightControlPointCoordinate(double z_right) { + if (m_twoPointControl) { + m_zRight = z_right; + } else { + throw CanteraError("StFlow::setRightControlPointCoordinate", + "Invalid operation: two-point control is not enabled."); + } +} + +void StFlow::enableTwoPointControl(bool twoPointControl) { + if (isStrained()){ + m_twoPointControl = twoPointControl; + } else { + throw CanteraError("StFlow::enableTwoPointControl", + "Invalid operation: two-point control can only be used" + "with axisymmetric flames."); + } +} + } // namespace diff --git a/test/python/test_onedim.py b/test/python/test_onedim.py index 3d8ce1c78..0e063abe8 100644 --- a/test/python/test_onedim.py +++ b/test/python/test_onedim.py @@ -155,6 +155,7 @@ class TestOnedim(utilities.CanteraTest): with pytest.raises(ValueError, match="mutually exclusive"): sim = cls(gas, grid=[0, 0.1, 0.2], width=0.4) + class TestFreeFlame(utilities.CanteraTest): tol_ss = [1.0e-5, 1.0e-14] # [rtol atol] for steady-state problem tol_ts = [1.0e-4, 1.0e-11] # [rtol atol] for time stepping @@ -1315,11 +1316,6 @@ class TestDiffusionFlame(utilities.CanteraTest): with pytest.raises(ct.CanteraError, match="right control point location is not set"): sim.right_control_point_temperature - - - - - class TestCounterflowPremixedFlame(utilities.CanteraTest): # Note: to re-create the reference file: # (1) set PYTHONPATH to build/python.