A&A 420, 411-422 (2004)
DOI: 10.1051/0004-6361:20035721

The cellular burning regime in type Ia supernova explosions

I. Flame propagation into quiescent fuel

F. K. Röpke1 - W. Hillebrandt1 - J. C. Niemeyer2


1 - Max-Planck-Institut für Astrophysik, Karl-Schwarzschild-Str. 1, 85741 Garching, Germany
2 - Universität Würzburg, Am Hubland, 97074 Würzburg, Germany

Received 21 November 2003 / Accepted 4 March 2004

Abstract
We present a numerical investigation of the cellular burning regime in Type Ia supernova explosions. This regime holds at small scales (i.e. below the Gibson scale), which are unresolved in large-scale Type Ia supernova simulations. The fundamental effects that dominate the flame evolution here are the Landau-Darrieus instability and its nonlinear stabilization, leading to a stabilization of the flame in a cellular shape. The flame propagation into quiescent fuel is investigated addressing the dependence of the simulation results on the specific parameters of the numerical setup. Furthermore, we investigate the flame stability at a range of fuel densities. This is directly connected to the questions of active turbulent combustion (a mechanism of flame destabilization and subsequent self-turbulization) and a deflagration-to-detonation transition of the flame. In our simulations we find no substantial destabilization of the flame when propagating into quiescent fuels of densities down to $\sim $ $10^7 ~{\rm g} ~{\rm cm}^{-3}$, corroborating fundamental assumptions of large-scale SN Ia explosion models. For these models, however, we suggest an increased lower cutoff for the flame propagation velocity to take the cellular burning regime into account.

Key words: stars: supernovae: general - hydrodynamics - instabilities

1 Introduction

Type Ia supernovae (SNe Ia) are usually attributed to thermonuclear explosions of white dwarf (WD) stars consisting of carbon and oxygen. Throughout this paper we will refer to the currently favored model in which the combustion proceeds in the form of a flame starting near the center of the WD and traveling outward. More specifically, the initial mode of flame propagation is assumed to be the so-called deflagration mechanism, in which the combustion wave is mediated by microphysical transport processes giving rise to a subsonic flame speed. Underlying our assumptions is the canonical single-degenerate Chandrasekhar-mass SN Ia model. For a review of SN Ia explosion models we refer to Hillebrandt & Niemeyer (2000).

Despite the success of recent attempts to numerically model Type Ia supernova explosions (Gamezo et al. 2003; Reinecke et al. 2002b), some basic questions regarding the explosion mechanism remain unanswered. This is mainly due to the fact that the whole range of scales relevant to this problem cannot be resolved in a single simulation in the foreseeable future. From the radius of the exploding white dwarf star down to the flame thickness it covers more than 11 orders of magnitude. Therefore, large-scale simulations provide an insight into the flame dynamics on the largest scales but have to rely on assumptions about the physics of flame propagation on smaller scales, which is poorly understood so far. However, the motivation for this approach is that the key feature of the SN Ia explosion is turbulent combustion. Turbulence is driven from large-scale instabilities - mainly the buoyant Rayleigh-Taylor instability and the Kelvin-Helmholtz (shear) instability - that give rise to the formation of a turbulent eddy cascade. The turbulent eddies wrinkle the flame and increase its surface. This enhances the net burning rate and accelerates the flame. In this way the energy production rate can reach values that are sufficient to power a SN Ia explosion. The underlying assumption of this model is that the flame evolution is dominated by the turbulent cascade originating from instabilities on large scales. Although there exist some theoretical ideas corroborating the assumption of flame stability on small scales, a justification by means of a hydrodynamical simulation has not provided definitive answers yet (Niemeyer & Hillebrandt 1995; Niemeyer & Woosley 1997). The results that will be reported in the present paper are, however, a contribution in this direction.

The work presented in the following stands in direct succession of Röpke et al. (2003), where the basic method was introduced and applied to describe flame propagation into quiescent fuel for a density of the unburnt material of $\rho_ {\rm u}= 5 \times 10^7 ~{\rm g} ~{\rm cm}^{-3}$. This example demonstrated that the applied method appropriately models the flame propagation on specific scales in the SN Ia explosion. Moreover, it confirmed for the first time by means of a full hydrodynamical simulation that the stabilization of the flame front in a cellular pattern known from chemical flames holds for thermonuclear flames under the conditions of SN Ia explosions. In the present paper we will extend the study of flame propagation to a wider range of fuel densities. The aim is to test the stability of the cellular pattern in dependence of this parameter. This will provide an overview over possible effects of flame propagation into quiescent fuel. We will additionally investigate the flow field resulting from the flame configuration in greater detail and measure the effective flame propagation velocities for the resulting flame structures.

The theoretical context of this study will be introduced in Sect. 2. After a brief description of the applied numerical methods in Sect. 3 we will present the results of numerical simulations of the flame evolution in Sects. 4 and 5. Finally, conclusions regarding the significance of the numerical investigations for large-scale SN Ia models will be drawn.

   
2 Theoretical background

Our study of flame evolution focuses on scales where flame propagation decouples from the turbulent cascade. On larger scales, this cascade dominates flame propagation by wrinkling the flame and enhancing its surface. But in order to contribute to the flame wrinkling, a turbulent eddy has to deform the flame significantly in a time shorter than the flame crossing time over the spatial extent of that eddy. Thus, equating the eddy turnover time to this time scale, one obtains the cutoff scale  $l_{\rm Gibs}$

 \begin{displaymath} v'(l_{\rm Gibs}) = s_{\rm l}, \end{displaymath} (1)

