Counterflow Diffusion Flame Validation

Docs/sphinx/manual/Validation.rst

@@ -375,3 +375,44 @@ subgrid-scale viscosity for the :math:`Re_{\tau}` = 934 case.
 
 Further work is ongoing to assess how to better handle large LES filter size changes in the context
 of AMR-LES
+
+Laminar counterflow diffusion flame
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+The laminar counterflow diffusion flame is a fundamental test case for non-premixed reactive flows. A 
+laminar counterflow diffusion flame setup is available in ``Exec/Production/CounterFlow``. 
+In the present configuration, the geometry consists of a 1.5 cm x 1.5 cm domain with inflow at the left 
+and right boundaries and outflow at the top and bottom boundaries. The left and right inflow boundaries 
+correspond to the oxidizer and fuel inlet sides respectively, with an oxidizer/fuel jet section of 
+diameter :math:`d_{jet} = 0.5 cm` (``prob.jet_radius`` of 0.0025) centered at the boundary. The
+remainder of the inflow boundary consists of inert inflow sections with composition of pure N2. 
+The oxidizer jet section composition is that of air, while the fuel jet section composition is pure 
+gaseous Dodecane. The inert inflow velocity on the left and right is set equal to the oxidizer and fuel 
+jet velocities, respectively, in order to eliminate shear interactions between the inert and jet 
+streams. The oxidizer and fuel jet velocities are 15.0 cm/s and 6.17 cm/s, respectively, corresponding 
+to a global strain rate of :math:`40 s^{-1}`. The initial condition corresponds to air composition 
+filling the left half of the boundary and gaseous Dodecane composition filling the right half of the 
+boundary with a hyperbolic tangent smoothly transitioning species mass fractions in the center of the 
+domain. Ignition is performed by patching a spherical kernel of elevated temperature (in this case 1000 K) 
+with radius of 0.3 cm in the center of the domain as to cause autoignition. Grid refinement criteria based 
+on the gradient of temperature is used and solution convergence is achieved with ``amr.max_level`` of 2.
+
+The `PeleLMeX` results obtained using the specified grid refinement criteria are compared to that of 
+Cantera in spatial and mixture fraction space in :numref:`CounterflowLMeXCantera_Spatial` and 
+:numref:`CounterflowLMeXCantera_Mixture`, respectively. Centerline profiles of major species as well as 
+temperature across the flame front at steady-state are displayed, with PeleLMeX results indicated by solid colored lines and Cantera results indicated by black ticked lines. The results show very good agreement between the PeleLMeX and Cantera centerline profiles, with only small differences visible when plotting versus spatial coordinate due to differences in geometry (1-D vs. 2-D). 
+
+
+.. figure:: images/validations/CounterflowFlame/Dodecane_Counterflow_Spatial_Dodecane_Added.png
+   :name: CounterflowLMeXCantera_Spatial
+   :align: center
+   :figwidth: 50%
+
+   : One-dimentional flame spatial profiles with PeleLMeX and Cantera (underlying black dashed lines)
+
+.. figure:: images/validations/CounterflowFlame/Dodecane_Counterflow_Mixture_Dodecane_Added.png
+   :name: CounterflowLMeXCantera_Mixture
+   :align: center
+   :figwidth: 50%
+
+   : One-dimentional flame profiles (mixture fraction space) with PeleLMeX and Cantera (underlying black dashed lines)

Exec/Production/CounterFlow/CMakeLists.txt

@@ -1,6 +1,6 @@
 target_include_directories(pelelmex PUBLIC $<BUILD_INTERFACE:${CMAKE_CURRENT_LIST_DIR}>)
 
-set(PELE_CHEMISTRY_MODEL drm19 PARENT_SCOPE)
+set(PELE_CHEMISTRY_MODEL dodecane_lu PARENT_SCOPE)
 set(PELE_EOS_MODEL       Fuego PARENT_SCOPE)
 set(PELE_TRANSPORT_MODEL Simple PARENT_SCOPE)
 set(PELE_DIM "2" PARENT_SCOPE)

Exec/Production/CounterFlow/GNUmakefile

@@ -29,9 +29,9 @@ THREAD_SANITIZER = FALSE
 USE_EFIELD = FALSE
 
 # PelePhysics
