Learning Flame Evolution Operator under Hybrid Darrieus Landau and Diffusive Thermal Instability

Recent advancements in the integration of artificial intelligence (AI) and machine learning (ML) with physical sciences have led to significant progress in addressing complex phenomena governed by nonlinear partial differential equations (PDE). This paper explores the application of novel operator learning methodologies to unravel the intricate dynamics of flame instability, particularly focusing on hybrid instabilities arising from the coexistence of Darrieus-Landau (DL) and Diffusive-Thermal (DT) mechanisms. Training datasets encompass a wide range of parameter configurations, enabling the learning of parametric solution advancement operators using techniques such as parametric Fourier Neural Operator (pFNO), and parametric convolutional neural networks (pCNN). Results demonstrate the efficacy of these methods in accurately predicting short-term and long-term flame evolution across diverse parameter regimes, capturing the characteristic behaviors of pure and blended instabilities. Comparative analyses reveal pFNO as the most accurate model for learning short-term solutions, while all models exhibit robust performance in capturing the nuanced dynamics of flame evolution. This research contributes to the development of robust modeling frameworks for understanding and controlling complex physical processes governed by nonlinear PDE.


I. INTRODUCTION
In recent years, the integration of artificial intelligence (AI) and machine learning (ML)   with the natural sciences and physical engineering has led to significant advancements, particularly in addressing the complexities of nonlinear partial differential equations (PDE).
These equations are fundamental in understanding various physical phenomena, ranging from turbulent fluid dynamics to complicate physico-chemical processes.Within the domain of nonlinear PDE systems lies a rich tapestry of intricate dynamics, including instabilities, multiscale interactions, and chaotic behaviors.To enhance predictive capabilities and design robust control strategies in engineering applications, computational methods are indispensable.These methods, often in the form of numerical solvers, enable the accurate simulation of PDE solutions across spatial and temporal domains.Implicit in these solvers is the concept of the functional-mapping operator, which could iteratively advances the PDE solution functions in time, providing a pathway to explore the evolution of physical systems over extended durations.A distinctive class of machine learning methods has emerged, capable of learning and replicating the behavior of these PDE operators.
These CNN-based approaches parameterize the PDE operator in a finite-dimensional space, enabling the mapping of discrete functions onto image-like representations.Building upon this foundation, recent strides have witnessed the development of neural operator methods 8,9 capable of learning operators of infinite dimensionality.Notable examples include Deep Operator Network 10 and the Fourier Neural Operator (FNO) 11 , both demonstrating remarkable proficiency across a diverse array of benchmark problems 12,13 .Furthermore, recent advancements have extended neural operators by amalgamating concepts from wavelet methods 14,15 and adapting approaches for complex domains 16 .
In our recent investigations 17,18 , we delved into the intricate dynamics of flame instability and nonlinear evolution, a canonical problem with profound implications for combustion science.Flames can undergo destabilization due to intrinsic instabilities, including the hydrodynamic Darrieus-Laudau (DL) mechanism 19,20 attributed to density gradients across a flame, and the Diffusive-Thermal (DT) mechanism 21,22 driven by heat and reactant diffusion disparities.Our previous work 17 primarily focused on DL flames, scrutinizing the evolution of unstable flame fronts within periodic channels of varying widths.Under DL instability, an initially planar flame morphs into a steady curved front; as the channel width increases, the curved front becomes sensitive to random noise, and small wrinkles start to emerge.At sufficiently large channels, DL flames give rise to complicated fractal fronts characterized by hierarchical cascading of cellular structures 23 .
The nonlinear evolution of DT flame development can be modeled by the Michelson-Sivashinsky equation 24 , while a more accurate but computationally expensive approach involves direct numerical simulation (DNS) of Navier-Stokes equations.Utilizing these two approaches to generate training datasets, our investigations 17 demonstrated that both CNN and FNO could effectively capture the evolution of DL flames, with FNO exhibiting superior performance in modeling complex flame geometries over longer durations.Subsequently, we embarked on developing parameterized learning methodologies capable of encapsulating dynamics across diverse parameter regimes within a single network framework.Through the introduction of pCNN and pFNO models 18 , we demonstrated their efficacy in replicating the behavior of DL flames across varying channel widths.Additionally, our methods have shown success in learning the parametric solutions of the Kuramoto-Sivashinsky equation 25 , which models unstable flame evolution due to the DT mechanism.However, a challenge remains in mitigating the tendency of these models to overestimate noise effects.
In this paper, we extend our research horizon to encompass the complexities arising from hybrid instabilities, specifically those arising from the coexistence of DL and DT mechanisms.These hybrid systems pose new challenges, as they embody a rich spectrum of behaviors stemming from the interplay of distinct instability modes.Leveraging our novel operator learning methodologies, we aim to unravel the nuanced dynamics underlying such hybrid instabilities, shedding light on their short-term evolution and long-term statistical properties.Furthermore, our endeavor holds promise for the development of robust modeling frameworks capable of capturing the intricate dynamics of real-world flame evolution scenarios.
The paper is organized as follows: first, we describe the problem setup for learning PDE operators, followed by brief descriptions of the two parametric learning methods to be used in this work.These methods will be compared in the context of learning parametric-dependent solution time-advance operators for the Sivashinsky Equation 26, which models unstable front evolution due to hybrid mechanisms of flame instability.Finally, we provide a summary and conclusion.