where $s_{\rm l}$ denotes the so-called laminar burning velocity of a deflagration flame, i.e. the propagation speed of a planar deflagration wave which is determined by microphysical transport. The turbulent velocity fluctuations at length scale l are represented by v'(l). Equation (1) defines a length scale $l_{\rm Gibs}$, which was identified with the Gibson scale for chemical flames (Peters 1986). We will make use of this term in the context of thermonuclear flames in SN Ia explosions, too.

In our model, we describe the flame as a discontinuity between burnt and unburnt states, ignoring any internal structure. This simplification is exploited by the numerical method of flame tracking applied in our simulations, resulting in an increased dynamical range compared to models that fully resolve the flame structure. To justify our thin flame approximation, we provide a rough estimate of the Gibson scale. The turbulent velocity fluctuations v' in Eq. (1) are given by the scaling law of the turbulent eddy cascade. Kolmogorov scaling, $v'(l) \propto l^{1/3}$, provides a reasonable approximation. Hence one obtains from Eq. (1)

\begin{displaymath}l_{\rm Gibs}=L\left( \frac{s_{\rm l}}{v'(L)} \right)^3\cdot \end{displaymath} (2)

The large scale instabilities give rise to the formation of ascending bubbles. A typical size of a Rayleigh-Taylor bubble is $L \sim 10^7 ~{\rm cm}$, which defines the integral scale of the turbulence. The rise of these bubbles invokes shear velocities v'(L) of about $10^7 ~{\rm cm} ~{\rm s}^{-1}$. Values for the laminar burning velocity have been determined by means of one-dimensional simulations resolving the flame structure by Timmes & Woosley (1992). The authors also provide the following fitting formula to the data which we will use in our simulations:

 \begin{displaymath} s_{\rm l} = 92.0\times 10^5 \left( \frac{\rho_{\rm u}}{2 \ti... ...{X(^{12}{\rm C})}{0.5}\right]^{0.889} ~{\rm cm} ~{\rm s}^{-1}. \end{displaymath} (3)

However, we note that the fit does not reproduce the simulation results as accurately as the authors claim for the fuel density range that will be investigated in the following. According to this formula, the laminar burning velocity at a fuel density of $5 \times 10^7 ~{\rm g} ~{\rm cm}^{-3}$ is of the order of $10^6 ~{\rm cm} ~{\rm s}^{-1}$. This yields a Gibson scale of $\sim $ $10^4 ~{\rm cm}$. Thus, keeping in mind the indicative character of our study, we will arrange our simulations at length scales around this value, which is still well-separated from the flame width. It is also worth noting that flame evolution at this scale will not be affected by gravity effects, which are expected to have drastic impact on larger scales. The reason for this is that the comparison of the burning speed of the flame with velocity fluctuations induced by the Rayleigh-Taylor instability sets a minimum scale at which the flame is affected by the latter (Timmes & Woosley 1992). This scale, however, is expected to be above the Gibson length, since the Rayleigh-Taylor induced velocity fluctuations scale with the square root of the length scale (Davies & Taylor 1950) while the Kolmogorov law predicts a cubic root scaling of the velocity fluctuations in the turbulent cascade with the length scale.

Below the Gibson scale, flame propagation does not interact with the turbulent cascade, as in the context of SN Ia explosion was first pointed out by Niemeyer & Hillebrandt (1995). Here, it would propagate as a laminar flame. However, it is well-known that laminar flames are subject to a hydrodynamical instability (Darrieus 1938; Landau 1944). The origin of this so-called Landau-Darrieus instability (LD instability henceforth) is the refraction of the streamlines of the flow on the density contrast across the flame. Mass flux conservation leads to a broadening of the flow tubes in the vicinity of bulges of the perturbation. Here the local fluid velocity decreases and becomes lower than the fluid velocities at $\pm$$\infty$. These velocities correspond to the laminar burning speeds with respect to fuel and ashes in the frame of reference comoving with the flame. Therefore the burning velocity $s_{\rm l}$ of the flame is higher than the corresponding local fluid velocity and this leads to an increment of the bulge. The opposite holds for recesses of the perturbed front. In this way the perturbation keeps growing.

This can be observed in our simulations. Figure 1a depicts the absolute value of the fluid velocity normalized to the laminar burning velocity of the flame in the linear regime of flame evolution in one of our simulations. It is well in agreement with the theoretical expectations and confirmes that our numerical method is capable of reproducing the LD instability. This is not trivial as can be seen from the failure of the passive implementation to develop the correct flow field (cf. Röpke et al. 2003).

  \begin{figure} \par\subfigure[]{\includegraphics[width=6.5cm,clip]{0721fig1.eps} } \subfigure[]{\includegraphics[width=6.5cm,clip]{0721fig2.eps} }\end{figure} Figure 1: Fluid velocity in the simulation with $\rho_ {\rm u}= 5 \times 10^7 ~{\rm g} ~{\rm cm}^{-3}$ and $100 \times 100$ cells resolution. The x-axis and the y-axis show position in cm; the z-axis represents $\vert\ensuremath\mbox{\boldmath{$v$ }}\vert/s_{\rm l}$ a) snapshot at $t=0.01 ~{\rm s}$, b) snapshot at $t=0.09 ~{\rm s}$.
Open with DEXTER

By means of a linear stability analysis, Landau (1944) derived a dispersion relation between the growth rate of the amplitude $\omega$ and the wavenumber k of the perturbation:

 \begin{displaymath} \omega_{\rm LD} = k s_{\rm l} \frac{\mu}{1 + \mu}\left(\sqrt{1 + \mu - \frac{1}{\mu}}-1\right), \end{displaymath} (4)

with $\mu = \rho_{\rm u}/\rho_{\rm b}$.

This instability, however, is not necessarily a contradiction to the assumption of flame stability on small scales. In fact, one observes stable flame propagation in chemical combustion experiments. The reason for this is that the flame stabilizes in a cellular pattern as soon as the growth of the perturbation enters the nonlinear regime. This effect has been explained by Zel'dovich (1966) in a geometrical approach. The propagation of an initially sinusoidally perturbed flame is followed by means of Huygens' principle (see also Fig. 1 in Röpke et al. 2003). At the intersections of two bulges of the front a cusp forms after a while. Here Huygens' principle is no longer applicable and the flame evolution enters the nonlinear regime. The propagation velocity at the cusp exceeds $s_{\rm l}$:

\begin{displaymath}v_{\rm cusp} = \frac{s_{\rm l}}{\cos \theta}, \end{displaymath} (5)

where $\theta$ denotes the inclination angle of the front adjacent to the cusp. This effect balances the perturbation growth due to the LD instability and leads to a stabilization of the flame in a cellular configuration. Figure 1b shows the absolute value of the velocity after a cusp has formed.

Assuming a parabolic shape of the cells, the temporal evolution of the amplitude of the perturbation $\hat{x}$ can now be written as (Zel'dovich 1966)

 \begin{displaymath} \frac{{\rm d} \hat{x}}{{\rm d} t} = \omega_{\rm LD} \hat{x} - \frac{2}{\pi^2} k^2 s_{\rm l} \hat{x}^2. \end{displaymath} (6)

The quadratic damping term accounting for the cellular stabilization introduces an extension to the linear stability analysis resulting in the dispersion relation (4). If one further assumes that the flame front advances with the constant laminar burning speed $s_{\rm l}$ independently of the flame shape, the effective speed  $v_{\rm eff}$ of the mean position of a flame wrinkled due to instabilities is given in terms of its surface A' compared to the surface of a planar flame A(0):

 \begin{displaymath} v_{\rm eff} = s_{\rm l} \frac{A'}{A_0}\cdot \end{displaymath} (7)

The flame thus accelerates with increasing surface area and the task to determine the effective flame propagation velocity reduces to the determination of the increase in flame surface area. From the assumption of a parabolic cell shape Zel'dovich et al. (1980) derived the acceleration of the flame corresponding to the increase in its surface:

 \begin{displaymath} \frac{v_{\rm cell}}{s_{\rm l}} = 1+ \epsilon, \qquad \epsilon = \frac{\pi^2}{24} \left(1-\frac{1}{\mu}\right)^2\cdot \end{displaymath} (8)

The cellular stabilization is the reason that large-scale SN Ia simulations assume flame stability on small scales. This is, however, not yet well confirmed. The cellular pattern may break up under certain conditions, e.g. density jumps over the flame front or interaction with turbulence. Blinnikov & Sasorov (1996) suspected that a fractalization of the cellular flame may drastically accelerate the flame at low fuel densities and Niemeyer & Hillebrandt (1995) reported on indications that the flame may destabilize at fuel densities of $\rho_{\rm u} \approx 5 \times 10^7 ~{\rm g} ~{\rm cm}^{-3}$. However, a failure of the cellular flame stabilization at $\rho_ {\rm u}= 5 \times 10^7 ~{\rm g} ~{\rm cm}^{-3}$could not be confirmed by Röpke et al. (2003). Nevertheless, the possible loss of flame stability at certain densities of unburnt material could have drastic impact on large-scale SN Ia models and has to be soundly tested. Besides the break-down of one fundamental assumption of these models, it could also have effects that seem appealing. Spectra produced by empirical one-dimensional SN Ia models fit the observations remarkably well if a transition of the flame propagation mode from subsonic deflagration to supersonic detonation is applied at low densities (Höflich & Khokhlov 1996; Iwamoto et al. 1999). The drawback of these delayed detonation models is that the mechanism providing the transition is unclear (Niemeyer 1999). As one possibility, an active production of turbulence and self-turbulization of the flame after a break-down of its cellular stabilization was proposed under the name active turbulent combustion by Niemeyer & Woosley (1997) based on an idea by Kerstein (1996). Following the reasoning provided by Niemeyer & Woosley (1997), this effect is directly connected to the burning regime we are aiming at in our simulations.

The above considerations point out the need for a thorough investigation of the cellular burning regime in SN Ia explosions at scales around or below the Gibson length. Numerical simulations addressing this question will be presented in Sects. 4 and 5.

   
3 Numerical methods and simulation setup

Since the present investigation focuses on the flame dynamics on small scales, it is performed in a simplified setup as compared to a realistic situation in a SN Ia explosion. The white dwarf structure is ignored in the sense that we assume pure carbon fuel and neglect gravity (cf. Sect. 2) as well as any density gradient. Though the flame propagation speed depends on the fuel density, a density gradient is unlikely to have a significant effect in our case. Its value is very small in late SN Ia explosion stages so that density changes over the scale of flame wrinkling are negligible.

The numerical methods we apply are based on the work by Reinecke et al. (1999) and have been described by Röpke et al. (2003). Therefore we will be brief in this and only mention some keywords here referring to Röpke et al. (2003) for the details.

Since our simulations aim at the range around the Gibson scale, which is well above the flame thickness for the fuel densities we are going to apply, it is justified to model the flame as a discontinuity between unburnt and burnt states. This flame description ignores any internal structure and therefore does not intrinsically provide the value of the laminar burning velocity. This also implies that, in a first approach, hydrodynamics and flame propagation can be modeled in separate steps.

The hydrodynamics is modeled applying the piecewise parabolic method (PPM) - a higher order Godunov scheme developed by Colella & Woodward (1984) - in the PROMETHEUS implementation (Fryxell et al. 1989). The corresponding part of the code treats only the fluid dynamics and does not account for the flame propagation. Since the flame is modeled in the discontinuity approximation, its propagation can be described applying the level-set method (Osher & Sethian 1988). Röpke et al. (2003) could show that the key feature that enables us to reproduce the theoretically anticipated flame evolution resulting from the LD instability is a specific implementation of flame/flow coupling, namely the in-cell reconstruction/flux-splitting technique proposed by Smiljanovski et al. (1997). It couples the flame propagation accurately to the flow and therefore allows the investigation of effects that originate from hydrodynamics, such as the LD instability. Our implementation corresponds to what was termed "complete implementation of the level-set method'' by Reinecke et al. (1999) with some minor changes (cf. Röpke et al. 2003). All simulations presented in the following were performed in two spatial dimensions. As our intention is to study flame dynamics rather than to provide a detailed nucleosynthetic description of burning in SNe Ia, we model the thermonuclear reactions in an extremely simplified way, converting the pure 12C fuel in a one-step reaction instantly to 56Ni and releasing $9.28667 \times 10^{17} ~ {\rm erg} ~ {\rm g}^{-1}$.

The basic approach is to study the flame evolution under the influence of a small perturbation (which can be introduced by noise or - as in our case - be initially imprinted on the flame structure). This perturbation should grow due to the LD instability. In the nonlinear regime, the formation of a cellular pattern should inhibit the further growth of the perturbation and stabilize the flame. In a previous publication (Röpke et al. 2003) we demonstrated that this mechanism holds for thermonuclear flames in degenerate matter. From theoretical considerations and semi-analytical models (resulting for instance from the Sivashinsky equation, e.g. Gutman & Sivashinsky 1990), it can be expected that the general features of the evolution of the flame front shape depend significantly on the numerical setup used in the simulations. Possible parameters that influence the flame evolution are the overall geometry of flame propagation, the width of the computational domain compared to the length scale of perturbations, the resolution, the boundary conditions, sources of numerical noise (which is likely to be different when the simulation is parallelized to a varying number of sub-processes), and the initial flame shape. Therefore care has to be taken in choosing the specific setup depending on the questions that are addressed by the simulations, and also in the conclusions drawn from the simulation results.

Regarding the overall flame geometry, two cases are commonly studied in the literature: an on average planar flame geometry and a flame that is on average circularly expanding. Although a naive approach would choose the second scenario for the supernova explosion, it is probably not the appropriate description of the flame propagation there. The flame evolution on scales of the star, where expansion effects are most pronounced, quickly becomes dominated by the Rayleigh-Taylor instability. This leads to a flame evolution completely diverging from a circular (or in three dimensions spherical) geometry (cf. Reinecke et al. 2002a). It rather proceeds in raising bubbles of burnt material. Nevertheless, it could be argued that these structures again partly resemble a spherically expanding geometry. Although the case of a circularly expanding flame reveals very interesting physical effects (as repeated mode splitting of the cells resulting from expansion effects and a possibly resulting fractalization of the flame front, cf. Blinnikov & Sasorov 1996), we will not follow this approach here for two reasons: First, in the scope of our implementation it is prohibitively expensive to follow a circular flame evolution for a sufficiently long time and, second, we aim at effects on scales around the Gibson length, at which global expansion effects are negligible. This has a technical advantage. It is much simpler to keep an overall planar flame in the center of the domain by choosing a comoving frame of reference. This enables us to study the long term flame evolution without requiring large computational domains. Thus it becomes possible to simulate the long-term flame evolution, as will be presented below. Note that the choice of an overall planar flame geometry determines a distinct flame evolution and possibly excludes the mechanism of repeated mode splitting.

Table 1: Setup values for the simulations of the flame evolution.

The influence of the width of the computational domain has already been discussed by Röpke et al. (2003). The simulations presented there lead to the conclusion that for sufficient numerical resolution the flame stabilizes in a single domain-filling cusp-like structure for periodic boundaries transverse to the direction of flame propagation. However, with increasing domain width (corresponding to a higher resolution in the simulations) this fundamental structure becomes superimposed by a short-wavelength cellular pattern. This is in accord with results from semi-analytical studies (Gutman & Sivashinsky 1990).

The discussion of the influence of the boundary conditions and the initial shape of flame perturbations will be postponed to the next section.

The "experimental setup'' applied our simulations is the following: The spatial extent of the computational domain was set to correspond to scales around the Gibson length and was fixed according to the discussion above. The flame was initialized in the center of the domain in an on average planar vertical shape with the unburnt material on the right hand side and the burnt material on the left hand side, so that in the laboratory frame of reference the flame would propagate to the right. The domain was set to be periodic in the y-direction. On the left boundary of the domain an outflow condition was enforced and on the right boundary we imposed an inflow condition with the unburnt material entering with the laminar burning velocity $s_{\rm l}$. This would yield a computational grid comoving with a planar flame. However, the LD instability leads to a growth of the perturbation and therefore increases the flame surface. According to Eq. (7), this causes an acceleration of the flame and therefore it is necessary to take additional measures to keep the mean position of the flame centered in the domain. One possibility is to detect the mean location of the flame and to simply shift the grid to keep the flame in the center. This method is consistent with the boundary conditions applied in x-direction. The described method allows the study of the long term flame evolution, so that a detailed investigation of the nonlinear stage of the LD instability becomes possible. To induce the development of perturbations, we usually imposed a sinusoidal perturbation on the initial flame shape.

The state variables were set up with values for the burnt and unburnt states obtained from (pseudo-)one-dimensional simulations performed with the "passive implementation'' of the level-set method (Reinecke et al. 1999), imposing a value for the density of the fuel. A compilation of the relevant setup values for different fuel densities is given in Table 1. Note that frequently in this paper fuel densities refer to a label of a specific set of values rather than giving the accurate fuel density. The laminar burning speed is calculated according to Eq. (3).

  \begin{figure} \par\includegraphics[width=14.4cm,clip]{0721fig3.eps} \end{figure} Figure 2: Flame propagation into quiescent fuel at $\rho_ {\rm u}= 5 \times 10^7 ~{\rm g} ~{\rm cm}^{-3}$; resolution: $200 \times 200$ cells, periodic boundary conditions; snapshots taken at a) 0, b) 5, c) 7.5, d) 10, e) 15, and f) 30 growth times $\tau _{\rm LD} = \omega _{\rm LD}^{-1}$ of a perturbation with $\lambda =200 \Delta x$. The vorticity $\zeta $ (cf. Eq. (9)) is color-coded and the flame position is indicated by solid white curves.
