Carsten Schilde
email: c.schilde@tu-bs.de

Hendrik Nolte

Christine Arlt

Arno Kwade

Effect of fluid-particle-interactions on dispersing nano-particles in epoxy resins using stirred-media-mills and three-roll-mills

Effect of fluid-particle-interactions on dispersing nano-particles in epoxy resins using stirred-media-mills

and three-roll-mills

as single primary particles but rather as particle collectives consisting out of several primary particles.

For the application in reinforced composites the particles must be suspended in epoxy resin as separately dispersed primary particles or in a certain aggregate size. Generally, the influencing parameters to break up the aggregates in a dispersion process can be divided into the stress mechanism, the intensity and the frequency of the dispersing machine itself, the properties of the dispersed particles (e.g. the particleparticle-interactions) the properties of the homogenous phase and the particle-resin-interactions. Besides the effect of the chosen dispersing machine the optimization of the dispersing process was investigated by applying modified particle surfaces and varying the fluid properties. The results show that the surface properties of the particles must fit to the epoxy resin properties and the attractive forces between the primary particles must be reduced or the stabilization improved, respectively. An indication for an improved stabilization and adjustment of the particles surface properties to the fluid properties can be obtained by measurements of the contact angle and the rheological properties. Generally, an increase of viscosity and mass fraction of the product leads to a higher energetic efficiency of the dispersion process in the stirred media mill and three-roll-mill. KEYWORDS: Nano composites, nano particles, mechanical properties, dynamic light scattering, dispersing, formulation

INTRODUCTION:

The idea of producing nano-sized particles is to reach certain properties, which are principally determined by the high ratio of surface to volume. This ratio strongly influences characteristics like attractive forces between particles, their sizes themselves or the transparency of the solution. By embedding nanoparticles into the epoxy matrix, the mechanical properties of the composites as well as the electrical and thermal characteristics can significantly be improved. In many cases nano-sized particles are produced by pyrolysis or in precipitation processes, often followed by drying processes. In most cases nano-sized particles are produced not as single primary particles but rather as particle collectives consisting out of several primary particles, like agglomerates or more compact aggregates.

For most applications, the particles must be available in a liquid as separately dispersed primary particles or in a certain agglomerate size. Generally, the influencing parameters of producing nanoscaled dispersions can be divided into three interacting groups (see Figure 1): The stress mechanism, the stress intensity and the stress frequency as well as the specific energy supply are appointed by the dispersing machine itself. The particle properties on the other hand specify the resistance against fragmentation of aggregates in a dispersing process. Basically this contains the material properties, the surface modification, the porosity, the structure, the particle size distribution and the particle-particleinteractions. Especially for nano-particulate agglomerates and aggregates the particle-particle interactions are very strong and can be described according to the Derjaguin-Landau-Verwey-Overbeektheory [START_REF] Wang | Dispersion investigation of TiO2 nanoparticles coated by pulsed RF plasma polymer[END_REF]. Additionally the dispersion results are strongly influenced by the formulation. This includes the rheological properties of the liquid phase, the ion concentration, the additives and especially the wetting behaviour.

For the deagglomeration of silicia and alumina nano-particles the selection of the liquid phase and the variation of their rheological properties (i.e. temperature, wetting behaviour) as well as a modification of the particle surface by functional groups can optimize a dispersion process. The modification of the formulation, e. g. the particle-particle-interactions and the stability can strongly influence the dispersing result and is independent on the stress mechanism of the dispersing machine.

EXPERIMENTAL SET-UP

MATERIAL AND ANALYTICAL METHODS:

In the experiments a pyrogenic, hydrophilic alumina Alu C (Aeroxide® Alu C, evonik), a modified hydrophobic alumina (Aeroxide® Alu C 805, evonik), a hydrophilic silicon oxide Aerosil 200 V and a hydrophobic silicon oxide Aerosil R 104 were dispersed in different fluids at various fluid conditions. Alumina Alu C 805 represents a modified alumina Alu C with octylsilane-groups on the particle surface.

Aerosil R 104 represents a modified Aerosil 200 V with octamethylcyclotetrasiloxane-groups on the particle surface. The particles have a median primary particle size of 13 nm (alumina) and 12 nm (silica). Figure 2 shows a schematic representation of the functional groups of the surface modifications of the different dispersed particles.

As liquid phase fluids of high viscosity like glycerol and two different epoxy resins were used.

Glycerol shows a polar behaviour, whereas the epoxy resins show a more non-polar behaviour.

Additionally, hardening of the deagglomerated samples for additional scanning electron microscopy (SEM) and measuring the mechanical properties of the enhanced composites, is possible. In order to avoid hardening through high temperatures and long residence times in the dispersing machines, resins of two components like the epoxy resin Araldite LY 556 (CIBA-Geigy) and L135i (Lange + Ritter) were used. As indication of the dispersing result the particle size distributions at different dispersing times were measured via dynamic light scattering (Nanophox, Sympatec GmbH) or rotating disc centrifuge (CPS, cps instruments). In order to determine the particle size distribution correctly, a well stabilized sample is necessary Using dynamic light scattering and the disc centrifuge, the different particle-fluidsystems must be prepared using different solvents, i. e. butanon, butanol or destilled water and different stabilizers (Byk Additives and Instruments). In case of dispersing the unmodified, hydrophilic alumina in glycerol, the particle size was measured in distilled water without using any stabilizer. . For reproducible results each measurement of the particle size was averaged from four single measurements.

However, each single measurement includes 100.000 scans and a standard deviation of the median particle size of approximately 3 nanometers could be determined.. The wetting of the particles in the different fluids was detected via sorption measurements in the tensiometer K11 (Krüss GmbH). The contact angle between the particles and the fluid was calculated using the Washburn equation:

         2 cos 2 r t l (1)
The compression of the bulk commodity has to be well defined to compare the results among each other. Therefore, the nano-sized powders were sieved and compacted in a stamp volumeter. The capillary radius was appointed by a calibration run with an ideal wetting liquid (n-heptane). Afterwards the powders were analyzed with the actual carrier fluid. For reproducible results each measurement of the contact angle was averaged from five individual measurements.

Additionally the apparent viscosities of the suspensions at low shear rates were measured using a rotating viscosimeter (CVO 120 HR, Bohlin Instruments), which are a good measure for the particleparticle interactions in the suspension.

In order to characterise the efficiency of the dispersion SEM analyses were carried out. For this purpose the samples consisting of epoxy resin LY 556 were cured by adding 50% of an anhydride hardener HY 917 related to the mass of the pure epoxy resin.

DISPERSING SYSTEMS:

The efficiency of the dispersing process can be characterized plotting the median particle size versus the specific energy. The specific energy is defined as the quotient of the energy input and the mass of the solid product particles [START_REF] Mende | Mechanical production and stabilization of submicron particles in stirred media mills[END_REF][START_REF] Schilde | Dispersing and Grinding of Alumina Nano Particles by Different Stress Mechanisms, ceramic forum internatinal[END_REF]:

           solid solid m m dt T T n π m dt P P E ) ( 2 ) ( 0 0 (2)
For the dispersing experiments two dispersing machines with different stress mechanisms, a stirred media mills (Dispermat® SL, VMA Getzmann GmbH) and a three-roll-mill Exakt ES 80 (Exakt Vertriebs GmbH) were used. Experiments for dispersing alumina and silica were carried out using the stirred media mill Dispermat® SL (see Figure 3) with a grinding chamber volume of 0.125 l and a separation gap wide of 300 µm at operational temperatures up to 85°C. The small closed grinding chamber, the small dead volume within the pipes as well as the integrated temperature control system makes the mill suitable for problematical and irritant fluids like epoxy resins. An integrated air pressure system helps discharging the grinding chamber.

The stress intensity in a stirred media mill is proportional to the kinetic energy of the grinding media and inversely proportional to the stressed particle mass [START_REF] Kwade | A stressing model for the description and optimization of grinding processes[END_REF]. Consequently, i in case of stressing single particles at constant operating parameters of the mill the stress intensity is inversely proportional to the third power of the particle size. Besides the particle size, the kinetic energy of the grinding media E kin , strongly influences the stress intensity in stirred media mills and is proportional to the grinding media diameter, d GM , the tip speed, v t , and the grinding media density, ρ GM (see Eq. 3). Due to frictional losses, energy dissipation, deformation energy and resistance against fluid displacement during the approach of 2 grinding media the stress intensity is not equal but proportional to the kinetic energy of the grinding media:

3 GM 2 t 3 GM 3 kin kin ρ v d E E x x m SI      (3)
Three-roll-mills are frequently used to process high concentrated suspensions. Here the experiments were performed in a three-roll-mill Exakt ES 80 (see Figure 4): The gap width between the first and second rolls is set to 5 µm, while the second gap is set to 10 µm. During the dispersing process the adjacent rollers rotate at different speeds and opposite directions. The dispersing effect results from the high shear stresses generated in the gap between the rolls. The stress intensity can be increased by reducing the width of the gap between the rollers and increasing the rotational speed of the rollers and/or the mass concentration of the suspension. Assuming the shear stress acting on particles is only caused by the fluid, the stress intensity is independent on the particle size. This is described by the model of Rumpf and Raasch for high viscous laminar flow as shown in equation 4 [START_REF] Rumpf | Desagglomeration in Strömungen[END_REF][START_REF] Raasch | Die Bewegung und Beanspruchung kugelförmiger und zylindrischer Teilchen in zähen Scherströmungen[END_REF][START_REF] Raasch | Beanspruchung und Verhalten suspendierter Feststoffteilchen in Scherströmungen hoher Zähigkeit[END_REF]. For high viscous laminar flow and spherical agglomerates or aggregates, the stress intensity SI is proportional to the shear stress τ and the product of the dynamic viscosity η and the shear rate γ:

(sphere) 2,5 SI η γ τ     (4)
However, additional stresses can develop due to interactions between the particles by applying shear stress between the surfaces of the rollers. Thus, the total stress intensity might increase with decreasing particle diameter [START_REF] Rumpf | Desagglomeration in Strömungen[END_REF][START_REF] Raasch | Die Bewegung und Beanspruchung kugelförmiger und zylindrischer Teilchen in zähen Scherströmungen[END_REF][START_REF] Raasch | Beanspruchung und Verhalten suspendierter Feststoffteilchen in Scherströmungen hoher Zähigkeit[END_REF]. Table 1 shows the operating parameters of the three-roll-mill and the stirred media mill. Besides the mechanism, intensity and frequency of stresses, the particle properties (especially the surface modification) as well as the rheological properties of the liquid phase (like the viscosity as function of the temperature and the wetting behaviour) influence the dispersion result. Especially, dispersing at high solid concentrations in the three-roll-mill strongly depends on the fluid viscosity and the wetting behaviour. At high solid concentrations and viscosities, wall gliding of the suspension can occur at the wall of the rolls depending on the wetting behaviour between the rollers and the liquid phase. In order to verify the influence of the viscosity and wetting behaviour of different fluids, the phenomenon of wall gliding must be eliminated.

Low temperatures and, thus, high viscosities deliver a product of low particle fineness as seen in Figure 5. In this case wall gliding occurred for dispersing alumina in the epoxy resin L135i using a three-roll-mill at high shear rates. The occurrence of wall gliding can be analysed by comparing theoretical values with experimental data: The theoretical shear rate can be calculated with help of the gap size and the peripheral velocity. In case of wall gliding (see Figure 5), the effective shear stress is much lower as the calculated shear rate.

Because of the wall gliding effect the product fineness increases with increasing temperature or decreasing viscosity, respectively. Normally an increase of the temperature and, therewith, a decrease of the viscosity leads to a decrease of the stress intensity acting on the particle surface [START_REF] Rumpf | Desagglomeration in Strömungen[END_REF][START_REF] Raasch | Die Bewegung und Beanspruchung kugelförmiger und zylindrischer Teilchen in zähen Scherströmungen[END_REF][START_REF] Raasch | Beanspruchung und Verhalten suspendierter Feststoffteilchen in Scherströmungen hoher Zähigkeit[END_REF] and, thus, to coarser particle sizes. Additionally, the dispersing process gets energetically inefficient.

In order to determine the effect of the wetting behaviour of different fluids and of the attractive and repulsive forces between the primary particles in different fluids on the dispersion process, the viscosity of the fluids should be equal to eliminate viscosity effects. Therefore, the viscosities of the fluids at different temperatures must be measured. L135i (60 °C) shows angels of 90° or higher. Besides the mass fraction of the suspension and viscosity of the solvent, the dispersion efficiency strongly depends on the wetting behaviour of the particles in the different solvents. Figure 8 shows the contact angle of alumina the unmodified alumina in the different solvents.

Figure 9 shows the effect of different surface modifications (unmodified alumina and modified alumina, unmodified silica and modified silica) on the dispersing result achieved in the three-roll-mill using epoxy resin LY 556. Different surface modifications have an influence on the particle-particle interactions and the wetting behaviour. Figure 2 shows a schematic representation of the functional groups of the surface modifications of the different particles. In case of the modified silica, the octamethylcyclotetrasiloxane-groups on the particle surface represent a good adjustment to the epoxy resin and the attractive forces between the primary particles are relatively low. Compared to the modified silica the attractive interaction forces of the unmodified silica are even lower. However, the affinity of the unmodified silica to the epoxy resins is comparatively bad. Comparing the two hydrophilic versions the unmodified silica and unmodified alumina, lower repulsive electrostatically interaction forces can be found for the unmodified silica. This can be a consequence of the hydrophobic siloxan groups (see Figure 2) on the particle surface of the unmodified silica. Additionally, zeta potential measurements in distilled water were performed using a DT 1200 (Dispersing Technology). The zeta potential of the unmodified silica can be found within the range of -18 and -30 whereas the zeta potential of the unmodified alumina shows values around 50. Concluding, the dispersion process of the modified silica in epoxy resin LY 556 is more efficient and leads to smaller median particle sizes as the unmodified silica.

In case of the unmodified and modified alumina (Figure 9), the adjustment of the surface modification to the epoxy resin is not optimal. The major difference in this case are the lower attractive forces of the unmodified alumina through positive charges on the particle surface compared to high attractive van der waals forces of the modified alumina through the long hydrophobic octylsilane-groups (see Figure 2).

For this reason the dispersion process of the unmodified alumina in epoxy resin LY 556 is more efficient and leads to smaller median particle sizes than the one of the modified alumina. These considerations can be confirmed by viscosity measurements of the different particle-fluid-systems at low shear rates (see Figure 10). The values of the dynamic viscosity at low shear rates give an indication for the particle-particle interactions. In case of hydrophobic alumina (modified) and silica (modified) the viscosity at low shear rates is much higher as the viscosity of the referring hydrophilic particles. The explanation is similar as the one given above. Attractive interactions between hydrophobic particles cause high viscosities at low shear rates. Strong repulsive electro statically interactions can especially be remarked in case of the unmodified alumina. Comparing the modified silica with the unmodified, the modified silica again delivers higher particle-particle interactions. On the other hand, the octamethylcyclotetrasiloxane-groups represent a better adjustment to the epoxy matrix. In this case the adjustment to the matrix has more influence on the dispersion and is more efficient [9-11].

For comparison Figure 11 shows the dispersing of aluminia and silica with a hydrophilic as well as hydrophobic surface in epoxy resin LY 556 using a stirred media mill. Here the dependency of the product fineness on the specific energy is shown. Generally, the dispersion of the two different alumina particles shows similar behaviour for the stirred media mill (Figure 11) and the three-roll mill (Figure 9).

As seen before the adjustment of the surface modifications of both alumina particles to epoxy resin is not optimal. The difference are the lower attractive van der waals interaction forces and high repulsive electrostatical interaction forces of the unmodified alumina with positive charges on the particle surface compared to higher attractive van der waals forces of the modified alumina with the long hydrophobic octylsilane-groups. For this reason the dispersion process of the unmodified alumina in epoxy resin LY 556 is more efficient and leads to smaller median particle sizes than the one with the modified alumina.

Anyway, Figure 11 shows higher fluctuations of the median particle size in the stirred media mill during the dispersion process as found for the three-roll-mill. The basic problem of operating the stirred media mill Dispermat Sl is the exact tempering of suspension during the experiments. This leads to fluctuating temperature values in the range of 5°C and, thus, to a fluctuating median particle size. Especially for dispersing the modified alumina, the tempering of the suspension is difficult and leads to different particle sizes especially for small specific energies. Generally, the dispersion processes in a stirred media mill and the three-roll-mill differ due to the different stress mechanism of the dispersing machine.

Moreover, the particles -grinding media -interactions change the dispersing behaviour in a stirred media mill. Especially dispersing silica particles in the stirred media mill results in a differing dispersion process. Figure 11 shows an unexpected dispersing behaviour for the unmodified silica. Dispersing the unmodified silica in the stirred media mill gets more efficient at high specific energies compared to dispersing the modified silica.

Figure 12 shows the SEM pictures of the hardened samples of the unmodified silica (a), modified silica (b), unmodified alumina (c) and modified alumina (d) in epoxy resin LY 556.

The hydrophilic particle-epoxy resin LY 556-systems of the unmodified alumina (c) and silica (a) after curing leads to a homogenous distribution of the particles within the matrix. The hydrophobic particleepoxy resin-systems, especially the one containing the modified alumina (d), are reagglomerated. For alumina the unmodified alumina (c) there are many small particles below 50 nm and a few large aggregates. The SEM pictures prove the efficient dispersing of the unmodified alumina obtained by the dispersing experiments. Additionally, the determined median particle sizes for silica are similar to those ones measured via dynamic light scattering.

For an increase of the Young's modulus of composites a high solid content of nanparticles in the matrix are advantageous. Thus, test specimen for the measurement of the mechanical properties were produced using suspensions, which were manufactured using the three-roll-mill with high solids concentrations and with small small particle sizes. Due to the strong agglomeration of the modified alumina during the curing process only mechanical tests of the unmodified alumina are carried out.

CONCLUSION:

Fibre reinforced composites are indispensable in the field of modern lightweight structures, such as used in aerospace, automotive industry or in wind power plants. Those materials provide high weight savings and significantly increase the efficiency of a structure. By embedding nanoparticles into the epoxy matrix, the mechanical properties as well as the electrical and thermal characteristics can significantly be improved. Besides the characteristics of the dispersing device, the formulation of the suspension significantly influences the energy efficiency and the obtained product fineness of the dispersing process. Besides the mechanism, the intensity and the frequency of the stress transferred by the dispersing machines, the surface modification, the fluid properties and at least the stabilization method have a great influence on the energy efficiency and the obtained product fineness of the dispersing process. Generally, an increase of the viscosity and mass fraction of the product leads to an increase of the energetically efficiency of the dispersion process in the stirred media mill and three-rollmill. In the three-roll-mill the effect of wall gliding of the suspension between the rolls leads to lower stress intensities and thus to coarser product quality. An optimization of the dispersing process can be achieved if the surface properties are adapted to the fluid properties. Additionally, the attractive forces between the primary particles have to be reduced or stabilized. The level of optimization can be indicated by measuring the contact angle and the rheological properties. For dispersing the same particle-liquid system in a stirred media mill instead of in a three-roll-mill the dispersion result can differ strongly. This is related to the different stress mechanisms and additional particle-grinding media interactions.

With increasing alumina solid content in the cured alumina epoxy resin composites the Young's modulus and the fracture strength rise, while the fracture strain decreases. At high solid contents the ceramic nanoparticles lead to an embrittlement of the elastic resin matrix.

ACKNOWLEDGEMENT:

This study was carried out in cooperation with company Exakt GmbH.

Figure 6

 6 shows the dependency of the fluid viscosity at a shear rate of 1000 1/s on the temperature for glycerol, epoxy resin L135i and epoxy resin LY 556. The viscosity of all fluids shows Newtonian flow behaviour. Using Figure 6 the operating temperature for further experiments in the three-roll-mill for the different fluids is determined. In the following the unmodified alumina was dispersed with the three-roll-mill in different fluids at different temperatures but approximately equal fluid viscosity of 0,72 Pas. As shown in Figure 7 the dispersion of the unmodified alumina in glycerol is more efficient than the dispersion in the epoxy resins. Dispersing the unmodified alumina in the epoxy resin L135i requires more specific energy than dispersing it in epoxy resin LY 556 and glycerol to achieve similar product fineness. In accordance with the dispersing experiments the results of the contact angle measurements confirm the tendency of the dispersing experiments: The contact angle of the unmodified alumina in glycerol (60 °C) was determined to approximately 40°, epoxy resin LY 556 (85 °C) was in the range of 80° and epoxy resin

Figure 13

 13 Figure13shows the percentaged change in the Young's modulus, the fracture strength and the fracture

Figure 1 :Figure 2 :Figure 3 :Figure 4 :

 1234 Figure 1: Influencing factors on the dispersing result

Figure 5 :Figure 6 :Figure 7 :Figure 8 :

 5678 Figure 5: Product fineness as function of the specific energy for dispersing nano-sized alumina in

Table 1 :

 1 Constant operating parameters of the three-roll-mill and the stirred media mill

	three-roll-mill
	gap 1	g 1	5 µm
	gap 2	g 2	10 µm
	relative peripheral velocity	v rel	0,186; 0,372 m/s
	mass fraction	c m	0,25
	stirred media mill
	diameter of the grinding media	d GM	1,1 mm (ZY 11)
	peripheral velocity	v	5,2 m/s
	mass fraction	c m	0,09

3. RESULTS AND DISCUSSION-EFFECT OF FLUID-PARTICLE-INTERACTIONS:

  The authors gratefully acknowledge for the useful discussions held with Prof. Dr. Garnweitner from the Institute for Particle Sun, Y., Zhang, Z., Wong, C.P. Study on mono-dispersed nano-size silica by surface modification for underfill applications, J. Colloid Interf. Sci., 2005, 292, 436-444 [10] Kang, S. et al. Preparation and characterization of epoxy composites filled with functionalized nanosilica particles obtained via sol-gel process,Polymer, 2001, 42, 879-887 [11] Jang, H. S., Park, D.-W., Shim, S. E., Measurement of dispersion stability of surface-modified nanosized carbon black in various liquids, J. Nanosci. Nanotechnol., 2007, 3827-3829 [12] Moloney, A. C., Kausch, H. H., Kaiser, T., Beer, H. R. Parameters determining the strength and toughness of particulate filled epoxide resins, Journal of Material Science, 1987, 22, 381-393 [13] Ng, C. B., Schadler, L. S., Siegel, R. W. Synthesis and mechanical properties of TiO2-epoxy nanocomposites, Nanostructured Materials, 1999, 12, 507-510

	[9]
	Technology
	SYMBOL INDEX
	c m	[-]	mass fraction
	d GM [mm]	grinding media diameter
	E m	[kJ/kg]	mass specific energy
	g	[µm]	gap range between the rollers
	L	[m]	wetting level
	m solid [kg]	mass of the solid
	P	[W]	power input
	P0	[W]	ideal power
	r	[µm]	capillary radius
	t	[s]	dispersing time
	T	[Nm]	torque
	T 0	[Nm]	ideal torque
	T	[°C1]	temperature
	v	[m/s]	tip speed
		[1/s]	shear rate
	η	[Pa s]	dynamic viscosity
	σ	[N/m]	surface tension
	Θ	[°]	contact angle


Tobias Ziegler

Achim Neubrand
email: achim.neubrand@iwm.fraunhofer.de

Romana Piat

Multiscale Homogenization Models for the Elastic Behaviour of Metal/Ceramic Composites with Lamellar Domains

Keywords: Metal-matrix composites (MMCs) (A), Mechanical properties (B), Finite element analysis (FEA) (C), Modelling (C)

This study predicts the elastic properties of an innovative metal-ceramic composite with statistically oriented domains of parallel ceramic platelets embedded in a eutectic Al-Si-alloy. For this purpose, a two-step homogenisation procedure was employed by finite element-and micromechanical modelling. In a first step, the microstructure of the specimen was divided in domains with the same orientation of lamellae and the elastic properties of single domains were calculated while a precise representation of the shape of the lamellae was attempted. In a second step the elastic characteristics of a large specimen consisting of many domains were computed both by finite element and micromechanical modelling. The experimental Young's modulus of such polydomain specimens was determined by an acoustic resonance method and was lower than predicted. The differences can be explained by microcracks caused by large residual microstresses produced in these materials when they are cooled from the manufacturing temperature.

Introduction

Most crystalline materials exhibit anisotropy at the microscopic scale, and homogenization is necessary to predict the response of large polycrystalline specimens from single crystal data. A common homogenization scheme is the meanfield-approach (MFA) by Hill [START_REF] Hill | Elastic properties of reinforced solids: some theoretical principles[END_REF] where the stress and strain fields are averaged over a representative volume element (RVE). Alternatively, a perturbation method can be used to compute the elastic constants of different arrays of rigid spheres in a homogeneous isotropic elastic medium [START_REF] Nunan | Effective elasticity tensor of a periodic composite[END_REF]. These methods were compared by Hori and Nemat-Nasser [START_REF] Hori | On two micromechanics theories for determining micromacro relations in heterogeneous solids[END_REF]. They found the first order terms in the expanded stress and strain fields to be the same and used this correspondence to formulate a hybrid approach. Counts and co-authors [START_REF] Counts | A comparison of polycrystalline elastic properties computed by analytic homogenization schemes and FEM[END_REF] investigated a polycrystalline MgLi alloy and an analytical self-consistency approach as well as finite element modelling (FEM) were employed to obtain the homogenized elastic properties. They favoured the analytical model for its simplicity. There are different micromechanical approximations which can be used for multiphase microstructure modelling. An overview of such methods can be found in monographs by Nemat-Nasser and Hori [START_REF] Nemat-Nasser | Micromechanics: overall properties of heterogeneous materials[END_REF] and Kachanov et al. [START_REF] Kachanov | Handbook of elasticity solutions[END_REF].

A comparison of FEM unit cell models with a Mori-Tanaka MFA was made by Duschlbauer et al. [START_REF] Duschlbauer | Computational simulation of Composites Reinforced by Planar Random Fibers: Homogenization and Localization by Unit Cell and Mean Field Approaches[END_REF] for fibre reinforced composites with a planar random orientation of the fibres. They found an excellent agreement between both methods and concluded that the Mori-Tanaka estimates are fully adequate for predicting the linear responses.

For materials with hierarchical microstructure the multistep homogenization procedure can be used. Such procedures are derived from the Benveniste method [START_REF] Benveniste | A new approach to the application of Mori-Tanaka's theory in composite materials[END_REF][START_REF] Benveniste | Stress fields in composites with coated inclusions[END_REF]. Examples of the application of such methods for different materials can be found for functional graded materials in [START_REF] Genin | Micromechanics and structural response of functionally graded, particulate-matrix, fiber-reinforced composites[END_REF], for polycrystals in [START_REF] Iwakuma | An estimate of average elastic moduli of composites and polycrystals[END_REF] and for modelling of the non-elastic repose of the elasto-plastic materials with hierarchical microstructure in [START_REF] Doghri | Micromechanical modelling and computation of elasto-plastic materials reinforced with distributed-orientation fibers[END_REF][START_REF] Friebel | General mean-field homogenization schemes for viscoelastic composites containing multiple phases of coated inclusions[END_REF]. The general multi-step homogenization method is proposed for shale-type rock materials in [START_REF] Ulm | Experimental microporomechanics[END_REF], for bone in [START_REF] Ch | Microelasticity of bone[END_REF] and for rock-like composites in [START_REF] Huynh | Application of results on Eshelby tensor to the determination of effective poroelastic properties of anisotropic rocks-like composites[END_REF][START_REF] Huynh | Effective poroelastic properties of transversely isotropic rock-like composites with arbitrarily oriented ellipsoidal inclusions[END_REF].

In this paper we will use such a homogenization procedure for modelling of the microstructure by semi-analytical and FE-methods.

The material under investigation in the present paper is a new type of metal/ceramic composite with a hierarchical structure with randomly orientated regions in one cross section in which all ceramic lamellae are parallel (hitherto called domains, cf. Figure 1 a) ). It is produced by infiltrating porous alumina preforms with molten AlSi12 via squeeze-casting. The lamellar microstructure of these preforms is produced by directional solidification along the Z-direction during freezing of a ceramic slurry. The size and internal structure of the domains can be controlled via the freeze casting parameters [START_REF] Waschkies | Tailor-Made Ceramic Preforms by Freeze-Casting[END_REF][START_REF] Fukasawa | Synthesis of porous ceramics with complex pore structure by freeze-dry processing[END_REF]. These lamellar preforms have a high permeability and are thus specifically suited as preform material for the production of metal-ceramic composites via squeeze-casting or even die-casting [START_REF] Ghomashchi | Squeeze casting: an overview[END_REF][START_REF] Clyne | The squeeze infiltration process for fabrication of metal-matrix composites[END_REF] with reinforcement content in the range from 30 to 70 vol.% [START_REF] Deville | Freezing as a part to build complex composites[END_REF][START_REF] Deville | Ice-templated porous alumina structures[END_REF]. Due to their elevated ceramic contents such composites are expected to have a high stiffness and wear resistance and are thus of particular interest for applications in the automotive industry.