-Chemistry_Model = drm19
+Chemistry_Model = dodecane_lu
 Eos_Model = Fuego
 Transport_Model = Simple
 
 PELE_HOME ?= ../../..
-include $(PELE_HOME)/Exec/Make.PeleLMeX
+include $(PELE_HOME)/Exec/Make.PeleLMeX

Exec/Production/CounterFlow/PeleLMeX_PatchFlowVariables.cpp

@@ -0,0 +1,68 @@
+
+#include "PeleLMeX.H"
+using namespace amrex;
+
+void
+patchFlowVariables(
+  const amrex::Geometry& geom, ProbParm const& lprobparm, amrex::MultiFab& a_mf)
+{
+
+  amrex::Print() << "\nPatching flow variables..";
+
+  const amrex::Real* prob_lo = geom.ProbLo();
+  const amrex::Real* prob_hi = geom.ProbHi();
+  const amrex::Real* dx = geom.CellSize();
+
+  for (amrex::MFIter mfi(a_mf, amrex::TilingIfNotGPU()); mfi.isValid(); ++mfi) {
+    const amrex::Box& bx = mfi.tilebox();
+
+    auto const& rho_arr = a_mf.array(mfi, DENSITY);
+    auto const& rhoY_arr = a_mf.array(mfi, FIRSTSPEC);
+    auto const& rhoH_arr = a_mf.array(mfi, RHOH);
+    auto const& temp_arr = a_mf.array(mfi, TEMP);
+    Real massfrac[NUM_SPECIES] = {0.0};
+
+    amrex::ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept {
+      auto eos = pele::physics::PhysicsType::eos();
+      amrex::Real massfrac[NUM_SPECIES] = {0.0};
+      amrex::Real sumYs = 0.0;
+
+      for (int n = 0; n < NUM_SPECIES; n++) {
+        massfrac[n] = rhoY_arr(i, j, k, n);
+#ifdef N2_ID
+        if (n != N2_ID) {
+          sumYs += massfrac[n];
+        }
+#endif
+      }
+#ifdef N2_ID
+      massfrac[N2_ID] = 1.0 - sumYs;
+#endif
+
+      AMREX_D_TERM(const amrex::Real x = prob_lo[0] + (i + 0.5) * dx[0];
+                   , const amrex::Real y = prob_lo[1] + (j + 0.5) * dx[1];
+                   , const amrex::Real z = prob_lo[2] + (k + 0.5) * dx[2];);
+
+      AMREX_D_TERM(const amrex::Real Lx = prob_hi[0] - prob_lo[0];
+                   , const amrex::Real Ly = prob_hi[1] - prob_lo[1];
+                   , const amrex::Real Lz = prob_hi[2] - prob_lo[2]);
+
+      AMREX_D_TERM(const amrex::Real xc = prob_lo[0] + Lx / 2.0;
+                   , const amrex::Real yc = prob_lo[1] + Ly / 2.0;
+                   , const amrex::Real zc = prob_lo[2] + Lz / 2.0);
+
+      amrex::Real radiusSq = AMREX_D_TERM(
+        (x - xc) * (x - xc), +(y - yc) * (y - yc), +(z - zc) * (z - zc));
+
+      amrex::Real radius = std::sqrt(radiusSq);
+      amrex::Real mixingWidth = 0.1 * lprobparm.ignitSphereRad;
+      amrex::Real mixingFunction =
+        0.5 *
+        (1.0 + std::tanh((lprobparm.ignitSphereRad - radius) / mixingWidth));
+      temp_arr(i, j, k) = mixingFunction * lprobparm.ignitT +
+                          (1.0 - mixingFunction) * lprobparm.T_inert;
+    });
+  }
+
+  amrex::Print() << "Done\n";
+}

Exec/Production/CounterFlow/README.md