Open with DEXTER

   
4 General features of flame evolution

The flame evolution for the exemplary case of a fuel density of $5 \times 10^7 ~{\rm g} ~{\rm cm}^{-3}$ was described by Röpke et al. (2003). The results obtained in this reference were

In this section we extend this investigation by presenting two simulations of the long-term flame evolution, focusing on the nonlinear regime. The two simulations were carried out on a regular Cartesian grid with $200 \times 200$ cells with a cell width of $\Delta x = 160 ~{\rm cm}$. The fuel density was $5 \times 10^7 ~{\rm g} ~{\rm cm}^{-3}$. In both cases we perturbed the flame front initially in a sinusoidal way with a wavelength of 1/6 of the domain width, but we applied different boundaries transverse to the flame propagation.

Figure 2 illustrates the flame evolution for periodic boundary conditions in the y-direction. The first snapshot shows the initial setup with the sinusoidal perturbation imposed on the flame. In the second snapshot the flame has already reached the nonlinear regime and exhibits a cellular shape. Here, the onset of an irregularity of the pattern is apparent. This effect leads to a merging of cells observable in the next three snapshots. Finally, the flame reaches a steady state which is a single domain-filling cusp. Note that the snapshots in Fig. 2 were not taken at equal time intervals. The final flame configuration was followed for more than 10 growth times $\tau _{\rm LD} = \omega _{\rm LD}^{-1}$corresponding to a perturbation of the width of the domain in order to ensure its stability (cf. Fig. 4a).

