Alexander Koller
email: koller@mmci.uni-saarland.de

Jörg Hoffmann
email: joerg.hoffmann@inria.fr

Waking Up a Sleeping Rabbit: On Natural-Language Sentence Generation with FF

We present a planning domain that encodes the problem of generating natural language sentences. This domain has a number of features that provoke fairly unusual behavior in planners. In particular, hitherto no existing automated planner was sufficiently effective to be of practical value in this application. We analyze in detail the reasons for ineffectiveness in FF, resulting in a few minor implementation fixes in FF's preprocessor, and in a basic reconfiguration of its search options. The performance of the modified FF is up to several orders of magnitude better than that of the original FF, and for the first time makes automated planners a practical possibility for this application. Beside thus highlighting the importance of preprocessing and automated configuration techniques, we show that the domain still poses several interesting challenges to the development of search heuristics.

Introduction

Natural language generation (NLG) -the problem of computing a sentence or text that communicates a given piece of information -has a long-standing connection to automated planning [START_REF] Perrault | A plan-based analysis of indirect speech acts[END_REF]. In the past few years, the idea has been picked up again in the NLG community, both to use planning as an approach for modeling language problems (e.g., [START_REF] Steedman | Planning dialog actions[END_REF]) and to tackle the search problems in NLG by using techniques from planning [START_REF] Koller | Sentence generation as planning[END_REF].

Existing planners are relatively good at solving modest instances of the planning problems that arise in this context. However, [START_REF] Koller | Experiences with planning for natural language generation[END_REF] recently showed that as the planning problem instances from [START_REF] Koller | Sentence generation as planning[END_REF] scale up, off-the-shelf planners become slow -indeed, too slow to be practically useful in NLG applications. They showed this for SGPlan and several variants of FF [START_REF] Hoffmann | The FF planning system: Fast plan generation through heuristic search[END_REF]; we show it herein also for LAMA [START_REF] Richter | Landmarks revisited[END_REF]. Koller and Petrick noted that there are major sources of inefficiency both in the search itself and in the commonly used preprocesses; they identified several features of the domain (e.g. comparatively large predicate and operator arities) that diverge from most IPC benchmarks. In this paper, we revisit this situation, and analyze the sources of inefficiency in detail for FF. Based on the analysis, we modify FF and obtain dramatically better performance that, for the first time, makes automated planners a realistic possibility for NLG applications. To achieve this quantum leap, all that is needed are a few minor implementation fixes in FF's preprocessor, and basic reconfiguration of its search options. While these changes certainly are no contribution to planning in their own right, we consider it a valuable -and slightly unsettling -lesson learned that such banalities can make the difference between success and failure. They can hold up progress in planning applications for years, unless users of a planner have direct access to its developer. We interpret this as a strong call for more research on automated configuration techniques, e.g. along the lines of [START_REF] Hutter | Automatic algorithm configuration based on local search[END_REF]. The detrimental effect of preprocessing details also sheds a light on the importance of this often neglected part of the planning process.

Apart from these insights, our domain remains an interesting challenge for the development of better heuristics. There still are realistic instances of the NLG planning domain on which our modified FF is not practical. Overcoming this weakness would entail overcoming fundamental weaknesses of FF's goal ordering techniques and the relaxed plan heuristic -variants of which are used in many planners today.

Sentence generation as planning

We consider a domain pertaining to the generation of natural language sentences based on tree-adjoining grammar (TAG; [START_REF] Joshi | Tree-Adjoining Grammars[END_REF]). This is a standard problem in computational linguistics, which we sketch by example; see [START_REF] Stone | Microplanning with communicative intentions: The SPUD system[END_REF][START_REF] Koller | Sentence generation as planning[END_REF] Figure 2: Actions for generating "The white rabbit sleeps." of Fig. 1 shows three trees for "the rabbit", "sleeps", and "white". Trees can be combined using the operations of substitution (purple dashed arrow in the figure) and adjunction (blue dotted arrows) to form larger trees. The end result of a grammatically correct derivation is a tree all of whose leaves are labeled with words, as on the right of Fig. 1. We can read off a sentence from such a tree from left to right.

To use TAG grammars for generation [START_REF] Stone | Microplanning with communicative intentions: The SPUD system[END_REF], we assume a set of ground atoms expressing the information we want the sentence to convey, such as {sleep(e, r 1 )}, and a knowledge base that contains all ground atoms we know to be true; say, {sleep(e, r 1 ), rabbit(r 1 ), rabbit(r 2 ), white(r 1 ), brown(r 2 )}. We then add variables ranging over constants from the knowledge base to the nodes of each elementary tree, and equip every elementary tree with a set of atoms over these variables to encode the meaning this tree expresses (printed in italics in the figure). Fig. 1 already shows specific instances of the elementary trees from the grammar, in which constants have been substituted for these variables. The derivation in the figure conveys the intended meaning -in particular, that r 1 sleeps. Crucially, it also describes uniquely who does the sleeping: The sentence "the rabbit sleeps" would not have done this, as "the rabbit" could be understood either as r 1 or r 2 . We say that r 2 is a distractor for the subject of the sentence, which is removed by adjoining the tree for "white". In short, the sentence generation problem consists in computing a grammatically correct TAG derivation from instances of the elementary trees in the grammar that conveys the intended meaning and describes all referents uniquely.

Sentence generation as a planning problem

The TAG sentence generation problem can be encoded as a planning problem [START_REF] Koller | Sentence generation as planning[END_REF]. The key idea is that each operator encodes the addition of an elementary tree to the derivation; the syntactic and semantic requirements and effects of doing this are captured in the operator's precondition and effect.

More precisely, each operator has a parameter u representing the syntax node into which the elementary tree is substituted or adjoined, and a parameter u i for each substitution node. There are also parameters x 0 , . . . , x k for the variables in the semantic representation of the elementary tree;

x 0 is the variable that occurs at the root of the tree. The planning state encodes a list of possible node names: current(u) expresses that u is the next node name that should be picked when a new substitution node is created, and next(u, v) says that the next node after u is v. These atoms are used to select names for the substitution nodes that are introduced by adding an elementary tree; the parameter u n represents the node that will be current after the addition.

The atom subst(A, u) expresses that we need to still substitute something into the substitution node u with label A; the initial state contains an atom subst(S, root) where S stands for "sentence" and root is an arbitrary name for the root of the derivation. The mapping from nodes to semantic constants is maintained using the ref atoms; the initial state for generating a sentence about the event e contains an atom ref(root, e). Finally, we keep track of the uniqueness of referring expressions in the distractor atoms: distractor(u, a) expresses that the expression at syntax node u could still be misunderstood as a.

The example generation problem above translates into a planning problem whose domain is shown in Fig. 2. The initial state of the problem encodes the generation problem's knowledge base; it contains atoms rabbit(r 1 ), sleep(e, r 1 ), etc. The goal requires syntactic completeness as ∀u∀A¬subst(A, u) and unique reference as ∀u∀x¬distractor(u, x); it also specifies the semantic representation we want to express in an atom expressed(sleep, e, r 1 ).

A plan for the example problem is sleeps(root, n 1 , n 2 , e, r 1 ); rabbit(n 1 , r 1 ); white(n 1 , r 1 ). This plan can be automatically decoded into the derivation shown in Fig. 1. The first two steps of the plan alone would not be a correct plan because they leave an atom distractor(n 1 , r 2 ) in the state, which contradicts the goal that all distractors have been eliminated.

FF Pitfalls and Fixes

The basic observation of FF's inefficiency in the NLG domain was made before [START_REF] Koller | Experiences with planning for natural language generation[END_REF]; here we examine its causes in detail. For some of the observed pitfalls, we are able to provide simple fixes. Others pose serious challenges to planning research.

Preprocessing

While preprocessing the ADL formulas, FF has a subroutine that checks the grounded formulas for tautologies. This involves a loop checking all pairs of sub-formulas within every conjunction and disjunction, testing whether they are identical/identical modulo their sign. For example, if we find φ and ¬φ within a disjunction, then the disjunction is replaced by TRUE. We haven't tested whether this subroutine has a notable effect in other benchmarks. In our domain, the only effect it has is to take a huge amount of runtime on the long goal conjunctions stating that we do not want any open substitution nodes or distractors. The fix is to switch this subroutine off. This trivial operation alone can render previously infeasible instances (preprocessing time > 15 minutes) harmless (instance solved in 6 seconds total time). FF's preprocessor contains a very small imprudence that apparently didn't surface in other domains, but that does in ours. FF distinguishes operators into "easy" (DNF precondition) and "hard" (non-DNF precondition), which go through different instantiation procedures. The procedure for hard operators uses a data structure encompassing an integer for every possible instantiation of all predicates. Since this does not take into account reachability, it is prohibitively huge in our domain. Now, we actually have no hard operators. But FF builds the data structure anyway, which exhausts memory and causes a segmentation fault on serious instances.

Search Techniques

FF's preprocessing pitfalls are certainly baffling, but easily fixed. We now turn our attention to some of FF's search techniques, variations of which are used in many planners, and which are very much not easily fixed. Our solution for the moment are simple configuration changes.

First, we consider [START_REF] Koehler | On reasonable and forced goal orderings and their use in an agenda-driven planning algorithm[END_REF] "reasonable goal orderings", variants and extensions of which are in prominent use today, e.g. in LAMA [START_REF] Richter | Landmarks revisited[END_REF]. FF as used in the 2000 competition approximates reasonable orders as a preprocess and then partitions the goal set into a goal agenda, a sequence of goal subsets. Each subset is posed to the planner in isolation. The final plan is the concatenation of the plans for all entries.

The goal agenda is completely off-target in our domain. The most striking point concerns the interaction between the goal to create a sentence (¬subst(S, root) in the rabbit example), and the goal to express a certain meaning (expressed(sleep, e, r 1 )). The goal agenda orders the former before the latter because it detects that, once a meaning is expressed, the sentence creation has started and it is not possible to start anew anymore. Unless the inverse also holds,1 subst(S, root) is alone in the first goal agenda entry. Consequently, the planner decides to generate any sentence, with some arbitrary meaning. If it happens to be the wrong sentence, the planner is stuck in a dead-end. It is not clear to us, at this point, how goal ordering techniques could be modified so that they return more sensible results.

In our improved FF, the goal agenda is simply switched off. Having thus disposed of goal orderings, let us turn our attention to the relaxed plan heuristic. This meets its Waterloo in the unusual life-cycle of facts encoding substitution nodes and distractors. These facts are required to be false in the goal, and most of them are initially false. However, any valid plan must make many of them true at some intermediate time point (e.g. distractors arise depending on the decisions taken by the planner). Hence, from the perspective of the relaxed plan heuristic, from one state to the next whole classes of goals suddenly appear "out of the blue". Consider again the rabbit example. In the initial state, the relaxed plan has length 1 because we only need to achieve ¬subst(S, root) and expressed(sleep, e, r 1 ), which is done by the single action "sleeps". However, that action has the effect ∀y.y = x 1 → distractor(u 1 , y), which introduces new distractors. Thus the relaxed plan for the second state contains two actions. In general (if there are n rabbits, say), the relaxed plan may become arbitrarily long. For FF's particular search algorithm, enforced hill-climbing (EHC), this means that FF turns into a blind breadth-first search: EHC starts at h value 1 and thereafter needs to solve the whole problem in order to find a (goal) state with lower h value.

