**4. Conclusions**

**REIT** **A. 2001m5–2017m7**

Beta R-squared Prob value (beta)

**B. 2001m5–2005m12**

Beta R-squared Prob value (beta)

**C. 2006m1–2009m1**

Beta R-squared Prob value (beta)

**D. 2009m1–2013m12**

Beta R-squared Prob value

0.00

0.00

0.00

0.00

0.00

0.00

0.00

0.00

0.00

0.10

0.00

(beta)

**E. 2014m1–2017m7**

Beta R-squared Prob value (beta)

**Table 6.**

Comparison of beta values for different sub-periods.

0.00

0.03

0.98

0.02

0.00 \*\*\*Shows that beta coefficient is significant at 1% level, \*\*shows beta coefficient is significant at 5% level, \*shows the significance at 10% level.

0.00

0.01

0.30

0.79

0.00

0.00

0.66

0.57

0.00

0.50

0.12

0.00

0.12

0.47

0.29

0.18

0.03

0.00

0.40

0.24

0.00

0.01

0.25

0.68\*\*\*

0.32\*\*

−0.01

0.62\*\*

0.90\*\*\*

0.55\*\*\*

0.50\*\*\* 0.35

0.07

0.75\*\*\*

0.91\*\*\*

0.15

−0.12

1.05\*\*\*

0.64

0.26

0.28

0.40

0.49

0.52

0.23

0.73

0.60

0.05

0.57

0.96\*\*\*

0.96\*\*\*

0.83\*\*\*

0.96\*\*\*

0.83\*\*\*

1.04\*\*\*

1.11\*\*\*

1.56\*\*\*

1.23\*\*\*

0.46\*

1.12\*\*\*

0.00

0.00

0.00

0.77

0.46

0.41

0.97\*\*\*

0.60\*\*\*

0.83\*\*\*

1.09\*\*\*

0.62 0.00

0.00

0.32

0.85\*\*\*

0.80\*\*\*

0.35 0.00

0.00

0.23

1.11\*\*\*

0.00

0.79

0.98\*\*\*

0.80\*\*\*

0.49 0.00

1.04\*\*\*

0.70 0.00

0.00

0.52

0.95\*\*\*

0.78\*\*\*

0.36 0.00

0.00

0.38

1.23\*\*\*

296 Financial Management from an Emerging Market Perspective

0.00

0.00

0.82

0.00

0.00

0.00

0.00

0.00

0.00

0.00

0.00

0.00

0.64

0.00

0.73

0.26

0.00

0.38

0.42

0.28

0.61

0.19

0.31

0.58

0.43

0.18

0.01

0.35

0.97\*\*\*

0.70\*\*\*

0.04

0.83\*\*\*

0.93\*\*\*

0.49\*\*\*

0.98\*\*\*

0.81\*\*\*

0.94\*\*\*

1.25\*\*\*

1.12\*\*\*

0.70\*\*\*

−0.08

1.18\*\*\*

**BIST REIT**

**AKMGY**

**AKSGYO**

**ALGYO**

**EKGYO**

**HLGYO**

**ISGYO**

**KLGYO**

**NUGYO**

**SNGYO**

**TRGYO**

**VKGYO**

**YGGYO**

**PEGYO**

The main purpose of the study is to analyse the performance stability of REIT index and individual REITs over different sub-periods. The sample for the study contains 12 individual REITs that have weights larger than 1% in the REIT index and the REIT index that survived for the whole REIT history in Turkish stock market, thus, in total, 13 individual REITs and REIT index itself. The performance of the REITs is compared to mainly that of BIST 100 Index which is

3 The present study differs at many aspects from the study by Altınsoy et al. [1]. This study uses monthly returns of both REIT index and individual REITs and uses simple CAPM over a longer period. The study carried out by Altınsoy et al. [1] uses weekly and daily returns for only REIT index for the period 2002–2009 and calculates time varying betas by using three different methods. It does not analyses the individual REITs.

4 The author has also run regressions with several macroeconomic factors including GDP. Only for the second period, the financial crisis period, 50% of the REIT returns has significant negative relation with GDP. For all other periods, the relationship between REIT returns and GDP is not consistently significant. The best can be said about the relationship is being unclear. These multifactor regressions are also consistent with the declared conclusions on alpha and beta behaviour. The results are not presented here for the sake of brevity and integrity.

accepted to be the reference portfolio in the study. Data employed in this study is the monthly returns for a 14-year period and Treasury bill rates are used as a proxy for the risk-free rate.

For analysing the performance stability, three different risk-adjusted measures, namely, Sharpe ratio, Treynor ratio and Jensen's alpha, are employed for four different periods. These periods are determined with respect to important regulatory changes in the REIT market and also to economic states of the country. The results show that although Treynor and Sharpe ratios use different risk measures, volatility and beta, they rank the REITs consistently for the high-growth periods. However, the rankings are not that consistent in low-growth periods and even they may contrast significantly. In terms of stability, whether the performance of the REITs changes over time, it is observed that the rankings of the REITs in the last period virtually are shuffled up with respect to the previous period. These results basically show us that regulatory changes almost have no impact on the performances of the REITs. Otherwise, we would observe stable improvements in the performances of the REITs throughout the years. In consistent with the previous literature, it is possible that state of the economy (whether it is a low-growth or high-growth state) is more important in terms of the performances of REITs. However, time-varying behaviour of betas does not tell us the same story. The betas do not consistently decrease from the high-growth to low-growth period. Thus, it is not possible to easily relate decreasing betas to economic states as it is done in the previous studies. There may be other important reasons for time varying betas. As a future study, it looks important to analyse why the betas of REITs are decreasing over time.


**Fernando Luiz Sacomano Filho 1,\*, Louis Dressler 2, Arash Hosseinzadeh 2, Amsini Sadiki 2 and Guenther Carlos Krieger Filho 1**

Received: 30 June 2019; Accepted: 29 October 2019; Published: 31 October 2019

**Abstract:** Evaporative cooling effects and turbulence flame interaction are analyzed in the large eddy simulation (LES) context for an ethanol turbulent spray flame. Investigations are conducted with the artificially thickened flame (ATF) approach coupled with an extension of the mixture adaptive thickening procedure to account for variations of enthalpy. Droplets are tracked in a Euler–Lagrangian framework, in which an evaporation model accounting for the inter-phase non-equilibrium is applied. The chemistry is tabulated following the flamelet generated manifold (FGM) method. Enthalpy variations are incorporated in the resulting FGM database in a universal fashion, which is not limited to the heat losses caused by evaporative cooling effects. The relevance of the evaporative cooling is evaluated with a typically applied setting for a flame surface wrinkling model. Using one of the resulting cases from the evaporative cooling analysis as a reference, the importance of the flame wrinkling modeling is studied. Besides its novelty, the completeness of the proposed modeling strategy allows a significant contribution to the understanding of the most relevant phenomena for the turbulent spray combustion modeling.

**Keywords:** spray combustion; evaporative cooling; flame surface wrinkling modeling; thickened flame; flamelet generated manifold

## **1. Introduction**

Turbulent spray combustion features the unsteady interaction of many strongly coupled physical phenomena varying in a broad range of time and length scales. The complexity resulting from such interacting phenomena can be appraised with the studies presented for instance by Jenny et al. [1], Gutheil [2], and Sacomano Filho [3]. Heat and mass transfer between phases and the effects of the turbulence on a flame, as well as their reciprocal, stand out among these coupled phenomena. Although both interactions are strongly coupled, they are not always comprehensively addressed in numerical simulations. The cooling down of the gas mixture caused by the evaporation process, namely the evaporative cooling (EC), is often neglected [4–8] or simplified [9,10] in the context of tabulated chemistry. Herein, the consideration of heat losses is not trivial and typically involves high computational efforts [11–13]. Commonly, the low volume fraction of diluted sprays are used to justify simulations without heat losses [14,15]. Nevertheless, as analyzed in [16] for one-dimensional flames propagating in droplet mists, the inclusion of the evaporative cooling is quite relevant in the spray combustion modeling. Yet the importance of the proper modeling of the turbulence-flame interaction is

a well known topic in the spray combustion community [1,4,6–10,13,17–21]. In view of the application of the artificially thickened flame (ATF) approach, many works have been using a similar value for the exponent of the power-law function used to model the flame surface wrinkling (FSW), namely 0.5. Recent investigations presented in [3,21] show that the characterization of the turbulence-spray flame interaction can be significantly improved according to the value choice of such a model parameter. However, this has not been consistently analyzed with spray flames accounting for the evaporative cooling up to the present study.

An initial approach to account for evaporative cooling effects is to consider heat losses in the mixture without including modifications in the chemistry table [9,10]. According to this formalism, influences of heat losses concentrate on the computation of droplet evaporation rates which indirectly interfere with the combustion process. As evaporation rates change, the mixture composition arriving at the reaction zone also changes. However, the influences of heat losses are not strictly accounted for in the chemistry. Knudsen et al. [13] proposed a low computational cost strategy to include heat losses exclusively caused by droplets in the flamelet-progress variable approach. They compute freely propagating and counter-flow flamelets to generate two tables (one based on freely propagating and another based on counter-flow flamelets) which would be chosen during the simulation evolution by means of a flame index. Instead of considering flamelets in diverse enthalpy levels, this variable is constrained with the mixture fraction by assuming that all the fuel existing in the fresh mixture comes from liquid evaporation. As a result, no enthalpy must be transported. This simplification can be applied in cases where the unitary Lewis number (*Le*) approach is adopted. The validity of this method for the description of the reaction evolution in mixtures richer than the mixture fraction corresponding to saturated vapor is not comprehensively investigated. As mentioned by van Oijen and de Goey [11], freely propagating flamelets become unphysical below a certain level of enthalpy. Therefore, care must be taken in order to suitably describe rich reactions with this method. Recently, Olguin and Gutheil [22] presented a new procedure to tabulate the chemistry involving spray flames. It does not only account for the effects of heat losses, but also to the presence of droplets in the reaction zone. Hu et al. [23] show that this method is quite promising to address turbulent spray flames. However, some obstacles can also be encountered there. For instance, the increase of table dimensions that request more storage and computational efforts, as well as the difficulty to find correspondence of the droplet trajectories in the counter-flow setup used to construct such a table [16] with the other found in a multi-dimensional flame. Similar to [22], Luo et al. [24] and Franzelli et al. [25] propose a method based on spray flamelets to characterize spray flames in view of tabulated chemistry. Despite being considered, effects of the evaporative cooling in [24] are limited since droplets are injected at constant (boiling point) temperature. As in [22], such limitation is not found in [25]. However, this last study considers a monodisperse droplet cloud, which is a more restrictive condition for the dispersed phase when compared with [22,24].