In the snapshots of the flame evolution depicted in Fig. 2 the vorticity of the flow field, $\ensuremath\mbox{\boldmath {$\zeta$ }} = \ensuremath\mbox{\boldmath {$\nabla$ }} \times \ensuremath\mbox{\boldmath {$v$ }}$, is color-coded. For our two-dimensional simulations the (scalar) vorticity reads

 \begin{displaymath} \zeta = \frac{\partial v_y}{\partial x} - \frac{\partial v_x}{\partial y}\cdot \end{displaymath} (9)

Note that we color-code the logarithm of the absolute value of the vorticity and therefore even slight deviations from zero vorticity are visible. The flame produces vorticity in the burnt material. This effect is especially strong behind cusps. Analytical approximations of the flame evolution frequently make the assumption of a potential flow ahead of the flame. In contrast to that, we observe the development of vorticity in the fuel region. This, however, is not surprising, since the flame propagates subsonically and the regions downstream and upstream of the flame are causally connected.

The merging proceeds such that some of the cells grow and the smaller cells disappear in the cusp between larger cells, as can be clearly seen in Fig. 4a. This effect forms the basis of flame stability as explained by Zel'dovich et al. (1980). There, the authors addressed the phenomenon that flames in experiments are often observed to stabilize in large cells. This, however, seems to be in contradiction with the result from Landau's linear stability analysis (Eq. (4)), because it predicts higher growth rates for perturbations with smaller wavelengths, which thus should dominate the flame structure. Even taking into account a modified dispersion relation resulting from a finite flame thickness (Markstein 1951) cannot explain this phenomenon. Zel'dovich et al. (1980) argue that the stable long-wavelength flame pattern is a result of the flow that develops upstream of the flame. A velocity component tangential to the flame is directed toward the cusp and advects small perturbations into this direction. A gradient in this tangential velocity stretches the perturbation wavelength and thereby retards its growth. Finally, the small perturbation disappears in the cusp (originally Zel'dovich et al. (1980) analyzed the flame propagation in tubes, but they point out the similarity of this configuration with a cellular flame). Zel'dovich et al. (1980) give an analytical description of the effect applying a WKB-like approximation.

