5TH EUROPEAN CONFERENCE FOR AERONAUTICS AND SPACE SCIENCES (EUCASS)
Large Eddy Simulation of Kerosene/LOx Supercritical
Mixing Characteristics of a Swirl Injector using Various
Cubic Equation of State
Junyoung Heo*, Kukjin Kim*, and Hong-Gye Sung*
*Korea Aerospace University
Goyang Gyeonggi, 412-791, South Korea
Abstract
The turbulent mixing of a kerosene and liquid oxygen coaxial swirl injector under supercritical
pressures has been numerically investigated in three dimensional system. Turbulent numerical model
is based on LES (Large Eddy Simulation) with real-fluid transport and thermodynamics over the entire
pressure range; Cubic Equations of State (SRK, PR, and RK-PR), Chung’s model for
viscosity/conductivity, and Fuller’s theorem for diffusivity to take account Takahashi’s compressible
effect. The effects of real-gas equation of state and swirl intensity on flow structures in the swirl
injector are investigated and power spectral densities of pressure fluctuations and mixing efficiency of
the injector are analyzed.
Nomenclature
D
E
f
p
q
t
u
x
Z
=
=
=
=
=
=
=
=
=
Diffusivity
Specific total energy
Mixture fraction
Pressure
Heat flux
Time
Velocity components
Spatial coordinate
Pseudo-time variable vector
=
=
=
=
Kronecker delta
Preconditioning matrix or circulation
Density
Viscous stress tensor
Greek




Superscripts
sgs
= Resolved-scale
= Favre-averaged resolved-scale
= Subgrid-scale
1. Introduction
An advanced liquid rocket engine is operating at high pressure in order to increase combustion efficiency and
thrust, thus the propellants are exposed to supercritical conditions in the combustion chamber. The supercritical
conditions are the state that arbitrary substance is enclosed by the environment over its critical point. In such
Copyright  2013 by Heo, J. Y., Kim, K. J., and Sung, H. G. Published by the EUCASS association with permission.
Heo, J. Y., Kim, K. J., and Sung, H. G.
conditions, the propellants experience the transition of thermophysical properties different from phase change
process at subcritical conditions. After the liquid propellant is introduced into the injector and combustion chamber,
its properties are varied continuously and the distinction between gas and liquid is ambiguous. At near the critical
point, propellant properties exhibit gas-like diffusivities and liquid-like densities so that the engine performance can
be achieved with high energy density. Hence, the study for the mixing and combustion process of the propellants at
supercritical conditions is an essential subject. Because of the peculiarities of substances at supercritical conditions,
the conventional methods for property analysis are inappropriate so that the applicable thermodynamic relationships
for a wide range of pressure must be employed.
For understanding of the atomization, mixing, and combustion phenomena in advanced liquid rocket engines,
experimental and numerical studies are performed in the wide range of pressure and temperature including the
supercritical conditions. In research of numerical analysis, the mixing and combustion processes in shear coaxial
injectors using propellants as methane/liquid oxygen and hydrogen/liquid oxygen at supercritical conditions have
been investigated. The thermophysical characteristics of hydrogen and liquid oxygen flame in a shear coaxial injector
at supercritical pressure have been investigated by Oefelein using direct numerical simulation [1]. Zong and Yang
applied the large eddy simulation for the identification of cryogenic fluid dynamics in a swirl oxidizer injector
operating at supercritical pressure [2]. Giorgi and Leuzzi performed the numerical analysis of mixing and
combustion in methane and liquid oxygen shear coaxial injector at supercritical conditions [3]. Recently, Giorgi et al.
simulated combustion characteristics of methane and liquid oxygen and compared between the results in accordance
with combustion models, reaction mechanisms [4], and equations of state at supercritical conditions. Sierra et al.
numerically investigated reaction mechanisms and combustion characteristics for the recess number using hydrogen
and liquid oxygen as propellants [5]. Hiroshi et al. analysed nitrogen injection characteristics and compared between
the numerical results to pressure and temperature in combustor [6]. H. Huo et al. simulated combustion
characteristics in a kerosene and liquid oxygen supercritical injector applying the flamelet library technique as a
combustion model [7].
In this study, on the basis of previous study [8], the turbulent mixing of a kerosene and liquid oxygen coaxial swirl
injector under supercritical pressures has been numerically investigated in three dimensional system. Turbulent
numerical model is based on LES with real-fluid transport and thermodynamics; SRK (Soave modification of
Redlich-Kwong) [9], PR (Peng-Robinson) [10], RK-PR (Redlich-Kwong-Peng-Robinson) [11] equations of state,
Chung’s model for viscosity/conductivity, and Fuller’s theorem for diffusivity to take account Takahashi’s
compressible effect. The effect of real-gas equation of state and swirl intensity are investigated and power spectral
densities of pressure fluctuations and mixing efficiency are analyzed.
2. Numerical Method
2.1 Filtered Transport Equations
The theoretical formulation is based on the filtered Favre averaged mass, momentum, energy, and mixture fraction
conservation equations in Cartesian coordinates. Turbulent closure is achieved by means of a LES technique in
which large scale turbulent structures are directly computed and the unresolved small scale structures are treated by
using the analytic or empirical modeling. The governing equations can be written as:
  u j