II. PROBLEM SETUP FOR LEARNING PDE OPERATORS
In this section, we delineate the problem setup for learning a parametric PDE operator, along with a description of recurrent training methods.
Consider a system governed by PDE, typically involving multiple functions and mappings between them.Our focus here is on a parametric operator mapping, denoted as where γ ∈ R dγ represents a set of parameters.The input function v(x), where x ∈ D, resides in a functional space V(D, R dv ), with domain D ⊂ R d and codomain R dv , while the output Our primary interest lies in the solution time advancement operator with parametric dependence, given by Ĝ : (ϕ(x, t), γ) → ϕ(x; t + 1) where ϕ(x; t) denotes the solution to a PDE under parameters γ, and t = t/∆ t represents normalized time with a small time increment ∆ t .For simplicity, we assume identical domain and codomain for both input and output functions, i.e., D ′ = D, V ′ = V, d ′ = d, and d ′ v = d v , with periodic boundary conditions on D.
To approximate the mapping Ĝ using neural network methods, let Θ denote the space of trainable parameters in the network.A neural network can be defined as where Θ represents the space of network parameters.Training the neural network involves finding an optimal choice of parameters θ * ∈ Θ such that G θ * approximates Ĝ.
Starting with an initial solution function ϕ(x; t 0 ) under fixed parameter values γ, the recurrent application of the operator G θ,γ := G θ (•, γ) can roll out predicted solutions of arbitrary length by iteratively updating the input function with its output from the previous prediction.Note, while the learned operator is expected to make accurate short-term predictions; its long-term prediction might be allowed to deviate if the ground truth PDE admits chaotic solutions.On the other hand, it is still desirable that the learned operator can reproduce the correct statistics in the long-term solutions.
Following previous studies 17,18 , our training approach adopts a one-to-many setup where the recurrent network is trained to make multiple successive predictions from a single input function.Such setup ensures numerical stability in the learned solution advancement operator, a crucial consideration highlighted in the prior work 17,18 .More specifically, let where v ∼ χ ′ and γ ∼ χ ′′ are randomly drawn according to independent probability measures of χ ′ and χ ′′ respectively.The cost function C :

III. PARAMETRIC OPERATOR LEARNING METHODS
In this section, we present a concise overview of two methods capable of learning the parametric operator Ĝ.Further details about these methods can be found in paper 18 .