The ATF model stands out among methods based on deterministic approaches to characterize the turbulence-flame interaction in turbulent spray flames [9,10,17,20,21]. This model showed a great capability to predict the behavior of turbulent spray flames in the Eulerian–Eulerian two-phase flow approach (e.g., [9,17]) as well as in the Euler–Lagrangian framework (e.g., [10,20,21]). Initially proposed to model laminar premixed flames [26], the ATF was gradually expanded to address turbulent premixed [27,28] and stratified flames [29]. In view of turbulent spray combustion, a further extension to include effects of the dispersed phase in the Euler–Lagrangian context are recently presented in [21]. A common characteristic to all of these preceding works involving turbulent flames (i.e., [9,10,20,21]) refers to the usage of a global and presumed definition of the exponent of the power-law function (*β*) proposed by Charlette et al. [28] to model the FSW. According to the analysis conducted for premixed flames in [30], no universal value exists for this parameter, which may be determined by adjusting with experimental data or by some modeling approach (e.g., the dynamic method proposed by Charlette et al. [30]). Preliminary studies accounting for the dynamic modeling of such a model parameter (see [3]) showed promising results for the characterization of turbulence-spray flame interactions. However, further investigations are still necessary to turn this approach more consistent for general turbulent spray flame simulations. On the other hand, analysis of the influence of the model parameter *β* are needed for spray flames, which have not been comprehensively addressed in cases where heat losses are considered.

In the present work, the effectiveness of the evaporative cooling and the FSW modeling to characterize a diluted turbulent spray flame is analyzed. The chemistry tabulation method proposed by van Oijen and de Goey [11] and comprehensively tested for spray flames in [16] is employed here. Namely, effects of the evaporative cooling interfere with reaction rates, and the consideration of heat losses are not limited to the evaporative cooling in the flamelet generated manifold (FGM) table. Simulations are conducted following the large eddy simulation (LES) approach coupled with the ATF method. Besides the novelty associated with such a consistent consideration of heat losses in a turbulent spray flame, the mixture adaptive thickening approach [21] is extended to include effects of enthalpy variations. The flame EtF5 of the Sydney diluted spray flame burner [31] is selected as a benchmark to the subsequent analysis. Detailed chemistry effects are included by 57 species and 379 intermediary reactions mechanism of ethanol-air proposed by Marinov [32] by means of the FGM methodology. The unsteadiness arising from the turbulent dispersion of evaporating droplets are captured by a Euler–Lagrangian spray module relying on the LES approach. An evaporation model accounting for the inter-phase non-equilibrium is applied to describe the droplet evaporation process. A parametric study has been performed to evaluate the effectiveness of the evaporative cooling and the turbulence flame interaction modeling consideration.

The remaining structure of this paper is divided into three parts. An overview of the theoretical background is described in Section 2. In Section 3, simulation results are presented. Analyses of effects of the evaporative cooling and its coupling with the flame surface wrinkling modeling are systematically addressed. In the last part, final remarks and the main conclusions are summarized.

#### **2. Materials and Methods**

A Euler–Lagrangian approach was adopted to represent the two-phase flow. Herein, a full inter-phase two-way coupling was accounted for. Carrier gas-phase quantities were interpolated into droplets positions, while influences of the dispersed phase were introduced through source terms in computational cells.

#### *2.1. Gas Phase*

The turbulent motions of the carrier phase are described in the LES context following a variable-density low Mach number formulation. According to this approach, mass and momentum equations are described by

$$\frac{\partial \overline{\overline{\rho}}}{\partial t} + \frac{\partial \overline{\overline{\rho}} \widetilde{u}\_j}{\partial x\_j} = \overline{S\_m} \tag{1}$$

$$\frac{\partial \overline{\rho} \overline{u}\_{i}}{\partial t} + \frac{\partial \overline{\rho} \overline{u}\_{i} \overline{u}\_{j}}{\partial \mathbf{x}\_{j}} = \frac{\partial}{\partial \mathbf{x}\_{j}} \left( 2 \overline{\mu} \overline{S}\_{i\bar{j}} - \frac{2}{3} \overline{\mu} \frac{\partial \widetilde{u}\_{k}}{\partial \mathbf{x}\_{k}} \delta\_{\bar{j}\bar{j}} - \overline{\rho} \tau\_{i\bar{j}}^{\text{sys}} \right) - \frac{\partial \overline{p}}{\partial \mathbf{x}\_{i}} + \overline{\rho} g\_{i} + \overline{S\_{u,i}}.\tag{2}$$

The dependent filtered variables were obtained from spatial filtering as *ψ* = *ψ*- + *ψ*" with *ψ*- = *ρψ*/*ρ*¯. Over-bars and tildes express spatially filtered and density-weighted filtered values with a filter width Δ*mesh*, respectively, while double prime represents sub-grid scale (SGS) fluctuations. *ρ* is the mixture density , *t* time, *uj* components of velocity in *j* (*j* = 1, 2, 3) direction, *p* pressure, *xj* Cartesian coordinate in *j* direction, *μ* the dynamic viscosity, *gi* the component of gravitational acceleration, *δij* the Kronecker's delta, and *Sij* the strain rate. The term *Sm* corresponds to the introduction of mass from the droplets into the carrier phase, while *Su*,*i* is the source term of momentum due to the presence of the dispersed phase. Both follow the implementations presented by Chrigui et al. [33]. The SGS stress tensor *τsgs ij* is closed by means of the Smagorinsky model with the dynamic procedure of Germano et al. [34]. More details about the mathematical treatment given here to mass and momentum equations can be found in Sacomano Filho et al. [21].

#### 2.1.1. Mixture Formation and Combustion Modeling

In order to account for the evaporative cooling and general heat exchanges, three scalar quantities are used to characterize the mixture following the FGM method: the mixture fraction *Z*, the reaction progress variable *Ypv*, and the absolute enthalpy of the gas mixture *h*. The transport equation for them can be written in terms of a general variable *ψ* within the ATF modeling as

$$\frac{\partial(\overline{\rho}\widetilde{\psi})}{\partial t} + \frac{\partial(\overline{\rho}\widetilde{\psi}\widetilde{u}\_{j})}{\partial \mathbf{x}\_{j}} = \frac{\partial}{\partial \mathbf{x}\_{j}} \left[ \left( FE^{\*}\_{\Delta} \frac{\overline{\mu}}{\sigma\_{\Psi}} + (1 - \Omega) \frac{\mu\_{t}}{\sigma\_{t,\Psi}} \right) \frac{\partial \widetilde{\psi}}{\partial \mathbf{x}\_{j}} \right] + \frac{E^{\*}\_{\Delta}}{F} \dot{\omega}\_{\Psi} + \overline{S\_{\Psi}}.\tag{3}$$

where *μt* is the turbulent viscosity. For the transport equations of *Z* and *Ypv*, *σψ* and *σt*,*ψ* respectively represent the laminar *Sc* and the turbulent *Sct* Schmidt numbers (*Sc* = *Sct* = 0.7), while for *h* both respectively correspond to the laminar *Pr* and the turbulent *Prt* Prandtl numbers (*Pr* = *Prt* = 0.7). It is important to highlight that, with these values for *Sc* and *Pr* the unitary *Le* = *Pr*/*Sc* approach is maintained. The quantity *F* corresponds to the thickening factor, *E*∗ Δ to the efficiency function, and Ω denotes the flame sensor. Details about these quantities are addressed below.

The term *ω*˙ *ψ* = *ω*˙ *ψ Z*-,*Y pv*, *h* corresponds to the reaction rate for the *Ypv* whereas it is set to zero for the mixture fraction and absolute enthalpy equation. Similarly to the reaction progress variable source term, *ρ* = *ρ Z*-,*Y pv*, *h* and *μ* = *μ Z*-,*Y pv*, *h* are obtained from the employed FGM table. Therefore, these are expressed as functions of the transported scalar quantities *Z*, *Ypv*, and *h*. The source term *Sψ* consists of the source of vapor introduced by the dispersed phase in the transport equation for *Z*, specifically *SZ* = *Sm*. Considering that the mass fraction of fuel is not present in the combination used to define the reaction progress variable (see Equation (15)), as well as no isolated droplet burning model is included in the employed approach, *Sψ* does not contribute (therefore it is set to zero) for the transport of *Ypv*. For the transport equation of *h*, *Sh* is given by

$$\overline{S\_h} = \sum\_{p=1}^{N} \frac{N\_p}{V} \left[ m\_{d,p} \left( \int\_{T\_{nf}}^{T\_d^\ddagger} c\_l dT - \int\_{T\_{nf}}^{T\_d^{t+\Delta t}} c\_l dT \right) + \dot{m}\_{d,p} \left( h\_f - L\_v \right) \right], \tag{4}$$

where *m*˙ *d*,*p* denotes the mass of vapor released by the parcel *p* into the control volume *V*, *md*,*p* is the droplet mass, *Np* the number of real droplets in parcel *p*, *N* the total number of tracked parcels, *cl* the specific sensible heat of liquid, *Tref* a reference temperature (298 K), *Tτ d* the droplet temperature at time step *τ*, *hf* the formations enthalpy, and *Lv* the heat of vaporization. More description about how the phase coupling source terms are computed in the ATF context can be found in [3,21].

Following the ATF method, the flame thickening is performed by means of a dynamic procedure. Accordingly, only the flame region is thickened and no interferences of the ATF with the pre-vaporization zone occur. The flame sensor Ω used by Aschmoneit [35] to simulate partially premixed flames is employed here. The quantity *ω*˙ *pv*,*max Z*-, *h* in Equation (5) is associated to the maximum value of the source term *ω*˙ *pv Z*-,*Y pv*, *h* at the same mixture composition and enthalpy level.

$$
\Omega = \min \left( 1.0, \left( 0.25 \Omega\_{\varepsilon} + 0.75 \times 2.0 \frac{\dot{\omega}\_{pv} \left( \widetilde{Z}, \widetilde{Y}\_{pv}^{-} \widetilde{h} \right)}{\dot{\omega}\_{pv, \max} \left( \widetilde{Z}, \widetilde{h} \right)} \right) \right), \tag{5}
$$

where Ω*c* is the flame sensor proposed by Durand and Polifke [36] as implemented in [29]. In this context, a progress variable in its normalized form is needed, which is defined as *c* = *Ypv*/*Yeq pv*. Herein, *Yeq pv* denotes the equilibrium value of *Ypv* for a specific mixture composition.

*Fluids* **2019**, *4*, 187

In order to facilitate the access of *ω*˙ *pv*,*max* during the simulation time, an extra look-up table is stored with its values mapped on *Z* and *h*. According to the dynamic method, the thickening factor *F* is defined as

$$F = 1 + (F\_{\text{max}} - 1)\Omega.\tag{6}$$

in which *F*max is defined according to an extension of the mixture adaptive thickening proposed in [21] to account for enthalpy variations as

$$F\_{\text{max}} = \max\left(1, \frac{\Delta\_{m\text{sol}}}{\Delta\_{\text{x,max}}}\right) = \max\left(1, \frac{V^{\frac{1}{3}}}{\Delta\_{\text{x,max}}}\right) \tag{7}$$

