CFD Online Logo CFD Online URL
Home > Wiki > Combustion


From CFD-Wiki

(Difference between revisions)
Jump to: navigation, search
(Flamelets in turbulent combustion)
(Flamelets in turbulent combustion)
Line 366: Line 366:
\widetilde{Y}_k= \int Y_F(\widetilde{Z},\widetilde{\chi}) P(Z) P(\chi) dZ d\chi
\widetilde{Y}_k= \int Y_F(\widetilde{Z},\widetilde{\chi}) P(Z) P(\chi) dZ d\chi
In [[Large eddy simulation (LES)]] context (see [[#LES equations]] for reacting  flow),
the [[probability density function]] is replaced by a [[subgrid PDF]] <math> \widetilde{P}</math>.
The same equation hold by replacing averaged values with filtered values.
\widetilde{Y}_k= \int Y_F(\widetilde{Z},\widetilde{\chi}) \widetilde{P}(Z) \widetilde{P}(\chi) dZ d\chi
The assumptions made regarding the shapes of the PDFs are still justified. In LES combustion the
[[subgrid variance]] is smaller than RANS counterpart (part of the large-scale fluctuations are solved)
and therefore the modelled PDFs are thinner.
====== Unsteady flamelets ======
==== Conditional Moment Closure (CMC)====
==== Conditional Moment Closure (CMC)====

Revision as of 11:17, 15 November 2005


What is combustion -- Physics versus modelling

Combustion phenomena consists of many physical and chemical processes with broad range of time scales. Mathematical description of combustion is not always trivial. Analytical solutions exists only for basic situations of laminar flame and because of its assumptions it is often restricted to few problems solved usually in zero or one-dimensional space.

Problems solved today concern mainly turbulent flows, gas as well as liquid fuels, pollution issues (products of combustion as well as for example noise pollution). These problems require not only extensive experimental work, but also numerical modelling. All combustion models must be validated against the experiments as each one has its own drawbacks and limits. However here the modelling part will be mainly addressed.

Reaction mechanisms

The combustion is mainly chemical process and although we can, to some extend, describe flame without any chemistry informations, for modelling of flame propagation we need to know the speed of reactions, product concentrations, temperature and other parameters. Therefore more or less detailed information about reaction kinetics is essential for any combustion model. Mixture will generally combust, if the reaction of fuel and oxidiser is fast enough to maintain until all of the mixture is burned into products. If the reaction is too slow, the flame will extinguish, if too fast, explosion or even detonation will occur. The reaction rate of typical combustion reaction is influenced mainly by concentration of reactants, temperature and pressure.

A stoichiometric equation of an arbitrary equation can be written as:

\sum_{j=1}^{n}\nu' (M_j) = \sum_{j=1}^{n}\nu'' (M_j),

where $\nu$ is the stoichiometric coefficient, M_j is arbitrary species. One prime specifies the reactants and double prime products of the reaction.

Reaction rate, expressing the rate of disappearance of reactant i of such a reaction, is defined as:

RR_i = k \, \prod_{j=1}^{n}(M_j)^{\nu'},

in which k is the specific reaction rate constant. Arrhenius found that this constant is a function only of temperature and this function is defined as:

k= A T^{\beta} \, exp \left( \frac{-E}{RT}\right)

where A is pre--exponential factor, E is activation energy and \beta is temperature exponent. These constants for given reactions can be found in literature. The reaction mechanism can be given from experiments for every reaction resolved, it could be also constructed numerically by automatic generation method (see [Griffiths (1994)] for review on reaction mechanisms). For simple hydrocarbon tens to hundreds of reactions are involved. By analysis and systematic reduction of reaction mechanisms global reaction (from one to five step reactions) can be found (see [Westbrook (1984)]).

Governing Equations for Reacting Flows

Together with the usual Navier-Stokes for compresible flows (See Governing Equations), additional equations are needed in reacting flows. The mass fraction transport equation for k-th species  Y_k is:

\frac{\partial}{\partial t} \left( \rho Y_k \right) +
\frac{\partial}{\partial x_j} \left( \rho u_j Y_k\right) = 
\frac{\partial}{\partial x_j} \left( \rho D_k \frac{\partial Y_k}{\partial x_j}\right)+ \dot \omega_k

where Ficks law is assumed for scalar diffusion with  D_k , the species difussion coefficient and  \dot \omega_k is the species reaction rate. A non-reactive scalar (like the mixture fraction  Z ) had the following transport equation:

\frac{\partial}{\partial t} \left( \rho Z \right) +
\frac{\partial}{\partial x_j} \left( \rho u_j Z \right) = 
\frac{\partial}{\partial x_j} \left( \rho D \frac{\partial Z}{\partial x_j}\right)

where  D is the diffusion coefficient of the passive scalar.

RANS equations

In turbulent flows, Favre averaging is often used and the mass fraction transport equation is transformed to

\frac{\partial \overline{\rho} \widetilde{Y}_k }{\partial t} +
\frac{\partial \overline{\rho} \widetilde{u}_j \widetilde{Y}_k}{\partial x_j}=
\frac{\partial} {\partial x_j} \left( \overline{\rho D_k  \frac{\partial Y_k} {\partial x_j} } -
 \overline{\rho} \widetilde{u''_i Y''_k } \right)
+ \overline{\dot \omega_k}

where the turbulent fluxes  \widetilde{u''_i Y''_k} and reaction terms   \overline{\dot \omega_k} needs to be closed.

The passive scalar turbulent transport equation is

\frac{\partial \overline{\rho} \widetilde{Z} }{\partial t} +
\frac{\partial \overline{\rho} \widetilde{u}_j \widetilde{Z} }{\partial x_j}=
\frac{\partial} {\partial x_j} \left( \overline{\rho D  \frac{\partial Z} {\partial x_j} } -
 \overline{\rho} \widetilde{u''_i Z'' } \right)

where similar to the mass fraction equation,  \widetilde{u''_i Z''} needs modelling.

In addition to the mean passive scalar equation, an equation for the Favre variance  \widetilde{Z''^2} is often employed

\frac{\partial \overline{\rho} \widetilde{Z''^2} }{\partial t} +
\frac{\partial \overline{\rho} \widetilde{u}_j \widetilde{Z''^2} }{\partial x_j}=
\frac{\partial}{\partial x_j}
\left(  \overline{\rho} \widetilde{u''_i Z''^2}  \right) -
 2 \overline{\rho} \widetilde{u''_i Z'' }
- \overline{\rho} \widetilde{\chi}

where  \widetilde{\chi} is the mean scalar dissipation rate defined as  \widetilde{\chi} =
2 D \widetilde{\left| \frac{\partial Z''}{\partial x_j} \right|^2 } This term and the variance diffusion fluxes needs to be modelled.

LES equations

The Large eddy simulation (LES) equation for reactive flows introduces equations for the filtered species mass fractions to the compressible flow field. Similar to #RANS equations, but using Favre filtering instead of Favre averaging. The filtered mass fraction transport equation is

\frac{\partial \overline{\rho} \widetilde{Y}_k }{\partial t} +
\frac{\partial \overline{\rho} \widetilde{u}_j \widetilde{Y}_k}{\partial x_j}=
\frac{\partial} {\partial x_j} \left( \overline{\rho D_k  \frac{\partial Y_k} {\partial x_j} } -
J_j \right)
+ \overline{\dot \omega_k}

where  J_j is the transport of subgrid fluctuations of mass fraction

J_j = \widetilde{u_jY_k} - \widetilde{u}_j \widetilde{Y}_k

and has to be modelled. Fluctuations of diffusion coefficients are often ignored and their contributions much smaller than apparent turbulent diffusion due to transport of subgrid fluctuations. The first term on the right hand side is then

\frac{\partial} {\partial x_j}  
\overline{ \rho D_k  \frac{\partial Y_k}  {\partial x_j} }
\frac{\partial} {\partial x_j}  
\overline{\rho} D_k  \frac{\partial \widetilde{Y}_k} {\partial x_j}  

Infinitely fast chemistry

All combustion models can be divided into two main groups according to the assumptions on the reaction kinetics. We can either assume the reactions to be infinitely fast - compared to e.g. mixing of the species, or of the comparable time scale of the mixing process. The simpler approach assuming chemistry fast enough, that the limiting process is mixing of the species is historically older approach and even today can be appropriate approach. It is simpler to solve then #Finite rate chemistry models, but introduces errors to the solution which may or may not be important.

Premixed Combustion

Premixed flame occurs in mixtures of fuel and oxidiser, homogeneously premixed prior to the flame. These flames are not limited only to gas fuels, but also to the pre-vaporised fuels. Typical example of premixed laminar flame is bunsen burner, where the air enters the fuel stream. The mixture burns in the wake of the riser tube walls forming nice stable flame. The premixed flames has many advantages in terms of control of temperature and products and pollution concentration, but introduce also some dangers like the autoignition (in the supply system).

Turbulent flame speed model

Eddy Break-Up model

The Eddy Break-Up model is the typical example of mixed-is-burnt combustion model. It is based on the work of Magnussen and Hjertager, and Spalding and can be found in all CFD packages. The model assumes the reactions to be completed in the moment of mixing, so that the reaction rate is completely controlled by turbulent mixing. The combustion is described by a single step global chemical reaction:

F + \nu_s O \rightarrow (1+\nu_s) P

in which F stands for fuel, O for oxidiser and P for products of the reaction. Alternativelly we can have multistep scheme, where each reaction has its own mean reaction rate. The mean reaction rate is given by:

\bar{\dot\omega}_F=A_{EB} \frac{\epsilon}{k} 

\bar{C} denotes mean concentrations for fuel, oxidiser and products respectively, A and B are model constants with typical values of 0.5 and 4.0 respectively. The values of these constants are fitted according to the experimental results and they are suitable for most of the general cases. Still they are just constants based on experimental fitting and they need not be suitable for all the situations. Care must be taken especially in highly strained regions, where the ratio of k to \epsilon is large (flame-holder wakes, walls ...). In those regions a positive reaction rate occurs and an artificial flame can be observed. CFD codes usually has some remedies to overcome this problem.

This model largely over-predicts temperatures and concentrations of species like CO and other species. Still this model is quite popular for its simplicity and relatively easy convergence and implementation.

Bray-Moss-Libby Model

Non premixed combustion

Conserved scalar equilibrium models

Finite rate chemistry

Premixed Combustion

Coherent Flame Model

Flamelets based on G equation

Non-premixed Combustion

Flamelets based on conserved scalar

Peters (2000) define Flamelets as "thin diffusion layers embedded in a turbulent non-reactive flow field". If the chemistry is fast enough, the chemistry is active within a thin region where the chemistry conditions are in (or close to) stoichiometric conditions, the "flame" surface. This thin region is assumed to be smaller than Kolmogorov length scale and therefore the region is locally laminar. The flame surface is defined as an iso-surface of a certain scalar  Z , mixture fraction in non-premixed combustion.

The reactive problems is therefore split into two parts: First, the mixing , which consists of the location of the flame surface which is a non-reactive problem concerning the propagation of a passive scalar. And second, the flame structure , which deals with the distribution of the reactive species inside the flamelet.

To obtain the distribution inside the flame front we assume it is locally one-dimensional and depends only on time and the scalar coodinate.

Using the following chain rules for the time

\frac{\partial Y_k}{\partial t} = \frac{\partial Z}{\partial t}\frac{\partial Y_k}{\partial Z}

and spatial coordinate

\frac{\partial Y_k}{\partial x_j} = \frac{\partial Z}{\partial x_j}\frac{\partial Y_k}{\partial Z}

to the species transport equation (see #Governing Equations for Reacting Flows) and re-arranging, we obtain

\rho \frac{\partial Y_k}{\partial t} + Y_k \left[  
\frac{\partial \rho}{\partial t} + \frac{\partial \rho u_j}{\partial x_j}
+ \frac{\partial Y_k}{\partial Z} \left[
 \rho \frac{\partial Z}{\partial t} + \rho u_j \frac{\partial Z}{\partial x_j} -
\frac{\partial}{\partial x_j}\left( \rho D \frac{\partial Z}{\partial x_j} \right)
 \rho D \left( \frac{\partial Z}{\partial x_j} \frac{\partial Z}{\partial x_j} \right) + \dot \omega_k

The second and third term in the LHS cancel due to continuity and mixture fraction transport, the equation therefore boils down to

\frac{\partial Y_k}{\partial t} = \frac{\chi}{2} \frac{\partial ^2 Y_k}{\partial Z^2} + \dot \omega_k

where  \chi = 2 D  \left( \frac{\partial Z}{\partial x_j} \right)^2 is called the scalar dissipation and controls the mixing, providing the interaction between the flow and the chemistry.

If the flame dependence on time is dropped, even though he field  Z  still depends on it.

 \dot \omega_k= -\frac{\chi}{2} \frac{\partial ^2 Y_k}{\partial Z^2}

This approach is called the Stationary Laminar Flamelet Model (SLFM) and has the advantage that libraries of  Y_k=f(Z,\chi) can be pre-computed and stored in look-up tables with all the required complex chemistry.

Flamelet Computation and Flamelet Libraries

The computation of non-premixed turbulent flames based on laminar-flamelet models is generally based on two-dimensional or three-dimensional CFD codes that employ standard models for fluid-mechanical closure of the govening equations. In many cases, for that purpose standard models such as the k-epsilon model are used, but occasionally more sophisticated models such as Reynolds-Stress models are also employed.

Chemical-source-term closure is a different matter. To this end, the CFD codes carry out suitable averaging procedures, such as pdf-avaraging on the basis of a beta function or a clipped Gaussian distribution. The quantities to be averaged are laminar-flamelet profiles, i.e., results from laminar-flamelet computations. Generally, these flamelet computations are carried oout a-priori, i.e, they are performed separately and prior to the turbulent-combustion simulation with the CFD code. Depending on the specific laminar-flamelet model used for the turbulent-combustion simulation, one or several parameters are varied in the laminat computatations. For instance, if the computations are based on

\frac{\partial Y_k}{\partial t} = 
\frac{\chi}{2} \frac{\partial ^2 Y_k}{\partial Z^2} + \dot w_k,

then the variable parameter is the scalar dissipation rate  \chi . The flamelet profiles for the various parameter values are stored in a dataset or file which is called a "flamelet library". For the generation of such libraries ready to use software is avalable such as Softpredict's Combustion Simulation Laboratory COSILAB [1] with its relevant solver RUN1DL, which can be used for a variety of relevant geometries; see various publications that are available for download.

Flamelets in turbulent combustion

In turbulent flames the interest is  \widetilde{Y}_k . In flamelets, the flame thickness is assumed to be much smaller than Kolmogorov scale and obviously is much smaller than the grid size. It is therefore needed a distribution of the passive scalar within the cell.  \widetilde{Y}_k cannot be obtained directly from the flamelets library  \widetilde{Y}_k \neq Y_F(Z,\chi) , where  Y_F(Z,\chi) corresponds to the value obtained from the flamelets libraries. A generic solution can be expressed as

\widetilde{Y}_k= \int  Y_F( \widetilde{Z},\widetilde{\chi}) P(Z,\chi) dZ d\chi

where  P(Z,\chi) is the joint Probability Density Function (PDF) of the mixture fraction and scalar dissipation which account for the scalar distribution inside the cell and "a priori" depends on time and space.

The most simple assumption is to use a constant distribution of the scalar dissipation within the cell and the above equation reduces to

\widetilde{Y}_k= \int Y_F(\widetilde{Z},\widetilde{\chi}) P(Z) dZ

 P(Z) is the PDF of the mixture fraction scalar and simple models (such as Gaussian or a beta PDF) can be build depending only on two moments of the scalar mean and variance, \widetilde{Z},Z''.

If the mixture fraction and scalar dissipation are consider independent variables,  P(Z,\chi) can be written as  P(Z) P(\chi). The PDF of the scalar dissipation is assumed to be log-normal with variance unity.

\widetilde{Y}_k= \int Y_F(\widetilde{Z},\widetilde{\chi}) P(Z) P(\chi) dZ d\chi

In Large eddy simulation (LES) context (see #LES equations for reacting flow), the probability density function is replaced by a subgrid PDF  \widetilde{P}. The same equation hold by replacing averaged values with filtered values.

\widetilde{Y}_k= \int Y_F(\widetilde{Z},\widetilde{\chi}) \widetilde{P}(Z) \widetilde{P}(\chi) dZ d\chi

The assumptions made regarding the shapes of the PDFs are still justified. In LES combustion the subgrid variance is smaller than RANS counterpart (part of the large-scale fluctuations are solved) and therefore the modelled PDFs are thinner.

Unsteady flamelets

Conditional Moment Closure (CMC)

In Conditional Moment Closure (CMC) methods we assume that the species mass fractions are all correlated with the mixture fraction (in non premixed combustion).

From Probability density function we have

\overline{Y_k}= \int <Y_k|\eta> P(\eta) d\eta

where  \eta is the sample space for  Z .

CMC consists of providing a set of transport equations for the conditional moments which define the flame structure.

Experimentally, it has been observed that temperature and chemical radicals are strong non-linear functions of mixture fraction. For a given species mass fraction we can decomposed it into a mean and a fluctuation:

Y_k= \overline{Y_k} + Y'_k

The fluctuations  Y_k' are usually very strong in time and space which makes the closure of  \overline{\omega_k} very difficult. However, the alternative decomposition

Y_k=  <Y_k|\eta> + y'_k

where  y'_k is the fluctuation around the conditional mean or the "conditional fluctuation". Experimentally, it is observed that  y'_k<< Y'_k , which forms the basic assumption of the CMC method. Closures. Due to this property better closure methods can be used reducing the non-linearity of the mass fraction equations.

The Derivation of the CMC equations produces the following CMC transport equation where  Q \equiv <Y_k|\eta> for simplicity.

\frac{ \partial Q}{\partial t} + <u_j|\eta>  \frac{\partial Q}{\partial x_j} =
\frac{<\chi|\eta> }{2} \frac{\partial ^2 Q}{\partial \eta^2} + 
\frac{ < \dot \omega_k|\eta> }{ <\rho| \eta >}

In this equation, high order terms in Reynolds number have been neglected. (See Derivation of the CMC equations for the complete series of terms).

It is well known that closure of the unconditional source term   \overline {\dot \omega_k} as a function of the mean temperature and species ( \overline{Y}, \overline{T}) will give rise to large errors. However, in CMC the conditional averaged mass fractions contain more information and fluctuations around the mean are much smaller. The first order closure 
 < \dot \omega_k|\eta> \approx \dot \omega_k \left( Q, <T|\eta> \right)
is a good approximation in zones which are not close to extinction.

Second order closure

A second order closure can be obtained if conditional fluctuations are taken into account. For a chemical source term in the form  \dot \omega_k = k Y_A Y_B with the rate constant in Arrhenius form  k=A_0 T^\beta exp [-Ta/T] the second order closure is (Klimenko and Bilger 1999)

< \dot \omega_k|\eta> \approx  < \dot \omega_k|\eta >^{FO}

\left[1+ \frac{< Y''_A Y''_B |\eta>}{Q_A Q_B}+ \left( \beta + T_a/Q_T \right)
\frac{< Y''_A T'' |\eta>}{Q_AQ_T} + \frac{< Y''_B T'' |\eta>}{Q_BQ_T}
\right) + ...


where  < \dot \omega_k|\eta >^{FO} is the first order CMC closure and  Q_T \equiv <T|\eta> . When the temperature exponent  \beta or  T_a/Q_T are large the error of taking the first order approximation increases. Improvement of small pollutant predictions can be obtained using the above reaction rate for selected species like CO and NO.

Double conditioning

Close to extinction and reignition. The conditional fluctuations can be very large and the primary closure of CMC of "small" fluctuations is not longer valid. A second variable  h can be chosen to define a double conditioned mass fraction

Q(x,t;\eta,\psi) \equiv <Y_i(x,t) |Z=\eta,h=\psi  >

Due to the strong dependence on chemical reactions to temperature,  h is advised to be a temperature related variable (Kronenburg 2004). Scalar dissipation is not a good choice, due to its log-normal behaviour (smaller scales give highest dissipation). A must better choice is the sensible enthalpy or a progress variable. Double conditional variables have much smaller conditional fluctuations and allow the existence of points with the same chemical composition which can be fully burning (high temperature) or just mixing (low temperature). The range of applicability is greatly increased and allows non-premixed and premixed problems to be treated without ad-hoc distinctions. The main problem is the closure of the new terms involving cross scalar transport.

The double conditional CMC equation is obtained in a similar manner than the conventional CMC equations

LES modelling

In a LES context a conditional filtering operator can be defined and  Q therefore represents a conditionally filtered reactive scalar.

Multiple Mapping Closure (MMC)

Linear Eddy Model

The Linear Eddy Model (LEM) was first developed by Kerstein(1988). It is an one-dimensional model for representing the flame structure in turbulent flows.

In every computational cell a molecular, diffusion and chemical model is defined as

\frac{\partial}{\partial t} \left( \rho Y_k \right)  =
\frac{\partial}{\partial \eta} \left( \rho D_k \frac{\partial Y_k}{\partial \eta }\right)+ \dot \omega_k

where  \eta is a spatial coordinate. The scalar distribution obtained can be seen as a one-dimensional reference field between Kolmogorov scale and grid scales.

In a second stage a series of re-arranging stochastic event take place. These events represent the effects of a certain turbulent structure of size  l , smaller than the grid size at a location  \eta_0  within the one-dimensional domain. This vortex distort the  \eta field obtain by the one-dimensional equation, creating new maxima and minima in the interval  (\eta_0, \eta + \eta_0) . The vortex size  l is chosen randomly based on the inertial scale range while  \eta_0  is obtained from a uniform distribution in  \eta  . The number of events is chosen to match the turbulent diffusivity of the flow.

PDF transport models

Probability Density Function (PDF) methods are not exclusive to combustion, although they are particularly attractive to them. They provided more information than moment closures and they are used to compute inhomegenous turbulent flows.

PDF methods are based on the transport equation of the joint-PDF of the scalars. Denoting  P \equiv P(\underline{\psi}; x, t) where  \underline{\psi} = ( \psi_1,\psi_2 ... \psi_N) is the phase space for the reactive scalars  \underline{Y} = ( Y_1,Y_2 ... Y_N) . The transport equation of the joint PDF is:

\frac{\partial <\rho | \underline{Y}=\underline{\psi}> P }{\partial t} + \frac{  
\partial <\rho u_j | \underline{Y}=\underline{\psi}> P }{\partial x_j} =
\sum^N_\alpha \frac{\partial}{\partial \psi_\alpha}\left[ \rho \dot{\omega}_\alpha P \right]
- \sum^N_\alpha \sum^N_\beta \frac{\partial^2}{\partial \psi_\alpha \psi_\beta}
\left[ <D \frac{\partial Y_\alpha}{\partial x_i} \frac{\partial Y_\beta}{\partial x_i} | \underline{Y}=\underline{\psi}> \right] P

where the chemical source term is closed. Another term appeared on the right hand side which accounts for the effects of the molecular mixing on the PDF, is the so called "micro-mixing " term. Equal diffusivities are used for simplicity  D_k = D

A more general approach is the velocity-composition joint-PDF with  P \equiv P(\underline{V},\underline{\psi}; x, t) , where  \underline{V} is the sample space of the velocity field  u,v,w . This approach has the advantage of avoiding gradient-diffusion modelling. A similar equation to the above is obtained combining the momentum and scalar transport equation.

The PDF transport equation can be solved in two ways: through a Lagrangian approach using stochastic methods or in a Eulerian ways using stochastic fields.


The main idea of Lagrangian methods is that the flow can be represented by an ensemble of fluid particles. Central to this approach is the stochastic differential equations and in particular the Langevin equation.


Instead of stochastic particles, smooth stochastic fields can be used to represent the probability density function (PDF) of a scalar (or joint PDF) involved in transport (convection), diffusion and chemical reaction (Valino 1998). This method is purely Eulerian and offers implementations advantages compared to Lagrangian or semi-Eulerian methods. Transport equations for scalars are often easy to programme and normal CFD algorithms can be used (see Discretisation of convective term)

A new set of  N_s scalar variables (the stochastic field  \xi ) is used to represent the PDF

P (\underline{\psi}; x,t) = \frac{1}{N} \sum^{N_s}_{j=1} \prod^{N}_{k=1} \delta \left[\psi_k -\xi_k^j(x,t) \right]


  • Griffiths, J. F. (1994), "Reduced Kinetic Models and Their Application to Practical Combustion Systems", Prog. in Energy and Combustion Science,Vol. 21, pp. 25-107.
  • Peters, N. (2000), Turbulent Combustion, Cambridge University Press.
  • Poinsot, T.,Veynante, D. (2001), Theoretical and Numerical Combustion, ISBN 1-930217-05-6, R. T Edwards.
  • Westbrook, Ch. K., Dryer,F. L., (1984), "Chemical Kinetic Modeling of Hydrocarbon Combustion", Prog. in Energy and Combustion Science,Vol. 10, pp. 1-57.

External links and sources

My wiki