# Design and Formulation of a Hydromechanical Fin for AUV Operations and Wave
Parameter Estimation

Brendan M. Unikewicz
Department of Ocean Engineering
University of Rhode Island
Narragansett, RI, 02882
Current affiliation: Department of Mechanical Engineering, Massachusetts
Institute of Technology, 77 Massachusetts Avenue Cambridge, MA 02139; please
contact<EMAIL_ADDRESS>for any inquiries.

(May 1, 2023)

###### Abstract

Ocean dynamics play a crucial role in global climate, ecosystems, and human
activities, necessitating accurate and efficient methods to characterize ocean
currents and waves. This paper presents the development of a novel bio-
inspired hydrofoil system for detecting and characterizing ocean currents and
waves based on a dolphin’s flipper foil design. The prototype’s performance
was assessed through controlled experiments, demonstrating the system’s
ability to quickly and accurately orientate itself in the direction of flow.
Wave mechanics and a geometric proof were applied to estimate wave parameters
such as wave height and current from the hydrofoil’s positioning. The proposed
hydrofoil system shows potential for use in various marine applications,
including oceanographic research, environmental monitoring, and navigation.
The bio-inspired hydrofoil system offers a promising approach to ocean current
and wave characterization, with the potential to significantly impact our
understanding and monitoring of ocean dynamics without greatly impacting
vehicle performance or increasing power utilization in operation.

_Keywords_ Hydrofoil $\cdot$ Ocean Waves $\cdot$ Autonomous Underwater Vehicle
(AUV) $\cdot$ Particle Trajectory $\cdot$ Wave Theory

##  1 Introduction

Hydrofoils have been used as propulsion methods for numerous undersea vehicles
[Sahoo et al.]. One of the more glaring benefits to utilizing hydrofoils is
that they grant increased vehicle mobility, providing added freedom and
degrees of movement towards undersea applications [Jaya and Kartidjo].
Additionally, hydrofoils are paralleled with species of Testudines. We have
seen research journals mimicking sea turtle flipper motion that address the
benefits of biomimicry; however, research has only gone as far as exploring
thrust capabilities and mobility characteristics within AUV and ROV systems
[Liu et al.].

Hydrofoils have the potential to grant increased range of operation within
AUVs via decreasing power consumption as well as offer the potential benefit
of being a data acquisition (DAQ) system for ocean-surface conditions.
Undersea vehicles have a lot to adapt from biomimicry and this project takes
an approach towards the potential integration of rheotaxis within an undersea
vehicle’s propulsion method, i.e. hydrofoils.

Rheotaxis is a behavioral trait seen in many aquatic vertebrate that will
allow them to orientate themselves in the direction of a fluid current through
a system of sensory organs found along the lateral line system [Jiang et al.].
The major unit of functionality of the lateral line is the neuromast. There
are two main varieties of neuromasts located in animals, canal neuromasts and
superficial or freestanding neuromasts. Superficial neuromasts are located
externally on the surface of the body, while canal neuromasts are located
along the lateral lines in subdermal, fluid filled canals [Webb et al.].
Through these discrete sensors along the body, the aquatic vertebrates may
detect movement, vibration, and pressure gradients and choose to either
maintain their current position or align themselves parallel to the direction
of flow. We will be modeling our work after the superficial neuromast, a
discrete sensor located flush with the boundary of the foil.

If we treat the aquatic vertebrate as a bluff body, i.e. an ellipsoid shape,
and the vertebrate chooses to align itself parallel to the direction of flow
then the lift force exerted on the body is a minimum, reducing the amount of
power required to move in a singular direction. This rationale of thinking may
be directly translated towards AUV systems and we can expect to see better
mission run positioning, tracking, and decreased power consumption through
multiple operations.

This does offer the question as to why we do not explore the integration of
discrete sensors along the body of an AUV as opposed to the intended
hydrofoil. This is done because we have the unique opportunity to not only
align ourselves in the direction of flow, but potentially use the same flow-
orientation scheme to measure ocean surface conditions. Through the smart
hydrofoil proposed in this paper, only one foil is required for decreasing
power requirements in AUVs with the additional advantage of characterizing
surface conditions at or near the upper ocean layer.

This paper reviews the theory as to why flow-orientation schemes provide
decreased power consumption in AUVs, validation of an in-air prototype for
potential wave tank testing, and the mathematical formulation for and design
of a hydrofoil capable of acting as a standalone DAQ system to characterize
currents and particle trajectories induced by gravity waves.

##  2 Design of a Foil Prototype

###  2.1 Effective Drag in AUV Systems

Previous works have investigated relationships between power, range, and speed
for AUVs [Singh et al.] that co-align with the goals of this research and a
prevalent assumption is that the vehicle drag is largely related to the
effective drag on the body and the relative speed of the body in the ocean
currents [Granville]. Using simplified assumptions, our aim is to decrease the
effective drag coefficient as much as possible, subsequenly increasing vehicle
range in AUV operations. This is done because power is widely considered to be
the most limiting factor in AUV design.

The smart hydrofoil’s sole purpose is to be able to accurately detect and
align itself with the dominant direction of flow. It is assumed, that a full
AUV system, upon determining the dominant direction of flow, will align itself
through additional on-board processing. Determining optimum shape parameters
for an AUV as well as how each on-board componen (i.e. propulsion methods,
depth sensors, etc.) draws power is a tedious task. It is much simpler to
consider the AUV as a bluff body with the intention of decreasing the form
drag. This drag force is a result of not only the shape of the body, but also
the angle of attack, described by basic forms of the drag force and force due
to lift [Bishop and Hassan].

###  2.2 Theoretical Approach to Reducing Effective Drag

For determining the dominant direction of flow in the prototype, we can ignore
the inertial and Froude-Krylov contributions towards Morison’s equation [Sumer
and Fredsøe]. This is because we are largely interested in reducing the
coefficient of lift, $C_{L}$, which is a function of the angle of attack, and
in doing so reducing the system drag as it hovers in the water column.

$F(t)=\rho C_{m}V_{b}\frac{du}{dt}+\frac{1}{2}C_{d}A_{b}u|u|$ (1)

Where, $F(t)$ is the total inline force on the body, $\frac{du}{dt}$, is the
flow acceleration, $V_{b}$, is the volume of the body, and the Froude-Krylov
force is $\rho V_{b}\frac{du}{dt}$. Multiplying, the Froude-Krylov force by
the inertia coefficient, $C_{m}$, where the inertia coefficient is simply
$1+C_{a}$, i.e. the added mass coefficient, we yield the total inertial force
which has been omitted for our considerations. Further, our area of interest,
the force due to drag, can be seen in the simplified version of the drag
equation, $\frac{1}{2}\rho C_{d}A_{b}u|u|$, where the coefficient of drag,
$C_{d}$, and the cross-sectional area of the body perpendicular to the
direction of the flow is $A_{b}$.

Additionally, the added forces in Morison’s equation for a moving body in an
oscillatory flow only provide an estimate of the variation about the mean
value that we are concerned with. At an angle of attack of zero degrees (and
depending on the foil parameters) the coefficient of lift will be nearly zero.
If we were to use the lift force along the body, a nonlinear function of the
current velocity, the minimum value would indicate the foil is aligned in the
dominant direction of flow and maximally reduce the lift force; however, we
can use even more simplified assumptions.

Assuming we have an approximately constant flow field for the top and bottom
of our fin platform, in the vertical plane (x- and z-axis), and apply mass
conservation laws, we know that the amount of fluid passing over the top of
the structure is equivalent to the amount of fluid passing underneath. In
potential flow theory, this applies to all angle–of–attacks, and you could
calculate the lift force; however, we want to utilize the resulting disparity
in flow speed to reorient ourselves in the mean, dominant direction of flow as
well as characterize particle trajectories at the upper ocean layer. The
application of these assumptions allows the avoidance of incorporating
dynamically varying pressure sensors, which may have brought complications for
near-surface operations due to the presence of the dynamic free surface
boundary condition.

Figure 1: CAD model of the symmetric dolphin-inspired foil. The cross-shaped
exclusion indicates where the servo motor is fitted, in-line with the center
of gravity. Hot-wire anemometers fit in each rectangular slot on the fin
platform and wires are funneled through the trailing edge.