The described mechanism acts also in our numerical model, as was discussed by Röpke et al. (2003) in connection with the flame structure superposed by a short-wavelength cellular pattern for high-resolved simulations. A close-up of the flow field in the vicinity of a cusp taken from one of our simulations (see Fig. 5) demonstrates the ability of our implementation to reproduce the effect. The arrows represent the velocity field after the formation of a cusp. Note that the streamlines converge toward the cusp ahead of the front. This leads to the formation of a layer upstream of the front in which the fluid velocity is directed toward the cusp. Thus the required flow field providing stability of long-wavelength patterns establishes itself in our simulations. In the downstream region, the flow diverges from the cusp.

The flame evolution in for reflecting boundary conditions in the y-direction is shown in Fig. 3. The mechanism of cellular stabilization acts similar to the case of periodic boundaries. However, the alignment of the steady-state of the flame shape is different. Whereas the periodic case developed a cusp in the center of the domain, in the reflecting case the crest of the pattern centers in the domain and the recesses of the front form a structure similar to a "half-cusp'' at the boundaries. This case is analogous to the configuration studied by Zel'dovich et al. (1980).

  \begin{figure} \par\includegraphics[width=14.4cm,clip]{0721fig4.eps} \end{figure} Figure 3: Flame propagation into quiescent fuel at $\rho_ {\rm u}= 5 \times 10^7 ~{\rm g} ~{\rm cm}^{-3}$; resolution: $200 \times 200$ cells, reflecting boundary conditions; snapshots taken at a) 0, b) 5, c) 7.5, d) 10, e) 15, and f) 30 growth times $\tau _{\rm LD} = \omega _{\rm LD}^{-1}$ of a perturbation with $\lambda =200 \Delta x$. The vorticity $\zeta $ (cf. Eq. (9)) is color-coded and the flame position is indicated by solid white curves.