0
t
x j
(1)
sgs
sgs
 ui  (  ui u j  p ij )  ( ij   ij  Dij )


t
x j
x j
(2)
 E  (  E  p)ui 



 qi  u j ij  Qisgs  Hisgs   ijsgs 
t
xi
xi
(3)

 f  (  ui f )
 
f


  isgs 
D
t
xi
xi 
xi

(4)
The unclosed sub-grid scale terms are defined
2
5TH EUROPEAN CONFERENCE FOR AERONAUTICS AND SPACE SCIENCES (EUCASS)
 ij sgs    ui u j  ui u j 
(5)
Dij sgs   ij   ijj 
(6)
Qi sgs   qij  qijj 
(7)

Hi sgs    Eui  Eui   pui  pui


(8)
 i sgs  u j ij  u j ijj

(9)


(10)
i sgs   ui f  ui f
 ij sgs , Dij sgs , Qij sgs , H ij sgs ,  ij sgs and  ij sgs are the sgs stress, nonlinearity of viscous stress term, heat flux, energy
flux, viscous work, and conserved scalar flux, respectively. The dynamic Smagorinsky model is employed to close
those sgs terms [12].
2.2 General-Fluid Thermodynamics and Transport
Intermolecular forces and the molecular body volume are considered as important factors in the real-fluid. In order
to understand the real-fluid phenomenon, effects such as compressibility, van der Waals force, molecular
dissociation, and non-equilibrium thermodynamic quantity must be taken into account. An equation of state capable
of handling dynamics of real-fluid at high pressure and low temperature regime is required for thermophysical
property evaluation of real-fluid. The most common equations of state utilized for calculating high pressure fluid
properties are the cubic equations of state, which include SRK, PR, and RK-PR equations of state.
SRK equation of state:
p
 Ru T
a
2
M w  M w  b 
(11)
a 2
 M w  b  M w  b    2b  M w
(12)
a 2
 M w  1b  M w   2b 
(13)
 M w  b 

PR equation of state:
p
 RuT
 M w  b 

RK-PR equation of state:
p
 Ru T
 M w  b 

2.3 Numerical Scheme
The governing equations are numerically solved by means of a finite volume method. This method allows for the
treatment of arbitrary geometry. The spatial discretization employs a fourth order central differencing scheme in
generalized coordinates. The temporal discretization is obtained by second order backward differencing scheme
using a fourth order Runge-Kutta scheme for integration of the real time term. In the regions of low Mach number
flows, the energy and momentum equations are practically decoupled and the system of conservation equations
becomes stiff. Pressure decomposition and preconditioning techniques with dual time stepping are applied to
circumvent the round-off errors associated with the calculation of the pressure gradient and the convergence
difficulties for the low Mach number flows in the momentum equation. First, a rescaled pressure term is used in the
momentum equation to circumvent the singular behavior of pressure at low Mach numbers. Second, a
3
Heo, J. Y., Kim, K. J., and Sung, H. G.
preconditioning technique with dual time-stepping integration procedure is established. A unified treatment of
general fluid thermodynamics is incorporated into a preconditioning scheme. The pseudo-time derivative may be
chosen to optimize the convergence of the inner iterations through the use of an appropriate preconditioning matrix
that is tuned to rescale the eigenvalues to render the same order of magnitude to maximize convergence. To unify the
conserved flux variables, a pseudo-time derivative of the form Z /  can be added to the conservation equation.
Since the pseudo-time derivative term disappears as the solution is converged, a certain amount of liberty can be
taken in choosing the variable, Z. We take advantage of this by introducing a pressure p as the pseudo-time
derivative term in the continuity equation. The code is paralleled using an MPI library for more effective calculation.
3. Model Description
Figure 1 shows a coaxial swirl injector geometry and the swirl numbers of liquid oxygen and kerosene are 1.598
and 6.548, respectively. The injector includes three major parts: a tangential inlet, a vortex chamber, and a discharge
nozzle. The kerosene (outer) and liquid oxygen (inner) injects into the swirl injector through the tangential passage.
In case of quasi-3D simulation, the radial and azimuthal velocities are determined from the tangential inlet port angle
and these velocity components are factor to set the swirl strength and mass flow rate, respectively.
LOx
Fuel
dt,f or dt,o
Lr
dn,o
So 
4d s , o d n , o
 nn,o dt ,o