Due to symmetry about the chord, when our hydrofoil is aligned in the
direction of flow, the flow over the top will be equivalent to the flow
underneath the hydrofoil. For in-air applications, a strong choice to measure
flow speeds are anemometers. The anemometers of choice utilize a constant
temperature circuitry that are resistor dependent. This was seen as ideal
because resistors could be exchanged as needed.

Constant temperature anemometers contain two fixed resistors, $R_{1}$ and
$R_{2}$, a temperature dependent resistor, $R_{w}$, and lastly a variable
resistor, $R_{3}$, to bring the bridge back into balance. As flow passes over
the temperature dependent resistor the unbalancing of the resistor bridge
results in a voltage difference, which is compensated through the current in
the system [Perry and Morrison]. The anemometers that were used required
3.3VDC and 20-40mA, depending on the wind velocity recorded.

$R_{1}R_{3}=R_{2}R_{w}$ (2)

In order to move the position of the hydrofoil, in accordance with the
resultant anemometer readings, a servo motor was used. With the anemometers
and a single servo motor, we can now describe code, in C/C++, to archive data
as well as reorient our hydrofoil.

###  2.3 DAQ Software Design

The DAQ software for the hydrofoil begins by moving the fin to a zero position
of 90°. A data file is established within the SD card and both a timer and
real-time clock is initialized. Once approximately 10 seconds has passed,
allowing the temperature dependent resistors to stabilize, the anemometers,
located symmetrically on opposite sides of the fin platform, are read into the
micro-controller as an analog value and a percent difference is taken. For our
validation, if the percent difference is below some threshold value, i.e. 5%,
then we are currently aligned in the dominant direction of flow and the sample
rate remains a constant value of 2 Hz. The data (servo position, anemometer
values, and percent difference) is then archived and this process is repeated.
This lower sample rate is introduced to minimize on-board processing and the
archiving of data.

If the percent difference is above the indicated threshold, the foil moves
towards the anemometer of higher value, introduces a higher sample rate
(20Hz), and archives the data once again. This process should introduce more
wind velocity to the lower valued anemometer until both anemometers are within
the threshold once again. Once the percent difference is within the threshold
once more, the sampling rate is reduced to the previous value.

This simplistic approach [Fig. 6] has room for improvement, i.e. releasing the
processor from waiting for peripheral data, improve data transfer throughput,
reducing the rate at which interrupts are generated, and lastly, mechanical
design considerations such as St. Venant’s principle, etc.; however, it’s a
strong validation approach when considering the simplicity of integration and
number of computations required. Additionally, the variable sampling rate is a
key factor in this system due to the availability of power and the ability of
the foil to characterize multi-frequency waves. In an ocean environment waves
are not monochromatic, and at a minimum the Nyquist frequency must be
satisfied to characterize the highest frequency encountered [Robinson and
Clark]; however, it will be necessary to increase the sampling rate further,
10-20x the Nyquist frequency, in order to fully characterize the particle
trajectory motions.

###  2.4 Complete System Design

For the first prototype system a laptop was used as the power source and
programming system. Once the electrical system was fully defined it was time
to design and additively manufacture the prototype fin.

The parameters for the fin were loosely based on a dolphin’s flipper foil
[Fig. 1]; however, the chord and chamber lines are the same length, 0.234
meters indicating a symmetric foil. Additionally, the maximum thickness of the
hydrofoil is 0.059 meters and the span is 0.15 meters, resulting in an aspect
ratio of 0.641. The servo motor was placed at the center of gravity, along the
chamber, which was 0.111 meters.

The first iteration fin is capable of detecting unidirectional flows,
operating in the x-y coordinate plane and is oriented orthogonally to such
plane. This was done to forgo the moment arm that would occur if the foil was
parallel to the x-y coordinate plane. Additionally, the foil is hollowed out
for decreased system mass; in water these considerations may be less of a
concern due to Archimedes’ principle, and a parallel orientation with the
ocean surface will be necessary to characterize the sea state above. The foil
was additively manufactured with a cartesian based 3-D printer using PLA
filament and mounted to a wooden board for stability, simulating a rigid body,
i.e. an AUV/ROV [Fig. 5].

##  3 Assessing Prototype Capabilities

###  3.1 Experimental Setup

Testing of the fin as a data acquisition system included a low-power fan
(120VAC, 0.25A, 60Hz) capable of in-line flow at the source of approximately
3.5 m/s and controlled angles at which the flow would reach the leading edge.
The fan is moved along the path line of a circle, equidistant to the center of
gravity for the hydrofoil, i.e. the servo motor position, 25cm away. The
hydrofoil begins at a resting position of 90°before any data is collected. The
fan faces the leading edge of the foil and is moved to angular positions of
90°, 140°, 105°, 60°, 150°, 105°, 150°, and 130°over the course of
approximately 50 seconds. The testing is designed to assess the ability of the
foil to accurately relay position, over time, to the micro-controller for
post-processing and future testing in an ocean environment.

Figure 2: Experimental testing setup; top-down view.

###  3.2 Experimental Results

After calibrating the anemometers, we see they exhibit a percent difference,
with no flow speed, of approximately 1.72% with anemometer one having a
resting voltage of 0.58 volts and anemometer two with 0.59 volts. The first 10
seconds of the data set were devoted to allowing the temperature dependent
resistor to stabilize. Afterwards the fan was turned on at the 140°location.
The anemometer on the left side of the fin, anemometer one, increased in value
and when the percent difference threshold was exceeded, the foil began to
converge towards that sensor, decreasing the percent difference measurement
until it was within an acceptable range [Fig. 2]. The position of the
hydrofoil depicts a clear story of what was occurring within the system,
moving to positions of approximately 90°, 140°, 105°, 60°, 150°, 105°, 150°,
and 130°throughout time. The foil was able to quickly and consistently move to
the new location and maintain its position until a new flow direction was
encountered [Fig. 3]. We can see that the rise and fall time of the foil
heading is approximately forty degrees per second with the servo motor used.
The sign of the percent differences indicate what side the flow was
originating, i.e. left or right, and we can further assess the hydrofoil’s
ability to accurately orientate itself in the direction of flow.

Figure 3: Output data from assessing the prototype capabilities featuring: the calibrated anemometer readings (top), the heading position in degrees (middle), and the percent difference between the two anemometer readings (bottom). Table 1:  Input flow speed angles versus measured hydrofoil angles at eight different positions. | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8
---|---|---|---|---|---|---|---|---
Intended Angle | 90.0° | 140.0° | 105.0° | 60.0° | 150.0° | 105.0° | 150.0° | 130.0°
Measured Angle | 90.0° | 138.6° | 106.2° | 61.2° | 151.2° | 104.4° | 147.6° | 131.4°
% Difference Angle | 0.00% | -1.00% | 1.14% | 2.00% | 0.80% | -0.40% | -1.60% | 1.08%

Without centering the heading of the foil with the direction of flow, there
are no percent differences larger than 2.00% and no percent differences
smaller than 0.40% (excluding the resting position). Improvement to this data
may take the form of decreasing the threshold that was used or swapping out
the hardware for more precise components. Limitations do exist with the
equipment utilized, and there may further be a need to incorporate Direct-
Memory-Access protocols into our software, working only from memory to
peripheral. When looking at the variable sampling rate data, the need for
Direct-Memory-Access becomes more apparent. We are able to satisfy the
threshold of 5% at approximately 28.65 seconds; however, the interrupt is not
generated and registered until another two samples are taken, increasing the
variability between measurements [Fig. 7]. This variable sampling rate will
help reduce the length of the outputted file, required processing, and
ultimately increasing the total range and power in the system; however, more
work must be done to increase efficiency. The benefits of a system like the
one proposed become more prevalent with its ease of operation and ability to
process the collected data. This is because the preliminary and final system
only has one purpose: orientate in and detect the direction of flow, quickly
and accurately.

##  4 Formulations for Wave Parameter Estimation

###  4.1 Application of Wave Mechanics

To obtain hydrodynamic parameters such as wave height and current,
small–amplitude wave theory, i.e. linearized wave theory, must be used. We can
begin by assuming that potential flow theory applies and we take our fluid to
be incompressible, inviscid, and irrotational.

$u=\nabla\phi$ (3) $\nabla^{2}\phi=0$ (4) $\frac{\partial\phi}{\partial
t}+\frac{1}{2}|u|^{2}+gz+\frac{p}{\rho}=0$ (5)