More generally, from this example we immediately see that, in the terms of [START_REF] Hoffmann | Where 'ignoring delete lists' works: Local search topology in planning benchmarks[END_REF], the exit distance from local minima under h + is unbounded. Also, we have unrecognized dead-ends, i.e., dead-end states that have a relaxed plan. An example is the initial state of a task that it unsolvable because not all distractors can be removed; its h value will still be 1. So the domain is in the most difficult class of [START_REF] Hoffmann | Where 'ignoring delete lists' works: Local search topology in planning benchmarks[END_REF] planning domain taxonomy. While it shares this property with several other challenging domains, turning EHC into a breadth-first search is unheard of.

Can we design heuristics that are better at dealing with the life-cycle of substitution nodes and distractors? It is easy to see that the problem of not noticing them in the initial state is present for most known heuristics [START_REF] Helmert | A planning heuristic based on causal graph analysis[END_REF][START_REF] Richter | Landmarks revisited[END_REF][START_REF] Karpas | Cost-optimal planning with landmarks[END_REF][START_REF] Helmert | Landmarks, critical paths and abstractions: What's the difference anyway?[END_REF][START_REF] Cai | Enhancing the context-enhanced additive heuristic with precedence constraints[END_REF]. The two notable exceptions are the h m family, and abstraction heuristics. We ran a few tests with a variant of the example problem involving 10 rabbits. The Graphplan estimate for the initial state is 4, i.e., h 2 is better than h + here (although still far from the real goal distance of 14). With merge&shrink [START_REF] Helmert | Flexible abstraction heuristics for optimal sequential planning[END_REF], we obtained the estimate 11 for the initial state, which is much better. Still the search space with that latter heuristic remained fairly large (5200 nodes vs. 13800 nodes for blind search), and it is not clear how well the technique will scale to larger instances and instances with a more complex grammar. Certainly, interesting challenges remain.

The only FF technique that appears to be clearly useful here is its action pruning technique. While the relaxed plan is often much too short, it appears to always contain at least one of the actions that actually are useful. Given the above observations, we modified FF to directly run best-first search with helpful actions (as opposed to original FF, which switches to best-first without helpful actions if EHC fails).

Experiments

To evaluate the effect of our changes on the efficiency of the planner on practical inputs, we ran a series of experiments using the XTAG Grammar (XTAG Research Group 2001), a large-scale tree-adjoining grammar of English, from which we selected all lexicon entries for seven words we used in the experiments. We set up inputs to generate sentences of the form "S 1 and S 2 and . . . and S n ", i.e. a conjunction of n sentences. Each sentence S i is of the form "X admires Y". X and Y are unique referring expressions, such as "the rich businessman". For each experiment, we have a parameter d that specifies how many adjectives are needed to distinguish the referent uniquely -in particular, d = 0 means we can say "the businessman", and d = 2 means we must say "the rich sad businessman". Consider first the results for d = 0, shown in Fig. 3 (a); the lines end where the planner exceeded its five-minute timeout. Note that the vertical axis is logarithmic. Clearly, "Old FF" (the original FF) is much too slow for practical use. The same is true -although to a less dramatic extent -for LAMA, which already takes 3 seconds to generate the relatively short sentence at n = 3 and times out for larger instances. "New preprocessor" is FF with only our fixes to FF's preprocess, "New FF" applies all fixes. The preprocessor fixes alone lead to a speed-up of at least 3 orders of magnitude, leading to practically useful runtime -e.g., the 24-word sentence at n = 4 is generated in under two seconds. By contrast, the search configuration does not make much of a difference. This changes in our next experiment.

Results for d = 2 are shown in Fig. 3 (b). Here the old FF preprocessor crashed with a segmentation fault even at n = 1. LAMA performs better, but times out beyond n > 2. Best-first search with helpful actions is clearly the best-performing search strategy here. It takes about 14 seconds for n = 3, but this sentence is already quite long at 29 words, and EHC takes 3.5 minutes on the same input. Best-first search without helpful actions timed out at n = 3, illustrating the importance of that pruning technique.

We remark that the relative behavior of the different search methods is not uniform in our domain, which exhibits a huge amount of structural variance due to the many different forms that the input grammar may take. We are far from an exhaustive exploration of all the possible parameter settings. For instance, with particular inputs (no distractors), EHC can be more effective than best-first search.

Conclusion

Generation of natural language sentences is interesting as a planning problem, both in its application value and in the challenges it poses. After the changes to FF described in this paper, for the first time an automated planner is a practical option for this application. The simplicity of the changes required is striking, and strongly suggests to look more closely at automatic planner configuration techniques.

The domain is far from "solved" in all its aspects. As we pointed out, several features of the domain are very hard to capture with existing search heuristics. In Fig. 3, this shows in the fact that no planner version scales beyond n > 4even in case (a) where no adjectives at all are required. It remains an interesting challenge to address this.

Figure 1 :

 1 Figure 1: Derivation of "The white rabbit sleeps."

Figure 3 :

 3 Figure 3: Runtimes for d = 0 (a) and d = 2 (b).

The inverse holds if, in the given grammar, the only way to communicate a meaning is in the main clause -i.e., there are no subordinate clauses like "the sleeping rabbit".

Acknowledgements. We thank Konstantina Garoufi for driving sentence generation with FF to its limits, Ron Petrick for fruitful discussions, Silvia Richter for making LAMA available, and the reviewers for their helpful comments.


Séverine Anne
email: severine.anne@cnrs-orleans.fr

Olivier Rozenbaum
email: olivier.rozenbaum@cnrs-orleans.fr

Pascal Andreazza
email: pascal.andreazza@univ-orleans.fr

Jean-Louis Rouet
email: jean-louis.rouet@univ-orleans.fr

Evidence of a bacterial carbonate coating on plaster samples subjected to the Calcite Bioconcept biomineralization technique

Keywords: Coating, bioremediation, carbonate production, plaster, grazing incidence, X-ray diffraction, bacterial surface treatment, weathering, SEM, microprobe mapping

Degradation of historical buildings is mainly due to the intrusion of water which is the main vector of pollutants. Different types of surface treatment have been proposed to avoid or limit this effect. One alternative to chemical treatments is the use of the carbonatogenesis property of some bacteria. This bacterial production has been evidenced on concrete and on limestone samples in an aqueous environment. However, the carbonate production was measured indirectly and the experimental protocol was far from real conditions of use. In this paper, with the same protocol as an industrial one, and using a surface selective investigation method, grazing incidence X-ray diffraction, we show the structural and morphological evolution of the carbonate coating produced on model plaster samples. This substrate was chosen in order to unambiguously detect the bacterial carbonate production.

Introduction

Weathering of buildings is a widespread problem encountered in most countries around the world, and concerns buildings made of concrete as well as historic buildings made of stone or brick. These porous materials are subjected to deterioration due to the action of external environmental agents (physical, chemical and biological) [START_REF] Camuffo | Physical weathering of stones[END_REF]2]. In all cases, water transfer within the whole volume of the porous media is the common point to weathering. To prevent this transfer, several surface treatments have been proposed in the past. The conventional ones involved water repellents that were also used for the chemical reinforcement of stone [START_REF] Amoroso | Stone decay and conservation : Atmospheric Pollution, Cleaning, Protection[END_REF]. In the short term, they give rather good results but can prove harmful in the long term due to their irreversible features. Indeed, these resins polymerize, resulting in a film on the pore surface, but in some cases, depending on the stone, the resin polymerizes inside the pores, completely filling up the porosity [START_REF] Perrier | Les roches ornementales[END_REF]. In the latter case, although water transfer 1 from outside toward the inside of the stones treated by water repellent is totally stopped, water (in its gaseous form) inside the pores can no longer escape and remains near the polymerized resin. And finally, water repellent induced desquamations [START_REF] Camaiti | Prodotti consolidanti impiegati nelle operazioni di restauro[END_REF][START_REF] Amoroso | Stone decay and conservation : Atmospheric Pollution, Cleaning, Protection[END_REF]. Moreover, it is clear that while it is necessary to prevent the intrusion of water into the stone, it is also important to maintain a gaseous exchange between the stone and its environment [START_REF] Černý | Methods for evaluation of water-proofness quality and diffusion properties of coating materials[END_REF]. To sum up, as complete stone waterproofing is too radical, human action should be confined to limiting water penetration from the exterior, while allowing gaseous transfer in both directions. Furthermore, the current trend is to preserve the original stone without any chemical additives and with a minimum human action. To meet these various criteria, a bio-mineralisation treatment was developed. The initial method, developed by Aldolphe, Castanier and others used the ability of bacteria to induce calcite precipitation [START_REF]Biosynthèse de calcite par association bactérienne aérobie[END_REF][START_REF] Castanier | Participation bactérienne à la production du carbonate[END_REF][START_REF] Billy | Isolement des constituants d'une association bactérienne productrice de calcite[END_REF][START_REF] Kahoui | Identification of Early Bacillus-induced Crystals in vitro Using Fourier Transform Infrared Spectroscopy[END_REF]. This process has been improved and used industrially by the Calcite Bioconcept company. Other ways to produce calcite have also been investigated. The crystallisation of struvite and calcite by dead cells and cellular fraction of myxococcus xanthus were studied by Rodriguez-Navarro et al. [START_REF] Rodriguez-Navarro | Conservation of Ornamental Stone by Myxococcus xanthus-Induced Carbonate Biomineralization[END_REF]. The Biobrush consortium selected bacteria for their ability to firstly remove black crusts and then to precipitate calcite [START_REF]Biobrush research monograph : novel approaches to conserve our European heritage[END_REF], while the Bioreinforce Consortium recommends the application of organic matrix molecules to reinforce the stones. Recently Tiano applied the latter method in situ on marble [START_REF] Tiano | Biomediated reinforcement of weathered calcareous stones[END_REF]. A review of these bio-techniques can be found in Webster and more recently in De Muynck [START_REF] Webster | Bioremediation of weathered-building stone surfaces[END_REF][START_REF] Muynck | Microbial carbonate precipitation in construction materials: A review[END_REF].