In a preceding paper, the anisotropic elastic behaviour of single domains in which all lamellae are parallel were investigated by FEM, a micromechanical model and experimentally [START_REF] Ziegler | Elastic constants of metal/ceramic composites with lamellar microstructures: Finite element modelling and ultrasonic experiments[END_REF]. In this paper finite element modelling is used to predict the elastic response of a real-world material sample consisting of many domains. The results are compared to micromechanical models established for ideal laminates and composites with ellipsoidal inclusions (Mori-Tanaka and inverse Mori-Tanaka) and to experimental results.

Microstructure characterisation and mechanical testing

The processing of the studied AlSi12/Al 2 O 3 composites was described in detail earlier [START_REF] Ziegler | Elastic constants of metal/ceramic composites with lamellar microstructures: Finite element modelling and ultrasonic experiments[END_REF][START_REF] Waschkies | Control of lamellae spacing during freeze casting of ceramics using double-sided cooling as a novel processing route[END_REF]. The studied materials were produced using melt infiltration of the ceramic preforms, which were produced by freeze-casting at -10°C. The preform architecture is induced through the geometry of the growing ice crystals, and results in a material microstructure with translation symmetry in the direction of ice crystal growth which can be adequately modelled in 2D [START_REF] Roy | Metal/Ceramic Composites from Freeze-Cast Ceramic Preforms: Domain Structure and Elastic Properties[END_REF]. The composites under consideration contained between 36 and 42 vol. % of ceramics, which is arranged in domains of parallel lamellae (the ceramic content of the composites was deduced from density measurements of the preform). The typical lamella spacing is about 50 µm (Figure 1 a). Rectangular poly-domain samples with dimensions 8x8x60mm³ were cut from cast blocks of 10x40x60 mm³ size. A schematic of such a poly-domain sample is shown in Figure 1 b). The edges of the cut poly-domain samples were always parallel to the edges of the original composite plate, i.e. one edge was always parallel to the freezing direction. Mass and dimensions of each sample were measured using a high precision laboratory balance and a micrometer scale, and their densities were determined there from.