Small surface elevations and fluid velocities, constant atmospheric pressure,
and water depth at the survey site permit the use of the linear dispersion
relationship for our data acquisition system.

$\frac{\sigma^{2}}{g}=ktanh(kh)$ (6)

###  4.2 Geometric Proof for Particle Trajectory Estimation

It’s possible to derive flow at the stagnation point from one of the regions
along the streamlined body; however, that would be computationally intensive,
most likely requiring numerical modeling to relate the measured current to the
actual current at or near the stagnation point or even the integration of
additional sensors, which would increase the operational power required. Since
we have our future motor with a known angle, $\Theta$, located at a known
distance from the stagnation point, $G$, we can use simple geometry and a
time-scale to resolve the vertical position, $r$, and vertical component of
ocean currents, $w$, induced by gravity waves [Fig. 4]. Small angle theory is
used, which means these operations must not occur too close to the ocean
surface, increasing the vertical distance the foil travels, voiding theory.

Figure 4: Proposed wave estimation diagram for a flow-orienting fin mounted to
an AUV for deep-water ocean waves.

The following proof applies to any known mounting distance from the stagnation
point. The center of gravity was not used for this proof because center of
gravity mounting for hydrofoils may inhibit their ability to produce thrust.
Since we are applying small angle theory we are assuming that the resultant
lift and drag force from small pitching movements in the foil are too little
to contribute greatly towards lifting the rigid body of which it is attached.
If the small angle theory is voided we can reasonably anticipate phase
differences in the particle trajectories and inaccurate estimates of wave
height. Additionally, for these operations, the hydrofoil must be placed
parallel to the ocean surface.

$\theta=cos(\frac{G}{r})$ (7) $w=\frac{d}{dt}r$ (8)

The deep and shallow water regimes were chosen because they provide stark
contrasts in wave mechanics. Deep-water waves experience dispersion, i.e.
individual wave celerity is greater than group celerity while shallow-water
waves experience no dispersion because all waves travel at the same speed,
including the group celerity [Newman]. This contrasting of dispersion provides
a wide frequency spectrum for the hydrofoil to be future tested. Depending on
the depth we lie, we can use the following relationships from solutions of the
linear dispersion relationship to achieve wave height and particle trajectory
information from our measurements.

$\frac{h}{L}>\frac{1}{2}$ (9) $w=\frac{g}{c}\frac{H}{2}e^{kz}sin(kx-wt)$ (10)
$\frac{h}{L}<\frac{1}{20}$ (11)
$w=\frac{g}{c}\frac{kz(cosh(kh)+khcosh(kz))}{cosh(kh)}sin(kx-wt)$ (12)

Where, $H$, is the wave height in meters, $c$, is the wave celerity in meters
per second, $k$, is the wave number, $L$, is the wavelength and, $h$, is the
depth in meters. It’s important to note that the only other component required
in order to complete this formulation is the integration of a depth sensor,
which is commonly found in AUV systems. These formulations can be processed,
real-time, on AUV systems through on-board processing and have the potential
to greatly increase the performance of existing fin-type propulsion systems
and/or systems seeking to use novel sensing mechanisms for wave parameter
estimation.

##  5 Conclusion

This paper presented the development of a novel bio-inspired hydrofoil system
for the purpose of detecting and characterizing ocean currents and waves. The
prototype was designed based on a dolphin’s flipper foil, which provided an
efficient and responsive platform for data acquisition. The system’s
capabilities were tested using a controlled experimental setup with varying
flow directions, and the results demonstrated the hydrofoil’s ability to
quickly and accurately orientate itself in the direction of flow for future
translation underwater.

The application of wave mechanics for particle trajectory estimation allowed
for the derivation of wave height and flow information to be formulated based
on the smart foil’s positioning and sensing principles. This method has the
potential to be utilized in various oceanic environments, from deep to shallow
water, providing valuable data for various marine applications such as
oceanographic research, environmental monitoring, and navigation.

While the initial prototype showed promising results, some areas of
improvement were identified, such as refining the variable sampling rate and
implementing Direct-Memory-Access protocols to increase the efficiency of the
data acquisition process. Additionally, further testing in real ocean
environments is necessary to validate the system’s performance under more
complex and dynamic conditions.

In conclusion, the bio-inspired hydrofoil system presents a promising approach
for ocean current and wave characterization, with the potential to
significantly impact marine applications. By continuing to refine the design
and testing the system in real-world conditions, this technology could become
a valuable tool for understanding and monitoring ocean dynamics.

##  6 Acknowledgements

The author wishes to acknowledge faculty from the University of Rhode Island’s
Ocean Engineering Department for approving this research as satisfactory for
the partial fulfillment of a Master’s of Science degree in Ocean Engineering
and further permitting its contents for electronic distribution to arXiv.

## References

  * Sahoo et al. [2019] Avilash Sahoo, Santosha K. Dwivedy, and P.S. Robi.  Advancements in the field of autonomous underwater vehicle.  _Ocean Engineering_ , 181:145–160, 2019.  ISSN 0029-8018.  doi:https://doi.org/10.1016/j.oceaneng.2019.04.011.  URL https://www.sciencedirect.com/science/article/pii/S0029801819301623.
  * Jaya and Kartidjo [2022] Arie Sukma Jaya and Muljo Widodo Kartidjo.  Thrust and efficiency enhancement scheme of the fin propulsion of the biomimetic autonomous underwater vehicle model in low-speed flow regime.  _Ocean Engineering_ , 243:110090, 2022.  ISSN 0029-8018.  doi:https://doi.org/10.1016/j.oceaneng.2021.110090.  URL https://www.sciencedirect.com/science/article/pii/S0029801821014177.
  * Liu et al. [2022] Yuhong Liu, Houcun Li, Shihan Deng, Shuxin Wang, Shuhe Liu, and Zhijie Wang.  Biomimetic robotic sea lion foreflippers: Design, modeling, and experimentation.  _IEEE/ASME Transactions on Mechatronics_ , 27(6):5679–5689, 2022.  doi:10.1109/TMECH.2022.3187014.
  * Jiang et al. [2019] Yonggang Jiang, Zhiqiang Ma, and Deyuan Zhang.  Flow field perception based on the fish lateral line system.  _Bioinspiration & biomimetics_, 14(4):041001, 2019.
  * Webb et al. [2021] Jacqueline F Webb, Karen P Maruska, Julie M Butler, and Margot AB Schwalbe.  The mechanosensory lateral line system of cichlid fishes: from anatomy to behavior.  _The behavior, ecology and evolution of cichlid fishes_ , pages 401–442, 2021.
  * Singh et al. [1997] Hanumant Singh, Dana Yoerger, and Albert Bradley.  Issues in auv design and deployment for oceanographic research.  In _Proceedings of International Conference on Robotics and Automation_ , volume 3, pages 1857–1862. IEEE, 1997.
  * Granville [1976] Paul S Granville.  Elements of the drag of underwater bodies.  Technical report, DAVID W TAYLOR NAVAL SHIP RESEARCH AND DEVELOPMENT CENTER BETHESDA MD SHIP …, 1976.
  * Bishop and Hassan [1964] Richard Evelyn Donohue Bishop and AY Hassan.  The lift and drag forces on a circular cylinder oscillating in a flowing fluid.  _Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences_ , 277(1368):51–75, 1964\.
  * Sumer and Fredsøe [2006] B. Mutlu Sumer and Jørgen Fredsøe.  Hydrodynamics around cylindrical structures.  2006\.
  * Perry and Morrison [1971] A. E. Perry and G. L. Morrison.  A study of the constant-temperature hot-wire anemometer.  _Journal of Fluid Mechanics_ , 47(3):577–599, 1971.  doi:10.1017/S0022112071001241.
  * Robinson and Clark [1991] Enders Robinson and Dean Clark.  Sampling and the nyquist frequency.  _The Leading Edge_ , 10(3):51–53, 1991.
  * Newman [2018] John Nicholas Newman.  _Marine hydrodynamics_.  The MIT press, 2018.

##  7 Appendix

Further information regarding the testing and assembly may be seen, below, in
the Appendix of this paper.

