CFD Simulations AC211
DelftJetinHotCoflow (DJHC) burner
Application Challenge AC211 © copyright ERCOFTAC 2023
CFD Simulations
Introduction
The DJHC flames have been studied in the literature using several different approaches. For a recent overview we refer to Perpignan et al. (2018). Here we present results from a comparative study of three different turbulencechemistry approaches also discussed in Perpignan et al. (2018). These three approaches are: Conditional Sourceterm Estimation (CSE) model, a new variant of the Eddy Dissipation Concept (EDC) model and a new variant of the Flamelet Generated Manifolds (FGM) / presumed PDF model.
Before reviewing all details of the model setup, the three modelling approaches for turbulence chemistry interaction are described in some detail.
Turbulencechemistry Interaction Approaches
Eddy Dissipation Concept
General introduction
The Eddy Dissipation Concept (Magnussen, 1981) has been frequently employed in simulations related to the Flameless combustion regime due to its availability in several CFD codes, the fact that it contains no assumptions related to the combustion regime, and the possibility of directly employing a chemical reaction mechanism of choice.
In many simulations of JetinHot Coflow flames, the standard EDC model was shown to predict ignition far upstream, leading to a clear peak in the radial mean temperature profile not present in the experiments. This is due to the fact that the EDC by its definition in no way represents ignition by separate local events (ignition kernels) responsible for the ignition in many JHC flames. Nevertheless, several authors attempt to improve the model by modifying the model constants (Christo & Dally, 2005; Sarras et al., 2014; De & Dongre, 2015; De et al., 2011; Aminian et al., 2012; Frassoldati et al., 2010; Mardani, 2017). Indeed, in the EDC model, two model parameters appear. The key model features are the volume fraction of the reaction zones (γ) and the residence time in these zones (t), which are dependent on two constants, Ct and Cγ respectively. For example, Aminian et al., 2012 proposed adjusting the value of Ct, while Evans et al. (2015) investigated a range of values for the model constants and recommended values for CH_{4}/H_{2} and C_{2}H_{4}/H_{2} flames. Both works were based on the Adelaide burner experiments and showed improvement in relation to the standard constants.
The proposed changes of model constants in the EDC model were not only purely empirical but were also motivated by the argument that the flameless combustion regime is different from the conditions for which the EDC was formulated originally. The modification of EDC models can potentially take into account the increased volumes caused by the different shapes of the reaction zones in the FC regime. This was in corroboration with earlier findings, also pointed in the conclusions of De & Dongre, 2015.
EDC model with local parameters (EDCLP)
Using the EDC, Bao, (2017) presented a formulation with locally determined constants (EDCLP), which was compared to the standard EDC formulation. The laminar and turbulent flame speeds are calculated in order to evaluate the constants, leading an expression for C_{γ} proportional to Da*^{3/4}, where Da* is the Damk”hler number based on the Kolmogorov scale (i.e. the inverse of the more commonly used Karlovitz number). Simulations were performed with the standard kε model and with a εequation based Reynolds stress model (RSM) using ANSYS Fluent®. with linear pressurestrain model and Lien and Leschziner model for turbulent diffusion (RSM). Simulations were performed in ANSYS Fluent®.
According to the paper of Ertesvåg & Magnussen (2000) in which the energy cascade model is discussed, it holds that the supply of mechanical energy (w) and the energy dissipation rate at a given energy cascade level (q) of the nth eddy are given by Eqs. 1 and 2, where ω is the strain rate, u is the velocity, ν is the viscosity, and C_{D1} and C_{D2} are model constants.
Eq. 1  
Eq. 2 
By assuming quasisteadiness, from energy conservation:
Eq. 3 
Considering the smallest eddy (*):
Eq. 4 
Ertesvåg & Magnussen (2000) showed that Therefore:
Eq. 5  
Eq. 6 
Instead of imposing the model constant values, Bao, (2017) used the ratio of the two constant (obtained by combining Eqs. 5 and 6) to find expressions for the model constants that can be locally defined, a formulation here named EDCLP (EDC with local paramenters).
Eq. 7 
Considering the velocity to be the turbulent flame speed, which can in turn be approximated by a relation containing the laminar flame speed and the turbulent Reynolds number (), the following relation is obtained:
Eq. 8 
Having the chemical timescale and the Damköhler number of the Kolmogorov scale :
Eq. 9 
Using the turbulent flame speed again, and the definition of the Damköhler number of the Kolmogorov scale , where is the Kolmogorov timescale, Eq. 5 can be rewritten as:
Eq. 10 
From this an expression for the constant can be written:
Eq. 11 
Using Eqs. 9 and 11, the original EDC constant values can be locally calculated:
Eq. 12  
Eq. 13 
Conditional Sourceterm Estimation
General introduction
The Conditional Sourceterm Estimation (CSE) model (Labahn et al, 2015;Labahn et al, 2016) is not as widely employed as the EDC, possibly due to its complexity, the fact that it is relatively recent, and its absence in commercial codes.
The CSE shares similarities with the Conditional Moment Closure (CMC) model in its formulation. Both CSE and CMC use conditional averages to compute the chemical source term, while the conditional fluctuations are assumed to be negligible (conditioned to mixture fraction). In the CMC, transport equations for the conditional means of species and temperature are solved, while in the CSE these conditional means are obtained by inverting the integral equation that relates the unconditional to the conditional mean in a relation involving the (presumed) PDF. The model's absence of assumptions related to the combustion regime and its representation of turbulent fluctuations make it attractive for simulations involving the Flameless regime.
CSE as applied to the DJHC burner
Labahn et al. (2015, 2016) applied the CSE to model the DJHC flames using RANS and LES. Their approach employed two mixture fractions as conditioning variables: the first was the conventional mixture fraction and the second was a variable representing at the same time the oxygen profile and temperature (or enthalpy deficit profile) at the inflow boundary. Handling oxygen and temperature profile with only one variable induces an approximation but was found to be effective. LES calculations were performed using the standard Smagorinsky SGS model. RANS calculations were performed with both the standard and the realizable kε models and were found to give only slightly different results. The results herein shown are for the standard kε model. Chemistry was tabulated with an approach called Trajectory Generated Low Dimension Manifold (TGLDM), explained in (Labahn, 2015), in which chemical states and reaction rates are stored as a function of the two mixture fractions and mass fractions of CO_{2} and H_{2}O. Simulations were performed using OpenFOAM (release not mentioned).
The CSE shares similarities with the CMC (Conditional Moment Closure) approach (Klimenko & Bilger, 1999), as both employ conditional averages to calculate the chemical source term. The underlying assumption is that the fluctuation of a given scalar can be associated to the fluctuations of a key quantity. Labahn et al. (2015) represent the conditional chemical source term of species k as the following function:
Eq. 14 
In this equation, η represents the value in sample space of the conditioning variable (e.g. mixture fraction). And for any stochastic variable A, 〈Aη〉 is the expected value conditional on the conditioning variable taking the value . The conditional average of the source term depends on the conditional averages of mass fraction (Y_{k}), temperature (T) and density (ρ). In contrast with CMC, these quantities do not have transport equations in the CSE. Instead, transport equations only for the Favreaveraged species mass fractions are included. Therefore, the mean unconditional reaction rates are determined with:
Eq. 15 
In which the spatial coordinate(x_{j}), time (t), and Favreaveraged PDF () are present. Labahn et al. (2015) assumed a βPDF distribution. In order to calculate , the integral that equals the Favreaveraged mass fraction (Eq. 16) is inverted to find .
Eq. 16 
The averages were conditioned to values of mixture fraction (Z) and of a variable able to capture the influence of the air stream in the DJHC problem. This variable (W) is defined by Eq. 17 and is equal to 0 at the location with the minimum amount of oxygen in the coflow, and is equal to 1 where the oxygen concentration is maximum (air).
Eq. 17 
The integral to be inverted was then:
Eq. 18 
Having the conditional mean after the inversion, is determined by consulting the chemical lookup table.
Flamelet Generated Manifolds
General introduction
First presented by van Oijen & de Goey (2000), the Flamelet Generated Manifolds (FGM) approach is based on tabulatedchemistry, in which laminar flame structures are considered and tabulated in a precomputation step, in which detailed chemistry is adopted. The important thermochemical properties and the source terms of a selected progress variable are tabulated as functions of independent variables (mixture fraction and/or progress variable, enthalpy), which are solved for in the CFD calculation. For turbulent flows usually an assumed PDF approached is used for the independent variables.
The relatively low computational costs make the FGM approach attractive and many works have been performed with the approach. Earlier works applying an FGM approach to the DJHC are (Sarras et al., 2014), (Abtahizadeh et al., 2015) and (Abtahizadeh et al., 2017). The approach presented here is a new variant presented in (Huang et al., 2017) and (Huang, 2018)
FGM with diluted air (DAFGM)
The application of FGM herein shown was developed to represent the dilution of reactants with burnt gases. Developed by Huang et al. (2017, 2018), the DAFGM (Diluted Air FGM) are generated having four control variables: a mixture fraction, a progress variable, enthalpy loss, and air dilution level. The progress variable was defined as the linear combination of the CO_{2}, H_{2}O, CO and H_{2} mass fractions. Therefore, the counterflow diffusion flamelets were calculated and tabulated based on this fourdimension set of variables. The RANS simulations were performed with the standard kε model and the LES used the standard Smagorinsky model. Simulations were performed using OpenFOAM, release 2.3.1, extended with a new storage and retrieval routine.
The precalculated flamelets were parametrized by four variables: a mixture fraction (ξ), the progress variable mentioned above (C), enthalpy loss (η) and dilution level (γ). The mixture fraction is defined by Eq. 19, where Z is the conventional mixture fraction based on fuel and air and ? represents the mass fraction of diluent (combustion products of stoichiometric combustion) in the mixture (Eq. 20).
Eq. 19  
Eq. 20 
In Eq. 20, Y_{d} is the dilution variable, defined by the sum of products mass fractions (CO_{2} and H_{2}O). The superscript Dil denotes the dilution variable value in combustion products of stoichiometric combustion, i.e. at the stoichiometric mixture fraction (Z_{st}). The dilution level of the air in the counterflow flames underlying the FGM then is given by:
Eq. 21 
The enthalpy loss is defined by Eqs. 22, in which is the enthalpy, is the adiabatic enthalpy defined by Eq. 23, while and are the enthalpies of the diluent for maximum and minimum enthalpy loss, respectively.
Eq. 22  
Eq. 23 
In the calculation of the DJHC flame radiative heat transfer can be neglected but the temperature variation at the coflow inlet leads to an enthalpy variation that should be taken into account. Accordingly, was taken zero and from the experimental data on the temperature, was taken as 2x10^{6} J/kg.
The four variables employed to parametrize the manifolds have values at every location in space from the transport equations of mixture fraction, progress variable, enthalpy and dilution variable that are solved throughout the domain. In the framework of RANS and LES transport equations are solved respectively for mean and for resolved variables. And for mixture fraction also the variance is calculated in order to determine the assumed βPDF. Therefore a fivedimensional FGM lookuptable is used. For further details we refer to (Huang, 2018).
Calculation Domain and Grids
The characteristics of the different modelling approaches for which calculations are compared are listed in Table 3. The selection of models and parameter values is explained following the rows in the table, starting with domain and grid size.
The computational domain extends from the fuel nozzle tip in axial direction beyond the last measurement station and in radial direction to a distance much wider than the width of the coflow. The latter is necessary because the entrainment of laboratory air into the coflow stream was shown to be crucial for accurate predictions and has to be included in the simulation.
RANSEDC  RANSCSE  RANSFGM  LESCSE  LESFGM  