The frequency of the first bending resonance was determined for bending in direction Z (axis convention in Figure 1 b). The elastic modulus of these specimens was determined from the weight, the dimensions of the specimen and the resonance frequency according to [27].

Modelling

The microstructural modelling of the studied materials was accomplished using two methods: FE-modelling of the real microstructure and micromechanical homogenisation. Both methods use a two step homogenization procedure:  The polydomain is divided in single domains and the material properties of polydomain (mesoscopic) sample are calculated by averaging over different oriented single domains.

 The anisotropic elastic material parameters of single domains with varying ceramic content are calculated and used as input for the polydomain model.

In the following we show this two-step procedure for FE and micromechanical modelling.

Finite Element (FE) modelling of the polydomain

For FE-modelling the micrographs of three polydomain specimens with an area of 60 x 40mm² each were used and divided into domains of conform lamella orientation (an example is shown in Figure 2a)). Each domain was then colour-coded with a colour representing the orientation of the lamellae  in steps of 15°. The resulting domain representation can be seen in Figure 2b). The different orientation distributions of the domains for three studied microstructures and their average (effective microstructure) in relation to a random domain distribution are presented in Figure 3. For all studied microstructures, some preferred orientations in 30-45° and 90° directions can be observed, but the preference is not very significant (max 8% deviation from random). Insignificance of the influence of such small microstructural variations for the effective material properties of polydomains will be demonstrated later by FEM and micromechanical modelling. The distribution of the lamellae orientation in this material map was used in the FEmodel and micromechanical model. The coloured image was imported into the software OOF2® [START_REF] Langer | OOF2: Object-Oriented Finite Element Analysis of Material Microstructures developed[END_REF] and a FE mesh was adjusted to it where each colour was treated as a different element group. The final mesh was then extruded in the third direction (Z-axis) using Patran® [29] and its elastic behaviour was modelled with the commercially available ABAQUS 6.8 FE-code [START_REF] Abaqus | [END_REF]. In order to determine the unknown components of the stiffness tensor C ijkl of the polydomain, displacements were applied to different sides/faces of the mesh while keeping the nodes on the other sides/faces fixed according to the associated mode of deformation. The components of the stiffness tensor were determined from the resulting stress on the displaced side for different single domains with varying ceramic contents. In the following we also use the Voigt notation of the stiffness tensors.

FE-modelling of the single domain

The anisotropic elastic constants for the single domains were taken from the previous modelling of single domain samples with different ceramic contents [START_REF] Ziegler | Elastic constants of metal/ceramic composites with lamellar microstructures: Finite element modelling and ultrasonic experiments[END_REF]. Real microstructure FE meshes were generated from microstructure images of entire surfaces of single domain specimens perpendicular to the freezing axis Z (= 1-axis in the single domain model) using the freeware program OOF2 [START_REF] Langer | OOF2: Object-Oriented Finite Element Analysis of Material Microstructures developed[END_REF]. The 2D-meshes were then extruded along the freezing direction (1-axis) in a subsequent step using Patran (see Figure 4). These meshes were used for obtaining all components of the stiffness tensor assuming an orthotropic symmetry. For alumina an isotropic Young's modulus of 400 GPa and a Poisson's ratio of 0.25 were taken as plausible values from the literature [31]. For the AlSi12 alloy the Young´s modulus was measured separately to 77.6 GPa, a Poisson's ratio of 0.3 was used for the AlSi12 (for more details we refer to [START_REF] Ziegler | Elastic constants of metal/ceramic composites with lamellar microstructures: Finite element modelling and ultrasonic experiments[END_REF]).

Micromechanical modelling of polydomains

Another possibility for calculation of the effective material properties of the studied materials is micromechanical homogenisation. In our studies we have assumed that the microstructure does not change along the freeze casting direction and the problem can be treated as a 2D-problem.

For micromechanical modelling the same mesoscopic distribution of the different domain orientations as in the FE-model was used (Figure 2 b). For calculation of the effective properties of the microstructure, which consists of differently oriented domains, the following homogenization procedures can be used:

Voigt approximation, which gives stiffer parameters than real ones:

   N f 1    C C ; (1)