Figure 5: Complete in-air prototype hydrofoil for unidirectional flow-sensing
featuring the electronic components utilized and 3-D printed CAD model. Figure
6: Procedural flow diagram outlining data acquisition and controls execution.
Figure 7: Variable sampling rate during data acquisition.


# NOx emissions trends in hydrogen lean premixed flamelets at high strain rate

A. Porcarelli∗, B. Kruljevic and I. Langella
TU Delft, Faculty of Aerospace Engineering, Flight Performance & Propulsion
group
Kluyverweg 1, 2629 HS Delft, Netherlands
<EMAIL_ADDRESS>

###### Abstract

NOx formation in lean premixed and highly-strained pure hydrogen-air flamelets
is investigated numerically. Lean conditions are established at an equivalence
ratio of 0.7. Detailed-chemistry, one-dimensional simulations are performed on
a reactants-to-products counter-flow configuration with an applied strain rate
ranging from $a=100\,{\rm s}^{-1}$ to $a=10000\,{\rm s}^{-1}$ and the GRI3.0
mechanism. Following a similar setup, two-dimensional direct numerical
simulations are also conducted for representative strain rates of $2000\,{\rm
s}^{-1}$ and $5000\,{\rm s}^{-1}$. Both solutions show a decreasing NOx trend
as the applied strain rate is increased. This decreasing emission outcome is
highlighted for the first time in this study for lean pure-hydrogen flamelets.
A deep analysis of the 2D solution underlines that there is no production of
NOx in the second dimension, thus proving that the emission trend is not a
result of a setup preconditioning, but is instead a direct physical effect of
stretch on the flame. Furthermore, a detailed analysis of the NOx formation
pathways at $a=2000\,{\rm s}^{-1}$ and $a=5000\,{\rm s}^{-1}$ is performed.
Thermal NOx and NNH pathways are shown to both contribute significantly to the
total NOx production. While the NNH route contribution is roughly constant at
different strain rates, a significant decrease is observed along the thermal
NOx route. Overall, results show that lean and highly-strained hydrogen flames
experience a significant decrease of NOx. This property is discussed and
analysed in the paper.

## Introduction

The yearly increasing energy demand is currently substantially met by fossil
fuels, with consequent emissions of CO2 and other pollutants. Hydrogen has
been identified as one of the possible alternatives to satisfy this demand and
reduce green-house emissions simultaneously. In fact, hydrogen is carbon-free,
and can be produced with electrolysis using renewable energy [1, 2]. In the
last decades, research efforts have focused on the possibility to burn
hydrogen in lean premixed conditions, where the lower adiabatic flame
temperature allows to decrease NOx formation via the thermal route.
Furthermore, hydrogen’s strong reactivity and high lower-heating-value allows
to reach ultra-lean conditions without approaching lean blow-off. Earlier
studies have discussed the influence of hydrogen addition on lean blow-off,
showing that a minimal hydrogen enrichment can consistently decrease the blow-
off equivalence ratio and thus broaden the possible burning regimes [3, 4, 5,
6]. However, many technological challenges are involved with controlling
hydrogen combustion in lean turbulent conditions. In fact, flashback and
uncontrolled flame propagation that are typical of this regime [7, 8] are
amplified by the combination of hydrogen’s high reactivity and so its higher
flame speed [9], and its ability to auto-ignite [10, 11] and quickly diffuse
[12].
One interesting property of hydrogen flames that has not been fully understood
yet is its behaviour under strain. Previous studies have highlighted the
peculiar performance of hydrogen-enriched laminar flames with varying strain
rates, suggesting that high strain regimes can potentially be exploited for
hydrogen combustion practical applications. Hydrogen addition has been shown
to delay the extinction strain rate [3, 4], proving that hydrogen is able to
sustain very high strain rates. It has also been shown for syngas fuels that
hydrogen percentage influences the lean flame response with strain in terms of
flame temperature and NOx emissions [13], and flame temperature and
consumption speed [9]. In particular, in lean conditions the consumption speed
is shown to increase with strain, and this trend is opposite to the one
observed for hydrocarbon flames [14]. Furthermore, the hydrogen-enhanced
differential diffusion effects and the strain sensitivity have been proved to
be interdependent and to have a combined influence on the flame response [15,
16]. Even more interestingly, while the mass burning rate decreases with
strain for pure methane flames, a certain amount of hydrogen addition inverts
this trend, particularly in lean conditions [17]. Similar considerations hold
for the heat release rate, which is found to increase with strain in lean
hydrogen enriched flames [18]. However, the effect of these distinctive
hydrogen burning features at high strain rates on NOx emission for purely
hydrogen lean laminar flames is still an open question.
In this study we numerically investigate counter-flow premixed hydrogen flames
to shed light on their behaviour under intensive level of strain, in
particular the effect on NOx emissions. This is, to best of the authors’
knowledge, the first time that such an investigation is conducted. The
increase in heat release rate with strain discussed in the previous paragraph
would suggest a corresponding increase of temperature and so of NOx emissions.
Counter-intuitively, this study highlights that NOx emissions do not increase
with strain, and conversely show a decreasing trend in lean conditions,
particularly at very high strain rates. Nevertheless, despite the discussed
advantages and potentialities, additional unexplored control challenges are
potentially introduced at high strain regimes, such as complex vortex
dynamics, turbulence-flame interaction and local flame extinctions. Further
modelling challenges are associated with the prediction of these phenomena
distinctive of the regime.
As reported in Figure 1, there are two premixed strained flamelet
configurations documented in literature.

(a)

(b)

Figure 1: Sketch of the reactants-to-reactsants (a) and reactants-to-products
(b) counter-flow premixed strained flame configurations.

The first one is represented by a back-to-back or reactants-to-reactants
counter-flow configuration, where two flames are stabilized symmetrically with
respect to the flow stagnation plane (Figure 1(a)). Several studies are
available where the effects of hydrogen enrichment and strain on lean blow-
off, extinction strain, mass burning rate, and NOx emissions are investigated
with this setup [4, 3, 17, 13]. Specifically on NOx emissions trends, a
technical report from Xie and Wang [19] underlines a decreasing trend with
strain with this configuration in rich pure-hydrogen flames, attributing to
the NNH pathway the predominant contribution to NOx formation at high strain
levels. However, the limitation of this configuration is represented by the
proximity of the flames to the stagnation plane, particularly at the high
strain rates achievable by purely hydrogen flames. In these conditions, in
fact, the combustion has no space to complete, as it is clearly visible in
previous studies. Minor species are shown not be fully burnt within the
stagnation plane in Figure 5 and 6 of Jackson et al. [3]. Similarly, the
reaction rates of NOx pathways achieve zero within the stagnation plane for
lower strain rates in Ning et al. [13] (see Figures 6-8 of their work), but
not at higher strain rates. The same happens to the water rate of production
in Xie and Wang [19], Figure 2. Hence, this configuration appears less
suitable for investigations on emissions, and ultimately not practical for the
development of flamelet databases for turbulent LES simulations.
The second premixed strained flamelet configuration is represented by a
reactants-to-products counter-flow, where a single flame stabilizes on the
reactants side of the domain (Figure 1(b)). On the one hand, the single flame
setup allows the combustion reactions to complete even at very high strain
rates, as the fuel and the radicals have space to be burnt completely. This
evidence can be found for instance in Marzouk et al. [18], Figure 4-6. This
consideration suggests that the configuration is more appropriate for emission
analyses in highly-strained premixed flamelets. On the other hand, the
presence of complete combustion hot products on one of the boundaries can
precondition the problem, particularly considering the products temperature.
Despite this, the reactants-to-products configuration has been employed for
physical analysis of strained syngas and hydrogen-enriched flames [9, 18], and
widely discussed for methane flamelets [20, 21], particularly in the
application of LES-FGM models with strained flamelets [22, 23]. Therefore, the
reactants-to-products configuration is employed in the present study.
The purpose of the present work is to investigate lean hydrogen strained
flamelets, with a particular focus on NOx emissions trends. Detailed-chemistry
one-dimensional and two-dimensional DNS analyses are performed to achieve a
deeper understanding of the employed counter-flow strained flamelet
configuration. It will be shown for the first time that NOx emissions of
purely-hydrogen lean laminar flames display a decreasing trend with strain. A
detailed chemical analysis of NOx formation pathways is performed to support
this evidence, further showing that thermal NOx mechanism is predominant. This
physical behaviour of laminar flames is the first step to shed light on
potential features of novel combustor systems, where NOx emissions are
controlled by stabilising the flame against intensive strain.
This paper is organised as follows. The flamelet equations along with the
modelling choices performed for the 1D and 2D cases are introduced first,
followed by an overview of the 1D and 2D numerical setups. Emission results
are discussed next. Finally, relevant conclusions are drawn on the influence
of strain on NOx emissions.