Domain (mm)  2D  2D  2D  3D cylindrical  3D hexahedral  
Axialradial  225 x 80  225 x 80  225 x 80  225 x 80  225 x 80  
Grid size (#)  190 x 145
=27175 
31250  300 x 215
=64000 
1.5×10^{6}  2.5×10^{6}  
Platform  ANSYS  OpenFOAM  OpenFOAM  OpenFOAM  OpenFOAM  
Turbulence  RANS  RANS  RANS  LES  LES  
Closure model  kε or RSM  Standard kε  Standard kε  Smagorinsky  Smagorinsky  
Kinetic scheme  GRI 1.2  GRI 2.11  GRI 3.0  GRI 2.11  GRI 3.0  
Reduction  DRM19  TGLDM  DAFGM:
Le=1 flamelets 
TGDLM  DAFGM:
Le=1 flamelets  
Turbcheminteraction  EDCLP
variable C_{T} 
Multistream CSE (two mixture fractions, βpdf)  Assumed PDF (, βpdf for mixture fraction and for scaled progress variable)  Multistream CSE (two mixture fractions, βpdf))  Assumed PDF (, βpdf for mixture fraction and for scaled progress variable)  
Thermochemical Scalar equations  Mean of species
Mean of enthalpy 
Mean and variance of two mixture fractions  Resolved mixture fraction and resolved progress variable
Resolved enthalpy 
Resolved mixture fractions and SGS variances
Resolved enthalpy 
Resolved mixture fraction and resolved progress variable
Resolved enthalpy  
Radiation  Not included  Not included  Not included  Optically thin.
TRI not included. 
Not included  