Reuss approximation, which gives the lower bound for material parameters:

     N f 1 1 1 ) (    C C . ( 2
)
 f is the surface fraction of the domains with orientation  and  C is the stiffness of these domains in the global coordinate system.

For calculation of the effective properties of the polydomain, the stiffness matrix of the single domain  C is needed.

Micromechanical modelling of a single domain

For the micromechanical modelling, the Mori-Tanaka and inverse Mori-Tanaka models were applied in which the lamellae (inclusions) with irregular boundaries were approximated as 2D ellipsoids. More about this procedure can be found in [START_REF] Piat | Material modelling of the CVI-infiltrated C-felt. II. Statistical study of the microstructure, numerical analysis and experimental validation[END_REF].

The inverse Mori-Tanaka method gives better estimations than the Mori-Tanaka method for large volume fractions of inclusions (in our case lamellae) and in our case the material had about 58-64 vol. % of inclusions (metal alloy). Mori-Tanaka and inverse Mori-Tanaka correspond to upper and lower estimates for material properties of the domain. In our modelling, we have supposed that the inclusions can be approximated as elliptic cylinders. For this type of inclusions, the components of the Eshelby tensor c S can be presented as functions of the half-axes of the ellipses and

Poisson ratio of alumina [START_REF] Nemat-Nasser | Micromechanics: overall properties of heterogeneous materials[END_REF][START_REF] Kachanov | Handbook of elasticity solutions[END_REF].

Then, the Mori-Tanaka model gives the effective stiffness tensor of the domain according to the relation:

m MT c m m c dom f ) ( : ) ( A C C C C    ; ( 3
)
1 1 4 ) ( )] ( : ) ( : [      c m c c c m m MT f f C C C S I A .
Where 4 I is the symmetric fourth-order unit tensor:

) ( 2 1 4 jk il jl ik ijkl I       ;

Elastic properties of single domains

The properties of single domains with parallel lamellae orientation and different ceramic contents were required as input data for the homogenisation. The components of the stiffness tensor, calculated as described in chapters 3. (perpendicular to the lamellae). The calculated shearing constant C 12 is higher than the other two constants C 13 and C 23 which differ only slightly (around 2% for the FEM).

Elastic properties of polydomains

The elastic constants for large specimens consisting of many domains (as shown in The constant C ZZ along the freezing direction is the highest for all models.

Differences between the three models are small (up to 2%) for this direction. The constants C XX and C YY (perpendicular to the freezing direction) are close but not equal for all three models, and the difference between all models is small. The shearing constants C XY , C XZ , and C YZ differ only slightly with the constant for shearing in the XY-plane being the smallest. All elastic constants increase monotonously with the ceramic content for both models as expected. The measured Young's modulus in the X-direction of polydomain specimens are compared to the FE-modelling results in Figure 7. The experimental Young's modulus data show a large scatter and are generally lower than the FEM results. No increase with the ceramic content is observable in the experimental results. The upper and lower Hashin-Shtrikman (HS) bounds [START_REF] Hashin | A variational approach to the theory of the elastic behaviour of multiphase materials[END_REF] and the Reuss isostress bound [START_REF] Reuss | Berechnung der Fließgrenze von Mischkristallen auf Grund der Plastizitätsbedingung für Einkristalle[END_REF] are given for comparison. Most experimental data points are outside of the HSbounds and some are so low they fall under the Reuss isostress bound. The HSbounds are given in terms of the bulk and shear moduli as:

m m m m c c m lower G K f K K f K K 4 3 3 1 *      ; c c c c m m c upper G K f K K f K K 4 3 3 1 *      ; ) 4 3 ( 5 ) 2 ( 6 1 * m m c m m m m c c m lower G K G f G K G G f G G       ; ) 4 3 ( 5 ) 2 ( 6 1
* c c c c c c c m m c upper G K G f G K G G f G G       .
where f m denotes the volume fraction of the metal in the composite and f c denotes the volume fraction of the ceramic accordingly.

Discussion

Elastic constants of single domains

The FEM and micromechanical results for the single domains in Figure 5 show the highest stiffness along the freezing direction in which a perfectly continuous microstructure is assumed. For a more detailed discussion of the results for normal modes of deformation and a comparison to experimental data on single domain specimens we refer to [START_REF] Ulm | Experimental microporomechanics[END_REF]. The shearing constants C 12 are the highest because for this deformation the ceramic lamellae are sheared lengthwise. C 13 and C 23 are smaller and have similar values. All constants increase with rising ceramic content of the samples. As real microstructures (an example is depicted in Figure 4) were used for obtaining these results small deviations from a monotonous relationship occur.

Comparison of micromechanical and FE modelling for polydomains

The Finite element modelling requires more effort as compared to the analytical model, but allows a better insight into the stress and strain distributions within the microstructure. An example is shown in Figure 9. It can be seen that domains with lamellae oriented along the loading direction show high tensile stresses (solid white arrows) and that domains perpendicular to the loading direction show lower tensile or even compressive stresses (dotted white arrows). This stress distribution is advantageous, because experiments on small single domain samples showed that the strength is highest along this direction. The FEM model can be extended to include plasticity and account for the thermal residual stresses that result from the infiltration process.

Comparison of modelling results with experimental values

Figure 7 shows a discrepancy between experimental and modelling results for the elastic modulus. This could be due to the resonance method not taking the anisotropy of the material into account. Calculation of the resonance modes on the 3D FE-model showed only a minor influence (less than 1%) of the occurring anisotropy on the resonance frequency of the material. Hence, errors in experimental methodology can be excluded as a reason for the discrepancy. Nevertheless, the experimental values are typically considerably lower than predicted by the FE-and HS models. This indicates that the assumptions made in these models (dense material, no microcracks or debonded interfaces) do not apply to the material as produced. Indeed, radiographic images of most composite specimens show some residual porosity (areas which are not completely infiltrated with metal). Density measurements indicate that this porosity cannot be very large (less than 3% on average). Both, estimates using microstructural models and our previous experience with infiltrated MMC indicate that this porosity could only account for a moderate modulus reduction of less than 10 GPa as compared to a dense specimen. However, radiographic images of some specimens showed cracks perpendicular to the image plane that were large enough to be visible. As the metal-ceramic-interfaces of the lamellar specimens lie in this plane it is very likely that the cracks run along or very near this interface. We believe that the aforementioned relatively large cracks can best account for the very low modulus observed in some specimens and its large scatter. We believe that these cracks may have formed during cooling from the infiltration temperature at the interfaces between domains of different lamella orientation. Elastoplastic FE calculations show a very strong anisotropy of thermal contraction, which can differ by a factor of more than three for different lamella orientations. Elastoplastic modelling of polydomain specimens shows high thermal residual stresses at domain interfaces which could well account for the observed crack formation.

Conclusion

The elastic constants of a novel MMC microstructure consisting of domains with parallel ceramic lamellae could be predicted with multiscale FE and micromechanical models. The FE model was based on the real microstructure and represented the irregular shape of the ceramic lamellae on the microscale whereas the micromechanical model was based on the Mori-Tanaka approach and approximated the lamellae as 2D ellipsoids. The results achieved with the two approaches are similar and the largest differences in predicted elastic constants of a polydomain specimen were only a few GPa. Both models have their advantages. The micromechanical model is easily applicable for different microstructures and could be used as a step for more complicated calculations. The FE model was able to reveal large differences in the stress state of different domains under external tensile loading -the largest stresses occurred in domains with lamella orientation parallel to the applied stress.

Figure 1 :

 1 Figure 1: a) Polarised light microscopy micrograph of the microstructure and single domain

Figure 2 :

 2 Figure 2: a) Micrograph of a typical polydomain specimen divided into domains of parallel

Figure 3 :

 3 Figure 3: Deviation of the domain distribution from random for three studied microstructures

Figure 4 :

 4 Figure 4: a) Micrograph of a typical single domain specimen viewed from the top (cross

  f are stiffness tensors and volume fractions of ceramic and metal phases correspondly.The stiffness tensor C of the single domain with orientation  to global coordinates can be calculated by transformating dom C into global coordinate system. Inverse Mori-Tanaka model is exactly the same method as Mori-Tanaka, however, the mechanical properties of the matrix and inclusions (i.e. of alumina and aluminium) are permuted.

Figure 5 :

 5 Figure 5: Elastic constants for single domain samples with different ceramic contents. C 11 , C 22

  2 and 3.4 using FE and micromechanical modelling are shown in Figure 5. The elastic constants C 11 along the freezing direction are predicted to be the highest, C 22 is slightly lower than C 11 (around 7% for the FEM), and the lowest values are expected for C 33

Figure 2 )

 2 Figure2) were calculated for different ceramic contents using FE-and

Figure 6 :

 6 Figure 6: Elastic constants of the mesoscopic models for polydomain specimens. The elastic

Figure 7 :

 7 Figure 7: Experimental Young's modulus and FE-modelling data for mesoscopic material are

5 . 3 . 8 Figure 8 :

 5388 Figure 8: Elastic constants of the virtual microstructures with random and average domain

Figure 9 :

 9 Figure 9: left: Detail of the  22 (in Pa) stress distribution in the microstructure for a vertical

Acknowledgements

We thank A. Wanner and S. Roy of the Institut für Werkstoffkunde, University of Karlsruhe for providing micrographs of single domain specimens, and M.J. Hoffmann, R. Oberacker, and T. Waschkies of the Institute for Ceramics in Mechanical


O Onokpe

P S Cumber
email: p.s.cumber@hw.ac.uk

MODELLING LIFTED HYDROGEN JET FIRES USING THE BOUNDARY LAYER EQUATIONS

The focus of this paper is the simulation of lifted hydrogen jet fires. The computational modelling of rim stabilised fires is mature. This is not the case for lifted jets where computational studies have focussed on understanding the mechanism for the location of the combusting flame base and good agreement between predicted and measured flow properties is not universal.

The simulation of hydrogen jet fires is an interesting and current area of research due to the sensitivity to global warming and the potential to address this problem with "The Hydrogen Economy" concept. To utilize hydrogen successfully requires the development of robust and accurate models to investigate new techniques for assessing the safety of operations involving hydrogen to ensure the inherent hazards of using hydrogen do not negate the benefit of reduced CO 2 emissions.

The paper presents preliminary findings of an investigation into modelling lifted hydrogen jet fires using the boundary layer equations. The use of the boundary layer equations means that any model calibration is rigorous confident that the numerical error is negligible. A number of lift-off models based on the flamelet quenching concept and the turbulence time-scale are implemented and evaluated using available experimental data taken from the open literature. The lift-off models based on flamelet quenching use multi-flamelet libraries for the state relationships used to calculate the mean flow properties. Initial results suggest that a lift-off model based on the flamelet quenching concept incorporating the small-scale strain rate gives reasonable agreement with the measured lift-off compared to the other modelling approaches considered. This is an interesting result as it contrasts with earlier studies of methane and propane lifted jets, where the large-scale strain rate gave better agreement between observation and theory. This would suggest the appropriate strain-rate model for a particular fuel and jet could be related to the residence time of the jet or the bulk strain rate.

INTRODUCTION