@@ -0,0 +1,13 @@
+## Counterflow Diffusion Flame
+
+Sample test case for running 2-D counterflow diffusion flame. The initial solution is obtained by first
+running a "coolflow" case using the input file ``input_coolflow.2d-regt`` without grid refinement to reduce
+computational cost (Note ``prob.do_ignition`` is set to 0). The domain will fill with oxidizer and fuel until
+a steady-state is reached when the oxidizer and fuel meet at approx. the center of the domain (Visualization
+software such as Paraview can be used to determine when steady-state is reached).
+
+The counterflow flame is then ignited using an ignition kernel by setting ``prob.do_ignition = 1`` and 
+``peleLM.initDataPlt_patch_flow_variables = true`` in the input file ``input.2d-regt``, which restores 
+the steady-state coolflow plotfile specified by ``amr.initDataPlt``. The ignition kernel parameters 
+``prob.ignition_SphT`` and ``prob.ignition_SphRad`` are used in ``PeleLMeX_PatchFlowVariables.cpp`` to 
+specify the patched kernel temperature and radius, respectively. 

Exec/Production/CounterFlow/input.2d-regt

@@ -1,8 +1,8 @@
 #----------------------DOMAIN DEFINITION------------------------
 geometry.is_periodic = 0 0               # For each dir, 0: non-perio, 1: periodic
-geometry.coord_sys   = 0                  # 0 => cart, 1 => RZ
-geometry.prob_lo     =  -0.01 -0.01 0.0        # x_lo y_lo (z_lo)
-geometry.prob_hi     =   0.01  0.01 0.016        # x_hi y_hi (z_hi)
+geometry.coord_sys   = 0                 # 0 => cart, 1 => RZ
+geometry.prob_lo     =  -0.0075 -0.0075 0.0		# x_lo y_lo (z_lo)
+geometry.prob_hi     =   0.0075  0.0075 0.0		# x_hi y_hi (z_hi)
 
 # >>>>>>>>>>>>>  BC FLAGS <<<<<<<<<<<<<<<<
 # Interior, Inflow, Outflow, Symmetry,
@@ -11,51 +11,72 @@ peleLM.lo_bc = Inflow Outflow
 peleLM.hi_bc = Inflow Outflow
 
 
-#-------------------------AMR CONTROL----------------------------
-amr.n_cell          = 64 64 32      # Level 0 number of cells in each direction   
+#-------------------------AMR CONTROL----------------------------  
+amr.n_cell          = 64 64 32	  	   # Level 0 number of cells in each direction
 amr.v               = 1                # AMR verbose
-amr.max_level       = 0                # maximum level number allowed
+amr.max_level       = 2                # maximum level number allowed
 amr.ref_ratio       = 2 2 2 2          # refinement ratio
 amr.regrid_int      = 5                # how often to regrid
 amr.n_error_buf     = 2 2 2 2          # number of buffer cells in error est
 amr.grid_eff        = 0.7              # what constitutes an efficient grid
-amr.blocking_factor = 16               # block factor in grid generation (min box size)
+amr.blocking_factor = 4
 amr.max_grid_size   = 64               # max box size
 
-
 #--------------------------- Problem -------------------------------
-prob.P_mean = 101325.0
-prob.T_ox   = 298.0
-prob.T_fuel = 298.0
-prob.massflow = 1.0
-prob.jet_radius = 0.005
-prob.inert_radius = 0.0075
-prob.inert_velocity = 0.2
-prob.pertmag = 0.000
-prob.pmf_datafile = "drm_CH4Air_stoich.dat"
+#prob.P_mean = 101325.0
+prob.P_mean = 792897.5 #7.82529*one_atm 
+
+prob.T_ox   = 450.0
+prob.T_fuel = 450.0
+prob.T_inert = 450.0
+prob.Y_O2_ox = 0.233
+prob.Y_N2_ox = 0.767
+prob.Y_N2_fuel = 0.0
+prob.Y_fuel = 1.0
+
+prob.massflow_ox = 0.9171
+prob.massflow_fuel = 2.228
+prob.jet_radius = 0.0025
+prob.inert_velocity_ox = 0.1500
+prob.inert_velocity_fuel = 0.0617
+
 prob.do_ignition = 1
-prob.ignition_SphT = 2200.0
+prob.ignition_SphT = 1000.0
+prob.ignition_SphRad = 0.0015
 
 #-------------------------PeleLM CONTROL----------------------------