The possibility of using this capacity in situ, in building conservation practice, was first studied by Levrel et al. [START_REF] Métayer-Levrel | Applications of bacterial carbonatogenesis to the protection and regeneration of limestones in building and historic patrimony[END_REF]. The bio-coating was successfully used for the first time on a surface area of 50 m 2 on the south-east tower of the Thouars church (France) in 1993 [START_REF] Métayer-Levrel | From Carbonatogenesis concepts to bacterial regeneration of limestones (microbial lifting). in: IAS-ASF-IGCP 380 International Workshop on "Microbial mediation in carbonate diagenesis[END_REF]. Since then, this treatment has been applied on numerous buildings, monuments and houses in Paris and all around France. The properties of the coating were determined by water absorption measurements (Karsten pipe), colorimetry, roughness measurement, macroscopic and microscopic observations of the stone surface and counting of the bacterial population [START_REF] Métayer-Levrel | Applications of bacterial carbonatogenesis to the protection and regeneration of limestones in building and historic patrimony[END_REF]. Two sets of measurements were done on Thouars church, six months and one year after the treatment. Stone powdering and roughness were decreased, no color change was observed and water absorption was decreased by a factor varying between four to eight times. The colorimetry test is important as a prerequisite for the restoration of historic buildings is that a treatment must not alter the color of the building [START_REF] Métayer-Levrel | Applications of bacterial carbonatogenesis to the protection and regeneration of limestones in building and historic patrimony[END_REF]. Furthermore, carbonate crystallization has been highlighted by other indirect methods, using different bacteria colonies [START_REF] De Muynck | Bacterial carbonate precipitation as an alternative surface treatment for concrete[END_REF][START_REF] Métayer-Levrel | Microbiogéologie du carbonate de calcium: applications industrielles, implications géologiques[END_REF]. Levrel sprayed limestone cubes with a bacteria solution and daily fed them (during five days) by spraying a nutrient solution onto the stones. The amount of carbonate produced was estimated by weighing the cubes before and after treatment. However, as pointed out by Levrel, the bacteria production was weighed with other crystallized salts present in the nutrient solution, making the absolute calcite quantification rather difficult. De Muynck et al. immersed mortar cubes and concrete cylinders for 24 hours in a stock culture of bacteria or covered them with ureolytic culture sludge prior to submersion in different compositions of culture media. The crystallographic form of the crystallized calcium carbonate by X-ray diffraction (XRD) revealed the presence of calcite and also small quantities of vaterite, which is strongly dependent on the composition of the culture media. Similar to many papers on biomineralization on limestone (see Rodriguez-Navarro et al.

[11], Jimenez-Lopez et et al. [START_REF] Jimenez-Lopez | Consolidation of degraded ornamental porous limestone stone by calcium carbonate precipitation induced by the microbiota inhabiting the stone[END_REF]) De Muynck et al. performed biomineralization experiments under optimal conditions in the lab (i.e. in aqueous environment) as proof-of-concept. Furthermore, their procedure improves the amount of produced carbonate rather than the bacterial colony growth.and needs to be adapted before it can be applied in practice. Castanier et al. used a bacterial sludge on plaster and pulverized a nutrient solution daily during 5 days. They observed a coating using optical microscopy and scanning electron microscopy, and showed that a modification of the surface was obtained [START_REF] Orial | Première approche de traitement de surface du plâtre par biominéralisation[END_REF].

The goal of this paper is to characterize and analyze the phase mineralogy produced by bio-treatment using the recommended industrial process already applied to Thouars church. In this process, Bacillus cereus was mainly selected for its rapid growth and important carbonatogenesis properties, as well as for its non-hazardous character for human health and for the environment. In order to limit the development of fungi present on the stone, a fungicide (Natamycine) was added to the nutrient solution (0.125g/l). Laboratory experiments showed that bacillus cereus was resistant to this amount of fungicide and was still able to produce carbonates while a larger quantity affected the bacteria development. Moreover, experimental tests showed that bacillus cereus was resistant to heavy metals [START_REF] Loubière | [END_REF]. In addition, the precipitate does not produce any stone color change [START_REF] Métayer-Levrel | Microbiogéologie du carbonate de calcium: applications industrielles, implications géologiques[END_REF].

The paper is organized as follows. In section 2, the choice of the substrate is discussed and then described as well as the bio-treatment process and the experimental techniques which were used to analyze the coating. The results are given in section 3. Discussion and conclusion follow in section 4.

Material and Methods

The bio-treatment used in this study was applied on plaster. This material is commonly used in building construction and is composed of gypsum crystals. These are tabular and prismatic with a typical size close to 10 micrometers. The plaster came from Soisy-sous-Montmorency and was marketed as "neige n • 2". XRD analysis indicates only the presence of gypsum without any additive. This material was chosen so as to enable the calcite production of the bacteria to be detected unambiguously. As the bio-treatment is expected to produce calcite, it would be difficult to distinguish deposited calcite from a substrate made of calcite such as limestone.

The bio-treatment used in this work, under patent by the Calcite Bioconcept firm, involved Bacillus cereus, a bacterium particularly suitable for limestones [START_REF] Métayer-Levrel | Microbiogéologie du carbonate de calcium: applications industrielles, implications géologiques[END_REF]. Since one of the aims of this technique is to produce a coating of the same nature as the stone substrate, this technique has often been applied on limestones such as tuffeau [START_REF] Dessandier | Relationships between mineralogy and porous medium of the craie tuffeau[END_REF][START_REF] Rozenbaum | 2-D image analysis: A complementary tool for characterizing quarry and weathered building limestone[END_REF].

For practical reasons, the bacteria were lyophilized in order to be easily used on a building site. The freeze-dried bacteria were first re-hydrated with the nutritive solution recommended by the Calcite Bioconcept firm. Fifteen hours later, the culture medium was sprayed onto the building facade (about 1 liter per m 2 ). The bacteria were fed with a nutritive solution 24, 32, 48 and 72 hours after spraying of the culture medium. During these three days the bacteria colony exponentially increased. This bio-mineralisation treatment was optimized on limestone so as to be completed within one week, which is one of the restorers' requirements.

Surface samples were firstly observed by using scanning electron microscopy (SEM); in addition a fresh fracture was also scanned in order to obtain depth information. All the specimens were coated with gold prior to observation using secondary electron emission.

X-ray diffraction is a physico-chemical analysis method commonly used to determine the reticular plane distance and the crystalline structure of samples. More particularly, Grazing Incidence X-ray Diffraction (GIXD) is well adapted to the structural analysis of coatings or surfaces of bulk samples [START_REF] Andreazza | In-depth structural X-ray investigation of PECVD grown diamond films on titanium alloys[END_REF]. The GIXD system used was configured with a standard Cu-Kα laboratory (λ = 0.15418 secondary parallel plate collimator followed by a flat graphite secondary monochromatic and proportional detector [START_REF] Andreazza | Buffer layer effect in nanostructured tin electrodeposition on insulating and conducting substrates[END_REF]. Diffraction spectra were obtained at fixed angles Ω (incidence angles through 2θ scan). The experiments presented in this paper were performed with a grazing incidence angle of θ = 0.3 • and θ = 1 • and a glancing incidence of θ = 15 • to modulate the analysis depth from the sample surface. For this experiment on massive samples, the plaster samples were moulded in order to obtain pellets (diameter: 20 mm thickness: 5 mm)

The following chemical elements were mapped with a micro-probe Camebax SX 50 on the surface samples: Calcium (Ca), Sulphur (S), Sodium (Na), Potassium (K), Phosphorus (P), Magnesium (Mg) and Chlorine (Cl). Each map consists of points regularly spaced on 512 rows and columns; for each point, the counting time was 100 ms, for a voltage of 15 kV and a current intensity of 50 nA. The following 7 samples were treated (1) and then analyzed:

• a reference sample (untreated: S1) that was treated after measurement (S1t),

• a treated sample (S2t), that was broken into two parts after measurement. The first part was treated a second time (S2tt) and the second part washed (S2tw),

• a treated sample for which the centre of the surface sample (S3t) and the border of the same surface sample (S3tb) were analyzed.

3 Results

SEM analysis

SEM micrographs for untreated plaster are given on figures 1(a) and (b) for magnifications of 1 000 and 5 000 respectively. They showed tangled gypsum needles giving a rather rough surface. After bio-treatment this surface was strongly modified and far less rough, with crusts covering the gypsum needles, producing a smooth surface plaster 2(a)). Furthermore, ovoid bacteria bodies (2 µm long and 0.5 µm wide) due to the mineralization of bacteria membranes (specifically carbonate mineralization [START_REF] Castanier | Ca-carbonates precipitation and limestone genesis -the microbiogeologist point of view[END_REF]) created reliefs on this smooth surface. Moreover, crystal growth was observed after bio-treatment as well as numerous cracks (Fig. 2(b)). These cracks are a common feature of biofilms observed by SEM, and have also been observed on limestone samples. However, this feature is probably an artifact due to the high vacuum needed in the SEM chamber.

GIXD analysis

As shown by SEM micrographs (Fig. 3), the treatment produced a coating. Consequently GIXD was used in order to obtain information about the structure modifications of the plaster material surface and especially the crystalline nature of the deposit layer. The samples were mounted on the sample holder in exactly the same position and orientation before and after bio-treatment. This precaution made it possible to obtain an identical substrate diffraction contribution before and after bio-treatment and thus to shed new light on the coating composition. This precaution was taken for all successive analyses on the same sample.

Figure 4 shows the GIXD spectra of an untreated (S1) and bio-treated (S1t) sample. The GIXD spectrum of the untreated plaster (S1) showed only typical gypsum peaks without any calcite contribution, whereas the treated sample also shows calcium carbonate peaks (here vaterite), halite and sylvite peaks. In addition, figure 4 exhibits the decrease in the gypsum surface contributions with respect to the untreated sample as a consequence of the formation of a coating by the treatment. Similar results were obtained for all the other analyzed samples, except that calcite was obtained instead of vaterite which was no longer detected.

To emphasize the effect of the treatment, it was applied on sample S2t a second time. The GIXD spectrum of this twice treated sample (S2tt) shows an increase in the amplitude of the calcite peaks, compared to the sample treated once only (S2t) (Fig. 5(a)). In addition the amplitude of the halite peaks also increases , those of the sylvite peaks remains almost the same, while the response of the substrate decreases as already observed for samples S1 and S1t.

High proportions of salts within the porous lattice of a building material can lead to severe damage due to the local force of the crystallization pressure (e.g. [START_REF] Benavente | Salt Influence on evaporation from porous building rocks[END_REF][START_REF] Scherer | Stress from crystallization of salt[END_REF][START_REF] Cardell | Weathering of limestone building material by mixed sulfate solutions. Characterization of stone microstructure, reaction products and decay forms[END_REF]). As rain is thought to discharge the salts from the surface, the coating was therefore washed in order to simulate the effect of rain. Sample S2t was immersed in distilled water for 30 seconds. The water was changed and the sample immersed again for 30 seconds. The sample was then dried in an oven during 12 days at 35 • C and then analyzed. Its GIXD spectra (Fig. 5(b)) showed lower halite peaks, in comparison with S2t spectra while sylvite peaks were completely removed. This can obviate the salt problems due to bio-treatment. Furthermore, calcite peaks remained, indicating a rather good durability of the coating.

In order to have information about both the surface homogeneity of the treatment and its time evolution, a sample was analyzed only 7 days after bio-treatment on two distinct places: at the centre of the surface sample (S3t) and at its border (S3tb) (Fig. 6). We observed that the peak intensity of the salts is higher at the border than in the center. This is in agreement with a drying of the sample from the border to the center, implying both a preferential salt crystallization on the border and diffusion of the concentration of the saline solution from the center to the border. This effect does not Figure 3: SEM image of a plaster sample section. Under the surface, needles of gypsum appear clearly while the coating is visible (on the top of the image). appear for the calcite for which the peak sizes are almost the same for the two places, showing that the calcite deposit is constant. It is also interesting to note that even after 7 days, we observed the same peaks (sample S3t) as for an older treated sample (S2t), indicating that the whole mineralogy is established at the very beginning of the process. Other samples, such as sample S4t, 17 days aging, yielded the same results.