2
 1.598 ,
doo
Sf 
dn,f
4  d n, f  do,o  d s , f
 n f dt , f 2
ds,f or ds,o
 6.548
Figure 1: Injector geometry and swirl number
Table 1: Operation conditions
Injector
Chamber Pressure
Fuel
Oxidizer
Mass Flow Rate
Fuel
Coaxial Swirl
100 bar
Kerosene; 350 K
Liquid Oxygen; 103 K
0.084 kg/s
Oxidizer
0.232 kg/s
Table 2: Computation conditions
Case
1
2
3
4
5
4
Equation of State
SRK
PR
RK-PR
RK-PR
RK-PR
Dimension
Quasi-3D
Quasi-3D
Quasi-3D
Quasi-3D
Full 3D
Tangential Velocity
x1
x1
x1
x2
Determined by Inlet Geometry
5TH EUROPEAN CONFERENCE FOR AERONAUTICS AND SPACE SCIENCES (EUCASS)
Figure 2: Computational domain (Quasi-3D)
Figure 3: Computational domain (Full 3D)
Operation and computation conditions are listed in Table 1 and 2, respectively. Quasi-3D flow fields with periodic
boundaries is treated herein because of the enormous computational time for simulating of full 3D using the large
eddy simulation. The tangential inlets which are circular ports connected to the injector are approximated with a thin
slit. Figure 2 shows a computational domain divided into 45 blocks, with each calculated on a single processor of a
distributed computing facility. The total grid points are about 765 000 nodes. The computational domain of full 3D
calculation is shown as Fig. 3 and grid points are 6 561 000 nodes which is divided into 146 blocks. Fuel and Oxygen
are supplied through the projected faces like Fig. 3. The disturbances at all of calculation conditions are generated by
a Gaussian random-number generator with an intensity of 5% of the mean quantity.
4. Results and Discussion
SRK and PR equations of state usually have been used to predict the thermophysical properties of real-gas. The
SRK equation of state has high accuracy for low molecular weight compounds while the PR equation of state does
for high weight. Considering kerosene and liquid oxygen used as propellants, the SRK equation of state can predict
the supercritical oxygen properties properly and the PR equation of state is accurate for kerosene properties at
supercritical conditions. However, the RK-PR equation of state can be properly applied to simulate the properties of
both kerosene and liquid oxygen [13-14]. The fuel and oxidizer having the only tangential and azimuthal velocity
components are supplied for the quasi-3D calculation of the swirl coaxial injector. In that calculation, the tangential
and azimuthal velocity determine the swirl strength and mass flow rate, respectively. The swirl strength is controlled
by changing the tangential velocity in this study.
5
Heo, J. Y., Kim, K. J., and Sung, H. G.
4.1 The Effect of Cubic Equation of State
(a) Kerosene
(b) Oxygen
Figure 4: Density distributions for different equation of state at 100 bar
(a) SRK
(b) PR
(c) RK-PR
Figure 5: Instantaneous density contours according to real-gas equation of state
Figure 4 represents density distributions to different equation of state as compared to that of the principle of ECS
(Extended Corresponding States) at 100 bar. The principal of ECS is that the properties of any interest fluid are
represented by scaling that of the reference fluid and may also be functions of density and temperature. The principle
of ECS is a very accurate method for prediction of thermophysical properties, applicable to the entire range of fluid
states, from dense liquid to dilute gas, as well as to the supercritical fluid conditions [15-17]. As shown in the results,
the RK-PR equation of state provides the most suitable results for kerosene and oxygen in the wide range of
temperature.
Figure 5 shows comparison results of density fields to the defferent equations of state with the same inlet
conditions of fuel and oxidizer. The SRK and PR equations of state predict lower kerosene density and higher
oxygen density respectively. On the other hand, the RK-PR equation of state has very high accuracy in the prediction
results of both kerosene and liquid oxygen. This study also shows that the SRK and PR equations of state have lower
density for kerosene and higher density for liquid oxygen, respectively, than the results by the RK-PR equation of
state. As a result, it is investigated that the liquid oxygen core sustains its length as long and the CTRZ (Central
Toroidal Recirculation Zone) is formed in the rear region compared with the other results in the results applied the
PR equation of state due to the high density of liquid oxygen. Also, the mixing layer of the fuel/oxidizer gets longer
and the CTRZ possesses the wider area since the kerosene density results by the RK-PR equations of state are higher
than those by the SRK equation of state.
6
5TH EUROPEAN CONFERENCE FOR AERONAUTICS AND SPACE SCIENCES (EUCASS)
4.2 Swirl Intensity Effect
Two tangential velocities are applied to investigate the swirl strength effect on mixing characteristics. Figure 6
represents the density fields for each swirl strength. The formed liquid oxygen film is closer to the injector wall since
a centrifugal force of swirl flow becomes stronger as the tangential velocity increases. As the swirl strength is
stronger, the gas phase core is formed in the central region of oxidizer injector which affects the film boundary area
enlarged in the injector and facilitates vorticity breakup.
Figure 6: Instantaneous density contour according to swirl intensity
Figure 7: Mixing efficiency at different swirl intensity
To quantify the level of mixing, mixing efficiency is considered using the following formula:
mix
 Y