## Model

Both one-dimensional and two-dimensional simulations are performed. For the
two-dimensional simulations, the reactingFoam solver in OpenFOAM [24] is
employed. The reacting Navier Stokes equations [25] are solved for mass,
momentum, absolute enthalpy and $N$ species with detailed chemistry. The
equation for the generic species $k$ is

$\frac{\partial(\rho Y_{k})}{\partial t}+\frac{\partial(\rho
u_{i}Y_{k})}{\partial x_{i}}=-\frac{\partial(\rho{V}_{k,i}Y_{k})}{\partial
x_{i}}+W_{k}\dot{\omega}_{k},$ (1)

where subscripts $i$ denotes direction $i$, $W_{k}$ is the molar mass, $\rho$
the mixture density, $V_{k,i}$ the diffusion velocity vector, and
$\dot{\omega}_{k}$ the molar rate of production of species $k$. A low Mach
number approximation is used in this study. Radiation, body forces, and
viscous dissipation effects are neglected. The Dufour effect on the heat flux
is also neglected.
The ideal gas law and the caloric equation of state are used as thermodynamic
model, where the species heat capacities are obtained using the JANAF
polynomials. Only laminar conditions are considered in this study, while the
effect of the turbulence-chemistry interaction will be investigated in a
future study. Detailed kinetic data of reactions from GRI3.0 mechanism [26]
are used to obtain consumption and production rate of species. In order to
speed up the simulations, the TDACChemistryModel is adopted, consisting of the
combination of the in situ adaptive tabulation (ISAT) algorithm with the
dynamic adaptive chemistry (DAC) reduction scheme [27]. A mixture-averaged
diffusion model is used [25] to account for the low Lewis number of the
hydrogen fuel as follows:

$V_{k,i}=-\frac{D_{k}^{M}}{X_{i}}\frac{\partial X_{k}}{\partial
x_{i}}\approx\frac{D_{k}^{M}}{Y_{i}}\frac{\partial Y_{k}}{\partial x_{i}}$
(2a) $D_{k}^{M}=\frac{1-Y_{k}}{\sum_{l,k\neq l}^{N}{X_{l}}/{D_{kl}}}$ (2b)

where the validity of the approximation performed in Eq. (2a) is verified in
post-processing. The binary diffusion coefficients $D_{kl}$ for the species
involved in the reactions are found with the Chapman-Enskog correlation [28,
29].
One-dimensional simulations are run with CHEM1D [30]. In one dimension and for
a flat reactants-to-products counter-flow flame, the set of conservation
equations solved in CHEM1D is re-arranged as follows [31]:

$\frac{\partial\rho}{\partial t}+\frac{\partial(\rho u_{x})}{\partial x}=-\rho
K$ (3a) $\frac{\partial(\rho Y_{k})}{\partial t}+\frac{\partial(\rho
u_{x}Y_{k})}{\partial x}=-\frac{\partial(\rho V_{x,k}Y_{k})}{\partial
x}+W_{k}\dot{\omega}_{k}-\rho KY_{k}$ (3b) $\frac{\partial\rho K}{\partial
t}+\frac{\partial\rho u_{x}K}{\partial x}=\frac{\partial}{\partial
x}\left(\mu\frac{\partial K}{\partial x}\right)+\rho_{\rm p}a^{2}-\rho K^{2}$
(3c) $\frac{\partial\rho h}{\partial t}+\frac{\partial\rho u_{x}h}{\partial
x}=\frac{\partial q}{\partial x}-\rho Kh$ (3d)

where the density of the products mixture $\rho_{\rm p}$, the applied strain
rate $a$, the local stretch rate $K$, along with Newton’s law for viscous
stresses have been introduced. The applied strain rate is a setup parameter
and is defined as the velocity gradient at the products boundary:

$a=-\left(\frac{du_{x}}{dx}\right)_{\rm p}.$ (4)

The influence of the $y$-component of the flow on the transport of scalars is
taken into account with the introduction of the local stretch rate [31]:

$K(x)=\frac{\partial u_{y}}{\partial y}.$ (5)

As a consequence of the continuity equation, the relation between the two
parameters reads $K(x\rightarrow\infty)=a$. Similarly to two-dimensional
simulations, GRI3.0 is used a chemical mechanism, along with a mixture-
averaged diffusion model (Eq. (2b)), while the binary diffusion coefficients
in CHEM1D are computed using molecular potentials and are tabulated as a
function of temperature in polynomial form [32].

## Computational setup

The combustion of hydrogen is evaluated at atmospheric conditions. Temperature
and species boundary conditions at the reactants and products boundary are
prescribed as follows for both 1D and 2D cases. An equivalence ratio of
$\phi=0.7$ is imposed for both streams. The reactants temperature is $T_{\rm
r}=300\,{\rm K}$, while the products temperature is set to the adiabatic flame
temperature of an unstrained hydrogen flamelet with the same equivalence ratio
computed with CHEM1D, $T_{\rm p}=2021\,{\rm K}$. Mass fractions at the
products boundary are imposed from complete combustion. A summary of the
temperature and species boundary conditions for both 1D and 2D simulations are
reported in Table 1.

Table 1: Temperature and species boundary conditions for both 1D and 2D simulations. Quantity | Left boundary (reactants) | Right boundary (products)
---|---|---
T [K] | 300 | 2021
$Y_{\rm H_{2}}$ [-] | 0.02 | 0
$Y_{\rm O_{2}}$ [-] | 0.228 | 0.068
$Y_{\rm H_{2}O}$ [-] | 0 | 0.18
$Y_{\rm N_{2}}$ [-] | 0.752 | 0.752

### 1D setup

For the one-dimensional simulations, a wide computational domain of $L_{\rm
1D}=20$ cm is chosen. CHEM1D uses an adaptive mesh algorithm, thus ensuring
that a proper mesh refinement in the flame region. The adaptive grid accounts
for 200 points in total. An exponential differential scheme is used for the
spatial discretization, and a second-order time integration of the
differential equations is performed by the stationary solver. The time step is
adjusted automatically by the numerical tool to achieve convergence, and
ranges between $10^{-6}\,{\rm s}$ and $10^{-8}\,{\rm s}$. The only input
parameter for the velocity field is represented by the applied strain rate
$a$, defined in Eq. (4). According to the applied strain rate definition, the
higher is $a$, the higher will the boundary velocities be and the stretch rate
experienced by the flame. A broad range of applied strain rates is
investigated for 1D simulations, from $100\,{\rm s}^{-1}$ up to $10000\,{\rm
s}^{-1}$ to eventually observe the occurrence of flame extinction. The
simulations wall clock time ranges from a few minutes for very high strain
rates up to a few hours for lower strain rates.

### 2D setup

The PISO time integration solver [33], with an implicit second-order Euler-
backward discretization scheme, is used for the 2D simulations. Similarly, a
second-order central scheme is used for the convective term of all resolved
quantities. A variable time step is used ensuring a maximum Courant number of
0.5 for all domain. The simulations are run for a physical time of
$t=0.5\,{\rm s}$ for both the applied strain rates investigated, corresponding
to more than 100 times the maximum estimated flow-through time inside the
domain. A typical computation required 225 and 895 CPU-hours respectively for
lower and higher strain rate cases.
The higher computational cost of two-dimensional DNS simulations led to the
choice of decreasing the dimension of the domain to $L_{\rm 2D}=2$ cm in both
the horizontal and the vertical direction. Nevertheless, this domain was
sufficient to observe a flame position enough far from the boundary for all
simulations. The parent mesh is made of 100 uniform quadrilaterals in the
x-direction and 40 in the y-direction, resulting into $\Delta x=0.2$ mm and
$\Delta y=0.5$ mm. Hence, the quads are refined progressively by a factor of 2
for three times in the region in the proximity of the flame. An example for
$a=5000\,{\rm s}^{-1}$ is provided in Figure 2. The flame position is
determined for every applied strain rate combining preliminary coarse 2D
simulations and 1D simulations, and is then double-checked in the refined
simulation solution a posteriori. After the refinement, the minimum spacing
results into $\Delta x_{\rm min}=0.2/2^{3}=0.025$ mm, which ensures at least
12 cells are present within the flame thickness. The number of cells enclosed
within the flame thickness is similar to that found in previous DNS studies on
reacting flows [34, 35, 36].