It is also possible to obtain depth information by varying the grazing incidence. Increasing grazing incidence will also increase probe depth. Sample S2t was analyzed again with an incidence angle of θ = 15 • . Taking into account the X-ray absorption of gypsum, an analytical calculation gives an estimation of the probe depth of about 1 µm for θ = 1 • and 10 µm for θ = 15 • . The diffraction lines shown on figure 7 correspond to the contributions of phases in the 1 micrometer layer and 10 micrometer layer from the surface. We estimated the depth evolution of the calcite and salts phases from the variations in integrated intensities. The (h,k,l) peaks of calcite (113), sylvite (200), halite (200) and gypsum (022) were chosen because they are isolated without overlapping with other phases. The intensity ratios I calcite /I gypsum , I sylvite /I gypsum and I halite /I gypsum for both incidence angles were calculated. We obtained:

I calcite I gypsum 15 • = 0.3 I calcite I gypsum 1 • = 1.1 I sylvite I gypsum 15 • = 10
I sylvite I gypsum 1 • = 16
I halite I gypsum 15 • = 0.5 I halite I gypsum 1 • = 3.
The calcite/gypsum and salts/gypsum ratios decrease with depth, demonstrating that salt and calcite are present only in the superficial layer. As already mentioned, the validity of these ratios holds only for the few first 10 micrometers from the surface. Moreover, the decrease is extremely strong in the halite salt case compared to the sylvite salt case. The results suggest that the halite salts were located on the carbonate coating and not within the porosity of the gypsum substrate, whereas the sylvite is rather in an intermediate position layer.

Microprobe analysis

Microprobe analysis was performed on treated and untreated samples on sulphur and calcium. Figures 8(a) and (b) give microprobe cartographies of calcium and sulphur respectively for a polished section of a treated sample. Figure 8(a) reveals the presence of a rather homogeneous calcium layer, clearly visible on the surface (right side of 8(a)), with a thickness of roughly 10 µm. The comparison between figure 8(a) and (b) shows that this layer corresponds to an area in which almost no sulfur is present, meaning that this calcite layer was formed on the surface of the sample. However, this new calcite layer has a low intensity compared to values obtained inside the sample, and it is difficult to infer the depth penetration of the treatment from this information alone . The untreated sample does not present this feature and the distribution of both sulphur and calcium is homogeneous as expected. Although microprobe analysis does not give information of a mineralogical nature, it is a complementary tool to GIXD because it provides a depth description. However, the microprobe analysis confirms the surface results obtained by GIXD.

2 2 2 4 2 6 2 8 3 0 3 2 3 4 3 6 3 8 4 0 4 2 2 F ( °) ( K Q C u ) A r b i t r a r y i n t e n s i t y S 2 t S 2 t w K ( 2 0 0 ) C ( 1 0 4 ) N a ( 2 0 0 ) * C ( 0 1 2 ) C ( 1 1 0 ) N a ( 1 1 1 )

Water absorption measurements

Water absorption measurements were approached by measuring the absorption of a drop by the surface of a treated and untreated sample. 0.2 ml of distilled water were placed on the surface of the stone and the time necessary for complete water absorption was measured. These experiments were done on five points for each sample and the average time is as follows: for the untreated stone, absorption was almost immediate (less than 2 s) while for the treated stone the drop keeped its shape for about 20 s (like a water drop on a hydrophobic surface) and was completely absorbed after 30 s. It is clear that the water properties were modified by the treatment.

Discussion and conclusion

A bio-treatment was performed in the laboratory on plaster samples following an industrial process. The coating produced was observed by SEM, and analysed by microprobe X and GIXD. First of all, Bacillus cereus bacteria (and the established industrial protocol) produced calcite coatings on plaster surfaces. Compared to the sample treated once only (S2t), the sample treated twice (S2tt) presented more calcite coating, showing that calcite was created at each treatment. Moreover, the GIXD technique enabled the crystallographic composition of the coating to be clearly displayed, making it possible to analyze thin surface layers. Despite the fact that the bio-treated samples were analyzed with different aging times, the crystal mineralogy remained the same. Indeed, except in one case, the only polymorphic state of carbonate calcium observed was calcite. No transitory state such as vaterite, which has been recognized to be a transient polymorph of calcite, was observed [START_REF] Spanos | The transformation of vaterite to calcite: effect of the conditions of the solutions in contact with the mineral phase[END_REF][START_REF] Katsifaras | Effect of inorganic phosphate ions on the spontaneous precipitation of vaterite and on the transformation of vaterite to calcite[END_REF][START_REF] Wolf | Thermodynamic aspects of the vaterite-calcite phase transition[END_REF]. In one case, vaterite was found (aging: 53 days) but was not found again on other samples. A possible explanation for this could be that the environmental conditions of treatment were insufficiently controlled (humidity, temperature . . . ) in the laboratory, leading to vaterite formation. However these variations could also be expected to occur in an outdoor environment and consequently to produce different carbonate polymorphs on treated building frontages.

The main problem caused by the process was salt deposit (coming from the feeding solution) on and within porous materials. This study shows that simply washing the sample (rain in the real case of buildings) sufficed to remove these salts to a large extent. Nevertheless the calcite coating was still present and was not removed. Then, once washed, the bio-treated surface was almost salt-free and the surface porosity should be reduced due to the calcite deposit. However durability tests should be carried out in order to control the efficiency of the treatment over several years.

GIXD spectra obtained on different areas of the same sample (S3 and S3b) show a higher concentration of halite and sylvite on the edges than in the center. This could be explained by the fact that the sample dried more quickly at the edges than in the center, creating favorable thermodynamic conditions for salt precipitation. This border effect, due to the shape of the building blocks could also appear in situ.

Initial measurements on water absorption indicate that water penetration was reduced for a treated sample, suggesting that the surface porosity was reduced, though this point needs to be controlled by porosity measurements. However, as the coating is 10 µm thick, intrusive porosity measurements (e.g. mercury porosimetry) would require extracting 10 µm thick slices in order to follow the depth porosity change. This is obviously not possible (without any resin induration) and an efficient way to measure this porosity can be obtained by image analysis. This requires access to X-ray tomography with high resolution (0.25 µm) only accessible nowadays on synchrotrons. This study will be the subject of future work and of a forthcoming publication.

Finally, as water absorption is reduced, this conservation technique appears to be perfectly efficient on plaster and this study should be extended to other porous building materials. Indeed, as this industrial protocol was recommended for limestones, future studies will be carried out on these building materials, which will involve addressing the difficulty of distinguishing bio-calcite from calcite limestone. A preliminary SEM study on bio-treated tuffeau showed the same kind of surface coating as the plaster coating. In addition, this study proves that grazing incidence X-ray diffraction, commonly used for low roughness (< 10 nm) material surfaces, is well adapted to highly porous and granular mineral material.

Figure 1 :

 1 Figure 1: SEM image (secondary electron) of a plaster surface before bio-treatment (a) 1000×magnification, (b) 5000×magnification.

  Figure 2: SEM image (secondary electron) of a plaster surface after bio-treatment (a) 1000×magnification, (b) 5000×magnification.

Figure 4 :

 4 Figure 4: GIXD spectra for θ = 0.3 • of sample 1. Continuous line untreated sample (S1) and dotted continuous line bio-treated (S1t). Gypsum, Vaterite, Halite and Sylvite peaks are labeled with stars (*) for the former and capital letters V, Na and K respectively for the others. The associated numbers in parentheses correspond to diffraction (h, k, l) lines.

Figure 5 :

 5 Figure 5: (a) GIXD spectra for θ = 0.3 • of sample 2 treated once (S2t) and treated twice (S2tt). (b) GIXD spectra for θ = 1 • of sample 2 treated once (S2t), then washed (S2tw). The peaks are labelled as in figure 4 and letter C stands for Calcite.

Figure 6 :

 6 Figure 6: GIXD spectra for θ = 0.3 • of sample S3t and S3tb analyzed on the center (dotted continuous line) and on the border (continuous line) respectively. The dotted line gives the GIXD spectra of the untreated sample S1 as a reference. The encapsulated picture is a zoom of the spectra. The peaks are labelled as in figure 5.

Figure 7 :

 7 Figure 7: GIXD spectra for θ = 1 • (continuous line) and θ = 15 • (dotted continuous line) of sample 2 treated once (S2t). The peaks are labelled as in figure 5.

Figure 8 :

 8 Figure 8: Microprobe cartographies of calcium (a) and sulphur (b) of a section of a treated plaster sample. The treated surface corresponds to the right side of the pictures, indicated by arrows. The microprobe measured intensity is given next to the color scale of the pictures.

Table 1 :

 1 Sample treatments and labels.

	Samples Bio-treatment	Elapsed days before analysis Label
	1	none	-	S1
		treated once	50	S1t
	2	treated once	110	S2t
		S2t treated once again	40	S2tt
		S2t washed	10	S2tw
	3	treated once	7	S3t
		treated once	10	S3tb

Acknowledgments

We thank J.-F. Loubière from Calcite Bio Concept firm for his support and for supplying the biological material as well as for interesting discussions. We also thank F. Bousta and G. Orial from the LRMH (Laboratoire de Recherche sur les Monuments Historiques) for their valuable comments. SEM images were taken at the "Centre de Microscopie électronique" of the university of Orléans.


Ronan Vicquelin
email: ronan.vicquelin@centraliens.net

Benoit Fiorina
email: benoit.fiorina@em2c.ecp.fr

Sandra Payet
email: sandra.payet@gdfsuez.com

Nasser Darabiha
email: nasser.darabiha@em2c.ecp.fr

Olivier Gicquel
email: olivier.gicquel@em2c.ecp.fr

Coupling tabulated chemistry with compressible CFD solvers

Keywords: tabulated chemistry, compressible flow, characteristic boundary conditions

Introduction

Large Eddy Simulation (LES) has become an affordable tool. This major evolution is due to the important improvement of the available computational power during the last decade. Using parallel computational resources, Boileau et al. [1] were able to simulate the ignition of a full combustion chamber. Similar simulations of realistic systems were also performed by [2,3]. All these works are of great interest to study turbulent flame dynamics, combustion instabilities or mixing characterization but they are limited in terms of chemistry description.

Indeed, because of large number of species, detailed chemistry simulations remain too expensive in terms of CPU time to be achieved in such configurations. Some direct numerical simulations have been done using detailed chemistry mechanisms for hydrogen-air flames [START_REF] Mizobuchi | Proceedings of the Combustion Institute[END_REF] or methane-air flames [START_REF] Hawkes | [END_REF] but they are limited to very small configurations that are far from industrial needs.

A classical approach used to take into account detailed chemistry effects in realistic LES for a low CPU cost is to use tabulated chemistry. Tabulated chemistry methods assume that chemical evolutions in the composition space can be parameterized and tabulated by a reduced set of variables (ψ 1 , ..., ψ n ) where n is the number of chemical database coordinates. In general, ψ l are combinations of species mass fractions. Knowing these variables in a simulation, all thermo-chemical variables ϕ can then be estimated by using the chemical database ϕ tab (ψ 1 , ..., ψ n ). Among these tabulation techniques one can mention ILDM [6], FPI [START_REF] Gicquel | Proceedings of the Combustion Institute[END_REF], FGM [START_REF] Van Oijen | [END_REF], REDIM [9] or ICE-PIC [10] methods. In these methods, instead of solving one balance equation for each chemical species involved in the detailed elementary reactions, only few equations for the reduced set of variables are solved.