A. Parametric Convolutional Neural Network (pCNN)
The operator G θ,γ can be regarded as an image-to-image map when applied for the temporal advancement of discretized solutions.Deep Convolutional Neural Networks (CNN) have demonstrated effectiveness in image learning tasks.The network architecture suitable for learning operators resembles a convolutional auto-encoder similar to that in U-Net 6 and ConvPDE 7 .This network comprises an encoder block and a decoder block, with input data passing through a series of transformations via convolutional layers.Additionally, the method incorporates side networks to handle additional parameter inputs.The pCNN model is outlined in Fig. 1.The output data channels c l for each convolution layer are indicated within brackets.
Let e + 0 denote the input function v(x j ) represented on an x-mesh.The encoder block follows an iterative update procedure: e + l → (e l+1 , e * l+1 ) → e + l+1 .This iteration occurs over the level sequence l = 0, 1, ..., L − 1. Denote the last decoding output as e ′ L = e + L , a subsequent decoding procedure is applied (e ′ l+1 , e + l ) → e ′ l through reversing level l.
Here, e l , e * l , e + l , e ′ l ∈ R c l ×N l represent four data sequences, each with c l channels and size of N l .The data size is halved as N l+1 = N l /2 for l ≥ 1.The first-stage encoder contains two sub-maps of e + l → e l+1 and e + l → e * l+1 , both are implemented using vanilla-stacked convolution layers (with a filter size of 3, stride 1, periodic padding, and ReLU activation).Some layers are replaced by Inception layers for improved performance.Additionally, a size-2 max-pooling layer is prepended to halve the image size for l ≥ 1.The second-stage encoder map is implemented as e + l+1 = e l+1 + e * l+1 • D l (γ).Here, D l is a simple function (a two-layer perceptron) that converts the PDE parameters γ into a scaling ratio.The decoder update (e ′ l+1 , e + l ) → e ′ l involves concatenating e ′ l+1 (but up-sample it to double its size) with e + l along the channel dimension.The final output is obtained as v ′ (x j ) = e ′ 1 .

B. Parametric Fourier Neural Operator (pFNO)
The parametric Fourier Neural Operator (pFNO) 18 was developed based on the original FNO method 11 , wherein learning for the infinite-dimensional operator is achieved by parameterizing the integral kernel operators in Fourier Space.The pFNO adopts an architecture of maps-composition as lifting map P, a sequence of hidden maps H l for l = 1, 2, . . ., L, and a projection map Q.
The first map C : the parameters γ to the co-dimension of input function v(x), yielding v c (x).The second , lifts the input to a higher-dimensional functional space projects back to low dimension functional space, finally yielding v ′ (x).
Both P and Q are implemented using simple multilayer perceptrons(MLP).The hidden maps H l+1 are implemented as parametric Fourier layers: where W l ∈ R dε×dε and b l ∈ R dε are learnable weights and biases respectively, and σ is a ReLU activation function.Here, F and F −1 represent the Fourier Transform and its inverse respectively.The function R * l : C κ max ×dε × R dγ → C κ max ×dε acts on the truncated Fourier modes, transforming them as : where R l , R * l ∈ C κ max ×dε×dε are two learnable weight tensors, and D * l : R dγ → R κ max is a function converting the parameters γ into κ max -number of scaling ratios.This function consists of a two-stage map γ → D l (γ) → D * l (γ), with D l (γ) ∈ R N D outputting N D scaling ratios and implemented as a MLP.The second map hierarchically redistribute these ratio across the wave numbers.In one dimension(d = 1), the distribution map reads: One might observes that we can deactivate the second weight tensor R * l in Eq. ( 6) by enforcing the map D * l (γ) to output only zeros.This modification still enables learning of the parameter operator due to the concatenation map C. Such a modified method can viewed as a simple tweak to the baseline method of FNO 11 and will be referred as pFNO* in later section.