Open with DEXTER

Another aspect of these two simulations should be noted. Röpke et al. (2003) studied the evolution of a flame that was initially perturbed by only one single domain-filling mode. This was done to be more efficient in the investigation of the linear growth due to the LD instability, and it was argued that a single domain-filling structure develops independently of the exact shape of the initial perturbation. This is reconfirmed by the presented simulations in accord with expectations from numerical simulations of the Sivashinsky equation (e.g. Gutman & Sivashinsky 1990) and from the pole-decomposition solutions of that equation (Thual et al. 1985). In both the simulation with periodic boundaries (cf. Fig. 2) and the simulation with reflecting boundary conditions (see Fig. 3) the cells merge. The steady state of the flame resulting from this process is a single domain-filling cell. Even a non-regular initial perturbation is not expected to alter this result. The regular sinusoidal initial perturbation was chosen in order to have a well-defined perturbation wavelength which can then be compared to theoretically predicted growth rates.

   
5 Flame stability at different fuel densities

Further numerical experiments addressed the flame stability at different fuel densities. With lower fuel density, the Atwood number - defined as ${\it At} = \vert\rho_{\rm u} - \rho_{\rm b}\vert/(\rho_{\rm u}+\rho_{\rm b})$ - and, equivalently, the density contrast $\mu$ over the flame front increase (cf. Table 1). At the same time the laminar burning velocity of the flame decreases. In the following, tests of the flame stability at a variety of fuel densities will be presented. These use the setup introduced in Sect. 3 with periodic boundaries in the direction transverse to the flame propagation. We again refer to Table 1 for a compilation of the setup values. The resolution chosen for this study was $200 \times 200$ cells. From the simulations presented in the previous section, it follows that an initial flame perturbation with a wavelength corresponding to the width of the computational domain is sufficient to see the flame stabilize in a universal steady state. In the case under consideration one expects a single centered domain-filling cusp-like structure. Although the restriction to an initial domain-filling perturbation wavelength may appear rather artificial, it avoids the time-consuming cell merging process before the final steady-state flame structure is reached.

  \begin{figure} \par\subfigure[]{\includegraphics[width=8.2cm,clip]{0721fig5.eps} } \subfigure[]{\includegraphics[width=8.2cm,clip]{0721fig6.eps} }\end{figure} Figure 4: Evolution of the flame front for $\rho_ {\rm u}= 5 \times 10^7 ~{\rm g} ~{\rm cm}^{-3}$ a) with periodic boundaries in the y-direction, b) with reflecting boundaries in the y-direction. Each contour marks a time step of 0.5 $\tau_{\rm LD}(\lambda_{\rm pert}=3.2 \times 10^4~{\rm cm})$. The contours are artificially shifted for better visibility and do not reflect the actual flame propagation velocity.
Open with DEXTER


  \begin{figure} \par\includegraphics[width=8.2cm,clip]{0721fig7.eps} \end{figure} Figure 5: Velocity field in a simulation with a fuel density of $\rho_ {\rm u}= 5 \times 10^7 ~{\rm g} ~{\rm cm}^{-3}$ and $100 \times 100$ cells resolution at $t=0.09 ~{\rm s}$; solid curve indicates flame position.
Open with DEXTER

Figure 6 shows the temporal evolution of the flame front for fuel densities of $2.5\times 10^7 ~{\rm g}~{\rm cm}^{-3}$ and $5 \times 10^7 ~{\rm g} ~{\rm cm}^{-3}$. In both cases the flame stabilizes in a single domain-filling cusp-like structure. Note that the evolution time has been normalized to the growth time of the LD instability $\tau _{\rm LD} = \omega _{\rm LD}^{-1}$corresponding to the initial perturbation wavenumber. This ensures comparability between the different simulations. As in Fig. 4, the spacing between individual contours is, however, artificial since the original simulations were performed in the frame of reference comoving with the flame. The spacing has been chosen in a way that the contours fill the plot window, and does not reflect flame propagation.