in which *V* denotes the cell volume and Δ*x*,*max* is the maximum cell size necessary to capture the laminar flame speed with less than 15 % error in one-dimensional simulations [29]. In agreement with Kuenne et al. [29], Δ*x*,*max* ≈ 0.3*δ*0 *l* (*Z*-, *h*). Additionally, two safety factors are also considered to this maximum cell size definition. The first of 1/ √3 is included to account for the worst scenario of flame propagating in cell diagonal direction for three-dimensional cases, the second 1/*α* is used to account for the non-uniformity of cells in *F*. Applying both safety factors, the maximum cell size is defined here as

$$
\Delta\_{\rm x,max} = \frac{0.3}{\alpha \sqrt{3}} \delta\_l^0(\vec{Z}, \vec{h}),
\tag{8}
$$

with *α* = 1.0, since cubic cells are employed in the jet core region (see Section 2.3).

To recover the modifications introduced by the ATF on the turbulence-flame interaction while including the effects of the unresolved FSW, the efficiency function *E*Δ follows the fractal model. It is written as

$$E\_{\Lambda} = \left(\frac{\Delta}{\delta\_l^0}\right)^{\beta},\tag{9}$$

where Δ = *F* · *δl* is the combustion LES filter width, *δ*0 *l* is the laminar flame thickness defined from the temperature gradient, *β* a model exponent, and *δl* an equivalent laminar flame thickness in terms of Gaussian LES filter as proposed in [37] (see more details in the sequel). This is a reasonable approach in the simulations presented here, since we could notice that

$$
\Gamma \frac{u\_\Lambda'}{s\_l} > \max \left( \frac{\Delta}{\delta\_l^0} - 1, 0 \right),
\tag{10}
$$

already for flames presenting a smaller Reynolds number (*Re*) in [3]. Considering the minimum operator in

$$E\_{\Lambda} = \left( 1 + \min\left[ \max\left( \frac{\Lambda}{\delta\_l^0} - 1, 0 \right), \Gamma \frac{u\_{\Lambda}'}{s\_l} \right] \right)^{\beta} \tag{11}$$

this inequality shows that the saturation holds.

Equation (9) is connected with the efficiency function *E*∗ Δ by means of

$$E^\*\_{\Lambda} = 1 + (E\_{\Lambda} - 1)\Omega. \tag{12}$$

Equation (12) is used to enforce that *E*∗ Δ → 1 outside the flame region, since an efficiency function is no longer necessary in this part of the domain. Hence, the term *FE*∗ Δ *μ*/*σψ* in Equation (3) reduces to *μ*/*σψ*, whereas (1 − Ω)*μt*/*σt*,*ψ* becomes *μt*/*σt*,*ψ*. Since *σ* = *σt* = 0.7, both terms are merged in which the molecular viscosity is added with the turbulent one and the effective viscosity *μeff* = *μ* + *μt* is obtained. In contrast to this scenario, within the flame region Ω → 1 and the term (1 − Ω)*μt*/*σt*,*ψ* is canceled, where molecular and turbulent transport in the flame are modeled by *FE*Δ *μ*/*σψ* .

The determination of Δ follows a formulation that can be directly applied for future implementations accounting for the dynamic modeling of *β*. Accordingly, the equivalent laminar flame thickness used to define the LES combustion filter is given by *δl* = *ζδ*0 *l* (in agreement with [37]), in which *ζ* is an adjusting factor to enforce that the efficiency function goes to one for laminar flames. Values of *ζ* are supposed to vary according to the mixture state, however only the variations associated with the composition are accounted for in this work. As a result, the actual formulation of the mixture adaptive thickening which accounts for enthalpy variations is considered to achieve

$$
\Delta = \delta\_l F\_{\text{max}} = \zeta \delta\_l^0 F\_{\text{max}} = \zeta \delta\_l^0 (\tilde{Z}, \tilde{h}) \left[ \max \left( 1.0, a \frac{\sqrt{3}}{0.3} \frac{V^{\frac{1}{3}}}{\delta\_l^0 (\tilde{Z}, \tilde{h})} \right) \right]. \tag{13}
$$

In regions where the cell sizes deliver *F* greater than 1, Equation (13) reduces to

$$
\Delta = \zeta \left( a \frac{\sqrt{3}}{0.3} V^{\frac{1}{3}} \right). \tag{14}
$$

#### 2.1.2. Chemistry

To construct our non-adiabatic FGM table, the strategy proposed by van Oijen and de Goey [11] is applied. Particularly, the reference table described in [16] is employed. It is based on the combination of freely propagating adiabatic flamelets for a mixture composition span of *Z* ∈ [0.050, 0.194] with fresh gas temperature varying from 250 K to 900 K with burner stabilized and extrapolated flamelets. The maximum value of this temperature determines the upper enthalpy value, while its minimum defines the lowest temperature (*Tlow*) allowed in the computational domain. Interpolations are done [38] between the limiting values of *Z* up to pure air and ethanol at 300 K, which corresponds to the ambient temperature of the investigations performed in this work. To compute both kinds of flamelets, the chemical mechanism proposed by Marinov [32] is employed. It represents the oxidation of ethanol in the air by means of 57 species and 379 intermediate reactions.

The inclusion of burner-stabilized and extrapolated flamelets in our table allows a more universal description of the heat losses. Burner-stabilized flamelets mimic premixed flames burning attached to a porous medium with controlled temperature. By keeping the temperature of this medium constant and adjusting the inflow velocity, heat losses can be introduced in the domain resulting in different enthalpy levels for the flame. In this particular step, the temperature is set to the same value as for the unburned gas of the coldest adiabatic flamelet. When the lowest enthalpy levels with the burner stabilized flamelets is reached, the burnt gas is still far away from *Tlow*. To cover the remaining range of physical states, step-wise extrapolations of the thermo-chemical data are done to the lowest enthalpy level. This is computed by assuming the equilibrium mixture composition at *Tlow*. Finally, once a manifold is generated, the look-up table is build up over the controlling variables *Z*, *Ypv*, and *h*. Following the strategy of Ketelheun et al. [12], an additional table is created to facilitate the establishment of the boundary conditions in terms of temperature.

The definition of *Ypv* follows the combination of mass fractions presented in [16]. This was empirically defined in view of a valid progress variable for the full range of enthalpy levels encountered in the employed table. Namely, a monotonically evolving variable from fresh to burnt gases along all flamelets used to construct the employed FGM table. The *Ypv* can be written as

$$Y\_{\rm pv} = \frac{1}{M\_{\rm CO\_2}} Y\_{\rm CO\_2} + \frac{1}{2.5 M\_{\rm H\_2O}} Y\_{\rm H\_2O} + \frac{1}{1.5 M\_{\rm CO}} Y\_{\rm CO\_2} \tag{15}$$

where *Mk* is the molar mass of species *k*. For more details about the definition of a progress variable in FGM context, the reader is referred to [11]).

#### *2.2. Liquid Phase*

The computation of the motion of non-rotating droplets considers only drag and buoyancy forces. Once that the density ratio of liquid ethanol and the gaseous mixture has an order of 103, complementary forces are neglected.

Heat and mass exchanges are computed following the model proposed by Miller et al. [39] with the correction procedure derived in [21] to include effects of the flame thickening onto the dispersed phase. Both phenomena are described by

$$\frac{dT\_p}{dt} = \frac{1}{F} \left[ \frac{f\_2 \text{Nu}}{3Pr} \left( \frac{\theta\_1}{\tau\_p} \right) (T - T\_p) + \left( \frac{L\_V}{c\_l} \right) \frac{\dot{m}\_p}{m\_p} \right] \tag{16}$$

and

$$\frac{dm\_p}{dt} = -\frac{1}{F} \frac{Sh}{3Sc} \left(\frac{m\_p}{\tau\_p}\right) H\_{M\prime} \tag{17}$$

with *T* the temperature, *Nu* the Nusselt number, *f*2 a correction factor due to evaporation computed as in [39], *Pr* the Prandtl number, *θ*1 the ratio of the gaseous and liquid specific sensible heat (*cp*/*cl*), *τp* = *ρpd*2 *p*/18*μ* expresses the particle relaxation time, *m*˙ *p* = *dmp*/*dt*, *Sh* the Sherwood number, and *HM* represents the specific driving potential for mass transfer according to [39]. It is important to highlight that vapor pressures, necessary to determine *HM* are computed here with the Wagner equation as in [40].

For the vast majority of droplets tracked in our simulations, the relaxation time exceeded the scales within the subgrid. As a consequence, SGS turbulent dispersion and micro-mixing modeling for the evaporation process had little effect.

#### *2.3. Experimental Configuration and Numerical Setup*

The flame EtF5 of the Sydney diluted spray flame burner [31] is used as a benchmark for the performed LES. This flame belongs to a series of flames (EtF1, EtF2, EtF5, and EtF7), where the premixed flame mode is dominant. This is an important feature in view of the chosen modeling approach, which has been originally developed for premixed flames. EtF5 was produced by a round jet burner with a diameter D of 10.5 mm fueled with ethanol. The flame stabilization was achieved with a series of pilot flames disposed of in a concentric arrange to the main jet nozzle with an annular diameter of 25 mm. A co-flow surrounded both streams with an annular diameter of 104 mm and a velocity of 4.5 m/s, which corresponds to the air velocity of the wind tunnel where the burner is mounted. The bulk velocities of main jet and pilot were 48 m/s and 11.6 m/s, respectively. At the nozzle exit, a mixture of droplets, air, and fuel vapor was formed with a *Re* of 38,200 and an equivalence ratio *φ* = 0.15 for the gaseous phase, which is found below the lower flammability limit (*φlow*).

The computational domain is constructed with a Cartesian grid in order to keep the cell growth rate smooth and to avoid large values of the factors *F* and *α* (see Equations (3) and (14), respectively). Specifically, mean values of *F* are about 7.0 but no bigger than approximately 9.0, while *α* = 1.0 due to the cubic cell arrangement in the jet core region. When compared to the "O-grid"arrangement (as employed in [3]), which would be an alternative to run in the used block-structured code, the Cartesian arrangement avoids the high cell diagonal growth rate. Figure 1 presents a sketch of the employed grid, which is composed of approximately 5.4 million cells distributed over 468 blocks. In this figure black lines represent block limits, whereas gray lines are cell sides. The dashed orange box indicates the jet core region, in which cells have uniform size with a side of 0.75 mm. A contour plot of the instantaneous value of *ω*˙ *pv* is included to show that the main reaction zone is entirely found within the jet core. The computational domain starts at the exit of the experimental burner and extends up to 70 D (735 mm). The most upstream part of this domain has a diameter of 10 D that linearly grows up its maximum with 20 D. The employed grid can be seen as an improvement of the coarse grid used in [3], since it preserves similar characteristics of the previous one (e.g., longitudinal and

radial cell sizes are the same, however with higher homogeneity here). In [3] a comprehensive grid study has been conducted for reactive and non-reactive flows.

Regarding the boundary conditions, attention is given to the description of the pilot flame. This flame is fueled by a mixture of acetylene and hydrogen at the hydrogen to carbon proportion of the stoichiometric reaction of ethanol. Under adiabatic assumptions, such a boundary condition can be easily represented by adjusting the mixture composition to the stoichiometry. This procedure is previously adopted for simulations presented in [3] and it is commonly employed by diverse research groups (e.g., [7,33,41]) which also used data from the Sydney diluted spray flame burner for modeling validation. Concerning a non-adiabatic formulation, the introduction of the enthalpy transport allows a more accurate representation of the pilot flame. It is important to notice that, despite the same proportion of hydrogen and carbon elements, the presence of the hydroxyl (OH) into the ethanol molecules does not allow the mixture of acetylene and hydrogen to produce the same composition of products. This characteristic becomes clear comparing the stoichiometric reaction of ethanol and the corresponding mixture of acetylene and hydrogen with air

$$\rm C\_2H\_5OH + 3\ (O\_2 + 3.76N\_2) \longrightarrow 2CO\_2 + 3H\_2O + 11.28N\_2 \,\tag{18}$$

$$\rm C\_2H\_2 + 2H\_2 + 3.5 \left(O\_2 + 3.76 N\_2\right) \longrightarrow 2CO\_2 + 3H\_2O + 13.06 N\_2. \tag{19}$$

Since the main function of a pilot flame consists in the production of hot products to stabilize the main flame, the mixture fraction is adjusted to match the products ratio generated by Equation (19). Namely, the ratio formed by the sum *Y*CO2 + *Y*H2O considering only products. As a result, the following reaction is used to define the mixture composition of the pilot flame for the non-adiabatic simulations performed here

$$0.89\text{C}\_{2}\text{H}\_{5}\text{OH} + 3\left(\text{O}\_{2} + 3.76\text{N}\_{2}\right) \longrightarrow 1.78\text{CO}\_{2} + 2.67\text{H}\_{2}\text{O} + 11.28\text{N}\_{2} + 0.33\text{O}\_{2} \,. \tag{20}$$

**Figure 1.** Computational domain. The dashed orange box indicates the jet core.

Once that the absolute enthalpy of the reaction given by Equation (20) does not match with that obtained from Equation (19), the reactants are lightly pre-heated in order to pair up with the desired value. It is important to highlight that with a tabulation method where heat losses are exclusively linked to the evaporation process as in [13,22], such a modeling feature would not be realizable. The relevance of the specification for this boundary condition has been analyzed in [3] for the flame EtF2.

The turbulent boundary condition follows the synthetic turbulence generator proposed by Klein et al. [42]. To account for the pre-vaporized mass of fuel exiting the injection's tube, the mixture fraction receives the explicit value of 0.0165 (corresponding to *φ*over = 0.15). Effects of the heat transfer throughout the injection tube are considered. For that, the mixture composed by dispersed fuel droplets and air entering at 300 K into the burner exchanges heat with the injections tube

wall at 300 K, which results in a mixture exiting the nozzle with 298 K. This arises from our assumption that the mixture of vapor and air with *φ* = 0.15 and *Re* = 38,200 results from the evaporation of liquid in an air stream where both were initially at 300 K. This mixture was then heated up by tube walls at 300 K assuming *Re* = 38,200. Considering the remaining boundary conditions for the enthalpy, walls are treated as adiabatic and incoming streams of air are kept at 300 K.

With respect to the boundary conditions of the diluted phase (i.e., droplet size distributions, conditioned velocities on size, injected mass, and injection position), these have been extracted from experimental measurements presented in [31]. An exception occured for the radial component of droplets velocities, for which the procedure suggested by De and Kim [6] was employed here. More details about this procedure are presented in Section 3.1. The initial droplet temperature has been set to 298 K, namely droplets and the pre-vaporized mass of fuel exiting the injection's tube had the same temperature.

Calculations for this research were conducted with the coupled version of the academic software FASTEST (see [3]). Different from previous works, some improvements in numerical methods and algorithms are included. For instance, the ordinary differential equations of the evaporation modeling are integrated in time here with the fifth-order Runge–Kutta scheme of [43].

Simulations are initialized from reactive single-phase computations obtained after 0.055 s, which coincides with one flow-through-time (FTT). The FTT of this configuration corresponds to the time necessary for the gas flow to fill a volume defined from the droplets injection position to 30 D downstream the burner exit, where the last measurement plane is set. Namely, a cylinder with a diameter of 10 D (to account for the burner coflow) and a length of 30 D. Thereafter, parcels start to be injected and a stable number of them is achieved after approximately 1 FTT. Once a stable number of parcels is reached, simulations run for an additional FTT until the averaging begins. Statistics (mean and rms) are collected for 2 FTT for all cases, which amounts circa 35,300 core-hours on 80 CPUs and corresponds to a reduced computational time (RCT) of ≈ 7.11 × 10−8 s [RCT = elapsed time × number of CPU/(number of cells × number of iterations)].

To perform the time integration, a second-order, three steps Runge–Kutta scheme (see details in [44]) has been employed. The time-step is maintained constant at 2.4 × 10−6 s. Particles are coupled at each three-time steps. As a result, particles were tracked in a coupling time step of 7.2 × 10−6 s. Despite this discrepancy, the corresponding Courant–Friedrichs–Lewy (CFL) number is kept below 1 respectively for droplets and the gas phase. This coupling strategy alleviates the computational costs of the particle tracking, which amounts in our calculations approximately 50% of the entire computational time.

#### **3. Results**

The evaluation of the proposed modeling strategies starts with the establishment of a reference scenario. This corresponds to simulation performed with a fixed exponent *β* (i.e., *β* = 0.5) of Equation (9) and the inclusion of evaporative cooling effects. The value 0.5 resembles the one used in the first investigations of Charlette et al. [28] and is a typical value assumed by different authors (e.g., [9,20,29,44]), therefore it is chosen as a reference here. Altogether three cases are studied, besides the mentioned reference case, one neglecting effects of the evaporative cooling and another considering a parametric variation of the model parameter *β* while including EC effects. In this last case, the value of *β* is defined as 0.2.

The influence of the evaporative cooling onto the flow is expressed in terms of the heat necessary to vaporize an amount of liquid (i.e., that associated with the latent heat of vaporization) and the heat-up of droplets. In order to investigate such an influence, EC effects are neglected by removing the

terms associated with the latent heat of vaporization and the heat-up of droplets from Equation (4). As a result, only the term including the vapor formation's enthalpy is preserved

$$\overline{S\_{\hbar}} = \sum\_{p=1}^{N} \frac{N\_p}{V} \left( \dot{m}\_{d,p} \mathbf{h}\_f \right). \tag{21}$$

The source term presented in Equation (21) is necessary for our approach since the energy equation is expressed in terms of the absolute enthalpy. In cases where sensible enthalpy is transported *Sh* = 0.

The remainder of this section is divided into two parts. In the first one, simulations results are compared with available experimental data following a validation process. In the second part, analysis of the simulated flames characteristics are conducted in a qualitative fashion.

#### *3.1. Validation of Modeling Strategies*

Radial profiles of mean and RMS values of the resolved longitudinal and radial components of droplet velocities are compared with experimental data in Figure 2. As a Cartesian coordinate system is employed the radial direction is represented by the coordinate Y, while the longitudinal by X. Focusing on the first column, one can observe that the proposed strategy is able to recover the mean values of droplets longitudinal velocities reasonably well. Deviations from experimental data are more pronounced downstream of the positions 25 X/D for all cases. Comparing the different simulated cases, further deviations are noticed between positions 5 X/D to 15 X/D. Similar to the observations in [3], it can be clearly noticed that the FSW modeling affects the estimation of mean droplet velocities. Specifically, when *β* = 0.2, results better agree with experimental data. Considering this aspect and the fact that such deviations are located in a radial region far from the jet centerline (i.e., Y/D > 0.5), it is possible to state that deviations occur in regions where droplets interact with reaction zones. Nevertheless, due to the strong coupling among the diverse modeling approaches it is not straightforward to draw an exact explanation for it. Influences of velocity fluctuations, thickening factor, and mixture stratification are also supposed to contribute to it, as pointed out throughout the discussions presented in this section. Concerning both cases computed with the same value for the model parameter *β*, marginal deviations are observed for Ud between them. Close to the jet core (i.e., Y/D < 0.5) both models deliver the same results. Influences are supposed to be caused by the variations of the flame speed caused by heat losses in the flame front as marginal deviations are noticed when Y/D > 0.5 between positions 5 X/D to 15 X/D.

Considering the second column of Figure 2, it is important to highlight that despite the absence of an SGS turbulent dispersion model, droplet velocity fluctuations agree reasonably well with the experimental data. The discrepancies noticed for the positions 10–15 X/D are in great parts assigned to the flame-dispersed phase interaction, which directly influences the droplet tracking. This aspect is supposed to justify the inflection observed when 0.5 < Y/D < 1.0 in the positions 10 X/D and 15 X/D. By comparing the different cases, it is again noticed that the flame wrinkling model also contributes to deviations.

According to the radial profiles of Vd, the experimental values are different from zero at Y/D ≈ 0.0, i.e., the jet centerline. As already discussed in [3], the measurements of this velocity component expose a problem on the experimental configuration. For vertical jets, the time averaged radial velocity along the centerline is always zero. However, this is not observed for all comparison planes, including the injection position. This result possibly indicates that the experimental spray is not strictly vertical. Due to this characteristic, the procedure suggested by De and Kim [6] is used here to define the radial velocity of droplets. Namely, the radial velocities of the carrier phase are used to define the starting conditions for droplets. This option justifies the disparity observed in the boundary condition's plane (Figure 2) for Vd. Altogether, simulations results follow a similar trend as the measured data. The existence of non zero values of Vd in the proximity of the jet centerline motivates the addition of a constant offset for all points (for which a zero velocity is attained at the centerline) of the experimental

data. Considering them as a reference, a better agreement is noticed for the computed sprays and even better for the case with *β* = 0.2. However, because of the anomalies perceived in the experimental measurements (see e.g., spray jet anomaly in [45,46]), further discussions about the radial components of mean droplets' velocity are difficult to be conducted. The same is valid for the last column of Figure 2, where the RMS values of the fluctuations of such velocity components are exhibited. These data are rather used here to give support to the discussions about the different simulated cases.

**Figure 2.** Radial profiles of mean (capital letter) and rms of the fluctuations (lower case letters) of the longitudinal (Ud and u d) and radial (Vd and v d) components of droplets velocity. Squares-original experimental data, triangles-offset, blue line-evaporative cooling (EC) with *β* = 0.5, red line-no EC with *β* = 0.5, black line-EC with *β* = 0.2.

An important quantity to observe the mass transfer between phases is the liquid volumetric flux (LVF), whose radial profiles at different longitudinal positions are presented in the first column of Figure 3. Both cases computed with *β* = 0.5 show lower values for LVF downstream of the position 10 X/D when compared with the case with *β* = 0.2. This is a straightforward outcome caused by the increased flame surface wrinkling obtained with *β* = 0.5. As the total FSW increases, the consumption flame speed also increases and the overall flame becomes shorter in the longitudinal direction while broader in the radial direction. Consequently, the spray inner core faces a gas flow with higher temperatures than the case computed with *β* = 0.2 (see Figure 4). Such higher temperatures increase the evaporation rate and reduce the corresponding LVF. The hot mixture that intensifies the evaporation process is an outcome of the combustion process, which is found in reaction and oxidation zones. As droplets cross a reaction zone, interference with the reaction progress is expected to occur. The amount of droplets interacting with reaction zones has not been quantified in the present work. Nevertheless, studies conducted in [9,16,21] show that FGM tables constructed with premixed flamelets are able to recover the effects of variations of the mixture fraction throughout a flame with

good accuracy.The effects of the temperature reduction of the jet inner core and the consequent increase of the simulated LVF are also noticed between both cases computed with *β* = 0.5. The inclusion of EC effects indeed decreases the mass transfer between phases as observed by the higher values of LVF downstream of the position 15 X/D. Sacomano Filho et al. [16] indicated by means of flames propagating in droplet mists that the evaporative cooling effects also contribute to the reduction of the laminar flame speed. Bearing this aspect in mind, the increase of LVF and decrease of gas temperature (see Figure 4) when including EC are combined features from the reduction of the inner core temperature and the flame speed. However, the contribution given by the inclusion of EC over the LVF is not so pronounced as that given by the reduced value of the model parameter *β*. It is important to highlight that the case computed with *β* = 0.2 agrees much better with the experimental data than the other cases. Again, this gives an indication of the relevance of the turbulence-flame interaction modeling.

Considering the droplet size distribution, reasonably good agreement with the experimental data can be perceived for the Sauter mean diameter (SMD-*D*32) in the inner part of the jet (i.e., Y/D < 0.5) and up to 25 X/D. In the farther downstream regions, the SMD becomes slightly smaller than the measured values in all radial positions. A probable cause for it is the penetration of droplets in higher temperature zones, which also agrees with the observed evolution of LVF radial profiles. The interaction of the dispersed phase with reaction zones may justify the deviations noticed for the SMD for radial positions far from the jet centerline (i.e., Y/D > 0.5) from 5 X/D to 20 X/D.

**Figure 3.** Radial profiles of mean liquid volumetric fluxes, the Sauter mean diameter *D*32, as well as mean (U) and the rms of the fluctuations (u') of the longitudinal component of the gas velocity. Marks-experiment, blue line-EC with *β* = 0.5, red line-no EC with *β* = 0.5, black line-EC with *β* = 0.2.

In the last two columns in Figure 3, radial profiles of the mean and the RMS values of the resolved longitudinal component of the carrier phase velocity are shown. Experimental data acquired for droplets smaller than 10 μm are used as a reference. In fact, both quantities are not strictly comparable, but the experimental data is used here to guide the subsequent discussions. Many aspects observed for the dispersed phase in Figure 2 are noticed in the gas flow. As for the spray, deviations occur close to reaction zones and in a similar fashion among the different cases. The correct prediction of the flame topology is supposed to contribute to these observed deviations. Concerning the velocity fluctuations, it is important to bear in mind that comparisons performed here do not take into account the contribution of SGS quantities since these are not available. Veynante and Knikker [47] showed that this lack of information does not modify mean but RMS values. As a consequence, the inclusion of these SGS contributions should modify the scenario regarding the rms values.

**Figure 4.** Radial profiles of time averaged and rms of the fluctuation of the equivalence ratio and the temperature of the carrier phase. Marks-experiment, blue line-EC with *β* = 0.5, red line-no EC with *β* = 0.5, black line-EC with *β* = 0.2.

To conclude the validation process, radial profiles of the mean gas temperature are compared with the available experimental data in Figure 4. Results obtained with *β* = 0.2 better agree with experiments when compared with those from the remaining cases. Effects of the evaporative cooling are noticed in all positions, but these are more pronounced as the distance from the burner nozzle increases. This is expected since the heat removed from the gas phase increases due to the cumulative impact of the evaporation. However, the effects of FSW modeling are more pronounced than the sole consideration of EC while keeping *β* = 0.5. Specifically, deviations between cases computed with *β* = 0.2 and *β* = 0.5 are more pronounced at 10 X/D and 20 X/D. At 30 X/D, similar results can be seen between cases including EC. This overall behavior regarding FSW modeling is explained considering the flame structures presented in the next section. The FSW modeling is strictly connected with the source term *ω*˙ *pv* (see Equation (3)). According to the qualitative contour plots presented in Figure 8, source terms of the progress variable are predominantly found upstream 25 X/D. Hence, the effects of the FSW modeling are more pronounced upstream of this region. Downstream of it, unburned combustion products are mixed together with fuel vapor released from droplets that are able to cross this main reaction zone. The resulting mixture burns further in a diffusion flame mode, which is less affected by the FSW model. Therefore, different approaches for the computation of *β* deliver similar results at 30 X/D.

Analyzing the temperature results with more details, simulations computed with *β* = 0.5 show broader regions with high-temperature values than the results obtained with *β* = 0.2 at 10 X/D. In fact, simulations computed with *β* = 0.2 agree quite good with the experimental data at this position. Contributions from the inclusion of EC are seen at the region where the temperature peaks, but these are not high. As previously mentioned, effects of the EC are better seen at 20 X/D and 30 X/D throughout the coordinate Y. At 20 X/D, the ascending part of the temperature profile in Y direction for the case in which *β* = 0.2 agrees better with experiments than the other cases. The performance of this case at these two positions explain the more accurate prediction of LVF in Figure 3.

#### *3.2. Qualitative Analysis of Simulated Cases*

To facilitate the interpretation of the mixture composition, Figure 4 presents time-averaged and RMS values of the equivalence ratio (obtained from the resolved fields of mixture fraction as in [48]). The time averaged values of *φ* point out how strong the FSW model and the inclusion of EC affect the flame structure. Particularly, this variable is connected with the LVF seen in Figure 3, as the vapor released from droplets directly affects the mixture composition. Overall deviations among cases are more pronounced between different definitions of *β* than the two treatments reserved for the EC. Deviations between both cases computed with *β* = 0.5 occur at the position 20 X/D. Observing the time-averaged Lagrangian source terms of mass *Sm* and enthalpy *Sh* respectively in Figures 5 and 6, we see that such effect is a direct outcome of the inclusion of EC. Both cases show similar *Sm* profiles at 20 X/D and upstream this position, while clearly different profiles of *Sh*. Therefore, the inclusion of EC effects clearly interferes with the mixture composition, as it is a product of higher transported temperatures of the gas phase (see Figure 4) through the higher droplet evaporation rates observed in terms of LVF profiles in Figure 3. With respect to the FSW modeling, the case computed with *β* = 0.2 already interferes with the mixture composition at 10 X/D. As seen in Figure 4 and reinforced in Figures 5 and 6, the main contribution to this aspect stems from the region where the reaction zone is found (Figure 8 points out the position of reaction zones by means of the time-averaged *ω*˙ *pv*), i.e., 0.5 < Y/D < 1.0. Figures 5 and 6 indicate that small impact of this modeling approach is found in the jet core region, namely Y/D < 0.5. As the flame evolves, more interactions with the spray at the jet inner core occur at a downstream position as seen through the Lagrangian source terms at 20 X/D and 30 X/D. The topology of the flame is directly related to this observation. As seen in Figures 7 and 8 upstream of 20 X/D both flames computed with *β* = 0.5 depict similar topology. Downstream of this position small flame interactions occur with the spray, whereas droplets that are able to cross the flame tip interact with combustion products coming out of the main flame. The Lagrangian source term profiles indicate that differences between models are also happening at 30 X/D. Nevertheless, these are less intense as no significant impact is noticed in time-averaged profiles of *φ* in Figure 4 at this position.

**Figure 5.** Radial profiles of the Lagrangian source term of mass *Sm*. Blue line-EC with *β* = 0.5, red line-no EC with *β* = 0.5, black line-EC with *β* = 0.2.

**Figure 6.** Radial profiles of the Lagrangian source term of enthalpy *Sh*. Blue line-EC with *β* = 0.5, red line-no EC with *β* = 0.5, black line-EC with *β* = 0.2.

Considering the RMS of *φ*, higher values are found in the region that coincides with the descending part of the temperature profiles. Such an outcome can be seen as a combined effect of the interaction of droplets with high-temperature combustion products and radial flow coming out of the main jet core as seen through the radial profiles of Vd in Figure 2. Deviations among cases are analogous to those found in the time-averaged profiles of the gas temperature. These are significant between FSW modeling approaches at 20 X/D, whereas between EC treatments at 30 X/D.

**Figure 7.** Contour plots of instantaneous values of the source term *ω*˙ *pv* [kmol/m3s] in a plane of size 10 D × 40 D at Z/D = 0 for: (**a**) EC with *β* = 0.5; (**b**) no EC with *β* = 0.5; (**c**) EC with *β* = 0.2.

Regarding the RMS of gas temperature, similar trends are noticed among cases, except at 20 X/D. At this position, higher values are found for the case with *β* = 0.2 in the proximities of the jet centerline. This characteristic can be explained considering the instantaneous plots of the source term *ω*˙ *pv* presented in Figure 7. Different from the cases computed with *β* = 0.5, the case with *β* = 0.2 deliver a longer flame in which pockets are frequently formed and released close to the flame tip. This occurs close to 20 X/D, where high fluctuations appear for *β* = 0.2 in Figure 4.

**Figure 8.** Contour plots of the time averaged source term *ω*˙ *pv* [kmol/m3s] in a plane of size 10 D × 40 D at Z/D = 0 for: (**a**) EC with *β* = 0.5; (**b**) no EC with *β* = 0.5; (**c**) EC with *β* = 0.2.

An additional aspect that can be extracted from the contour plots of *ω*˙ *pv* presented in Figures 7 and 8 refers to the mixture stratification. Such a feature can be observed in terms of the different intensities of *ω*˙ *pv* throughout the flame. The stratification is also noticed in *φ* radial profiles of Figure 4. The variation of the mixture composition within the computational domain is a direct outcome of the evaporation of liquid droplets. Due to the free jet configuration, the amount of released vapor increases as the spray penetrates into the domain. The opposite behavior is noticed in the liquid phase by means of the LVF profiles in Figure 3.

Contour plots of time-averaged values of the gas temperature give a more complete overview of this quantity for the different cases in Figure 9 than the radial profiles shown in Figure 4. The higher penetration length of the flame computed with *β* = 0.2 is also observed here considering the penetration of the cold jet core. Additionally, the higher temperatures observed for the case where EC effects are neglected (using the other case computed with *β*=0.5 as a reference) are also noticed. The extension of the slice plane used in Figure 9 allows the visualization of another import aspect. Considering the jet centerline (i.e., Y/D = 0.0), as the distance from the nozzle increases, the temperature initially grows from the cold-core level up to a first pick which coincides with the main flame tip position. Downstream of it, the temperature slightly reduces and grows up again for positions above 30 X/D. Particularly, at the region located between 20 X/D and 35 X/D for cases computed with *β* = 0.5, whereas between 24 X/D and 38 X/D for the one computed with *β* = 0.2, a local cold-core is identified. A similar structure is also perceived, but with much less contrast, in Figure 8 in terms of time-averaged values of the source term *ω*˙ *pv*. The appearance of it is a result of the combustion of unburnt products mixed with the vapor released from droplets that are able to cross the main reaction zone. In Figure 7 isolated flame structures are found downstream of the main reaction zone. Besides being quite sparse these structures are preferentially located on the jet borders, where the mixture of products and vapor diffuses with the co-flowing air stream. Accordingly, the occurrence of these secondary reactions resemble processes controlled by the flow mixing. This characteristic aligns with the definition of diffusion flames, namely mixing controlled combustion reactions [49]. It is important to highlight that such a multi-regime flame differs from those observed in [10], where the main premixed flame is enveloped by a diffusion-reaction layer. Such an overall behavior is supposed to do not significantly

affect the conducted analysis in the present work since the focus has been given to positions upstream of the second cold-core region. Nevertheless, the inclusion of the method recently proposed in [40] to improve the evaporation modeling of droplets interacting with combustion products, as well as the adoption of a multi-regime combustion model (e.g., [50]), may improve the characterization of this post main flame phenomena which may not be reserved to this configuration.

**Figure 9.** Contour plots of the time averaged values of the gas phase temperature [K] in a plane of size 10 D × 40 D at Z/D = 0 for: (**a**) EC with *β* = 0.5; (**b**) no EC with *β* = 0.5; (**c**) EC with *β* = 0.2.

#### **4. Conclusions**

Investigations of the effects caused by the consistent consideration of the evaporative cooling and the turbulence flame interaction demonstrated that both phenomena do interfere not only with the estimation of the flame structure, but also with the prediction of spray properties. As previously pointed out in [3] for a similar flame but with lower Reynolds number, such an outcome reveals the strong interaction among turbulence, flame, and dispersed phase. The necessity of an appropriate computation of each of these phenomena is again reinforced in the present study.

The combination of evaporative cooling effects through the FGM method and the flame surface wrinkling model considering *β* = 0.2 demonstrated to have the best performance in the validation process. The corresponding simulated case showed the best overall agreement with available experimental data and could recover the main phenomena observed experimentally. For the validation process, the experimentally measured flame EtF5 of the Sydney diluted spray flame burner [31] was used as a benchmark. The other two cases have been used to allow a parametric study of the influence of the EC effects and FSW modeling. The raised hypothesis in that EC effects could be quite important [16] was accessed by comparing two cases, with and without consideration of EC effects at a rather standard of the FSW model parameter *β*. Although EC effects were shown to be important for the correct flame prediction, the influence of the turbulence-flame interaction modeling demonstrated to be more pronounced.

Thanks to the phenomena explicitness that LES and ATF approaches allow, the main features observed in the simulation results could be recalled by analyzing further computed quantities. Such a clear interpretation of different phenomena could not be easily done if the proposed modeling extensions to account for heat losses were not consistently done.

Places for improvements are found as well as pointed out throughout the result presentation and discussion. It is important to highlight that the proposed modeling strategy is still limited since part of the underlying phenomena to turbulent spray combustion is not accounted for. Indeed, some effects were not comprehensively studied or even understood to allow to infer on their relevance in the simulations. Two of them are the inclusion of actual mixture composition on droplet evaporation (see e.g., [40]) and effects of differential diffusion in spray flames (see e.g., [16]). Additional effects that may improve the predicting ability of the proposed modeling strategy is, on the one hand, the inclusion of models capable to represent multi-regime combustion, as studied by Knudsen et al. [13] and Hu and Kurose [50], and the incorporation of single droplet burning effects (see [51]) on the other. Nevertheless, it worths to be noticed that the proposed methods following the inclusion of the evaporative cooling effects in the chemistry and the selected FSW modeling converge to an optimum combination of strategies which better approaches the experimental observations. Models able to represent the above-listed phenomena are desired and should improve the capability of the presented module.

**Author Contributions:** Conceptualization and formal analysis, F.L.S.F.; methodology, F.L.S.F., L.D. and A.H.; software, F.L.S.F. and A.H.; resources, G.C.K.F., and A.S.; writing—original draft preparation, F.L.S.F. and L.D.; supervision, G.C.K.F., and A.S.; funding acquisition, F.L.S.F., G.C.K.F., and A.S.

**Funding:** This research was funded by São Paulo Research Foundation (FAPESP) grant number 2017/06815-7, as well as by the German Research Council (DFG) grant number TRR 150.

**Acknowledgments:** We gratefully acknowledge the computational resources by Lichtenberg HPC of the TU Darmstadt and A. Masri to make the experimental data available.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **References**

c 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).


**4.3. Identification of AMY and GLL1 using RP-HPLC for protein**

Nowadays, identification of proteins usually employs techniques with sophisticated and delicate instrumentation e.g. peptide-mass finger print mass spectrometry (PMF-MS) (24). Unfortunately, this technique is rather pricey and furthermore, requires convenient access to the protein database and to the amino acid sequence of the protein. Particularly in the lack of the latter, which is a very common situation in the early stage of protein works, the separation of AMY from GLL1 was confirmed by a more simple and robust technique, exploiting the use of an RP-HPLC system (25). The analysis using RP-HPLC is based on the fragmentation profile of a protein after proteolytic digestion. Proteolytic digestion of AMY was expected to result in fragmentation profile that differs to that of GLL1, as reflected from their chromatographic profiles.

276 Chromatography – The Most Versatile Method of Chemical Analysis

isolation of purified AMY and GLL1 was unanimous.

(ml)

Volume

glucoamylase.

chromatography (RP-HPLC) column. In the Spheron 300 LC column, α-amylase (strain KZ) was eluted at an ammonium sulfate concentration of 10% whilst glucoamylase was at 5%. The separation profile, however, is strikingly similar to that of AMY and GLL1 with butyl-Toyopearl column. However, separation of AMY from GLL1 in the butyl-Toyopearl column has better resolution, suggesting variations in the surface properties of *S. fibuligera* αamylase and glucoamylase from different strains, or the butyl-Toyopearl resin serves for better separation because there is no interaction between the proteins and the supporting matrix. Nevertheless, these results emphasize the power of HIC to separate α-amylase from

**4.2. Subsequent purification of AMY or GLL1 in chromatography columns**

Total activity (Units)

**Table 1.** The final purification scheme of AMY. AS stands for ammonium sulfate.

**4.3. Identification of AMY and GLL1 using RP-HPLC for protein**

Crude 2000 2640000 16500 160 100 25% AS 2200 2244000 13046 172 85 Butyl Toyopearl<PHONE_NUMBER> 753 1438 41 95% AS 35 765600 213 3595 29 DEAE-Toyopearl 60 580800 100 5808 22

Nowadays, identification of proteins usually employs techniques with sophisticated and delicate instrumentation e.g. peptide-mass finger print mass spectrometry (PMF-MS) (24). Unfortunately, this technique is rather pricey and furthermore, requires convenient access to the protein database and to the amino acid sequence of the protein. Particularly in the lack of the latter, which is a very common situation in the early stage of protein works, the separation of AMY from GLL1 was confirmed by a more simple and robust technique, exploiting the use of an RP-HPLC system (25). The analysis using RP-HPLC is based on the fragmentation profile of a protein after proteolytic digestion. Proteolytic digestion of AMY

Total protein (mg)

Specific activity (U/mg)

Yield (%)

Intended for their characterization, AMY and GLL1 were independently collected and subjected to subsequent purification with DEAE-Toyopearl 650M AEX column. An SDS PAGE analysis showed that these purification steps resulted in pure AMY (5) and GLL1 (16). The final purification scheme of AMY is summarized in Table 1. Further, whilst the presence of two and three types of α-amylase and glucoamylase, respectively, were reported upon purification of these amylolytic enzymes from strain KZ by a Mono-Q anion exchanger column (22), AMY and GLL1 appeared to be the only amylolytic enzyme species from strain R64. Nevertheless, additional analysis was performed to confirm that the The analysis of AMY and GLL1 by RP-HPLC was carried out essentially following Soedjanaatmadja and co-workers on the identification of pseudo-hevein from the latex of rubber tree (26). Firstly, purified AMY was incubated for four hours at 37oC with chymotrypsin (EC. <IP_ADDRESS>), at an AMY to chymotrypsin mass ratio of 100:1, in 200 mM ammonium bicarbonate buffer, pH 8.0. The reaction was stopped by an addition of 10 mM hydrochloric acid, to lower the pH of the solution. Preparation of GLL1 sample was done in the same way, including the substrate to chymotrypsin mass ratio.

The reaction mixture was immediately applied to an RP-HPLC nucleosil 10 C18 column (30 x 0.45 cm) and the separation was performed for 60 minutes, at a flow rate of 1 ml/min, using 0-70% acetonitrile gradient in 0.1% trifluoroacetic acid (TFA). Elution of fragments was monitored at λ 214 nm, which is specific for detection of peptide bonds. The fragments are eluted according to their hydrophobicity, where more hydrophobic fragments are retained longer in the RP-HPLC column.

Chymotrypsin cleavage takes place on peptide bonds at the C-terminal part of preferably tyrosine, phenylalanine, tryptophan, and leucine residues, and with (much) less extent of methionine, valine, isoleucine, histidine, glycine and alanine residues (27). Due to this broad specificity, the cleavage may result in many fragments varying in size and hydrophobicity, as observed in the chromatographic profiles of proteolytically digested AMY and GLL1 (Fig. 3).

**Figure 3.** Chromatographic profile of fragments from proteolytic digestion AMY (A) or GLL1 (B) as monitored at λ 214 nm.

The fragmentation profile of AMY was significantly different from that of GLL1. Six sets of resolved peaks were recovered in both cases (Fig. 3) but their retention time, intensities, and

peak distribution were different. Intensity of the peaks suggested much less materials were recovered from GLL1 digested samples than from AMY. This situation is likely contributed from highly hydrophobic or negative charged fragments that were not eluted from the RP-HPLC column due to their poor solubility in the solvent used (28). The amino acid sequences of AMY and GLL1 (see 3.5) showed that the charged amino acid distribution in the amino acid sequence of GLL1 (100 charged amino acid residues out of 494, ~20.2%) is higher than that of AMY (86 out of 468, ~18.4%), thus AMY contains more hydrophobic and non-charged amino acids. Unfortunately, correlation between the result from RP-HPLC and hydrophobic amino acid distribution cannot firmly be established because the nature of chymotrypsin digestion was unclear. The results may indeed indicate that hydrophobic amino acids in GLL1 are likely more clustered than in AMY, resulting in highly hydrophobic peptide fragments in GLL1. This hypothesis may be related to the structure of GLL1 (Fig. 8), which consists of one globular molecule that opposes two separable domains of AMY. However, possibility for the existence of highly negative charged peptide fragments from GLL1 can also not be excluded. Nevertheless, of the eluted fragments (Fig. 3), majority of AMY fragments are localized at peaks 5 and 6 whilst GLL1 are at peaks 3 and 4, suggesting different fragmentation had occurred. Thus, the RP-HPLC profiles of AMY and GLL1 indicated successful separation of the two enzymes.

Chromatography as the Major Tool in the Identification and

the Structure-Function Relationship Study of Amylolytic Enzymes from Saccharomycopsis Fibuligera R64 279

mobile phase, which was the mixture of butanol : ethanol : water (at a ratio of 5:5:3, v/v/v). The plate was retrieved from the TLC separation chamber after the mobile phase migration reached three quarters of the length of the plate and then immediately short-submersed (few seconds) in a mixture of water : methanol : sulphuric acid (at a ratio of 45:45:10, v/v/v). The plate was then heated at 120oC for 15 minutes on a hot plate for visualization of the sugars.

The TLC profile (Fig. 4) shows that the product of starch hydrolysis by GLL1 was solely monosaccharide i.e. glucose (G1). The intensity of the G1 spot was increasing from 5 to 45 minutes of incubation, indicating that glucose was produced over time. The spots at the sample application points were from the un-hydrolyzed soluble starch substrate. Apparently, most of the starch molecules were hydrolyzed within 15-45 minutes of incubation. Most importantly, no higher oligomeric sugar forms were detected, even during the first 5-10 minutes of incubation, suggesting the absence of α-amylase. This result provided the undisputed proof that GLL1 had successfully been separated from AMY despite of α-amylase activity was detected in the GLL1 fraction. Thus, the detected α-

**Figure 4.** Analysis of starch hydrolysis products of GLL1 on a TLC plate (16). Lane 1-4 is hydrolysis products after 5, 10, 15, and 45 minutes, respectively. Lane 5 is the mono- (G1), di-(G2), and oligo-

Since the purity of AMY and GLL1 was firmly established, each enzyme could now be characterized. AMY was found active at a pH range of 5.0-7.5 and a temperature range of 30-60oC, with an optimum at 5.5 and 50oC, respectively (5). GLL1 was also found active a pH range of 4.6-6.8 and a temperature range of 30-65oC, with an optimum at 5.6-6.4 and 50oC, respectively (16). These results suggest that the two enzymes are active at similar conditions.

Recently, the amino acid sequences of AMY (GenBank accession code HQ172905) and GLL1 (HQ415729) were successfully elucidated as derived from their chromosomal DNA (16, 33).

Moreover, their characters are similar to most other *S. fibuligera* amylases.

amylase activity was not originated from AMY.

saccharide (G3-G7) markers.

**4.5. Properties of AMY and GLL1**

## **4.4. Analysis of starch hydrolysis products to discriminate AMY and GLL1**

Glucoamylase hydrolyzes starch at the non-reducing end of amylose or amylopectin to result in glucose therefore the enzyme activity assay is normally based on the release of reducing sugar (21, 29). Unfortunately, α-amylase random digestion of starch may also result in reducing sugar i.e. maltose (30) therefore detection of glucoamylase activity in an α-amylase preparation can be anticipated. Glucoamylase can also, at lesser extent, act on 1,6 α-glycosidic bond of amylopectin (31), although its efficiency diverse greatly depending on the source organisms. Hydrolysis of the 1,6-α- bond may result in a less integrated starch molecule, hampering the formation of iodine-starch complex, which is the basis of standard α-amylase activity assay (20). Therefore, cross-detection of the two enzyme activity assays is inevitable. This phenomenon was notorious upon separation of AMY from GLL1 at pH 5.8 (5), where α-amylase activity was detected in GLL1 fractions.

One approach to discriminate α-amylase and glucoamylase is *via* the evaluation of their hydrolysis product. α-Amylase action results in various kind of oligomeric sugars (32) whilst glucoamylase hydrolysis product is glucose (monomeric). Based on this principle, successful separation of GLL1 from AMY was assessed through their hydrolysis products, as analyzed with Thin Layer Chromatography (TLC) (16), using purified GLL1 that was obtained from the purification procedure carried out at pH 5.8.

Purified GLL1 was incubated with soluble starch substrate at 37oC and samples were taken after 5, 10, 15 and 45 minutes of incubation. Each sample was then applied onto a silica gel 60 TLC plate (20 cm x 20 cm, Merck, Darmstadt, Germany) using capillary glass tube. The plate was then placed in a TLC separation chamber that had been equilibrated with the mobile phase, which was the mixture of butanol : ethanol : water (at a ratio of 5:5:3, v/v/v). The plate was retrieved from the TLC separation chamber after the mobile phase migration reached three quarters of the length of the plate and then immediately short-submersed (few seconds) in a mixture of water : methanol : sulphuric acid (at a ratio of 45:45:10, v/v/v). The plate was then heated at 120oC for 15 minutes on a hot plate for visualization of the sugars.

The TLC profile (Fig. 4) shows that the product of starch hydrolysis by GLL1 was solely monosaccharide i.e. glucose (G1). The intensity of the G1 spot was increasing from 5 to 45 minutes of incubation, indicating that glucose was produced over time. The spots at the sample application points were from the un-hydrolyzed soluble starch substrate. Apparently, most of the starch molecules were hydrolyzed within 15-45 minutes of incubation. Most importantly, no higher oligomeric sugar forms were detected, even during the first 5-10 minutes of incubation, suggesting the absence of α-amylase. This result provided the undisputed proof that GLL1 had successfully been separated from AMY despite of α-amylase activity was detected in the GLL1 fraction. Thus, the detected αamylase activity was not originated from AMY.

**Figure 4.** Analysis of starch hydrolysis products of GLL1 on a TLC plate (16). Lane 1-4 is hydrolysis products after 5, 10, 15, and 45 minutes, respectively. Lane 5 is the mono- (G1), di-(G2), and oligosaccharide (G3-G7) markers.

## **4.5. Properties of AMY and GLL1**

278 Chromatography – The Most Versatile Method of Chemical Analysis

and GLL1 indicated successful separation of the two enzymes.

(5), where α-amylase activity was detected in GLL1 fractions.

obtained from the purification procedure carried out at pH 5.8.

**4.4. Analysis of starch hydrolysis products to discriminate AMY and GLL1**

Glucoamylase hydrolyzes starch at the non-reducing end of amylose or amylopectin to result in glucose therefore the enzyme activity assay is normally based on the release of reducing sugar (21, 29). Unfortunately, α-amylase random digestion of starch may also result in reducing sugar i.e. maltose (30) therefore detection of glucoamylase activity in an α-amylase preparation can be anticipated. Glucoamylase can also, at lesser extent, act on 1,6 α-glycosidic bond of amylopectin (31), although its efficiency diverse greatly depending on the source organisms. Hydrolysis of the 1,6-α- bond may result in a less integrated starch molecule, hampering the formation of iodine-starch complex, which is the basis of standard α-amylase activity assay (20). Therefore, cross-detection of the two enzyme activity assays is inevitable. This phenomenon was notorious upon separation of AMY from GLL1 at pH 5.8

One approach to discriminate α-amylase and glucoamylase is *via* the evaluation of their hydrolysis product. α-Amylase action results in various kind of oligomeric sugars (32) whilst glucoamylase hydrolysis product is glucose (monomeric). Based on this principle, successful separation of GLL1 from AMY was assessed through their hydrolysis products, as analyzed with Thin Layer Chromatography (TLC) (16), using purified GLL1 that was

Purified GLL1 was incubated with soluble starch substrate at 37oC and samples were taken after 5, 10, 15 and 45 minutes of incubation. Each sample was then applied onto a silica gel 60 TLC plate (20 cm x 20 cm, Merck, Darmstadt, Germany) using capillary glass tube. The plate was then placed in a TLC separation chamber that had been equilibrated with the

peak distribution were different. Intensity of the peaks suggested much less materials were recovered from GLL1 digested samples than from AMY. This situation is likely contributed from highly hydrophobic or negative charged fragments that were not eluted from the RP-HPLC column due to their poor solubility in the solvent used (28). The amino acid sequences of AMY and GLL1 (see 3.5) showed that the charged amino acid distribution in the amino acid sequence of GLL1 (100 charged amino acid residues out of 494, ~20.2%) is higher than that of AMY (86 out of 468, ~18.4%), thus AMY contains more hydrophobic and non-charged amino acids. Unfortunately, correlation between the result from RP-HPLC and hydrophobic amino acid distribution cannot firmly be established because the nature of chymotrypsin digestion was unclear. The results may indeed indicate that hydrophobic amino acids in GLL1 are likely more clustered than in AMY, resulting in highly hydrophobic peptide fragments in GLL1. This hypothesis may be related to the structure of GLL1 (Fig. 8), which consists of one globular molecule that opposes two separable domains of AMY. However, possibility for the existence of highly negative charged peptide fragments from GLL1 can also not be excluded. Nevertheless, of the eluted fragments (Fig. 3), majority of AMY fragments are localized at peaks 5 and 6 whilst GLL1 are at peaks 3 and 4, suggesting different fragmentation had occurred. Thus, the RP-HPLC profiles of AMY

> Since the purity of AMY and GLL1 was firmly established, each enzyme could now be characterized. AMY was found active at a pH range of 5.0-7.5 and a temperature range of 30-60oC, with an optimum at 5.5 and 50oC, respectively (5). GLL1 was also found active a pH range of 4.6-6.8 and a temperature range of 30-65oC, with an optimum at 5.6-6.4 and 50oC, respectively (16). These results suggest that the two enzymes are active at similar conditions. Moreover, their characters are similar to most other *S. fibuligera* amylases.

> Recently, the amino acid sequences of AMY (GenBank accession code HQ172905) and GLL1 (HQ415729) were successfully elucidated as derived from their chromosomal DNA (16, 33).

The amino acid sequence of AMY is similar to that of published earlier by Itoh and coworkers (34), having mutations at six amino acid residues (Asp>Asn153, Ile>Val159, Ser>Asn190, Ser>Xaa216, Asp>Asn239, and Ser>Thr295). The six deviating residues in AMY comprise an additional predicted glycosylation site (Asn153), which according to a structural modelling, is highly plausible because it resides in a long surface loop (33).

Chromatography as the Major Tool in the Identification and

the Structure-Function Relationship Study of Amylolytic Enzymes from Saccharomycopsis Fibuligera R64 281

The hydrophobicity profile of AMY and GLL1 was analyzed based on the amino acid sequence of the mature enzymes (38), using an online analysis service at ExPASy (http://web.expasy.org/protscale/) (35). Another handful online service is also available at http://www.vivo.colostate.edu/molkit/hydropathy/index.html. The window size values for the frame normalization were scanned from 3 to 21 and compared to determine the significance of the regions that represent the hydrophobic character. The profile at window frame 9 (recommended for hydrophilic protein) and 19 (for hydrophobic protein) are presented in Fig. 5. The prediction points were linearly weighted with 100% relativity at the window edges. The score for hydrophobicity is ranged from -4.5 (hydrophilic) to 4.5 (hydrophobic), where the curve above the midpoint (zero) is interpreted as regions with hydrophobic character. As shown in Fig. 5, regions with hydrophobic character in AMY are at the residues 50-75, 210-230, 310-330, 375-410, and 425-440, whilst for GLL1 are at 225-230, 375-380, 425-445, and 475-480. This *ab initio* results suggest that AMY has more fragments with hydrophobic character, which may also be due to the presence of two separable domains in AMY that each has its own hydrophobic core. However, this result may not be corresponded with RP-HPLC profile as the fragmentation of AMY and GLL1 upon digestion with chymotrypsin occurs randomly, where

more hydrophobic fragments could be produced from GLL1.

**Figure 5.** Hydrophobicity profile of AMY and GLL1 as calculated according to their amino acid sequences (38). The upper graphs are produced with window frame of 9 whilst the lower was with

that of 19.

The amino acid sequence of GLL1 is similar to glucoamylases from strain HUT 7212 (GLU) and KZ (GLA) (16). GLU and GLA share high sequence identity (4), with only seven amino acid residues different. However, these seven residues are responsible for differences in their characteristics, where GLA exhibits potential thermal stability and GLU has higher affinity towards substrate. The amino acid sequence of GLL1 also differs to both GLU and GLA at precisely those particular seven residues, which four are being identical to GLU and three to GLA. Interestingly, these mutations result in GLL1 to adopt the potential thermal stability and higher affinity towards substrate. Thus, GLL1 behaves as a hybrid of GLU and GLA.

The calculated theoretical pI values (35) of AMY and GLL1 (based on the amino acid sequences of their mature protein sequences: AMY sequence starts at residue Glu27 of the full-length protein containing signal and pre-peptides as encoded by *AMY* gene whilst GLL1 sequence starts at residue Ala1 as reported in the data base) were 4.4 and 4.3, respectively. The isoelectric point (pI) value of AMY was confirmed experimentally, being 4.6±0.2. This finding supports the early observation that AEX or CEX columns are unable to separate the two enzymes. Similarly, the calculated theoretical molecular masses of AMY and GLL1 were also similar, being 51.7 kDa and 54.6 kDa, respectively. Both values were confirmed experimentally, being 54±2 kDa for AMY and 56±2 kDa for GLL1. Additional ±2 kDa of AMY and GLL1 masses is likely contributed from glycosylation. This result also confirms the inability of SEC column to separate them.

The activity of AMY is diminished in the presence of ethylene diamine tetra acetate (EDTA), a chelating agent. Inactivation of α-amylase by chelating agent is well known, as the enzyme requires calcium ion for its activity and integrity (36). This inactivation by chelating agent was not observed with GLL1. However, the activities of both AMY and GLL1 were increased in the presence of calcium or magnesium ions (data not shown). Further, AMY demonstrated lower activity in phosphate-citrate buffer. This phenomenon may arise from the citrate, which is reported to interact with calcium ion (37). The activity of AMY decreased concomitantly with the increase of the citrate buffer concentration but was fully recovered upon back-titration with calcium chloride (data not shown). Based on this observation, citrate buffers (1, 5) should only be used with considerable reservation and the use of tris buffer is recommended. This recommendation is in line with the necessity to perform purification at basic pH for HIC. In addition, this finding may also explain higher glucoamylase activity detected in AMY fractions when the purification was carried out using phosphate citrate buffer, pH 5.8 (5). Since the buffer does not influence the glucoamylase activity as it does to AMY, the disparity between the two enzyme activities was much less pronounced in comparison to the one presented in Fig. 2.

confirms the inability of SEC column to separate them.

was much less pronounced in comparison to the one presented in Fig. 2.

GLA.

The amino acid sequence of AMY is similar to that of published earlier by Itoh and coworkers (34), having mutations at six amino acid residues (Asp>Asn153, Ile>Val159, Ser>Asn190, Ser>Xaa216, Asp>Asn239, and Ser>Thr295). The six deviating residues in AMY comprise an additional predicted glycosylation site (Asn153), which according to a

The amino acid sequence of GLL1 is similar to glucoamylases from strain HUT 7212 (GLU) and KZ (GLA) (16). GLU and GLA share high sequence identity (4), with only seven amino acid residues different. However, these seven residues are responsible for differences in their characteristics, where GLA exhibits potential thermal stability and GLU has higher affinity towards substrate. The amino acid sequence of GLL1 also differs to both GLU and GLA at precisely those particular seven residues, which four are being identical to GLU and three to GLA. Interestingly, these mutations result in GLL1 to adopt the potential thermal stability and higher affinity towards substrate. Thus, GLL1 behaves as a hybrid of GLU and

The calculated theoretical pI values (35) of AMY and GLL1 (based on the amino acid sequences of their mature protein sequences: AMY sequence starts at residue Glu27 of the full-length protein containing signal and pre-peptides as encoded by *AMY* gene whilst GLL1 sequence starts at residue Ala1 as reported in the data base) were 4.4 and 4.3, respectively. The isoelectric point (pI) value of AMY was confirmed experimentally, being 4.6±0.2. This finding supports the early observation that AEX or CEX columns are unable to separate the two enzymes. Similarly, the calculated theoretical molecular masses of AMY and GLL1 were also similar, being 51.7 kDa and 54.6 kDa, respectively. Both values were confirmed experimentally, being 54±2 kDa for AMY and 56±2 kDa for GLL1. Additional ±2 kDa of AMY and GLL1 masses is likely contributed from glycosylation. This result also

The activity of AMY is diminished in the presence of ethylene diamine tetra acetate (EDTA), a chelating agent. Inactivation of α-amylase by chelating agent is well known, as the enzyme requires calcium ion for its activity and integrity (36). This inactivation by chelating agent was not observed with GLL1. However, the activities of both AMY and GLL1 were increased in the presence of calcium or magnesium ions (data not shown). Further, AMY demonstrated lower activity in phosphate-citrate buffer. This phenomenon may arise from the citrate, which is reported to interact with calcium ion (37). The activity of AMY decreased concomitantly with the increase of the citrate buffer concentration but was fully recovered upon back-titration with calcium chloride (data not shown). Based on this observation, citrate buffers (1, 5) should only be used with considerable reservation and the use of tris buffer is recommended. This recommendation is in line with the necessity to perform purification at basic pH for HIC. In addition, this finding may also explain higher glucoamylase activity detected in AMY fractions when the purification was carried out using phosphate citrate buffer, pH 5.8 (5). Since the buffer does not influence the glucoamylase activity as it does to AMY, the disparity between the two enzyme activities

structural modelling, is highly plausible because it resides in a long surface loop (33).

The hydrophobicity profile of AMY and GLL1 was analyzed based on the amino acid sequence of the mature enzymes (38), using an online analysis service at ExPASy (http://web.expasy.org/protscale/) (35). Another handful online service is also available at http://www.vivo.colostate.edu/molkit/hydropathy/index.html. The window size values for the frame normalization were scanned from 3 to 21 and compared to determine the significance of the regions that represent the hydrophobic character. The profile at window frame 9 (recommended for hydrophilic protein) and 19 (for hydrophobic protein) are presented in Fig. 5. The prediction points were linearly weighted with 100% relativity at the window edges. The score for hydrophobicity is ranged from -4.5 (hydrophilic) to 4.5 (hydrophobic), where the curve above the midpoint (zero) is interpreted as regions with hydrophobic character. As shown in Fig. 5, regions with hydrophobic character in AMY are at the residues 50-75, 210-230, 310-330, 375-410, and 425-440, whilst for GLL1 are at 225-230, 375-380, 425-445, and 475-480. This *ab initio* results suggest that AMY has more fragments with hydrophobic character, which may also be due to the presence of two separable domains in AMY that each has its own hydrophobic core. However, this result may not be corresponded with RP-HPLC profile as the fragmentation of AMY and GLL1 upon digestion with chymotrypsin occurs randomly, where more hydrophobic fragments could be produced from GLL1.

**Figure 5.** Hydrophobicity profile of AMY and GLL1 as calculated according to their amino acid sequences (38). The upper graphs are produced with window frame of 9 whilst the lower was with that of 19.

Another computational study was performed using the on-line hydrophobic cluster analysis (HCA) program (39). This approach has previously been done to compare the hydrophobic clusters in α-amylases (40). The sequence of AMY and GLL1 were submitted to the drawcha server (http://bioserv.impmc.jussieu.fr/hca-form.html) and the resulted profiles were analyzed manually following Gaboriaud *et al.* (41). The profiles were manually mismatched with the structures of AMY (*in silico*, generated from its homolog, see section 5) and GLU, correspondingly. AMY was predicted to have 34 hydrophobic clusters whilst GLL1 was 17. However, only 3 hydrophobic clusters of AMY are exposed on the protein surface as oppose to 7 of GLL1. The analysis suggests the presence of more hydrophobic patches on the surface of GLL1 than on that of AMY. Thus the result corresponded with the experimental finding with HIC column.

Chromatography as the Major Tool in the Identification and

the Structure-Function Relationship Study of Amylolytic Enzymes from Saccharomycopsis Fibuligera R64 283

X-ray crystallography. Nevertheless, based on the result of overexpression, the attempt for heterelogous expression of AMY in *E. coli* was withdrawn for the use of expression system

*S. cerevisiae* has been a popular choice for heterelogous expression of proteins (48) and tested to express α-amylase and glucoamylase from *S. fibuligera* strain HUT 7212 (49, 50). *S. cerevisiae* expression system for expression of *S. fibuligera* α-amylase and glucoamylase has also been improved, i.e. the enzymes expression was drastically increased under the control of the *S. cerevisiae* constitutive phosphoglycerate kinase (PGK) promoter (51). Motivated from that success, AMY was overexpressed in *S. cerevisiae* INVSc1, using yEP-secretex vector and galactose as the inducer. However, although AMY recombinant demonstrated similar activity to AMY, its molecular mass is substantially higher, being 67±2 kDa. Unfortunately, no such information was provided from the expression of *S. fibuligera* α-amylase from strains HUT 7212 (50) and Eksteen (51). Since no polypeptide chain or protein tag was added, different glycosylation profile was likely the reason for the increase in the molecular mass of the recombinant protein. The finding of a plausible additional glycosylation site in

from higher organism i.e. yeast *S. cerevisiae*.

the amino acid sequence of AMY applauds this proposal.

respective buffer in the cold room (~4oC).

**5.1. AMY recombinant behavior upon purification in AEX columns**

Likewise AMY, AMY recombinant was secreted into overproduction medium therefore it was harvested by cold centrifugation at 6000 *g* for 30 minutes to remove the yeast cells. The enzyme was captured from the medium by fractionation with ammonium sulfate at a saturation degree of 0-100%. The ammonium sulfate protein precipitate was dissolved in 50 mM Tris-HCl buffer, pH 8.0 (~5 ml) and then dialyzed overnight against one litre of that

The dialyzed AMY recombinant was then loaded onto a resource-Q AEX column (GE Healthcare Europe GmbH, Diegem, Belgium), which had been equilibrated with that respective buffer. Purification was performed in a cold cabinet (~7oC) using a fast protein liquid chromatography (FPLC) Äkta system (GE Healthcare Europe GmbH, Diegem, Belgium) and monitored on-line with the Unicorn program. The enzyme was recovered from the column upon an elution with an increasing gradient of sodium chloride 0-1 M. As

Purification of AMY in the resource-Q column showed that minor contaminants were still present in the purified sample (Fig. 6). The major contaminant has, however, higher absorbance at λ 260 nm, suggesting that it may not be protein. This contaminant appeared yellowish in colour, which may be originated from the overproduction medium component that was co-purified in HIC, AEX, and SEC columns. Further, the elution profile of AMY recombinant was very similar to that of AMY, except for an additional large protein peak upon sample application and washing step prior to the sodium chloride salt gradient. These additional peaks unmistakably are originated from other proteins and components of overproduction medium because the AMY recombinant sample applied was not purified prior to this column. AMY recombinant was not pre-purified with the HIC because no

the control, purified AMY was also subjected to purification with the same column.