IV. NUMERICAL EXPERIMENTS AND RESULT DISCUSSIONS
In this section, we employ the pFNO and pCNN methods to learn flame evolution under hybrid instabilities arising from both Darrieus-Landau (DL) 19,20 and Diffusive-Thermal where Γ : ψ → −H(ψ x) is a linear singular non-local operator defined using the Hilbert transform H, or equivalently written as Γ : ψ → F −1 (|κ|F κ (ψ)) using the spatial Fourier transform F κ (ψ) and its inverse F −1 .
In Eq. (7), Ω is the density ratio between burned product and fresh reactant; Le * is a ratio (positive or negative) depending on the Lewis number of deficient reactant and another critical Lewis number.Introduce three constants (a, b, c) for variable transformation on time t = a 2 t, space x = b −2 x and the displacement function ψ(x, t) = c 2 ϕ(x, t), then equation (7)   can be rewritten with In this work we consider the flame front solution ϕ(x, t) of Eq. ( 8) in a channel domain subjected to periodic boundary condition, i.e. x ∈ D = (−π, π].One might notice that Eq. ( 8) admits a zero equilibrium solution being a flat flame (i.e.ϕ * (x, t) = 0), a perturbation analysis around this zero solution yields a linear dispersion relation with the perturbed solution being ϕ(x, t) = κ φκ (t)e iκx + φ * κ (t)e −iκx (superscript * denotes complex conjugate) and the Fourier mode of perturbation evolving as φκ (t) ≈ φκ (0) • e ω(κ)•t .Equations ( 8) and ( 9) present a straightforward approach to hybridize two flame instabilities of DT and DL mechanisms.This strategy is accomplished by specifying two parameters, ρ and β, while the remaining parameters (µ, ν, and τ ) can be determined by additional constraints outlined below.Initially, the parameter ρ (in between 0 and 1) is defined to allow for the continuous blending of these two instabilities.When ρ = 1, the Sivashinsky equation ( 8) yields a pure DL instability described by the Michelson-Sivashinsky (MS) equation 24 : whereas, at the other end (ρ = 0), it recovers the pure DT instability as described by the Kuramoto-Sivashinsky (KS) equation 25 : Secondly, the parameter β is determined as the largest value for which the dispersion relation of Eq. ( 9) equals zero (i.e.ω(β) = 0).This definition yields ν = µ−ρ.Consequently, we can prescribe β to establish the largest unstable wave number.To mitigate variability in the remaining parameters, a third constraint is imposed: the maximum value of ω(κ) over the interval 0 < κ < β must be 1/4.Furthermore, τ = ρβ/10 + (1 − ρ) is employed to better accommodate the timescales attributed to the various hybrid instabilities.This strategy allows for the determination of all remaining parameters given the values of ρ and β.This is illustrated in Figure 2, which presents dispersion relation plots and associated parameters.
Before proceeding further, it may be worthwhile to mention a few well-known results.
The KS equation ( 11) is often utilized as a benchmark example for PDE learning studies and is renowned for exhibiting chaotic solutions at large β.On the other hand, the MS equation (10), although less familiar outside the flame instability community, can be precisely solved using a pole-decomposition technique 27 , transforming it into a set of ODEs with finite freedoms.Moreover, at large β, the MS equation admits a stable solution in the form of a giant cusp front.However, at smaller β, the equation becomes susceptible to noise, resulting in unstable solutions characterized by persistent small wrinkles atop a giant cusp.Additional details about known theory can be found in references 23,[28][29][30][31][32][33][34][35] .

B. Training dataset
Equation ( 8) is tackled using a pseudo-spectral approach combined with a Runge-Kutta

C. Result analysis
The training datasets described in the previous section are utilized to train parametric solution advancement operators, denoted as Ĝ(γ) : ϕ(x; t) → ϕ(x; t + ∆ t ) with γ := (ρ, β) As a reminder, one ending value of ρ= 0 enables the pure DT instability while the other ending value of ρ = 1 active the pure DL instability.
In this study, three models-pFNO, pFNO*, and pCNN-described in Section II are employed to learn the two-parameter dependent operator Ĝ(ρ,β) .As explained in the last paragraph in Section II, pFNO* is a simple variant of the baseline FNO method 11 that includes the parameters in the codomain of the input function.On the other hand, pCNN has shown poor performance in learning the full operator Ĝ(ρ,β) , with a high training error exceeding 3 percent, see Table I.Therefore, we resort to two slightly restricted models (pCNN10 and pCNN40) which learn the single parameter (ρ) dependent operators Ĝ(ρ,β=10) and Ĝ(ρ,β=40) , with each model being trained using one third of total dataset at ρ=10 and 40 respectively.
The learned operator at given parameters is expected to make recurrent predictions of solutions over an extended period.The training for such operators aims not only for accurate short-term predictions but also for robust predictions of long-term solutions with statistics similar to the ground truth.As demonstrated in previous studies 17,18 , achieving this involves organizing the training data in a 1-to-20 pair, as expressed in Eq. ( 4), optimized for accurately predicting 20 successive steps of outputs from a single input over a range of parameter values.
] are compared in Fig. 8 for the normalized total front length ( (ϕ 2 x + 1) 1/2 dx/(2π)), in Fig. 9 for the model errors accumulated through recurrent predictions, and in Fig. 10 for the long-term auto-correlation function.This auto-correlation function characterizes the long-term recurrently predicted solutions: where ϕ * (x) denotes the predicted solutions obtained after a sufficiently long time.Numerical calculation for the expectation E in Eq. ( 12) is implemented by averaging over seven randomly initialized sequences of model predictions for a time duration 1000 < t/∆ t < 4000.
Moreover, for each of the learned models G θ,(β,ρ) ≈ Ĝ(β,ρ) , we compute an approximated dispersion relation: where J is the operator Jacobian This Jacobian is computed using the automatic differential tool (e.g., torch.autograd.functional.jacobian in PyTorch).Fig. 3 compares the dispersion relations by all models with the reference ones.Additionally, Fig. 3 shows one example of learned operator Jacobian, which is clearly diagonal dominant.