( x) 
 Y
Fuel
u dA
Fuel
udA
1

 
1
  /
global
 local
( local  global )
( local  global )
(14)
7
Heo, J. Y., Kim, K. J., and Sung, H. G.
Where ρ and YFuel are density and mass fraction of kerosene, respectively. The u is the velocity component normal
to the area, dA , local and global are local and global equivalence ratio, respectively. The parameter  is decided by
local equivalence ratio. This parameter means that if the local equivalence ratio is equal or less than global
equivalence ratio, the mixing is regarded as completed. On the other hand, local efficiency is less than 1 when local
equivalence ratio is large. Therefore, the mixing efficiency represents how much fuel spread into combustion
chamber, in other words, the efficiency becomes 1 when local equivalence ratio is equal to global equivalence ratio
through whole cross sectional area. Figure 7 shows mixing efficiencies on the time averaged fields. The efficiency
increases as moves to downstream after the oxidizer injector post where the mixing starts. As shown in Fig. 7, the
stronger swirl intensity is, the earlier the mixing efficiency reaches 1, meaning that most regions become lean burn
environment because the mixing completely takes place not much far from the oxidizer injector post. Then, when
increasing the swirl strength, the film is formed inside the oxidizer injector by the strong swirl flow, rapid mixing of
fuel and the oxidizer is performed.
Figure 8: Probe location
Figures 8 shows the locations of pressure measurement: inner side of injector (probe 1, 2), near the oxidizer
injector post (probe 3), phase-change region (probe 4), mixing layer (probe 5), and front of the CTRZ (probe 6).
Figure 9 represents the power spectral densities of the pressure fluctuations at the locations. The Fast Fourier
Transform (FFT) technique was implemented for the spectral analysis.
(a) Case 3
(b) Case 4
Figure 9: Power spectral densities of pressure fluctuations according to swirl intensity
As shown in case 3, a frequency of 1.05 kHz is dominantly estimated inside the oxidizer injector filled with high
density fluid and its amplitude is higher than the other locations. Frequencies of 1.05 and 2.75 kHz are observed at
the other locations too. The oscillation at 2.75 kHz is same as a cycle of the vortex shedding at the oxidizer injector
exit. Also, high amplitude of the pressure oscillation at 2.75 kHz is observed in kerosene injector due to the Doppler
8
5TH EUROPEAN CONFERENCE FOR AERONAUTICS AND SPACE SCIENCES (EUCASS)
effects. In case 4, a frequency of 2.23 kHz which is same as a cycle of the vortex shedding at the tip of oxidizer
injector post is dominantly estimated for all locations while the low frequency observed inside the oxidizer injector
did not occurred as compared with those of case 3. The initial turbulence flow at the tangential inlet generates a
succession of wavy structures and is convected downstream through the dense liquid oxygen. In a swirl injector,
those waves are generated by the presence of oscillations in the chamber or injector post. When the swirl intensity is
strong, the core boundary (film) is formed in propellant feed-line. Also, gasification of the liquid oxygen is promoted
from the core boundary and the mixing is rapidly performed in order of the vortex formed at the interface of the
liquid phase. The pressure wave of the low frequency is not generated inside oxidizer injector due to the film.
4.3 3D Effect
The swirl injector flow structures computed in full 3D are presented in Fig. 10. Liquid oxygen and kerosene are
introduced into the injector through slots. Since the inlet ports consist of the tangential slit, the strong radial pressure
gradient is formed by centrifugal force. Because the chamber pressure is higher than the critical pressure, there is no
obvious boundary between the injected fluid and the ambience. Near the oxidizer injector post, a slowly swirling
gaseous region is formed. The adverse pressure gradient in this region gives rise to recirculation flow and finally
produces a central toroidal recirculation zone (CTRZ), resulting from vortex breakdown, which plays an important
role in determining the flame stabilization characteristics. The CTRZ accelerates the earlier burst of vortices near the
fuel/oxygen mixing zone. In the CTRZ, these vortical structures move downstream and upstream, and are rapidly
dissipated by the turbulent diffusion and viscous damping. In the full 3D results, unlike the quasi-3D calculation, the
core of the gas phase is formed in the oxidizer injector.
Pressure oscillations in the chamber are reflected in the gas phase region inside injector interacted acoustically
with the chamber, which deforms the boundary of the liquid and gaseous phases. As shown in Fig. 10, gasification of
the liquid oxygen is promoted by the pressure wave transmitted to the inside of the injector, the dynamic flow of the
injection is amplified in the formation of vortices at the film boundary. It may cause a significant effect on the
mixing structure in the chamber.
Figure 10: Flow structure of coaxial swirl injector
Figure 11 represents the temporal evolution of pressure and temperature field over one cycle of the dominant
pressure oscillation in oxidizer injector. As the propellants are discharged into the chamber, sinuous KelvinHelmholtz waves develop along the outer wall of the injector. As these waves move downstream, a hairpin vortice is
generated and grows up. The oscillation is same as a cycle of the vortex shedding caused by the hairpin vortices at
the LOx post.
The pressure fluctuations from rotating in the tangential direction by the full 3D effect, these pressure waves affect
the change of film thickness in the tangential direction. Similar to the quasi-3D results, fuel - and oxidizer injectors
9
Heo, J. Y., Kim, K. J., and Sung, H. G.
have a dominant frequency about 1.4 kHz and 3.4 kHz, respectively, i.e. the dominant frequency in the fuel injector
is higher than a frequency in the oxidizer injector due to the thin flow passage.
θ=0˚
θ = 120 ˚
θ = 240 ˚
(a) Pressure
(b) Temperature
Figure 10: Temporal evolution of pressure and temperature field over one cycle of the dominant pressure oscillation
in the oxygen injector
5. Conclusions
The mixing characteristics of co-axial swirl injector are investigated according to real-gas equation of state and
swirl intensity under supercritical condition. The RK-PR equation of state provides more reasonable flow and mixing
structures than the SRK, PR equation of state. In order to compare the mixing efficiency according to the swirl
intensity, inlet conditions are controlled by different tangential velocities. The case increasing tangential velocity up
to twice reaches mixing efficiency 1 about 12% ahead of a baseline case. As the swirl intensity increases, the gaseous
core is obviously formed inside oxidizer injector. Also, gasification of the liquid oxygen is facilitated from the core
boundary so that the low frequency, 1.05~1.4 kHz, inside oxidizer injector disappears due to the active oscillation on
the film boundary. Full 3D simulation has been calculated in order to apply accurate tangential inlet ports. The gasphase core is formed in the oxidizer injector unlike the quasi-3D calculation. The vertical structures are amplified at
the core interface, causing a significant effect on the mixing structure in the chamber.
10
5TH EUROPEAN CONFERENCE FOR AERONAUTICS AND SPACE SCIENCES (EUCASS)
Acknowledgement
This research was supported by National Space Laboratory (NSL) Program (No. 2008-2006287) through the
National Research Foundation of Korea (NRF) funded by the Ministry of Education, Science and Technology.
References
[1] Oefelein, J. C. 2005. Thermophysical Characteristics of Shear Coaxial LOx-H2 Flames at Supercritical Pressure.
Proceedings of the Combustion Institute. Vol. 30. Issue 2. 2929-2937.
[2] Zong, N. and Yang, V. 2008. Cryogenic Fluid Dynamics of Pressure Swirl Injectors at Supercritical Conditions.
Physics of Fluids. Vol. 20. Issue 5. 056103-1-14.
[3] Giorgi, M. G. D. and Leuzzi, A. 2009. CFD Simulation of Mixing and Combustion in LOx/CH4 Spray Under
Supercritical Conditions. 39th AIAA Fluid Dynamics Conference.
[4] De Giorgi, M., G., Sciolti, A., Ficarella, A. 2011. Comparisons between Different Combustion Models for High
Pressure LOx/CH4 Jet Flame. AIAA Fluid Dynamics Conference and Exhibit.
[5] Sierra, P., Masquelet, M., Menon, S. 2013. Large-Eddy Simulation of a Reactive Shear Coaxial Injector
Configuration. AIAA Aerospace Sciences Meeting.
[6] Hiroshi, T., Mitsuo, K. 2013. Characterization of Cryogenic Nitrogen Jet Mixings under Supercritical Pressures.
AIAA Aerospace Sciences Meeting.
[7] Huo, H., Yang, V. 2013. Large-Eddy Simulation of Supercritical Combustion of Liquid Oxygen and Kerosene
of a Bi-Swirl Injector. AIAA Aerospace Sciences Meeting.
[8] Heo, J. Y., Kim, K. J., Sung, H. G. 2011. Numerical Study on Kerosene/LOx Supercritical Mixing
Characteristics of a Swirl Injector. AIAA Aerospace Sciences Meeting.
[9] Soave, G. 1972. Equilibrium Constants from a Modified Redlich-Kwong Equation of State. Chemical
Engineering Science. Vol. 27. Issue 6. 1197-1203.
[10] Peng, D. and Robinson D. B. 1976. A New Two-Constant Equation of State. Industry Engineering Chemical
Fundamentals. Vol. 15. No. 1. 59-64.
[11] Cismondi, M. and Mollerup, J. 2005. Development and Application of a Three-Parameter RK-PR Equation of
State. Fluid Phase Equilibria. Vol. 232. Issue 1-2. 74-89.
[12] Germano, M. 1991. A Dynamic Subgrid-scale Eddy Viscosity Model. Phys. Fluids A. Vol. 3. No. 7. 1760-1765.
[13] Kim, S. K., Choi, H. S., Kim, Y. M. 2012. Thermodynamic Modeling Based on a Generalized Cubic Equation
of State for Kerosene/LOx Rocket Combustion. Combustion and Flame. Vol. 159. Issue 3. 1351-1365.
[14] Kim, K. J., Heo, J. Y., Sung, H. G. 2012. Study on Thermophysical Property Characteristics of Kerosene
Surrogates in a Swirl Injector at Supercritical Pressure Condition. The Korean Society of Propulsion Engineers
Spring Conference.
[15] Ely, J. F. and Hanley, H. J. M. 1981. Prediction of Transport Properties. 1. Viscosity of Fluids and Mixtures.
Industrial and Engineering Chemistry Fundamentals. Vol. 20. No. 4. 323-332.
[16] Ely, J. F. and Hanley, H. J. M. 1983. Prediction of Transport Properties. 2. Thermal Conductivity of Pure Fluids
and Mixtures. Industrial and Engineering Chemistry Fundamentals. Vol. 22. No. 1. 90-97.
[17] Yang, J. C., Vazquez, I., Boyer, C. I., Huber, M. L., and Weber, L. 1997. Measured and Predicted
Thermodynamic Properties of Selected Halon Alternative/Nitrogen Mixtures. International Journal of
Refrigeration. Vol. 20. No. 2. 96-105.
11