Vladimir A Sabelnikov

Andrei N Lipatnikov

Nilanjan Chakraborty

Shinnosuke Nishiki

Tatsuya Hasegawa

A Balance Equation for the Mean Rate of Product Creation in Premixed Turbulent Flames

Keywords: premixed turbulent combustion, mean reaction rate, flame stretch, DNS, modeling

Transport equations for reaction rate W and its Favre-averaged value W are derived from first principle in the case of premixed turbulent combustion. The assumptions made for derivation hold for unity Lewis number premixed flames at least in the flamelet regime of turbulent burning. Analysis of the latter equation shows that it involves two dominant terms, but the difference between them vanishes if reaction zones retain the structure of the zone in the unperturbed laminar flame. However, in such a case, turbulent burning velocity cannot grow with time during interaction of an initially laminar flame with a turbulent flow. Therefore, the analysis indicates a vital role played by local perturbations of reaction zone structure in premixed turbulent combustion. The dominance of these two terms and the important role played by the difference between them are confirmed by analyzing three DNS databases associated with both the corrugated flamelets and thin reaction zones regimes of premixed turbulent burning. Moreover, the DNS data show that perturbations of local displacement speed due to perturbations of local flamelet structure are also of paramount importance for modeling transport of flame surface density even in weakly turbulent flows. Finally, by simulating curved and/or strained 1 laminar premixed flames and integrating the transport equation for W across the flames, the integral is shown to depend linearly on the stretch rate even in highly perturbed flames, with results obtained from variously stretched flames being close to each other. Based on this finding, the difference between the two dominant terms in the transport equation for the mean rate W is hypothesized to depend linearly on the stretch rate conditioned to the reaction zone. Application of this hypothesis to the DNS data associated with the corrugated flamelets combustion regime yields encouraging results, thus, confirming a crucial role played by local perturbations of reaction zone structure even in weakly turbulent flames.

Introduction

As reviewed elsewhere [START_REF] Peters | Turbulent Combustion[END_REF][START_REF] Lipatnikov | [END_REF][3][START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF][START_REF] Lipatnikov | [END_REF][6][START_REF] Chakraborty | Turbulent Premixed Flames[END_REF][START_REF] Lipatnikov | Fundamentals of Premixed Turbulent Combustion[END_REF], a variety of models have been proposed to evaluate Favre-averaged reaction rate ρ W ≡ ρW in premixed turbulent flames in the case of a single-step chemistry. Some of these models invoke an extra balance equation for (i) the mean scalar dissipation rate (SDR) ρ χ ≡ 2ρD∇c • ∇c [START_REF] Chakraborty | Turbulent Premixed Flames[END_REF][START_REF] Mantel | [END_REF], (ii) mean Flame Surface Density (FSD) Σ [10][11][12][13], or (iii) variance ρc 2 used to close presumed Probability Density Function (PDF) P (c, c, ρc 2 /ρ) [14,15]. Here, ρ is the density, c is the combustion progress variable, D is the molecular diffusivity, W is the rate of product creation, and c ≡ c -c. Mantel and Borghi [START_REF] Mantel | [END_REF] converted their mean SDR transport equation to a transport equation of ρW by (i) utilizing the following expression ρW = ρ χ/(2c m -1) derived by Bray [16] for high Damköhler number (i.e. Da 1) flames and (ii) considering c m = 1 0 ρW cf (c)dc/ 1 0 ρW f (c)dc to be a constant. Here, f (c) is the burning mode probability density function [16]. It is worth noting, however, that the latter transport equation is a model equation, which is restricted to Da 1 and does not include e.g. dilatation effects [START_REF] Chakraborty | Turbulent Premixed Flames[END_REF] or spatial variations of c m , which are well pronounced even at high Da [17]. The goal of the present work is to derive and analyze an exact transport of W .

In the next section, an unclosed transport equation for W is derived and discussed. In the third section, numerical simulations are summarized. In the fourth section, results are reported, followed by conclusions.

Transport Equation

Derivation

Let us assume that the state of the mixture is characterized with a single combustion progress variable c, e.g. ρ = ρ(c), i.e. the flame is adiabatic, the Mach number is low, the Lewis number Le = 1, etc. [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]. Let us also assume that the rate W depends only on c, i.e. W = W (c). Then, using the continuity

∂ρ ∂t + ∇ • (ρu) = 0, (1)
and c-transport

ρ ∂c ∂t + ρu • ∇c = ∇ • (ρD∇c) + ρW (2)
equations, we arrive at

∂ ∂t (ρW ) + ∇ • (ρuW ) = ρ ∂W ∂t + ρu • ∇W = dW dc ρ ∂c ∂t + ρu • ∇c = dW dc ∇ • (ρD∇c) +ρW dW dc = ∇ • ρD dW dc ∇c -ρD∇c • ∇ dW dc + ρW dW dc = ∇ • (ρD∇W ) -ρN d 2 W dc 2 + ρW dW dc , (3)
where t is time, u is the flow velocity vector, and N ≡ D∇c•∇c. Equation (3) has a standard structure, i.e. unsteady and convection terms on the Left Hand Side (LHS) and diffusion and sink/source terms on the Right Hand Side (RHS). These sink/source terms change their sign within a flame, because the

sign of W = d 2 W/dc 2 or W = dW/dc changes with c. ∂ ∂t ρ W + ∇ • ρuW = ∇ • ρD∇W -ρN W + ρW W (4)
Equation ( 3) can also be filtered for LES. Equation ( 4) has been derived invoking a single assumption that the rate W depends solely on c if Le = 1 and this assumption is expected to be valid for unity

Lewis number premixed flames at least in the flamelet regime of premixed turbulent burning. Such an assumption is widely used by the theory of stretched laminar flames [18,19] and by various models of premixed turbulent burning [START_REF] Lipatnikov | [END_REF][3][START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF][START_REF] Lipatnikov | [END_REF][6][START_REF] Chakraborty | Turbulent Premixed Flames[END_REF][START_REF] Lipatnikov | Fundamentals of Premixed Turbulent Combustion[END_REF][START_REF] Mantel | [END_REF][10][11][12][13][14][15].

Discussion

First, application of Eq. ( 3) to an unperturbed planar 1D laminar premixed flame, followed by integration along a normal to the flame results in

∞ -∞ ρN W -ρW W dx = 0. ( 5
)
Therefore, the source/sink terms exactly balance each other after integration along a normal to an unperturbed laminar premixed flame.

Second, application of Eq. ( 4) to a statistically planar, 1D turbulent premixed flame yields

∂ ∂t ρ W + ∂ ∂x ρuW T 1 = ∂ ∂x ρD ∂W ∂x T 2 -ρN W T 3 + ρW W T 4 . (6)
Integration of this equation along x results in

ρ u ∂U t ∂t = ∞ -∞ ρW W -ρN W dx. (7)
Here,

U t ≡ ρ -1 u ∞ -∞
ρ W dx is the turbulent burning velocity. Equation [START_REF] Chakraborty | Turbulent Premixed Flames[END_REF] proves that the integrated terms T 3 and T 4 should not balance each other in order for turbulent burning velocity to develop, e.g. starting from U t = S L at t = 0. Accordingly, the sum of T 3 and T 4 should not vanish in a developing turbulent flame. In a fully-developed flame, T 3 and T 4 are unlikely to locally balance each other either.

Indeed, if (T 3 + T 4 ) = 0, then, the stationary Eq. ( 7) reduces to ρuW = ρD ∂W ∂x , but a non-trivial solution to the latter equation is hardly expected if e.g. ρD =const (this simplification is widely used in theoretical studies of premixed flames [18,19]). As will be discussed later local imbalance between terms T 3 and T 4 is associated with perturbations of the local structure of reaction zones by turbulent stretching.

Third, in the corrugated flamelets regime associated with large Damköhler Da = τ t /τ c and low Karlovitz Ka = (δ L /η) 2 numbers, the probability γ of finding intermediate (0 < c < 1) states of the reacting mixture can be assumed to be low everywhere within the mean flame brush and the following bimodal PDF [16] P

(x, t, c) ≈ [1 -c(x, t)]δ(c) +c(x, t)δ(1 -c) + γ(x, t)P f (x, t, c) (8)
can be invoked. Here, τ t = L/u and τ c = D u /S 2 L are turbulence and flame time scales, respectively, u is the rms turbulent velocity, L and η = LRe -3/4 t are integral and Kolmogorov length scales, respectively, δ L = D u /S L is the laminar flame thickness, Re t = u L/ν u is the turbulent Reynolds number, ν u is the kinematic viscosity of unburned gas, δ(1 -c) and δ(c) are Dirac delta functions. For brevity, dependencies of P , c, γ, etc. on spatial coordinates x and time t will not be specified in the following.

Accordingly,

ρ W = ρW = 1 0 ρ(c)W (c)P (c)dc = γ 1 0 ρ(c)W (c)P f (c)dc = γρ u W * , (9)
where W * is a scale for the local rate W (c). In order for the term ρ W to play a role in the ensembleaveraged Eq. ( 2), this source term should be comparable with the LHS terms, i.e. the product of γW * should be on the order of U t /δ t , where δ t is the mean flame brush thickness. The use of the same W * and Eq. ( 8) in order to estimate various terms in Eq. ( 6) shows that T 3 and T 4 are on the order

of γρ u W 2 * ∝ γ -1 ρ u (U t /δ t ) 2 ρ u (U t /δ t ) 2
, whereas other terms are on the order of ρ u (U t /δ t ) 2 or less.