The main problem when coupling this procedure with compressible Navier-Stokes equations is that perturbations due to compressibility effects are not considered during the database generation. This limitation is not a problem if the CFD solver is based on a low Machnumber assumption where the introduction of tabulated chemistry is straightforward [11,12,13]. But introducing tabulated chemistry in a compressible CFD solver is more challenging. In this case, temperature can not be directly read in the database because the tabulated temperature does not take into account acoustic perturbations. Galpin et al. [14] have proposed to solve, in addition to the reduced set of variables, balance equations for some chemical species selected to estimate the temperature from energy. For methane-air combustion, the number of supplemental species required to have a fair estimate of the temperature is nine [14]. Therefore this technique leads to a large increase of the number of equations. In addition, a divergence between the additional transported species and the tabulated ones is frequently observed and requires specific treatment. A second problem in compressible codes is the interaction between acoustics and boundaries [15,16]. When using tabulated chemistry, boundary conditions for the reduced set of variables should therefore take into account perturbations due to acoustics.

In the present work, a new technique called TTC (Tabulated Thermochemistry for Compressible flows) is developed to introduce tabulated chemistry strategies in compressible solvers. This approach only requires to solve additional balance equations for the database coordinates. First, the governing equations of reactive flow described by tabulated chemistry are briefly presented and a method is proposed to take into account temperature deviation due to acoustics. Then characteristic boundary condition treatment is detailed for tabulated chemistry. Finally, validation tests are presented. Reformulation of temperature computation and boundary treatment are validated with 1-D tests by comparison with multi-component simulations. Further validations are conducted by performing a non-reactive 3-D large eddy simulation of a round jet and a reactive case where the FPI tabulation method is chosen to reproduce a 1-D laminar flame.

Governing equations for tabulated chemistry

Conservative variables are gathered in the vector U = (ρ, ρu i , ρe t , ρψ 1 , ..., ρψ n ) ⊤ where ρ is the mass density, u i are the velocity components, e t = u i u i /2 + e is the total energy, sum of kinetic and internal energy, and ψ l is the l th coordinate used to describe the thermochemical database. In multi-component formulation, mixture internal energy e is given as e = N k=1 e k Y k where N is the number of species and Y k , e k are the k th species mass fraction and energy, respectively. At a given temperature T , e k reads:

e k = T T 0 C vk (T )dT - RT 0 W k + ∆h 0 f,k (1)
where R = 8.314 J/mol/K is the ideal gas constant, C vk is the species specific heat capacity at constant volume, W k is the k th species molar mass, ∆h 0 f,k is the species enthalpy of formation and T 0 is the reference temperature (T 0 = 298 K).

When using tabulated chemistry, governing equations for adiabatic reactive flows can be written as:

∂U ∂t + ∂F j ∂x j + ∂F j d ∂x j = S (2)
where the source term vector is S = (0, 0, 0, 0, 0, ωψ 1 , ..., ωψ n ) ⊤ with ωψ l the chemical mass production rate of the variable ψ l . F j and F j d are respectively the Eulerian and diffusion fluxes in direction j:

F j =                                       ρu j ρu j u 1 + δ 1 j P ρu j u 2 + δ 2 j P ρu j u 3 + δ 3 j P u j (ρe t + P) ρu j ψ 1
where δ i j is the Kronecker operator: δ i j = 1 for i = j and δ i j = 0 for i j. P is the thermodynamic pressure given by the ideal gas law, P = ρRT/W, where

W = N k=1 Y k /W k -1
is the mixture molar mass. τ i j is the viscous tensor, q j is the heat flux and ψ l follows a Fick diffusion law with diffusivity coefficient D l . The fluid is considered as Newtonian. The heat flux vector is given using Fourier law and assuming unity Lewis number for species diffusion:

q j = - λ C p ∂h ∂x j (4)
where h = e + P/ρ is the mixture enthalpy, λ is the mixture thermal conductivity. C p and C v are the mixture specific heat capacities at constant pressure and constant volume respectively. The ratio of heat capacities γ = C p /C v is introduced. Diffusivity coefficients D l are computed by assuming unity Lewis number: D l = D = λ ρC p . The coupling strategy developed in this article is not affected by this assumption.

In tabulated chemistry framework, C p , C v , λ and D are stored in a look-up table as a function of (ψ 1 , ..., ψ n ). A technique adapted to tabulated chemistry to estimate the temperature from the transported energy is proposed in the following section.

Temperature computation with tabulated chemistry

Conservative variables ρ, ρu i , ρe t and ρψ l are solved from their respective balance equations. The remaining difficulty is the computation of temperature. In low Mach-number approach, as pressure fluctuations are assumed to be small, the temperature is directly extracted from the database: T = T tab (ψ 1 , ..., ψ n ). This assumption is not accurate when compressibility effects are taken into account. Indeed, in this case, the transported values of e and T are different from the tabulated ones: e tab and T tab . A solution to compute the temperature is to approximate the fluid composition by a reduced sample of species as proposed by Galpin et al. [14]. It prevents a prohibitive storage in the database of all species thermodynamical properties and mass fractions. However, as mentioned previously, this method presents several drawbacks. Another solution is to tabulate the energy as a function of database coordinates and temperature. An additional coordinate (the temperature) needs therefore to be added to the chemical table and an iterative algorithm is required to deduce the temperature from energy. Such an approach is memory space and CPU time consuming. A less expensive alternative is here proposed.

The first order truncated Taylor expansion of e around T = T tab leads to the following linear approximation to estimate the temperature in TTC formalism:

T = T tab (ψ 1 , ..., ψ n ) + e -e tab (ψ 1 , ..., ψ n ) C tab v (ψ 1 , ..., ψ n ) (5)
The "compressible" temperature T can therefore be approximated from the tabulated energy and temperature e tab and T tab , respectively and the transported energy e = e t -u i u i /2. Note that this assumption is valid in the case of moderate acoustic perturbations which induce small temperature variations. It is justified when combustion operates at constant pressure such as in gas turbines, furnaces, unconfined laboratory flames, etc... Other realistic situations exist where pressure variations are sufficiently large to affect the chemistry (internal combustion engines, detonation waves, ...) and extracoordinates have to be added to the chemical database.

Characteristic boundary conditions for tabulated chemistry

Characteristic wave treatment

In compressible CFD solvers, boundary conditions are provided using characteristic boundary con-ditions. Navier-Stokes Characteristic Boundary Conditons (NSCBC) have been first derived for singlecomponent flows [15], and later for multi-component flows [16]. As N species transport equations have been replaced by n < N equations for variables ψ l , characteristic wave decomposition must be expressed in terms of n characteristic waves corresponding to ψ l and not N species waves anymore in TTC formalism. For that purpose, Eq. ( 2) is written introducing the vector of nonconservative variables V = (ρ, u i , P, ψ 1 , ..., ψ n ) ⊤ :

∂V ∂t + A j ∂V ∂x j + J -1 ∂F j d ∂x j = J -1 S (6)
where A j is a (n + 5) × (n + 5) matrix:

A j =                                       u j ρδ 1 j ρδ 2 j ρδ 3 j 0 0 ••• 0 0 u j 0 0 1 ρ δ 1 j 0 ••• 0 0 0 u j 0 1 ρ δ 2 j 0 ••• 0 0 0 0 u j 1 ρ δ 3 j 0 ••• 0 0 ρa 2 δ 1 j ρa 2 δ 2 j ρa 2 δ 3 j u j 0 ••• 0 0 0 0 0 0 u j ••• 0 . . . . . . . . . . . . . . . . . . . . . . . . 0 0 0 0 0 0 ••• u j                                      
. a = (γrT ) 1/2 is the speed of sound and r = R/W. The jacobian matrix J = ∂U ∂V reads:

J =                                       1 0 0 0 0 0 ••• 0 u 1 ρ 0 0 0 0 ••• 0 u 2 0 ρ 0 0 0 ••• 0 u 3 0 0 ρ 0 0 ••• 0 e t -rT β ρu 1 ρu 2 ρu 3 1 β -ρ β ϑ ψ 1 ••• -ρ β ϑ ψn ψ 1 0 0 0 0 ρ ••• 0 . . . . . . . . . . . . . . . . . . . . . . . . ψ n 0 0 0 0 0 ••• ρ                                      
where β = γ -1 and ϑ ψ l is defined as:

ϑ ψ l = N k=1 RT W k -βe k ∂Y k ∂ψ l (7)
Assuming that the boundary normal direction is along the x 1 axis, Eq. ( 6) is written as:

∂V ∂t + A 1 ∂V ∂x 1 = RHS ( 8
)
where RHS= J -1 S -J -1 ∂F j d ∂x j -A 2 ∂V ∂x 2 -A 3 ∂V ∂x 3 . Matrix A 1 has 5 + n eigenvalues λ m which are given here with one example of right eigenvectors r m :

λ 1 = u 1 + a, r 1 = ( ρ 2a , 1 2 , 0, 0, ρa 2 , 0, ..., 0) ⊤ λ 2 = u 1 , r 2 =
(1, 0, 0, 0, 0, 0, ..., 0) ⊤ λ 3 = u 1 , r 3 = (0, 0, 1, 0, 0, 0, ..., 0) ⊤ λ 4 = u 1 , r 4 = (0, 0, 0, 1, 0, 0, ..., 0) ⊤

λ 5 = u 1 -a, r 5 = ( ρ 2a , - 1 2
, 0, 0, ρa 2 , 0, ..., 0) ⊤ λ 6 = u 1 , r 6 = (0, 0, 0, 0, 0, 1, ..., 0) ⊤ ... λ 5+n = u 1 , r 5+n = (0, 0, 0, 0, 0, 0, ..., 1) ⊤

The right eigenvectors r m are gathered in the transformation matrix R in columns. Inverting R gives the matrix L whose lines are the left eigenvectors l m . Amplitude time variations of characteristic waves are defined as L m = λ m l m ∂V ∂x 1 [17] and are gathered in the vector L :

L =                                           λ 1 ∂u 1 ∂x 1 + 1 ρa ∂P ∂x 1 λ 2 ∂ρ ∂x 1 -1 a 2 ∂P ∂x 1 λ 3 ∂u 2 ∂x 1 λ 4 ∂u 3 ∂x 1 λ 5 -∂u 1 ∂x 1 + 1 ρa ∂P ∂x 1 λ 6 ∂ψ 1 ∂x 1 ... λ 5+n ∂ψ n ∂V ∂t + RL = RHS (10)
Eq. ( 10) specifies boundary conditions on the characteristic waves amplitudes L m . Indeed, if transverse, diffusion and source terms are neglected, i.e. RHS = 0, the Local One Dimensional Inviscid system [15] is retrieved. Once L m are known, conservation equations are expressed to provide Navier-Stokes Characteristic Boundary Conditions (NSCBC) at the boundary face:

∂U ∂t + J RL + ∂F 2 ∂x 2 + ∂F 3 ∂x 3 + ∂F j d ∂x j = S (11)

Tabulation of ϑ ψ l

The terms defined by Eq. ( 7) need to be computed at the boundary but can not be directly tabulated since they depend on temperature and energy which are both sensitive to compressible effects. Consequently, as done for temperature correction in Eq. ( 5), ϑ ψ l is compared to its corresponding tabulated value. The difference

∆ϑ ψ i = ϑ ψ i -ϑ tab ψ i (ψ 1 , ..., ψ n
) is approximated following a linear approximation by:

∆ϑ ψ l ≈ N k=1 R W k -βC vk ∂Y k ∂ψ l T -T tab (ψ 1 , ..., ψ n ) (12) Therefore, ϑ tab ψ l and σ tab ψ l = N k=1 (R/W k -βC vk ) ∂Y k
∂ψ l are first stored in a look-up table as a function of (ψ 1 , ..., ψ n ). Then to introduce compressible effects in the simulation, ϑ ψ l is computed as:

ϑ ψ l = ϑ tab ψ l (ψ 1 , ..., ψ n ) + σ tab ψ l (ψ 1 , ..., ψ n )(T -T tab ) (13)

Validation tests

The present method to introduce tabulated chemistry into a compressible CFD solver is validated. For that purpose, the temperature and characteristic boundary condition corrections presented in sections 3 and 4 are tested. First, the temperature correction method is validated by performing a one-dimensional simulation in the case of an acoustic wave traveling a gas mixture at rest. Then, the NSCBC terms corrections are validated by sending acoustic and entropic waves through boundaries. Finally, a turbulent jet and a laminar premixed flame are simulated to illustrate the model performances. For all these validations two simulations are performed. The first one uses multi-component transport formulation for the four species present in the mixture (reference simulation), while the other one uses tabulated chemistry following the TTC formalism. All simulations are done with the compressible CFD solver AVBP [19] using a third-order numerical scheme [20].

In the database used for non-reactive test cases, all thermo-chemical quantities are stored in a look-up table in term of a unique coordinate: n = 1 and ψ 1 = z where z is the mixture fraction. This database represents a rich methane-air mixture (z = 1 with equivalence ratio φ = 4.4) injected in a hot vitiated air (z = 0) coflow produced by the lean combustion of a hydrogen-air mixture (φ = 0.4). Table 1 details species mass fractions and temperature on the fuel side (Y k f and T f ) and on the oxidizer side (Y k ox and T ox ). Species mass fractions and temperature are stored in the database as:

Y k (z) = (Y k f -Y k ox )z + Y k ox ( 14
)
and

T (z) = (T f -T ox )z + T ox ( 15
)
In the present case only one balance equation for the mixture fraction z is added to the mass, momentum and energy balance equations.

Temperature correction test

The compressible deviation correction of temperature introduced by Eq. ( 5) is validated by computing an acoustic wave traveling across a periodic domain of length L = 0.005 m in a pure oxidizer mixture (z 0 = 0, Y k 0 , T 0 , P 0 = 1 atm). The initial solution is given by:

u = ±A exp - (x -x 0 ) 2 d 2 , P = P 0 + ρ 0 a 0 u (16) ρ = ρ 0 + ρ 0 a 0 u , T = P ρr
where ρ 0 = 0.24 kg/m 3 and a 0 = 734.6 m/s are the mass density and sound speed in the unperturbed initial solution. The sign of the velocity is chosen so that the wave travels in the domain towards positive x values. Periodic boundary conditions are prescribed. Temperature time evolution at x = 0 is plotted in Fig. 1. The same temperature field is predicted by both tabulated chemistry and multi-component simulations. Similar conclusions are observed for pressure, velocity and mass density time evolutions.

Characteristic boundary conditions validation

The NSCBC modifications for tabulated chemistry described in section 4 are tested by injecting acoustic and entropic waves through inflow and outflow boundary conditions. The wave amplitudes L m are imposed following the procedure described in Ref. [16].

In the proposed test cases, the initial mixture fraction is given by:

z = z 0 + z ′ , with z 0 = 1 - x L ( 17
)
where z ′ is a perturbation. When the flow is not perturbed by acoustic, species and temperature profiles are deduced using Eqs. ( 14) and ( 15), respectively. Figure 2 shows species and temperature profiles in this situation and without mixture fraction perturbation (z ′ = 0).

Acoustic wave through subsonic non-reflective outlet (right).

An acoustic wave defined by Eq. ( 16) with parameters: A = 0.01a 0 , d = 2.5 10 -3 m and x 0 = 0.025 m is superimposed. In this first case, the mixture fraction field is not perturbed: z ′ = 0. Fig. 3 shows velocity and pressure profiles passing through the right boundary. Identical solutions are observed in both multicomponent and tabulated chemistry simulations.

Entropic wave through subsonic non-reflective outlet (right). The flow velocity is now uniformly initialized: u = 20 m.s -1 and P = 1 atm. The mixture fraction is perturbed by:

z ′ (x, t = 0) = A exp - (x -x 0 ) 2 d 2 (18)
where A = -0.5, d = 2.5 10 -3 m and x 0 = 0.025 m. Fig. 4 shows methane mass fraction, temperature and density at different instants. Species mass fractions are transported in the multi-component simulation and extracted from the look-up table in the tabulated chemistry simulation. The entropic wave exits properly the computational domain without introducing numerical artifacts in both multi-component and tabulated chemistry simulations.

3-D large eddy simulation of a round jet in a coflow

Validation of the TTC method is addressed here on a three-dimensional simulation of a round jet impinging in a vitiated coflow. The configuration is identical to Ref. [18]. The fuel jet diameter d is 4.57 mm and the bulk velocity V jet is 100 m/s (Reynolds number is 28,000). The coflow velocity is 5.4 m/s. Composition of the fuel jet and the coflow are the same as specified in Tab. 1. Large Eddy Simulation (LES) is performed with multi-component and tabulated chemistry formulations without taking combustion into account. Therefore, only the mixing phenomenon is studied here. The subgrid LES flux is modeled using the WALE model [21]. The domain of length L = 0.7 m and of diameter 0.35 m is discretized on a non-uniform mesh composed of 30 million tetrahedra.

An instantaneous mixture fraction isosurface obtained by the tabulated chemistry method is represented in Fig. 5. Instantaneous data are averaged during a physical time equal to one convective time L/V jet to compute mean and root-mean-square (rms) quantities. Figure 6(a) shows that mean axial temperature and rms profiles obtained by both simulations are similar. Same behavior are observed in Fig. 6(b) for radial profiles of mean species mass fraction.

Error in temperature estimation due to the first order truncation of the energy Taylor expansion in Eq. ( 5) is determined. The bias B = (T T TC -T re f )/T re f is introduced. T T TC and T re f are both computed from the transported energy e of the tabulated chemistry simulation. T T TC is estimated from Eq. (5) while T re f is computed using a more expensive iterative algorithm without truncation assumption. In the present 3D LES, the maximum error in the whole computational domain remains very low: Max(B)=0.3%.

One-dimensional laminar premixed flame with tabulated chemistry

A simulation with combustion is now considered. The TTC method is compatible with any tabulated chemistry model. It is here illustrated with the FPI method [START_REF] Gicquel | Proceedings of the Combustion Institute[END_REF]13]. The objective is to simulate a onedimensional laminar premixed flame with tabulated chemistry method in a compressible CFD code. The reference solution with detailed chemistry is computed with the PREMIX code [22] for a stoichiometric methane-air mixture. In PREMIX, pressure is fixed to 1 atm and fresh gas temperature is 298K. Transport is described with Fick law and unity Lewis numbers for all species. The detailed mechanism GRI 3.0 [23] which involves 53 species is used. The flame structure is tabulated following the FPI methodology [START_REF] Gicquel | Proceedings of the Combustion Institute[END_REF]. All thermochemical quantities are stored in a look-up table in term of a unique coordinate: n = 1 and ψ 1 = c where the progress variable c is defined as the coordinate of the chemical database:

c = Y CO 2 Y CO 2eq (19)
where Y CO 2eq is the equilibrium value of CO 2 in burnt gases. AVBP initial solution is obtained by interpolating the detailed chemistry solution profile of c on a mesh with a length L = 0.04 m. At the beginning of the computation, an acoustic bump is created and exits the domain because the isobaric CHEMKIN solution is not compatible with a fully compressible formulation.

A convective time is defined as τ conv = L/S l . The tabulated chemistry simulation runs over a physical time equal to 0.2 τ conv = 27 ms. Comparisons between PRE-MIX and AVBP solutions are shown in Fig. 7. Species mass fractions and temperature match perfectly, demonstrating the performance of the temperature and boundary condition corrections proposed. Moreover, CPU saving is important since 53 species transport equations are replaced by one balance equation for progress variable.

In this reactive case, the maximum temperature estimation error is Max(B)=0.15%. Figure 7(a) shows the pressure profile accross the flame front. The pressure drop ∆P due to the variation of momentum flux ∆(ρu 2 ) across the deflagration wave is clearly visible.

Conclusion

In the present study, we have developed a new technique called TTC (Tabulated Thermochemistry for Compressible flows) to accurately couple tabulated chemistry with compressible flow solvers. A corrective term is proposed to take into account the deviation between the temperature stored in the chemical database and the one calculated in the CFD solver, due to compressible effects. NSCBC boundary conditions are adapted to tabulated chemistry formalism. It is shown that acoustic and entropic waves propagations are very well predicted. A LES of a turbulent round jet shows that the mixing is well described using this new coupling technique with the same accuracy as the multi-component formalism. Finally, a laminar flame simulation shows that the method perfectly reproduces premixed combustion phenomena. This new procedure to efficiently introduce tabulated chemistry in compressible flow solvers enables to take into account detailed chemistry effects in realistic LES configurations.

Fuel side (z = 1) Oxidizer side (z = 0)

T f = 320 K T ox = 1350 K Y O2 f = 0.195 Y O2 ox = 0.142 Y N2 f = 0.591 Y N2 ox = 0.758 Y H2O f = 0.0 Y H2O ox = 0.1 Y CH4 f = 0.214 Y CH4 ox = 0.0 φ CH 4 = 4.4 φ H 2 = 0.4

Figure 1 :Figure 2 :

 12 Figure 1: Time variation of temperature at x = 0 in the acoustic wave propagation test case. The initial gaussian profile is parametrized by A = 0.01a 0 , d = 5 10 -4 m and x 0 = 0.0025 m. Symbols: multicomponent solution. Line: tabulated chemistry solution.

Figure 3 :

 3 Figure 3: Acoustic wave passing through a non-reflective boundary. Pressure (left) and velocity (right) profiles are plotted for multicomponent (symbols) and tabulation (line) formalism at different times t + = tL/a 0 (x = L/2) in the simulation. (a) Initial solution, t + = 0 ; (b) t + = 0.4 ; (c) t + = 0.51.

Figure 4 :

 4 Figure 4: Composition wave passing through a reflective boundary. Methane mass fraction, temparature [K] and mass density [kg/m 3 ] profiles are plotted for multi-component (symbols) and tabulation (line) formalism at different times t + = tL/u 0 in the simulation. (a) Initial solution, t + = 0 ; (b) t + = 0.27 ; (c) t + = 0.45 ; (d) t + = 0.54 ; (e) t + = 0.63.

Figure 5 :Figure 6 :Figure 7 :

 567 Figure 5: Mixture fraction isosurface z = 0.1 colored by longitudinal velocity.

Table 1 :

 1 Temperature, species mass fractions and equivalence ratio φ on fuel and oxidizer sides in the Cabra burner[18].

Acknowledgment

This work was supported by the project ANR-PANH and performed using HPC resources from GENCI-IDRIS (Grant 2009-x2009020164).