D. Findings
Overall learning: Our study underscores the robust learning capabilities of pFNO and pCNN methodologies in capturing the nuanced dynamics of flame front evolution, modulated by varying DL and DT instabilities blends.Both pFNO and pFNO* demonstrate good performance in learning the full two-parameter front evolution operator Ĝρ,β modulated by a ρ-varying blends of DL/DT instabilities as well as by a β-varying size of the largest unstable wavelength.While pCNN encounters difficulty in learning the full operator, the method still performs well in learning differently-hybrid instabilities when being restricted for the single-parameter operators Ĝρ,β=10 and Ĝρ,β=40 .
Short term learning: Across the board, all learned models (pFNO,pFNO* and pCNN10/40) demonstrate good accuracy in short-term predictions, with training/validation errors below 2 percent (Table I) and small accumulated errors (Fig. 9).This precision extends to various metrics, including front displacement, front slope, and normalized front length (Fig. t ≤ 50∆ t ), affirming the models' fidelity in capturing short-term dynamics.Moreover, pFNO demonstrates the smallest error and is the most accurate model for learning short term solutions.
Long term learning: Detailed analysis of reference solutions unveils distinct characteristics of isolated instabilities.At ρ = 1, DL fronts evolve toward a single giant cusp structure, either remaining stationary at small β = 10 (Fig. 4) or exhibiting noise-induced wrinkles at larger β = 40 (Fig. 5).Conversely, at ρ = 0, DT fronts adopt an overall flat shape interspersed with oscillatory wavy structures, with decreasing wavelength and amplitude as β increases from 10 to 40 (Figs. 4, 5).The slope plots for DT front evolution result in the typical zebra stripe pattern (Figs. 6, 7).Intermediate values of ρ showcase a gradual transition between these features, with the front structure blending wavy oscillations together with the cusp shape.
Long-term predictions by all learned models (pFNO*, pFNO, pCNN10/pCNN40) accurately replicate these characteristic behaviors across diverse parametric configurations, encompassing pure DL and DT instabilities as well as blended scenarios.Quantitative comparisons through auto-correlation functions (Fig. 10) and total front length (Fig. 8) confirm the models' proficiency in capturing long-term solutions.
Learning challenges: However, a common challenge across both pFNO and pCNN models lies in over-predicting the impact of noise-induced wrinkles, particularly noticeable at small β = 10 (Fig. 4, 6).This tendency leads to an overestimation of the total front area, Extra finding: It is particularly interesting to point out that the two models pFNO and pFNO* learn well on the parametric-dependent linear dispersion relations, as seen in Fig.

3.
Except for a moderate level of mismatch at a few parameter conditions toward large ρ at 25 and 40, pFNO and PFNO* reproduce the relations quite accurately.
Such learning performance is impressive considering the fact that the data effective for learning these linear relations (i.e., the initial near-zero solutions) is just a tiny portion of the total dataset.For pCNN-based models, Fig. 3 shows pCNN10 learns the dispersion quite accurately while pCNN40 learned relations show a more significant deviation than ones by pFNO.