Therefore, in the corrugated flamelets regime, T 3 and T 4 should dominate.

Fourth, in a general case, we can write

P (c) = αδ(c) + βδ(1 -c) +γ 1 P f,1 (c) + γ 2 P f,2 (c), (10)
where coefficients α, β, γ 1 , γ 2 and functions

P f,1 (c), P f,2 (c) are unknown, but P f,1 (c > c w ) = P f,2 (c < c w ) = 0
, with c w being associated with the boundary of the reaction zone, i.e. W (c < c w ) can be disregarded, but the rate W is substantial at c > c w . Recent experimental data reviewed elsewhere [START_REF] Lipatnikov | Fundamentals of Premixed Turbulent Combustion[END_REF]20] indicate that reaction zones are thin even in the thin reaction zones regime. In such a case, we can assume that γ 2 1 and arrive at

W = O(γ 2 W * ) = O(U t /δ t ), T 1 = O[ρ u (U t /δ t ) 2 ], T 2 = O[ρ u D u U t /δ 3 t ], T 3 = O[ρ u γ -1 2 (U t /δ t ) 2 ],
and

T 4 = O[ρ u γ -1 2 (U t /δ t ) 2 ]
. Therefore, we can expect that terms T 3 and T 4 dominate even in the thin reaction zones regime.

Fifth, if we come back to the corrugated flamelets regime and disregard perturbations of local flamelet structure, then, P f (c) in Eq. ( 8) can be modeled as follows [15,16], where ε 1 is sufficiently small in order for contributions of ranges of 0 < c < ε and 1 -ε < c < 1 to the mean terms to be negligible. Accordingly,

P f (ε < c < 1 -ε) ∝ (δ L |∇c|) -1
ρW W -ρN W ∝ γ 1-ε ε ρW W -ρN W dc δ L |∇c| = γ ξ 2 ξ 1 ρW W -ρN W dξ = 0 ( 11
)
by virtue of Eq. ( 5). Here, ξ is the non-dimensional spatial coordinate locally normal to a flamelet and integration is performed over a single flamelet.

However, the equality of (T 3 + T 4 ) = 0, resulting from Eq. ( 11), is wrong. Indeed, if (T 3 + T 4 ) = 0, then, ∂U t /∂t = 0 due to Eq. ( 7) and the burning velocity cannot grow even if the initial U t = S L .

Thus, an assumption that local flamelet structure perturbations can be neglected, which is equivalent to substitution of P f ∝ (δ L |∇c|) -1 into the BML PDF given by Eq. ( 8), results straightforwardly in Eq. ( 11) and, hence, in a wrong conclusion that U t cannot develop. The assumption appears to be wrong even in statistically stationary case, because, otherwise, Eq. ( 6) would involve solely transport terms, but would not involve a term that controls an increase in W from W = 0 in unburned gas.

Thus, local flamelet structure perturbations cannot be neglected, i.e. they play a vital role in premixed turbulent burning. This direct consequence of the introduced W -transport equation shows that it offers an opportunity to gain fundamental insights into the physics of turbulent flames.

At a first glance, the conclusion regarding an important role played by local flamelet structure perturbations in premixed turbulent burning could seem to be trivial. Indeed, if the effect of turbulence on premixed combustion is associated with stretching and wrinkling of flamelets (or thin reaction zones) by turbulent eddies, then, the aforementioned perturbations appear to be taken into account. However, this is not true, because many models that address the flamelet stretching and wrinkling place the focus of consideration on an increase in the flamelet surface area by turbulent stretching, but, at least in the case of sufficiently weak turbulence, neglect the flamelet structure perturbations caused by the same turbulent stretching.