The outcome from simulations with even lower fuel densities reveals a completely different flame evolution. Figure 7 shows the flame for fuel densities of $1\times 10^7 ~{\rm g}~{\rm cm}^{-3}$ and $1.25\times 10^7 ~{\rm g}~{\rm cm}^{-3}$. Although the initial flame evolution resembles that of $\rho_{\rm u} \approx 5 \times 10^7 ~{\rm g} ~{\rm cm}^{-3}$, the forming cusp lacks long-term stability. As can be seen from the plots, the stable structure breaks up from the cusps outward and the flame subsequently evolves in an irregular pattern. The origin of these effects is most likely numerical. For a fuel density of $1\times 10^7 ~{\rm g}~{\rm cm}^{-3}$ a simulation run with $600 \times 600$ cells resolution did not result in a disruption of the cusp. Thus, a high grid resolution is required in order to describe a stable cusp properly at low fuel densities ($\lesssim$ $1.25\times 10^7 ~{\rm g}~{\rm cm}^{-3}$).

A rather unexpected result is obtained from simulations for fuel densities around $1\times 10^8 ~{\rm g}~{\rm cm}^{-3}$. Here, the flame does not propagate in a stable manner in the first stages of the flame evolution. This is illustrated by Fig. 8. Small-wavelength perturbations are superposed on the initial flame shape shortly after the beginning of the simulation. These appear to dominate the flame evolution for a while (see in particular Fig. 8b) but then the flame stabilizes in the single-cell structure. This implies that (i) the flame structure is less stable against small-wavelength perturbations at these fuel densities and (ii) the mechanism for a stabilization in a preferred long-wavelength pattern as analytically predicted (see the discussion in Sect. 4) does finally stabilize the flame evolution.

  \begin{figure} \par\subfigure[]{\includegraphics[width=8.3cm,clip]{0721fig8.eps} } \subfigure[]{\includegraphics[width=8.3cm,clip]{0721fig9.eps} } \end{figure} Figure 6: Evolution of the flame front at a) $\rho_{\rm u}=2.5\times 10^7 ~{\rm g}~{\rm cm}^{-3}$ and b) $\rho_ {\rm u}= 5 \times 10^7 ~{\rm g} ~{\rm cm}^{-3}$. The contours mark evolution steps of 0.2 $\tau _{\rm LD}$. Every 10th contour is highlighted bold for better visibility.
Open with DEXTER


  \begin{figure} \par\subfigure[]{\includegraphics[width=8.2cm,clip]{0721fig10.eps} } \subfigure[]{\includegraphics[width=8.2cm,clip]{0721fig11.eps} } \end{figure} Figure 7: Evolution of the flame front at a) $\rho_{\rm u}=1\times 10^7 ~{\rm g}~{\rm cm}^{-3}$ and b) $\rho_{\rm u}=1.25 \times 10^7 ~{\rm g}~{\rm cm}^{-3}$. The contours mark evolution steps of 0.2 $\tau _{\rm LD}$. Every 10th contour is highlighted bold for better visibility.
Open with DEXTER


  \begin{figure} \par\subfigure[]{\includegraphics[width=8.2cm,clip]{0721fig12.eps} } \subfigure[]{\includegraphics[width=8.2cm,clip]{0721fig13.eps} } \end{figure} Figure 8: Evolution of the flame front at a) $\rho_{\rm u}=7.5\times 10^7 ~{\rm g}~{\rm cm}^{-3}$ and b) $\rho_{\rm u}=1\times 10^8 ~{\rm g}~{\rm cm}^{-3}$. The contours mark evolution steps of 0.2 $\tau _{\rm LD}$. Every 10th contour is highlighted bold for better visibility.
Open with DEXTER

The initial flame destabilization becomes less pronounced in flame evolution at even higher densities. Figure 9 shows an example with $\rho_{\rm u} = 1\times 10^9 ~{\rm g}~{\rm cm}^{-3}$.

  \begin{figure} \par\includegraphics[width=8.5cm,clip]{0721fig14.eps} \end{figure} Figure 9: Evolution of the flame front at $\rho_{\rm u} = 1\times 10^9 ~{\rm g}~{\rm cm}^{-3}$. The contours mark evolution steps of 0.2 $\tau _{\rm LD}$. Every 10th contour is highlighted bold for better visibility.
Open with DEXTER

We measured the flame acceleration in these simulations via the increase in flame surface area (cf. Eq. (7)). This method has been discussed by Röpke et al. (2003). Although there are some minor differences compared to measuring the effective flame propagation velocity directly (most importantly, Eq. (7) strictly holds only for a flame idealized as a discontinuity, which is fulfilled to a reasonable degree in our simulations, but probably not quite exactly, cf. Röpke et al. 2003), this approach can be expected to yield rather reliable results. Some examples of the temporal evolution of the flame surface area are plotted in Fig. 10. This figure illustrates the different cases of flame evolution. While the simulations of low fuel densities (Fig. 10a) fail to converge, the flame reaches a steady-state propagation velocity for higher fuel densities. The initial destabilization of the flame pattern and final stabilization for fuel densities around $10^8 ~{\rm g} ~{\rm cm}^{-3}$ is apparent in the evolution of the corresponding flame surface area (Fig. 10c).

Table 2: Flame propagation velocities in the cellular regime.

Table 2 compares the final steady-state velocities of the cellular flames to the theoretical predictions from Eqs. (8). These equations, of course, ignore any structure superimposed on the flame front and additionally assume a parabolic shape of the smooth cells. Both effects obscure the expected trend of decreasing $v_{\rm cell}/s_{\rm l}$ with increasing fuel density. This can be clearly seen from Figs. 6-9, where the deviations are most pronounced for fuel densities around $10^8 ~{\rm g} ~{\rm cm}^{-3}$. Therefore the measurements agree only in order of magnitude with the theoretical values. This indicates that ignoring the superimposed small-wavelength pattern is a too simple approach and underestimates the actual flame surface.