V. SUMMARY AND CONCLUSION
This paper delves into the potential of machine learning (ML) for understanding and predicting the behavior of flames experiencing hybrid instabilities.These instabilities arise Sivashinsky equation, we introduce two parameters: ρ and β.These parameters control the blending of DT and DL instabilities, as well as the cutoff wavelength for unstable flame behavior.
Our learning problem focuses on understanding the PDE solution time advancement operator under different parameter combinations.This operator, when repeatedly applied with its input solution as the output from the previous iteration, yields a time sequence of solutions of arbitrary length.We employ two recently developed operator-learning mod-   full two-parameter operator, enabling variation in both ρ and β.
Challenges: However, both pCNN and pFNO tend to overestimate noise-induced wrinkles associated with DL instability, leading to inaccurate predictions of the total flame area, especially at lower instability levels.
In conclusion, this work showcases the potential of operator learning methodologies for analyzing complex flame dynamics arising from hybrid instabilities.While challenges persist, particularly related to noise overestimation, these methods offer promising tools for understanding and predicting real-world flame behavior in combustion systems [36][37][38][39][40][41][42][43] .Future research directions may involve incorporating additional physical mechanisms or exploring

FIG. 1 .
FIG.1.The parametric CNN model adopted in this study is demonstrated for input function v(x i ) discretized at 1D mesh of 256 points, with L = 6 levels of encoding and decoding.Standard convolution layers are represented by gray rectangles, while Inception layers are depicted in magenta.

(
DT) 21,22 mechanisms.The dynamics of such unstable flame development are encapsulated by the Sivashinsky equation 26 .To facilitate parametric learning, we begin by reformulating the Sivashinsky equation, introducing two parameters that enable straightforward specification for blending the two instabilities and controlling the largest unstable wave numbers.By sampling across these parameters, we construct an extensive training dataset covering a range of relevant scenarios subjected to different DL/DT mixing.Subsequently, we present the results and compare the performance of the different methods in learning these hybrid instabilities. A. Governing equations Consider modeling the unstable development of a statistically planar flame front.Let t denote time and x represent the spatial coordinate along the normal direction to flame propagation.Introduce a displacement function ψ(x, t) : R × R → R describing the streamwise coordinate of a flame front undergoing intrinsic flame instabilities.Such evolution can be modeled by the Sivashinsky equation 26 :

FIG. 2 .
FIG. 2. Dispersion relations (left) and relevant parameter values at prescribed values of ρ and β.

( 4 , 5 )
figuration tuples (ρ, β), formed as the Cartesian product of three values for β in the range[10, 25, 40] and five values for ρ in the range [0, 1/4, 1/2, 3/4, 1].For each of the fifteen parametric configurations, we generate 250 sequences of short-duration solutions, as well as a single sequence of long-duration solutions.Each short solution sequence spans a time

FIG. 3 .
FIG. 3. Left three columns: Comparison of reference dispersion relations (Eq.9, solid lines) with those computed for the learned operators of all models (Eq.13, non-solid lines), where line colors indicate different ρ.Right column: Illustration of a learned operator Jacobian (Eq.14), with dark colors indicating small values.
especially pronounced at lower β values of 10 and 25 (Fig. 8 at ρ = 1).When learning for the hybrid DL and DT instabilities, excessive noisy wrinkles also show up in all the model predictions (at ρ = 3/4 in figs.6 and 7), however, the issue becomes less discernible toward smaller values of ρ when DT instability plays a larger role, as also evident by the front length in Fig. 8.

FIG. 7 .
FIG. 7. Comparison of front slopes predicted by pFNO* and pCNN40 at β = 40.Other details are the same as in Fig. 6.

TABLE I .
Relative L2 train/validation errors for all operator learning networks Model Parameter configurations (β, ρ) Train L2 Valid.L2 t = 0.15.Additionally, each sequence starts from random initial conditions ϕ 0 (x) sampled from a uniform distribution over the range [0, 0.03].The long sequence covers a time duration of 0 ≤ t ≤ 18, 750 and comprises 125,000 consecutive solutions outputted at the same interval ∆ t .A validation dataset is similarly created for all fifteen parameter tuples, but it contains only 10 percent of the data present in the training dataset.

Table I
17,18nts the relative training/validation errors for various models.The validation errors in TableIare consistent with those reported in our previous work17,18.Additional details on training and model hyper-parameters are provided in Appendix A.