For instance, FSD models [3,[START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF] deal with the following unclosed transport equation [10][11][12][13] ∂Σ ∂t

+ ∇ • uΣ = a t Σ -∇ • S d nΣ + S d Σ∇ • n (12)
for a mean FSD Σ = |∇c|. Here, a t = ∇ • unn : ∇u is the strain rate, n = -∇c/|∇c|, and

S d = 1 ρ|∇c| [∇ • (ρD∇c) + ρW ] ( 13
)
is the displacement speed. Equation ( 12) is exact and addresses both the increase in flamelet surface area by turbulent stretching, see the first term on the RHS, and the flamelet structure perturbations, because ρS d (x, t) = ρ u S L in a general case. However, at a first glance, Eq. ( 12) does not seem to necessitate ρS d (x, t) = ρ u S L , but appears to be consistent with a simplification1 of ρS d (x, t) = ρ u S L [21]. A similar simplification of (ρS d ) s = ρ u S L is widely used to close Eq. ( 12), as reviewed elsewhere [3,[START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]. Thus, contrary to Eq. ( 4) derived above, the FSD Eq. ( 12) or, to the best of the present authors' knowledge, another model of turbulent flames does not reveal a crucial role played by flamelet structure perturbations in premixed turbulent combustion.

In the next sections, we shall illustrate certain features of Eq. ( 4) by processing DNS data.

Numerical Simulations

DNS attributes

We evaluated various terms in Eq. ( 4) by processing three DNS databases computed by Rutland

and Cant [START_REF] Rutland | Proc. of 1994 Summer Program[END_REF], by Nishiki et al. [START_REF] Nishiki | [END_REF]24] in Nagoya University, and by Chakraborty et al. [25,26] in Newcastle University. In the following, these three studies will be called RC, Na, and Ne DNS, respectively. Because the DNS data were discussed in detail elsewhere [START_REF] Rutland | Proc. of 1994 Summer Program[END_REF][START_REF] Nishiki | [END_REF][24][25][26] and were already used in a number of investigations, see [27][28][29][30][31][32][33][34][35], [17,35,[36][37][38][39][40][41][42][43], [44][45][46][47][48] for the RC, Na, or Ne DNS, respectively, we will restrict ourselves to a brief summary of the simulations.

In all three cases, unsteady 3D balance equations for mass, momentum, energy, and mass fraction of the deficient reactant were numerically solved. The ideal gas state equation was used. Combustion chemistry was reduced to a single reaction. The Lewis and Prandtl numbers were equal to 1.0 and 0.7, respectively. Other flame characteristics are reported in Table 1, where Da th = LS L /(u δ th ), Rectangular computational domains (8×4×4 mm, 36δ th × 45δ th × 45δ th , and 36.2δ th × 24.1δ th × 24.1δ th in the Na, RC, and Ne DNS, respectively) were resolved using uniform Cartesian meshes of 512×128×128, 261×128×128, and 345×230×230 points, respectively. In all cases, the mean flow velocity was parallel to the x-axes and normal to the mean flame brush, with the periodic boundary conditions being set at the transverse sides.

Ka th = (u /S L ) 3/2 (L/δ th ) -1/2 , δ th = (T b -T u )/
In all cases, homogeneous isotropic turbulence was used to initialize velocity fluctuations and a single planar laminar flame was embedded into the computational domain at t = 0. In the RC and Ne DNSs, turbulence decayed with time and the mean inlet velocity U = S L was constant. Averaging was performed over transverse planes at t/τ t = 4 (RC flame), 2 (Ne flame D), 3 (Ne flames C and E), and 4.34 (Ne flame B). In the Na DNS, homogeneous isotropic turbulence was generated in a separate box, was injected into the computational domain at x = 0, and decayed along the direction x. Averaging was performed over transverse yz-planes and over time interval on the order of 10 ms (about 200 snapshots) during that both U t (t) and mean flame brush thickness δ t (t) oscillated around statistically steady values [42].

These earlier DNSs did not aim at studying Eq. ( 4), which is introduced here. Accordingly, the used meshes may not be sufficiently fine to properly resolve spatial variations in W W and N W . Therefore, in the present paper, the DNS data are analyzed to gain a qualitative insight into the behavior and relative magnitudes of various terms in Eq. ( 4), whereas target-directed DNS with a very fine mesh is necessary to investigate these terms quantitatively.

Simulations of perturbed laminar flames

To gain an insight into a link between Eq. ( 4) and local perturbations of reaction zone, we also simulated stretched laminar flames by numerically integrating the following unsteady 1D balance equations

∂ ∂t ( Φ) + gΦ + 1 r k ∂ ∂r r k vΦ = 1 r k ∂ ∂r r k 1-m d φ ∂Φ ∂r + S φ , (14)
as discussed in detail elsewhere [49]. Here, Φ = {1, g, c},

d φ = {0, D u , D u }, S φ = {0, -g 2 + J 2 / , W }, = ρ/ρ u = 1/[1 + (σ -1)c] is the normalized density, g is rate of strain, k = {0,
Here, r m corresponds to the outlet boundary of the computational domain and J is an input parameter required to vary the strain rate. The initial conditions described a small pocket of combustion products.

The flame stretch rate was evaluated as follows

ṡ = g + k R f dR f dt , (16)
with the flame coordinate R f being obtained from a constraint of c[r = R f (t), t] = c * , where c * = c m is associated with the peak W (c). Such a method of defining flame surface is widely used in theoretical [18,19] and numerical [50] studies of laminar flames. Because the laminar flame thickness is small and almost constant in expanding flames, computed results are weakly sensitive to the choice of c * . Typical results of such simulations are reported in Fig. 1, which supports the simplest linear fit

1 R k f ∞ 0 ρW W -ρN W r k dr = ρ u V ṡ ( 17
)
for various types and strengths of flame perturbations. The dimensional coefficient V can be called "Markstein velocity" via an analogy with Markstein lengths used to characterize response of a laminar flame speed to a weak stretch rate.

DNS Results and Discussion

Various terms in Eq. ( 6) are shown in Figs. 2 and3 in cases associated with the corrugated flamelets and thin reaction zones regimes, respectively. Due to averaging over reaction zones, the term T 3 (or T 4 ) assumes positive (or negative) sign within the flame brush, contrary to the instantaneous predecessor term in Eq. ( 3) or (14), which changes its sign following the sign of W (or W ). As the probability of finding burning fluid is smaller than finding unburned or burned gas, the signs of W and W towards the burned gas side (i.e. 0.8 < c < 1) ultimately determine the sign of T 3 and T 4 , because W and W remain negligible on the unburned gas side (0 < c < 0.5). In all cases, terms T vanishes under conditions of the Na DNS and is small in other five cases.

The fact that (T 3 +T 4 ) does not vanish, but plays a substantial role highlights flamelet perturbations, as discussed in Sect. 2.2. An important role played by the perturbations even in weakly turbulent flames is supported in Fig. 4, which compares dependencies of ρuW on c, extracted from the DNS and obtained by numerically integrating the following equation

∂ ∂x ρuW = ∂ ∂x ρD ∂W ∂x + V S L ρ W ṡ|c w,1 < c < c w,2 . (18)
Here, V was evaluated by simulating the counterpart perturbed laminar flames, as discussed in Sect.

3.2, while all other terms on the RHS were extracted from the DNS. This equation can be obtained if (i) ρN W and ρW W are averaged invoking the PDF given by Eq. ( 8) with

P f ∝ (δ L |∇c|) -1 , (ii)
the probability γ is modeled as follows γ = ρ W δ L /S L by applying the same PDFs to averaging W , and (iii) Eq. ( 17) is used to integrate ρW W -ρN W over flamelets. Because W , ρN W , and ρW W vanish beyond the reaction zone, the local stretch rate ṡ is conditioned to c-values associated with the Because (i) importance of local perturbations of reaction zones is recognized for highly turbulent flames and (ii) the model Eq. ( 18) relies on the flamelet paradigm, we will restrict ourselves to assessment of this model using weakly turbulent Na DNS data, which pertain [17] to the corrugated flamelets regime. Moreover, these DNS cases deal with statistically stationary flames, and, thus, the suggested joint closure relation for (T 3 + T 4 ) can be tested using a well-resolved quantity such as ρuW , thereby, circumventing the aforementioned problem of eventually insufficiently fine numerical resolution of spatial variations in ρN W and ρW W in the case of the Na DNS. For the thin reaction zones regime [START_REF] Peters | Turbulent Combustion[END_REF],

the problem of closing (T 3 + T 4 ) requires further analysis.

Figure 4 shows that Eq. ( 18) predicts ρuW (c) in a reasonably satisfactory manner, thus, indicating an important role played by local flamelet structure perturbations even in weakly turbulent flames. If W = (ρuW -ρu W )/(ρũ) is evaluated using the model ρuW , with other terms being extracted from An important role played by flamelet structure perturbations is consistent with the FSD Eq. ( 12).

To show this consistency, we evaluated uΣ by numerically integrating Eq. ( 12), with all terms on the RHS being extracted from the Na DNS. Subsequently, Σ = (uΣ-u Σ )/ū was calculated using the DNS data on u Σ and ū. Results obtained by substituting S d computed using Eq. ( 13) into Eq. ( 12) agree very well with Σ extracted straightforwardly from the DNS, cf. solid and dotted lines in Fig. 5, thus, validating the simulations.2 However, if the correct Eq. ( 13) is substituted with an approximation of S d (x, t) = ρ u S L /ρ when extracting the two last terms on the RHS of Eq. ( 12) from the same DNS, the mean FSD is substantially overestimated, see dashed lines. In particular, Σ does not vanish at c = 1 due to the use of this approximation. Moreover, the burning velocity S L Λx 0 |∇c|dx is significantly (about 50 %) overestimated if S d is assumed to be equal to ρ u S L /ρ. This test shows that the approximation is wrong, i.e. flamelet structure perturbations significantly affect FSD transport even in weakly turbulent premixed flames, in line with the above discussion.

The fact that terms T 3 and T 4 dominate in Eq. ( 6), with (T 3 + T 4 ) playing a substantial role, is a challenge for closing this transport equation, because both T 3 and T 4 should be modeled with a high precision in order to predict (T 3 +T 4 ). A similar problem arises for the closure of the transport equation for the mean scalar dissipation rate χ at high Reynolds numbers, but models that deal with the latter equation are successfully developed, as reviewed elsewhere [START_REF] Chakraborty | Turbulent Premixed Flames[END_REF]. Moreover, a solution could consist of studying (T 3 + T 4 ) instead of modeling each term separately.

Figure 4 indicates that Eq. ( 18) is a reasonable first step to development of such a joint closure in the corrugated flamelets regime (Ka < 1). Difference between the DNS and model results could be attributed to inability of simple problems addressed in Sect. 3.2 to represent variety of unsteady local burning structures that exist in a premixed turbulent flame [52]. Such effects are more pronounced at Ka > 1 and substantial difference between DNS data and model prediction has been found for Ne-database. The difference can also be attributed to numerical resolution of these DNS cases, because a mesh sufficiently fine to resolve terms in Eq. ( 2) could be too coarse to resolve rapidly varying terms in Eq. (3).

Furthermore, Eq. ( 18) and Fig. 4 explain why (T 3 + T 4 ) changes its sign from positive to negative when c is increased, see dotted lines in Fig. 3. The point is that positive (negative) stretch rates ṡ|c w,1 < c < c w,2 statistically dominate at the leading (trailing) part of a premixed turbulent flame brush, at least under conditions of the present DNS.

Conclusions

Transport equations for instantaneous reaction rate and its mean value have been derived for premixed mode of burning under the assumption which is valid at least for the flamelet regime of combustion. The latter equation reveals a vital role played by local perturbations of reaction zone structure even in weakly turbulent flames. 12), with all other terms in this equation, including S d , being extracted from the DNS data. Dashed lines show SL|∇c| obtained by integrating Eq. ( 12), with S d (x, t) = ρuSL/ρ and all other terms being extracted from the DNS data.

  max |dT /dx| is the thermal laminar flame thickness, and σ = ρ u /ρ b is the density ratio. The RC and Na DNSs address the corrugated flamelets regime, whereas the Ne DNS data are associated with the thin reaction zones regime[START_REF] Peters | Turbulent Combustion[END_REF].

  1, 2} for planar, cylindrical, and spherical flames, respectively (g = 0 if k = 2), m is a power exponent of the temperature dependence of D = D u / m . The boundary conditions were set as follows v(0, t) = ∂c ∂r (0, t) = ∂g ∂r (0, t) = 0, ∂c ∂r (r m , t) = 0, g(r m , t) = J.

Figure 1 :

 1 Figure 1: The LHS of Eq. (17) vs. normalized stretch rate τc ṡ. Simulations were performed with σ, m, and W (c) used in case H of Na DNS. 1 -expanding spherical flame (k = 2, J = 0), 2 -expanding cylindrical flame (k = 1, J = 0), 3steady strained cylindrical flame (k = 1, various J), 4 -steady strained planar flame (k = 0, various J).

  3 and T 4 dominate in line with the above discussion. Moreover, |T 3 | differs from |T 4 |, with the difference magnitude, see dotted lines, being much smaller than |T 3 | or |T 4 |, but comparable with |T 1 | or |T 2 |. The unsteady term

Figure 2 :

 2 Figure 2: Various terms in Eq. (6) obtained from flames associated with the corrugated flamelets regime.

Figure 3 :

 3 Figure 3: Various terms in Eq. (6) obtained from flames associated with the thin reaction zones regime.

Figure 4 :

 4 Figure 4: Normalized flux ρuW Du/(ρuS 3 L ) vs. c. Solid and dashed lines show results extracted from the DNS and obtained by integrating Eq. (18), respectively. cw,1 = 0.55 and cw,2 = 0.99.

Figure 5 :

 5 Figure 5: Mean FSD multiplied with the laminar flame speed, computed in Na DNS cases H, M, and L. Dotted lines show SL|∇c| extracted from the DNS data. Solid lines show SL|∇c| obtained by integrating Eq. (12), with all other terms in this equation, including S d , being extracted from the DNS data. Dashed lines show SL|∇c| obtained by integrating Eq. (12), with S d (x, t) = ρuSL/ρ and all other terms being extracted from the DNS data.

Fig. 1 .

 1 Fig.1. The LHS of Eq. (17) vs. normalized stretch rate τc ṡ. Simulations were performed with σ, m, and W (c) used in case H of Na DNS. 1 -expanding spherical flame (k = 2, J = 0), 2 -expanding cylindrical flame (k = 1, J = 0), 3steady strained cylindrical flame (k = 1, various J), 4 -steady strained planar flame (k = 0, various J).

Fig. 2 .

 2 Fig.2. Various terms in Eq. (6) obtained from flames associated with the corrugated flamelets regime.

Fig. 3 .

 3 Fig.3. Various terms in Eq. (6) obtained from flames associated with the thin reaction zones regime.

Fig. 4 .

 4 Fig. 4. Normalized flux ρuW Du/(ρuS 3 L ) vs. c. Solid and dashed lines show results extracted from the DNS and obtained by integrating Eq. (18), respectively. cw,1 = 0.55 and cw,2 = 0.99.

Fig. 5 .

 5 Fig. 5. Mean FSD multiplied with the laminar flame speed, computed in Na DNS cases H, M, and L. Dotted lines show SL|∇c| extracted from the DNS data. Solid lines show SL|∇c| obtained by integrating Eq. (12), with all other terms in this equation, including S d , being extracted from the DNS data. Dashed lines show SL|∇c| obtained by integrating Eq. (12), with S d (x, t) = ρuSL/ρ and all other terms being extracted from the DNS data.

Table 1 :

 1 Flame characteristics Case u /S L L/δ th σ Re t Da th Ka th

	RC	1.4	9.6	3.3 57	6.8	0.54
	Na-H	0.9	15.9	7.5 96	18.0	0.21
	Na-M	1.0	18.0	5.0 96	17.8	0.24
	Na-L	1.3	21.8	2.5 96	17.3	0.30
	Ne-B	6.3	1.4	5.5 24	0.23	13.0
	Ne-C	7.5	2.5	5.5 48	0.33	13.0
	Ne-D	9.0	4.3	5.5 100 0.48	13.0
	Ne-E	11.3	3.75	5.5 110 0.33	19.5

Later, we will show that the use of such a simplification can result in Σ(c = 1) = 0 if Σ(c = 0) = 0.

The fact that these DNS cases are well resolved to analyze the FSD transport is also supported by small magnitudes of residuals obtained by extracting each term of Eq. (12) from the DNS, as discussed in[51], where the relation between S d and SL is also assessed straightforwardly.

Acknowledgments

VS gratefully acknowledges the financial support by ONERA. AL gratefully acknowledges the financial support by CERC and Chalmers Transport Area of Advance. NC gratefully acknowledges the financial support provided by EPSRC, UK.


W Tigra

F Lucas

N Benjamin

G David

P Antonio

F Emerson

L Violaine

G Anthony

C Azevedo Coste

TOWARDS FES-ASSISTED GRASPING CONTROLLED BY RESIDUAL MUSCLE CONTRACTION AND MOVEMENT ON PERSONS WITH TETRAPLEGIA

Introduction/Background

Functional electrical stimulation (FES) can be used on individuals with upper motoneuron dysfunctions to restore grasping functions. To be functional and useful in daily tasks, the patient must be able to pilot the device by means of an interface. We have investigated two techniques that could address the situation of persons with tetraplegia: electromyography (EMG) and inertial measurement units (IMU), respectively measuring muscle activity and limb movements.

Material and Method

One group of 8 tetraplegic patients (EMG group) was equipped with two pairs of EMG electrodes located on muscles of the upper arm. Selected muscles could be voluntary activated in a comfortable way. Another group of 9 tetraplegic patients (IMU group) was equipped with one IMU (3-axis accelerometer and 3-axis gyroscope) located on the arm. In the EMG group a threshold detection algorithm was used to detect muscle contractions while in IMU group a classification algorithm was used to detect two distinct movements of the limb where the IMU was located. The detection algorithms outputs were used to trigger two postures of a robot hand. The same outputs were also used, for patients where FES of arm muscles lead to visible hand movements, to trigger the contractions of two muscles on their contralateral lower arm.

Results

All the patients were able to control muscle contraction or limb movements to trigger different actions. Patients were asked to control the robotic hand, following a randomized sequence of two pre-defined postures. In EMG group the success score was of 95% and 92% in IMU group. Whenever FES was used, patients could activate grasping movements and wrist extension. Three patients were able to perform functional tasks such as grasping and relocating objects.

Conclusion

Despite little training, all the patients were able to voluntarily control a robotic hand and, when it was possible to test, their own hand with FES.

Keywords:

Movement analysis Electromyography interface Inertial Measurement Unit interface


Claudia Colonnello

Miroslav Kramár

Dynamics of silo deformation under granular discharge

We use Topological Data Analysis to study the post buckling behavior of laboratory scale cylindrical silos under gravity driven granular discharges. Thin walled silos buckle during the discharge if the initial height of the granular column is large enough. The deformation of the silo is reversible as long as the filling height does not exceed a critical value, Lc. Beyond this threshold the deformation becomes permanent and the silo often collapses. We study the dynamics of reversible and irreversible deformation processes, varying the initial filling height around Lc. We find that all reversible processes exhibit striking similarities and they alternate between regimes of slow and fast dynamics. The patterns that occur at the beginning of irreversible deformation processes are topologically very similar to those that arise during reversible processes. However, the dynamics of reversible and irreversible processes is significantly different. In particular, the evolution of irreversible processes is much faster. This allows us to make an early prediction of the collapse of the silo based solely on observations of the deformation patterns.

I. INTRODUCTION

Thin walled cylindrical silos are used extensively as storage facilities for granular materials in industrial applications. Despite great efforts to develop appropriate design and construction protocols, their structural failure remains a widespread problem [? ? ? ], leading to important economic losses and risks for human personnel.

One of the main factors hindering the efficient prevention of silo failure is the difficulty in accounting for the interaction between the grains and the silo [? ], due to the complex stress patterns that arise within confined granular media (e.g [? ? ? ]). The model proposed by Janssen [? ] provides a simple description of the stress on a silo filled with a granular material and is widely used as a reference in engineering applications [? ? ]. It accurately describes the total stress on the walls and bottom of a silo if local contact shear forces exerted by the grains are in a state of maximal mobilization [? ? ? ]. However, the state of mobilization is strongly affected by the details of the granular packing and the filling procedure [? ? ]. During a granular discharge of the silo the maximal mobilization condition is attained spontaneously if it is not initially satisfied [? ]. Once the maximal mobilization condition is met, the total stress during the discharge is well described by the model [? ? ].

Another limitation of Janssen's model is that it fails to predict local stress fluctuations on the walls of the silo, as it is based on a description of the granular material as a continuum. Local stress measurements and simulations show that fluctuations can be very large, particularly at the onset of a granular discharge [? ? ]. Moreover, their distribution may not be symmetric in the angular direction [? ? ]. Such stress fluctuations are likely to induce geometric defects and load eccentricities on the silo wall, that can affect the stability of the structure. In fact, the critical buckling load of a thin empty shell exhibits a strong imperfection sensitivity [? ? ], which represents an enduring challenge for the accurate description of the shell's stability [? ? ? ].

The most common failure mode for silos with thin and flexible vertical walls is buckling due to the effective axial compression resulting from frictional forces exerted by the grains [? ]. In particular, if the initial filling height of the silo is large enough, its vertical walls buckle during the gravity driven discharge of the granular material [? ? ]. During the discharge, diamond shaped indentations develop at the lower end of the cylindrical shell and form intricate deformation patterns (see Fig. ??). The number and size of the indentations are strongly affected by local stress fluctuations due to granular collisions and they tend to increase with the initial filling height of the silo. If the filling height exceeds a critical value, L c , then the indentations grow large enough to produce a permanent deformation of the shell. This is usually followed by a catastrophic collapse of the structure. Otherwise, all indentations disappear completely by the end of the granular discharge. The scaling behavior of L c has been studied systematically in [? ], where it is successfully described by a model that considers the stability of an empty shell subject to an effective axial stress derived from Janssen's model. However, the evolution of the deformation patterns is far from being understood.

Local features of the stress distribution on the silo wall, such as fluctuations and friction mobilization, are expected to play a fundamental role in the deformation process [? ]. However, direct local stress measurements on thin walled silos can be challenging. A grid of strain sensors has been used on large scale experiments [? ? ? ], but such methods are difficult to implement at laboratory scale, and could affect significantly the response of very thin shells. These difficulties have been circumvented by considering the deformation patterns that develop on the shell as it interacts with the granular material [? ? ]. For small shell deformations, this allows for the reconstruction of the stress field on the shell.

In this paper we present a quantitative study of the dynamics of post buckling deformation patterns induced by granular discharges on thin walled cylindrical silos. We use images obtained from laboratory scale experiments using paper silos, and employ Topological Data Analysis (TDA) to study the patterns that develop on the silo while varying the initial filling height around the critical value, L c . We use persistent homology[? ? ] to represent a deformation process as a trajectory in the space of persistence diagrams. This space is well suited to capture important topological information about the pattern and quotient out the overwhelming geometrical variations [? ]. In this way, we compare different deformation processes and gain novel insights into their evolution, that contribute to our understanding of the conditions that lead to the failure of the structure.

We find that deformation patterns at early stages of irreversible processes leading to the collapse of the silo are very similar to the patterns exhibited by reversible processes. However, already at these early stages their dynamics is significantly different. In particular, irreversible processes evolve faster. This allows us to effectively predict the collapse at an early stage of the discharge. The existence of such early differences between reversible and irreversible processes indicates that the collapse is triggered by partially mobilized grain-towall friction, as proposed in [? ? ]. We also observe evidence of slow-fast dynamics in the evolution of the deformation and argue that it is caused by the intrinsic nonlinear behavior of the thin shell.

The paper is organized as follows. Section ?? gives a detailed account of our methods. In Sec. ?? we present our results. First, in Sec. ?? we introduce a set of qualitative criteria to classify a deformation process after the discharge. We define three classes according to the deformation sustained by the silo. In Secs. ?? and ?? we analyze the dynamics of deformation processes in each of these classes, and discuss the evidence of fast-slow dynamics as well as the reversibility of the trajectories. We compare the trajectories for all deformation processes in Sec. ?? and analyze quantitative differences and similarities between trajectories in each class. In Sec. ?? we discuss the early classification of the processes. Finally, in Sec. ?? we discuss the origin of the fast-slow dynamics. We conclude with some final remarks in Sec. ??.

II. METHODS

A. Experimental setup

Figure ??(a) shows a sketch of the experimental setup. Silos used in our experiments are built following the protocol described in [? ]. A cylindrical paper shell, with thickness h = (24 ± 5) µm and radius R = (4.0 ± 0.1) cm, is fitted tightly on a solid acrylic base fixed to a level surface. The upper boundary of the silo remains free. The solid base has a centered orifice of diameter a = (0.5 ± 0.1) cm. We start by closing the orifice with a foam plug. Then the silo is filled at constant flow rate with glass beads of diameter d = (1.50 ± 0.05) mm, through a fixed funnel centered above the silo. After filling the silo, the plug is removed and the gravity driven granular discharge starts.

The critical filling height for our silos (measured as described in [? ]) is L c = (11.5 ± 0.6) cm. To estimate the effect of fluctuating deformations induced by the granular flow, we perform one granular discharge with an initial filling height L f o = (7.0 ± 0.3) cm. This filling height is sufficiently below L c to guarantee that the silo wall only exhibits deformations due to collisions with the moving particles. The remaining granular discharges are performed with initial filling heights L o > 0.8L c , so that persistent deformation patterns develop as well on the silo wall. We study five different silos and a total of 27 granular discharges. We progressively increase L o in successive granular discharges of a single silo until an evident permanent deformation develops on its wall.

In order to record the configuration of the deformation patterns, the silo is obliquely illuminated by a strong light source from above. In this way, the typical diamond shaped indentations are clearly visible as pairs of approximately symmetric dark and bright features, see Fig. ??(a). A digital camera placed in front of the silo is used to capture a sequence of images during the granular discharge. The sampling starts simultaneously with the discharge and continues for 8.3 s. The duration of the sampling is set by limited hardware resources and does not encompass the whole granular discharge. However, this has only a minor impact on our results, as our main interest lies in the emergence and early development of the deformation. Moreover, if a permanent deformation occurs it is always before 8.3 s and we truncate the sequence shortly after we detect it by visual inspection.

Grayscale digital images with size 440 × 400 pix are collected at a fixed rate of 30 frames per second, i.e. we take an image at every time t i = i/30 s for 0 ≤ i ≤ 250. Images capture a fixed area of the shell spanning roughly 100 o of its circumference, as depicted in Fig. ??(a).

B. Image processing

To study the topological features of the deformation patterns and their evolution during a granular discharge, we construct a time series of dimensionless scalar fields {ω i } for each discharge. A piecewise constant scalar field ωi is obtained from the digital image of the shell taken at time t i (see Sec. ??). In the rest of this section we will describe how to construct ω at a fixed time from an image of the shell. Then we will show that ω retains relevant information about the topology of the deformation pattern exhibited by the shell at this time.

First, the intensity field of the digital image is normal-ized by a median filtered (size=50pix) reference image taken right before the onset of the discharge. This compensates for inhomogeneities in the illumination of the shell partly caused by its curvature. The reference image may exhibit some incipient deformations produced during the filling of the silo but their size is marginal and their effect on the median filtered intensity field is negligible. The normalized intensity field defined on the pixels of the original image is denoted by I. If a point (x, y) in the image belongs to a dark feature, corresponding to the upper part of a diamond shaped indentation, then I(x, y) < 1. Conversely, for I(x, y) > 1 the point corresponds to the bottom part of the indentation. We retain only information about the upper part of the indentations to simplify the procedure and define ω(x, y) = min{I(x, y), 1} -1.

Thanks to the characteristic symmetry of the indentations, considering only their upper part is sufficient, as we will show below. Finally, we apply an erosion filter on ω(x, y) to remove small noisy features produced by the texture of the paper. An example of an original image of the pattern and the corresponding field ω resulting from this procedure are shown in Figs. ??(a) and (b), respectively.

The shape of the deformed shell is completely described by the radial displacement from its original cylindrical state, ω. The fields ω and ω are far from being identical. However, we will show that their similarities within the set Ω0 = {(x, y) : ω(x, y) < 0} make ω a good descriptor of the topology of the deformation pattern.

The definition of ω and the illumination of the silo ensure that for a paper sheet folded to recreate the Yoshimura pattern [? ], there exists a monotone function f : R → R such that |f (ω(x, y)) -ω(x, y)| ≤ 0.07f (ω(x, y)) for (x, y) ∈ Ω0 . Figures ??(c)-(d) show the fields ω and ω on a cross-section of an indentation from the Yoshimura pattern. The folded shell is made with the same paper and has the same radius as the silos used in our experiments. The resulting regular indentations are very similar in size and shape to those that develop during the discharge. Therefore, we expect that in our experiments the differences between ω and f (ω) have the same bound as in the case of the Yoshimura pattern.

Finally, we show that ω is a good descriptor of the topology of the deformation pattern. In particular, consider the set of indentations given by {(x, y) : ω(x, y) < 0}. The geometry of this set differs considerably from the geometry of Ω0 but their topology is closely related. For every diamond shaped indentation, its upper part is contained in Ω0 , so for every indentation there is a connected component in Ω0 . Moreover, if two indentations merge they always do so along an edge, as the vertices exhibit an increased stiffness thanks to their increased curvature [? ]. Thus, the connected components in Ω0 corresponding to these indentations merge into a single one as well.

These observations support the existence of a close relation between the topology of the sub-level set Ωγ = {(x, y) : ω(x, y) < γ} and the set {(x, y) : ω(x, y) < γ}. This justifies our choice to study the topological features of the pattern represented by ω.

C. Dynamics in the Space of Persistence Diagrams

Persistent homology [? ? ] is used to encode topological features of the pattern given by a scalar field ω. In particular, two persistence diagrams PD 0 (ω) and PD 1 (ω) are assigned to ω . These diagrams describe changes in topology of the sub-level sets Ωγ , for all γ ∈ R. Every point (b, d) ∈ PD 0 (ω) (PD 1 (ω)) encodes the appearance of a connected component (loop) of Ωγ at γ = b and its disappearance at γ = d. The lifespan of (b, d), defined by d -b, is often interpreted as a measure of prominence of the topological feature represented by (b, d). For an in-depth presentation describing the connection between the persistence diagrams and the underlying pattern we refer the reader to [? ].

More classical measures, such as number and maximal depth of indentations, are readily available from the persistence diagrams. Due to the characteristic diamond shape of buckling indentations, even the total area they cover can be estimated. Conversely, some of our results can be obtained by analyzing the evolution of the area covered by indentations during the discharge.

There are well defined notions of distance between pairs of persistence diagrams [? ]. We will only use the so called Wasserstein distance, d W 2 . This distance measures the overall difference between the persistence diagrams, encoding the difference between the patterns they represent. As such, it measures changes in both the number and depth of the indentations. Larger changes in depth contribute more to the distance between the corresponding diagrams.

The set of all persistence diagrams equipped with the d W 2 distance is a complete metric space denoted by Per. This, combined with the fact that the persistence diagrams are good descriptors of patterns, makes Per a good choice of an observation space. For a given discharge, a trajectory X : [0, ∞) → Per describes the evolution of the observed deformation pattern. By using the methods presented above we produce a sample of this trajectory

X(t i ) = (PD 0 (ω i ), PD 1 (ω i )), at times t i = i/30 s, for 0 ≤ i ≤ 250.

III. RESULTS AND DISCUSSION

A. Qualitative classification

Observation of the pattern through the granular discharge allows for a qualitative classification of the deformation processes into three classes: (i) Far-from-collapse (FFC), (ii) Close-to-collapse (CTC), and (iii) Collapse (C) processes. A collapse process is one that results in a permanent deformation of the shell, detectable by visual inspection after the granular discharge has concluded. These processes often lead to a catastrophic collapse of the silo. They occur if L o > L c and always produce indentations that grow until they exceed the elastic limit of the shell.

Conversely, if L o < L c a noncollapse deformation process occurs, i.e., all visible deformations completely disappear by the end of the discharge. We say a noncollapse process is far-from-collapse if all indentations remain comparable in size to the fluctuating deformations induced by the granular flow. Since these deformations are produced by collisions of the grains against the silo wall they exhibit small values of the radial displacement of the shell (ω 10 -2 R). Nevertheless, large clusters of localizations may develop during a FFC process. If a noncollapse process exhibits indentations that attain a larger radial displacement before receding, we say it is a close-to-collapse process. Figure ?? shows characteristic examples of the deformation pattern for each class. In the rest of this paper we will demonstrate quantitative differences and similarities between the three classes.

B. Dynamics of individual deformation processes

In this section we consider the trajectories (in Per) of typical deformation processes in each of the qualitative classes introduced in Sec. ??. For the sake of clarity we only consider trajectories corresponding to three successive discharges of the same silo, with filling heights L a o = (11.2 ± 0.3) cm (FFC), L b o = (11.5 ± 0.3) cm (CTC), and L c o = (11.6 ± 0.3) cm (C). In Sec. ??, we show that trajectories of processes in the same class exhibit strong quantitative similarities. To simplify the notation, we denote the trajectories of the representative FFC, CTC, and C processes by X a , X b , and X c , respectively.

We start by studying the deviation of a trajectory X from its initial condition X(0), which represents the marginal deformation pattern sustained by the shell at the onset of the discharge. Figure ??(a) shows that for all three trajectories the distance from the initial condition, Dev(X(t)) := d W 2 (X(0), X(t)), initially increases. This reflects the growth and proliferation of indentations at the beginning of the discharge.

As expected, the deviation from the initial condition is smallest for the FFC trajectory X a and largest for the C trajectory X c . Notice that both Dev(X a (t)) and Dev(X b (t)) stop increasing around t ≈ 4 s. In noncollapse processes all indentations stop growing before the elastic limit of the shell is exceeded and they disappear completely by the end of the granular discharge. So, the trajectories X a and X b eventually return close to their initial condition. Indeed, Dev(X a (t)) starts decreasing at t = 4.86 s as the trajectory starts returning towards X a (0). The restricted size of our sample does not allow us to observe the decrease in Dev(X b (t)) but its value starts fluctuating around its maximum as indentations stop growing. Conversely, Dev(X c (t)) exhibits a sustained increase, reflecting the ongoing proliferation and growth of indentations that leads to the collapse of the silo.

Another important property of the trajectory is how fast it evolves. In order to asses this, we define the magnitude of the displacement on the trajectory X at time t and time scale τ by

∆ τ (X(t)) = d W 2 (X(t), X(t + τ )). Figures ??(b)-(d)
show that for a short time scale, τ = 0.03 s, the values of ∆ τ (X a (t)), ∆ τ (X b (t)), and ∆ τ (X c (t)) are very similar. Moreover, they are comparable with fluctuating deformations of the shell induced by collisions during the granular flow. To estimate the magnitude of these fluctuations, we consider the trajectory X f ⊂ Per of a discharge with a low initial filling height, L f o , for which only deformations induced by collisions are present (see Sec. ??). The maximum value of ∆ 0.03 (X f (t)) is 0.1 and we interpret it as the magnitude of deviations caused by flow induced fluctuations along the trajectory of any deformation process. Once larger persistent indentations appear on the shell they are likely to favor collisions and thus increase the level of fluctuations. However, interactions between persistent and fluctuating deformations are expected to be complex [? ] and their analysis is out of the scope of the present paper.

For time scales longer than τ ≈ 0.40 s, ∆ τ (X) exhibits a more interesting behavior. Remarkably, the noncollapse trajectories X a and X b tend to evolve slower than X c , as corroborated by ∆ 0.43 (X(t)) shown in Figs. ??(b)-(d). After t ≈ 4 s, the values of ∆ 0.43 (X a,b (t)) and ∆ 0.03 (X a,b (t)) are comparable, indicating that the evolution of the systems is dominated by fluctuations [? ]. This does not come as a surprise considering the behavior of Dev(X a,b (t)) for t > 4 s.

The curve ∆ 0.43 (X c (t)), depicted in Fig. ??(d), indicates that X c evolves fastest soon after the onset of the discharge and slows down around t ≈ 1 s. This slowdown causes a drop in the slope of Dev(X c (t)). The slope of Dev(X c (t)) changes once again around t ≈ 4 s but with no corresponding change in ∆ 0.43 (X c (t)). This implies that, instead of slowing down, at this time the trajectory X c changes the direction in which it evolves. Visual inspection of the pattern shows that this coincides with a significant change in the dynamics of the deformation. Namely, the growth of individual indentations almost stops and the subsequent changes of the pattern are predominantly due to the merging of neighboring indentations. Merging constitutes an important growth mechanism for indentations, which individually do not overcome a characteristic size determined by the properties of the shell [? ]. In this way, very large indentations that cause the collapse of a silo can develop.

Figures ??(c)-(d) show that for t < 4 s the curves ∆ 0.43 (X a,b (t)) exhibit a sequence of pronounced maxima and minima. Thus, at the time scale τ = 0.43 s, both X a and X b alternate between regimes of slow and fast dynamics. In the rest of this section we will study this phenomenon in detail. To show that it is not an artifact caused by a particular choice of timescale, we introduce a distance matrix that provides a comprehensive description of the dynamics at all available time scales. We will also show how to detect the transitions between slow and fast regimes using the distance matrix. The results obtained will prove fundamental for our analysis of the origin of these transitions described in Sec .??.

Let X(t) in Per be a trajectory sampled at times t i , for 0 ≤ i ≤ n -1, then the distance matrix of X is an n × n matrix defined by

D(i, j) = d W 2 (X(t i ), X(t j )).
If there is no possibility of confusion, we denote D(i, j) by D(t i , t j ) to indicate the time corresponding to the samples i, j. In this notation D(t i , t j ) is the d W 2 distance between points X(t i ) and X(t j ). It quantifies the magnitude of the displacement on the trajectory between times t i and t j . So, ∆ τ (X(t i )) = D(t i , t i + τ ) and Dev(X(t i )) = D(0, t i ).

Figures ??(a)-(c) show distance matrices D a , D b and D c for the representative trajectories X a , X b and X c , respectively. These matrices exhibit different patterns, revealing interesting features of the dynamics of the underlying trajectories. We stress that similar patterns are observed in all matrices corresponding to the same class. Now we analyze the dynamics of the trajectory X a by using its distance matrix. In Fig. ??(a) one distinguishes a sequence of distinct dark blocks along the diagonal of D a . Inside of each block the evolution of X a is slow. To be more precise, if the times t i and t j fall in a single block, then the distance D a (t i , t j ) is of the order of the estimated value of fluctuations. Note that the separation between distinct dark blocks is rather sharp. This means that during a short time interval around the boundary between two successive blocks the dynamics suddenly becomes fast.

To detect transition points along a trajectory X, represented by a distance matrix D, we consider the following ratio inspired by Fisher's discriminant [? ],

ρ w (D, t i ) = w n,m=0 D(t i-n , t i+m ) w n,m=0 (D(t i-n , t i-m ) + D(t i+n , t i+m ))
, defined on a sliding window with length 2w. If the time interval between two consecutive fast transitions of X is longer than w∆t, and the duration of a single fast transition is short with respect to w∆t, then ρ w (D, t) attains prominent local maxima around the times at which these transitions take place. Thus, these time scales determine the choice of w for an optimal detection of the fast transitions. One must also consider that if w is too small fluctuations might compromise the results. Accordingly, we choose w ∈ [10,25], which corresponds to time scales between 0.3 s and 0.8 s. We find that the results are almost identical within this interval, thus we report the results only for w = 10.

For the FFC trajectory X a , Fig. ??(d) shows sharp peaks of ρ 10 (D a , t) with well defined maxima indicating abrupt transitions. The positions of the maxima are used to obtain the boundaries of the blocks, drawn in Fig. ??(a), separating time intervals during which the evolution of the orbit is slow. Also CTC trajectories exhibit such prominent transitions, as illustrated in Figs. ??(b) and ??(e) for X b . However, CTC trajectories exhibit a faster evolution than FFC trajectories at the onset of the discharge (see Figs. ??(a)-(c)). As a result, the peaks exhibited by ρ 10 (D b , t) are wider and more rugged, as shorter time scales are required to differentiate some transitions.

The curve ρ 10 (D c , t) for the C trajectory X c , shown in Fig. ??(f), is different from the curves for noncollapse processes. Its values tend to be larger and its peaks are less prominent, with a single exception soon after the onset of the discharge. The initial peak of ρ 10 (D c , t) is caused by the rapid slowdown of the system soon after the onset of the discharge. Therefore, at time scales larger than 0.3 s there is no evidence of an alternation between slow and fast dynamics in C processes. This is in accord with the behavior of ∆ 0.43 (X c (t)) shown in Fig. ??(d). Nevertheless, slow-fast dynamics may be present at shorter time scales. Indeed, ρ 5 (D c , t) exhibits a succession of peaks in the first 0.7 s of the discharge, when the C process evolves fastest. However, due to the temporal resolution of our sampling, it is hard to guarantee that these are not just artifacts caused by fluctuations of the system.

We will consider the origin of the fast-slow dynamics exhibited by the deformation processes in Sec. ??. Before that, we will discuss the reversibility of noncollapse trajectories and develop necessary machinery for a quantitative comparison of different trajectories.

C. Reversibility of noncollapse trajectories

At the end of a noncollapse process all visible deformation of the shell disappears. Thus the trajectories of noncollapse processes eventually return close to their initial condition. This is corroborated by the fact that distances between the initial conditions of all trajectories are smaller than the size of fluctuations. Since, hardware limitations restrict the duration of our sampling, we cannot observe the whole trajectory. However, we will demonstrate in this section how to estimate the time when trajectories of noncollapse processes start returning. This will prove useful for comparing different trajectories in Secs. ?? and ??.

We start by considering the FFC trajectory X a . Small values of D a in the rectangular region marked by a dotted line in Fig. ??(a) indicate that the states X a (t) for t ∈ [6.00, 8.30] s are close to states visited previously, for t ∈ [3.63, 4.23] s. Thus, the trajectory is retracing itself towards its initial condition. This is not surprising as it reflects the reversible nature of the deformation. Nevertheless, it is remarkable that we can directly observe it by considering only a small portion of the deformation pattern.

The last two transitions of X a in Figs. ??(a) and (d) occur while the trajectory is returning to its initial con-dition. Accordingly, we call them backward transitions. All other observed transitions occur as the deviation from the initial condition increases and we call them forward transitions. Note that X a evolves slower as it returns. Also, note that the backward transitions are less abrupt than forward transitions and correspond to relatively small peaks of ρ 10 (D a , t), see Fig. ??(d). Clearly, the trajectory X a starts returning to X a (0) at some time t ∈ [4. 27, 5.97] s, between the last forward and the first backward transitions. This is consistent with the fact that Dev(X a (t)) reaches its maximum at t = 4.84 s. In general, fluctuations make it difficult to determine the precise time when a trajectory starts returning. However, the last forward and the first backward transitions can be identified with certainty for all FFC processes. The fact that the trajectories move very little between these transitions implies that the turning point of a trajectory X is contained in a small ball around X(t r ), where t r is the time of the last forward transition. In the following, we refer to t r as the time of onset of return.

In the case of CTC processes our sample time is too short to capture the first backward transition. For the trajectory X b , Figs. ??(b) and ??(e) show that after the transition at t ≈ 4 s the system merely fluctuates during the rest of the sampled time. Given the total duration of the discharge, and the time the system spends fluctuating after the last observed transition, we expect that the trajectory does not evolve much further before it turns around. Thus, in analogy with the FFC case, we refer to the time of the last detected transition of a CTC trajectory as its time of onset of return, t r .

We expect the deformation to start receding after the effective axial stress on the shell starts decreasing. In our experiments the granular discharge is triggered right after the filling of the silo. Under these conditions, one expects that local frictional forces exerted by the grains will not satisfy the maximal mobilization condition [? ? ? ? ]. Thus, the effective axial stress will increase at the beginning of the discharge as a result of a spontaneous mobilization process [? ? ]. Once the maximum possible stress component is transmitted to the shell, Janssen's saturating profile [? ] is realized, and the axial stress starts decreasing. One can thus interpret t r as a rough estimate for the duration of the mobilization process.

We remark that since t r is defined as the time of the last forward transition it represents only a lower bound on the time when the trajectory starts returning. A more precise determination of the return time would require the observation of the whole pattern. It is therefore not surprising that we find a large dispersion in the values of t r . For FFC and CTC trajectories t r ∈ [2, 4] s and there is no evidence of a dependence on their class or initial filling height. The values of t r are larger but of the same order of magnitude as the duration of the mobilization regime measured by [? ] for a system with comparable dimensions, except for the shell thickness, and different boundary conditions. On the other hand, in [? ] it has been observed that the duration of the mobilization process can be substantially longer than the values of t r , as it is associated to the details of the granular flow pattern inside the silo.

D. Comparing the trajectories

So far we studied the dynamics of individual trajectories. In this section, we consider distances between all the sampled states to investigate relations between the trajectories for different discharges. We show that the trajectories of all FFC processes lie in a narrow tube in the space of persistence diagrams. Therefore, the evolution of the deformation patterns for all FFC processes is topologically similar, although the geometry of the patterns can vary significantly. Trajectories of CTC and C processes also start inside the narrow tube containing the FFC processes but exit it a short time after the onset of the discharge. Nevertheless, during this short time they traverse a large portion of the tube.

To construct the tube encompassing the FFC trajectories we use the point to set d W 2 distance in Per, defined as follows. Let S ⊂ Per be a set and x a point in Per, then the distance between x and S is given by

d(x, S) = min {d W 2 (x, y) : y ∈ S} .
For every FFC trajectory X we define a tubular neighborhood with radius q by O q (X) = {x ∈ Per : d(x, X([0, t r ])) < q}.

Note that O q (X) is a neighborhood of the forward moving part of X, i.e., X([0, t r ]) = {X(t) : t ∈ [0, t r ]}. Let q X be the smallest radius such that O q X (X) contains all FFC trajectories. We choose the FFC trajectory X R for which q X R = min{q X : X ∈ FFC }. This gives q X R = 0.17, which is of the same order as the estimated lower bound on the fluctuations in the system. To simplify the notation we denote O q X R (X R ) by O. We estimate the length of the tube O using the arc length of the forward moving part of the trajectory X R . We could define the arc length as the sum of distances between the consecutive sample points. However, to mitigate the artificial elongation of the trajectory caused by fluctuations, we consider displacements over a longer time scale τ = 0.43 s (corresponding to 13 time steps) and parametrize the trajectory by

(X(t i )) = 1 13 i j=0 ∆ 0.43 (X(t j )).
Using this definition, the arc length of the forward moving part of the trajectory X R is (X R (t r )) = 1.37 > 8q X R . We stress that all observed FFC trajectories are contained in the relatively narrow tubular neighborhood O centered around X R . Thus, all FFC trajectories closely follow X R and it is reasonable to use X R as a reference trajectory. Figure ?? shows the deviation of all trajectories from the forward moving part of X R . Note that all the trajectories start close to X R . By construction, every FFC trajectory X stays in O and d(X(t i ), X R [0, t r ]) < 0.17 for all t i . However, C trajectories leave the set O in less than 1.5 s and they keep moving away from X R until the collapse of the silo. CTC trajectories tend to stay in O longer, some of them exit only after 4 s, and their distance to X R eventually stops increasing as the pattern stabilizes.

We have shown that C and CTC trajectories start within the tube O but they exit it very fast. Finally, we identify what fraction of O is traversed by C and CTC trajectories before they exit. To do this, we define the arc length coordinate for any point in the tube as follows. Let x ∈ O, then (x) = (X R (t i (x))), where t i (x) = argmin tj ∈[0,tr] {d w 2 (x, X R (t j ))}. Suppose that a trajectory X exits the tube O at exit time t e , i.e. X(t) ∈ O for 0 ≤ t ≤ t e and X(t) ∈ O for t > t e . Then (X(t e )) measures the part of O traversed by X before exiting.

For all C and CTC trajectories (X(t e )) > 0.99. This means that they follow closely a large part of X R , so that initially they are similar to the FFC trajectories. On the other hand, the fact that the C processes exit O in such a short time indicates that, already at an early stage of the discharge, their dynamics is very different from FFC processes. Note that t e for all C processes is smaller than the minimum value of t r , our estimated lower bound for the duration of the mobilization process. This supports the idea, proposed in [? ? ], that the collapse of the silo is triggered by partially mobilized friction. Moreover, the evolution of all C processes is fast compared to FFC processes. In a quarter of the time required by X R they traverse more than 70% of its length. This suggests that comparing how fast the trajectories evolve at the onset of the discharge may allow for an early classification of the deformation processes. We explore this idea in Sec. ??.

E. Early collapse detection and the initial condition

Early prediction of silo collapse is an important problem in many industrial applications. As a toy example, we demonstrate a method that uses TDA for predicting the collapse in the case of our paper silos. The silo is known to collapse when its initial filling height exceeds the critical value L c [? ? ]. Thus one can expect that knowledge of L o and L c is enough to predict the outcome of the granular discharge. Unfortunately, as visible in Fig. ??(a), the error in L c is quite large and predicting the class of a deformation process based solely on L o leads to misclassifications.

In Sec. ?? we showed that C processes evolve faster than noncollapse processes at the onset of the discharge. Thus, we propose a classifier based on the value of the initial displacement rate ∆ τ (X(0)). For small values of τ the classifier does not yield satisfactory results because ∆ τ (X(0)) is dominated by noise. However, for time scales τ > 0.4 sec, long enough to overcome noise, it performs consistently well. In particular, the value of ∆ 0.43 (X(0)) allows us to distinguish collapse processes better than the value of L o .

As shown in Fig. ??(a), for all C processes ∆ 0.43 (X(0)) > 0.28, while for noncollapse processes ∆ 0.43 (X(0)) < 0.28. It does not come as a surprise that for some CTC processes ∆ 0.43 (X(0)) is very close to the threshold between these two classes, since according to our qualitative criteria CTC processes come very close to causing the collapse of the silo. We stress that by using the value of ∆ 0.43 (X(0)) we distinguish collapse and noncollapse processes already 0.43 s after the onset of the discharge. At this time all trajectories in our sample are inside the neighborhood O and thus indistinguishable from an FFC trajectory through instantaneous observations.

Figure ??(a) suggests that there exists a relationship between the initial rate of evolution, ∆ 0.43 (X(0)), and L o . Although the dispersion of the results is fairly large, it is reduced when considering each silo separately. This suggests that some intrinsic properties of individual shells affect the dynamics under the discharge. It is possible that microscopic defects of the paper sheet used to build a silo and small asymmetries introduced during its construction are not detected by visual inspection. Such defects can make a shell more prone to develop initial deformations during filling. These are known to have an important effect on the resistance of the shell to axial compression [? ? ] and are thus likely to alter the dynamics of the deformation during the granular discharge. In particular, we expect ∆ 0.43 (X(0)) to be largely determined by the initial deformations, which must be reflected in the initial condition X(0).

To assess the effect of the initial condition, for each trajectory X we consider the value of (X(0)). Figure ??(b) shows that ∆ 0.43 (X(0)) tends to increase with (X(0)). So, a larger initial deformation leads to a faster evolution of the deformation pattern at the beginning of the discharge. We expect that the large spread of the points in Fig. ??(b) is caused by the fact that we consider only a small portion of the deformation pattern. We also note that for the FFC processes (X(0)) is smaller than for C and CTC processes. This suggests that the dynamics of the deformation pattern can be predicted from the initial condition. Thus, it is likely that by considering the initial deformation pattern on the whole surface of the silo one could distinguish collapse and noncollapse processes even before the discharge begins.

F. Nature of the fast transitions

In this section we discuss the mechanism causing the fast-slow dynamics described in Sec. ??. At first glance, the observed abrupt transitions might look surprising because the effective axial stress on the shell is expected to evolve rather smoothly [? ? ]. Naturally, fluctuations induced by the granular flow may affect the effective stress on the shell [? ]. Moreover, both simulations [? ? ] and direct measurements on large scale silos [? ? ] exhibit very large local stress fluctuations, particularly at the onset of the granular discharge. Thus, one possible explanation is that fast transitions in the evolution of the deformation are triggered by fluctuations induced by the flow. However, in what follows we provide evidence that stress fluctuations are not the main source of these abrupt transitions. We observe that transitions do not occur at random places along the tube O, as one would expect if they were predominately caused by fluctuations. Thus, we argue that the slow-fast dynamics is intrinsic to the response of the shell during the granular discharge.

We concentrate only on the FFC trajectories because they are confined in the narrow tube O around the trajectory X R . Thus, the positions of their transitions can be meaningfully compared using the arc length coordinate of points in O defined in Sec. ??. We use the methods described in Sec. ?? to detect the transitions of FFC trajectories. As already mentioned, the last forward transitions can be identified with certainty for all FFC processes. Figure ?? shows their positions along the tube O. It is clear that they typically fall close to some transition of the trajectory X R . This observation holds for 85% of the transitions detected by our methods. In the remaining cases, visual inspection of the distance matrices shows the detection of the transition time is not optimal, either due to fluctuations or because our sampling time is too long with respect to the time between two successive transitions. Moreover, for 50% of the FFC trajectories all consecutive transitions fall close to consecutive transitions of X R .

The moderate amount of misaligned transitions suggests that the slow-fast dynamics is intrinsic to the behavior of the shell. This behavior is indeed reminiscent of the behavior of an empty cylindrical shell subjected to axial compression. It has been experimentally observed [? ] that, as the axial stress is increased, the empty shell suffers a sequence of successive transitions towards metastable states with increasing curvature. These metastable states are localized with respect to the length of the cylinder and periodic around its circumference. As suggested by [? ], the circumferential wave number can be predicted from the cylinder length. The lack of the characteristic axial symmetry in our experiments can be explained by a stabilizing effect of the flowing granular material [? ]. In fact, one observes similar patterns when a close fitting rigid mandrel is fitted inside a thin cylindrical shell subject to axial compression [? ? ], in this case the mandrel constrains the depth of any developing indentations. However, the precise nature of any stabilization mechanism induced by the granular flow remains obscure and out of the scope of the present paper. Now we return to Fig. ?? and discuss the evident sep-aration of the last forward transitions into two distinct groups. For trajectories with ∆ 0.43 (X(0)) < 0.10 the value of (X(t r )) is generally smaller than 0.70 and their last forward transition corresponds to one of the first four transitions of X R . On the contrary, for all trajectories with ∆ 0.43 (X(0)) > 0.10 the value of (X(t r )) is larger than 0.99 and their last forward transition reaches at least the fifth transition of X R . This suggests the presence of a bifurcation along the control parameter ∆ 0.43 (X(0)), which describes how fast the deformation evolves at the onset of the discharge. Thus, we hypothesize that the transition exhibited by X R at (X R (t r )) = 1.01 can only be triggered if the initial evolution of a deformation process is fast enough. Namely, if ∆ 0.43 (X(0)) is larger than the critical value 0.10. Moreover, we propose that there is a second bifurcation along the parameter ∆ 0.43 (X(0)) around the value 0.28, that accounts for the separation between collapse and noncollapse trajectories observed in Fig. ??.

We remark that we do not observe any bifurcation along L o and (X(0)). We attribute this to the fact that they fail to capture some important aspect of the initial condition. As described in Sec. ??, the effect of any initial defects of the shell cannot be accounted for by L o . We also recall that we only observe a fraction of the silo, so (X(0)) does not provide complete information about the initial condition. The parameter ∆ 0.43 (X(0)) is also obtained only from a fraction of the shell, but the observed initial evolution of the deformation seems to be affected by the global state of the shell. Thus, it yields a more comprehensive description of the initial state of the system.

IV. CONCLUSION

In this paper we have demonstrated that TDA constitutes an effective framework to investigate the complex deformation patterns sustained by thin walled silos during gravity driven granular discharges. Using TDA, a deformation process is associated with a trajectory in the space of persistence diagrams. This allows for a meaningful description of the time evolution of the deformation and a quantitative comparison of different processes.

We found that a short time after the onset of the granular discharge, significant dynamical differences arise between those deformation processes that are reversible and those that lead to the collapse of the silo. By this time, the trajectories representing all the deformation processes are very close, but the deformation develops much faster in processes that lead to the collapse of the silo. Based on these observations, we proposed an early classification criterion that distinguishes processes that lead to a collapse while the deformation of the shell is still incipient. This is remarkable considering the low temporal and spatial resolution of our sample, and the fact that we observe only a small portion of the pattern. The efficacy of this toy-model prediction method for the collapse demonstrates the strong potential of our methods for designing early-warning signals. Moreover, recent progress in computational homology makes it possible to detect the differences in the dynamics of the processes in real time, even for considerably larger data sets than considered in this paper.

We also showed that reversible deformation processes are characterized by slow-fast dynamics. Our results suggest that this is associated to the intrinsic nonlinear behavior of the shell, evocative of the sequence of post buckling deformation patterns observed on empty cylindrical shells subject to an increasing axial compression. However, additional experimental observations are needed to study this phenomenon. In particular, we found that in some cases our sampling rate was too slow to detect accurately the transitions between intervals of slow evolution. Specially for deformation processes that lead to the collapse of the silo, for which our sampling rate does not allow to conclusively decide if slow-fast dynamics is present.

FIG. 1 .

 1 FIG. 1. Snapshots of deformation patterns observed in our experiments. Using the qualitative criteria of Sec. ??, the corresponding deformation processes are classified as (a) farfrom-collapse, (b) close-to-collapse, and (c) collapse processes. (a) and (b) show the shell when the deformation stops growing, right before it starts returning to its original undeformed state. (c) shows a shell after catastrophic failure. The upper section of the shell has collapsed against a protective glass cylinder surrounding the silo.

FIG. 2 .

 2 FIG. 2. (Color online) (a) Schematic cross-section of the experimental setup indicating its dimensions. L indicates the height of the granular column. A typical image used in our analysis is shown, marking the area of the shell it captures. Due to oblique illumination from above, the upper (lower) parts of diamond shaped indentations produce dark (bright) features in the image. (b) Field ω obtained from processing the image of the deformation pattern shown in (a). (c) Schematic representation of ω(θ, z), for fixed θ, of the crosssection at the center of a single diamond shaped indentation in the Yoshimura pattern. Referential directions corresponding to cylindrical coordinates (r, θ, z) on the shell are indicated. Here θ = θ(x, y) and z = z(x, y), with (x, y) Cartesian coordinates on the image, as indicated in (d). (d) Schematic representation of ω(x, y), corresponding to the cross-section of the indentation shown in (c).

  FIG. 3. (Color online) (a) Deviation from the initial condition along typical trajectories representing FFC, CTC and C deformation processes X a (black line), X b (blue line), and X c (pink line), respectively. (b)-(d) Magnitude of the displacement for X a , X b and X c at two different time scales, τ = 0.03 s (gray line) and τ = 0.43 s (lines colored as in (a)).

1 FIG. 4 .

 14 FIG. 4. (Color online) Distance matrices (a) D a , (b) D b and (c) D cfor the trajectories X a (FFC), X b (CTC) and X c (C), respectively. Ratio (d) ρ10(D a , t), (e) ρ10(D b , t), and (f) ρ10(D c , t). In (d), vertical lines mark the maxima of ρ10(D a , t), corresponding to fast transitions in the evolution of X a . The green line marks the last forward transition of X a , i.e., its onset of return time, tr. Dark and light gray lines correspond to the remaining forward and backward transitions of X a . In (a) boundaries of diagonal blocks exhibiting slow dynamics are identified as transition times marked in (d), with corresponding line coloring. The dotted gray rectangular line marks an off diagonal block within D a with small values, which signals the onset of the return of X a towards its initial condition.

FIG. 5 .

 5 FIG.5. Distance d W 2 (X(t), X R ) between X(t) and the reference trajectory X R for all FFC (black lines), CTC (blue lines) and C (pink lines) trajectories in our sample. The black dashed line marks the radius , q X R of O.

FIG. 6 .

 6 FIG.6. For each trajectory X in our sample, ∆0.43(X(0)) versus: (a) Initial filling height of the corresponding granular discharge. The continuous (dotted) vertical gray line marks the critical filling height, Lc (and its error bounds). The horizontal line marks a classification threshold for C trajectories. (b) (X(0)). All trajectories X resulting from successive discharges of a single silo are represented with the same symbol, they are colored according to the class of X: FFC (black), CTC (blue) or C (pink).

FIG. 7 .

 7 FIG.7. Arc length coordinate of the points X(tr) for each FFC trajectory X versus ∆0.43(X(0)). Vertical lines mark the coordinates of transition points of the trajectory X R .

V. ACKNOWLEDGEMENTS

The authors would like to thank Konstantin Mischaikow for many fruitful discussions and the use of his cluster computer Conley3. They also express their gratitude to Lou Kondic for very useful comments on the manuscript. CC would like to acknowledge the late Gustavo Gutiérrez, who contributed with many valuable ideas and the laboratory where experiments were performed. Efficient algorithms implemented in the GUDHI library were used to compute the persistence diagrams. MK was supported by ERC project GUDHI (Geometric Understanding in Higher Dimensions).