6 Conclusions

The simulations presented in Sects. 4 and 5 provide a deeper insight into the dynamics of thermonuclear flame fronts subject to the conditions of SN Ia explosions. One fundamental assumption was that the flame propagates into quiescent fuel.

  \begin{figure} \par\includegraphics[width=14.2cm,clip]{0721fig15.eps} \end{figure} Figure 10: Increase in flame surface area at fuel densities of a) $\rho_{\rm u}=1.25 \times 10^7 ~{\rm g}~{\rm cm}^{-3}$, b) $\rho_ {\rm u}= 5 \times 10^7 ~{\rm g} ~{\rm cm}^{-3}$, c) $\rho_{\rm u}=1\times 10^8 ~{\rm g}~{\rm cm}^{-3}$, d) $\rho_{\rm u} = 1\times 10^9 ~{\rm g}~{\rm cm}^{-3}$.
Open with DEXTER

In Sect. 4 we addressed the question of the influence of the boundary conditions and the shape of the initial flame perturbation on the long-term flame evolution. The conclusions that can be drawn from our simulations are:

1.
In the given setup, the initial perturbation shape is not retained. Although at first it develops into a cellular structure of the same wavelength in the nonlinear regime, the flame later shows the tendency to align in the shape of a single domain-filling cell. The transition to this steady-state pattern takes place due to the growth of preferred cells from the initial perturbation and the disappearance of smaller cells in newly formed cusps. This effect gives rise to a "merging process'' of cells, consistent with numerical solutions (Gutman & Sivashinsky 1990) and with pole decomposition analysis (Thual et al. 1985) of the Sivashinsky equation (Sivashinsky 1977).
2.
The merging of the small cells proceeds smoothly without drastic effects on the flame shape. In particular it does not lead to a loss of flame stability.
3.
The cases of reflecting and periodic boundary conditions differ in the alignment of the ultimately developed cell. In the case of periodic boundary conditions the cusp develops at the center of the domain and the crest tends toward the boundary. In principle no specific alignment should be preferred for periodic boundaries. This is not true for reflecting boundary conditions. Here, the crest centers in the domain, similar to the simulation with adiabatic boundaries presented by Gutman & Sivashinsky (1990).
These results support the generality of the simulations presented by Röpke et al. (2003) and were used to set up the study of the influence of the fuel density on the flame stability that was presented in Sect. 5. Here, the general features of the flame evolution resemble those of Röpke et al. (2003) and Sect. 4 of the present paper. However, some peculiarities were observed. At very low fuel densities ($\lesssim$ $1.25\times 10^7 ~{\rm g}~{\rm cm}^{-3}$), the cusp-like flame shape breaks up. It was argued that this is a numerical effect. Nevertheless, if the reason for the disruption of the stable flame pattern is numerical noise, then one could conclude that flames at lower fuel densities are likely to become more sensitive to velocity fluctuations in the fuel region. These, however, can be expected to be present in realistic scenarios of SN Ia explosions. Therefore we will address that question in a subsequent study.

A destabilization of the flame structure shortly after initialization was observed at densities around $10^8 ~{\rm g} ~{\rm cm}^{-3}$. However, after $\sim $ $2 \tau_{\rm LD}$ the flame stabilized again forming the cusp-like steady-state structure. Thus, we conclude that the flame tends to stabilize in a cellular pattern in the parameter space explored in our simulations.

The implications of the presented investigation for large-scale SN Ia models can be summarized as follows:

(i)
The flame is expected to propagate in a stable way at least down to fuel densities of around $10^7 ~{\rm g} ~{\rm cm}^{-3}$. This corroborates fundamental assumptions of large-scale SN Ia models (e.g. Reinecke et al. 2002b) which rely on a smooth flame evolution on unresolved scales. We disagree with previous simulations by Niemeyer & Hillebrandt (1995), who reported on a possible loss of flame stability at a fuel density of $5 \times 10^7 ~{\rm g} ~{\rm cm}^{-3}$. Therefore, a self-turbulization and subsequent drastic acceleration of the flame seems unlikely - at least in the case of propagation into quiescent fuel. No convincing indications for active turbulent combustion (Kerstein 1996; Niemeyer & Woosley 1997) were observed in our simulations. However, flame propagation into a turbulent flow field may reveal completely different flame structures, since one expects strong enough turbulent eddies to disrupt the cellular flame pattern. The arguments above provide strong hints into this direction.

(ii)
Large-scale SN Ia simulations as performed by Reinecke et al. (2002b) apply a subgrid-scale model to determine the contribution from unresolved scales to the effective flame propagation velocity. However, as a lower cut-off they use the laminar burning velocity $s_{\rm l}$. On the basis of the performed simulations we suggest a slightly higher value in order to take into account the cellular structure of the flame below the Gibson scale. This may have impact on the nucleosynthesis yields from SN Ia models and will be addressed in a forthcoming study.

Acknowledgements
This work was supported in part by the European Research Training Network "The Physics of type Ia Supernova Explosions'' under contract HPRN-CT-2002-00303 and by the DFG Priority Research Program "Analysis and Numerics for Conservation Laws'' under contract HI 534/3. A pleasant atmosphere to prepare this publication was provided at the workshop "Thermonuclear Supernovae and Cosmology'' at the ECT*, Trento, Italy. We would like to thank M. Reinecke, S. Blinnikov, and W. Schmidt for stimulating discussions. The numerical simulations were performed on an IBM Regatta system at the computer center of the Max Planck Society in Garching.

References


Copyright ESO 2004