Pierre Auzillon
email: pierre.auzillon@em2c.ecp.fr

Benoit Fiorina

Ronan Vicquelin

Nasser Darabiha

Olivier Gicquel

Denis Veynante

Modeling chemical flame structure and combustion dynamics in LES

Keywords: Large Eddy Simulation, Turbulent premixed combustion, Tabulated chemistry

Introduction

Large Eddy Simulation is an attractive tool to predict flame dynamics and pollutant emissions in industrial combustion chambers [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]. The primary challenge is then to capture simultaneously the flame front propagation speed, the thermal expansion and the chemical species formation.

These phenomena are first driven by chemical kinetics and require therefore an accurate description of combustion chemistry. A widely used approach to describe fluid/chemistry interactions at a reduced computational cost is to tabulate chemical species mass fractions and/or reaction rates as functions of a reduced set of variables. A large variety of techniques has been developed to construct chemical databases. Among them, the ILDM (Intrinsic Low Dimensional Manifold) is based on a direct mathematical analysis of the dynamic behavior of the nonlinear chemical system response [START_REF] Maas | Proc. Combust. Inst[END_REF]. An alternative strategy is to map the chemical flame structure to a reduced phase subspace from elementary combustion configurations [START_REF] Gicquel | Proc. Combust. Inst[END_REF][START_REF] Van Oijen | [END_REF]. For instance, in the Flame Prolongation of ILDM method [START_REF] Gicquel | Proc. Combust. Inst[END_REF]5], the chemical subspace of a turbulent premixed flame can be approximated by a collection of 1-D laminar flames. In such simple situations, all thermochemical quantities are related to a single progress variable related to the temperature or to a linear combination of chemical species. This tabulation chemistry method has been found very efficient and successful in laminar premixed and partially-premixed flames simulations [START_REF] Gicquel | Proc. Combust. Inst[END_REF]6,7] and to compute turbulent premixed and partially-premixed flames in both RANS (Reynolds Averaged Navier-Stokes equations) [8,9,10] and LES contexts [11,12]. As recently shown by Nguyen et al. [13], very complex flame structures such as triple flames can also be accurately recovered by tabulated chemistry techniques.

Unfortunately, an accurate description of the combustion chemistry is not sufficient to predict the correct flame front behavior in Large Eddy Simulations. The primary recurrent problem is that the flame thickness is typically thinner than the LES grid size. To overcome this issue, different strategies have been developed. The level-set or G-equation approach consists in tracking the inner layer by solving a propagation equation. Initially developed for RANS [START_REF] Peters | Turbulent Combust[END_REF], the G-equation has been reformulated for LES [START_REF] Menon | [END_REF]16,17]. However as level-set techniques provide information only on the thin reaction zone position and not on the filtered flame structure, the coupling with the flow equations remains challenging. In particular, an additional modeling effort is needed to estimate the temperature field and thermal expansion [18]. Moreover, the coupling between G-equation formalism and tabulated chemistry is a difficult task.

A common way to couple tabulated chemistry with turbulence in RANS [8,9,10] and in LES [19] is to presume the shape of the sub-grid-scale probability density function (PDF). Unfortunately, when extended to LES, this formulation does not guarantee a proper description of regimes where sub-grid-scale flame wrinkling vanishes [12]. These situations, observed when the sub-grid fluctuations are lower than the laminar flame speed, are encountered in practical LES of premixed combustion [12,18,19]. Additionally LES should tend toward DNS when the filter size becomes lower than the Kolmogorov scale.

A solution to propagate a flame front on a coarse grid is to artificially thicken the flame front (Thickened Flame model for LES) and to model subgrid scale wrinkling [20]. In premixed combustion, as this technique can be applied to the progress variable balance equation, coupling TFLES and tabulated chemistry is straightforward [21]. An alternative is to introduce a filter larger than the mesh size in order to resolve the chemical flame structure. Initially developed for infinitely thin flame front and for single-step chemistry [22,23,24], this strategy has recently been extended to tabulated chemistry in the Filtered Tabulated Chemistry for LES (F-TACLES) model [12,25]. As the flame front is artificially spread in TFLES and F-TACLES approaches, flame structure and flow field may be affected.

The objective of this present work is to compare the effects of flame front thickening and filtering on both flame structure and dynamics. To avoid uncertainties related to the modeling of sub-grid-scale turbulence and flame interactions, the comparison is conducted in situations where the flame front is not wrinkled at the subgrid-scale level. TFLES and F-TACLES mathematical formalism are first briefly described in situations where flame wrinkling is fully resolved on the LES meshes (section 2). A criterion based on reaction zone thickness is proposed in section 3 to achieve quantitative comparison between the two competitive approaches. Finally, TFLES and F-TACLES models are applied to the simulation of an unsteady pulsed 2-D laminar Bunsen flame for various thickening factors and flame filter sizes, respectively (section 4). Results are compared to a filtered DNS solution serving as a reference solution.

Coupling tabulated chemistry with LES

Problem formulation

The FPI technique is retained here to determine lowdimensional trajectories covered by combustion chemistry [START_REF] Gicquel | Proc. Combust. Inst[END_REF]5]. For premixed combustion, one-dimensional freely propagating flames are first computed using detailed chemical schemes. Thermodynamical and chemical quantities are assumed to depend on a unique monotonic progress variable c, where c = 0 corresponds to fresh gases and c = 1 to fully burnt gases respectively. For instance, progress variable reaction rate is directly given by a look-up table

: ωc = ωTab c [c].
The chemical database is then coupled to the flow field by adding a balance equation for the progress variable to the Navier-Stokes equations. For LES, under unity Lewis Number assumption, this equation is filtered leading to:

∂ ρ c ∂t +∇•( ρ u c) = ∇• ρD∇c -∇• ρ uc-ρ u c + ωc (1)
where ρ is the density, u the velocity vector, D is the diffusivity coefficient and ωc is the progress variable reaction rate. The overbar denotes the spatial filtering operation:

φ(x) = R 3 G(x -x ′ )φ(x ′ )dx ′ (2)
where φ represents reactive flow variables and velocity components and G the filtering operator, here a Gaussian function:

G(x) = 6 π∆ 2 3/2 exp - 6x 2 ∆ 2 (3)
where ∆ is the filter size. The tilde operator denotes the density-weighted filtering defined by ρ φ = ρφ.

The sub-grid-scale transport term -∇ • ( ρ uc -ρ u c), the filtered laminar diffusion term ρD∇c and the filtered source term ωc , require closure models. When fully compressible solvers are considered, the closure of the energy balance equation must also be carefully addressed. For the sake of clarity, only the closure of the progress variable balance equation is discussed here. The full description of the F-TACLES model that include the closure of the filtered momentum and energy balance equations is given in [12].

The two possible approaches TFLES and F-TACLES to couple tabulated chemistry with flow field through the progress variable are now presented.

Introduction of tabulated chemistry in the Thickened Flame model

A correct flame propagation on a coarse grid is ensured by artificially thickening the flame front [26]. Eq. ( 1) is then replaced by:

∂ ρ c ∂t +∇•( ρ u c) = ∇•(ρDF∇ c)+ ωTab c [ c] F ( 4
)
where F is the thickening factor. As both molecular diffusion and chemical source terms have been multiplied and divided by F respectively, the flame front propagates at the same laminar flame speed S 0 l than the original flame but is thickened by a factor of F: c(x) = c(x/F).

Filtered Tabulated Chemistry model (F-TACLES)

The F-TACLES approach solves the filtered progress variable balance equation ( 1) where unclosed terms are modeled from filtered one-dimensional laminar premixed flames. The flame structure in the direction n normal to the flame front is assumed identical to the structure of a planar 1-D freely-propagating premixed laminar flame, computed with a detailed chemical mechanism involving N s species and filtered using a Gaussian filter (Eq. 3) of size ∆. The filtered source term for the filtered progress variable is then directly stored as a function of c and ∆ in a filtered chemical database:

ωc = ωFTab c [ c, ∆] (5)
As discussed in [12], the filtered molecular diffusion term in Eq. ( 1) is modeled by:

∇ • (ρD ∇c) = ∇ • α FT ab c [ c, ∆] ρD ∇ c (6)
where the correction factor α FT ab c [ c, ∆] is defined as:

α FT ab c [ c, ∆] = ρD ∂c * ∂x * ρD ∂ c * ∂x * -1 (7)
and estimated from filtered 1-D flame solutions and stored in a 2-D look-up table [ c, ∆]. The superscript * corresponds to extraction from one-dimensional solution. Finally, assuming that the flame is not wrinkled by unresolved flow motions, the unresolved convection term can be estimated as [12]:

-∇ • ρ uc -ρ u c = -ρ 0 S 0 l ∂c * ∂x * - ∂ c * ∂x * = Ω FT ab c [ c, ∆] (8)
where Ω FT ab c [ c, ∆] is also estimated from the filtered chemical database. The filtered progress variable balance equation for the model F-TACLES, neglecting subgrid-scale flame front wrinkling is therefore written:

∂ ρ c ∂t +∇•( ρ u c) = ∇• α FT ab c [ c, ∆]ρD ∇ c +Ω FT ab c [ c, ∆]+ ωc [ c, ∆] (9)
The filtered energy balance equation is closed similarly when fully compressible flow are considered [12]. Situations where unresolved flame front wrinkling should be considered are also discussed in reference [12].

Comparison of filtered and thickened flames

TFLES and F-TACLES models described in Section 2 are implemented into the AVBP code [27]. The coupling of the chemical table and the compressible solver is ensured following Tabulated Thermo-chemistry for Compressible flows (TTC) formalism [28]. Planar onedimensional premixed methane/air flames are computed using GRI 3.0 mechanism [29] for an equivalence ratio equal to 1.05 to generate chemical databases. The ability of TFLES and F-TACLES to exactly reproduce planar laminar flame propagation speed has been already shown elsewhere [START_REF] Poinsot | Theoretical and Numerical Combustion[END_REF]12], but the present parametric study focuses on numerical grid requirements.

Two thicknesses are now introduced to characterize laminar premixed flames: the reaction zone thickness δ 0 r , defined here as the full width at half maximum (FWHM) of the reaction rate ωc , and the thermal thickness δ 0 T , defined as the inverse of the maximum progress variable gradient (δ 0 r < δ 0 T ). Corresponding values for TFLES and F-TACLES flames are respectively ( δ r = Fδ 0 r , δ T = Fδ 0 T ) and (δ r , δ T ). Note that, when the filter size ∆ is larger than the flame front thickness, as assumed in F-TACLES, δ r ≈ δ T while for TFLES δ r < δ T . Then, the lower progress variable scales to be resolved in the simulation on a numerical grid of mesh size, ∆ x , are the reaction thicknesses δ r and δ r in the TFLES and F-TACLES approaches respectively.