-peleLM.v = 1
-peleLM.incompressible = 0
-peleLM.rho = 1.17
-peleLM.mu = 0.0
-peleLM.use_wbar = 1
-peleLM.sdc_iterMax = 2
-peleLM.floor_species = 1
-peleLM.deltaT_verbose = 0
-
-#amr.restart = chk01500
-#amr.regrid_on_restart = 1
-amr.check_int = 100
-amr.plot_int = 20
-amr.max_step = 100
+peleLM.v = 3 						# [OPT, DEF=0] Verbose
+peleLM.incompressible = 0 			# [OPT, DEF=0] Enable to run fully incompressible, scalar advance is bypassed
+peleLM.rho = 1.17					# [OPT, DEF=-1] If incompressible, density value [MKS]
+peleLM.mu = 0.0						# [OPT, DEF=-1] If incompressible, kinematic visc. value [MKS]
+peleLM.use_wbar = 1 				# Include Wbar term in species diffusion fluxes
+peleLM.sdc_iterMax = 2				# Number of SDC iterations
+peleLM.floor_species = 1			# [OPT, DEF=0] Crudely enforce mass fraction positivity
+peleLM.deltaT_verbose = 0			# [OPT, DEF=0] Verbose of the deltaT iterative solve algorithm
+
+#amr.restart = chk2000
+amr.initDataPlt = plt02656_coolflow
+peleLM.initDataPlt_patch_flow_variables = true
+amr.regrid_on_restart = 1
+
+amr.check_int = 2000
+#amr.plot_int = 200
+
+amr.plot_per = 0.005  				#Plot every t=5ms
+amr.plot_per_exact = 1            	# [OPT, DEF=0] Flag to enforce exactly plt_per by shortening dt
+
+#amr.max_step = 1000
 amr.dt_shrink = 0.01
-amr.stop_time = 1.1
-#amr.stop_time = 1.00
-amr.cfl = 0.25
-amr.derive_plot_vars = avg_pressure mag_vort mass_fractions
+amr.init_dt = 1.0e-6
+amr.stop_time = 0.150
+amr.cfl = 0.1
+amr.derive_plot_vars = avg_pressure mag_vort mass_fractions mixture_fraction
+
+peleLM.fuel_name = NC12H26
+peleLM.mixtureFraction.format = Cantera
+peleLM.mixtureFraction.type   = mass
+peleLM.mixtureFraction.oxidTank = O2:0.233 N2:0.767
+peleLM.mixtureFraction.fuelTank = NC12H26:1.0
 
 peleLM.chem_integrator = "ReactorCvode"
 peleLM.use_typ_vals_chem = 1          # Use species/temp typical values in CVODE
@@ -64,39 +85,21 @@ ode.atol = 1.0e-5                     # Absolute tolerance factor applied on typ
 cvode.solve_type = denseAJ_direct               # CVODE Linear solve type (for Newton direction) 
 cvode.max_order  = 4                  # CVODE max BDF order. 
 
-godunov.use_ppm = 0
-godunov.use_forceInTrans = 0
-
 nodal_proj.verbose = 0
 mac_proj.verbose = 0
-#diffusion.verbose = 2
+diffusion.verbose = 0
+mac_proj.rtol = 1.0e-10
+nodal_proj.rtol = 1.0e-10
 
 peleLM.do_temporals = 1
 peleLM.do_mass_balance = 1
 
 #--------------------REFINEMENT CONTROL------------------------
-amr.refinement_indicators = temp
-amr.temp.max_level     = 2
-amr.temp.value_greater = 1500
-amr.temp.field_name    = temp
-
-#amr.refinement_indicators = magVort
-#amr.magVort.max_level     = 1
-#amr.magVort.value_greater = 500.0
-#amr.magVort.field_name    = mag_vort
-
-#amr.refinement_indicators = yH_Crse yH_Fine CH2O
-#amr.yH_Crse.max_level     = 1
-#amr.yH_Crse.value_greater = 1.50e-4
-#amr.yH_Crse.field_name    = Y(H)
-#
-#amr.yH_Fine.max_level     = 3
-#amr.yH_Fine.value_greater = 2.00e-4
-#amr.yH_Fine.field_name    = Y(H)
-#
-#amr.CH2O.max_level     = 4
-#amr.CH2O.value_greater = 1.00e-3
-#amr.CH2O.field_name    = Y(CH2O)
+amr.refinement_indicators = gradT
+
+amr.gradT.max_level = 2
+amr.gradT.adjacent_difference_greater = 100
+amr.gradT.field_name = temp
 
 #amrex.fpe_trap_invalid = 1
 #amrex.fpe_trap_zero = 1