The current interest in hydrogen as an alternative fuel to fossil fuels as a strategy for CO 2 reduction, the so called 'hydrogen economy' is a driver for research into hydrogen safety such that safety systems can be designed and hazards assessed. For example any vehicle fuelled by hydrogen is likely to have a pressure relief system in case of accident. Should the relief valve lift then a jet of hydrogen will form. If the jet of hydrogen finds a source of ignition it will form a lifted jet fire based on the typical valve diameters and source velocities encountered. Therefore the understanding of lifted jet fires is of interest to industrial combustion engineers. Turbulent jet flames occur in many areas of industry either by accident such as following the ignition of a leak from a high pressure plant processing hydrocarbons or in a controlled environment such as a furnace or industrial burner. In a turbulent jet the flow is characterised by the diameter of the nozzle or pipe, the fuel composition and the speed of the fuel at the nozzle exit. For a given fuel and nozzle diameter, if the jet speed is sufficiently low the flame is attached or rim-stabilised. If the jet speed is gradually increased then at some critical velocity the jet lifts off and the flame structure fundamentally changes. Immediately downstream of the nozzle the flow is non-reacting, further downstream combustion is initiated at a fluctuating flame base followed by the main body of the jet flame. Figure 1 shows a schematic diagram of the lift-off region of a lifted flame with the mean lift-off height labelled. A wide range of small-scale, Kalghati (1984) and [START_REF] Pitts | Importance of isothermal mixing processes to the understanding of liftoff and blowout of turbulent jet diffusion flames[END_REF] and large-scale experimental studies of lifted fires, Birch and Hargreave (1988) exist in the open literature. There are also a number of theoretical studies, [START_REF] Sanders | Modeling and calculation of turbulent lifted diffusion flames[END_REF], [START_REF] Cumber | Modelling lifted methane jet fires using the boundary layer equations[END_REF] being examples.

The computational modelling of rim-stabilised fires is mature with good agreement between the predicted and measured mean temperature and major chemical species established, [START_REF] Gore | Spectral and total radiation properties of turbulent hydrogen/ air diffusion flames[END_REF] for example. This is not the case for lifted jets where computational studies have focussed on understanding the mechanism for the location of the combusting flame base and good agreement between predicted and measured flow properties

U 0 z L Fluctuating Flame Base

Stoichiometric

Contour is not universal. The possible mechanisms for the location of the lift-off height will be considered further below. The mechanism for flame lift-off is a competition between chemical and mixing processes characterised by a chemical reaction time and a flow residence time.

The ratio of these time-scales is defined as the Damköhler number. As the jet velocity increases the flow residence time decreases. At the critical velocity or Damköhler number there is insufficient time for the chemical reaction to take place and lift-off occurs, [START_REF] Sanders | Modeling and calculation of turbulent lifted diffusion flames[END_REF]. This is the accepted physical explanation for flame lift-off but the mechanism involved in flame stabilisation at the lift-off height after over three decades of research is still debated vigorously. The early accepted theory was the lift-off height was determined by the downstream location on the stoichiometric contour where the mean jet velocity equals the turbulent burning velocity. The location of the flame base is governed by the premixed nature of the fuel air mixture. Computational models using this criterion were able to predict lift-off heights for the jet fires considered [START_REF] Vanquickenborne | The stabilisation mechanism of lifted diffusion flames[END_REF]. Experimentalists were also able to correlate measured lift-off heights for a range of fuels using the premixed assumption and the turbulent burning velocity, Kalghati (1984). [START_REF] Peters | Lift-off characteristics of turbulent jet flames[END_REF] were the first to question this theory suggesting that insufficient premixing could occur for the premixed theory to be valid in the lifted jet fires they considered. Peters and Williams proposed another mechanism, based on laminar flamelet quenching that relies on an analysis of laminar flamelets distorted by the turbulence field. At the fluctuating flame base the lifted flame can be considered as an ensemble of laminar flamelets in the instantaneous field; the laminar flamelets varying from air-fuel mixtures at or near stoichiometric conditions to fully entrained air. Each of the flamelets can be considered to be burning or non-reacting depending on the local instantaneous composition and level of turbulent distortion or stretch the flamelet is subjected to. The location of the flame base is determined by the proportion of burning to non-reacting flamelets. A candidate property for characterising the degree of stretch of the flamelets is the scalar dissipation rate, [START_REF] Peters | Laminar diffusion flamelet models in non-premixed turbulent combustion[END_REF], although [START_REF] Sanders | Modeling and calculation of turbulent lifted diffusion flames[END_REF] present a convincing argument against the scalar dissipation rate and suggest the strain rate as the appropriate parameter. [START_REF] Muller | Partially premixed turbulent flame propagation in jet flames[END_REF] use a model for the strain rate of the largest eddies to predict the lift-off height, whereas Sanders and Lamers favour the small-scale strain rate. Sanders and Lamers base their conclusions on a calibration of lift-off height using isothermal simulations, the argument being the lift-off height is primarily determined by the flow upstream of the flame base. As well as the two theories discussed above [START_REF] Pitts | Importance of isothermal mixing processes to the understanding of liftoff and blowout of turbulent jet diffusion flames[END_REF] analysed a number of lifted jet fires and suggested that the location of the lift-off height is determined by large-scale rather than smallscale structures in the jet, and isothermal mixing upstream of the lift-off region is important in determining the location of the flame base. Even though Pitts presents a compelling case for his favoured mechanism he still states that models based on a premixed theory or small-scale turbulent structures can not be discounted. More recently [START_REF] Upatnicks | Lift-off turbulent jet flames-Assessment of edge flame and other concepts using cinema-PIV[END_REF] using a high speed photography technique, Cinema PIV identified another mechanism for the location of the lift-off height based on an analysis of the edge of the flame. One reason for the lack of clarity as to the correct mechanism is all theories show some degree of agreement with observation, and in many experiments there is considerable overlap in prerequisite conditions for each theory to apply.

In this paper a model originally proposed by [START_REF] Sanders | Modeling and calculation of turbulent lifted diffusion flames[END_REF] and modified by [START_REF] Cumber | Modelling lifted methane jet fires using the boundary layer equations[END_REF] based on the laminar flamelet quenching concept is extended to predict lifted hydrogen jet fires. The original model was applied to lifted methane jet fires, [START_REF] Cumber | Modelling lifted methane jet fires using the boundary layer equations[END_REF] demonstrating that for lifted methane jet fires the large scale strain rate, the eddy strain rate was the important parameter in determining the lift-off height rather than the strain rate at the smaller scales nearer the molecular scale. The lifted jet fire methodology presented in [START_REF] Cumber | Modelling lifted methane jet fires using the boundary layer equations[END_REF] is based on a parabolic flow model and as such the use of laminar flamelet quenching is a convenient lift-off model. The methodology presented with minor modification could be implemented with other lift-off models. The overall computational framework is a parabolic flow solver, utilising the boundary layer equations, similar to GENMIX, [START_REF] Spalding | GENMIX: A general computer program for two-dimensional parabolic phenomena[END_REF], extended to include elliptic features found in lifted fires. In lifted jet fires the important elliptic features are the feedback mechanisms downstream of the lift-off region that determines the lift-off height and local flame structure. In the context of the lift-off model implemented here the elliptic feedback mechanism is introduced through the probability of burning field. The details of implementation are considered further below. As jet fires have a dominant flow direction it is common for the flame structure of a rim-stabilised jet fire to be calculated in this way, [START_REF] Gore | Spectral and total radiation properties of turbulent hydrogen/ air diffusion flames[END_REF]. The computational advantages of a parabolic flow model over a fully elliptic solver are clear, as no iteration cycle of a pressure correction algorithm is required in its solution as the dominant flow direction has a time like quality allowing a marching procedure in the dominant flow direction. This makes it possible to calculate complex flame structures using readily available computational resources. For example [START_REF] Wang | PDF modelling of turbulent non-premixed combustion with detailed chemistry[END_REF] have reported the computation of a laboratory scale rim-stabilised flame, the simulation is calculated using a Probability Density Function (PDF) Transport Model with detailed chemistry, including 53 species and 32 elemental reactions and a multi-time-scale k-ε turbulence model. Either of these modelling approaches in isolation would make the jet fire calculation prohibitively computationally intensive when combined with an elliptic solver on a standard PC. Using a flow solver based on the boundary layer equations Wang and Chen report a run-time of the order of 1 week on a relatively high specification PC. The use of a parabolic solver is not only an issue of convenience as it makes it relatively simple to demonstrate mesh independence in any predicted flow fields presented. In addition any model calibration is purely dependent on the quality of the experimental data used in the calibration. The issue of robust calibration is a particularly important one as without it, model development looses rigour. [START_REF] Cumber | Modelling lifted methane jet fires using the boundary layer equations[END_REF] showed that [START_REF] Sanders | Modeling and calculation of turbulent lifted diffusion flames[END_REF] adoption of the small-scale strain rate model for lifted methane jet fires was incorrect as the large scale strain rate model gave far closer agreement with observed lift-off heights once both models were calibrated appropriately. In the next section the mathematical model is presented.

MATHEMATICAL MODEL

The basis of the flow equations is the parabolized Favre averaged Navier Stokes equations in an axisymmetric coordinate system. The system is closed using a variant of the k-ε turbulence model.

( ) ε ρ σ µ ρ ρ - + ∂ ∂ ∂ ∂ = ∂ ∂ + ∂ ∂ k k eff P r k r r r r k V r r z k U 1 1 ~ ( ) ε ε ρ ε σ µ ε ρ ε ρ ε 2 1 1 1 ~C P C k r r r r r V r r z U k eff - + ∂ ∂ ∂ ∂ = ∂ ∂ + ∂ ∂ (1) ~ ∂ ∂ = r U P eff k µ ε ρ µ µ µ 2 k C l eff + = 3 1 4 . 3 4 . 1 cl dz U d k C - = ε C µ =0.09, C 2 =1.84, σ k =1, and σ ε =1.3
The version of the k-ε turbulence model given above is a modification of the typically implemented variant of the k-ε turbulence model (C 1 =1.44, C 2 =1.92) to take account of the round jet/ plane jet anomaly, [START_REF] Pope | An explanation of turbulent round-jet/ plane-jet anomaly[END_REF] where the spreading rate of round jets tends to be overpredicted by the 'standard' version of the turbulence model. This is a well known limitation of most two equation turbulence models and 2 nd moment closure models unless some modification is introduced to account for the reduced spreading rate. The modification to the k-ε model introduced above is due to [START_REF] Morse | Axisymmetric turbulent shear flows with and without swirl[END_REF], and has been used successfully in previous rim-stabilised fire simulations, [START_REF] Moss | Flow field modelling of soot formation at elevated pressure[END_REF]. The axisymmetric correction is used here as it gives an appropriate balance between model complexity and predictive capability.