Figure 2: Sketch of the two-dimensional setup with mesh refinement for
$a=5000\,{\rm s}^{-1}$.

In contrast to CHEM1D, where the applied strain rate $a$ was the only input
parameter for the velocity field, a uniform velocity at the two inlets is
prescribed in OpenFOAM’s setup. Given the extended-domain one-dimensional
simulation at a given applied strain rate, the following horizontal velocity
boundary conditions are thus applied in the 2D case:

$\left(U_{\rm r,2D}\right)_{a=\rm const}=\left(U_{\rm
1D}(x=-1\,\rm{cm})\right)_{a=\rm const},$ (6a) $\left(U_{\rm
p,2D}\right)_{a=\rm const}=\left(U_{\rm 1D}(x=1\,\rm{cm})\right)_{a=\rm
const}.$ (6b)

Hence, at the upper and lower outlets, a zero gradient boundary condition is
prescribed. The discrepancy in the inlet boundary conditions definition
precludes the possibility to directly compare 1D and 2D simulations at a given
applied strain rate, because the resulting 2D horizontal velocity profile does
not perfectly overlap with the one in the 1D solution. However, as it will be
further discussed in the Results section, the purpose of two-dimensional
simulations is solely to show that the same NOx emission trends with strain
are observed with respect to 1D simulations. Therefore, this setup discrepancy
is accepted for the objectives of the present study. Furthermore, only two
high applied strain rates configurations are investigated in the 2D case,
$a=2000\,{\rm s}^{-1}$ and $a=5000\,{\rm s}^{-1}$.

## Results

### NOx emissions trends

NOx generally refers to nitrogen oxides, including both NO and NO2. However,
NO2 mass fraction peaks are shown to be at least three orders of magnitude
smaller than the ones of NO for any setup investigated. Therefore, similarly
to previous studies [13], NO emissions can alone be considered representative
of the overall NOx emissions and are thus examined for the purposes of the
present study. The behaviour of $Y_{\rm NO}$ in the longitudinal direction at
different applied strain rates predicted by CHEM1D in the proximity of the
stagnation plane is shown in Figure 3. It is immediately visible that both the
peak and the area under the curve decrease consistently as the applied strain
is increased. In addition, no extinction is observed at the highest strain
rates achieved. As one could expect considering previous studies involving
hydrogen-enriched fuels [3, 4], this evidence confirms that pure hydrogen
flames can potentially sustain very high strain levels. In Figure 4(a) the
decrease of the peaks of $Y_{\rm NO}$ in one-dimensional and the center line
of two-dimensional simulations is further highlighted. The discrepancies in
the exact values of $Y_{\rm NO}$ peaks at the same applied strain can be
pointed to the different definition of the inlet boundary conditions discussed
in the 2D setup section. Nevertheless, the graph shows that two-dimensional
simulations already confirm the trend of decreasing NOx emissions with
increasing strain rate. Similar considerations hold for the trend of NO flux
reported in the left y-axis of Figure 4(b), which is calculated as

$\Phi_{\rm NO,1D}=\int_{-L/2}^{L/2}\rho Y_{\rm NO}dx.$ (7)

Figure 3: NO mass fraction as a function space for different applied strain
rates.

(a)

(b)

Figure 4: Variation of NO mass fraction (a) and NO flux (b) with strain.

One might object that NO is a slow species which would form downstream of the
flame in the hot products region, and this formation in the products side is
prevented by the right boundary condition, where $Y_{\rm NO}(x=L/2)=0$ is
imposed. In fact, in the present configuration $Y_{\rm NO}$ peaks at the
stagnation plane to then suddenly drop to fulfill the boundary condition, as
clearly visible in Figure 3. Nevertheless, the investigation of the 2D
simulations, where the NO is free to form in the $y$-direction (tangential to
the flame in Figure 1), suggests this not to be the case, i.e. NO is
suppressed at high strain by some other means. Hence, the focus is shifted to
the understanding of the emission phenomena in the flame-tangential direction
of the two-dimensional simulations. One can compute the two-dimensional NO
flux for the two setups investigated by integrating on the 2D-domain surface
as follows:

$\Phi_{\rm NO,2D}=\int_{A}\rho Y_{\rm NO}dA.$ (8)

This flux is reported on the right y-axis in Figure 4(b). It can be
immediately observed that the same decreasing NO flux trend is observed for
the 1D and 2D cases. Furthermore, from the analysis of the 2D data, it is
highlighted that

$\Phi_{\rm NO,2D}\approx\Phi_{\rm NO,1D}\cdot L_{\rm 2D},$ (9)

i.e. the 1D flux of NO calculated across the center line multiplied by the
vertical dimension is approximately equal to the 2D flux. These considerations
raise confidence on the fact that the decrease of NO emissions observed in 1D
simulations and in the 2D domain centreline as strain is increased is not
compensated by any additional NO formation in the vertical direction. This
evidence can be further proved following the same streamline at the two
different strain rate cases investigated with OpenFOAM. The streamline
investigated originates at $y=-0.2$ mm, such that it crosses the flame and
thus travels in the products and exits the domain without clashing with the
stagnation plane. As shown in Figure 5(b), the particle residence time inside
the $a=5000\,{\rm s}^{-1}$ domain is estimated around $\Delta
t=1.23\,\rm{ms}$. The path that the particle originating in the same position
has travelled in the same time interval in the $a=2000\,{\rm s}^{-1}$ domain
is reported in Figure 5(a). Hence, a domain cut is taken at the $y$-location
the particle has reached at $\Delta t$ in the two cases (green line in Figure
5), and the NO fluxes computed along the line cuts are reported in Table 2.

(a) $a=2000\,{\rm s}^{-1}$.

(b) $a=5000\,{\rm s}^{-1}$.

Figure 5: Travelling path of the a streamline originating at $y=-0.2$ mm on a contour plot of the progress variable $c$ based on $Y_{\rm H_{2}O}$. Table 2: Flux of NO after the same particle has travelled for $\Delta t=1.23$ ms at two different applied strain rates. $a$ [$\,{\rm s}^{-1}$] |  $y(\Delta t)$ [mm] |  $\Phi(y)$ [g/cm2]
---|---|---
2000 | -0.53 | 9.51 $\cdot$ 10-10
5000 | -10 | 2.56 $\cdot$ 10-10

Not only do the data show a still lower NO flux at higher strain rates, but it
can be also observed that the values are very close to the ones at the center
line (reported in the graph in Figure 4(b)). Therefore, it is further proved
that there is no additional NO formation in the vertical direction at higher
strain rates. It can be concluded that, at least within the studied framework,
the observed NO emissions decrease with increasing strain is not a consequence
of the specific counter-flow configuration chosen and of the boundary-
preconditioned numerical setup, but is instead a direct effect of the
increased tangential velocity gradients on the hydrogen flame combustion and
dissociation reactions.

### NO formation pathways

The focus is now shifted to the evaluation of the NO formation pathways, with
the goal of further understanding the physical reasons behind the observed NO
decrease with strain. As carbon species are not involved in hydrogen
combustion, the prompt NOx formation pathway is not considered, and the only
three pathways evaluated are thermal NO, NNH-NO and N2O-NO. CHEM1D solutions
data are used for the purposes of the present NO formation routes
investigations, since the trends of NO emissions with strain predicted by one-
dimensional and two-dimensional simulations have been shown to compare with
good approximation (see Figure 4(b)). Two high applied strain rates setups of
$a=2000\,{\rm s}^{-1}$ and $a=5000\,{\rm s}^{-1}$ are analysed. For a sample
reaction $r$ involved in a pathway
$\nu_{A}A+\nu_{B}B\leftrightarrow\nu_{C}C+\nu_{D}D$, the forward reaction
coefficient is obtained with Arrenhius law, and thus the reverse reaction
coefficient can be found [25]:

$K_{f,r}(T)=A_{r}T^{\beta_{r}}e^{\left(\frac{-E_{a,r}}{R_{0}T}\right)},\hskip
28.45274ptK_{r,r}(T)=\frac{K_{f,r}}{\left(\frac{p_{a}}{R_{0}T}\right)^{\sum_{s}\nu_{s}}e^{\left(\frac{\Delta
S_{r}^{0}}{R_{0}}-\frac{\Delta H_{r}^{0}}{R_{0}T}\right)}},$ (10)

where $A_{r}$, $\beta_{r}$ and $E_{a,r}$ are found in the GRI3.0 mechanism
reaction web page [26], and $\Delta S_{r}^{0}$ and $\Delta H_{r}^{0}$ are
respectively the reaction enthropy and enthalpy that can be found as a
function of temperature with JANAF polynomials. Thus, the forward and reverse
reaction rates are found as a function of $T(x)$ introducing respectively the
reactants and products concentrations:

$\dot{\omega}_{f,r}=K_{f,r}(x)\cdot[A(x)]^{\nu_{A}}[B(x)]^{\nu_{B}},\hskip
28.45274pt\dot{\omega}_{r,r}=K_{r,r}(x)\cdot[C(x)]^{\nu_{C}}[D(x)]^{\nu_{D}}.$
(11)

Three-body and pressure-dependent reactions have been treated with dedicated
formulas [37]. Hence, following the formulation of Ning et al. [13], the
overall reaction rate ($ORR$) of the reaction $r$ is obtained by integrating
the reaction rate across the longitudinal (flame-normal) direction:

$ORR_{r}=\int_{-L/2}^{L/2}\dot{\omega}_{r}(x)dx.$ (12)

Summing the contribution of the single reactions $ORR$ along a specific route,
a chart diagram showing the NO formation pathways at $a=2000\,{\rm s}^{-1}$
and $a=5000\,{\rm s}^{-1}$ is obtained and is reported in Figure 6(a). While
all the reaction rates involved in the pathways show a decrease at higher
strain, it can already be observed that this decrease is remarkably more
consistent considering the reactions involved in the thermal NO pathway.
Hence, the $ORR$ of the reactions directly producing NO in each pathway are
summed and the bar chart in Figure 6(b) is obtained. The histogram shows that
the main NO production pathways involved in lean hydrogen flames at high
strain rates are both the thermal and the NNH. Considering syngas fuel
mixtures and twin counter-flow configuration, Ning et al. [13] have shown with
a similar analysis that the NNH pathway was by itself predominant, even with a
high hydrogen content. The reason of the discrepancy with the present study is
probably related to the lower equivalence ratio of their investigation
($\phi=0.5$), which may decrease the maximum flame temperature and thus the
weight of thermal NO pathway on the total NO formation process. In addition,
they conclude that with a high hydrogen content in the fuel and high strain
rates the NO production is exacerbated. Once again, the fuel mixture and the
equivalence ratio, along with the maximum strain rates achieved and the
counter-flow configuration chosen do not allow for a direct comparison even in
the trends with the present study, which instead investigates on completely
novel setup with 100% hydrogen fuel and very high strain rates.

(a)

(b)

Figure 6: (a) Chart diagram showing the overall reaction rates
[mol/(cm${}^{2}\cdot$s)] along the NO formation pathways. Thermal NO pathway
is highlighted in red, NNH-NO pathway in green, NNH-NO-HNO intermediate
pathway in black, and N2O-NO pathway in blue. (b) Contribution of each NO
formation pathway to the total NO formation. For strain level $a=2000\,{\rm
s}^{-1}$ and $a=5000\,{\rm s}^{-1}$.

Figure 6(b) confirms that the main decrease of NO production at higher strain
rates is associated to the thermal pathway. Speth and Ghoniem [9] have shown
that for syngas with a high content of H2, the flame temperature has an
increasing trend with strain within their limited strain range, i.e. up to
$a=500\,{\rm s}^{-1}$. Furthermore, extending the strain range would probably
result in a non-monotonic function, as their study shows for fuel mixtures
with less hydrogen content. NOx emissions can be expected to follow directly
the flame temperature function. In contrast, the present study shows that NOx
emissions decrease monotonically in a very broad range of strain rates, thus
indicating that the emissions function does not follow directly the
temperature function. Hence, the reason behind the observed NO decreasing
trend may be not only imputable to a change in the flame temperature. This
consideration suggests that a more extensive analysis should be performed on
the local temperature and concentration of the radicals involved in the
pathway reactions in order to further shed light on the physical explanation
behind the phenomenon of decreasing NOx with strain in lean premixed hydrogen
flames.

## Conclusion

Detailed-chemistry one-dimensional and two-dimensional DNS analyses are
conducted on pure hydrogen lean premixed and strained flamelets in a counter-
flow reactants-to-products configuration. A broad range of strain rates is
investigated, up to very high levels which have been rarely considered before
in the literature, and which hydrogen is shown to be able to sustain without
experiencing extinctions. Both 1D and 2D simulations show that NOx emissions
have a decreasing trend as the strain rate is increased. This is shown for the
first time for pure hydrogen fuel and lean flamelets. The analysis of the 2D
results also indicates that there is no compensating NOx formation in the
second dimension, and thus that the trend observed is not a consequence of any
possible setup preconditioning. Therefore, the greater tangential velocity
gradients may have a direct effect on the flame, somehow affecting combustion
efficiency and consequently the rate of reactions through which NOx is formed.
In-depth analyses of the NOx formation mechanisms and reaction rates further
underline that the main contribution to the decreasing trend arises from the
thermal NOx pathway. Future work will focus on a more extensive local
temperature and radical concentration study on the setup investigated and at
further conditions in order to achieve a deeper understanding of the physics
behind the phenomenon. Consumption speed trend with strain and the role of
differential diffusion effects are also hydrogen-peculiar features that may
play a role and that should be further explored. If the emission trend with
strain are confirmed in turbulent conditions, the low NOx outcomes of lean
premixed and highly-strained hydrogen flames can be exploited in future
technological applications of hydrogen combustion systems.

## Acknowledgments

AP and IL acknowledge the Dutch Ministry of Education and Science for
providing funding support to this project via the Sector Plan scheme. BK and
IL further acknowledge support from the Engineering and Physical Sciences
Research Council, grant n. EP/T028084/1.