Figure 1 displays the predicted F-TACLES (respectively TFLES) flame front propagation speed ratio S ∆ /S 0 l (respectively S F /S 0 l ) as a function of the length scale ratio δ r /∆ x (respectively δ r /∆ x ) for two values of the ratii δ r /δ 0 r and δ r /δ 0 r . The predicted filtered and thickened flame front propagation speeds, S ∆ and S F , remain very close to the reference laminar flame speed S 0 l as long as δ r /∆ x and δ r /∆ x ratii remain larger than 3. Below this limit, a departure between computed and reference flame speeds is observed because of the underresolution of reaction rate profiles. Therefore, both minimal filtered and thickened reaction zone resolutions on a grid of cell size ∆ x are given by: δmin r = δ min r = n∆ x (10) where n+1 is the minimal number of grid points needed to resolve the flame front. The minimal thickening factor for TFLES is then: F min = n∆ x /δ 0 r . For instance, present results show that flame simulations performed with AVBP using the third order TTGC numerical scheme [30] without artificial viscosity require n = 3.

Thickening and filtering operators as defined in Section 2 are applied to the detailed chemistry solution to determine δ r (∆) and δ r (F). The relation between F and ∆ such as δ r (∆) = δ r (F) is then deduced and displayed in Fig. 2. In practice, this relation may be used to determine ∆ min from δ r (∆ min ) = δ r (F min ). In Appendix A, an analytical relation between F and ∆ is derived modeling the progress variable reaction rate ωc by a Gaussian function and plotted in Fig. 2 (solid line).

Figure 3 shows the effect of the filtering and thickening operators on the flame structure when the resolved reaction zone thickness is preserved, i.e. when δ r = δ r . Two reaction zone thicknesses corresponding to F = 2.5 (∆ = 3δ 0 r ) and F = 5.7 (∆ = 7.9δ 0 r ) are considered. Fig. 3(a) confirms that both filtered and thickened progress variable reaction rates have the same thickness but the thermal thickness, defined as the inverse of the progress variable gradient, is larger for TFLES than with F-TACLES. Figure 4 compares F-TACLES (δ T ) and TFLES ( δ T ) thermal flame thicknesses plotted as a function of the thickening factor F = δ r /δ 0 r when reaction layer thicknesses are identical (δ r = δ r ). The thermal flame thickness is considerably more affected by TFLES than F-TACLES. As expected, δ T ≈ δ r , while δ T ≈ 2.5 δ r (i.e. δ T ≈ 2.5δ T ). This behavior is also visible in Figs. 3(b) and3(c), where temperature and velocity profiles are shown for both models. The response of a flame to a flow motion of size r is primarily controlled by the ratio r/δ 0 T . It has been shown in Ref. [20] that when the flame front is artificially broadened, then the flame response decreases. Therefore as δ T > δ T under similar mesh conditions, thickened flames will be always less affected by turbulence than filtered and original flames. In the next section, a 2-D unsteady laminar flame simulation is performed to verify these effects.

Unsteady 2-D Bunsen flame simulation

TFLES and F-TACLES formalisms are now tested against direct numerical simulations (DNS) to investigate the behavior of a bunsen laminar flame submitted to hydrodynamic excitation. Methane-air mixture at an equivalence ratio of 1.05 is issued from a 22 mm exit diameter nozzle burner. The inlet flow velocity is modulated harmonically according to v = 0.97 m.s -1 , v RMS = 0.19 m.s -1 and f = 62.5 Hz where v, v RMS and f are the mean inflow velocity, the inflow velocity RMS and the excitation frequency. This configuration is a 2-D simplification of the Schuller et al. [31] experimental set-up.

One-dimensional laminar methane-air flames are computed using the PREMIX software [32] and the GRI 3.0 chemical mechanism [29] to generate chemical look-up tables as a function of the progress variable defined from carbon monoxide and dioxide mass fractions

as c = (Y CO + Y CO 2 )/(Y eq CO + Y eq CO 2 )
where eq superscript denotes equilibrium state.

The objective of the simulations is to evaluate the ability of F-TACLES and TFLES formalisms to reproduce unsteady flame behavior when flow motions are fully resolved (no sub-grid scale flame front wrinkling) but not the internal flame structure. A direct numerical simulation (DNS) is first performed on a 1.1 million grid elements to fully resolve the laminar flame structure. The mesh has been designed to ensure that at least 15 nodes and 6 nodes are present in the thermal and reaction zones, respectively. The smallest unsteady flow structure is resolved over approximately 20 nodes. From the DNS solution, a reference "LES" field is then obtained by filtering the progress variable field with the operator defined in Eq. (3).

LES have been performed with both F-TACLES and TFLES formalisms on two coarse grids of 0.6 and 0.3 millions of elements, respectively. On each grid, the reaction zone thicknesses of both models are set identical. The meshing criterion (Eq. 10) is respected (four grid points across δ r ). Model parameters are summarized in Table 1. Results are analyzed after ten excitation cycles to ensure the independence of the solution from initial conditions.

Figures 5 and6 show consecutive sequences of the pulsed Bunsen flame front position during an excitation cycle. F-TACLES (dashed-dotted lines) and TFLES (dotted lines) are compared to the reference DNS solution plotted with solid lines. For moderate values of filter size relative to the reaction thickness i.e. ∆ = 3 δ 0 r , a very good agreement is observed between F-TACLES and the DNS (Fig. 5) while the corresponding thickened flame solution obtained with F = 2.5 slightly differs from the reference solution. The heat release rate integrated over the computational domain, displayed in Fig. 7 along an excitation cycle, is very well predicted by both approaches.

When the filter size increases up to ∆ = 7.9 δ 0 r , the flame dynamics remains fairly reproduced by F-TACLES simulation while the TFLES solution obtained with F = 5.7 fails to predict the flame patterns (Fig. 6). The heat release prediction is strongly impacted as shown in Fig. 8. F-TACLES predictions remain in good agreement even if a slight phase shift corresponding to a time lag of 0.1/ f is noticed. On the other hand, TFLES misses most of the unsteady behavior of the global reaction rate. As explained previously, for a given reaction layer thickness, the thermal thicknesses are differently affected in TFLES or F-TACLES approaches. As indicated in Table 1, the thickened thermal structure is 2.5 times larger than the reactive layer while the F-TACLES thermal layer remains of the same order than the filtered reaction zone.

It has been shown in Ref. [31] that a flame behaves as a low-pass filter. The thermal layer thickness broadening tends to decrease the filter cut-off frequency, increasing the flame response time. As discuss previously, the problem has three main length scales: r, the size of flow motions, δ T and δ T , the thermal resolved flame thicknesses when using F-TACLES and and TFLES respectively ( δ T > δ T ). When r >> δ T > δ T , the flame wrinkling will be fully resolved and both models are expected to provide similar results. When δ T > r > δ T . F-TACLES is still expected to provide good results while TFLES will damp a part of the flame wrinkling. As discussed in section 3, for similar meshes leading to similar reaction thickness, the thermal flame thickness is larger with TFLES than with F TACLES. Therefore more flow structures will be damped by TFLES than F-TACLES. F-TACLES is then able to explicitly resolve a larger part of the flame turbulence interactions than TFLES.

Conclusion

The ability of TFLES and F-TACLES approaches to reproduce flame structure and dynamics has been investigated. Relevant comparison between both models are performed for identical filtered and thickened reaction layer thicknesses ( δ r = δ r ) that were found to be the thinnest zone to be resolved on a given numerical grid. The flame structure and consequently the flow field are more affected by the thickening operation (TFLES) than by the filtering procedure (F-TACLES). In fact, the thermal thickness of the resolved flame front is of the same order of the reaction zone thickness when using F-TACLES (δ T ≈ δ r ) while it is larger, here by a factor of about 2.5, for TFLES. Simulations of a twodimensional Bunsen flame submitted to acoustic excitations show that the spreading of the thermal thickness strongly impacts the flame dynamics and might affect the prediction of combustion instabilities, a point to be investigated in the future. However, the F-TACLES formalism is expected to be more efficient than TFLES to describe flame dynamics and chemical flame structure in turbulent combustion regimes.

Table 1 :Figure 1 :

 11 Figure 1: Predicted flame speed ratii as a function of δ r /∆ x for δ r = δ r = 2.5δ 0 r (top) and δ r = δ r = 5.7δ 0 r (bottom). : F-TACLES S ∆ /S 0 l , : TFLES S F /S 0 l .

Figure 2 :

 2 Figure 2: Relation between ∆/δ 0 r and F when δ r = δ r . Symbols: Computations. Solid line: Analytical model (Appendix A).

Figure 4 :

 4 Figure 4: Evolution of the resolved thermal thickness as a function of the thickening factor F, : δ T , •: δ T

Figure 7 :

 7 Figure 7: Heat release rate integrated over the computational domain as a function of time during a cycle. δ r = δ r = 2.5δ 0 r ; -DNS; • -• F-TACLES; ---TFLES.

Figure 3 :Figure 5 :

 35 Figure 3: Effect of filtering and thickening operators on the flame structure. Solid lines: F-TACLES. Dashed lines: TFLES. Without symbol: δ r = δ r = 2.5δ 0 r . Symbols: δ r = δ r = 5.7δ 0 r . a) Progress variable source term, b) Temperature, c) Velocity

Figure 6 :

 6 Figure 6: Iso-line c = 0.5 during a cycle for δ r = δ r = 5.7δ 0 r ; -DNS filtered with ∆ = 7.9δ 0 r ; • -• F-TACLES with ∆ = 7.9δ 0 r ; ---TFLES with F = 5.7.

Figure 8 :

 8 Figure 8: Heat release rate integrated over the computational domain as a function of time during a cycle. δ r = δ r = 5.7δ 0 r ; -DNS; • -• F-TACLES; ---TFLES.

Acknowledgments

This work was supported by the ANR-07-CIS7-008-04 Grant of the French Ministry of Research and accessed to the high performance computing resources of IDRIS and CINES under the allocation 2009-i2009020164 made by GENCI (Grand Equipement National de Calcul Intensif). The authors also gratefully acknowledge Drs. T. Schuller and D. Durox for helpful discussions about their experiments.

Appendix A. Analytical relation between the filter size ∆ and the thickening factor F

The laminar reaction rate is modeled as a Gaussian function: ωc (x) = ρ 0 S l G(e r , x) with e r = 3 2 ln(2) 1/2 δ 0 r (A.1) where δ 0 r is the full width at half maximum of the progress variable reaction rate. In the expression e r is chosen so that the integral of the progress variable source term and its full width at half maximum are respectively equal to those obtained by detailed chemistry 1-D calculations done with PREMIX code. G(a, x) is the Gaussian function of size a defined as :

In TFLES approach, the thickened progress variable reaction rate reads:

The F-TACLES source term is written, for a filter size ∆:

remembering that filtering a Gaussian function of thickness e r with a Gaussian filter of size ∆ corresponds to a Gaussian function of size (e 2 r + ∆ 2 ) 1/2 . According to Eqs. (A.3) and (A.4), the condition δ r = δ r gives :

List of Figures

1 Predicted flame speed ratii as a function of δ r /∆ x for δ r = δ r = 2.5δ 0 r (top) and δ r = δ r = 5.7δ 0 r (bottom). : F-TACLES S ∆ /S 0 l , : TFLES S F /S 0 l . . . 9 2

Relation between ∆/δ 0 r and F when δ r = δ r . Symbols: Computations. Solid line: Analytical model (Appendix A). . . . . 9 4

Evolution of the resolved thermal thickness as a function of the thickening factor F, : δ T , •: δ T . . . . . . . . . . . 9 7

Heat release rate integrated over the computational domain as a function of time during a cycle. δ r = δ r = 2.5δ 0 r ; -DNS; • -• F-TACLES; ---TFLES. .