Turbulent Combustion Model

The turbulent combustion model is a laminar flamelet combustion model, with two flamelet libraries, one for combustion and the other for isothermal mixing. Combustion is assumed to be infinitely fast with a prescribed probability density function, a β function, [START_REF] Fairweather | Predictions of a turbulent reacting jet in a cross-flow[END_REF]. The shape of the β function at any spatial location is determined by a conserved scalar, the mixture fraction f and its variance f ''2 , which are calculated using modelled transport equations, [START_REF] Fairweather | Predictions of a turbulent reacting jet in a cross-flow[END_REF]. The combusting flamelet is derived by calculating the adiabatic flame temperature at stoichiometric conditions, based on the sum of the chemical and sensible enthalpies are conserved during the reaction.

( ) ( ) ( ) ( ) ( ) (
)

0 1 1 = - + - - = = prod react N i ref i adia i i N i ref i amb i i fu T h T h x T h T h x H
(2)

The flamelet for equivalence ratios other than the stoichiometric ratio is specified by fitting temperature measurements in laminar flames and instantaneous temperature measurements in turbulent flames, [START_REF] Gore | Spectral and total radiation properties of turbulent hydrogen/ air diffusion flames[END_REF]. This gives a flamelet similar to one calculated by assuming equilibrium apart from the measurements used in the curve fitting will have a radiative heat loss associated with them. This is desirable as the flames considered in this paper also have a radiative heat loss. The equilibrium assumption for the hydrogen reaction is a reasonable assumption given the reactivity of hydrogen in air. Any mean property can be calculated as a weighted average of the burning and isothermal mixing flamelet weighted by the pdf and integrated over instantaneous mixture fraction,

( ) ( ) ( ) ( ) ( ) - + = 1 0 1 0 1 ~df f f P P P df f f P P P m b d b b d ϕ ϕ ϕ (3)
where P b is a probability of burning defined below and P d is a probability that the axial location is above the fluctuating flame base.

( )

∞ > = 0 L L base d dz z z P P (4)
Sanders and Lamers prescribe the pdf for the location of the instantaneous flame base, z L to have a triangular shape, the apex located at the mean lift-off height and the base of the triangle is taken to be five diameters. It should be noted that this is an assumption of convenience rather than one based on observation, however a sensitivity study has shown that the overall flame structure is insensitive to this aspect of the model. The mean density and mean adiabatic temperature are given by the relations, ( ) ( ) ( ) ( )

1 1 0 1 0 1 - - + = df f f P P P df f f P P P m b d b b d ρ ρ ρ (5) ( ) ( ) ( ) ( ) ( ) ( ) ( ) - + = 1 0 1 0 1 df f f T f P P P df f f T f P P P T m m m b d b b b b d adia ρ ρ ρ ρ
To close the system the probability of burning and the mean lift-off height must be modelled. The probability of burning is given by,

( ) = q s b ds s P P 0 (6)
s is the strain rate, s q is the quenching strain rate and P(s) is a quasi Gaussian pdf for the strain rate. [START_REF] Sanders | Modeling and calculation of turbulent lifted diffusion flames[END_REF] give a quench strain rate of 565 sec -1 derived from an analysis of diluted methane-air counter flow diffusion flames. However a sensitivity study indicates that lift-off predictions are fairly insensitive to the quench value used, [START_REF] Cumber | Modelling lifted methane jet fires using the boundary layer equations[END_REF]. The strain rate can either be taken to be the strain rate of the small-scale turbulence,

2 / 1 1 , 2 = l s C s ν ε (7)
Or the strain rate of the large-scale turbulence.

k C s s ε 2 , = (8)
Where C s,1 and C s,2 are calibration constants. The laminar dynamic viscosity is given by Sutherland's law. The mean lift-off height is prescribed using percolation theory, [START_REF] Peters | Lift-off characteristics of turbulent jet flames[END_REF] and the probability of burning on the stoichiometric contour, The turbulence time-scale is evaluated on the jet axis. For methane jets t fresh is 5x10 -3 seconds. In this study it will be fixed by calibrating against one of Kalghati's (1984) jet fires. This lift-off model is considered here as it is consistent with a parabolic flow model without modification and has been used in the past to predict the lift-off height of methane jet fires in a cross-flow, [START_REF] Fairweather | Predictions of a turbulent reacting jet in a cross-flow[END_REF]. In this model the transition from no reaction to combustion is instantaneous at the lift-off height.

Boundary Conditions and Numerical Parameters

In all of the simulations presented below the bulk inlet conditions are given by the nozzle diameter and the average source velocity. The mean stream-wise velocity distribution and radial velocity distribution are taken to be consistent with fully developed pipe flow that is a 1 /7th power law is prescribed for the stream-wise velocity distribution and zero for the radial velocity component. Nozzle exit turbulence profiles are prescribed to be given by equilibrium profiles for fully developed turbulent pipe flow. For some of the jets considered the source velocity is above the local sound speed. To accommodate the compressibility of the flow the source condition is modelled as a nozzle and the jet expanded isentropically to atmospheric pressure. For all simulations in this article 80 control volumes in the radial co-ordinate direction spanning the jet radius are used to calculate the flame structure, with a maximum fractional step in the axial direction of less than 1% of the radial control volume spacing. A number of simulations using 160 control volumes were also completed to confirm that the predictions presented are independent of further mesh refinement. It is estimated that the predicted lift-off heights are within 2% of the fully mesh converged values.

DESCRIPTION OF THE NUMERICAL ALGORITHM

When considering how a lifted flame differs from a rim-stabilised flame in the context of the mathematical model described above, the key difference is the composite probability of burning field,

P burn = P b P d (11)
In a rim-stabilised flame P burn is one everywhere, whereas in a lifted flame this takes a value between zero and one. Zero below the lift-off region, between zero and one in the liftoff region and one above it. Therefore an algorithm based on a "guess and correct" approach for the composite probability of burning field suggests itself as one way of extending a parabolic flow model suitable for simulating rim-stabilised jet fires to simulate lifted jet fires. An overview of the algorithm in the form of a flow chart is given in Figure 2. To initialise the process some estimate of the composite probability of burning field must be prescribed. There are a number of possibilities, for example the fire could start off as being rimstabilised, hence the initial composite probability of burning field could be set to one everywhere. An alternative choice is to prescribe the jet to be isothermal, with a composite probability of burning of zero everywhere. Both options have some appeal as in the initial rim-stabilised fire approach (P burn,0 =1), this is how a lifted jet fire would initiate, alternatively the analysis of [START_REF] Pitts | Importance of isothermal mixing processes to the understanding of liftoff and blowout of turbulent jet diffusion flames[END_REF] work on isothermal jets indicates that the non-reacting isothermal region upstream of the flame base is important in the flame stabilisation process, favouring the second approach(P burn,0 =0). In [START_REF] Cumber | Modelling lifted methane jet fires using the boundary layer equations[END_REF] both approaches were evaluated and taking the initial lift-off height to be zero i.e. a rim stabilised fire gave a marginal improvement in convergence rate compared to starting with a non-combusting isothermal simulation of the jet. Once the composite probability of burning field is prescribed, an estimate of the flame structure can be calculated by running the parabolic flow model. As the flame structure is calculated flow fields required to estimate the lift-off height and the composite probability of burning field, such as the mean mixture fraction, mean temperature, turbulence kinetic energy, and the dissipation rate of the turbulence kinetic energy, are stored. Once the flame structure is complete, a new estimate of the lift-off height is calculated and the composite probability of burning field recalculated. If the change in liftoff height is small the algorithm is terminated otherwise the flame structure is recalculated and the process continues.

Large and Small Scale Strain Rate

An iterative procedure poses the important question of 'how fast does the algorithm converge and what is an appropriate convergence criterion?

To illustrate the convergence behaviour of the algorithm a typical jet fire is used. The jet fire simulated is the same in each case, (d=2.95mm, U 0 =1190 m/s). Figure 3a shows the convergence history for two simulations using the two different strain rate models. The relative change in mean lift-off height is plotted against the iteration number. Considering the simulation using the large-scale strain rate sub-model convergence is monotonic for the first twenty iterations after which round-off error prevents further convergence, however the differences in predicted lift-off heights for successive iterations using the large scale strain rate model is less than 10 -5 % and changes of less than 1% is sufficient for the lift-off height to be converged to the visual resolution of Figure 3b. In Figure 3b the lift-off height as a function of iteration number is shown for the two simulations. Figure 3a and Figure 3b taken together suggest for the large scale strain rate model the algorithm takes around 5-6 iterations to converge.

Considering the convergence history for the small-scale strain rate simulation, the situation is less satisfactory. The convergence rate is slower than the large scale strain model simulation and round off error has a larger effect on the convergence history. Oscillations in the relative difference occur with peak relative differences of the order of 10 -3 % although a relative difference of 0.5% is sufficient for the methodology to have converged to the resolution of Figure 3b. The slower convergence means that for the small scale strain rate model the flame structure takes of the order of 10 iterations to converge. The relatively poor convergence behaviour of the simulations using the small-scale strain rate sub-model is believed to be due to its dependence on the dissipation rate of turbulence and the kinematic viscosity. These two properties are weakly correlated and when the lift-off height changes these two properties respond differently to the change in the mean temperature field. This convergence problem is lessened when more control volumes are introduced in the radial coordinate direction, with of course an increase in computational cost. The large-scale strain rate sub-model does not suffer from this problem as its dependency on turbulence parameters alone means that it responds in a consistent way to changes in the mean temperature field brought about by changes in the predicted mean lift-off height.