## References

  * Turner [2004] John A Turner.  Sustainable hydrogen production.  _Science_ , 305(5686):972–974, 2004.
  * Jensen et al. [2007] Søren H Jensen, Peter H Larsen, and Mogens Mogensen.  Hydrogen and synthetic fuel production from renewable energy sources.  _International Journal of Hydrogen Energy_ , 32(15):3253–3257, 2007.
  * Jackson et al. [2003] Gregory S Jackson, Roxanne Sai, Joseph M Plaia, Christopher M Boggs, and Kenneth T Kiger.  Influence of h2 on the response of lean premixed ch4 flames to high strained flows.  _Combustion and Flame_ , 132(3):503–511, 2003\.
  * Cho and Chung [2009] Eun-Seong Cho and Suk Ho Chung.  Improvement of flame stability and nox reduction in hydrogen-added ultra lean premixed combustion.  _Journal of Mechanical Science and Technology_ , 23(3):650–658, 2009.
  * Barbosa et al. [2007] Séverine Barbosa, Marta de La Cruz Garcia, Sébastien Ducruix, Bernard Labegorre, and François Lacas.  Control of combustion instabilities by local injection of hydrogen.  _Proceedings of the combustion Institute_ , 31(2):3207–3214, 2007.
  * Lipatnikov and Chomiak [2005] AN Lipatnikov and Jerzy Chomiak.  Molecular transport effects on turbulent flame propagation and structure.  _Progress in energy and combustion science_ , 31(1):1–73, 2005.
  * Huang and Yang [2009] Ying Huang and Vigor Yang.  Dynamics and stability of lean-premixed swirl-stabilized combustion.  _Progress in energy and combustion science_ , 35(4):293–364, 2009.
  * Ducruix et al. [2003] Sebastien Ducruix, Thierry Schuller, Daniel Durox, and Sébastien Candel.  Combustion dynamics and instabilities: Elementary coupling and driving mechanisms.  _Journal of propulsion and power_ , 19(5):722–734, 2003.
  * Speth and Ghoniem [2009] Raymond L Speth and Ahmed F Ghoniem.  Using a strained flame model to collapse dynamic mode data in a swirl-stabilized syngas combustor.  _Proceedings of the Combustion Institute_ , 32(2):2993–3000, 2009.
  * Hu and Kurose [2019] Yong Hu and Ryoichi Kurose.  Large-eddy simulation of turbulent autoigniting hydrogen lifted jet flame with a multi-regime flamelet approach.  _International Journal of Hydrogen Energy_ , 44(12):6313–6324, 2019.
  * Johannessen et al. [2015] Birgitte Johannessen, Andrew North, Robert Dibble, and Terese Løvås.  Experimental studies of autoignition events in unsteady hydrogen–air flames.  _Combustion and Flame_ , 162(9):3210–3219, 2015\.
  * Dinesh et al. [2016] KKJ Ranga Dinesh, H Shalaby, KH Luo, JA Van Oijen, and D Thévenin.  High hydrogen content syngas fuel burning in lean premixed spherical flames at elevated pressures: Effects of preferential diffusion.  _International Journal of Hydrogen Energy_ , 41(40):18231–18249, 2016.
  * Ning et al. [2017] Daoguan Ning, Aiwu Fan, and Hong Yao.  Effects of fuel composition and strain rate on no emission of premixed counter-flow h2/co/air flames.  _international journal of hydrogen energy_ , 42(15):10466–10474, 2017.
  * Liang et al. [2017] Wenkai Liang, Fujia Wu, and Chung K Law.  Extrapolation of laminar flame speeds from stretched flames: Role of finite flame thickness.  _Proceedings of the Combustion Institute_ , 36(1):1137–1143, 2017.
  * Abbasi-Atibeh and Bergthorson [2019] Ehsan Abbasi-Atibeh and Jeffrey M Bergthorson.  Differential diffusion effects in counter-flow premixed hydrogen-enriched methane and propane flames.  _Proceedings of the Combustion Institute_ , 37(2):2399–2406, 2019.
  * Vance et al. [2022] Faizan H Vance, Y Shoshin, LPH de Goey, and Jeroen A van Oijen.  Quantifying the impact of heat loss, stretch and preferential diffusion effects to the anchoring of bluff body stabilized premixed flames.  _Combustion and Flame_ , 237:111729, 2022.
  * Van Oijen et al. [2016] JA Van Oijen, A Donini, RJM Bastiaans, JHM ten Thije Boonkkamp, and LPH De Goey.  State-of-the-art in premixed combustion modeling using flamelet generated manifolds.  _Progress in Energy and Combustion Science_ , 57:30–74, 2016.
  * Marzouk et al. [2005] Youssef M Marzouk, Raymond L Speth, and Ahmed F Ghoniem.  Combined effects of curvature and strain on hydrogen enriched lean methane flames.  Technical report, Sandia National Lab.(SNL-CA), Livermore, CA (United States), 2005.
  * Xie and Wang [2013] Tian Fang Xie and Pei Yong Wang.  Analysis of no formation in counterflow premixed hydrogen-air flame.  In _Applied Mechanics and Materials_ , volume 284, pages 601–607. Trans Tech Publ, 2013.
  * Libby and Williams [1982] Paul A Libby and Forman A Williams.  Structure of laminar flamelets in premixed turbulent flames.  _Combustion and flame_ , 44(1-3):287–303, 1982\.
  * Darabiha et al. [1986] N Darabiha, SM Candel, and FE Marble.  The effect of strain rate on a premixed laminar flame.  _Combustion and flame_ , 64(2):203–217, 1986\.
  * Kolla and Swaminathan [2010] H Kolla and N Swaminathan.  Strained flamelets for turbulent premixed flames, i: Formulation and planar flame results.  _Combustion and Flame_ , 157(5):943–954, 2010\.
  * Langella and Swaminathan [2016] Ivan Langella and Nedunchezhian Swaminathan.  Unstrained and strained flamelets for les of premixed combustion.  _Combustion Theory and Modelling_ , 20(3):410–440, 2016.
  * Weller et al. [1998] Henry G Weller, Gavin Tabor, Hrvoje Jasak, and Christer Fureby.  A tensorial approach to computational continuum mechanics using object-oriented techniques.  _Computers in physics_ , 12(6):620–631, 1998\.
  * Poinsot and Veynante [2005] Thierry Poinsot and Denis Veynante.  _Theoretical and numerical combustion_.  RT Edwards, Inc., 2005.
  * Gregory P. Smith, David M. Golden, Michael Frenklach, Nigel W. Moriarty, Boris Eiteneer, Mikhail Goldenberg, C. Thomas Bowman, Ronald K. Hanson, Soonho Song, William C. Gardiner, Jr., Vitali V. Lissianski, and Zhiwei Qin [2000] Gregory P. Smith, David M. Golden, Michael Frenklach, Nigel W. Moriarty, Boris Eiteneer, Mikhail Goldenberg, C. Thomas Bowman, Ronald K. Hanson, Soonho Song, William C. Gardiner, Jr., Vitali V. Lissianski, and Zhiwei Qin.  GRI-MECH 3.0, 2000.  URL http://combustion.berkeley.edu/gri-mech/version30/text30.html.
  * Contino et al. [2011] Francesco Contino, Hervé Jeanmart, Tommaso Lucchini, and Gianluca D’Errico.  Coupling of in situ adaptive tabulation and dynamic adaptive chemistry: An effective method for solving combustion in engine simulations.  _Proceedings of the Combustion Institute_ , 33:3057–3064, 1 2011.  ISSN 1540-7489.  doi: 10.1016/J.PROCI.2010.08.002.
  * Chapman and Cowling [1990] Sydney Chapman and Thomas George Cowling.  _The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases_.  Cambridge university press, 1990.
  * Novaresio et al. [2012] Valerio Novaresio, María García-Camprubí, Salvador Izquierdo, Pietro Asinari, and Norberto Fueyo.  An open-source library for the numerical modeling of mass-transfer in solid oxide fuel cells.  _Computer Physics Communications_ , 183:125–146, 1 2012\.  ISSN 00104655.  doi: 10.1016/j.cpc.2011.08.003.
  * Eindhoven University of Technology [2021] Eindhoven University of Technology.  CHEM1D. A one dimensional laminar flame code, 2021.  URL https://github.com/thijsa93400/TUe_chem1d.
  * Ramaekers [2011] WJS Ramaekers.  _Development of flamelet generated manifolds for partially-premixed flame simulations_.  PhD thesis, Technische Universiteit Eindhoven, 2011.
  * Somers [1994] Bart Somers.  _The simulation of flat flames with detailed and reduced chemical models_.  PhD thesis, Gastec NV NL, 1994.
  * Issa [1986] Raad I Issa.  Solution of the implicitly discretised fluid flow equations by operator-splitting.  _Journal of computational physics_ , 62(1):40–65, 1986.
  * Klein et al. [2020] M Klein, A Herbert, H Kosaka, B Böhm, A Dreizler, N Chakraborty, V Papapostolou, Hong G Im, and J Hasslberger.  Evaluation of flame area based on detailed chemistry dns of premixed turbulent hydrogen-air flames in different regimes of combustion.  _Flow, Turbulence and Combustion_ , 104(2):403–419, 2020.
  * Lee et al. [2022] HsuChew Lee, Peng Dai, Minping Wan, and Andrei N Lipatnikov.  A dns study of extreme and leading points in lean hydrogen-air turbulent flames-part ii: local velocity field and flame topology.  _Combustion and Flame_ , 235:111712, 2022.
  * Berger et al. [2020] Lukas Berger, Raik Hesse, Konstantin Kleinheinz, Michael J Hegetschweiler, Antonio Attili, Joachim Beeckmann, Gregory T Linteris, and Heinz Pitsch.  A dns study of the impact of gravity on spherically expanding laminar premixed flames.  _Combustion and Flame_ , 216:412–425, 2020.
  * Kee et al. [1989] Robert J Kee, Fran M Rupley, and James A Miller.  Chemkin-ii: A fortran chemical kinetics package for the analysis of gas-phase chemical kinetics.  Technical report, Sandia National Lab.(SNL-CA), Livermore, CA (United States), 9 1989.  URL https://www.osti.gov/biblio/5681118.