Figure 9: Twodimensional grid employed by Bao (2017). Axial direction from left to right. Radial direction from bottom to top. 
All simulations considered three inlets: fuel inlet through the jet nozzle, coflow inlet, and air inlet. An example of a twodimensional computation domain and grid used in RANS is given in Figure 9. The grid is refined close to the nozzle, near the centre and in the regions where a shear layer is expected (as indicated by the dashed rectangles). A grid independence study for RANSEDC simulations has been reported in (De et al., 2011) and concluded that a grid of 180x125=22500 cells is sufficiently fine. Details of the grids used in this study are given in Table 3.
Software and solvers
EDC
Steady RANS simulations are performed using ANSYSFluent, version 14.5. The EDCLP model is implemented using userdefined functions. (Bao, 2017)
CSE
OpenFOAM is used for both RANS and LES. (Labahn, 2015, 2016)
DAFGM
OpenFOAM (version 2.3.1) is used for both RANS and LES. The momentum equations are solved using a pressurebased solver. The transport equations for the scalar variables are solved In every iteration after the velocity prediction, and the controlling parameters needed for table lookup are derived and properties are retrieved from the table. (Huang, 2017, 2018)
Turbulence model
The turbulence chemistry interaction approaches have been used with different models for the statistical description of turbulence. The three considered approaches have been used with RANS (ReynoldsAveraged NavierStokes) using the standard kε model (Bao, 2017) (Labahn, 2015, 2016) (Huang, 2017, 2018).
The simulations using EDC were also done using a standard Reynolds stress model (RSM) more precisely an εequation based RSM with linear pressurestrain model and Lien and Leschziner model for turbulent diffusion (ANSYSFluent, 2013) (Bao, 2017).
CSE and DAFGM have also been combined with LES (Labahn, 2015, 2016) (Huang, 2018). The standard Smagorinsky model was used for closure of the subgrid stress.
Boundary Conditions
The boundary conditions at the inflow of fuel and coflow are set using information available for measurements at z=3mm, but differences between the different approaches are mentioned below.
Neumann boundary conditions are imposed at the outflow. At the side boundary on the air side a slip wall boundary condition was assumed. In the 2D RANS simulations, symmetry is imposed at the axis of symmetry. LES simulations were performed in 3D and then the symmetry axis is not a boundary.
EDC
The experimental data measured at the axial position of 3 mm are employed to determine the jet and coflow inlet boundary conditions. Equilibrium is assumed in order to estimate the species distribution at the boundary based on the available temperature and O_{2} profiles. Bao (2017) imposed the turbulence kinetic energy based on the available profiles of the Reynolds stresses, with the relation shown in Eq. 24, while the dissipation rate followed the relation in Eq. 25. In these equations, are the axial, radial, azimuthal, and shear Reynolds stresses, respectively. The mean axial velocity is represented by .
Eq. 24  
Eq. 25 
The outer boundary is defined as a slip wall, having zero shear stress. A constant profile for mean axial velocity equal to 0.3 m/s is imposed at the air inlet.
CSE
The experimental data measured at the axial position of 3 mm was employed to determine the jet and coflow inlet boundary conditions. Equilibrium was assumed in order to determine the thermochemical variables at the boundary based on the information on mean temperature and O_{2} profiles. The oxygen mass fraction profile is assumed to have an analytical relation with the mean temperature profile. See below the thermochemical scalar equations for more details. The resulting oxygen mass fraction profile is shown in Figure 2.
The fuel composition was assumed to be 85% CH4 and 15% N2.
For the LES simulations, turbulence at the inlet is generated using the synthetic turbulence method of Kornev & Hassel (2007) based on prescribed mean velocity profiles, the corresponding Reynolds stress tensor and integral length scales. Turbulent fluctuations for temperature, mixture fraction, oxidiser split or species mass fractions are assumed zero at the inlet.
The outer boundary is defined as a slip wall. A constant profile for mean axial velocity equal to 0.5 m/s is imposed at the air inlet. The choice of this value is further discussed in the Best Practice Advice section.
DAFGM
For the RANS simulations the inlet boundary conditions used the available experimental data at 3 mm, except for the oxygen concentration. The turbulent kinetic energy and the dissipation rate were obtained in the same manner as in the EDC simulations explained before.
Equilibrium was assumed in order to estimate the scalar variables at the boundary based on the available temperature and a modified O_{2} profiles. In order to exactly satisfy the constraint that the fuel/air ratio agrees with the fuel and air flow rates to the coflow burner measured by mass flow controllers, the mixture fraction boundary condition was modified and, as can be seen in Figure 10, the resulting assumed O_{2} – profile has a lower value in the inner region r < 25mm than measured with the sampling probe.
For the LES simulations, turbulence at the inlet is generated using the synthetic turbulence method of Kornev & Hassel (2007) based on prescribed mean velocity profiles, the corresponding Reynolds stress tensor and integral length scales. No turbulent fluctuations for enthalpy and dilution variable were included at the inlet.
Figure 10: O_{2} mass % in experiment (Oldenhof et al, 2011) and in simulations (CSE and DAFGM) 
Kinetic scheme and chemistry reduction
All simulations are based on versions of the GRImech detailed mechanism, but they significantly differ in the chemistry reduction method.
EDC
In EDC transport equations are solved for mean values of all species. The use of the DRM19 mechanism, which is derived from GRImech 1.2 allows significant reduction in computational cost.
CSE
Employing GRImech 2.11, a lower dimensional manifold for chemical evolution is derived using the trajectory generated low dimensional manifold (TGLDM) method. GRImech 2.11 is a mechanism for 49 species, but a lower dimensional manifold in composition space is derived from a set of trajectories in this space. Mixing variables and a limited number of reacting coordinates determine the dimension of the manifold (Labahn, 2015). See below under scalar variables for more information.
DAFGM
Employing GRImech 3.1, a lower dimensional manifold was derived using a new variant of the Flamelet Generated Manifold (FGM) method. GRImech is a mechanism for 53 species, but a lower dimensional manifold for species evolotion is derived from steady laminar counterflow flames between fuel and diluted and cooled air for strain rates from low value to extinction strain rate and also including an extinguishing flamelet. The flamelets are calculated assuming unity Lewis number for all species (Huang, 2017, 2018) The method has been explained in more detail above in the section on turbulencechemistry interaction.
Thermochemical scalar equations
EDC
In EDC transport equations are solved for mean values of species mass fractions and enthalpy. Since the DRM19 has 19 species in total 20 equations are solved.
CSE
To describe the local conditions two nonreacting scalars and two reacting scalars (H_{2}O and CO_{2} mass fractions) are used.
The first nonreacting scalar is mixture fraction Z based on fuel and air. The second nonreacting scalar W represents the enthalpy variation at the coflow inlet. It is defined as a normalised temperature: for r > 2.25mm for r < 2.25mm. The value of temperature for radial distance larger than the last experimental point, i.e. from 35.1 mm to 44.6 mm, was obtained by linear extrapolation of the known profile. The mass fraction of oxygen at the inlet boundary is assumed to vary linearly with W between minimal value in the middle of the coflow temperature maximum (W=0) and maximal value in the air (W=0). Mean (or resolved) values of these variables are solved (4 equations) and equations for the variance of the nonreacting scalars. In total 6 equations are solved. In addition to that, unconditional and conditional averages of all species mass fractions are calculated through the CSE algorithm (Labahn, 2015). In the LES, the mean enthalpy equation was also solved and the radiation heat loss was taken into account via an optically thin model (emission only).
DAFGM
To describe the local conditions, four scalar quantities are needed: mixture fraction, progress variable, dilution variable and enthalpy. In turbulent flow equations the mean (or resolved) values of these quantities are solved (4 equations) and equations for the variance of mixture fraction. In total 5 equations are solved. Fluctuations in progress variable are not taken into account, but it was checked that including them has negligible influence on the results. (Huang, 2018)
References
 Abtahizadeh E, de Goey P, van Oijen J. Development of a novel flameletbased model to include preferential diffusion effects in autoignition of CH4/H2 flames. Combustion and Flame 2015; 162:43584369.
 Abtahizadeh E, de Goey P, van Oijen J. LES of Delft JetinHot Coflow burner to investigate the effect of preferential diffusion on autoignition of CH4/H2 flames. Fuel 2017;191:3645.
 Aminian J, Galletti C, Shahhosseini S, Tognotti L. Numerical investigation of a MILD combustion burner: analysis of mixing field, chemical kinetics and turbulencechemistry interaction. Flow, Turbulence and Combustion 2012; 88:597623.
 Ansys Fluent, Ansys fluent theory guide 15.0, Inc, Canonsburg, PA (2013).
 Bao H. Development and validation of a new Eddy Dissipation Concept (EDC) model for MILD combustion, MSc Thesis, Delft University of Technology, 2017.
 Christo FC, Dally BB. Modeling turbulent reacting jets issuing into a hot and diluted coflow. Combustion and Flame 2005;142:117129.
 De A, Dongre A. Assessment of turbulencechemistry interaction models in MILD combustion regime. Flow, Turbulence and Combustion 2015;94:439478.
 De A, Oldenhof E, Sathiah P, Roekaerts DJEM. Numerical simulation of DelftJetinHotCoflow (DJHC) flames using the Eddy Dissipation Concept model for turbulencechemistry interaction. Flow, Turbulence and Combustion 2011;87:537567.
 Ertesvåg IS, Magnussen BF. The eddy dissipation turbulence energy cascade model. Combustion science and technology 2000;159:213235.
 Evans MJ, Medwell PR, Tian ZF. Modeling lifted jet flames in a heated coflow using an optimized eddy dissipation concept model. Combustion Science and Technology 2015;187:10931109.
 Frassoldati A, Sharma P, Cuoci A, Faravelli T, Ranzi E. Kinetic and fluid dynamics modeling of methane/hydrogen jet flames in diluted coflow. Applied Thermal Engineering 2010;30:376383.
 Huang X, Tummers MJ, Roekaerts DJEM. Experimental and numerical study of MILD combustion in a labscale furnace. Energy Procedia 2017;120:395402.
 Huang X. PhD Thesis, Delft University of Technology, 2018
 Klimenko AY, Bilger RW. Conditional moment closure for turbulent combustion. Progress in Energy and Combustion Science 1999;25:595687.
 Kornev N, Hassel E. Method of random spots for generation of synthetic inhomogeneous turbulent fields with prescribed autocorrelation functions. Communications in Numerical Methods in Engineering 2007;23:3543.
 Labahn JW, Devaud CB. Large Eddy Simulations (LES) including Conditional Sourceterm Estimation (CSE) applied to two DelftJetinHotCoflow (DJHC) flames. Combustion and Flame 2016;164:6884.
 Labahn JW, Dovizio D, Devaud CB. Numerical simulation of the DelftJetinHotCoflow (DJHC) flame using conditional sourceterm estimation. Proceedings of the Combustion Institute 2015;35:35473555.
 Mardani A. Optimization of the Eddy Dissipation Concept (EDC) model for turbulence chemistry interactions under hot diluted combustion of CH4/H2. Fuel 2017;191:114129.
 Perpignan A.A.V., Gangoli Rao A., Roekaerts, D.J.E.M., Flameless Combustion and its Potential Towards Gas Turbines, Progress in Energy and Combustion Science 2018; 69: 2862
 Sarras G, Mahmoudi Y, Arteaga Mendez LD, van Veen EH, Tummers MJ, Roekaerts, DJEM. Modeling of turbulent natural gas and biogas flames of the Delft JetinHotCoflow burner: effects of coflow temperature, fuel temperature and fuel composition on the flame liftoff height. Flow, Turbulence and Combustion 2014;93:607635.
 van Oijen JA, de Goey PH. Modelling of premixed laminar flames using flameletgenerated manifolds. Combustion Science and Technology 2000;161:113137.
Contributed by: André Perpignan, Dirk Roekaerts, E. Oldenhof, E.H. van Veen, M.J. Tummers, Hesheng Bao, Xu Huang — TU Delft
Contributed by: Jeffrey W. Labahn — Stanford University
© copyright ERCOFTAC 2018