Lift-off Height Model Calibration

The three lift-off models are calibrated for lifted hydrogen jet fires, using one of Kalahari's experiments, d=2.95 mm and U 0 =1190 m/s. This particular jet was chosen as it had one of the highest source velocities that was still subsonic. This has the advantage that the isentropic nozzle approximation does not influence the calibration. Fitting the models to this data gives values of C s,1 =0.0023 and C s,2 =0.0073 and t fresh =7.9x10 -5 s (12)

for the calibration constants in the strain rate sub-models and the turbulence time-scale submodel respectively. predicted lift-off height.

MODEL PREDICTIONS

Rim-stabilised Jet Fires

To gain confidence in the model two rim-stabilised hydrogen fires are considered first as for the lifted jets the lift-off height alone is used for validation. The two jet fires considered are presented in [START_REF] Gore | Spectral and total radiation properties of turbulent hydrogen/ air diffusion flames[END_REF], for both the nozzle diameter is 5mm and the source velocity is 66.3 m/s (Re=3,000) and 108.4 m/s (Re=5,722) . Figure 4 shows the predicted and measured axial temperature distribution on the jet axis for the two jet fires. For the Re=3,000 jet fire the agreement between the model and the measured temperature distribution is very good. For the higher Reynolds number, Re=5722, there is a tendency to under predict the temperature field by approximately 10%. This level of agreement is marginally superior to that reported by [START_REF] Gore | Spectral and total radiation properties of turbulent hydrogen/ air diffusion flames[END_REF]. This is not surprising as for rim-stabilised fires the models are similar, differing in the mesh resolution, a finer mesh is used here, the turbulence model implemented and the temperature flamelet. All of these differences are points of detail that do not radically change the performance of the model.

Lifted Jet Fires

Before presenting the validation of the model for lifted jets it is of interest to explore how the lift-off models influence the local flame structure in the lift-off region. Figure 5 shows the predicted temperature field in the lift-off region. Three simulations are shown, each differing in the sub-model used to predict lift-off. This jet fire (d=2.95mm, U 0 =1190 m/s) is used for illustrative purposes as it is the jet fire used for calibrating the models and so in each simulation the lift-off height is the same. For each temperature field 5 equi-spaced contours are shown, starting at 600K increasing in increments of 300K, the lowest temperature nominally defining the boundary of the jet. To aid interpretation an arrow on each of the predicted temperature fields indicates where the mean temperature is increasing. The main points of note are that the temperature fields predicted using the two strain rate models are remarkably similar. This is different to previous simulations of lifted methane jets, [START_REF] Cumber | Modelling lifted methane jet fires using the boundary layer equations[END_REF]. Another point of note is the flame structure in both simulations is physically plausible as the combusting region is initiated off axis where the stoichiometric contour is located and the mean temperature increases relatively smoothly with increasing axial distance. Considering the predicted temperature field calculated using the turbulence time scale to predict lift-off the predictions are not physically plausible as combustion is initiated as a step change in flow conditions at the lift-off height across the width of the jet. Taking this analysis one stage further would result in an over prediction of the received radiation field in the near field. However considering all three simulations someway downstream of the lift-off region, the predicted temperature field is insensitive to how the liftoff region is modelled.

Figure 6 shows the lift-off height as a function of source velocity for 4 different nozzle diameters. The lift-off height for a given fuel is primarily dependent on the source velocity with the nozzle diameter having a much weaker influence on the lift-off height. As discussed above, to make a fair comparison of the three lift-off models they were all calibrated against a single datum, d=2.95mm, U 0 =1190 m/s, see Figure 6c. Considering the small scale strain rate model predictions of lift-off height, this model gives the overall best agreement with the measured lift-off height measurements as the lift-off measurements for nozzle diameters of 1.74mm, 2.95mm Is underpredicted making it a poor choice for application to jet fires outside the envelope of Kalghati's experiments. The turbulence-timescale model exhibits the worst performance of all three lift-off models, as for each nozzle diameter it predicts the lift-off height for one source velocity. Overall all models tend to underpredict the lift-off height for small diameters and tend to overpredict the lift-off height for larger diameter nozzles.

The superiority of the small-scale strain rate model compared to the large scale strain rate for lifted hydrogen jet fires is particularly interesting given previous experience of simulating lifted methane jet fires, [START_REF] Cumber | Modelling lifted methane jet fires using the boundary layer equations[END_REF] where the opposite was the case. The reason for this change in behaviour is being investigated by comparing the strain rate fields for methane and hydrogen jet fires.

Considering how the lift-off height predictions could be improved, the most obvious weakness is the isentropic nozzle approximation used to account for compressibility effects for sonic source conditions. For the smaller diameter nozzles where the higher source velocities are possible, combining the flame structure model with a compressible flow solver, Cleaver et al. (2001) to resolve the shock structure downstream of the nozzle should improve model performance. As well as resolving the shocks this will also introduce compressibility effects on the turbulence field as turbulent mixing tends to be reduced in high Mach number flows. The mechanism for this is believed to be shocks in the instantaneous flow fields, so called 'eddy shocklets'. This effect would tend to increase the predicted lift-off height. Note this would improve the agreement for the smaller diameter nozzle jets for source velocities above 1200 m/s. For the larger nozzle jets most measured lift-off heights are for source velocities below Mach numbers where significant compressibility effects on turbulent mixing would occur.

CONCLUSION

In this article an algorithm for calculating the flame structure of lifted hydrogen jet fires is presented. The algorithm is based on an extension of a fire model used to simulate rimstabilised flames. An important characteristic of the lifted jet fire methodology is the flame structure is calculated using the boundary layer approximation. The advantage of this is finite volume mesh independent predictions of the mean flow fields can be calculated on readily available computer resources. The simulations presented here typically required of the order of 10 minutes run-time on a PC with a 3 GHz Pentium IV processor. This has allowed a rigorous calibration of the lift-off models considered.

Three lift-off models were implemented and evaluated. The simplest uses a turbulence time-scale threshold to determine the lift-off height. This model is consistent with the boundary layer equations and can be implemented easily, but leads to unrealistic predictions of the mean temperature in the lift-off region. The other two lift-off models were formulated by [START_REF] Sanders | Modeling and calculation of turbulent lifted diffusion flames[END_REF] and are based on laminar flamelet quenching for the flame stabilisation mechanism at the flame base. The two lift-off models are differentiated by the strain rate sub-model implemented. One is a small-scale strain rate sub-model the other is a large-scale strain rate sub-model. The lift-off models were evaluated by comparing predicted lift-off predictions and of Kalghati(1984).

In this study the small-scale strain rate sub-model proved to be more accurate of the two strain rate sub-models. This is in contrast to previous work on lifted methane jets, [START_REF] Cumber | Modelling lifted methane jet fires using the boundary layer equations[END_REF] were the large scale strain rate model was superior.

The next step in the development of a hydrogen jet fire model for safety assessment of hydrogen processing operations is to validate the lift-off models using a wider range of measurements, both other lift-off data and detailed measurements of the mean temperature and chemical species in the lift-off region. This task is in hand and will be reported in a subsequent publication.

NOMENCLATURE

Figure 1 -

 1 Figure 1 -Schematic diagram of a lifted flame in the lift-off region.

  is an approximation to the error function and implicitly gives the percolation threshold, P c , see[START_REF] Cumber | Modelling lifted methane jet fires using the boundary layer equations[END_REF] for more details of the lift-off model.A third lift-off model is considered, based on a correlation for lift-off height using a turbulence time-scale threshold.

StartFigure 2 -

 2 Figure 2-Flow chart of the lifted jet fire methodology.

Figure 3

 3 Figure 3 Typical convergence histories of a lifted jet fire simulation a) Relative difference, b)predicted lift-off height.

Figure 4 -

 4 Figure 4 -Axial temperature distributions for two rim-stabilised jet fires, a) Re=3,000, b) Re=5,722.

Figure 5 -

 5 Figure 5 -Predicted temperature field in the lift-off region, (U 0 =1190 m/s, d=2.95mm) a) Small scale strain rate, b) Large scale strain rate and c) turbulence time-scale model. (Arrow indicates increasing temperature.)

Figure 6 -

 6 Figure 6 -Lift-off height vs. source velocity for a) d=1.08mm, b) d=1.74mm. c) d=2.95mm and d) d=4.06mm.

  amb Ambient value(-) q Quench value (-) st Property at stoichiometric conditions (-) 0Initial condition or source condition(-)

  C 1 , C 2 , C µ Parameters in the k-ε turbulence model (-) C s,1 , C s,2 Calibration constants in strain rate models(-)

	d	Nozzle diameter (mm)
	H fu	Heat of combustion (J/Mol)
	h i k	Specific enthalpy of specifies i (J/Mol) Turbulence kinetic energy (m 2 /s 2 )
	P	Probability density function (-)
	P b	Probability of burning (-)
	P burn	Composite probability of burning (-)
	P c	Percolation threshold (-)
	P d P k	Probability of burning relating to the location of the fluctuating flame base (-) Turbulence production source term (kg/ m s 3 )
	Re	Reynolds Number (-)
	r s	Radial co-ordinate (m) Strain rate (s -1 )
	t fresh