close

Вход

Забыли?

вход по аккаунту

?

170.Журнал Сибирского федерального университета. Сер. Техника и технологии №1 2012

код для вставкиСкачать
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Æóðíàë Ñèáèðñêîãî ôåäåðàëüíîãî óíèâåðñèòåòà
2012
Journal of Siberian Federal University
5 (1)
Òåõíèêà è òåõíîëîãèè
Engineering & Technologies
Редакционный совет
академик РАН Е.А.Ваганов
академик РАН И.И.Гительзон
академик РАН А.Г.Дегерменджи
академик РАН В.Ф.Шабанов
чл.-к. РАН, д-р физ.-мат. наук
В.Л.Миронов
чл.-к. РАН, д-р техн. наук
Г.Л.Пашков
чл.-к. РАН, д-р физ.-мат. наук
В.В.Шайдуров
член-корр. РАН, д-р физ.-мат. наук
В.В. Зуев
Editorial Advisory Board
Chairman:
Eugene A. Vaganov
Members:
Josef J. Gitelzon
Vasily F. Shabanov
Andrey G. Degermendzhy
Valery L. Mironov
Gennady L. Pashkov
Vladimir V. Shaidurov
Vladimir V. Zuev
Editorial Board:
Editor-in-Chief:
Mikhail I. Gladyshev
Founding Editor:
Vladimir I. Kolmakov
Managing Editor:
Olga F. Alexandrova
Executive Editor for Engineering &
Technologies:
Vladimir A. Kulagin
CONTENTS / ÑÎÄÅÐÆÀÍÈÅ
Robert Möckel, Jens Götze, Sergey A. Sergeev,
Igor N. Kapitonov, Elena V. Adamskaya,
Nikolay A. Goltsin and Torsten Vennemann
Trace-Element Analysis by Laser Ablation Inductively Coupled
Plasma Mass Spectrometry (LA-ICP-MS): a Case Study for
Agates from Nowy Kościoł, Poland
–3–
Sergey B. Sidelnikov, Vladimir N. Baranov,
Ekaterina S. Lopatina and Leonid P. Trifonenkov
Study of Effect of Nickel and Rare-Earth Metals on the Structure
and Properties of Cast and Deformed Semi-Finished Products
Made of Aluminum Alloys
– 19 –
Viktor N. Timofeev, Àlexander I. Korchagin,
Evgeniy A. Pavlov and Nikolay V. Timofeev
Control of Convective Flows in Liquid Metal in Fluenced by
Electromagnetic Forces
– 28 –
Ivan V. Shulgin,
Aleksey A. Gerasimenko and Zhou Su Quan
Stochastic Simulation of Covariance Matrix and Power Load
Curves in Electric Distribution Networks
– 39 –
Â.À. Êóëàãèí, Ò.À. Ïüÿíûõ
h““лед%"=…,е
*=",2=ц,%……/.
м=2ем=2,че“*%г% м%дел,!%"=…,
2ече…,L
“!ед“2"=м,
– 57 –
Редактор И.А. Вейсиг Корректор Т.Е. Бастрыгина
Компьютерная верстка Е.В. Гревцовой
Подписано в печать 20.02.2012 г. Формат 84x108/16. Усл. печ. л. 9,8.
Уч.-изд. л. 9,3. Бумага тип. Печать офсетная. Тираж 1000 экз. Заказ 7868.
Отпечатано в ПЦ БИК СФУ. 660041 Красноярск, пр. Свободный, 82a.
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Editorial board for Engineering &
Technologies:
Vladimir A. Kulagin
Yury D. Alashkevich
Viktor G. Anopchenko
S. T. Batmunkh
Yury B. Galerkin
Gennadiy I. Gritsko
Georg Guggenberger
Carsten Drebenstedt
Lev V. Endjievsky
Sergey V. Kaverzin
Feng-Chen Li
Vladimir А. Makarov
Alexander V. Mineev
Vladimir V. Moskvichev
Bernard Nacke
Oleksandr F. Nemchin
Valeriy A. Nikulin
Oleg Ostrovski
Harald A. Oye
Vasiliy I. Panteleev
Sergey P. Pan’ko
Peter V. Polyakov
Anatoli M. Sazonov
Viktor N. Timofeev
Ibragim Khisameev
Anatoly Z. Shvidenko
Galina A. Chiganova
Свидетельство о регистрации СМИ
ПИ № ФС77-28-722 от 29.06.2007 г.
Серия включена в «Перечень ведущих рецензируемых научных журналов и изданий, в которых должны
быть опубликованы основные научные результаты диссертации на
соискание ученой степени доктора и
кандидата наук» (редакция 2010 г.)
Igor A. Shimanskiy,
Vladimir G. Babkin, Vladimir K. Frizorger,
Alexander S. Samoylo and Alexey B. Nabiulin
Influence of Initial Components Mechanical Activation on the
Properties of Protective Enamel Coating on the Cast-Iron
Surface
– 63 –
Olga M. Sharonova, Natalia A. Oreshkina,
Larisa I. Kurteeva and Alexandre G. Anshits
Aerodynamic Separation of Fly Ashes of Selective Sampling
from Pulverized Combustion of Coals of Different Ranks
– 72 –
Irina Yu. Botvich, Alexander F. Sidko,
Tamara I. Pisman and Anatoly P. Shevyrnogov
Determination of Chlorophyll Photosynthetic Potential in
Vegetation Using Ground-Based and Satellite Methods
– 87 –
Vladimir A. Kodnyanko
Static Characteristics of Active Hydrostatic Two-Row Radial
Bearing with Restriction of Output Lubricant Flow
– 98 –
Vladimir O. Andreev, Sergey E. Tinykov,
Oksana P. Ovchinnikova and Gennady P. Parahin
Extreme Value Theory and Peaks Over Threshold Model in the
Russian Stock Market
– 111 –
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Journal of Siberian Federal University. Engineering & Technologies 1 (2012 5) 3-18
~~~
УДК 621.315.592
Trace-Element Analysis
by Laser Ablation Inductively Coupled
Plasma Mass Spectrometry (LA-ICP-MS):
a Case Study for Agates
from Nowy Kościoł, Poland
Robert Möckela,
Jens Götzea, Sergey A. Sergeevb*,
Igor N. Kapitonovb, Elena V. Adamskayab,
Nikolay A. Goltsinb and Torsten Vennemannc
a
TU Bergakademie Freiberg, Institute of Mineralogy,
14 Brennhausgasse, D-09596 Freiberg, Germany
b
A.P. Karpinsky All Russian
Geological Research Institute (VSEGEI),
74 Sredny Prospekt, St.-Petersburg 199106, Russia
c
Institut de Minéralogie et Géochimie,
Université de Lausanne,
UNIL-BFSH2, CH-1015 Lausanne, Switzerland 1
Received 6.02.2012, received in revised form 13.02.2012, accepted 20.02.2012
Laser ablation inductively coupled plasma mass spectrometry (LA-ICP-MS) was applied to detect
trace elements in agate from Permian volcanics (Nowy Kościoł, Poland) in low concentrations and
with high spatial resolution. The used LA-ICP-MS system consists of a DUV 193 laser ablation system
linked to a Thermo Finnigan Element 2 mass spectrometer. The use of a 193 nm ArF excimer laser
(50-200 mJ energy output) and the standards NIST 611 and NIST 612 enables to produce and analyse
small crater diameters down to 5 μm.
Trace-element profiles have been analyzed for the elements Ti, Ge, Al, Fe, Mn, U, Th, Ba, Sr, Rb, Cs,
and Y in the ppm- and sub-ppm level. The concentrations of the REE are sometimes below the detection
limit of the method. Almost all elements (except Cu) display higher contents in chalcedony than in the
macrocrystalline quartz. Fe, for instance, shows a 100 times higher concentration in agate bands
compared to quartz, which may be due to finely distributed iron oxide particles in the chalcedony
which probably act as colour pigments.
The trace elements in agate are released simultaneously with Si during alteration of the surrounding
volcanic rocks. Oxygen isotope data indicate that silica accumulation and agate formation took
place at temperatures below 120°C. The characteristic trace-element distribution patterns in agate
result from a “self-purification” process during crystallization of chalcedony and quartz from a
silica gel.
*
1
Corresponding author E-mail address: vsegei@vsegei.sp.ru
© Siberian Federal University. All rights reserved
#3#
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Robert Möckel, Jens Götze…Trace-Element Analysis by Laser Ablation Inductively Coupled Plasma Mass Spectrometry…
Keywords: agate, quartz, trace elements, LA-ICP-MS, oxygen isotopes, CL, Nowy Kościόł.
Introduction
Agates are banded forms of microcrystalline α-quartz (chalcedony) and often form spectacular,
multicoloured geodes or veins. These forms of SiO2 can be found all over the world [1]. Despite the
intense research concerning the conditions of agate formation during the last 150 years, detailed
theories and suggestions are widespread and still controversely discussed (e.g. [2,3]). During the last
decades, modern analytical methods such as fluid inclusion studies, oxygen isotope and trace-element
analysis or luminescence microscopy and spectroscopy have been more involved into agate research
(e.g. [4,2,5,6,7,8,9,10,11]). Nevertheless, many of the genetic problems are still unresolved and specific
aspects have not been studied yet.
For instance, up to now only few information exist about the spatial distribution of trace elements
in agates. Pilot studies by Heaney & Davis [12] or Götze et al. [10] revealed that the trace element
patterns of chalcedony and macrocrystalline quartz can be different in the same agate. On the other
hand, rhythmic banding of agate layers can result in variations of the trace-element distribution. The
detailed knowledge about the trace-element behaviour may provide further information about the
physico-chemical conditions during agate formation.
The method of laser ablation inductively coupled plasma mass spectroscopy (LA-ICP-MS) is
available since the mid-1980s [13] and seems to suit perfectly to determine a spatial resolved traceelement distribution in agates. During the last two decades, the method was brought to perfection with
detection limits down to a few ppb, depending on analytical conditions (e.g. spot size) and sample
material [14]. Recently, LA-ICP-MS is a widely accepted powerful analytical method in science and
technology, including geoscientific research. An overview about requirements and recent technical
developments is given in Flem et al. [15], Heinrich et al. [16] or Günther & Hattendorf [17]. In the
present study, LA-ICP-MS was firstly applied to the study of agate to answer the questions: (1) is LAICP-MS applicable for very low trace-element contents expected in quartz (agate), and (2) is it possible
to detect trace-element variations with high spatial resolution along a profile perpendicular to the agate
banding.
1. Materials and Methods
The agate material investigated in the present study originates from occurrences in the vicinity of
Nowy Kościόł (Lower Silesia, Poland). The agate outcrops occur at the bottom part of Permian rhyolites
and rhyodacites close to the contact with Middle Rotliegend tuffs and Upper Permian claystones [18].
Porphyric nodules are mostly spherical shaped and can reach sizes of 2–70 cm in diameter. The agates
from Nowy Kościoł show both chalcedony bandings and macrocrystalline quartz zones. They are
typically multicoloured fortification agates with shades of white, reddish, green, brown and yellow
(Fig. 1). Several accessoric phases have been detected in the agates including hematite, sulphide
minerals and monazite [18,19]. A specific feature is the occurrence of organic matter (aliphatic and
subordinately aromatic compounds), which has been related to the postgenetic alteration of the volcanic
rocks [19].
The phase composition of the volcanic host rocks was determined by X-ray diffraction (XRD)
studies. Therefore, material from the surrounding adjacent host rock of the agates was separated and
#4#
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Robert Möckel, Jens Götze…Trace-Element Analysis by Laser Ablation Inductively Coupled Plasma Mass Spectrometry…
Fig. 1. Multi-coloured volcanic agates from Novy Kościόł (Poland) with typical fortification banding
grinded to <20 μm grain size. An URD 6 XRD device (Seifert/Freiberger Präzisionsmechanik) with
CoKα radiation was used and the irradiated area was kept constant at 15 mm². Samples were scanned
with 2θ-steps of 0.03° in the range from 5–80° at step times of 5 seconds per step. Quantification of
powder diffraction patterns was carried out using the Rietveld algorithm BGMN [20].
Polished thin sections of agates and surrounding host rocks were prepared for microscopic
studies (ZEISS Axio Imager A1m) to get information about internal agate structures. In addition, CL
measurements on carbon coated thin sections were performed using a “hot cathode” CL microscope
HC1-LM [21]. The system was operated at 14 kV accelerating voltage and a current density of about
10 μA/mm 2. Luminescence images were captured “on-line” during CL operations using a peltier
cooled digital video-camera (KAPPA 961-1138 CF 20 DXC). CL spectra in the wavelength range
380 to 900 nm were recorded with an Acton Research SP-2356 digital triple-grating spectrograph
with a Princeton Spec-10 CCD detector that was attached to the CL microscope by a silica-glass
fibre guide.
Oxygen isotopic studies on macrocrystalline quartz and chalcedony were performed by a
microanalytical laser-based method at the University of Lausanne, Switzerland. Small samples free of
visible microinclusions were placed on a nickel sample holder and were heated with a 20 W CO2 laser,
while fluorinating due to BrF5 atmosphere. The released oxygen was converted to CO2 by combustion
with a hot carbon rod after passing over a cold trap and through a hot mercury fluorine-getter and then
#5#
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Robert Möckel, Jens Götze…Trace-Element Analysis by Laser Ablation Inductively Coupled Plasma Mass Spectrometry…
Fig. 2. Selected agate from Novy Kościόł (Poland) for detailed LA-ICP-MS investigations; the rectangle represents the prepared area for the trace-element profile
Fig. 3. Schematic sketch of the LA-ICP-MS device
admitted on-line to the mass spectrometer. The precision and accuracy of the method is ± 0.1 ‰ for <
100 μg samples [22].
For trace-element analysis by LA-ICP-MS, a section perpendicular to the agate banding was cut,
polished and fixed with epoxy resin in a 1 inch diameter ring, which fits in the sampling cell of the
LA-ICP-MS device (see marked area in Fig. 2). Measurements were performed at the VSEGEI (St.
Petersburg, Russia) on a system with a Thermo Finnigan Element 2 mass spectrometer, linked to a DUV
193 laser ablation system including optics and viewing by New Wave Research/Merchantek Products.
A principle sketch of a LA-ICP-MS system is given in Figure 3. The ArF excimer laser operates with
a wavelength of 193 nm and an energy output of 50–200 mJ. This constellation enables to produce and
analyse small crater diameters down to 5 μm. The two accredited standards NIST 611 and NIST 612
were used for calibration procedures [23]. Calibration procedures were made at the beginning and at
the end of the measuring cycle and after analysing five spots, mainly in order to minimize memory
effects. The parameters of the system are given in Table 1.
Trace elements were classified concerning possible overlapping of their mass/charge ratio
and first measured with low resolution (LR, R=300) and medium resolution mode (MR, R=4000),
respectively. Table 2 shows the modes for the measured isotopes. Isotopes were selected considering
natural abundances. However, the higher resolution mode is opposite to intensity and therefore, to the
#6#
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Robert Möckel, Jens Götze…Trace-Element Analysis by Laser Ablation Inductively Coupled Plasma Mass Spectrometry…
Table 1. Technical parameters of the LA ICP-MS system
Excimer Laser Compex
Energy output
50–200 mJ
pulse duration/repetition frequency
15 ns/ 1–20 Hz
beam diameter
5–300 μm
carrier gas flow (ablation cell)
approx. 1 l/min (optimized automatically)
Thermo Finnigan Element 2
type
sector field mass spectrometer
auxiliary gas flow
approx. 1 l/min (optimized automatically)
cooling gas flow
approx. 15 l/min (optimized automatically)
power
approx. 1 kW
Table 2. Selected isotopes for LA-ICP-MS measurements (LR low resolution mode, MR medium resolution
mode). Distributions of trace elements printed in bold letters are displayed in Figs. 7 and 8.
LR
Li, 27Al, 28Si, 29Si, 44Ca, 60Ni, 69Ga, 71Ga, 72Ge, 74Ge, 85Rb, 88Sr, 89Y, 133Cs, 138Ba, 139La, 140Ce, 141Pr,
Nd, 146Nd, 147Sm, 151Eu, 153Eu, 155Gd, 157Gd, 159Tb, 161Dy, 163Dy, 165Ho, 166Er, 167Er, 169Tm, 172Yb,
173
Yb, 175Lu, 232Th, 238U
Na12 , 29Si, 30Si, 47Ti, 55Mn, 56Fe, 57Fe, 59Co, 60Ni, 62Ni, 63Cu, 65Cu
7
142
MR
detection limit, which is around 0.01 ppm and 1 ppm for LR and MR, respectively. The given values are
approximate values, since the detection limit is typically three times the standard deviation of the noise
level (blank). Strictly speaking, the detection limit has to be calculated for each element and analyzing
spot, resulting in varying detection limits depending on the specific trace element concentration.
Because of this limitation, the approximate values as stated above were set after all measurements
and proved to be adequate. The same strategy was applied for the error of measurements, which is
estimated to be below 5–10% for all values given below. In order to lower the analytical error, the
measurements of all elements was repeated 20 times on each spot.
As can be seen in Figure 4, a typical signal of LA-ICP-MS measurement for a single element
during ablation of material consists of three phases. In the first phase, the blank from the carrier gas
is recorded, where no laser action takes place and no sample material is introduced into the ICP-MS.
The initiation of the laser is marked by rapid increase of the signal, followed by a stable plateau phase
during which sample material is constantly carried to the ICP-MS (phase two). The proper value of
intensity is averaged from this stable plateau. Of course the length of this plateau depends on the
volume of ablated material and therefore, on laser parameters (laser pulse, repetition rate etc.) and not
least on the material. Phase three is characterized by a steady decrease of the intensity down to the
blank level, after ablation ended.
Measuring points were chosen to be located in both agate and macrocrystaline quartz areas, with
an average distance of 870μm. Since the amount of ablated material from one “shot” is insufficient for
operating 20 analyses, there had to be done several “laser shots” arranged in a row in order to avoid
fractionating processes in the crater [17,24,25].
#7#
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Robert Möckel, Jens Götze…Trace-Element Analysis by Laser Ablation Inductively Coupled Plasma Mass Spectrometry…
Fig. 4. A simplified diagram of a typical LA-ICP-MS signal for one element displaying three different phases
2. Results and discussion
Agate formation and structure
The XRD analyses provided first results concerning the mineral composition of the adjacent host
rocks of the agates (Table 3). The fresh material is dominated by quartz and K-feldspar and represents
a rhyolitic composition. Alteration processes resulted in a loss of amorphous material (glassy matrix)
and feldspar, respectively and the formation of clay minerals (especially illite). The illitization of the
volcanic host rocks can be one of the processes providing silica for the agate formation.
For the estimation of formation temperatures, δ18O values were measured for the silica matrix.
The δ18OSMOW data of the analyzed agate range from 29.36 ‰ (quartz) to 30.25 ‰ (chalcedony)
(Table 4). Because of the unknown origin of the participating fluids, the temperatures were
calculated for different sources according to Méheut et al. [26] (see Table 4). The calculations
point to relatively low temperatures between 20°C and 120°C. These results are in accordance with
measured homogenization temperatures of primary fluid inclusions in quartz, which scatter between
100 and 215°C [27].
The agates itself are typically banded with alternate bedding of chalcedony bands and
macrocrystalline quartz layers (Fig. 5a). Isolated spherulitic silica growth is common (Fig. 5b), and
areas with macrocrystalline quartz have formed in the central and marginal parts of the agates (Fig. 5c).
In some cases, the transition from quartz to chalcedony is sharp. Iron oxides are often interbedded
along the crystallization front of chalcedony, probably resulting from a “self-purification” process of
originally iron-rich silica.
CL imaging revealed several differently coloured zones within the agates. The chalcedony layers
display reddish, rose and orange colours (Fig. 5). Phanerocrystalline quartz may exhibit similar CL
colours as the surrounding chalcedony, but sometimes more bluish CL colours and distinct complex
zoning (both sector and concentric zoning) (Fig. 5). The comparison of the differently luminescent
agate bands and macrocrystalline quartz (Fig. 6) illustrates that the CL spectra generally consist of
the same emission bands. The different CL colours and intensities, respectively, are a result of varying
intensity ratios of the main CL emission bands. The most common and most intense CL emission
#8#
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Robert Möckel, Jens Götze…Trace-Element Analysis by Laser Ablation Inductively Coupled Plasma Mass Spectrometry…
Table 3. Phase composition of adjacent host rocks of the investigated agates determines by XRD analysis (contents in weight-%)
Fresh host rock
(greyish brown)
35 ± 2
Altered host rock
(greenish)
35 ± 1
Monoclinic K-feldspar
40 ± 2
38 ± 2
Plagioclase
6±1
5±1
Quartz
Kaolinite
2±1
4±1
Illite
2±1
9±1
Amorphous
15 ± 4
10 ± 3
Table 4. Oxygen isotope data of agate and quartz with calculated temperatures according to [26]; the temperatures
are calculated for fractionation with meteoric water (-10‰), oceanic water (±0‰) and magmatic water (+8‰)
Sample
δ18OSMOW(‰)
Tmet (°C)
Toc (°C)
Tmag (°C)
quartz
29.36
31
70
115
agate
30.25
28
66
109
band occurs around 650 nm (Fig. 6) and is related to the so called non-bridging oxygen hole centre
(NBOHC [28]). This defect type may result from several precursors such as silanol groups, peroxy
linkage (oxygen-rich samples) or strained Si-O bonds [29]. In the case of agates we have high contents
of silanol groups (Si-OH), which are common in microcrystalline quartz [9,11,30]. The lower intensity
of the 650 nm emission in macrocrystalline quartz is a result of the lower density of this defect type.
Other emission bands (450 nm, 580 nm) appear only subordinately, partly as shoulder of the main
CL emission (Fig. 6). The yellow CL emission (580 nm), which can be related to the E’ centre (electron
defect on an oxygen vacancy), was up to now only observed in agates, silicified wood and hydrothermal
quartz [9,11,31]. The emission band around 450 nm, which is especially present in macrocrystalline
quartz, is due to the recombination of the so-called self-trapped exciton [29], an electron hole pair on
an oxygen vacancy and a peroxy linkage. The oxygen vacancies may indicate crystallization under
conditions with oxygen deficit.
The results of the microscopic studies indicate that the agate crystallized in a continuous process
from a silica-rich precursor. Crystallization started with spherulitic growth and the formation of
chalcedonic quartz, later changing into well developed crystals. Accordingly, the phanerocrystalline
quartz seems to be the prosecution of chalcedony “fibres”. The high defect density and the occurrence
of sector zoning in quartz may result from rapid crystallization under non-equilibrium conditions.
Except for the outermost layer of chalcedony, CL colours and spectra are almost the same in both quartz
and chalcedony and display a similar local defect structure, which implies the origin from a single
ongoing crystallization process. The sometimes observed repeated change between chalcedony and
macrocrystalline quartz layers may result from changes in the silica concentration of the mineralizing
fluids or from several mineralization stages.
#9#
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Robert Möckel, Jens Götze…Trace-Element Analysis by Laser Ablation Inductively Coupled Plasma Mass Spectrometry…
Fig. 5. Polarized light (Pol) and cathodoluminescence (CL) micrograph pairs of the investigated agates from Novy
Kościόł (Poland); (a) alternate bedding of chalcedony bands and macrocrystalline quartz layers showing varying
CL colours; (b) spherulitic chalcedony growth with distinct internal structure; (c) strongly zoned macrocrystalline quartz
Trace element studies by LA-ICP-MS
The results of the present study show that most of the measured trace-elements are present in
very low concentrations. Because of the specific analytical conditions, the limits of detection vary
for different elements (Table 5). Therefore, the strategy of the analytical procedure was to lower the
limits of detection by repeated “laser shots” (more sample material) and valuable standard materials.
Nevertheless, concentrations of some elements were below the detection limits. For example, the
contents of the rare earth elements (REE) did not provide important data in terms of REE distribution
patterns.
Problems during measurements appeared with Ca and Ni. Ni is part of the cone material of the
skimmer technology of the mass spectrometer. Ni was measured in LR as well as MR mode (see
# 10 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Fig. 6. Typical CL spectra from chalcedony and macrocrystalline quartz; the CL spectrum of quartz shows two
broad emission bands at 450-500 nm and ca. 650 nm; the CL of chalcedony is dominated by a strong emission
band at 650 nm and shows at least three additional bands (shoulders) at 580, 500 and 450 nm
Table 5. Average trace-element contents and variation limits (in ppm) of the investigated profile perpendiculat to
the agate banding (see Fig. 2) (quartz: measunring points 1–7; agate: measuring points 8–12)
quartz
average
agate
max
min
average
max
min
Li
1.01
2.77
0.03
2.73
4.69
0.84
Na
24.0
59.2
1.35
192.3
294.1
106.7
Al
13.7
47.2
2.71
325.4
640.8
58.0
Ge
0.63
1.06
0.32
12.12
20.42
7.42
Rb
0.12
0.17
0.07
1.60
3.05
0.66
Sr
0.07
0.09
0.03
3.29
7.58
0.48
Ba
1.40
0.37
0.03
5.00
6.38
0.19
Y
0.01
0.02
0.01
0.20
0.03
–
La
0.09
0.09
0.08
–
0.62
–
Ce
0.10
0.12
0.08
0.04
0.12
0.01
Nd
0.03
0.03
0.03
0.05
Sm
0.02
0.01
–
–
0.05
–
0.05
–
0.01
Eu
–
0.02
–
0.01
0.01
Gd
–
–
–
0.03
0.04
0.01
Dy
–
–
–
0.04
0.08
0.01
Ho
–
–
–
0.03
0.03
0.03
Er
–
–
–
0.03
0.04
0.01
Tm
–
–
–
0.01
0.01
0.01
Yb
–
–
–
0.06
0.06
0.06
Th
0.37
0.37
0.37
0.02
0.03
0.01
U
0.15
0.02
0.52
2.39
3.37
1.81
Ti
<1
<1
<1
3.40
<1
<1
<1
<1
1.20
Mn
0.61
1.39
<1
Fe
0.95
2.62
<1
90.4
119.8
69.3
(-) below detection limit
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Robert Möckel, Jens Götze…Trace-Element Analysis by Laser Ablation Inductively Coupled Plasma Mass Spectrometry…
Fig. 7. Distribution profiles of selected trace elements (Ba, Sr, Rb, Cs, Cu, U) in the investigated agate from Novy
Kościόł (Poland)
Table 2), mostly in order to avoid overlapping with other isotopes (e.g. 44Ca16O, 24Mg36Ar). However,
measured values differed sometimes more than 50% and therefore, have been not reliable. The 44Ca
values were mostly unreliable as well due to overlapping with 28Si16O, which is basic matrix material.
A higher resolution mode for Ca or the use of another inert gas instead of argon (40Ar overlaps
40
Ca), respectively, may prevent this problem. In the present study, the problem was tried to solve
by an afterwards calculation of a 28Si16O-blank. Although the measured Ca contents are similar to
concentrations (330–634 ppm) reported by [19], there was still some doubt on the reliability of these
corrected values and therefore, they will not be presented here.
A number of trace-elements could be analyzed with high precision. The results of the traceelement analyses are given in Table 5 (average values and limits of variation). The spatial distribution
of selected element contents is presented in Figures 7 and 8.
The comparison of the data from Nowy Kościόł with trace element contents of agates from
different regions [3,10] revealed that concentrations presented here are slightly below other published
values, but more or less in the same order of magnitude. In general, it is obvious that a number of
# 12 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Robert Möckel, Jens Götze…Trace-Element Analysis by Laser Ablation Inductively Coupled Plasma Mass Spectrometry…
Fig. 8. Distribution profile of selected trace elements (Al, Li, Na, Ge, Ga, Fe) in the investigated agate from Novy
Kościόł (Poland)
elements is accumulated in agate. Besides the common elements Al, Fe, and Na, also Ba, Sr, U and
Ge show relatively high concentrations. These elements may be mobilized together with Si during the
alteration of the volcanic host rocks. The decomposition of volcanic glass and feldspar provides high
amounts of Si, Al, Ca, Na and K. The mobility of U during the alteration of volcanic rocks was reported
by Zielinski [32] who observed the parallel accumulation of Si and U. In the case of Ge, the identical
geochemical character with Si promotes the common transport and the incorporation of Ge into the
quartz structure [33,34]. Therefore, high concentrations of Ge are typical for agates from acidic host
rocks [3].
One characteristic feature of the trace-element distributions is that almost all elements show higher
contents in chalcedony compared to macrocrystalline quartz (Fig. 7 and 8). Aluminium, for instance,
has an average concentration of 325.4 ppm in chalcedony and 13.7 ppm in quartz. Fe shows a hundred
times higher concentration in chalcedony parts (average 90.4 ppm) than in the quartz centre (average
0.95 ppm). The only exception for this is Cu, which is slightly enriched in quartz (average of 0.07 ppm)
compared to chalcedony (average 0.04 ppm). Moreover, elements which behave geochemically similar
# 13 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Robert Möckel, Jens Götze…Trace-Element Analysis by Laser Ablation Inductively Coupled Plasma Mass Spectrometry…
such as Ba and Sr or Rb and Cs display comparable patterns (Figure 7). These similar patterns suggest
a simultaneous mobilization and precipitation of the elements.
The differences in most element concentrations between chalcedony and quartz can probably be
explained by a “self-purification” process during the crystallization of the agate. The trace elements
are firstly precipitated together with the primary silica gel and are finely dispersed in the silica matrix.
Starting spherulitic growth of the chalcedony may initially incorporate high contents of trace elements.
During further growth, the concentration of impurities drops, which is macroscopically detectable in
the changing visible colour of the agate bands. In general, deeply coloured parts of the agate inherit
higher trace-element concentrations than colourless parts. For example, the outermost measuring
points are situated in a dark reddish band and reveal in almost all cases higher element contents,
whereas the colour gets more and more light towards the inner core.
Heaney [35] stated that the high trace-element contents in silica can cause the formation of frequent
defects, Brazil twinning and typical twisting in chalcedony. He concluded that chalcedony fibres in
agates may be the result of dislocation growth. This statement would explain why macrocrystalline
quartz can only grow after removing most of the impurities. The transition zone between chalcedony
and quartz is characterized by a characteristic drop of most trace-element contents. On the other hand,
the higher defect density in chalcedony was also detected by CL spectroscopy.
The kind of incorporation of various trace elements into quartz can be different [34]. Whereas
Al, Ti, Ge, Fe, Na or Li are often substitute for Si or placed at interstitial lattice positions, elements
such as Ba, Sr, Rb, Cs or U may be concentrated in fluid and mineral inclusions or at grain
boundaries. However, the incorporation into the quartz structure is limited as is visible for iron.
The occurrence of high iron contents and sometimes iron oxides within the agate matrix is an
indication for high iron contents of the silica-bearing mineralizing medium. The accumulation
of hematite spheres along the banding of the next chalcedony layers may be a result of the “selfpurification” of chalcedony during crystallization. The iron was removed from the silica during
crystallization, moved with the crystallization front and accumulated forming iron oxide spheres.
Especially the coloured agate bands exhibit high concentrations of iron emphasizing that iron
is mostly concentrated in fi ne iron oxide particles within the chalcedony layers acting as colour
pigments.
3. Conclusions
The results of the present study clearly demonstrate that LA-ICP-MS is a suitable analytical
method for the detection of trace elements in very low concentrations. Although the element contents
of most of the elements in quartz are in the ppm and sub-ppm level, the method enables to realize a
spatially resolved analysis. Therefore, the method provides valuable data for genetic interpretations in
samples with marked micro- and/or macro-textures. Only for a few elements (especially the REE), the
concentrations are often too low to be detected by LA-ICP-MS.
The trace-element analysis with high spatial resolution enabled to monitor the distribution of
specific elements during proceeding agate formation. One characteristic feature is the release of
elements (e.g., Al, Fe, Ge) together with Si from the volcanic host rocks due to the ongoing alteration,
and the subsequent accumulation within the silica matrix during agate formation. Elements showing
a similar geochemical behaviour (e.g., Ba-Sr or Rb-Cs) display identical spatial distribution pattern
# 14 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Robert Möckel, Jens Götze…Trace-Element Analysis by Laser Ablation Inductively Coupled Plasma Mass Spectrometry…
within the agate texture. Most of the impurities precipitate at the beginning of the agate formation,
whereas further growth and crystallization of the silica matrix results in lower impurity concentrations.
Within the last-stage macrocrystalline quartz in the agate centre, trace-element contents are very low
or even not detectable. These facts point to a “self-purification” process of the silica matrix during
ongoing crystallization, sometimes resulting in the admixture of impurity phases (e.g. iron oxides/hydroxides).
Acknowledgement
We gratefully acknowledge providing of agate samples from Novy Kościόł by Marek Stepisiewicz
(Warsaw). We also thank Dmitry Sergeev for helping a lot to overcome language barriers at the VSEGEI
St. Petersburg.
References
[1] Zenz, J. (2005) Agate. Bode Verlag, Haltern.
[2] Godovikov, A.A., Ripinen, I.I., Motorin, S.G. (1987) Agates. Izd. Nedra, Moscow (in
Russian).
[3] Blankenburg, H.-J. (1988) Agate. VEB deutscher Verlag für Grundstoffindustrie, Leipzig (in
German)
[4] Blankenburg, H.-J., Pilot, J., Werner, C.D. (1982) Erste Ergebnisse der Sauerstoffisotopenuntersuchungen an Vulkanitachaten und ihre genetische Interpretation. Chemie Erde 41:
213–217.
[5] Fallick, A. E., Jocelyn, J., Hamilton, P. J. (1987) Oxygen and Hydrogen stable isotope
systematics in brazilian agates. In: Rodriguez-Clemente, R. and Tardy, Y. (Eds.), Geochemistry and
mineral formation in the earth surface, 99–117.
[6] Harris, C. (1988) Oxygen isotope Geochemistry of a Quartz-Agate Geode from northwestern
Namibia. Communs geol. Surv. S.W. Africa Namibia 4: 43–44.
[7] Blankenburg, H.-J., Thomas, R., Klemm, W., Leeder, O. (1990) Interpretation der Ergebnisse von
Einschlußuntersuchungen an den Quarzinkrustaten aus Vulkanitachaten. Zeitschrift für geologische
Wissenschaften 18: 81–85.
[8] Strauch, G., Nitzsche, H.-M., Holzhey, G. (1994) Isotopenuntersuchungen an Rhyolithen und
Achatbildungen. N. Jb. Mineral., Abh. 165: 103–104.
[9] Götze, J., Plötze, M., Fuchs, H., Habermann, D. (1999) Defect structure and luminescence
behaviour of agate - results of electron paramagnetic resonance (EPR) and cathodoluminescence (CL)
studies. Mineralogical Magazine 63: 149–163.
[10] Götze, J., Tichomirowa, M., Fuchs, H., Pilot, J., Sharp, Z.D. (2001) Geochemistry of agates:
a trace element and stable isotope study. Chemical Geology 175: 523–541.
[11] Moxon, T. and Reed, S.J.B. (2006) Agate and chalcedony from igneous and sedimentary hosts
aged from 13 to 3480 Ma: a cathodoluminescence study. Mineralogical Magazine 70: 485–498
[12] Heaney, P. J. & Davis, A. M. (1995) Observation and origin of self-organized textures in
agates. Science 269: 1562–1565.
[13] Gray, A.L. (1985) Solid sample introduction by laser ablation for inductively coupled plasma
source mass spectroscopy. Analyst 110: 551–556
# 15 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Robert Möckel, Jens Götze…Trace-Element Analysis by Laser Ablation Inductively Coupled Plasma Mass Spectrometry…
[14] Günther, D., v. Quadt, A., Wirz, R., Cousin, H., Dietrich, V.J. (2001) Elemental analyses using
laser ablation-inductively coupled plasma-mass spectrometry (LA-ICP-MS) of geological samples
fused with Li2B4O7 and calibrated without matrix-matched standards. Mikrochimica Acta 136: 100–
107
[15] Flem, B., Larsen, R.B., Gromstvedt, A., Mansfeld, J. (2002) In situ analysis of trace elements
in quartz by using laser ablation inductively coupled plasma mass spectrometry. Chemical Geology
182: 237–247.
[16] Heinrich, C.A., Pettke, T. Halter, W.E. Aigner-Torres, M., Audétat, A., Günther, D., Hattendorf,
B., Bleiner, D., Guillong, M., Horn, I. (2003) Quantitative multi-element analysis of minerals, fluid
and melt inclusions by laser-ablation inductively-coupled-plasma mass-spectrometry. Geochimica et
Cosmochimica Acta 67(18): 3473–3496
[17] Günther, D. and Hattendorf (2005) Solid sample analysis using laser ablation inductively
coupled plasma mass spectrometry. Trends in Analytical Chemistry 24(3): 255–265
[18] Gaweda, A. and Rzymełka, J.A. (1992) Bituminous agates from rhyolites in the environs of
Nowy Kościół, Lower Silesia. Mineralogia Polonic 23(1): 72–84
[19] Dumańska-Słowik, M., Natkaniec-Nowak, L., Kotarba, M.J., Sikorska, M., Rzymetka, J.A.,
Loboda, A., Gaweł, A. (2008) Mineralogical and geochemical characterization of the „bituminous“
agates from Nowy Kościół (Lower Silesia, Poland). Neues Jahrbuch für Mineralogie, Abhandlungen
184(3): 255–268.
[20] Taut, T., Kleeberg, R., Bergmann, J. (1998) Seifert Software: The new Seifert Rietveld
program BGMN and its application to quantitative phase analysis. Materials Structure 5(1): 57–66.
[21] Neuser, R.D., Bruhn, F., Götze, J., Habermann, D., Richter, D.K. (1995) Kathodolumineszenz:
Methodik und Anwendung. Zentralblatt für Geologie und Paläontologie, Teil I, H. 1/2: 287–306.
[22] Sharp, Z.D. (1992) In situ laser microprobe techniques for stable isotope analysis. Chemical
Geology 101: 3–19.
[23] Hollocher and Ruiz (1995) Major and trace element determinations on NIST glass standard
reference materials 611, 612, 614 and 1834 by inductively coupled plasma-mass spectrometry.
Geostandards Newsletter 19(1): 27–34
[24] Fryer, B.J., Jackson, S.E., Longerich, H.P. (1995) The design, operation and role of the laserablation microprobe coupled with an inductively coupled plasma-mass spectrometer (LAM-ICP-MS)
in the earth sciences. The Canadian Mineralogist 33: 303–312
[25] Fernández, B., Claverie, F., Pécheyran, C., Donard, O. (2007) Direct analysis of solid
samples by fs-LA-ICP-MS. Trends in Analytical Chemistry 26(10): 951–966.
[26] Méheut, M., Lazzeri, M., Balan, E., Mauri, F. (2007) Equilibrium isotopic fractionation in the
kaolinite, quartz, water system: Prediction from first-principles density-functional theory. Geochimica
et Cosmochimica Acta 71: 3170–3181.
[27] Kozłowski, A. (1985) Studies of fluid inclusions in agates from Plόczki Górne and Nowy
Kościόł, Poland. Fluid Inclusion Research – Proceedings of COFFI, Michigan University Press, Ann
Arbor, USA, 18: 221–222.
[28] Siegel, G. H. and Marrone, M. J. (1981) Photoluminescence in as-drawn and irradiated silica
optical fibers: An assessment of the role of nonbridging oxygen defect centres. Journal of Noncrystalline
Solids 45: 235–247.
# 16 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Robert Möckel, Jens Götze…Trace-Element Analysis by Laser Ablation Inductively Coupled Plasma Mass Spectrometry…
[29] Stevens Kalceff, M.A., Phillips, M.R., Moon, A.R., Kalceff, W. (2000) Cathodoluminescence
microcharacterisation of silicon dioxide polymorphs. In: Pagel, M., Barbin, V., Blanc, P., Ohnenstetter,
D. (eds.) (2000) Cathodoluminescence in geosciences. Springer Verlag, Berlin Heidelberg New York,
193–224.
[30] Flörke, O.W., Köhler-Herbertz, B., Langer, K., Tönges, I. (1982) Water in microcrystalline
quartz of volcanic origin: agates. Contributions to Mineralogy and Petrology, 80, 324–333.
[31] Götze, J., Plötze, M., Habermann, D. (2001) Origin, spectral characteristics and practical
applications of the cathodoluminescence (CL) of quartz – a review. Mineralogy and Petrology 71:
225–250.
[32] Zielinski, R.A. (1979) Uranium mobility during interaction of rhyolitic obsidian, perlite and
felsite with alkaline carbonate solution: T=120°C, P=210 kg/cm2. Chemical Geology 27: 47–63.
[33] Walenzak, Z. (1969) Geochemistry of minor elements dispersed in quartz (Ge, Al, Ga, Fe, Ti,
Li and Be). Archiwum Mineralogiezne 28: 189–335.
[34] Götze, J., Plötze, M., Graupner, T., Hallbauer, D.K., Bray, C. (2004) Trace element incorporation
into quartz: a combined study by ICP-MS, electron spin resonance, cathodoluminescence, capillary
ion analysis and gas chromatography. Geochimica et Cosmochimica Acta 68: 3741–3759.
[35] Heaney, P.J. (1993) A proposed mechanism for the growth of chalcedony. Contributions to
Mineralogy and Petrology 115: 66 – 74.
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Robert Möckel, Jens Götze…Trace-Element Analysis by Laser Ablation Inductively Coupled Plasma Mass Spectrometry…
Трассировочный анализ кривых
с помощью метода лазерной абляции
с индуктивно связанной плазмой
масс-спектрометрией (LA-ICP-MS):
на примере агатов
из Нового Костела, Польша
Р. Мёкела, Й. Гетца, С. А. Сергеевб,
И. Н. Капитоновб, Е. В. Адамскаяб,
Н. А. Гольдзинб, Т. Веннеманв
а
Технический университет
Фрайбергская горная академия, Институт минералогии,
Германия, Фрайберг, D-09596, Brennhausgasse 14
б
Всероссийский научно-исследовательский
геологический институт им. А. П. Карпинского
Россия 199106, Санкт-Петербург, Средний проспект, 74
в
Институт минералогии и геохимии,
Университет Лозанны, UNIL-BFSH2,
Швейцария, Лозанна, СН-1015
Метод лазерной абляции с индуктивно связанной плазмой масс-спектрометрией (LA-ICP-MS)
был применен для обнаружения следов элементов в агатах из пермских вулканитов (Новый
Костел, Польша), имеющих низкие концентрации и высокие пространственные разрешения.
Используемая система LA-ICP-MS состоит из DUV 193 системы лазерной абляции, связанной с
термоэлементом Finnigan 2 масс-спектрометра. Применение 193 нм ArF эксимерного лазера (с
выходом 50-200 мДж) и стандарта NIST 611 и NIST 612 позволили получить и проанализировать
кратеры малого диаметра до 5 мкм.
Были проанализированы кривые следующих элементов Ti, Ge, Al, Fe, Mn, U, Th, Ba, Sr, Rb, Cs
и Y на уровнях < млн-1 и ниже. Концентрации редкоземельных элементов иногда были ниже
предела обнаружения метода. Почти все элементы (кроме меди) показали более высокие
концентрации в халцедоне, чем в крупнокристаллическом кварце, железо, например, имеет
в 100 раз более высокую концентрацию полос агата, чем кварц, это может быть связано с
тонким распределением частиц оксида железа в халцедоне, которые могут действовать в
качестве цветных пигментов.
Кривые элементов в агате освобождаются одновременно с кремнием при изменении
окружающих вулканических пород. Данные изотопов кислорода показывают, что накопление
кремнезема и формирование агата происходит при температурах ниже 120 °С. Характер
схемы распределения элементов в агате объясняется процессом «самоочищения» в ходе
кристаллизации халцедона и кварца из силикагеля.
Ключевые слова: агат, кварц, микроэлементы, LA-ICP-MS, изотопы кислорода, CL, Новый
Костел.
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Journal of Siberian Federal University. Engineering & Technologies 1 (2012 5) 19-27
~~~
УДК 671.03
Study of Effect of Nickel and Rare-Earth Metals
on the Structure and Properties
of Cast and Deformed Semi-Finished Products
Made of Aluminum Alloys
Sergey B. Sidelnikovа*, Vladimir N. Baranovа,
Ekaterina S. Lopatinaа and Leonid P. Trifonenkovb
а
Siberian Federal University,
79 Svobodny, Krasnoyarsk, 660041 Russia
b
LLC “RUSAL ETC”,
37/1 Pogranichnikov, Krasnoyarsk, 660111 Russia 1
Received 6.02.2012, received in revised form 13.02.2012, accepted 20.02.2012
Technology for making electrical aluminum alloys with different contents of rare-earth metals and
nickel was developed. Cast and deformed semi-finished products made of new alloys were obtained
as a result of experimental studies. The structure and mechanical properties of metal were studied
as well as electrical resistivity measurements of rods and wires were carried out. The effect of nickel
and rare-earth metals on these characteristics was found out and alloy compositions with the most
favorable combination of strength and electrical properties were determined.
Keywords: aluminum, nickel, rare earth metals, alloy, mechanical properties, electro conductivity,
structure, casting parameters.
The results of this study were obtained during realization of the project which is funded jointly by
the Ministry of Education and Science of Russian Federation and UC RUSAL and aims to organize an
industrial production of low-alloyed aluminum alloy rod for the subsequent manufacture of electrical
wires. The wires should have following properties: electrical resistivity similar to electrotechnical
grade aluminum A5E (0,0295 ohm*m/mm2) and tensile strength up to 250 Mpa with the increased
value of specific strength up to 92.6 N*m/g while it does not exceed 60 N*m/g for the best grades of
hard-drawn aluminum wires made of aluminum A5E. Specific strength can reach values up to 110
N*m/g for the electrical wires made of aluminum alloys with different contents of nickel and rareearth metals (REM) but at the same time conductivity is reduced to 86% of electrical conductivity of
aluminum.
The most productive way to obtain wire rod being used for electrical wires is combined
casting and rolling process in the casting-rolling units (CRU). However it is designed to produce
the wire rod made of technical (A5, A7) and low-alloyed aluminum so production of high-strength
*
1
Corresponding author E-mail address: sbs270359@yandex.ru
© Siberian Federal University. All rights reserved
# 19 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Sergey B. Sidelnikov,Vladimir N. Baranov… Study of Effect of Nickel and Rare-Earth Metals on the Structure...
a
b
Fig. 1. Diagrams of phase equilibrium in aluminum – nickel (a) and nickel – cerium (b) systems
alloys can lead to equipment failure and defective products. Therefore it is reasonable to use the
continuous technologies of combined casting, rolling and pressing (CCRP) for these alloys [1].
These technologies are characterized by low energy consumption and high productivity and allow
the fi nal product to be obtained in one or two technological steps, which is important in modern
manufacturing.
Thus the combination of the highest specific strength and electrical conductivity (minimal
electrical resistance equates to aluminum wire) of electrical wire rod made of low-alloyed aluminum
alloys with optimum complex alloying can be achieved through the use the author’s CCRP method
for producing of long-length semi-finished products and subsequent drawing with a given degree of
deformation and annealing.
The efficiency of hardening effects of the alloying at the expense of complicating the structure
of the hardening phases such as through additional alloying by transition metals was evaluated in this
study. These transition metals are slightly soluble in aluminum like REM otherwise they can form
intermetallic phases with aluminum and REM as well as with iron and silicon, which are the main
impurities in aluminum. This further reduces the content of such impurities in the solid solution and
leads to an increase in hardening while maintaining the electrical resistivity. Typical representative
of such transition metals is nickel. Diagram of phase equilibrium with aluminum (Fig. 1 a) shows the
possibility of such an influence of nickel on the phase state of electrical alloys.
According to the data provided by the phase equilibrium diagram the maximum solubility of
nickel at the eutectic temperature (640 ºC) is 0,05 % (atomic), but it sharply decreases to 0.01% at
the temperature of 500 ºC (Table 1). Nickel as well as REM forms intermetallic compound Al3Ni
which is in equilibrium with the aluminum solid solution, whose particles are practically insoluble in
aluminum at the temperatures below 500 ºC. The presence of nickel in aluminum alloy in conjunction
with REM can lead to the interaction between nickel and these components, as double diagram
indicates (Fig. 1 b).
According to [2] nickel in the presence of iron (in equal amounts) in aluminum-based alloy
interacts with it forming a ternary compound FeNiAl9. This phase is practically insoluble during the
heating and do not adversely affect the mechanical properties of the alloy.
# 20 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Sergey B. Sidelnikov,Vladimir N. Baranov… Study of Effect of Nickel and Rare-Earth Metals on the Structure...
Table 1. Features of dual diagrams Al-Ce, Al-La, Al-Ni
Metal solubility limit
in aluminum, mass %
(at the eutectic temperature)
0,05
0,05
0,05
Metal
Ce
La
Ni
Eutectic
concentration,
mass %
12,0
12,0
6,0
Eutectic
temperature, °С
Phase
637
642
640
Al4Ce (Al11Ce3)
(Al11La3) Al4La
Al3Ni
Table 2. Casting modes and chemical compositions of Al-REM and Al-REM-Ni alloys
Sample №
1
2
3
4
5
6
7
8
Ni concentration,
mass %
0,15-0,20
0,20-0,25
0,25-0,30
0,30-0,35
Relative content
of REM in alloy
Very low
Low
Medium
High
Very low
Low
Medium
High
Heat temperature,
°С
Pouring
temperature, °С
750-800
740±10
Analysis of phase equilibriums in alloys based on aluminum showed that it is possible to create
semi-finished products of the Al-REM-Ni alloy for electrical purpose with the correct iron and
silicon content in amounts corresponding to the content of these components in the commercial grade
aluminum.
In order to experimentally verify the analytical conclusions the new aluminum alloys with
different contents of rare-earth metals (cerium based) and nickel were developed in the laboratory of
Siberian State University as well as the technological parameters of their preparation and casting were
studied. Al-Ni and Al-REM waffle alloys made by KBM Company were used to prepare electrical
purpose alloys. Primary aluminum melt was prepared in graphite crucibles in a high-frequency melting
furnace and heated to a temperature of 750-800 °C. Then the preheated Al-Ni alloy with calculated
composition and mass (in amount of 0,15-0,35 mass%) was injected directly into the crucible and held
for 3-5 minutes and after that Al-REM alloy was added into the crucible. The melt was poured into
a specially prepared mold and technological sample to determine the chemical composition. Casting
modes and chemical compositions of studied alloys are shown in Table 2.
The macrostructure of the alloys, which was studied using Stemi 2000-C stereoscopic microscope,
mainly consist of columnar crystals zone located on the periphery of the ingot and the zone of large
equiaxed crystals in the center of the ingot. The adverse fan structure is observed in alloys with small,
medium and high REM content and the higher total concentration of alloying elements, the greater
area occupied by fan crystals (Fig. 2).
Research of microstructure of the alloys with different contents of rare-earth metals and Ni was
carried out using light and scanning electron microscope EVO 50 with energy-dispersion microanalyzer
Inca Energy 350. The microstructure of the alloys is shown in Fig. 3. It was found that the dendritic
cells size reduces with increasing concentrations of rare-earth elements and there are Al4Me eutectic
# 21 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Sergey B. Sidelnikov,Vladimir N. Baranov… Study of Effect of Nickel and Rare-Earth Metals on the Structure...
ʋ1
ʋ5
ʋ2
ʋ6
ʋ3
ʋ7
ʋ4
ʋ8
Fig. 2. Macrostructure of Al-REM alloys (samples №1-4) и Al-REM-Ni alloys (samples №5-8)
phases along the boundaries of dendritic cells in the low-alloyed alloys (sample №1). The amount of
eutectic component increases with increasing concentration of REM and eutectic colonies (α+Al4Ме)
form the volume of dendritic cells (samples № 2-4).
Alloying of Al-REM alloys with small additions of Ni from 0.15 to 0.35 % (samples № 5-8) has
no significant effect on the microstructure of the alloys.
Hot-deformed semi-products with circular cross-section in the form of 9 mm rod were produced
from these alloy ingots using combined casting, rolling and pressing process [1]. Then the wire 2 mm
in diameter was produced from this rod by cold-drawing and the samples obtained were researched in
relation to their structure, mechanical properties and electrical resistance.
Typical distribution of constituents throughout the cross-section of the samples occurs when
producing the rods from Al-REM and Al-REM-Ni alloys by CCRP process which leads to both
macro- and microheterogeneity. In low-alloyed alloys (samples № 1, № 5) there are some areas
# 22 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Sergey B. Sidelnikov,Vladimir N. Baranov… Study of Effect of Nickel and Rare-Earth Metals on the Structure...
ʋ1
ʋ5
ʋ2
ʋ6
ʋ3
ʋ7
ʋ4
ʋ8
Fig. 3. Microstructure of Al-REM alloys (samples № 1-4) and Al-REM-Ni alloys (samples № 5-8), ×200
with an increased number of eutectic component in the form of lines oriented in the direction of
deformation.
An increase in concentration of rare-earth metals in the alloys leads to significant heterogeneity
throughout the cross-section of the rod (samples № 2 – 8). Areas containing eutectic colonies of crystals
are located along the boundaries of deformed grains. Eutectic colonies generate the areas in the form of
clusters with a high content of alloying elements. Microheterogeneity of alloying elements throughout
the cross section of the rod is confirmed by the presence of large intermetalloids predominantly of
platelet shape in the eutectic structure.
Electron probe microanalysis of intermetalloids (samples № 3, 4) shows the presence of Al, Ce,
La and Si. The rods high-alloyed with rare-earth metals have Ce content higher than La and it reaches
4.5% in Al-REM alloys and 7.9% in Al-REM-Ni alloys. Microheterogeneity of the rods appears in the
form of areas with a low amount of eutectic constituent (Table 3, Fig. 4).
# 23 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Sergey B. Sidelnikov,Vladimir N. Baranov… Study of Effect of Nickel and Rare-Earth Metals on the Structure...
Table 3. Microhardness of deformed semi-finished products made of Al-REM and Al-REM-Ni alloys
Rod microhardness HV, kgf/mm 2
Wire microhardness HV, kgf/mm 2
№
Al-REM
№
Al-REM-Ni
№
Al-REM
№
Al-REM-Ni
1
2
3
4
38,4±2,7
42,7±1,8
41,7±1,6
45,5±1,0
5
6
7
8
42,2±1,3
40,7±3,4
46,7±2,0
45,7±2,0
1
2
3
4
52,7±2,3
54,0±2,6
54,4±1,0
61,6±2,3
5
6
7
8
57,0±1,1
51,1±1,2
61,5±1,5
67,4±2,1
ʋ1
ʋ5
ʋ2
ʋ6
ʋ3
ʋ7
ʋ4
ʋ8
Fig. 4. The microstructure of the rods made of Al-REM (Samples № 1-4) and Al-REM-Ni (samples № 5-8) alloys
produced by CCRP process, ×200
# 24 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Sergey B. Sidelnikov,Vladimir N. Baranov… Study of Effect of Nickel and Rare-Earth Metals on the Structure...
Wire drawing out of studied rods did not eliminate the macro- and microheterogeneity of the rods.
The structures of the rods and wires with different compositions are similar and differs by presence
more expressed streaked structure for cold deformation. It should be noted that plastic properties of
the developed alloys provide cold deformation without annealing and ruptures with a total degree of
deformation up to 70-80 %.
Studies of tensile strength of the test rods and wires were carried out using the universal
electromechanical machine LFM 400 (400 kN) which can record the main parameters of the process
on a computer. The research results are shown on Fig. 5-6. Results of specific electrical resistivity
measurements of obtained rods and wires are presented in Fig. 7-8.
Fig. 5. Tensile strength of the rods made of Al-REM and Al-REM-Ni alloys produced by CCRP process
Fig. 6. Tensile strength of the wires made of Al-REM and Al-REM-Ni alloys
# 25 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Sergey B. Sidelnikov,Vladimir N. Baranov… Study of Effect of Nickel and Rare-Earth Metals on the Structure...
Fig. 7. Specific electrical resistivity of the rods made of Al-REM and Al-REM-Ni alloys produced by CCRP
process
Fig. 8. Specific electrical resistivity of the wires made of Al-REM and Al-REM-Ni alloys
The results of studies of the mechanical properties of rods and wires of the alloys shows that
their strength properties reach values up to 150 MPa for rods and up to 200 MPa for the wire, what is
much higher than tensile strength according to the standard requirements for deformed semi-finished
products made of electrotechnical grade aluminum. But on the other hand resistivity has the same level
in case if the content of rare-earth metals and nickel in the alloy decreases.
The made studies shown that in principle it is possible to obtain a new aluminum alloys with
additions of rare-earth metals and nickel with the required level of mechanical and electrical properties
and fine structure. Developed technical solutions in addition to an increase in production efficiency of
long-length semi-finished products of aluminum alloys [3] enable to get most favorable combination of
specific strength and electrical conductivity due to the optimal complex alloying.
# 26 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Sergey B. Sidelnikov,Vladimir N. Baranov… Study of Effect of Nickel and Rare-Earth Metals on the Structure...
The obtained quantitative assessment of mechanical and electrophysical properties and
regularities of their changes depending on process parameters were used to create equipment design
and develop operating practices for producing of long-length products made of the investigated
aluminum alloys.
References
[1] Sidelnikov S.B., Dovzhenko N.N., Zagirov N.N. Combined methods of treatments of nonferrous metals and alloys: a monograph. M.: MAKS Press, 2005. 344 p.
[2] Mondolfo, L.F. Structure and properties of aluminum alloys / Trans. from English edited by F.
I. Kvasov, G.B. Stroganov, I.N. Friedlyander. M.: Metallurgy, 1979. 640 p.
[3] Belyaev, S.V., Dovzhenko N.N., Sidelnikov S.B. et. al. Increase of efficiency of aluminum alloy
profile production on the basis of management of thermal conditions during pressing // The Journal of
Siberian State University. Engineering & Technologies. 2009. Vol. 4. P. 418-426.
Исследование влияния никеля
и редкоземельных металлов на структуру
и свойства литых
и деформированных полуфабрикатов
из алюминиевых сплавов
С.Б. Сидельникова, В.Н. Баранова,
Е.С. Лопатинаа, Л.П. Трифоненковб,
a
Сибирский федеральный университет
Россия 660041, Красноярск, пр. Свободный,79
б
Компания «Русский алюминий»,
Россия 660111, Красноярск, Пограничников, 37/1
Разработана технология приготовления сплавов алюминия электротехнического назначения с
различным содержанием редкоземельных металлов и никеля. В результате экспериментальных
исследований получены литые и деформированные полуфабрикаты из новых сплавов.
Изучены структура, механические свойства металла, а также проведены замеры удельного
электросопротивления на образцах прутков и проволоки, произведенных из них. Установлено
влияние содержания никеля и РЗМ на указанные характеристики, и определены составы
сплавов с наиболее выгодным сочетанием прочностных и электрофизических свойств.
Ключевые слова: алюминий, никель, редкоземельные металлы, сплав, механические свойства,
удельное электросопротивление, структура, параметры литья.
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Journal of Siberian Federal University. Engineering & Technologies 1 (2012 5) 28-38
~~~
УДК 532: 538.4:536.24:621.36: 669.187
Control of Convective Flows
in Liquid Metal in Fluenced
by Electromagnetic Forces
Viktor N. Timofeev*, Аlexander I. Korchagin,
Evgeniy A. Pavlov and Nikolay V. Timofeev
Siberian Federal University
79 Svobodny, Krasnoyarsk, 660041 Russia 1
Received 6.02.2012, received in revised form 13.02.2012, accepted 20.02.2012
During preparation of metal alloys stirring of the melt plays an important role in order to even its
chemical composition and temperature in the whole bath volume. It is necessary to create specific
convectional streams of metal in the bath for effective stirring. Based on the known analysis of selforganization of convectional streams in the liquid heated from beneath (Rayleigh-Bernard effect), the
present article describes a MHD-stirrer which allows receive a similar pattern of streams in the molten
melt by using electromagnetic force. We have considered stirring processes in the rectangular bath
and a MHD-stirrer which allows increase efficiency and reduce energy costs on the alloy preparation
process.
Keywords: structure flows, liquid metal, MHD-stirrer, self-organization, control of flows.
Introduction
The end of the 20th century saw great development of high-intensive technologies, based on
self-organizing processes of dissipative structures. A system is called self-organizing if it develops
a special, temporary or functional structure without a specific outer influence. We understand outer
influence as some influence which is followed by the corresponding modifications of the system. For
example, when feeding energy to the system reaches some critical value and instability occurs, the
system develops a new state.
A laser has become a paradigm for self-organizing coherent processes because accordingly
interacting atoms in a laser lead to appearance of a coherent laser wave.
There is another example of appearance of a macroscopic structure as a result of self-organization
known in the hydrodynamics. When a layer of liquid is heated from the bottom, characteristic, spacetime structures – Bernard’s cells occur. When liquid heating is more intensive, new space-time
structures can occur – oscillations /1/.
The present article attempts to observe application of self-organizing effect in technological
processes of alloy preparation with usage of MHD-technology.
*
1
Corresponding author E-mail address: viktortim0807@mail.ru
© Siberian Federal University. All rights reserved
# 28 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Viktor N. Timofeev, Аlexander I. Korchagin… Control of Convective Flows in Liquid Metal in Fluenced...
As it is known, MHD-technology is based on the influence from an electromagnetic field on melts
with a special purpose. For example, in the baths of melting furnaces and out-of-furnace metal treatment
devices with liquid metal, forced stirring of the melt can be carried out in order to intensify the heat-mass
exchange by influencing the melt with certainly distributed and directed electromagnetic forces.
Under certain heating conditions of the liquid system, actuation of self-consistent dissipative
structures occurs with an over-critical value of the leading parameter which leads to drastic increase
of the heat-mass exchange intensity. It is possible to design casting-melting devices with MHD-stirrers
so that to take into account possibilities of self-organization in a liquid metal in order to reach the end
result with minimal energy demands and automation of the technological process.
Stirring of a melt in casting-melting devices
One of the main ways to improve melting processes and refining of casting alloys is intensification
of heat-mass exchange in a liquid bath of melting furnaces and out-of-furnace metal treatment devices
with the help of forced melt stirring. Stirring provides maximum temperature gradient and concentration
on the border of the molten metal with the interacting phases (metal, alloying and refining hard and
liquid reagents, atmosphere and vacuum). The gradient occurs thanks to supply of new portions of nonreacting metal to those borders and to decrease of the laminar border layer thickness, mass transfer,
which is determined by the coefficients of diffusion and heat conductivity, the values lower than the
corresponding coefficients with convectional heat transfer. In a number of cases the metal movement
induces stirring of the reagent and speed of the heat transfer which can become limiting on a certain
stage. Intensification of the process can also happen due to increase of the phase contact surface when
stirring by creating waves, shredding metal or reagent for drops e.t.c.
Possibilities to use mechanical stirring are quite limited by high melt temperatures and their
large chemical activity. Therefore in the last decades non-contact MHD-devices for force influence on
molten metal have been widely spread.
The staff of the chair “Electrotechnics and electrical engineering” of Siberian Federal University
(Krasnoyarsk) have designed and implemented different MHD-devices of metallurgical purposes in
industrial enterprises: MHD-stirrers in furnaces and mixers; MHD-stirrer in crucibles; MHD-stirrers
of a liquid centre of crystallizing ingots; Device for casting billets with a MHD-mould; MHD-rotators
of liquid metal in channels of induction channel furnaces /2/.
In all the devices force impact of pulsing, running or rotating magnetic field is used on molten
metal in order to receive certain technological effects. In the process of alloy preparation it is necessary
to reach evenness both in chemical composition and temperature in the whole melt volume.
In order to carry out the technological process of alloy preparation the energy is needed, both
to provide the necessary temperature rate and to provide stirring. A way and values of the heat and
electromagnetic energies imposed on the melt define the behavior of the system, enable it to find its
structure which corresponds to the distribution of heat, gravitation and electromagnetic fields as well
as to the geometry of melt bath.
Casting-melting equipment on metallurgical plants
At the present time metallurgical plants employ different furnaces and mixers for preparation of
alloys.
# 29 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Viktor N. Timofeev, Аlexander I. Korchagin… Control of Convective Flows in Liquid Metal in Fluenced...
Fig. 1 Electrical resistance furnace for preparation of aluminum alloys with a rectangular bath
The most popular ones are resistance furnaces, induction crucible and channel furnaces /7/. Fig. 1
shows a schematic of a rectangular resistance furnace for preparation of aluminum alloys. Here 1 –
molten metal, 2 – lining, 3 – steelframe, 4 –electric heaters installed under the furnace roof, 5 – electric
hetaers installed in the furnace hearth.
When electric heaters 4 are installed only under the furnace roof, the heat stream coming from
them reaches the melt mainly due to radiation and thus the upper melt layer has a larger temperature
than the lower layer. With the melt height h about 1 m the temperature difference ∆Т between the upper
and the lower melt layers can reach 1400С. The heat inside of the melt goes from the upper layer to the
lower one mainly due to heat conductivity. In order to even the temperature in the whole melt volume,
it is necessary to carry out its forced stirring. If electric heaters 5 are installed only in the furnace
hearth/8/, the heat will be delivered to the lower melt layer due to heat conductivity. With a certain
temperature difference ∆Т a convectional stream occurs and the heat goes to that bath both due to heat
conductivity and convection.
When electric heaters are installed both under the furnace roof and in the furnace hearth, heating
will be combined and we can control electric heating both from the bottom and from the upper melt
surface.
An example of a more wide spread casting-melting equipment with the melt heated from the lower
surface are induction crucible furnaces (ICF).
Fig. 2 shows a ICF with a cylindrical bath. Here 1 – melt; 2- refractory; 3 – steelworks; 4 –
inductor winding; 5- inductor core; 6 – channel part with molten metal. When the inductor winding is
connected to the varying voltage net, electric current occurs which heats the melt and the heat from the
bath lower part is emitted to its upper layers with heat conductivity and convection.
When using MHD-stirring both in rectangular furnaces and cylindrical ones electromagnetic
forces which significantly influence heat-mass transfer and convectional stream structure will
occur.
A convectional stream has a certain structure and depends on the temperature drop between the
upper and lower layers, liquid characteristics and bath geometry.
# 30 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Viktor N. Timofeev, Аlexander I. Korchagin… Control of Convective Flows in Liquid Metal in Fluenced...
Fig. 2 Induction channel furnace with a cylindrical bath
Analysis of self-organization in a layer
of liquid when heated from the bottom
All the liquids can be divided into liquids characterized by Prandtl number based on their influence
on heat exchange conditions /4/.
Prandtl number is the ratio
Pr
Q
,
k
(1)
where ν – cinematic viscosity; k- thermal conductivity.
In liquids with a big value Pr (all those are non-metallic dropping liquids) processes of movement
quantity transfer prove to be more significant that heat transfer processes.
In liquids where Pr << 1 (all molten metals) heat transfer processes are more intensive that
processes of movement quantity transfer.
The convection in horizontal liquid layers heated from the bottom has been studied since 1900. In
the following years many experiments were carried out which were named experiments of RayleighBernard with a high and low Prandtl number.
An example of self-organization in a drain ditch with liquid with Pr > 1, heated from the bottom
and cooled from a free surface is shown in Fig. /4/. Because of heating and cooling between the lower
and upper liquid surface, there is a certain temperature difference ∆Т = Т1- Т2 which occurs. As the
experiment shows, this temperature difference surprisingly influences the liquid behavior. If the
temperature difference ∆Т between the lower and upper surfaces lower than a critical value ∆Ткрит,
liquid is in a state of peace and heat is transferred only with heat conductivity. When the temperature
difference is larger than critical (Т1- Т2)› ∆Ткрит macroscopic liquid movement suddenly becomes
visible. For example, in the form of bowls, as in Fig.3, if the form of a drain ditch is rectangular or in
a form of regular hexagonal cells, such cells have received the name of Bernard’s cells. We can say
that the parameter ∆Т governs a macroscopic system behavior; therefore ∆Т is called a governing
parameter.
# 31 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Viktor N. Timofeev, Аlexander I. Korchagin… Control of Convective Flows in Liquid Metal in Fluenced...
z
y
Vz
Vy
ly
Vx
h
/y
x
Lx
Fig. 3 Convectional streams in a rectangular bath
q
Formationof
cells
Atrest
0
'Tɤɪɢɬ.
'T
Fig. 4 Dependence of heat stream on the lower surface to the upper one of the temperature difference ΔT
Dependence of the heat flow q per a period of time from the lower surface to the upper one
depending on ∆Т is presented in Fig. 4.
When the value of the temperature difference is critical, the mode of motionless liquid is unstable
(dashed line) and it is replaced with a stable mode characterized with presence of convectional cells.
With a large value of ∆Т fluid at rest is not able to handle transfer of the corresponding large
amount of heat and therefore a convectional mode is steadied which is more favorable for heat transfer.
In comparison with an even balanced distribution, the convection cells are a more high-organized
structure occurring as a result of cooperative movement of liquid molecules.
Physical explanations of the process are as follows. Because the liquid near the lower surface has a
lower density that the liquid near the upper surface due to thermal expansion, due to gravity force and
consequently buoyant Archimedean force, “light lower layer” and “ heavy” upper layer tend to replace
each other. As a result of the liquid viscosity there is no movement even with minor ∆Т.
# 32 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Viktor N. Timofeev, Аlexander I. Korchagin… Control of Convective Flows in Liquid Metal in Fluenced...
When the temperature gradient reaches a critical value, a convectional stream occurs which has a
special structure in the shape of cylinders or hexagonal cells.
A non-dimensional criteria which is proportional to the ratio of buoyant force to viscosity force
is a Rayleigh number
Ra
a ˜ g ˜ h 3 'T
,
Qk
(2)
where α – coefficient of bulk expansion, g – acceleration of gravity.
Following that we can say that when a critical value of Rayleigh number is reached, there occurs an
organized structure in the heated liquid as a result of cooperative movement of the liquid molecules.
The results of Rayleigh-Bernard convection (silicon oil) with a high Prandtl number in a rectangular
bath are presented in /5/. A layer of the observed liquid was placed between a lower polished copper
plate and an upper mono-crystal sapphire plate with high heat conductivity. The convectional structure
was visualized on the horizontal display with a shadow method.
If Rayleigh’s number Ra (proportionally applied to the temperature difference layer in a vertical
direction) exceeds critical value of Raс, there is a convectional movement in the liquid in the shape
of rotating cylinders. In the layer limited by the planes of a rectangular frame of axes (at that in a
horizontal plane Lx > Ly, and the width Ly is not too big in comparison with the depth h of the layer),
the convectional structure consists of straight cylinders with the axes perpendicular to a larger side Lx.
This periodic structure has a certain wave length Λ x, which is close to 2h neat the threshold Raс. The
measures showed that with Ra >Rac (but near Rac ) when we neglect edge effects, the speed changes in
space as in Fourier component
V
Vm cos
2Sx
,
/x
(6)
where х – a coordinate in the direction perpendicular the axes of cylinders.
Thus for a rectangular vessel with liquid, an orderly system of cylinders is the most preferable
convectional structure.
When Ra is increased to a certain value, the convectional stream stays a bidimensional one
(Vx, Vz, Vy = 0), however higher special harmonics V2 with a period Λ x/2 and V3 period Λ x/3 are added
to the main harmonics V with a period Λ x.
New bifurcation is higher than a certain threshold Ra when a third speed component Vy appears,
the convectional structure become a three-dimensional one.
The results of experimental study of the convection in the mercury-liquid with a low Prandtl
number (Pr = 0,031 при 273 K) are presented in /6/. The convection was investigated both in a cylindrical
vessel and in a parallelepiped with the dimensions:
D = 1,7 cm, h = 0,85 cm, Lx =1,6 cm, Ly =0,8 cm, h = 0,8 cm; here D – diameter, h – vessel height,
Lx – large horizontal size, Ly – smaller horizontal size.
Temperature difference on the layer thickness was created by heating on the lower copper block
with constant current.
When analyzing the breaks of the dependence of Nusselt number on Ra the following conclusions
are made. The first break with Ra ≈ Rac is connected with occurrence of convection. For liquids with
# 33 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Viktor N. Timofeev, Аlexander I. Korchagin… Control of Convective Flows in Liquid Metal in Fluenced...
low Prandtl number the main instability of the cylinder bidimensional system is oscillatory instability
which increases buoyant force. The second break is observed with Ra ≈ 2,5 Rac. It can be connected
with occurrence of the second frequency in a spectrum of temperature. In the marginal case of low
Prandtl numbers a mode occurs instead of balance between buoyant force and viscous dissipation
where buoyant force is balanced with inertial members. It happens with Ra ≈ 4 Rac.
Summing up the results of the described experiments we can make a conclusion that with certain
influences on molten metal self-organizing convectional streams can occur in it. Their structure depends
on the melt characteristics, bath geometry and temperature gradient value. Knowing the structure of
convectional streams we can design MHD-devices in such a way that the created electromagnetic
fields will intensify processes of heat-mass transfer or provided energy efficient stirring when there is
no temperature drop.
The convectional movement of molten metal in the shape of cylindrical bowls in a rectangular
bath can also be obtained with the help of electromagnetic forces.
Mathematical and experiment Results
Fig. 5 shows a rectangular bath 1 with molten metal 2. If, for example, a large side Lx≈5 h where
h – metal height, it is possible to install 3 coils. Three consequently connected coils inclined to the right
present one phase A and three consequently connected coils inclined to the left present another phase
B. When the phases A and B are connected to a two-phase power supply with a phase shift of 900,
alternating currents appear in the coils iA=Imsinωt; iВ=Imsin(ωt+900); where Im and ω – is an amplitude
and angle current frequency. We know from the electrotechnics that two coils space-shifted and fed by
the phase-shifted currents create magnetic fields.
As a result of the influence of rotator magnetic fields on the melt convectional streams in the form
of rotating cylinders occur.
A system of three pairs of coils can be replaced with an inductor containing a magnetic wire 3 and
six coils wound on the magnetic wire frame. Here each phase consists of three consequently connected
L
h
x
x
h
h
h
x
x
h
x
x
1
2
iB
iA
x
iA
x
x
iA
x
x
iB
iB
Fig. 5 Convectional MHD – movement in the rectangular bath
# 34 #
iB
x
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Viktor N. Timofeev, Аlexander I. Korchagin… Control of Convective Flows in Liquid Metal in Fluenced...
̌
b
d
c
Fig. 6. Pathway with a different frequency: а) f=50 Hz; b) f=800 Hz; c) f=1500 Г=Hz; d) f=10000 Гц
coils. The inductor can be installed under the bath with a melt. Such a construction has an advantage
as it does not include the coils located in the bath with a larger temperature.
A mathematic and physic models have been built to investigate a pattern of flows in a rectangular
bath when an inductor electromagnetic field influences the melt, installed under the bath.
Bath manufactured from stainless steel has the following dimensions (mm): Lx = 125; Ly = 40;
h = 25.
As liquid metal gallium was used with heat conductivity γGa ≈3, 106 Om-1·m-1 and the melting
temperature 38 °С.
Mutual solution of electrodynamical and hydrodynamical solutions was carried out with a program
package ANSYS Academic Research.
Fig. 6а, 6b, 6с and 6d show the results of numerical simulation for the molten metal convectional
streams with a different inductor current frequency. As it follows from Fig. 6а with f = 50 Hz metal
movement is presented by rotating cylinders. With an increase of frequency Fig. 6b (f = 800); Fig. 6с
(f = 2000) the cylinders begin to deteriorate, slowly the movement of distinct eddies becomes a complex
turbulent movement (Fig. 6d, f =10 000).
This conclusion is also confirmed by Fig. 7, which shows the graphs of distribution for the velocity
component Vx along the bath length у of the melt bottom (z≈0). When with frequency f<1000Hz
velocity curves are close to the sinusoid , with a period Т≈2h, then with the frequency f>1000Hz higher
# 35 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Viktor N. Timofeev, Аlexander I. Korchagin… Control of Convective Flows in Liquid Metal in Fluenced...
Dependence of Velocity u-comp in point x=0,0375 m
0,0070
0,0060
0,0050
0,0040
0,0030
0,0020
0,0010
0,0000
-0,0010
0
2 000
4 000
6 000
8 000
10 000
12 000
-0,0020
-0,0030
Fig 7. A graph of distribution of the velocity compound V on the bath length у of the melt bottom
Fig.8 Dependence of the velocity V on the frequency f
harmonic appears. Fig. 8 shows dependence Vx on the frequency, in the point: х =3/2h ; z≈0. Maximal
value of the velocity is reached with the frequency f≈900 Hz, with a further increase of the frequency,
the velocity first decreases and after changes a sign.
All the calculation have been carried out with a constant current value in the coils I = 5А. As it
follows from the above calculations, the current frequency governs the pattern of the convectional
streams and thus it can be called a governing parameter. Obviously we can govern the pattern of the
convectional streams with the current value I on the inductor coils with constant frequency.
# 36 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Viktor N. Timofeev, Аlexander I. Korchagin… Control of Convective Flows in Liquid Metal in Fluenced...
By an experimental investigation it was possible to prove existence of bowls at the inductor
winding feeding with the currents of 50Hz to 800Hz frequency. At that a propeller with flat planes of
stainless steel was inserted into the melt area, the propeller diameter is a few centimeters lower than
h. A fine wire was welded to one of the vanes. The propeller was lowered under the oxide film of the
molten gallium and the angle velocity was defined by the wire movement. The measurements of the
velocity were approximate as the wire passing the oxide film experiences large resistance.
Together with that we registered change of the rotation direction for bowls and increase of the
rotational speed with increase of frequency.
Thus by changing the current frequency of the inductor coil (governing the parameters) we can
obtain a different structure of the convectional streams in the melt leading to effective stirring.
Conclusions
1. One of the main directions in updating the technology of alloy preparation is intensification
of convectional heat-mass exchange in a melt. The most effective dynamics should embrace maximally
possible set of hydrodynamic structures for the corresponding macroscopic and microscopic stirring
along the whole bath volume.
2. Taking into account the melt characteristics, thermal strains, bath geometry we can design
MHD-devices which will strengthen the effect of self-organization of dissipative structures by the
influence of electromagnetic fields on the melt and will enable us to create energy-efficient castingmelting equipment.
References
[1] Haken H. Synergetics. Springer-Verbag Berlin Heidelberg New Yorke 1978, 405p.
[2] Timofeyev V., Khristinich R., Boyakov S., Pervukhin M., Golovenko E., Nadtochiy D.
Application of MHD-technology in alluminium alloy production. Proceedings of the International
Symposium on Heating by Electromagnetic Sources. HES-07, Padua, June 19-22, 2007, pp.145-154.
[3] Shmitz C. (ed.) Handbook of Aluminium recycling, vulkan verlag, 2006, 510p.
[4] Ebeling W. Strukturbildung bei irreversible prozessen. Elne Elnfukrung in die T heovie
dissipativer Strukturen. Bsb b. g. Teubner Verlagsgesellschaft 1976, 279p.
[5] Berge P. Rayleigh-Benard Convection in High Prandtl Number Fluid. – in: Chaos and Order
in Nature. Proc. Int. Symp. On Synergetics at Schloss Elmau, Bavaria, April 27-May 2, 1981/Ed. H.
Haken. – Berlin, Heidelberg, New York: Springer-Verlag, 1981, pp. 14-24.
[6] Fauve S., Libchaber A. Rayleigh-Benard Experiment in a Low Prandtl Number Fluid,
Mercury. – in: Chaos and Order in Nature. Proc. Int. Symp. On Synergetics at Schloss Elmau,
Bavaria, April 27-May 2, 1981/Ed. H. Haken. – Berlin, Heidelberg, New York: Springer-Verlag,
1981, pp. 25-35.
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Viktor N. Timofeev, Аlexander I. Korchagin… Control of Convective Flows in Liquid Metal in Fluenced...
Управление конвективными потоками
в жидком металле
под действием электромагнитных сил
В.Н. Тимофеев, А.И. Корчагин,
Е.А. Павлов, Н.В. Тимофеев
Сибирский федеральный университет
Россия 660041, Красноярск, пр. Свободный
В процессе приготовления металлических сплавов перемешивание играет важную роль в
выравнивании химического состава и температуры во всем объеме ванны. Для эффективного
перемешивания необходимо создать специфическую структуру конвективных потоков.
На основе анализа известного эффекта самоорганизации конвективных потоков в
жидкости, подогреваемой снизу (эффект Релея – Бенара), в данной статье описан МГДперемешиватель, который позволяет получить аналогичные потоки в жидком металле под
действием электромагнитных сил. Рассмотрены процессы перемешивания в прямоугольной
ванне, предложен МГД-перемешиватель, который увеличивает эффективность и снижает
энергопотребление в процессе приготовления сплава.
Ключевые слова: структура потоков;
самоорганизация; управление потоками.
жидкий
металл;
МГД-перемешиватель;
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Journal of Siberian Federal University. Engineering & Technologies 1 (2012 5) 39-56
~~~
УДК 621.316.11
Stochastic Simulation of Covariance Matrix
and Power Load Curves
in Electric Distribution Networks
Ivan V. Shulgina,
Aleksey A. Gerasimenko * and Zhou Su Quanb*
a
Siberian Federal University,
79 Svobodny, Krasnoyarsk, 660041 Russia
b
Harbin Institute of Technology
China 150001, Harbin, 92 West Dazhi Street,
Nan Gang District 1
a
Received 6.02.2012, received in revised form 13.02.2012, accepted 20.02.2012
An algorithm of stochastic simulation of covariance matrix and nodal power load curves is developed
for electric distribution networks based on factor analysis. Statistical stability of factor power load
model is confirmed. Application of this model is able to identify a general regularity of nodal power
changing, and to simplify the analysis of multivariate operating conditions in operational problems of
electric distribution networks and their optimization.
Keywords: stochastic simulation, electric power load, covariance matrix, electric distribution network,
principal component analysis, energy losses, operating condition optimization.
Introduction
The adoption of automated meter reading (AMR) systems in the industry makes it possible to store
statistical data about power transmission and consumption. Based on the above-mentioned systems
and modern mathematical methods, it is possible to solve a series of problems: multifactor simulation,
prediction and standardization of energy consumption and some integral characteristics of power
systems; production activity analysis and optimization of power system functioning; diagnostics of
electrical equipment in electric power supply systems etc. [1–3].
Power supply continuity and safety of electric power supply depend on a stability of a whole chain:
“electric power generation – transmission – distribution”. Electric distribution networks, which are the
master link in that chain, are the most problematic ones and outlay elements influence not only the electricity
tariffs, but also the economic efficiency. About half of the current power sector’s basic assets are related
to electric distribution networks, and most of the energy is lost just in these networks. However, the role
of the electric distribution networks is still often underestimated on the background of global construction
problems. Thus, dangerous and far reaching consequences, both economic and social, may arise [4].
*
1
Corresponding author E-mail address: gerasimenkoaa@yandex.ru
© Siberian Federal University. All rights reserved
# 39 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Ivan V. Shulgin, Aleksey A. Gerasimenko... Stochastic Simulation of Covariance Matrix and Power Load Curves...
Recently, taking into account the new computer technologies and the development of modern
control measuring systems, the models of power consumption have been mainly developed by
means of stochastic methods of component analysis, which include the principal component method
[1, 5–12]. The models of power consumption or, in other words, the models of power load curves
are able to decrease the volume of initial information needed for problem solving, and to simplify
analysis of multivariate operating conditions in electric distribution systems. The problems include
determination of power integral characteristics (energy losses, ranges of changing the operating
condition parameters in the electric nodes and between power systems etc.), and reactive power
compensation, both of which are very important when it comes to the complex optimization of
power system and energy saving.
Deterministic methods, statistical simulation of power operating conditions, and determination
of integral characteristics in power engineering have been advocated by different authors for some
years [1, 5–8, 10–12]. The related research studies faced some difficulties: large dimension of the nodal
power covariance matrix, large volume of information about power loads and operating condition
parameters, complicated processing of initial information as well as underdevelopment of measuring
systems, computers and programming. Therefore, the application area of stochastic analysis methods
was limited. More recently, taking into account the adoption of SCADA and AMR systems, the
above-mentioned disadvantages have been gradually disappearing, and the development of stochastic
methods, partially the methods based on the principles of factor analysis, is more promising [13].
Conceptual Description of Principal Component Method
A component analysis as a method was developed by Pearson1; he proposed a method of databank
compression which allocates a maximal variance. This method was also developed by Hotelling2.
The factor analysis is a multivariate analysis which researches an internal structure of the
covariance or correlation matrices. It is applied for the statistical research of a system of random
variables which have a correlation by means of stable random or nonrandom factors [9].
The principal component method is based on simple and ordinary conceptions, which depend on
the covariance matrix analysis and the matrix linear transformation.
The modeling, characterizing the behavior of a random variable, is implemented by different ways
of regression and factor analysis. In regression analysis the factors and model structure are entered a
priori; in factor analysis we assume that, factors exist while their exact number and the model structure
are only determined during the process of problem solving.
The principal component method is a dismemberment of a covariance matrix on the orthogonal
vectors (components) or directions corresponding to the number of variables. These vectors correspond
to the eigenvalues and the eigenvectors of the matrix. We agree that by a characteristic value we mean
the set of eigenvalues and eigenvectors of the matrix.
Based on this method, the characteristic values are formed in descending order, which is important
since only few components have to be used for the description of the initial data. The vectors are
pairwise orthogonal ones, and their components are uncorrelated. A few components can reflect most
1
2
Pearson, K. On lines and planes of closest fit to systems of points in space. – Phil. Mag. – 1901. #6, p. 559–572.
Hotelling, H. Analysis of complex of statistical variables into principal components. – Jep, #24, 1933. – p. 417–441,
498–520.
# 40 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Ivan V. Shulgin, Aleksey A. Gerasimenko... Stochastic Simulation of Covariance Matrix and Power Load Curves...
of the sum variance of initial variables; however, all components are required for accurate reproduction
of correlations between variables.
The principal component method is used for total simulation of initial random variables. However,
we do not need to put forward the hypotheses about variables because the variables do not even have
to be random variables. In practice the observations of random variables are samples from some
population.
In order to decrease the complexity of the statistical calculations we can replace an n-dimensional
random variable by k<n linear functions from the initial variables. The simulation is called a
reconstruction of function using a linear predictor3, which is implemented by means of eigenvectors of
the covariance matrix [12].
Consider the multivariate random variable X which is an n-dimension sample data
X
ª x11 x 21 ... x n1 º
« x x ... x »
n2 »
« 12 22
«... ... ... ... »
«
»
¬ x1m x 2 m ... x nm ¼
For the analysis of other random variables depending on X, it is necessary to determine the
mathematical expectations, for instance, sample mean or average values МХ1, МХ2, …, МХn, and
variations (changing) of initial random variables in the neighborhood of their average values ΔХ1, ΔХ2,
…, ΔХn. A characteristic of a random variable variation in the neighborhood of their average values is
a variance. It may be that the linear combination of initial random variables has the maximal variance
(certainly, we should compare only the normalized linear combination of random variables because
any random variable can be multiplied by a large number, so any large variance can be obtained).
The linear transformation of initial variables is implemented by means of uncorrelated and
normalized linear variables υ.
Linear Combination of Random Variables
with Maximal Variance
Consider the possible linear combinations of random variables X i
G
Xuȣ
where
>X ... X @u >X ...X @,
1
ª x11 º
«x »
« 12 »
« x1i » ; X k
« »
«... »
« x1n »
¬ ¼
X1
k
1
(1)
k
ª xk1 º
«x »
« k2 »
« x ki » ; X1
« »
«... »
« x kn »
¬ ¼
ªX11 º
« »
«X12 »
«X1 j » ; X k
« »
«... »
«X »
¬ 1m ¼
ªX k 1 º
« »
«X k 2 »
«X kj » ,
« »
«... »
«X »
¬ km ¼
taking into account a limitation which fulfills the following condition of normalization
m
¦X
2
ij
1,
i 1,2,..., k
(2)
j 1
3
Predictor is a superior system of variables
# 41 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Ivan V. Shulgin, Aleksey A. Gerasimenko... Stochastic Simulation of Covariance Matrix and Power Load Curves...
The variance of the linear combination (1) is determined by the following formula based on [7,
12]
ı 2G
DG
ªX11
«X
« 21
«...
«
¬X k1
X12 ... X1m º
ªX11 X 21 ... X k1 º
»
«X
X 22 ... X 2 m »
X 22 ... X k 2 »»
u K(X) u « 12
...
Xk 2
... ... »
»
... X km ¼
ȣ ɬ u K(X) u ȣ
«... ... ... ... »
«
»
¬X1m X 2 m ... X km ¼
ªO1
« 0.0
«
«...
«
¬ 0.0
(3)
0.0 ... 0.0 º
O2 ... 0.0 »»
,
... ... ...»
»
0.0 ... O k ¼
where K(X)=K – covariance matrix of initial random variables Х1, Х2, …, Х k;
k – rank of matrix K(X);
т – index of the transpose of a matrix;
m – total number of changing of the random variable Хi.
We assume that the matrix’s elements k(X i X j) are the estimations calculated from samples: xi1,
xi2, …, xim and xj1, xj2, …, xjm. Thus, the selection of a random factor having a maximal variance is
found through a minimum of function (3) satisfying condition (2). The optimization problem is solved
by a method of Lagrange multipliers. Introduce an auxiliary objective function, which is Lagrange
function
Ɏ
m
V 2 G l1 ¦ (X12j 1) ,
(4)
j 1
where l1 – the Lagrange multiplier.
An absolute minimum of the function (4) corresponds to the conditional minimum of function (3)
subject to condition (2). The function including all variables is differentiated; the minimum condition
is obtained as
wɎ
wX r
m
2¦ K jrX1 j 2l1X1r
m
0 ;r
1,2,..., k ;
j 1
¦X
2
1j
1.
(5)
j 1
The solutions of system (5) are all normalized eigenvectors of the matrix K(X). Every solution
determines an extreme point or specific point of the function. The coordinates of the eigenvector
corresponding to the maximal eigenvalue λ1 correspond to the global minimum.
In factor analysis the components of the vector G are new random variables which are the linear
combination of initial X or centered ΔХ random variables.
The eigenvalues λ and eigenvectors X of power covariance matrix have useful properties
which are applied in component analysis. The eigenvalues are real values and the eigenvectors can
be chosen as perpendicular to each other. The eigenvectors defi ne an undergoing pure tension or
compression direction of the linear transformation corresponding to the matrix K(X). These vectors
are also named the principal components of the matrix [12], and the eigenvalue λ is a coefficient
of the transformation. The variance of the i-th principal component is equal to the eigenvalue λi of
matrix K(X).
It is known that the eigenvalues λ and eigenvectors X of matrixes satisfy the following
expression
# 42 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Ivan V. Shulgin, Aleksey A. Gerasimenko... Stochastic Simulation of Covariance Matrix and Power Load Curves...
Kuȣ
ȣu Ȝ.
(6)
1
By multiplying both sides of the expression (6) on the left of the matrix ȣ , we will arrive at
Ȝ
ȣ 1 u K u ȣ .
(7)
Expression (7) is considerably simplified when the original matrix K is defined as positive, which
is the case for the power covariance matrix [7]. Thus, all eigenvectors can be made orthonormal ones,
i. e. satisfying the expressions
Xiɬ uX j
0 when i z j ; X i ɬ u X j
1 when i
j.
It is easy to verify that the inverse matrix ȣ 1 is equal to the conjugate one, and expression (7) can
be rewritten as (3)
Ȝ
ȣɬ u K u ȣ ,
(8)
where ȣ – an orthonormal matrix whose columns consist of eigenvectors X1 , X 2 , ..., X k . Condition (2) is
satisfied for every column of the orthonormal matrix.
Inverse expression between the original matrix K and the matrix λ
K
>X1
X 2 ... X k @u Ȝ u >X1 X 2 ... X k @ɬ
ªX11 X 21 ... X k1 º ªO1 0.0 ... 0.0 º ªX11
«X X ... X » «0.0 O ... 0.0 » «X
k2 »
2
« 12 22
» u « 21
u«
«... ... ... ... » «... ... ... ...» «...
«
» «
» «
¬X1m X 2 m ... X km ¼ ¬0.0 0.0 ... O k ¼ ¬X k 1
,
X12 ... X1m º
X 22 ... X 2 m »» .
...
Xk2
(9)
... ... »
»
... X km ¼
The method allows us to select the orthogonal factors among the factors-arguments, i. e.
statistically-independent components, which provide a linearity of the method and additive efficiency.
The above-mentioned properties of the eigenvectors show that the full totality of them is
equivalent to the original probabilistic model corresponding to the vector X. Moreover, both sets
of variables X and G define the same vector space.
However, from the set of vectors G, it is sufficient to select a small number of M principal
components (factors) explaining most of the relationships between all components of the initial vector of
random variables X. The main factors are not directly observed, but they characterize the change in the
original variables. Therefore, we can get the task of obtaining a linear predictor of dimension M (M <k).
For each k, the best predictor is the first M eigenvectors of K, corresponding to the maximal eigenvalues.
As a result of studying the internal structure of the matrix K, the selection of main factors is made
in such order that at the beginning the first of them makes the greatest contribution to the variance
of the variables, then the second one is the largest contribution to the variance of the variables
remaining after taking into account the main factor, etc. Ultimately, these vectors constitute a set
of linearly independent basis vectors, oriented in such a way that each of them makes the maximum
contribution to the variance of the original variables Х.
On the practical side, the factor model makes it possible to adequately estimate the covariance
structure between the relatively large number of observed variables by means of a smaller number
of common factors. Evaluation of the factor structure is carried by the required number of factors
explaining the correlations between variables and load factors in these variables. Component analysis
# 43 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Ivan V. Shulgin, Aleksey A. Gerasimenko... Stochastic Simulation of Covariance Matrix and Power Load Curves...
is most useful when all variables xi are measured in the same units. If not, the method is much more
difficult to validate [9].
Factorial or component methods of statistical analysis are used in the computation of operational
energy losses [13, 14], as well as for short-term forecasting and optimization [15, 16]. When solving
the problem of factorial simulation of electric loads via a stochastic approach, the information about
the characteristics of the random variable is approximately determined by a partial sample from the
general population. In the factor simulation of power loads [5–7, 10, 11] the curves of active and reactive
nodal powers are considered as a training sample having 2n-dimension. In operational computations
the scope of the method is limited to the modeling of daily power load curves of an unobservable
network.
The modeling of electric power loads on the basis of factor analysis allows us to:
– find hidden regularities, which are determined by many internal and external causes of load
changing;
– carry out the compression of information by describing all curves by means of the common
factors or principal components, whose number is much smaller than the number of initial
curves;
– identify the statistical dependence between the power load curves and the main factors;
– predict the random component curves based on the regression equation constructed on the
basis of factor analysis;
– simplify the methods for determining the integral characteristics of power systems.
Determination Methods of Principal Components
The problem of determination of the principal components is a classical problem of determining
the characteristic values from a covariance matrix of random variables, as which the nodal power
loads are considered. The determination of eigenvalues and eigenvectors of matrices in linear algebra
is called the problem of characteristic values, and it is a complicated task which is implemented in
several statistical software applications. The value λ is called the eigenvalue of K, if there is a nonzero
vector (eigenvector of K) satisfying the equation
(K O u E) u X
0,
(10)
where E – unity matrix; 0 – null vector.
The system (10) is a homogeneous system of liner equations because the free members of its
equations are zeros. It has nontrivial solutions if the determinant of matrix K O u E equals zero, i. e.
On E1On 1 E 2 On 2 ... E n 1O E n
0,
(11)
where E1 ,..., E n – coefficients of the characteristic polynomial.
The methods for determining the eigenvalues and eigenvectors can be divided into two groups [11]:
the first group includes iterative methods which often use a similarity transformation and solve a linear
system of equations (10); the second group includes the direct methods that calculate the characteristic
polynomial (11). The problems (10) and (11) have different conditionality, as the roots of polynomials
(11) are often highly sensitive to errors which are inevitably arising in the calculation of polynomial
coefficients. That was the main reason of the almost complete exclusion of direct methods.
# 44 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Ivan V. Shulgin, Aleksey A. Gerasimenko... Stochastic Simulation of Covariance Matrix and Power Load Curves...
The direct application of covariance matrix in various algorithms is greatly complicated by its
dimension. In order to compensate the above-mentioned disadvantage, the modeling of covariance
matrix is implemented by using modern computer software interactive systems, such as MATLAB,
MATCAD, C++, ANSYS, FORTRAN, etc.
The main criterion for the normalization of the eigenvectors in MATLAB consists of
ȣɬ u ȣ
E.
(12)
Small changes in matrix elements, such as rounding errors, can cause large changes in the
characteristic values. The power covariance matrix is a square matrix that is easier to use for matrix
transformations in comparison with other matrices.
Stochastic Model of Covariance Matrix
and Electric Power Loads
Most of the research [5–8] conducted in this field was aimed at the modeling of power loads and
its application for a daily time interval. This is due to the peculiarities of the energy business and
information support of power utilities during the development of this technique. Structural changes that
have occurred in the management of an integrated power grid led to the need of periodic computations
between separate power business entities. Today, the main period of the financial settlement is one
month. The modeling of power consumption on a monthly time interval was proposed for the first time
in [10, 11]. At the moment the technique and its possible application are implemented insufficiently and
therefore further research and elaboration is needed.
The simulation of the power covariance matrix is based on several properties of eigenvalues and
eigenvectors of a matrix which is expanded by 2n eigenvalues and eigenvectors; the first few M characteristic
values (M <<2n) accurately reflect the total variance of the initial power load curves [5–7, 10–12].
Statistical analysis of power operating conditions uses information about variances of power loads
2
V Pi , V 2 Qi and cross-covariance functions k ( Pi Pj ) , k ( Pi Q j ) , k (Qi Q j ) between random variables of
different power nodes
V 2 Pi
1 d
¦ ( Pim MPi ) 2 ; V 2 Qi
dm1
k ( Pi Q j )
1 d
¦ (Qim MQi ) 2 , i 1, n ;
dm1
1 d
¦ ( Pim MPi )(Q jm MQ j ) , i, j
dm1
k ( Pi Pj )
1 d
¦ ( Pim MPi )( Pjm MPj ) , i, j
dm1
k (Qi Q j )
1 d
¦ (Qim MQi )(Q jm MQ j ) , i, j
dm1
1, n ;
(13)
1, n , i z j ;
1, n , i z j ,
where i, j – indexes of nodes; m – index of every time interval for the T period; n – the number of nodes
of power distribution systems with known or simulated power load curves.
The elements (13) of a power covariance matrix characterize the degree of irregularity of power
load curves, which remains approximately constant over a long period. The degree may be determined
on the basis of daily measurements performed on different days in power system sectors. This is an
important advantage of the statistical method.
# 45 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Ivan V. Shulgin, Aleksey A. Gerasimenko... Stochastic Simulation of Covariance Matrix and Power Load Curves...
The variances and cross-covariance functions of power loads form the square covariance matrix
K= K(P,Q) as follows
ª K 11 K 12 º
«K K »
¬ 21 22 ¼
K(P, Q)
ª ªV 2 P1 k ( P1 P2 )......k ( P1 Pn ) º ªk ( P1Q1 ) k ( P1Q2 )......k ( P1Qn ) º
««
» «
»
2
« «k ( P2 P1 ) V P2 ......k ( P2 Pn ) » «k ( P2 Q1 ) k ( P2 Q2 )......k ( P2 Qn ) »
« «.......................................... » «.................................................»
««
» «
»
« «¬k ( Pn P1 ) k ( Pn P2 )......V 2 Pn »¼ ¬k ( Pn Q1 ) k ( Pn Q2 )......k ( Pn Qn ) ¼
«
« ªk (Q1 P1 ) k (Q1 P2 )......k (Q1 Pn ) º ªV 2 Q1 k (Q1Q2 )......k (Q1Qn ) º
»
««
»«
2
« «k (Q2 P1 ) k (Q2 P2 )......k (Q2 Pn ) » «k (Q2 Q1 ) V Q2 ......k (Q2 Qn ) »
« «..........................................
» «.......................................... »
»
««
»«
«¬ ¬k (Qn P1 ) k (Qn P2 )......k (Qn Pn )¼ «¬k (Qn Q1 ) k (Qn Q2 )......V 2 Qn »¼
º
»
»
»
»
»
»
»
»
»
»
»
»¼
(14)
The simulation of the covariance matrix is implemented corresponding with expression (9), which
may also be written in the usual form
M
K
¦O
i
uXi uXiɬ .
(15)
i 1
Every eigenvalue of covariance matrix corresponds to a general load diagram of power (GLD) Ƚ i ,
which is a linear combination of 2n initial nodal power load curves centered at expectation MPi, MQi
Ƚ 2n
>'P1 ...'Pn
'Q1 ...'Qn @u >ȣ 2n @
ɬ
>Ƚ
Ƚ 2 Ƚ 3 ... Ƚ 2 n @, i 1,..., n ,
1
(16)
where >ȣ2n @ɬ – the transposed matrix of eigenvectors obtained from the covariance matrix of statistical
sample for initial power loads, which has 2n×2n dimension;
ΔР1 , ΔQ1 – the deviations centered at expectation of active and reactive power in the node # 1 for
a certain time period T
ª'p11 º
ª'q11 º
«'p »
«'q »
« 12 »
« 12 »
'P1 «'p13 » , 'Q1 «'q13 » , m=1, 2,…,d .
«
«
»
»
«.... »
«.... »
«'p1m »
«'q1m »
¬
¬
¼
¼
The obtained GLDs can be considered as new independent centered random variables with zero
mathematical expectations. The GLDs have a property of orthogonality, i. e. cross-covariance functions
k ( Ƚ i Ƚ j ) , k ( Ƚ j Ƚ i ) are equal to zero. These new random variables (factors) are a suitable coordinate
system for accurate simulation of initial random variables Pi, Qi; therefore, using M of them
Ƚk
>Ƚ ....Ƚ @ Ƚ
1
k
2n
, k
1,..., M 2n ,
(17)
which corresponds to the maximal eigenvalues of power covariance matrix, allow us to simulate initial
load changing with a sufficient accuracy for a certain time interval T
S i | ONE m,1 u >MP1 ....MPn MQ1 ...MQn @ >'P1 ....'Pn 'Q1 ...'Qn @
ONE m ,1 u >MP1 ....MPn MQ1 ...MQn @ Ƚ k u >ȣ k @
ɬ
# 46 #
>P1 ...Pn
Q1 ...Qn @ , k
1,..., n ,
(18)
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Ivan V. Shulgin, Aleksey A. Gerasimenko... Stochastic Simulation of Covariance Matrix and Power Load Curves...
where ONEm,1 – a column vector consisting of units, and which have m = 1, 2,..., d rows.
>ȣk @ɬ – the transposed matrix of prime k eigenvectors X k corresponding to the first maximal
eigenvalues λk of the power covariance matrix K(P,Q) (14);
MQ1 – a mathematical expectation of reactive power curve in the node #1 for the accounting
period T.
Рi – possible variation of active power in i node for the accounting period T.
The simulation of the power loads allows us to track the variation of load parameters. It should be
noted that initial power load curves are fully simulated on formula (18) using all GLDs. The eigenvalues
λ of the initial covariance matrix K(P,Q) are variances of GLDs (3). Hence, any power load curve can
be represented as a linear combination of the GLDs which reflect the general regularities of power
expectation changing for the initial collection of power nodes.
The following daily/monthly samples of the power load curves have been analyzed:
1) 36 simple weekday and weekend power curves [24 hrs] of a 10 kV electric distribution network
(d=4); the unit of active power is kW; the unit of reactive power is kvar.
2) 18 atypical real daily active and reactive power curves of a 10 kV electric distribution network
(d=24) [11]; the unit of active power is kW; the unit of reactive power is kvar.
3) 42 typical daily power curves for different industries (d=12) [12]; the unit of power is relative
unit [r. u.].
4) 30 monthly active power load curves for 110–220 kV overhead lines for August–September in
2009 [11]. Number of intervals d of period T was reduced from 744 to 31 by calculating the average
power for each day. The unit of active power is MW.
Computational results of eigenvalues and GLDs (17) from the above-mentioned samples of data
#1–4 are presented in Table 1–3, and in Fig. 1–2.
The original power load curves characterize a different degree of irregularity, therefore every
sample data #1–4 has its own principal factors linking the power load curves to the system of
characteristic values. In all cases, the error of simulation of power load curves is, in a dozen times
and sometimes more, less than the error of the covariance matrix simulation (9) taking into account
Table 1. The Six Maximal Eigenvalues Obtained from Sample Data #1–4 for Simulation of Covariance Matrix
and Power Load Curves
Sample
Data
1)
2)
3)
4)
Contribution of principal
components to the total
variance of loads
λ
Θ, %
λ
Θ, %
λ
Θ, %
λ
Θ, %
Eigenvalues of the original sample data
in decreasing sequence
λ1
λ2
λ3
45806.73
25494.41
5039.49
60.00
33.40
6.60
17699.71
9289.16
53.28
27.96
1.06705
λ5
λ6
-
-
0.0
-
-
2247.51
1136.81
595.20
502.52
6.77
3.42
1.79
1.51
0.432602
0.104090
0.102429
60.10
24.36
5.86
5.77
1.60
0.885
1510.38
369.29
172.26
80.00
68.83
47.70
65.06
15.91
7.42
3.45
2.96
2.05
# 47 #
λ4
4.64·10
-12
0.0284292 0.0157218
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Fig. 1. Two Principal Daily GLDs Corresponding to the Maximal Variances for Sample Data #1 in the Same Units
as the Initial Variables
Fig. 2. Three Principal Daily GLDs Corresponding to the Maximal Variances Obtained from Sample Data #2 in
the Same Units as the Initial Variables
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Ivan V. Shulgin, Aleksey A. Gerasimenko... Stochastic Simulation of Covariance Matrix and Power Load Curves...
Table 2. Six Daily Nonnormalized GLDs Corresponding to the Maximal Variances Obtained from Typical Daily
Power Curves for Different Industries (Sample Data #3) [r. u.]
Ƚ1
Ƚ2
-1.58281
0.107604
2–4
-1.63002
0.111508
4–6
-1.00645
-0.0537509
6–8
0.826063
0.369511
8–10
0.965881
0.888987
10–12
0.186023
0.684827
12–14
0.671946
14–16
1.33107
t, hrs
0–2
Ƚ3
-0.184982
Ƚ4
Ƚ5
Ƚ6
0.0774817
-0.00315256
0.0536194
-0.193562
0.0402251
0.0274119
0.0268668
0.339956
-0.721240
0.201227
-0.130268
-0.435288
-0.464289
-0.383642
0.0144767
0.169168
0.111077
-0.0132479
-0.241056
0.589090
-0.0899582
-0.0842728
0.267227
0.990801
0.0328418
0.464211
0.1607393
0.00610122
-0.155875
-0.588556
-0.189836
0.303330
0.0769681
16–18
0.832963
-1.12614
0.168255
0.0598246
0.0587576
0.125265
18–20
0.687673
-0.991643
0.289502
0.0722762
-0.101288
-0.107512
20–22
-0.0217858
-0.704435
0.00610245
0.325773
-0.0908384
-0.0676617
22–24
-1.26054
-0.121393
-0.192526
0.314455
-0.0750239
-0.0240258
the same number of characteristic values and GLDs. The reason for this is that the elements of
the power covariance matrix are relatively small values compared with the values of the initial
power load curves. The maximal error of simulation of power load curves does not exceed 25%,
and the average one is 8.83% for sample data #3 taking into account only six maximal characteristic
values.
The obtained GLDs can also be used to determine the total normalized or weight average GLDs
which are used for the simulation of unknown power load curves.
By the means of MATLAB system using all GLDs and characteristic values for the sample data,
the power covariance matrix (9) and power load curves (18) are simulated with high accuracy. The
simulation on the basis of expression (18) is primarily designed for modern automated meter systems
and it uses operating condition information from them. The lack of the above-mentioned systems
in most distribution systems can be replaced by the power load simulation [7, 11] and so to use the
advantages of factor simulation.
Stability of Factor Power Load Model
Factor simulation of power loads has a practical application only if the estimations of factor values
(GLDs), which are obtained for different random processes of load changing, are statistically similar
ones, i. e. have statistical stability. A research study of the statistical stability of the factor model, which
was based on real data on power load curves for different power utilities with large statistical volume of
information, showed the presence of a collective and dynamic stability for daily, weekly and monthly
power load curves [4, 5, 7, 10, 11].
Statistical stability describes the possibility of using GLDs derived from one learning sample for
the simulation of powers which were not included in the learning sample. Dynamic stability describes
the comparison of different time realizations of the factor model for constant power utilities; collective
stability is the comparison of GLDs belonging to different power utilities.
# 49 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Ivan V. Shulgin, Aleksey A. Gerasimenko... Stochastic Simulation of Covariance Matrix and Power Load Curves...
Table 3. Six Monthly GLDs Corresponding to the Maximal Variances Obtained from Sample Data #4 [MW]
d
Ƚ1
Ƚ2
Ƚ3
Ƚ4
Ƚ5
Ƚ6
1
-61.0910
-19.4337
-5.91616
0.987755
-17.5551
2.85404
2
-49.9621
-17.3637
0.624555
6.63003
-10.2079
-7.32189
3
-56.1344
-11.5625
10.9850
9.11488
-9.88028
-11.2969
4
-44.3974
-8.20962
10.8664
5.96288
0.0544416
-17.9550
5
-29.6979
2.44294
-19.8471
8.92353
-2.01733
2.65349
6
7.69728
16.7571
-23.6115
8.33256
-4.15039
9.31310
7
50.0524
42.8409
0.189204
4.71245
-14.7875
-3.14985
8
51.1595
42.7496
-11.1642
-6.21158
-10.4222
-7.43516
9
54.5808
14.4466
-22.6867
-2.33433
6.02504
-11.5762
10
43.7955
-33.5056
-13.8767
4.15772
2.12694
-4.75286
11
42.7563
-28.2635
-14.5850
4.16292
2.03908
-4.62368
12
22.2479
0.475527
-4.00582
-13.7653
-16.9796
3.54289
13
-82.9028
22.9289
8.61478
-27.3533
-3.32166
-0.580486
14
-89.8760
32.9555
-14.5703
2.33812
20.4750
4.03684
15
-64.0213
8.74729
-18.2781
11.1739
4.56532
4.64756
16
24.1308
-4.13785
-10.5377
8.70718
8.82819
0.635262
17
7.18290
-15.8651
8.20768
-3.85646
7.90065
-5.18798
18
6.47624
-14.2352
3.61354
-9.05050
8.53147
-5.94042
19
11.5967
-4.96041
-1.97678
-11.5722
5.81598
-2.55565
20
13.6620
-8.57086
-3.04740
-12.5938
5.49477
-0.341653
21
9.67771
-13.5638
-3.17290
-11.0283
3.86244
2.55857
22
8.28646
-17.2128
1.59427
-6.97636
3.11115
2.62101
23
3.85626
1.39991
11.8803
0.615213
-1.67941
3.97740
24
11.9982
-10.1443
8.96102
4.82252
-3.59272
6.98522
25
14.4479
-15.7646
12.6704
5.13767
-0.285371
7.41964
26
14.5485
4.14165
10.7802
-1.09407
2.70395
8.33036
27
19.7817
15.1572
24.4290
9.36452
4.34313
-0.319466
28
19.8501
26.3028
27.7652
14.1172
2.67827
-0.136420
29
17.8695
5.40787
15.9846
1.12875
7.99281
1.50315
30
15.8722
5.87849
9.71566
-4.84532
5.26982
5.52796
31
6.55574
-19.8384
0.394553
0.291731
-6.93900
16.5671
The research studies of daily power load curves obtained by the AMR system in more than 100
points of head line sections in 6–110 kV distribution networks for 13 days identified a strong statistical
relationship between the first GLDs, i. e. proximity of variances of the power load curves and close
correlation dependence [11]. This allows us to conclude that the factor model of the covariance matrix
and power load curves have statistical stability, and it is possible to apply it for the simulation of power
consumption irregularity for the posterior or previous analogous time periods.
The results of multiple computations indicate a rather large contribution of the three principal
components (70–80%) to the total variance of the whole sample of original power load curves. The first
# 50 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Ivan V. Shulgin, Aleksey A. Gerasimenko... Stochastic Simulation of Covariance Matrix and Power Load Curves...
principal components of GLDs indicate the presence of common internal reasons for daily irregularity
in power load curves and intersystem power flows [6].
The collective stability of the factor model allows us to suggest that GLDs reflect the main reasons
for power changing without specific factors. Therefore, in order to determine the GLDs, we do not need
to analyze the power load curves in all nodes; it is enough to take only some modeling subset of power
load curves into account, for instance, combined diagrams of consumer groups. These are formed on
the basis of check measurements of power consumption, which is carried out by an inspectorate. The
computation results showed that sets of GLDs which were obtained on the basis of analysis of large
samples of power load curves are quite close to each other [6].
Computations for different samples of daily and monthly power load curves have also confirmed
the statistical stability of the factor model (see above). The contribution of the first principal component
to the total variance of power loads is over 50%, while the significant contribution of the first three
components was also confirmed (see Table 1).
Wide application of factor analysis for power load simulation offers a possibility to limit the
volume of the statistical sample to not more than 100 elements [7, 8], which allows to manipulate not
large covariance matrices. In this case, the requirements of representativeness are carried out, and the
obtained GLDs are statistically stable diagrams.
A Number of Principal Components
A recent review of the approaches to determine the required number M of eigenvalues and
eigenvectors for covariance matrix simulation identified the absence of a consensus among specialists
on the factor analysis [11]. One opinion is that this number is usually not higher than four. However, in
some cases depending on the accuracy of the covariance matrix simulation, it is required to take into
account a greater number of characteristic values and GLDs. Researches of power covariance matrix
and power load curves showed that the required number of characteristic values depends on sample
data properties and on the irregularity of power load curves. The maximal number is equal to the first
six maximal characteristic values in MATLAB system. Most likely, this number was chosen after
extensive analysis of covariance matrices and factor simulation.
A simple approach for the determination of a reasonable number of characteristic values is to
estimate the overall contribution of the principal components’ sequence Ƚ 1 , Ƚ 2 ,..., Ƚ k to the total
variance. If the summarized contribution to the total variance is 75–90%, we should stop on the value
k=1, 2,…, М [5, 6]. The total percentage contribution Θ to the variance for fixed M is calculated by the
following formula
M
Ĭ
¦O
k
k 1
2n
¦O
˜ 100% ; 75 d Ĭ d 90% ,
(19)
i
i 1
where
2n
¦O
i
– the sum of eigenvalues of original power covariance matrix, which is called spur of
i 1
matrix [7, 8]. Θ is a criterion for the accuracy of the covariance matrix simulation and the simulation of
original active and reactive power load curves. This criterion is sufficient to carry out the computations
of the power integral characteristics and power system optimization.
# 51 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Ivan V. Shulgin, Aleksey A. Gerasimenko... Stochastic Simulation of Covariance Matrix and Power Load Curves...
Component analysis is used to simulate both random and deterministic dependencies similar to
the regression method of approximating the functions dependent on the time. The most effective test
of statistical hypothesis and the determination of the required number of principal components is the
repeated application of component analysis for various samples of the same general population. If the
statistical hypothesis about the existence of common dominant trends in all random variables Хi is
true, then the statistical characteristics, at least for the first principal component selected on the basis
of different sample data, will be close to each other.
The statistical characteristics for principal components Гk not characterized by the properties
of one general population differ substantially from sample to sample. Repeated factor simulation for
different sample data is a universal method to select statistically stable principal components which
characterize the properties of the general population.
The factor simulation of a set of random variables is a useful instrument of statistical analysis
if the dimension of the model space of M initial random variables is a sufficiently small one. Such a
situation is typical for the simulation of nodal power loads. The application of factor analysis methods
allows us to simulate hundreds of power load curves by the means of 2–3 GLDs [7, 8].
For a reliable application of the method, it is necessary that the first eigenvalues of the K are
significantly different from each other. The matrix K, which corresponds to the nodal power load
curves, usually satisfies this condition [12].
The required number of GLDs for the simulation of original power load curves or required number
of characteristic values for the simulation of power covariance matrix depends on the properties of
the concrete set of research random variables. The research studies on different collections of power
load curves confirmed the hypothesis about high quality simulation of power covariance matrix. This
simulation is based on a small number of characteristic values reflecting the maximal part of the total
variance for the whole general population; hereby a small number means 2÷5 GLDs, the exact number
depending on the properties of the sample data, desired accuracy and purposes of simulation.
Stochastic Simulation of Covariance Matrix
and Power Load Curves’ Algorithm
Referring to the above-mentioned points we can formulate a simulation algorithm for covariance
matrix and power load curves:
1. The original power covariance matrix K(P,Q) is formulated on formulas (13), (14) based on
retrospective analysis of active and reactive power curves in distribution networks for a certain time
period T and control measuring data.
2. 2n eigenvalues and 2n eigenvectors of the K(P,Q) are determined by means of any suitable
software based on the principal component method. The characteristic values identify a degree of
statistical relationship between random deviations of powers Pi, Pj, Qi, Q j from their mathematical
expectations MPi, MPj, MQi, MQ j.
3. A selection of first M maximal eigenvalues of the K(P,Q) λ k < λ and eigenvectors υk < υ is
arranged in descending order. The number of eigenvalues and GLDs needed for covariance matrix
simulation is an acceptable one for practical computations of integral characteristics and system
optimization if condition (19) is satisfied. In other words, we should stop on such value of M where
the contribution of sum of the first diagonal elements of the covariance matrix to the total variance
# 52 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Ivan V. Shulgin, Aleksey A. Gerasimenko... Stochastic Simulation of Covariance Matrix and Power Load Curves...
of the general population of active and reactive powers is 75–90%. The recommended range of M is
2 ≤ M ≤ 5.
4. The simulation of the power covariance matrix K(P,Q) is carried out on formula (9) by the
means of M characteristic values. Computations like that can be performed only once a month based
on the analysis of a training sample of power load curves.
5. The GLDs are determined on formula (16) for the whole collection of the corresponding sample
data.
6. The first M maximal GLDs Гk are selected on (17), which correspond to the M eigenvectors
υk and eigenvalues λk of K(P,Q). The original active and reactive power curves are simulated on the
expression (18) by the means of Гk.
This method is effective when the condition М <<2n can be limited to take into account only the
first few υk and Гk. In addition, the GLDs Гk obtained for the different random process realizations
of power load changing must have statistical stability. The property of universality was confirmed by
means of computations of the GLDs Гk of statistically representative sample data for different power
utilities. In all cases, 3–4 Гk are usually sufficient for the reflection of up to 75–95 % of the total
variance of original power loads (see tab. 1–3).
Thus, the original power load curves can be presented in the form of some characteristics: the
mathematical expectations and the coefficients λk and υk. This is used for effective determination of
integral characteristics and optimization algorithm in power distribution systems.
In order to obtain a more accurate model of covariance matrix and power load curves it is necessary
to take into account a larger number of characteristic values and GLDs, and it is possible to use the
expansion in a Fourier series or similar methods [8]. However, the component analysis differs from
other statistical methods by offering a more economical and convenient way for system optimization
and a good form for the presentation of information. One advantage of the component analysis is
that it is determined in such way that the function package of simulation is not selected randomly, as
for instance, in Fourier analysis, but on the basis of analyzing the principal regularity of power load
changing. The indicated regularities obtained on the basis of the main factors, their number is less than
in other simulation methods, determine the possibility of their application for small samples containing
4–6 points of a daily curve.
The methodology of the integral characteristic determination is developed on the basis of
stochastic simulation of covariance matrix and power load curves particularly the load-dependent
energy loss determination [13, 14] and optimization algorithms of operating conditions on reactive
power in electric distribution networks [15, 16].
Conclusions
1. This paper proposes the factor model of power loads on daily and monthly time periods, which
allows us to identify the general regularity of nodal power changing in electric distribution networks.
The advantages and possibilities of the model application are described.
2. The computational results of the general load diagrams are derived for different sample data of
original active and reactive power curves for daily and monthly time periods T. The statistical stability
of factor power load model is confirmed for 6–110 kV electric power networks. The contribution of the
first principal component to the total variance of power loads is more than 50%.
# 53 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Ivan V. Shulgin, Aleksey A. Gerasimenko... Stochastic Simulation of Covariance Matrix and Power Load Curves...
3. The simulation algorithm for covariance matrix and nodal power load curves is formulated by
means of only a subset of the first main factors (2 ≤ М ≤ 5); it is possible to decrease the complexity
of computations in comparison to traditional computation of multiple operating conditions and to
simplify the determination of power integral characteristics.
4. The algorithms of energy loss determination and ranging of reactive power changing are
developed on the basis of the proposed stochastic model of covariance matrix and power load curves
[13–16]. This can be effectively used for solving the problems of energy saving and power system
optimization.
References
[1] Надтока И.И., Седов А.В., Холодков В.П. Применение методов компонентного анализа
для моделирования и классификации графиков электрической нагрузки // Известия высших
учебных заведений. Электромеханика. 1993. № 6. С. 21–29.
I. I. Nadtoka, A. V. Sedov, V. P. Holodkov. Application of Component Analysis Methods for
Simulation and Classification of Electric Power Loads// Proceedings of the Higher Education
Institutions “Elecromechanics”, 1993, #6, pp. 21–29.
[2] Haibin Wang, Noel N. Schulz. Using AMR data for load estimation for distribution system
analysis// Electric Power Systems Research, www.elsivier.com/locate/epsr. 76 (2006). Pp. 336–342.
[3] Juan Jose Gonzalez de la Rosa, Antonio Moreno Munoz. Aurora Gil de Castro, Victor Pallares
Lopez, Jose Antonio Sanchez Castillejo. A Web-based distributed measurement system for electrical
power quality assessment. www.elsevier.com/locate/measurement, Measurement 43 (2010), pp. 771–
780.
[4] Журавлев В., Грицай М., Артамонов И. Распределительные сети нового поколения//
Электрика. №2, 2007. С. 13–15.
V. Zhuravlev, M. Gricai, I. Artamonov, Distribution Networks of New Generation. #2, 2007. – pp.
13–15.
[5] Герасименко А.А., Липес А.В. Применение компонентного анализа для
определения интегральных характеристик электрических систем// Оптимизация режимов
электропотребления промышленных предприятий и районов: Межвузовский сборник.
Красноярск. КПИ, 1982. С. 101–110.
A.A. Gerasimenko, A.V. Lipes, Application of Component Analysis for Determination of Power
Integral Characteristics in Industry. Interuniversity Collection. Krasnoyarsk. KPU, 1982. – pp. 101–
110.
[6] Арзамасцев Д.А., Липес А.В., Ухалов В.А. Алгоритм статистического определения
интегральных характеристик установившихся режимов электроэнергетических систем/
Известия Академии наук СССР. Энергетика и транспорт. 1984. – С. 39–48.
D.A. Arzamastsev, A.V. Lipes, V.A. Uhalov, Algorithm of Statistical Determination of Power
Integral Characteristics in Power Engineering Systems/ Proceedings of the Academy of Sciences
“Energy and transport”, 1984. Pp. 39–48.
[7] Липес А.В. Применение методов математической статистики для решения
электроэнергетических задач: учеб. пособие. Свердловск, изд. УПИ им. С. М. Кирова, 1983,
88 с.
# 54 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Ivan V. Shulgin, Aleksey A. Gerasimenko... Stochastic Simulation of Covariance Matrix and Power Load Curves...
A.V. Lipes, Application of Mathematical Methods of Statistics for Solving Power Engineering
Problems/ Tutorial book. Sverdlovsk, Edition of Ural Polytechnic Institute, C. M. Kirov, 1983,
88 p.
[8] Арзамасцев Д.А., Липес А.В. Снижение технологического расхода электроэнергии. – М.:
Высшая школа, 1989. – 127 с.
D.A. Arzamastsev, A.V. Lipes, Reduction of Technological Energy Losses: Tutorial, Moscow,
“Higher School”, 1989, 127 p.
[9] D.N. Lawley, A.E. Maxwell. Factor Analysis as a Statistic Method, “World”, 1967. 144 p.
[10] Герасименко А.А., Тихонович А.В. Факторное моделирование нагрузок
распределительных сетей электроэнергетических систем/ Вестник Ассоциации выпускников
КГТУ. 2005. Выпуск 12. С. 147–156.
A.A. Gerasimenko, A.V. Tihonovich, Factor Simulation of Power Loads in Electric Distribution
Systems// Herald of the Graduate Student Association, Krasnoyarsk State Technical University, 2005,
#12, pp. 147–156.
[11] Тихонович А.В. Расчёт потерь электроэнергии в распределительных электрических
сетях на основе объединения детерминированного и стохастического методов и алгоритмов:
автореф. дис. … канд. техн. наук. Красноярск, 2008. 20 с.
A.V. Tihonovich, Calculation of Energy Losses Based on Deterministic and Stochastic Methods
and Algorithms in Electric Distribution Systems/ PhD Thesis Summary. Krasnoyarsk, 2008, 20 p.
[12] Герасименко А.А. Применение ЭЦВМ в электроэнергетических расчетах: учеб.
пособие. – Красноярск: Изд. КПИ, 1983. – 116 с.
A.A. Gerasimenko, EDC Application in Power Engineering Computations// Tutorial Book,
Krasnoyarsk, 1983, 116 p.
[13] Герасименко А.А., Нешатаев В.Б., Шульгин И.В. Расчёт потерь электроэнергии
в распределительных электрических сетях на основе вероятностно-статистического
моделирования нагрузок: Известия высших учебных заведений. Электромеханика, 2011, № 1,
С. 71–77.
A.A. Gerasimenko, V.B. Neshataev, I.V. Shulgin, Energy Loss Evaluation Based on Power Load
Simulation in Electric Distribution Networks/ Proceedings of the Higher Education Institutions
“Electromechanics”, 2011, #1, pp. 71–77.
[14] Герасименко А.А., Нешатаев В.Б., Шульгин И.В. Вероятностно-статистическое
определение потерь электроэнергии в задаче оптимальной компенсации реактивной мощности
в распределительных сетях // Энергетика в современном мире: материалы IV Всероссийской
научно-практической конф. Ч. 1. Чита: ЧитГУ, 2009. С. 214–221.
A.A. Gerasimenko, V.B. Neshataev, I.V. Shulgin, Stochastic Determination of Energy Losses in
Optimal Reactive Power Compensation Problem in Distribution Networks/ Power Engineering in the
Modern World/ Proceedings of IV All-Russian Scientific Practical Conference. Vol. 1. Chita: Chita
State University, 2009, pp. 214–221.
[15] Герасименко А.А., Липес А.В. Оптимизация режимов электрических систем на основе
метода приведенного градиента // Электричество, 1989, № 9. С. 1–7.
A.A. Gerasimenko, A.V. Lipes, Optimization of Power Operating Conditions Based on Reduced
Gradient Method// “Electricity” Journal, 1989, # 9, pp. 1–7.
# 55 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Ivan V. Shulgin, Aleksey A. Gerasimenko... Stochastic Simulation of Covariance Matrix and Power Load Curves...
[16] Герасименко А.А., Нешатаев В.Б., Шульгин И.В. Оптимальная компенсация реактивных
нагрузок в системах распределения электрической энергии// Известия высших учебных
заведений. Проблемы энергетики. 2008. №11–12/1. С. 81–88.
A.A. Gerasimenko, V.B. Neshataev, I.V. Shulgin, Optimal Compensation of Reactive Power
Loads in Energy Distribution Systems// Proceedings of the Higher Education Institutions. “Power
Engineering Problems” Journal # 11-12/1, 2008, pp. 81–88.
Стохастическое моделирование матрицы
корреляционных моментов и графиков нагрузок
мощностей узлов
распределительных электрических сетей
И.В. Шульгина,
А.А. Герасименкоа, Су-Чуан Джоуб
a
Сибирский федеральный университет
Россия 660041, Красноярск, пр. Свободный,79
б
Харбинский политехнический университет
КНР, Харбин
Разработан алгоритм стохастического моделирования матрицы корреляционных моментов и
графиков активной и реактивной мощностей узлов распределительных электрических сетей
на основе факторного анализа. Обоснована статистическая устойчивость факторной модели
электрических нагрузок. Применение данной модели позволяет выявить общие закономерности
изменения мощностей нагрузочных узлов сети и упростить, сделать эффективными методы
анализа и учёта многорежимности в задачах эксплуатации распределительных электрических
сетей и оптимизации их режимов.
Ключевые слова: вероятностно-статистическое моделирование электрических нагрузок,
матрица корреляционных моментов, распределительная электрическая сеть, метод главных
компонент, потери электроэнергии, оптимизация режимов.
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Journal of Siberian Federal University. Engineering & Technologies 1 (2012 5) 57-62
~~~
УДК 533.528
Исследование кавитационных течений
средствами математического моделирования
В.А. Кулагин*, Т.А. Пьяных
Сибирский федеральный университет,
Россия 660041, Красноярск, пр. Свободный, 79 1
Received 6.02.2012, received in revised form 13.02.2012, accepted 20.02.2012
Выполнен обзор современных методов математического моделирования кавитационных
течений, а также анализ их возможностей и недостатков.
Ключевые слова: кавитационные течения, CFD, математическая модель, баротропная модель,
бароклинная модель.
С расширением области приложения кавитационной технологии в различных отраслях
производства возрастает потребность в исследовании кавитационных течений. Однако одновременное существование граничной динамики, фазового перехода и сильного изменения
плотности значительно затрудняет эту задачу.
В начале 1990-х гг. в результате совершенствования CFD-методов появились кавитационные модели, включающие решение уравнений Навье-Стокса. Выделяют два основных подхода
при моделировании многофазных потоков: с взаимопроникающими средами и без взаимопроникновения.
В подходе, который рассматривает невзаимопроникающие среды, четко определяется поверхность раздела пара и жидкости. Этот подход был разработан для моделирования устойчивых
кавитационных течений, в которых известны начальная форма каверны и модель области замыкания (следа). Уравнения движения решаются только для жидкой фазы, а паровая фаза учитывается граничными условиями на поверхности раздела. Массовый поток через поверхность раздела
фаз не учитывается. Начальную форму каверны и ее область замыкания R. Hirschi, P. Dupont и F.
Avellan [1, 2] определяли с помощью уравнения Рэлея-Плессета, тогда как M. Deshpande, J. Feng
[3] и Y. Chen, S.D. Heister [4] применяли критерий статического давления пара. Однако модели,
предполагающие четкую поверхность раздела фаз, не нашли широкого применения.
Как показал анализ работ по исследованию кавитационных течений, наиболее часто встречается подход, основанный на взаимопроникающих средах. Здесь не предполагается поверхность раздела между двумя несмешивающимися жидкостями, поэтому объемная доля фазы
может изменяться от нуля до единицы в зависимости от занимаемого пространства в двухфазном потоке. Данный подход включает одно- и многожидкостные модели.
*
1
Corresponding author E-mail address: vak-sfu@mail.ru
© Siberian Federal University. All rights reserved
# 57 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
В.А. Кулагин, Т.А. Пьяных. Исследование кавитационных течений средствами математического моделирования
В многожидкостных моделях уравнения сохранения применяются для каждой фазы потока [5-7]. Следовательно, фазы в любой точке расчетной области могут характеризоваться различными скоростями и температурами, т.е. они не находятся в равновесии. Взаимодействие
между фазами учитывается введением источниковых членов в уравнения сохранения. Многожидкостные модели наиболее точно описывают физическую сущность кавитационных течений, но их использование связано со значительными трудностями в определении источниковых членов.
В одножидкостных моделях кавитационный поток рассматривается как гомогенная смесь
[8]. Эти модели основываются на предположении локального кинематического и термодинамического равновесия между фазами. В этом случае уравнения сохранения применяются
для смеси, также необходимо замыкающее уравнение для нахождения объемной доли фазы.
Одножидкостные модели классифицируются по типу уравнений, которые используются для
определения содержания доли фазы в потоке, а именно: алгебраические и дифференциальные.
Алгебраические модели предполагают мгновенное влияние локального давления на плотность гомогенной смеси. Поэтому они также известны как баротропные модели или модели,
основанные на уравнении состояния. Упрощая уравнение энергии, можно получить следующее баротропное уравнение:
dUm dP
,
(1)
dt
dt
где cm – скорость звука в смеси, м/c; ρm = αρv + (1 − α ) ρl – плотность смеси, кг/м3, ρv и ρl –
cm2
плотности пара и жидкости соответственно, кг/м3; α – объемная доля пара; P – давление, Па.
O. Coutier-Delgosha [9] учитывал скорость звука в смеси с помощью полиномиального
уравнения. D.P. Schmidt и др. [10] использовали переменную скорость звука, основываясь на
классической гомогенной модели равновесия Wallis:
1
cm2
wȡ m
wP
§ Į
1 Į ·
ȡm ¨ 2 .
2 ¸
© ȡ v cv ȡ lcl ¹
(2)
В результате интегрирования уравнения (2) можно получить зависимость между давлением и плотностью, которая дает начало алгебраическим моделям кавитации [11-13]:
P ȡm Pɧ ª§ ȡ c ·
ȡv ȡl
1
log «¨ m m ¸
2
2
ȡ vȡ l ȡ lcl ȡ v cv «¬© ȡlcl ¹
2
º
».
»¼
(3)
Однако в баротропной модели плотность зависит только от давления, вследствие чего не
учитываются некоторые особенности кавитационных течений. Это может быть легко замечено
при записи уравнения переноса завихренности:
wȦ
1
u ˜ ’Ȧ Ȧ ˜ ’u 2 ’Um u ’P viscous terms,
wt
ȡm
(4)
где u – абсолютная скорость, м/с; Ȧ ’ u u – завихренность, 1/с.
Второй член с правой стороны равенства при ȡ m ȡ m P обращается в ноль. В результате
исчезает значительный источник турбулентности, который является существенным в области
закрытия каверны [14].
# 58 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
В.А. Кулагин, Т.А. Пьяных. Исследование кавитационных течений средствами математического моделирования
Для более полного описания природы гидродинамической кавитации были разработаны
дифференциальные модели. К ним относятся модели, основанные на уравнении переноса объемной доли паровой фазы, а также модели, включающие уравнение Рэлея-Плессета.
Наиболее часто применяется уравнение переноса объемной доли паровой фазы, имеющее
следующий вид:
w
(5)
ȡ v Į ’ ȡ v Į u SĮ ,
wt
где SĮ – источниковый член, который зависит от свойств жидкости и процесса фазового
перехода (испарение или конденсация).
Основная трудность данного метода заключается в определении SĮ. A. Alajbegovic и др.
[15], а также W. Yuan и др. [16] предложили кавитационную модель, базирующуюся на уравнении Рэлея. C.L. Merkle [17], R. Kunz [18], V. Ahuja [19], N. Singhal [20] применяли полуаналитические уравнения для массового переноса. В качестве примера приведем выражение для
определения SĮ, полученное N. Singhal [20]. В данном случае источниковый член складывается из двух составляющих:
1/2
Re
Ce
§2 P P·
Vch
ȡ lȡ v ¨ ˜ ɧ
¸
ı
© 3 ȡv ¹
Rc
Cc
§ 2 P Pɧ ·
Vch
ȡ lȡ l ¨ ˜
¸ Į,
ı
ȡl ¹
©3
1 Į ,
ɩɪɢ P d Pɧ ;
1/2
ɩɪɢ P t Pɧ ,
(6)
(7)
где Vch – характеристическая скорость, м/с; Pн – давление насыщения, Па; σ – поверхностное
натяжение воды, Н/м; Ce = 0,02 и Cc = 0,01 – эмпирические константы.
I. Senocak и W. Shyy [14] получили полностью аналитическую кавитационную модель
переноса, основанную на уравнениях сохранения массы и импульса для поверхности раздела
пар – жидкость.
Модели, предложенные A. Alajbegovic и др. [15], а также W. Yuan и др. [16], описывают
предельный случай роста сферического пузырька при изменении давления окружающей жидкости. Однако данные модели плохо адаптированы для учета коллапса пузырька и пренебрегают многими эффектами, которые определяют поведение пузырьков. Влияние инерции, вязкости, поверхностного натяжения жидкости и сжимаемости парогазовой составляющей смеси
на структуру кавитационного потока моделируется более универсальным уравнением РэлеяПлессета [21]:
§
3 R 2 ·
ȡ¨ R u R
¸
2 ¹
©
Pv Pg Pl 4ȝ
R
ı
2 ,
R
R
(8)
где R – радиус кавитационного пузырька, м; Pv , Pg , Pl – давление пара, газа и жидкости
в кавитационной области соответственно, Па; μ – коэффициент динамической вязкости,
Па·с.
В течение долгого времени уравнение Рэлея-Плессета (8) решалось в одномерном виде.
Однако в связи с повышением мощности ЭВМ появились работы, отражающие расчетный анализ двух- и трехмерных нестационарных кавитационных потоков [22-25].
# 59 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
В.А. Кулагин, Т.А. Пьяных. Исследование кавитационных течений средствами математического моделирования
В настоящее время стало возможным исследование кавитации при помощи «тяжелых» в
вычислительном смысле моделей E. Giannadakis, M. Gavaises, H. Roth, C. Arcoumanis [25] и G.L.
Chahine [24], которые включают уравнение для динамики группы пузырей радиуса R и плотности их числа n, полученное A. Kubota и др. [22]:
3Ȗ
2ı 4ȝ R
Pv P Pg ,0 § R0 ·
¨ ¸ ȡ
ȡ © R¹
ȡR ȡ R
nR
2 R 2nRR 2 ,
2ʌ ˜ 'r 2 ˜ nR 2 R
3 R 2
Ru R
2
(9)
где Pg ,0 и R0 – первоначальные давление газа и радиус пузырька, γ – показатель политропы.
Последний член в этом уравнении моделирует динамику пузырьков на подсеточном масштабе
Δr. Эта особенность становится полезной с использованием грубых сеток, где давление
распределяется между несколькими пузырями в вычислительной ячейке.
Сравнение математических моделей (в таблице приведены свойства трех широко применяемых типов одножидкостных моделей) кавитационных потоков показывает, что наиболее
распространенными для инженерных расчетов являются одножидкостные модели.
Свойства математических моделей кавитации
Дифференциальные модели
Равновесие фаз
Hеравновесное
состояние фаз;
постоянная скорость
фазового перехода
Модели, основанные
на уравнении РэлеяПлессета
Неравновесное
состояние фаз;
переменная скорость
фазового перехода
Баротропная модель
Бароклинная модель
Бароклинная модель
Не учитываются
Могут быть учтены
Могут быть учтены
Пренебрегается или
учитывается для смеси
Учитывается для смеси
Учитывается для смеси
или жидкости
Алгебраические модели Модели, основанные на
уравнении переноса
Фазовый
переход
Зависимость плотности
от давления
Эффекты
кавитационных ядер
Турбулентность
Вывод
Алгебраические модели предполагают явное влияние давления на объемную долю паровой фазы, тем самым представляя природу потока баротропной. Модели, которые используют дифференциальные уравнения, описывают бароклинную природу кавитационного потока,
более соответствующую реальности. Модели кавитации, основанные на уравнении переноса
объемной доли паровой фазы, являются устойчивыми и популярными в CFD относительно
моделей, основанных на уравнении Рэлея-Плессета. Однако последние более пригодны для
моделирования нестационарных кавитационных течений, хотя и требуют больших вычислительных ресурсов.
# 60 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
В.А. Кулагин, Т.А. Пьяных. Исследование кавитационных течений средствами математического моделирования
Список литературы
[1] Hirschi R. Prediction par Modelisation Numerique Tridimensionnelle des Effets de la Cavitation a Poche dans les Turbomachines Hydrauliques. PhD thesis, Ecole Polytechnique Federale de
Lausanne. 1998. №1777.
[2] Hirschi R., Dupont P., Avellan F. Centrifugal pump performance drop due to leading edge
cavitation: Numerical predictions compared with model tests // Trans. of ASME, J. Fluids Eng., december 1997. №120. Р. 705-711.
[3] Deshpande M., Feng J. Cavity flow predicitions based on the euler equations // Trans. of
ASME, J. Fluids Eng., march 1994. №116. Р. 36-44.
[4] Chen Y., Heister S.D. A numerical treatment for attached cavitation // Trans. of ASME, J. Fluids Eng., september 1994. №116. Р. 613-618.
[5] Von Berg E., Edelbauer W., Alajbegovic A., et al. Coupled simulations of nozzle flow, primary
fuel jet breakup, and spray formation // J. Eng for Gas Turbines Power, 2005. №127(4). Р. 897-908.
[6] Chiavola O., Palmieri F. Modeling needle motion influence on nozzle flow in high pressure
injection system // SAE Paper, 2007. 2007-01-0250.
[7] Wang X., W.H. Su. A numerical study of cavitating flows in high-pressure diesel injection
nozzle holes using a two-fluid model // Chinese Science Bulletin, 2009. №54(10). Р. 1655-1662.
[8] Zhang X.B, Qiu L.M., Gao Y., Zhang X.J. Computational fluid dynamic study on cavitation in
liquid nitrogen // J. Cryogenics, 2008. №48 (9-10). Р. 432-448.
[9] Coutier-Delgosha O. Modelisation des Ecoulement Cavitants: Etude des Comportements Instationnaires et Application Tridimensionnelle aux Turbomachines. PhD thesis, LEGI-INPG, Grenoble,
France, Nov. 2001. №5519.
[10] Schmidt D.P., Rutland C.J., Corradini M.L. A numerical study of cavitating flow through
various nozzle shapes // SAE paper, 1997. 971597. Р. 10.
[11] Schmidt D.P., Corradini M.L. The internal flow of diesel fuel injector nozzles: a review // Int.
J. Engine Research, 2001. Vol. 2. №1. Р. 1-22.
[12] Qin J.R. Direct Calculations of Cavitating Flows by the Space-Time CE/SE Method // AIAA
conf. Paper, 2001. P. 324.
[13] Dumont N., Simonin O., Habchi C. Numerical simulation of cavitating flows in Diesel injectors by a homogeneous equilibrium modelling approach // 4th Int. Symposium on Cavitation: CAV2001,
2001. P. 44.
[14] Senocak I., Shyy W. Evaluation of cavitation models for Navier-Stokes computations // Proceedings of the 2002 ASME Fluids Engineering Division Summer Meeting, 2002. №31011.
[15] Alajbegovic A., Grogger H.A., Philipp H. Calculation of transient cavitation in nozzle using
the two-fluid model // Proc. ILASS-Americas’99 Annual Conf., 1999. P. 373-377.
[16] Yuan W., Schnerr G.H. Numerical simulation of two-phase flow in injection nozzles: interaction of cavitation and external jet formation // J. Fluids Eng, 2003. №125(6). P. 963-969.
[17] Merkle C.L., Feng J., Buelow P.E. Computationalmodeling of the dynamics of sheet cavitation //
Third international symposium on cavitation (CAV1998), 1998. P. 307-311.
[18] A preconditioned navier-stokes method for two-phase flows with application to cavitation prediction / R.F. Kunz, D.A. Boger, D.A. Stinebring, T.S. Chyczewski, H.J. Gibeling, S. Venkateswaran,
T.R. Govindan // Computers & Fluids, 2000. №29. P. 849-875.
# 61 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
В.А. Кулагин, Т.А. Пьяных. Исследование кавитационных течений средствами математического моделирования
[19] Ahuja V. Hosangadi A., Arunajatesan S. Simulations of Cavitating Flows Using Hybrid Unstructured Meshes // Trans. of ASME, J. Fluids Eng, June 2001. №123. P. 331-340.
[20] Singhal N.H., Athavale A.K., Li M., Jiang Y. Mathematical basis and validation of the full
cavitation model // Tr. ASME, J. of Fluids Engineering, 2002. Vol.124. P. 617-624.
[21] Brennen C. E. Cavitation and Bubble Dynamics. Oxford University Press, New-York, 1995.
P. 294.
[22] Kubota A., Kato H., Yamaguchi H. A new modelling of cavitating flows : A numerical study
of unsteady cavitation on a hydrofoil section // J. Fluid Mech, March 1992. №240. P. 59-96.
[23] Delale C.F. Thermal Damping in Cavitating Nozzle Flows // Proc. of the 4th Int. Symposium
on Cavitation: CAV2001, 2001. P. 7.
[24] Chahine G.L. Nuclei Effects on Cavitation Inception and Noise // 25th Symposium on Naval
Hydrodynamics. St. Johns, August 2004. №14. P. 8-13.
[25] Giannadakis E., Gavaises M., Roth H., Arcoumanis C. Cavitation Modelling in SingleHole Injector Based on Eulerian-Lagrangian Approach // Conference on Thermo- and Fluid Dynamic
Processes in Diesel Engines, 2004. P. 14.
Research of Cavitating Flows by Methods
of Mathematical Simulation
Vladimir A. Kulagin and Tatyana A. Pyanykh
Siberian Federal University
79 Svobodny, Krasnoyarsk, 660041 Russia
The purpose of given article – to execute the review of modern methods of mathematical modelling
cavitating flows, and also to realize the analysis of their possibilities and limitation.
Keywords: cavitating flows, CFD, mathematical model, barotropic model, baroclinic model.
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Journal of Siberian Federal University. Engineering & Technologies 1 (2012 5) 63-71
~~~
УДК 666.293.53
Influence of Initial Components
Mechanical Activation on the Properties
of Protective Enamel Coating
on the Cast-Iron Surface
Igor A. Shimanskiyа*,
Vladimir G. Babkinа, Vladimir K. Frizorgerb,
Evgeniy S. Goloskina and Alexey B. Nabiulinа
а
Siberian Federal University,
79 Svobodny, Krasnoyarsk, 660041 Russia
b
«Rus-engineering» Ltd.,
Krasnoyarsk, 660049 Russia 1
Received 6.02.2012, received in revised form 13.02.2012, accepted 20.02.2012
The influence of initial components mixture mechanical activation on the synthesis processes in
enamel mass and on the molten enamel adhesion to the metal was investigated. It’s established, that
the process of albite formation in the coating material intensifies by components mixture activation.
On the base of measured contact angles of cast iron with molten enamel was revealed, that adhesion
of the molten enamel to the metal noticeably increases by initial components activation.
Keywords: cast iron, gas manifold, skirts, high temperature gas corrosion, protective coatings,
mechanical activation.
Introduction
Electrolysis of aluminum from cryolite-alumina melt produces cell gases (СО, СО2, HF, SO2)
and bath evaporation products (NaAlF4). Soderberg cells are equipped with gas collection system the
most important element of which is the gas manifold consisting of cast iron skirts. The skirts are fixed
on the lower part on the anode casing. Together with alumina cover the manifold forms over the cell
perimeter a channel to collect and remove gas. At temperatures up to 1000 К the manifold is subject to
gas corrosion by oxygen, water steam and cell gases. This deteriorates the gas skirts in cell operation.
Disintegration products which are iron compounds get into the bath and then into aluminum. As a
result the grade of aluminum produced from the same raw material in a Soderberg cell is lower (and
accordingly, its price) than the aluminum made in cells with prebaked anode.
The methods to protect cast-iron castings from gas corrosion by applying protective coatings on
its surface, such as painting, tinning, enameling, electrochemical deposition, spraying etc. are known
[1]. General disadvantage of these methods is necessity of metal surface preliminary preparation
*
1
Corresponding author E-mail address: shim_wassup@mail.ru
© Siberian Federal University. All rights reserved
# 63 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Igor A. Shimanskiy, Vladimir G. Babkin… Influence of Initial Components Mechanical Activation on the Properties...
including degreasing, etching, machining etc. that defi nes process multistage and requires special
equipment.
Earlier, in papers [2,3] the composition of an albite-based protective enamel coating for
gas manifold and new one-step method of its application directly in the casting process were
proposed.
The aim of present work is investigation of initial components mixture mechanical activation
influence on the formation process and properties of protective enamel coating on the cast iron
surface.
Research methods
Element and phase composition of materials under study were studied with X-ray analysis. X-ray
patterns were made on automated X-ray equipment XRD-6000 (Shimadzu).
Microstructure of samples was studied with optical microscope Axio Observer А1 (Carl Zeiss)
with magnification from 200 to 500.
Dilatometric measures were carried out with device DIL 402 C (Netszch).
Thermal properties of coating were measured with LFA 457 Micro Flash (Netzsch).
Mechanical activation of materials was carried out with Rocklabs ring mill.
Determination of the contact angle of cast iron with molten enamel was carried out with lying drop
method using vacuum device “Drop” (Giredmet). Preliminary synthesized enamel material granules
were placed on the heated to the desired temperature substrate. After melting and reaching desired
temperature the drop was photographed.
Results
In paper [2] we had shown that inert with respect to aluminum cell aggressive gas environment
material is albite NaAlSi3O8, that confirmed by Gibbs energy change positive value of its reaction with
main corrosion source – hydrogenium fluoride.
NaAlSi3O8 + 4HF = NaF + AlF3 + 3SiO2,
∆Go973 K = +7,5, ∆Go1073 K = +36,46, kJ.
(1)
Albite-based protective coating was made by enamel ceramic mass, its composition is given in
Table 1.
Broken glass and alumina were the sources of sodium, silicon and aluminum oxides required to
form albite. Quartz sand was the filler, argil was the binder. Copper oxide and sodium fluoride acted
as surfactants and flux, respectively.
In paper [3] we had offered one-stage method of applying protective anti-corrosion coating on the
cast-iron casting surface directly in the casting process. The proposed method of forming corrosionresistant coatings on cast-iron castings is to deposit a layer of enamel ceramic mass on the casting mold
surface. Then the casting mold is filled with cast iron melt overheated to the temperature 150-200 К
higher than the melting temperature. Ceramic coating layer is heated due to heat conductivity of the
solidifying casting to the melting temperature and cools to turn into a compact enamel corrosionresistant coating.
# 64 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Igor A. Shimanskiy, Vladimir G. Babkin… Influence of Initial Components Mechanical Activation on the Properties...
Table 1. Composition of enamel ceramic mass for protective coating
Material
Content, mass %
Broken glass
60
Alumina
15
Copper oxide
1
Sodium fluoride
3
Quartz sand
17
Argil
4
In this work the influence of enamel mass initial components activation with a view to intensify
albite synthesis and melting processes and increase coating material adhesion to the metal was
studied.
Thus proceeded with the fact that joint activation of initial components allows to achieve their
even distribution in the mix material volume and leads to mechanical intrusion of one substance to the
surface layer of another with the process of chemical interaction between them.
Activation length was 1÷7 minutes. Upper limit is caused by the fact that this time length provided
maximum dispersion of the system, above which aggregation of particles was observed.
Enamel mass microphotographs in primary and activated conditions are presented on Fig. 1
It’s revealed, that the average particle size in non-activated sample is ~ 90 micron and large
particles the size up to 120 micron are present. After components mixture activation during 7 minutes
the average particle size decreases to ~ 20 micron and the maximum is ~ 40 micron.
Fig. 2 shows X-ray photograms of activated and non-activated enamel mass samples passed
annealing at temperature 1173 K during 40 minutes.
It’s established, that the character of course of phase interactions in investigated system is changed
under activation influence. Qualitative difference in the samples composition was not revealed.
But quantitative composition, presented in the Table 2 considerably differs at the same time.
As it follows from experiments results, albite content in the non-activated sample is ~ 35,6 mass
%. In the activated sample albite content increasing up to 59,7 mass % was revealed.
So, mechanical activation allows to increase yield of the reaction product during the synthesis
under the same conditions significantly, that caused by increasing particles specific surface area and
thus increasing their reactivity.
It confirms by the results of differential thermal analysis, presented in Fig. 3, from which follows
that the value of the thermal effect at a temperature of 1050 K referred to the albite forming reaction
decreases. This fact indicates that the activation energy of chemical transformation in enamel ceramic
mass is decreasing.
Also it’s revealed that annealing of non-activated components mixture leads to porous and
inhomogeneous sample structure. Activation allows to obtain dense and solid samples with more
homogeneous structure. Fig.4 shows samples microstructure after annealing.
Samples’ seeming density was defined by hydrostatic weighing method. Its value is 1,75 and 1,95
g/cm3 for non-activated and activated samples respectively.
# 65 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Igor A. Shimanskiy, Vladimir G. Babkin… Influence of Initial Components Mechanical Activation on the Properties...
ɚ)
ɛ)
Fig. 1. microphotographs of powders, ×250: а) non-activated sample; б) sample activated during 7 minutes
Enamel mass mechanical activation may influence not only on the course of phase transformations
during coating formation, but it also determines the interaction between enamel melt and cast-iron
casting surface, what is possible to judge by adhesion work value.
Adhesion work was calculated with using contact angle measures by Dupre-Young equation:
Wa = σL·(1 + сosθ)
(2)
where σL – enamel melt surface tension, θ – contact angle.
Enamel melt surface tension was determined by Pavlov-Popel equation:
V
V1 2000 lg
¦
k
Fx
(3)
l 1 1 1
# 66 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
1000
950
900
850
800
750
700
650
ɂɧɬɟɧɫɢɜɧɨɫɬɶ
600
550
500
450
400
350
300
250
200
150
100
50
0
6
8
10
12
14
16
18
20
22
24
26
28
30
32
34
36
38
ɍɝɨɥ (2Q)
40
42
44
46
40
42
44
46
48
50
52
54
56
58
60
62
64
66
68
ɚ)
1000
950
900
850
800
750
700
650
600
550
500
450
400
350
300
250
200
150
100
50
0
6
8
10
12
14
16
18
20
22
24
26
28
30
32
34
36
38
ɍɝɨɥ (2Q)
48
50
52
54
56
58
60
62
64
66
68
Fig. 2. Results of X-ray phase analysis: а) non-activated sample; б) sample activated during 7 minutes
Table 2. Quantitative composition of samples
Main phases
Albite
Content, mass %
Non-activated
Activated during 7 minutes
35,6
59,7
Diopside
19,4
16,0
Quartz
15,8
1,84
Nepheline
7,76
17,5
Al2O3
19,3
3,22
CuO
0,73
0,64
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
ɚ)
ɛ)
Fig. 3. Results of differential thermal analysis: а) non-activated sample; б) sample activated during 7 minutes
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Igor A. Shimanskiy, Vladimir G. Babkin… Influence of Initial Components Mechanical Activation on the Properties...
ɚ)
ɛ)
Fig. 4. Samples microstructure: а) non-activated sample; б) sample activated during 7 minutes
where σ1 – solvent surface tension (70 mole % SiO2, 30 mole % Na2O – it approximately matches with
the content of SiO2 and Na2O in broken glass+quartz sand) which is equal to
282; Fi – the parameter
characterizing component’s surface activity; k – quantity of components.
The following enamel composition was accepted in the calculation: (mole %): SiO2 (65), Na2O
(10,9), Al2O3 (10,9); CaO (5,7), MgO (2) and others (5,5). Corresponding values of constants Fi of
enamel melt components are: 1,02; 0,90; 0,65; 0,72; 0,70.
After substitution of σ1 и Fi values equation (3) takes the following form:
σL = 282 – 2000lg (1,02xSiO2 + 0,9xNa2O +
+ 0,65x Al2O3 + 0,72xCaO + 0,70xMgO + ...).
(4)
Value of enamel melt σ calculated by equation (4) is 386 mJ/m 2. Calculation results are presented
in Table 3.
On the base of this data we can conclude that enamel melt forms large contact angles on the castiron. It demonstrates rather large value of interfacial tension σSL because of the weak interaction between
melt and metal. Decreasing of θ was observed at temperatures upper 1373 K for the melt obtained from
# 69 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Igor A. Shimanskiy, Vladimir G. Babkin… Influence of Initial Components Mechanical Activation on the Properties...
Table 3. Mechanical activation influence on the contact angles and adhesion work in the enamel melt-solid cast
iron system
Enamel melt
characteristic
Melt made from nonactivated mix material
Melt made from activated
mix material
Temperature, К
θ, deg
Adhesion work Wa, mJ/m 2
1273
118
204,8
1293
117
210,7
1353
117
210,7
1373
111
247,7
1293
119
198,9
1353
108
266,7
1373
108
266,7
1408
58
590,5
activated mix material. According to this, adhesion work more than 2 times increases. Some (about
5 % on every 100 oC [6]) can’t noticeably affect on Wa increasing rate. Molten enamel transition from
unwettability to wettability can be associated only with σSL reduction, i.e. with interaction with metal
increasing.
Contained in the cast-iron cementite carbon can act like a reducing agent to enamel components
placed on the surface layer (Na2O, SiO2). Sodium oxide reduction reaction can be presented in following
form:
Na2OS + Fe3CS = 2NaG + 3FeS + СОG
(5)
Balance constant of this reaction depends on substances conditions. Thus the essential role is
played by enamel disperse components reducing interaction reaction between melt and metal activation
energy. Iron oxidation is possible in the enamel and metal contact area. Iron oxides transfering to
enamel melt is also promotes the adhesion work growth.
Conclusion
Enamel mass components mechanical activation leads to intensification of their interaction
process with albite formation, to increasing of the enamel melt adhesion to cast-iron casting surface,
and to the coating material porosity reduction and density increasing.
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Igor A. Shimanskiy, Vladimir G. Babkin… Influence of Initial Components Mechanical Activation on the Properties...
Влияние механоактивации исходных компонентов
на свойства защитного эмалевого покрытия
на поверхности чугунных изделий
И.А. Шиманскийа,
В.Г. Бабкина, В.К. Фризоргерб,
Е.С. Голоскинб, А.Б. Набиулина
а
Сибирский федеральный университет,
Россия 660041, Красноярск, пр. Свободный 79
б
ООО «Рус-инжиниринг», Красноярск, Россия,
Исследовано влияние механоактивации смеси исходных компонентов на процессы синтеза
в эмалевой массе, а также на адгезию материала покрытия к металлу. Установлено, что
при активировании исходной смеси интенсифицируется процесс образования альбита в
материале покрытия. На основании измеренных краевых углов смачивания чугуна расплавом
эмали выявлено, что при активировании исходных компонентов адгезия расплава эмали к
металлу заметно увеличивается.
Ключевые слова: чугун, газосборный колокол, секции, высокотемпературная газовая коррозия,
защитные покрытия, механоактивация.
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Journal of Siberian Federal University. Engineering & Technologies 1 (2012 5) 72-86
~~~
УДК 662.613+66-96
Aerodynamic Separation of Fly Ashes
of Selective Sampling from Pulverized Combustion
of Coals of Different Ranks
Olga M. Sharonovaa*, Natalia A. Oreshkinab,
Larisa I. Kurteevaa and Alexander G. Anshitsa,b
a
Institute of Chemistry and Chemical Technology,
Siberian Branch of Russian Academy of Sciences,
42 K.Marx Str., Krasnoyarsk, 660049 Russia
b
Siberian Federal University,
79 Svobodny, Krasnoyarsk, 660041 Russia 1
Received 6.02.2012, received in revised form 13.02.2012, accepted 20.02.2012
This paper reports on the results of the investigation of fly ashes from pulverized combustion of coals
of two ranks, i.e., the B2 (sub C) rank, which is taken from Fields 1−4 of the electrostatic precipitators
at the BSDPS -1 (the Bfa series), and the T (sa) rank, which is taken from Fields 1 and 2 of the
electrostatic precipitators at the MPP-22 (the Mfa series).
It has been shown that fly ashes of the Bfa series have a lower bulk density and a higher dispersity
and a unimodal particle size distribution. The size of particles increases along the gas-and-dust
flow, in particular from the Field 1 to the Field 2 of the electrostatic precipitators. Fly ashes of the
Mfa series are characterized by a higher density, a substantially lower dispersity, a bimodal particle
size distribution and have rather close size characteristics for ashes from the Fields 1 and 2 of the
electrostatic precipitators. The use of the aerodynamic separation has made it possible to obtain
three products from fly ashes of both series which differ significantly in the bulk density and size of
particles. Combination of selective sampling and aerodynamic separation of fly ashes from pulverized
combustion of coals are promising techniques for manufacturing ash products with a controlled
density and dispersity.
Keywords: Fly ash; Aerodynamic separation; Dispersity.
Introduction
In pulverized combustion of coal, the fly ash amounts to 80−95% of the ash and slag wastes of
coal-fired power plants and, consequently, represents the most significant (in volume) by-product of
coal combustion [1−3]. The main directions in the utilization of fly ashes in the world are provided by
industry of building materials, highway engineering, and reclamation of pits [1, 2, 4]. The improvement
of the quality of fly ashes has been achieved, in particular, by means of the design and development of
facilities for conditioning the fly ash with the use of air classification, sieving, drying, and grinding for
the purpose of producing ash fractions with controlled sizes and properties [5, 6]. The size-classified ash
*
1
Corresponding author E-mail address: shar@icct.ru
© Siberian Federal University. All rights reserved
# 72 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Olga M. Sharonova, Natalia A. Oreshkina… Aerodynamic Separation of Fly Ashes of Selective Sampling from Pulverized...
has been employed to manufacture dry binding mixes, cements with improved properties, pigments,
ceramic membrane filters, composite plastic materials, and industrial rubber goods [5−9].
Since coals of different ranks differ in reactivity, their combustion has been carried out using
boilers with different degrees of grinding and at different temperatures [10,11], what is the important
factor affecting the size of the ash particles. The fly ash is partially classified in an ash collection facility
at a coal-fired power plants. For example, ash products of different compositions and dispersities are
formed in different fields of mechanical and electrostatic precipitators [12−14]. For example, Lee et al.
[15] studied aluminosilicate ashes taken from different fields of electrostatic precipitators and showed
that, when changing over from the first field to the third field, the quartz content decreases, whereas the
content of the glass phase increases and that, in the composition of the glass phase, the Al2O3 content
increases, whereas the SiO2 content decreases. Shanthakumar et al. [16] also observed an increase in
the content of the glass phase (from 85.8 to 96.2%) in the fly ashes taken from six sampling points
and revealed a tendency toward a decrease in the ferrospinel content, whereas the Fe2O3 content in the
chemical composition of the fly ash taken from the third field increased with respect to the first and
second fields of the electrostatic precipitators. In our recent study [17], we showed that, when changing
over from the first field to the fourth field of the electrostatic precipitators, the contents of the glass
phase and calcium sulfate in high-calcium fly ashes increase by a factor of 1.4 and ~2.0, respectively,
whereas the contents of calcium aluminate, lime, quartz, and hematite decrease by a factor of 2.4, 2.8,
8.0, and 5.7, respectively.
The most substantial differences between the fly ashes taken from different sampling points are
observed in the particle sizes. When there is a general tendency toward a decrease in the particle size
along the gas-and-dust flow, the dispersities of fly ashes taken from different coal-fired power plants
can differ significantly. This is associated with the grinding fineness of coal, the presence of unburned
carbon−mineral particles, and other factors. For example, according to the data available in the literature
for different coal-fired power plants [14], the content of the fraction <63 μm in the fly ash from the first
field of the electrostatic precipitators amounts to 35−91%, whereas the corresponding content in the
same ashes from the third field of the electrostatic precipitators reaches 60−98%. According to the data
reported by Lee et al. [15], the average size of the ash particles decreases by a factor of 7 (from 28.1 to
3.8 μm) when changing over from the first field to the third field of the electrostatic precipitators.
In the case where the fly ashes are used as additives to cements and concretes, they should
satisfy the requirements for dispersity. In particular, the fly ashes should have a fi neness retained
80 μm sieve ≤15−20% according to State Standard [18]; and a fi neness retained 45 μm sieve ≤34%,
≤40%, ≤12%, and ≤10% according to the foreign standards ASTM 618-2005, BS EN450-1995,
BS3892-1996, and JIS A 6201-1999, respectively [19−22]. In order to satisfy these requirements, it is
common practice to use the mechanical activation of ashes [23], sieving [24], conditioning facilities
[5], and aerodynamic classification [25−27]. At the same time, the selective sampling of the fly ash
at different points of the ash collection facility, in many cases, makes it possible to satisfy these
requirements for dispersity.
Fly ashes of fi ner fractions < 10 or 20 μm (Microsit M10 and M20, respectively) have been used
for manufacturing high-strength concretes, plastic materials, and paints [5]. A fraction < 5 μm has
served as a filler in polymer products or as a component of road pavements with a high strength and
a high water resistance [28] or high-strength multifunctional cements [29]. It should be noted that,
# 73 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Olga M. Sharonova, Natalia A. Oreshkina… Aerodynamic Separation of Fly Ashes of Selective Sampling from Pulverized...
in practical applications of fly ashes, the positive effect (improvement of plasticity, regularity of
the structure, etc.) has been achieved due to a microspherical shape of ash particles. Therefore, the
grinding of fly ash is undesirable. A combination of selective sampling with aerodynamic separation
is a promising technique for producing ash products with a controlled dispersity without grinding
of the ash.
The purpose of this work was to study of the fly ashes sampled from different fields of electrostatic
precipitators (EPs) at two coal-fired power plants, which have burnt coals of different ranks, and the
products of aerodynamic separation of these ashes.
Objects of investigation and experimental techniques
The objects of our investigation are fly ashes from the pulverized combustion of brown coal from
the Berezovsky surface mine in the Kansk−Achinsk Basin at the Power Plant BSDPS-1 and those of
black coal from the Kuznetsk Basin at the Moscow Electric Power Plant No. 22 (MPP-22). The fly
ashes series are designated as the Bfa and Mfa, accordingly. The coal characteristics, boiler types,
coal combustion conditions, and sampling points of the fly ashes are presented in Table 1 [10, 30,31].
As follows from this table, the BSDPS-1 uses a low-rank coal, i.e., coal of the B2 (sub C) rank with
high moisture, high volatile matter content, low ash content, and low heat of combustion. This coal
is combusted in a P-67 boiler with the dry-ash removal at a temperature of 1350−1450°C. In contrast
to the BSDPS-1, the MPP-22 uses a high-rank coal, i.e., coal of the T (sa) rank with low moisture,
low volatile matter content, high ash content, and high heat of combustion. This coal is combusted
in a TPP-210A boiler with the slag-tap removal at a temperature of ≥1600°C. The difference in the
Table 1. Coal characteristics, boiler types, fly ash yields, and sampling points for the fly ashes of Bfa and Mfa
series [4,30,31]. (Wtr - the total moisture content of the fuel under operating conditions, Qidaf – the heat of combustion and Vdaf – the volatile substances in dry ash-free basis, Ad – ash content of fuel in the dry basis)
No.
Characteristic
1.
Power Plant
2.
Coal
(deposit)
Coal rank according to:
GOST 25543-88 [32]
ASTM D388-98a [33]
3.
4.
Wtr (%)
5.
Qi
6.
Vdaf (%)
7.
Ad (%)
8.
Slag removal type
daf
(MJ/kg)
9.
Boiler type
10.
Furnace temperature ( С)
11.
Fly ash yields (%)
12.
Sampling points (fields of
electrostatic precipitators)
o
Fly ash series
Bfa
Mfa
BSDPS-1
MPP-22
Berezovsky coal mine (KanskAchinsk Basin)
Kuznetsk Basin
B2
sub C
Т
sa
33.0
7.0
25.92
34.0
48.0
12.5
7
20
Dry-ash removal
Slag-tap removal
P-67
TPP-210А
1350-1450
>1600
≥95
80-90
1, 2, 3 and 4
1 and 2
# 74 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Olga M. Sharonova, Natalia A. Oreshkina… Aerodynamic Separation of Fly Ashes of Selective Sampling from Pulverized...
Table 2. Yields and characteristics of size distribution for initial fly ashes of the Bfa and MS series and the
products of their aerodynamic separation
Fly ashes and products
of separation
Bfa, Field 1
Yield
(%)
-
ρ bulk,
(g/cm3)
0.93
Average diameter
(μm)
9.3
d90
(μm)
41.8
Maximum
(μm)
21.3
- heavy
29 - 33
1.13
26.2
62.2
37.9
- medium
52 - 56
0.79
7.5
26.5
13.6
- light
10 - 12
0.58
2.8
10.2
4.3
-
0.69
3.1
11.8
5.3
0.67
2.8
9.1
4.0
Bfa, Field 2
Bfa, Field 3
Bfa, Field 4
-
0.63
2.4
7.9
3.4
- heavy
53 - 62
0.72
4.8
13.4
8.9
- medium
23 - 31
0.59
4.7
13.0
8.7
9
0.35
2.9
7.6
4.6
- light
Mfa, Field 1
-
1.31
16.0
73.6
58.7
- heavy
45 - 47
1.05
51.0
97.8
60.0
- medium
35 - 40
1.16
9.6
31.3
16.8
8-9
0.73
4.8
13.9
9.5
-
1.18
15.9
69.9
68.8
- light
Mfa, Field 2
temperature regimes is associated with the high reactivity of the B2 coal and the low reactivity of the
T coal in the course of combustion [10].
The aerodynamic separation was carried out on a laboratory setup, which was described in our
previous paper [27]. An fly ash sample 10 cm3 in volume was charged into the setup in one cycle of its
operation. The heavy product remained in the wind tunnel, the medium product accumulated in the
lower part of the cyclone, and the light product collected on the filter installed at the exit of the cyclone.
The fly ash was separated using the air flow provided by a compressor at linear velocities of 0.051,
0.060, and 0.068 m/s. Table 2 presents the relevant data obtained for the velocity of 0.051 m/s. As the
velocity increases, the separation parameters change toward an increase in the yield of the light and
medium products with increasing fraction of larger particles. The number of cycles was varied from
one to three, and the yield of products was determined in each cycle and was equal to 87-99 wt.%.
Then, the identical products were mixed and averaged, the samples were taken for analysis according
to the State Standards [32,33].
The granulometric compositions of the fly ashes from the first fields of EPs were determined
using dry sieving according to the State Standard [34] on a VP-S/220 vibrodrive equipped with a set
of sieves with meshes 0.400, 0.200, 0.160, 0.100, and 0.063 mm in size. The compositions of the major
components of the ashes and the losses on ignition (LOI) were determined by the chemical analysis
according to the State Standard [35], and the bulk density was determined by the method of measuring
containers. The particle size distribution was investigated on a Fritsch MicroTech 22 laser analyzer. In
this case, samples of Bfa and Mfa series were studied using dry and wet measuring cells, respectively.
The dispersity and morphology of the ash particles were examined on a Carl Zeiss Axioscop-40 optical
microscope equipped with a color digital video camera (Canon).
# 75 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Olga M. Sharonova, Natalia A. Oreshkina… Aerodynamic Separation of Fly Ashes of Selective Sampling from Pulverized...
Results and discussion
Fly ashes from different fields of electrostatic precipitators
The fly ashes studied in this work were obtained from the combustion of coals of different ranks,
which differ in the reactivity, temperature conditions of combustion in different boilers, and ash content.
In particular, they differ in the content of unburned carbon, which is lower than 1% in the fly ashes
from the combustion of the high-reactivity Berezovsky coal (the Bfa series) [17] and is considerably
higher, i.e., 9.2 and 10.3 wt %, in the fly ashes from the Fields 1 and 2 of the electrostatic precipitators,
respectively, in the case of the low-reactivity Kuznetsk coal (the Mfa series). Moreover, as follows
from Fig. 1, they have a substantially different composition of the mineral part and, according to the
classification by the chemical composition [36], belong to different types of ashes. Note also that fly
ashes from the combustion of the Berezovsky coal belong to high-calcium ashes of the Calsialic (CS)
type and have a high sulfur content (from 6.7 to 15.6 wt % SO3) [17]. In fly ashes of the Kuznetsk coal,
the aluminosilicate component dominates, which is characteristic of ashes of the Sialic (S) type.
One of the technical characteristics used in practice for a dispersed material is the bulk density.
It is a complex parameter, which depends on the composition, the shape of particles, and their size
distribution. As follows from the data presented in Table 2, the bulk density (ρbulk) of fly ashes of the
Bfa series (0.93−0.63 g/cm3) is substantially lower than that of fly ashes of the Mfa series (1.31 and 1.18
g/cm3). At that the decrease in the bulk density decrease along the gas-and-dust flow, which is more
significant for ashes of the Bfa series for the Fields 1 and 2 – from 0.93 to 0.69 g/cm3, respectively,
Fig. 1. Classification of the fly ashes according to the chemical composition [36]: 1 – fly ashes taken from Felds
1−4 of the electrostatic precipitators at the BSDPS-1 (Bfa series); 2 – fly ashes taken from Fields 1−2 of the electrostatic precipitators at the MPP-22 (Mfa series)
# 76 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Olga M. Sharonova, Natalia A. Oreshkina… Aerodynamic Separation of Fly Ashes of Selective Sampling from Pulverized...
Table 3. Granulometric composition of fly ashes for Bfa and Mfa series from the Fields 1 of electrostatic precipitators
Fraction (mm)
Content (wt.%)
Bfa
Mfa
> 0.40
0.03
0.22
0.2-0.4
0.07
0.7
0.16-0.20
0.2
0.72
0.10-0.16
0.7
4.72
0.063-0.1
3.6
13.0
<0.063
95.4
80.6
∑
100
100
and is less significant for ashes of the Mfa series – from 1.31 to 1.18 g/cm3 for the Fields 1 and 2,
respectively.
A comparison of the granulometric compositions of fly ashes of both types taken from the first
field of the electrostatic precipitators (Table 3) showed that the fly ash of the Bfa series differs by its
higher dispersity: the content of the fraction <0.063 mm in this ash is higher than 95%, whereas the
content of the same fraction in the fly ash of the Mfa series does not exceed 81%. The analysis of the
data available in the literature [3,14−16,37] has demonstrated that, according to the content of the
fraction <0.063 mm, the fly ash of the Bfa series from the Field 1 of electrostatic precipitators is close
only to the fly ashes from the first field of EPs from combustion of the Angren brown coal (Angren
SDPS) and the Estonian slate (Estonian SDPS) in TP boilers. The content of the fraction <0.063 μm in
the other ashes from the first fields of EPs amounts to 35−82%.
A similar conclusion regarding the difference between the sizes of particles in the fly ashes under
investigation from the first field of the electrostatic precipitators follows from the data obtained on the
laser analyzer (Table 2). In particular, the average diameter of particles in the Bfa fly ash from the Field
1 is equal to 9.3 μm and the quantity d90 (the size of particles for which the integral distribution curve
corresponds to 90%) is 41.8 μm, whereas these quantities for the Mfa ash amount to 16.0 and 73.6 μm,
respectively. Furthermore, they differ in the shape of the curve and in the position of the maximum:
a unimodal size distribution with the maximum at 21.3 μm and the asymmetry toward submicron
particles is observed for the fly ash of the Bfa series (Fig. 2), and a bimodal size distribution with the
main maximum at 58.7 μm, the secondary maximum at 16.0 μm, and the asymmetry toward micron
and submicron particles is observed for the fly ash of the Mfa series (Fig. 3).
For fly ashes from different fields of electrostatic precipitators, there exists a general tendency
toward a decrease in the size of particles along the gas-and-dust flow, and the regularities observed
in this case for different ash collection facilities are somewhat different. The dispersity can increase
monotonically from the first field to the subsequent fields of electrostatic precipitators [14,15], or the
size of particles in the fly ash from the second field can decrease rather abruptly with respect to that
from the first field [14,39]. There are also examples where the differences between the sizes of particles
in the fly ashes from the first and second fields are less significant compared to the third and fourth
fields [14].
# 77 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Fig. 2. Differential particle size distributions for fly ashes of the Bfa series from Fields 1−4 of the electrostatic
precipitators
Fig. 3. Differential particle size distributions for fly ashes of the Mfa series from Fields 1 and 2 of the electrostatic
precipitators.
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Olga M. Sharonova, Natalia A. Oreshkina… Aerodynamic Separation of Fly Ashes of Selective Sampling from Pulverized...
For ashes of the Bfa series, the average diameter of particles decreases from 9.3 μm to 2.4 μm
when changing over from the Field 1 to the Field 4 of the electrostatic precipitators, respectively, while
the maximum shifts from 21.3 μm to 3.4 μm and the quantity d90 changes from 41.8 μm to 7.9 μm,
respectively (Table 2). The contribution from the particles with a submicron size increases in the same
order (Fig. 2). It follows from these data that the largest decrease in the sizes of the particles is observed
when changing over from the first field to the second field of the electrostatic precipitators.
For fly ashes of the Mfa series, the dispersities of the ashes from the Fields 1 and 2 of the
electrostatic precipitators differ insignificantly: the average size of the particles changes from 16.0 to
15.9 μm, and the quantity d90 changes from 73.6 to 69.9 μm. The positions of the maxima differ more
significantly; i.e., the maxima are observed at 58.7 and 68.8 μm, respectively (Table 2), and become
more pronounced for the fly ash from the Field 2 (Fig. 3).
Therefore, fly ashes of the Bfa series sampled from different fields of electrostatic precipitators at
the BSDPS-1 differ significantly in the bulk density and size of particles, which is especially true for
the fly ashes from the first and second fields of the electrostatic precipitators. By contrast, fly ashes of
the Mfa series sampled from the first and second fields are rather close in size. Moreover, they have a
substantially higher bulk density and a lower dispersity as compared to fly ashes of the Bfa series of
the same sampling.
Fly ash products of aerodynamic separation
The aerodynamic separation makes it possible to more finely separate fly ashes, thus providing
a means for effectively isolating small and light particles from larger and heavier particles, and leads
to the formation of products with a stable composition. Kruger [6] described the air classification of
ashes of the aluminosilicate composition (class F) and obtained fractions <45 μm and <10 μm, which
have close chemical compositions but exhibit different properties. The use of the fraction <45 μm as
a 30% additive to the Portland cement does not lead to a deterioration of its strength characteristics,
and the ultrafine fraction <10 μm, which is characterized by a higher pozzolanic activity, makes it
possible to produce cements for high-strength concretes. Fisher et al. [38] used the aerodynamic
separation of fly ashes to prepare four fractions with average particle diameters of 20.0, 6.3, 3.2,
and 2.2 μm and demonstrated that, with a decrease in the size of the fraction, the contribution from
massive homogeneous microspheres to this fraction increases, whereas the contribution from particles
of a different morphology, especially from cenospheres, decreases. Fractions <5−10 μm, which were
prepared by the air classification of fly ashes, are intended for the use as fillers in thermoplastics,
rubbers, and asphalt, as well as for manufacturing paints [5−7,9,39].
The aerodynamic separation with the formation of three (heavy, medium, and light) products was
performed for the fly ash of the Bfa series from the Fields 1 and 4 of the electrostatic precipitators and
for the fly ash of the Mfa series from the Field 1 of the electrostatic precipitators. For the fly ash of
the Bfa series (Table 2), the bulk density substantially decreases when changing over from the heavy
product to the light product: from 1.13 g/cm3 to 0.79 and 0.58 g/cm3 for the Field 1 and from 0.72 g/
cm3 to 0.59 and 0.35 g/cm3 for the Field 4, respectively. It follows from the data obtained on the laser
analyzer that the fly ash of the Bfa series from the Field 1 of the electrostatic precipitators is effectively
separated and that the products differ significantly in the dispersity (Table 2, Fig. 4). In particular, the
average diameter of particles in the heavy product is equal to 26.2 μm and decreases to 7.5 and 2.8 μm
# 79 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Olga M. Sharonova, Natalia A. Oreshkina… Aerodynamic Separation of Fly Ashes of Selective Sampling from Pulverized...
Fig. 4. Differential particle size distributions for products of the aerodynamic separation of the fly ash of the Bfa
series taken from Field 1 of the electrostatic precipitators (at the top) and optical images of the (a) light and (b)
heavy products
in the medium and light products, respectively; in this case, the quantity d90 decreases from 62.2 μm
to 26.5 and 10.2 μm with a shift of the maximum from 37.9 μm to 13.6 and 4.3 μm, respectively.
A different conclusion follows from the separation of the fly ash of the Bfa series from the Field
4 of the electrostatic precipitators. The size characteristics for the heavy and medium products are
very close to each other, and those for the light product change to a lesser extent as compared to the
corresponding parameters for the fly ash from the first field of the electrostatic precipitators (Table 2,
Fig. 5). The fly ash from the Field 4 of the electrostatic precipitators is characterized by the highest
dispersity and contains micron and submicron particles that are prone to aggregation, which, in turn,
hampers their separation in air flow. Moreover, according to the size of particles, this fly ash without
an additional air separation satisfies requirements for many applications.
The fly ash of the Mfa series was also separated in our setup with a high efficiency. From the fly ash
of the Field 1 with a bimodal particle size distribution, we obtained three products with unimodal size
# 80 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Olga M. Sharonova, Natalia A. Oreshkina… Aerodynamic Separation of Fly Ashes of Selective Sampling from Pulverized...
Fig. 5. Differential particle size distributions for the products of the aerodynamic separation of the fly ash of the
Bfa series taken from Field 4 of the electrostatic precipitators (at the top) and optical images of the (a) light and
(b) heavy products
distributions; in this case, the average diameter of particles decreases from 51.0 μm in the heavy product
to 9.6 and 4.8 μm in the medium and light products, respectively (Table 2). Considerable variations are
observed in the quantity d90, which changes from 97.8 μm to 31.3 and 13.9 μm, respectively, and in
the position of the maximum in the differential distribution curve, which changes from 60.0 μm to 16.8
and 9.5 μm, respectively (Fig. 6). Unlike the fly ashes of the Bfa series, the bulk density ρbulk in this case
reaches a maximum value for the medium product.
The examination of the optical images of the heavy and light products obtained by the separation
of fly ashes of both series (examples of these images are shown in Figs. 4−6) confirms the laser analyzer
data on the differences in the dispersities of the products.
The analysis of the images obtained with an optical microscope in transmitted light allows us
to conclude that the heavy product formed from the fly ash of the Bfa series from the Field 1 (Fig. 4)
predominantly contains particles with sizes ranging from 20 to 60 μm, including transparent and opaque
# 81 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Olga M. Sharonova, Natalia A. Oreshkina… Aerodynamic Separation of Fly Ashes of Selective Sampling from Pulverized...
Fig. 6. Differential particle size distributions for the products of the aerodynamic separation of the fly ash of the
Mfa series taken from Field 1 of the electrostatic precipitators (at the top) and optical images of the (a) light and
(b) heavy products
angular large particles, black opaque magnetic microspheres, and spherical particles with different sizes
and transparencies. The light fraction predominantly contains particles with sizes smaller than 10 μm,
among which the number of microspheres, especially transparent species, increases and the number of
black opaque microspheres and angular particles is considerably smaller. Moreover, the fraction contains
a number of aggregates consisting of micron and submicron particles, which partially disintegrate during
the preparation of samples for microscopic investigations in an emulsion oil or glycerol.
The aggregation of particles is observed to a greater extent in products of the separation of fly
ashes of the Bfa series from the Field 4 of the electrostatic precipitators, which decreases the efficiency
of the separation process. It can be seen from the images shown in Fig. 5 that these products contain
even a larger number of microspherical particles and that the number of transparent microspheres
increases as compared to their number in the products obtained from fly ashes of the same series from
the Field 1.
# 82 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Olga M. Sharonova, Natalia A. Oreshkina… Aerodynamic Separation of Fly Ashes of Selective Sampling from Pulverized...
The products of separation of fly ashes of the Mfa series predominantly contain microspherical
particles of different hues and transparencies; moreover, they include a considerable number of opaque
angular particles of unburned carbon (Fig. 6). It should be noted that the size of particles in the heavy
products exceeds their size in the light product by a factor of more than 10. Moreover, the light products
contain a larger number of transparent microspheres.
Therefore, the performed investigations have demonstrated that the aerodynamic separation makes
it possible to prepare products from fly ashes that differ significantly in the density and size of particles.
Selective sampling of fly ashes from the ash collection facility and variations in the parameters of the
setup for aerodynamic separation allow one to prepare narrow fractions with a controlled density and a
controlled dispersity for the use as additives and fillers in manufacturing cements with improved properties,
concretes, high-quality plastic materials, rubbers, road pavements, paints, and ceramic goods.
Conclusions
A comparative study has been carried out of fly ashes from the pulverized combustion of coals of
two ranks, i.e., the B2 (sub C) rank, which is taken from Fields 1−4 of the electrostatic precipitators at
the BSDPS -1 (the Bfa series), and the T (sa) rank, which is taken from Fields 1 and 2 of the electrostatic
precipitators at the MPP-22 (the Mfa series).
It has been shown that fly ashes of the Bfa series exhibit a lower bulk density, a higher dispersity
and a unimodal particle size distribution. The latter parameter considerably increases along the gasand-dust flow from Field 1 to Field 4 of the electrostatic precipitators. In particular, the average diameter
of particles decreases from 9.3 μm to 2.4 μm and quantity d90 changes from 41.8 μm to 7.9 μm for fly
ashes from the Fields 1 and 4, respectively. The strongest change of density and size of particles occurs
at transition from 1 to 2 field of the EPs.
Fly ashes of the Mfa series are characterized by a higher density, a substantially lower dispersity,
and a bimodal particle size distribution and have rather close size characteristics for fly ashes from the
Fields 1 and 2 of electrostatic precipitators. In particular, the average diameter of particles changes
from 16.0 to 15.9 μm and the quantity d90 changes from 73.6 to 69.9 μm for fly ashes from the Fields
1 and 2, respectively.
The use of the aerodynamic separation has made it possible to obtain three products from fly ashes
of both series which differ significantly in the bulk density and size of particles. In particular, for fly
ashes of the Bfa series from the Field 1, the average diameter of particles decreases from 26.2 μm in
the heavy product to 7.5 and 2.8 μm in the medium and light products; in this case, the quantity d90
decreases from 62.2 μm to 26.5 and 10.2 μm, respectively. For fly ashes of the Mfa series from the Field
1, the average diameter of particles decreases from 51.0 μm in the heavy product to 9.6 and 4.8 μm in
the medium and light products, while the quantity d90 changes significantly from 97.8 μm to 31.3 and
13.9 μm, respectively. Thus, the selective sampling and aerodynamic separation of the fly ashes from
the pulverized combustion of coals are promising techniques for manufacturing ash products with a
controlled density and a controlled dispersity.
Acknowledgments
This study was supported in part by the Siberian Branch of the Russian Academy of Sciences
within the framework of the Interdisciplinary Integration Project No. 77.
# 83 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Olga M. Sharonova, Natalia A. Oreshkina… Aerodynamic Separation of Fly Ashes of Selective Sampling from Pulverized...
References
[1] European Coal Combustion Products Association, 2006. URL: http://www.ecoba.com.
[2] American Coal Ash Association, 2008. URL: http://www.acaa-usa.org.
[3] Sharonova O.M., Anshits N.N. , Orujeinikov A.I., Akimochkina G.V., Salanov A.N., Anshits
A.G. // Chemistry for Sustainable Development. 2003. Vol.11. N. 3. P.639.
[4] Yeremin I.V. , Bronovets T.M. Coal rank composition and rational utilization handbook,
Moscow; Nedra, 1994. p. 254 (in Russian).
[5] Sharonova O.M., Anshits N.N., Yumashev V.V., Anshits A.G. // Fuel. 2008. Vol. 87. N. 10-11. P.
1989.
[6] Ozasa K.// Advanced Clean Coal Technology, International Symposium 2001, Japan, Tokyo,
2001. P. 202
[7] Lewandowski W. , Feuerborn H.J. // Informational and Analytical Center «Ecology in Power
Engineering», 2010. URL: http://www.ecopower.ru.
[8] Kruger R.A.// Informational and Analytical Center «Ecology in Power Engineering», 2010.
URL: http://www.ecopower.ru.
[9] Kruger R.A. // Fuel. 1997. Vol. 76. N. P. 777.
[10] Jo Y.M., Huchison R., Raper J.A.// Waste Management and Research. 1996. Vol.14.
P. 281.
[11] Kruger R., Hovy M., Wardle D. // International Ash Utilization Symposium, 1999, Center for
Applied Energy Research, University of Kentucky, Paper № 72.
[12] Mardon S.M., Hower J.C. // Int. J. Coal Geol. 2004. Vol. 59. P. 153.
[13] Ovcharenko G.I. Coal Ashes of the Kansk-Achinsk Fuel and Energy Complex in Building
Materials, Krasnoyarsk University, Krasnoyarsk, Russia; 1992. p. 216 (in Russian).
[14] Panteleev V.G., Larina E.A., Melentiev V.A., Sergeeva T.E., Mokrushin A.R., Composition
and characteristics of ash and Slag of heat power plants handbook, Leningrad: Energoatomizdat; 1985.
pp. 288, (in Russian).
[15] Lee S.H., Kim K.D., Sakai E., Daimon M. // Journal of the Ceramic Society of Japan. 2003.
Vol. 111. P. 0011
[16] Shanthakumar S., Singh D.N., Phadke R.C. // The 12th International Conference of International
Association for Computer Methods and Advances in Geomechanics (IACMAG), Goa, India, October
1−6, 2008. P. 2261.
[17] Sharonova O.M., Solovyov L.A., Oreshkina N.A., Yumashev V.V., Anshits A.G. // Fuel
Processing Technology. 2010. Vol. 91. P.573.
[18] GOST 25818-91, Fly Ash of Heat Power Plants for Use in Concretes. 1991.
[19] ASTM C618-05. Standard Specification for Coal Fly Ash and Raw or Calcined Natural
Pozzolan for Use in Concrete. 2005.
[20] EN 450:1994, Fly ash for concrete. 1994.
[21] BS 3892:1997, Pulverized fuel ash − Part 1: Specifications for pulverized fuel ash for use with
Portland cement. 1997.
[22] JIS A 6201-1999, Standards for fly ash used for concrete. 1999.
[23] Blanco F., Garcia M.P., Ayala J. // Fuel. 2005. Vol. 84. P. 89.
[24] Poon C.S., Qiao X.C., Lin Z.S. // Cem. Concr. Res. 2003. Vol. 33. P. 1857.
# 84 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Olga M. Sharonova, Natalia A. Oreshkina… Aerodynamic Separation of Fly Ashes of Selective Sampling from Pulverized...
[25] Ronin V. Patent MX PA06012978 (2007). Mexica, http://worldwide.espacenet.com. (14
November 2011)
[26] Ghosal S., Self S.A. // Fuel. 1995. Vol. 74. P. 522.
[27] Vereshchagin S.N., Kurteeva L.I., Anshits A.G. // Chem. Sustainable Dev. 2008. Vol. 16.
P. 529.
[28] Hase Yuji, Saito Tomotaka. Patent 2001090012 (2001) Japan. http://worldwide.espacenet.
com. (14 November 2011)
[29] Long Cao. Patent 1095049 (1994)China, http://worldwide.espacenet.com. (14 November
2011)
[30] Vdovchenko V.S., Martynova M.I., Novitskii N.V., Yushina G.D., Power plant fuel in the USSR:
Handbook, Moscow : Energoatomizdat, 1991. pp. 184 (in Russian).
[31] Chernyshev E.V., Zuev O.G. , Borovkin B.A., Koptev A.S., Shanin M.P. // Elektr. Stn., 1999,
11 (in Russian).
[32] GOST 23148-98, Powders Used in Powder Metallurgy: Sampling. 1998.
[33] GOST 26565-85, Unmolded Refractories: Methods for Withdrawing and Preparing Samples.
1985.
[34] GOST 18318-94, Metallic Powders: Determination of the Particle Size Using Dry Sieving.
1994.
[35] GOST 5382-91, Concretes and Materials for Use of Concrete Production: Methods of
Chemical Analysis. 1991.
[36] Vassilev S., Vassileva C. // Fuel. 2007. Vol. 86. P.1490.
[37] Valentim B., Guedes A., Flores D., Ward C.R., Hover J.C. // Coal Combustion and Gasification
Products. 2009. Vol.1. P.14.
[38] Fisher G.L., Prentice B.A., Silberman D., Ondov J.M., Bierman A.H., Ragaini R.C., McFarland
A.R. // Environ. Sci. Technol. 1978. Vol. 12. P. 447.
[39] William S.R., Seyi A.K. Patent 6139960 (2000) USA, http://worldwide.espacenet.com. (14
November 2011)
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Olga M. Sharonova, Natalia A. Oreshkina… Aerodynamic Separation of Fly Ashes of Selective Sampling from Pulverized...
Аэродинамическое разделение летучей золы
селективного отбора от пылевидного сжигания
разных марок углей
О.М. Шароноваa, Н.А. Орешкинаб,
Л.И. Куртееваa, А.Г. Аншицa,б
a
Институт химии и химической технологии,
Сибирское отделение РАН,
Россия 660049, Красноярск, ул. К.Маркса, 42
б
Сибирский федеральный университет,
Россия 660041, Красноярск, пр. Свободный, 79
Изучены летучие золы от пылевидного сжигания двух марок углей – марки Б2, отобранные с
полей 1-4 электрофильтров на БГРЭС-1 (серия Bfa), и марки Т, отобранные c полей 1 и 2 ЭФ на
Московской ТЭЦ-22 (серия Mfa).
Показано, что летучие золы серии Bfa имеют пониженную насыпную плотность, более
высокую дисперсность и мономодальное распределение частиц по размеру. При этом размер
частиц значительно уменьшается по ходу газопылевого потока, в особенности от 1-го ко
2-му полю ЭФ. Золы серии Mfa отличаются более высокой плотностью, имеют существенно
меньшую дисперсность, бимодальное распределение частиц по размеру и достаточно близкие
размерные характеристики для зол поля 1 и 2 ЭФ. Аэродинамическое разделение позволило
получить из летучих зол обеих серий три продукта, существенно отличающиеся по плотности
и размеру частиц. Комбинирование селективного отбора и аэродинамического разделения
летучих зол от пылевидного сжигания углей является перспективным способом для получения
зольных продуктов заданной плотности и дисперсности.
Ключевые слова:
разделение.
летучая
зола,
электрофильтры,
дисперсность,
аэродинамическое
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Journal of Siberian Federal University. Engineering & Technologies 1 (2012 5) 87-97
~~~
УДК 577.355(571.51)
Determination of Chlorophyll Photosynthetic
Potential in Vegetation Using Ground-Based
and Satellite Methods
Irina Yu. Botvicha*, Alexander F. Sidkoa,
Tamara I. Pismana,b and Anatoly P. Shevyrnogova,b
a
Institute of Biophysics Siberian Branch of the RAS,
50-50 Akademgorodok, Krasnoyarsk, 660036 Russia
b
Siberian Federal University,
79 Svobodny, Krasnoyarsk, 660041 Russia 1
Received 6.02.2012, received in revised form 13.02.2012, accepted 20.02.2012
Potential productivity of agricultural crops can be evaluated using a new method based on chlorophyll
photosynthetic potential (CPSP) derived from satellite information. The chlorophyll photosynthetic
potential is determined by the amount of light absorbed by the plants in the region of the red band of
chlorophyll absorption in a definite growth period. The value of parameter CPSP is determined as the
area of the triangle whose vertex coordinates are exoatmospheric reflectance in the green, red, and
near infrared bands and the respective mean wavelengths of bands. CPSPs of the crops have been
examined using medium and high spatial resolution satellite (MODIS/Terra and Landsat7 ETM+)
data obtained during the growing season of the plants. CPSP and crop productivity have been found
to be interrelated. The obtained satellite evaluations are in good agreement with the ground-truth
observation data.
Keywords: chlorophyll photosynthetic potential; satellite data; agricultural plants
In our previous studies, in order to quantify the relationship between the spectral brightness
coefficients (SBC) of the crops and chlorophyll “a” concentration in the leaves of the upper plant layer,
we proposed using chlorophyll photosynthetic potential (CPSP) of plants. This optical remote sensing
technique of determining CPSP was based on the data of ground-truth measurements of SBC. CPSP
is the most comprehensive indicator of the relationship between SBC and plants’ physiological and
biometric parameters. The chlorophyll photosynthetic potential is determined by the amount of light
absorbed by the plants in the region of the red band of chlorophyll absorption in a definite growth
period.
In this paper we present a new method for productivity evaluation of agricultural crops, based
on chlorophyll photosynthetic potential derived from satellite data. The value of parameter CPSP is
determined as the area of the triangle whose vertex coordinates are exoatmospheric reflectance in the
green, red, and near infrared bands and the respective mean wavelengths of bands. CPSPs of the crops
*
1
Corresponding author E-mail address: irina.pugacheva@mail.ru
© Siberian Federal University. All rights reserved
# 87 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Irina Yu. Botvich, Alexander F. Sidko… Determination of Chlorophyll Photosynthetic Potential in Vegetation Using...
have been examined using medium and high spatial resolution satellite (MODIS/Terra and Landsat7
ETM+) data obtained during the growing season of the plants.
The values of chlorophyll photosynthetic potential obtained using MODIS/Terra satellite
data were found to be related to the degree of weed contamination of the crops on the territory
of the Krasnoyarskii Krai and the Republic of Khakasia (Russia). Chlorophyll photosynthetic
potential (based on Landsat 7 ETM+) is an indicator of crop productivity in the study areas.
Chlorophyll photosynthetic potential of the crops based on Landsat 7 ETM+ satellite data can
be used to determine crop species composition. CPSP and crop productivity have been found to
be interrelated. The obtained satellite evaluations are in good agreement with the ground-truth
observation data.
Introduction
The major objectives of satellite monitoring of terrestrial vegetation cover are to identify
agricultural areas, determine their species composition [1, 2] and morphophysiological properties
of plants [3, 4, and 5]. Space imaging affords ways to both enhance agricultural statistics gathering,
enabling more accurate, uniform, objective and frequent observations, and essentially improve
methods of global and regional online observation and monitoring of the crops and crop productivity
forecasting [6, 7].
The reflectance of vegetation cover is a parameter providing very much information on the
physiological state of the plants. For instance, based on intensity and spectral composition of the light
reflected by the vegetation cover, one can determine the pigment and species composition of plants
and their moisture content, the state of the surface, and the architectonics of phytoelements. The
most significant contribution to the formation of the spectral portrait of vegetation cover is made by
green pigments – chlorophylls. Up to 95% variations in spectral brightness coefficients (SBC) of the
vegetation cover within the visible spectrum result from changes in chlorophyll content. Carotenoids
absorb 10% to 20% of the total radiant energy taken up by plant pigments. Hence, it is important to
know how chlorophyll content is related to SBC of vegetation cover in order to be able to develop
remote sensing methods and accurately interpret the data obtained using these methods for evaluation
of plant productivity [8, 9].
In our previous studies, in order to quantify the relationship between the SBC of the crops and
chlorophyll “a” concentration in the leaves of the upper plant layer, we proposed using chlorophyll
photosynthetic potential (CPSP) of plants [9, 10]. The proposed optical remote sensing technique of
determining CPSP was based on the data of ground-truth measurements of SBC.
The purpose of this study was to assess the feasibility of evaluation of crop productivity based on
CPSP calculated from satellite data.
The major objectives of the study were as follows:
- to develop the method for calculating CPSP from satellite data, using the technique of
quantifying CPSP based on ground-truth spectrophotometric data;
- to study the CPSP of agricultural crops calculated from medium- and high-resolution satellite
data (MODIS/Terra and Landsat 7 ETM+) throughout plant growing season;
- to determine the relationship between the CPSP of crop communities and their
productivity.
# 88 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Irina Yu. Botvich, Alexander F. Sidko… Determination of Chlorophyll Photosynthetic Potential in Vegetation Using...
1. Methods and material
Previous ground-truth studies [8, 9, 11, 12, 13] of the SBC of crops showed that CPSP is determined
by the amount of the light in the chlorophyll “a” absorption region of the red band (λmax= 680 nm)
absorbed by plants during their growing season (Fig. 1). CPSP is the most comprehensive indicator
of the relationship between SBC and plants’ physiological and biometric parameters. CPSP (S) is
determined from the following formula:
730
S
90 ˜ U 730 t U550 t ³ U O , t dO ,
(1)
550
where ρ 550 and ρ 730 are the mean values of crop SBC at λ = 550 and λ = 730 nm, with λmax = 680 nm
(wavelength that reflects maximum of absorption chlorophyll); t – the limits of the time interval.
Wheat (Triticum aestivum L.), buckwheat (Fagopyrum esculentum) and oat (Avena sativa L.) fields
were studied regularly during growth period. Study sites were situated in the Minusinskii District of
the Krasnoyarskii Krai (wheat and buckwheat fields) and in the Altaiskii District of the Republic of
Khakasia (oat fields) (Fig. 2). Examination of the plant communities was performed on fixed plots of
total area 5300 ha, using conventional geobotanical methods, throughout the plant growing season of
2006 [14]. To quantify green aboveground biomass, 3 to 5 replicate samples were collected from 1×1 m2
plots and weighed. Study site coordinates were recorded using a GPS navigator.
CPSPs of wheat, oats, and buckwheat were studied using information from Landsat 7 ETM+ and
MODIS/Terra satellites. The study was based on the following satellite data:
• MODIS/Terra (MOD09GHK, MOD09GQK products) data in the visible (band 4: 545 – 565
nm), (band 1: 620 – 670 nm) and near infrared (band 2: 841 – 876 nm) spectral regions, of 250-m
Fig.1. Wheat crop SBC vs. wavelength
# 89 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Irina Yu. Botvich, Alexander F. Sidko… Determination of Chlorophyll Photosynthetic Potential in Vegetation Using...
Fig. 2 Schematic map of the study area location (c). The study area of wheat crops (№ 1-3), buckwheat crops (№ 4,
5) were situated in the Minusinskii District of the Krasnoyarskii Krai (a); the study area of oat crops (№ 6, 7) were
situated in the Altaiskii District of the Republic of Khakasia (b)
spatial resolution. A MODIS 250 m, daily composite imagery (MOD09) product developed by
the NASA MODIS MODLAND Science Team, is available to the user community (http://
modis.gsfc.nasa.gov/).
• Landsat 7 ETM+ data in the green (band 2: 525 – 605 nm), red (band 3: 630 – 690 nm), and near
infrared (band 4: 750 – 900 nm) spectral regions, of 30-m spatial resolution.
The crops were monitored throughout their growing season in the Krasnoyarskii Krai and the
Republic of Khakasia (Russia) using MODIS/Terra images for the time period from May 10 through
September 10, 2006. Evaluation of CPSP based on high spatial resolution satellite data and determination
of geographical coordinates of study areas were performed using satellite data obtained by Landsat 7
ETM+ on September 2, 2006.
Satellite data were processed using the ENVI software. The obtained data were statistically
processed and visualized (by graphing) using the Statistica software. Statistical processing involved
calculation of the following parameters: mean, minimum, maximum values and standard deviation
from the mean. The Pearson’s correlation analysis evaluating the association between CPSP and the
phytomasses of plants was carried out.
1.1 Satellite data processing
Preprocessing of Modis/Terra satellite data
MODIS/Terra satellite data were processed step by step. MOD09GQK and MOD09GHK image
projections were reprojected from Sinusoidal projection into Universal Transverse Mercator (UTM),
using MODIS Reprojection Tool, and, simultaneously, reformatted from HDF to GeoTIFF. Analysis
# 90 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Irina Yu. Botvich, Alexander F. Sidko… Determination of Chlorophyll Photosynthetic Potential in Vegetation Using...
of vegetation development using satellite data can only be based on time series of observation data
that would be free of such negative factors as cloud cover and cloud shadows. Cloud masks were
constructed using the “Maska” software program written in IDL. The program is based on the following
principles:
1. Observations made at the zenith angle above 45 degrees are excluded from the initial satellite
data [7];
2. Cloud, snow cover, and cloud shadow pixels are detected, using measurements of reflectance
in MODIS band 3 (459-479 nm) and band 6 (1628-1652 nm), which are normed to the cosine
of the solar zenith angle at the viewpoint, in order to compensate for lighting conditions. The
algorithm also uses Normalized Difference Snow Index (NDSI) [7];
3. Cloud shadow areas are excluded from further examination in the satellite image, based on the
data on the height of the cloud cover and information on the positions of the satellite and the
Sun (MODMGGAD product) [7].
Preprocessing of Landsat 7 ETM + satellite data
The image used was collected by the Landsat 7 ETM+ sensor on September 2, 2006 (Path 143,
Row 023). It is a standard L1T product. This image was radiometrically and geometrically corrected
in the EROS Data Center of the U.S. Geological Survey (USGS) using standard methods. Raw digital
numbers associated with the ETM+ pixels were converted to spectral Radiance using ENVI software
programs (TM Calibration Parameters Module).
Atmospheric correction was carried out using the ENVI FLAASH (Fast Line_of_sight
Atmospheric Analysis of Spectral Hypercubes) module. FLAASH was developed by Spectral Sciences,
Inc., a world leader in optical phenomenology research, under the sponsorship of the U.S. Air Force
Research Laboratory. The module is based on proven MODTRAN (MODerate resolution atmospheric
TRANsmission), code and algorithms, producing accurate results based on every image’s unique
parameters. Spectral Radiance was converted to exoatmospheric reflectance by Flaash.
The cloud cover assessment is detected using the algorithm proposed by Richard R. Irish (the
article “Landsat 7 Automatic Cloud Cover Assessment”). This approach is based on the two passes
through ETM+ data. In pass one, the reflective and thermal properties of scene features are used to
establish the presence or absence of clouds in a scene. If present, a scene-specific thermal profile for
clouds is established. In pass two, a unique thermal signature for clouds is developed and used to
identify the remaining clouds in a scene.
1.2. Thematic image processing
Every pigment (chlorophyll, carotenoid) is capable of selectively absorbing solar energy of different
spectral regions, essentially determining the coloration of plant communities during the course of their
growing season, when the composition and the amounts of the pigments in plants vary changeable [9,
13].
The relationship between chlorophyll content of plants and crop productivity was assessed using
exoatmospheric reflectance (R) of plants via S – chlorophyll photosynthetic potential. The chlorophyll
photosynthetic potential depends on the amount of light absorbed by the plants in the region of the red
band of chlorophyll absorption in a definite growth period. The value of parameter S is determined as
# 91 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Irina Yu. Botvich, Alexander F. Sidko… Determination of Chlorophyll Photosynthetic Potential in Vegetation Using...
Fig. 3. A schematic representation of the triangle with vertices (λi, r i) used to calculate CPSP
the area of the triangle whose vertex coordinates are exoatmospheric reflectance (ri) in the green, red,
and near infrared bands and the respective mean wavelengths (λi) of bands (Fig. 3).
The area of the triangle is calculated using geometric properties of the vector product from the
following formula:
S
O 2 O1 ˜ r3 r1 O3 O1 ˜ r1 r2 2,
(2)
where λi is the wavelength; r i is the R value; i=1,2,3.
Chlorophyll photosynthetic potential calculation based on satellite information uses:
• MODIS/Terra data for the green (545 – 565 nm), red (620 – 670 nm) and near infrared (841 –
876 nm) bands,
• Landsat 7 ETM + data for the green (525 – 625 nm), red (630 – 690 nm) and near infrared (750 –
900 nm) bands.
2. Results and Discussion
Ground-based field studies yielded the following results. Analysis of seasonal variations in the
CPSP (∑S(t)) determined from the crops’ SBC, dry and green biomasses, and chlorophyll concentration
(Cchl), using two-dimensional criterion for random values, showed a rather close positive relationship
between these parameters. The¦highest correlation coefficients (R = 0.9 – 0.95) were recorded between
∑S(t) and ∑Сchl(t) for the crop growing season. Correlation coefficients between ∑S(t) and dry and
green biomasses, however, were 0.75 – 0.85 [8]. Thus, crop CPSP can be used to evaluate potential crop
productivity. Our calculations also showed¦a close interrelationship
between ∑S(t) and productivity of
¦
one crop species (e.g., 30 wheat
¦ varieties), which did not depend on either the variety of the crops or
their growing and processing conditions. Correlation coefficients reached 0.85 – 0.90. The higher ∑S(t)
for the growing season, the greater the productivity in the study plots and fields [8].
Figure 4 shows the relationship of crop productivity (U) to ∑S(t). The graph was plotted based
on the data on the study crops collected over many years (over 110 variants), regardless of differences
# 92 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Irina Yu. Botvich, Alexander F. Sidko… Determination of Chlorophyll Photosynthetic Potential in Vegetation Using...
Fig. 4. Productivities of different crops vs. CPSP
in their sowing and pre-sowing treatment, the amounts of fertilizers (N, K, P) added per square
meter, or crop varieties (wheat, barley, or oat), which evidently had different ripening schedules and
productivities [8].
Below we present results of processing of medium spatial resolution satellite data (MODIS/Terra)
for the crops, collected during their growing season. Graphs of variations in CPSP and the phytomasses
of oat crops are shown in Figure 5. There is a clear difference between the CPSP curves. As tillage,
fertilizer application, seed sowing schedule, and watering were similar, this difference is mainly
accounted for by the extent to which the crops were contaminated with weeds. Weed contamination of
the oat crops in Plot № 2 was higher than in Plot № 1 throughout the oat growing season. Thus, CPSP
can be used as an indicator of weed contamination of the crops.
Table 1 presents correlation coefficients between CPSP and the mass of reproductive organs,
the vegetative mass, the mass of weeds, and the total phytomass of oat, wheat, and buckwheat crops.
For evaluation of weed contamination of crop communities, the table presents correlation coefficients
between CPSP and the phytomass of the crops with different degrees of weed contamination. If the
mass of weeds is greater than 20% of the total phytomass (100%), the crops are considered to be
“contaminated to a high degree”. If the mass of weeds is smaller than or equal to 20% of the total
phytomass, the crops are considered to be “contaminated to a low degree”. Processing of the MODIS/
Terra satellite data showed that for the fields with a high degree of weed contamination there is a good
correlation between CPSP and the weed mass (Table 1).
In addition to this, the data listed in Table 1 showed that the correlation coefficients between CPSP
and the mass of reproductive organs higher than 0.87 (for the crops with low weed contamination)
were indicative of the close relationship between these parameters. Analysis of the results suggests
that CPSP can be used to evaluate potential crop productivity based on both satellite information
and ground-truth spectrometric measurements. In both cases, correlation coefficients between crop
productivity and CPSP (for the crop with low weed contamination) amount to at least 0.85.
# 93 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Fig. 5. Variations in CPSP and the phytomasses of oat plant communities over the growing season:
mass of reproductive organs of oat field №1
mass of reproductive organs of oat field №2
vegetative mass of oat field №1
vegetative mass of oat field №2
weed mass of oat field №1
weed mass of oat field №2
total phytomass of oat field №1
total phytomass of oat field №2
CPSP of oat field №1
CPSP of oat field №2
Table 1 Results of Pearson’s correlation analysis evaluating the association between CPSP and the phytomasses
of wheat, oat, and buckwheat plants for the plant growing season of 2006, p<0.05
Pearson’s correlation between CPSP and phytomass
mass of reproductive
vegetative mass
organs
oat field No. 1
“contaminated to a low degree”
oat field No. 2
“contaminated to a high degree”
wheat field No. 1
“contaminated to a low degree”
wheat field No. 2
“contaminated to a high degree”
buckwheat field No. 1
“contaminated to a low degree”
buckwheat field No. 2
“contaminated to a high degree”
weed mass
total
phytomass
0,91
0,72
0,33
0,53
0,2
0,91
0,64
0,74
0,87
0,9
0,12
0,7
0,7
0,44
0,75
0,9
n/a
0,83
0,05
0,69
n/a
0,75
0,44
0,75
(n/a - not available data)
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Irina Yu. Botvich, Alexander F. Sidko… Determination of Chlorophyll Photosynthetic Potential in Vegetation Using...
Table 2 Results of Pearson’s correlation analysis evaluating the association between CPSP and the phytomasses
of wheat plants (September 2, 2006), p<0.05
Pearson ‘s correlation between CPSP and phytomass
mass of reproductive
organs
vegetative mass
weed mass
total phytomass
-0,9
-0,12
-0,6
-0,72
Fig. 6. Graphs of exoatmospheric reflectance of crops obtained using Landsat7 ETM+ satellite data (bands 2, 3,
4) of September 2, 2006
The data presented below are the results of processing of high spatial resolution satellite data
(Landsat 7 ETM+) of September 2, 2006. Figure 6 shows graphs of exoatmospheric reflectance of
wheat, buckwheat, and oat crops based on Landsat 7 ETM+ satellite data. The graphs were plotted
using the data of three bands (2, 3, and 4); the band numbers are recorded along the abscissa and
reflectances in each of the bands – along the ordinate. The points in the graph are connected by straight
lines. To improve the visual perception of the information, we constructed triangles by connecting
points (r2, b2) and (r4, b4) with a solid line (where b2, b4 are numbers of the bands and r2 , r4 are
reflectances in bands 2 and 4, respectively). Having examined the spatial arrangement of the triangles
within the (λ, r) coordinates, the areas of the triangles, and the angles of the triangles, we concluded
that the triangles were different. Thus, chlorophyll photosynthetic potential of the crops can be used to
determine crop species composition.
Table 2 presents correlation coefficients between CPSP and the mass of reproductive organs, the
vegetative mass, the mass of weeds, and the total phytomass of wheat crops. The CPSP for three study
plots were calculated using Landsat7 ETM+. Correlation coefficients between CPSP and wheat crop
phytomass show that CPSP is closely related to the mass of reproductive organs of the study crops
# 95 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Irina Yu. Botvich, Alexander F. Sidko… Determination of Chlorophyll Photosynthetic Potential in Vegetation Using...
(Table 2). The SBC curve constructed using autumn data flattens (the area of S is decreased and so is
CPSP). Thus, as the mass of reproductive organs increases, CPSP goes lower.
3. Conclusions
1. A method for calculating chlorophyll photosynthetic potential of agricultural crops using
MODIS/Terra and Landsat 7 ETM+ satellite data was proposed.
2. The values of chlorophyll photosynthetic potential obtained using MODIS/Terra satellite data
were found to be related to the degree of weed contamination of the crops on the territory of
the Krasnoyarskii Krai and the Republic of Khakasia (Russia).
3. Chlorophyll photosynthetic potential (based on Landsat 7 ETM+) is an indicator of crop
productivity in the study areas.
4. Chlorophyll photosynthetic potential of the crops based on Landsat 7 ETM+ satellite data can
be used to determine crop species composition.
Acknowledgments
The study was supported by the project “Development of methods for space monitoring of forests
in the Krasnoyarskii Krai in order to estimate the dynamics of biodiversity under natural and human
impacts”, Subprogram No. 23 “Diversity and monitoring of forest ecosystems in Russia”.
References
[1] Erol H., and Akdeniz F. A multispectral classification algorithm for classifying parcels in an
agricultural region// Int. J. Remote Sens. 1996. vol. 17. pp. 357–3371.
[2] Grignetti A., Salvatori R., Casacchia R., and Manes F. Mediterranean vegetation analysis by
multitemporal satellite sensor data// Int. J. Remote Sens. 1997. vol. 18. pp. 1307–1318.
[3] Roberts D. A., Smith M. O., and Adams J. B. Green vegetation, nonphotosynthetic vegetation,
and soils in AVIRIS data// Remote Sens. Environ. 1993. vol. 44. pp. 255–269.
[4] Pax-Lenney M., and Woodcock C. E. The effect of spatial resolution on the ability to monitor
the status of agricultural lands. Remote Sens. Environ. 1997. vol. 61. pp. 210–220.
[5] Clevers J. G. P. W. A simplified approach for yield prediction of sugar beet based on optical
remote sensing data// Remote Sens. Environ. 1997. vol. 61. pp. 221–228.
[6] Santhosh K. Seelan, Soizik Laguette, Grant M. Casady, George A. Seielstad Remote
sensing applications for precision agriculture: A learning community approach// Remote Sensing of
Environmen. 2003. vol. 88. pp. 157–169.
[7] Bartalev S.A., Lupyan E.A., Neishtadt I.A., Savin I.Y. Classification of some types of
agricultural crops in south Russia, using Modis satellite data// Issledovaniye Zemli iz kosmosa (Invest.
Earth Space). 2006. № 3. pp. 68–75 (in Russian).
[8] Sid’ko A.F., Shevyrnogov A.P. Seasonal dependence of the spectral brightness of agricultural
crops on plant chlorophyll content and physiological parameters// Earth.Obs.Rem.Sens, 2000. vol. 16.
pp. 487-500.
[9] Sidko A.F. A remote sensing technique of determining chlorophyll photosynthetic potential of
agricultural crops: wheat, barley, and oats // Izvestiya Akademii nauk. Ser. biol. nauk (Proc. Acad. Sci.
Biol. Ser.). 2004. № 5. pp. 547–555 (in Russian).
# 96 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Irina Yu. Botvich, Alexander F. Sidko… Determination of Chlorophyll Photosynthetic Potential in Vegetation Using...
[10] Sidko A.F., Shevyrnogov A.P. A study of seasonal relationship between the spectral brightness
of agricultural crops and chlorophyll content and physiological parameters of the plants // Issledovaniye
Zemli iz kosmosa (Invest. Earth Space). 1998. № 3. pp. 96–105 (in Russian).
[11] Sidko A.F., Shevyrnogov A.P. Spectral luminance factor of plants as a basis for remotesensing of agricultural cultures sowings// Doklady Akademii Nauk (Proc. Russian Acad. Sci.). 1997.
vol. 354. № 1. pp. 120–122.
[12] Sid`ko A. Remote Assay for Chlorophyll Photosynthetic Potential of Grops on the Example of
Wheat// Biology Bulletin of the Russian Academy of Sciences. 2004. vol. 31. № 5. pp. 450-456.
[13] Shevyrnogov A.P., Vysotskaya G., Sid`ko A., Dunaev K. Typification of natural seasonal
dynamics of vegetation to reveal impact of land surface change on environment (by satellite data) //
Adv. Space Res. 2000. vol. 26. № 7. pp. 1169-1172.
[14] Zhukova Yelena Yu., Shevyrnogov Anatoliy P., Zhukova Vera M., Zorkina Taisia M.,
Pugacheva Irina Yu. Seasonal dynamics of production of agrocoenosises of the south of the minusinsk
hollow// Vestnik Tomskogo Gosudarstvennogo universiteta (Vestnik of Tomsk State University). 2009.
vol. 323. pp. 354- 357 (in Russian).
Изучение хлорофилльного
фотосинтетического потенциала растительности
Юга Красноярского края и Республики Хакасия
наземными и спутниковыми методами
И.Ю. Ботвича, А.Ф. Сидькоа,
Т.И. Письмана,б, А.П. Шевырногова,б
а
Институт Биофизики СО РАН
Россия 660036, Красноярск, Академгородок, 50-50
б
Сибирский федеральный университет
Россия 660041, Красноярск, пр. Свободный, 79
Представлена
новая
методика
оценки
потенциальной
урожайности
посевов
сельскохозяйственных культур с помощью хлорофилльного фотосинтетического потенциала
(ХФСП) на основе спутниковой информации. ХФСП определяется количеством поглощенного
растениями света в области красной полосы спектра поглощения хлорофилла в определенные
периоды вегетации. Величина параметра ХФСП определяется как площадь треугольника,
координатами вершин которого являются значения СКО в зеленом, красном и ближнем
инфракрасном каналах и средние значения длин волн в соответствующих каналах. Проведено
изучение ХФСП сельскохозяйственных культур по спутниковым данным среднего и высокого
пространственного разрешения (MODIS/Terra и Landsat7 ETM+) в течение периода вегетации.
Определена степень взаимосвязи ХФСП и урожайности. Полученные спутниковые оценки
хорошо согласуются с наземными наблюдениями.
Ключевые слова: хлорофилльный фотосинтетический потенциал; спутниковые данные;
посевы сельскохозяйственных культур.
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Journal of Siberian Federal University. Engineering & Technologies 1 (2012 5) 98-110
~~~
УДК 621.9: 621.89
Static Characteristics of Active Hydrostatic
Two-Row Radial Bearing with Restriction
of Output Lubricant Flow
Vladimir A. Kodnyanko*
Siberian Federal University,
79 Svobodny, Krasnoyarsk, 660041 Russia 1
Received 6.02.2012, received in revised form 13.02.2012, accepted 20.02.2012
The design of active hydrostatic radial bearing with smooth cylindrical surfaces and lubricant output
flow restrictors in the form of movable rings on membrane suspension is presented.
The device is several times less power-consuming compared with known devices of flow control. The
bearing has a negative and zero compliance (infinite stiffness), and therefore can be used in machine
tools to suppress the negative influence of elastic system deformation on the accuracy of processing.
On the basis of two-dimensional model of lubricant flow developed a mathematical model, method
and procedure for calculating the bearing load capacity and flow rate. It is established, that the
calculation of static characteristics of bearing in the entire range of operating loads can be correctly
performed only on the base of two-dimensional model. For small eccentricities the characteristic of
zero and negative compliance can be calculated with sufficient accuracy by the simplified method,
based on one-dimensional motion of lubricant flow. Bearing of zero or negative compliance have
load capacity range, which is 20 – 50% more than conventional bearings of the same overall
dimensions. The setting of input throttling slits resistance decisive influence on the optimal static
characteristics of the bearing. The optimal values of its resistance for conventional and active
bearing are practically identical.
Keywords: energy-saving, hydrostatic bearing, zero compliance, negative compliance, infinite
stiffness, smooth cylindrical surface.
Introduction
In the radial gas and hydrostatic bearings (HB) with input flow regulators (AGH-IR), which are
able to exert an active influence on the performance of important characteristics, in particular the
significant decrease in compliance of the load-carrying lubricant film to zero and negative values
[1]. These bearings can be used in machine tools in order to diminish the negative impact of the
deformation of the elastic system on accuracy. The active bearing of this type characterized by two
major disadvantages – to maintain their working ability requires a considerable flow rate [2] and in this
connection bearing shafts may make limited movement only.
In [3], an opposite principle to control the flow of hydraulic fluid, which is used in bearings the
restrictors of output lubricant flow (AGH-OR). In comparison with the AGH-IR such designs have a
*
1
Corresponding author E-mail address: kowlad@rambler.ru
© Siberian Federal University. All rights reserved
# 98 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir A. Kodnyanko. Static Characteristics of Active Hydrostatic Two-Row Radial Bearing with Restriction...
much better static characteristics and have no disadvantages of AGH-IR [4]. In addition, for example,
not full scope radial or open thrust AGH-OR with little or no increase in energy consumption can make
a significant movement of mobile elements, in which the AGH-IR would obviously inoperable [5]. This
unique property of AGH-OR opens up the possibility of their to use not only in machine tools, but in
micro shock absorbers, energy-saving hovercrafts, low flow rate aerostatic suspensions of ground and
suspended trains and high-load machines.
In [4] presents the method of calculating the characteristics profiled AGH-OR on the basis of
one-dimensional flow model of working fluid in a thin lubricating gaps. In conventional radial bearing
more often used smooth cylindrical compounds, which are distinguished by simplicity of design and
manufacturing technology. For such GH one-dimensional models are not suitable because it does not
take into account the influence of circumferential lubricant overflows on the characteristics of bearings,
which leads to significant errors of calculations. In this paper presents design, mathematical model,
calculating procedure and results of study for radial hydrostatic bearing of this type with smooth
cylindrical working surfaces.
Principle of the bearing
In Fig. 1 shows a longitudinal section of the bearing. The design contains a shaft 2 and the housing
1 with a throttling slits 3, through which fluid from source under pressure pS enters into the bearing. At
the ends of housing there are circular protrusions on the inner side of which are a ring-type membranes
4, forming with the housing and protrusions the cavity 5, hydraulically connected with the slits 3. To
the membrane tightly attached rigid rings 6 with cylindrical surfaces forming with the housing axial
gaps 8. The inner surface of the housing 1 and the surface of shaft 2 form slot gap 9 of main loadcarrying lubricant film of thickness h, and the surface of the shaft and ring elements 6 – end load-
Fig. 1. Scheme of bearing
# 99 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir A. Kodnyanko. Static Characteristics of Active Hydrostatic Two-Row Radial Bearing with Restriction...
carrying lubricant film with thickness ht. Membrane hermetically attached to the outer edge of the
housing and by inner edge – to the rings.
Bearing operates as follows. The pumped lubricant, overcoming the hydraulic resistance of slits 3,
comes into a thin lubricating gap of cavities 5 and through the channels 8 in carrying lubricating gaps
9 and 7, and then flows from the bearing. Hydraulic forces which created by pressure on the membrane
in cavities 5 and gaps 7, counterbalanced by force of elastic membrane material deformation. The
integral reaction of the pressure forces on the rings 6 in gaps of the cavities 5 is always greater than the
same reaction from side of gaps 7, as can be verified by analyzing the pressure distribution diagrams of
lubrication in these gaps. Therefore rings 6 will always carry opposite direction against the direction
of external force vector f, thereby creating a greater obstruction for lubricant outflow from the bearing
and raising the hydraulic force of the impact on displaced under the load shaft. Depending on the
membrane flexibility, affecting the motion characteristics of rings 6, it provides decrease of capacity
to zero and negative values.
Mathematical modeling
In modeling of AGH-OR work expected that observed parallelism of axes housing, shaft and
movable rings (Fig. 1). Calculations were carried out using dimensionless quantities. For the scales
adopted: radius r0 of shaft – for linear dimensions, supply pressure pS – for pressures, πr02pS – for loads
and forces, πh03ps /6μ – for volumetric flow rates, h0 – for gaps, clearances and eccentricities of moving
parts, where h0 – thickness h of lubricant film in the coaxial arrangement of movable elements (at f = 0),
μ – viscosity of lubricant. Dimensionless variables are indicated by capital letters.
Due to the symmetry of bearing was considered his right half (Fig. 1). For the central (c-area) and
end (t-area) parts of the base gap and for gap of cavity 5 (m-area), introduced local coordinate systems.
Longitudinal coordinate Z is measured from the left edge of areas, and the circumferential coordinate
φ from the line in which indicated vector of external force f.
Function of dimensionless pressure P(Z, φ) in thin lubricating gaps for an incompressible lubricant
satisfies the stationary differential Reynolds equation [6]
w § 3 wP · w § 3 wP ·
¨ Hk
¸ 0,
¨ Hk
¸
wM ©
wM ¹ wZ ©
wZ ¹
(1)
where function of dimensionless thickness of the lubricant film for c-area H(φ) = 1 – ε Cos(φ), for t-area
Ht(φ) = Ht0 – εt Cos(φ), for the m-area Hm(Z, φ) = [Hm0 – εm Cos(φ)](1 – Z/L1); ε – eccentricity of shaft and
housing, εt – eccentricity of rings and shaft, εm – eccentricity of rings and housing.
For pressures are obvious conditions
wPc
(0, M )
wZ
0, Pt (L1 ,M )
0, Pc ,t ,m (Z,M )
Pc ,t ,m (Z,2S -M ) ,
(2)
where L – dimensionless half-length of bearing (summary length of c-area and t-area), L1 – dimensionless
length of rings (length of c-area and m-area). The first and last of them due to the symmetry of function
P with respect to the central longitudinal and transverse planes of bearing on the line of intersection of
which lies the external force vector f, the second defines a function of pressure at the outlet of lubricant
duct.
# 100 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir A. Kodnyanko. Static Characteristics of Active Hydrostatic Two-Row Radial Bearing with Restriction...
Since the cavity 5 is stagnant, then at any point φj on the circle of pairing membranes and bearing
housing, flow rate is absent
Qm ( L1 , M j )
M j G
Rm
wP
( L1 , M )dM
H m3 ( L1 , M )
2S M j³G
wZ
(3)
0,
where Rm – dimensionless outer radius of rings.
It is also clear that in the absence of hydraulic resistance of lubricant flow in the channel 8 at the
entrance to the lubricating gaps same functions of pressure are equal
Pc ( L2 , M )
Pt (0, M )
Pm (0, M )
Pn (M ),
(4)
аnd balance of the dimensionless flow rates of working fluid in this section is determined by the
equation
Qn (M j ) Qm (0, M j ) Qt (0, M j ) Qc ( L2 , M j ) 0.
(5)
Here index n refers to the channel 8, L2 = L–L1 – length of c-area.
In mathematical modeling also suggested that the width of the movable ring is small compared
to its diameter. This assumption is justified by the fact that the narrow ends of the bearings have high
carrying load, lower radial and angular compliance. The presence of the narrow ends allowed neglect
the influence of circumferential overflows in the lubrication of t-area and m-area, i. e. assume that for
them ∂P/∂φ = 0. Simplifies equation (1) and its solution for t-area, taking into (2) can be found as a
function
L Z
(6)
Pt ( Z , M ) Pn (M ) 1
.
L1
Relative to coordinate Z linear is also the general solution of the simplified equation (1) for marea, and (3) is equal
wPm ( L1 , M )
wz
0.
In accordance with condition (4) solution of simplified equation (1) for m-area obtained as
Pm ( Z , M )
Pn (M ).
(7)
At c-area mentioned assumption can’t be applied, so solution of equation (1) founded by numerical
finite-difference method [7]. To do this, the intervals [0, L2] and [0, π] broke an even number n and m
equal parts of length ν = L2/n and τ = π/m, respectively.
Split points and function P values in finite-difference grid nodes are
Zi
iQ (i
0,1, 2,..., n), M j
jW ( j
0,1, 2,..., m), Pi j
P( Z i , M j ).
Equation (1) presented in a more convenient form
3 dH wP w 2 P w 2 P
H dM wM w 2M w 2 Z
0.
For interior points of area derivatives changed on symmetric finite-difference relations [6]
# 101 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir A. Kodnyanko. Static Characteristics of Active Hydrostatic Two-Row Radial Bearing with Restriction...
D j Pi j 1 Pi j 1 Pi j 1 2 Pi j Pi j 1 E Pi j1 2 Pi j Pi j1 0,
(8)
where
Dj
3WH Sin(M j )
§W ·
¨ ¸ .
©Q ¹
2
2 ¬ª1 H Cos (M j ) ¼º
, E
The desired pressure in nodes presented by the vectors
Pi
P , P , P ,..., P
0
i
1
i
2
i
m 1
i
, Pi m ,
i
0,1, 2,..., n and by means of (8) described in matrix-vector form
Pi 1 APi Pi 1
0, (i 1, 2,3,..., n 1),
(9)
where
A
2
0
0
ª K
«1 D
K
1 D1
0
1
«
«
0
1
D
K
1
D2
1
2
«
..
...
...
E « ...
« 0
0
0
0
«
0
0
0
¬ 0
...
...
0
0
0
0
...
0
0
...
...
...
... 1 D m 1 K
...
0
2
º
»
»
0 »
» ,K
... »
1 D m 1 »
»
K ¼
0
0
2(1 E ).
System of linear equations (9) is solved by matrix sweep method [7], using the recurrence
formulas
Pi 1
X i Pi Yi , (i 1, 2,3,..., n 1),
(10)
where Xi, Yi – sweep matrixes and vectors.
The system of equations (9) supplemented by vector equation
P0 4 P1 3P2
(11)
0,
which is a finite-difference analogue of the first boundary condition (2) for Z = 0.
In system of two linear matrix equations (11) and (10) for i = 1, after elimination of P2 found
P0
1 ·
§
¨ 2 E A ¸ P1 ,
2 ¹
©
(12)
where E – identity matrix.
Comparing (12), (10) for i = 1, found
X1
1 ·
§
¨ 2 E A ¸ , Y1
2 ¹
©
(13)
0.
After substituting (10) into (9) found the recurrence formulas
X i 1
A X i , Yi 1
1
X i 1Yi , (i 1, 2,..., n 1).
# 102 #
(14)
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir A. Kodnyanko. Static Characteristics of Active Hydrostatic Two-Row Radial Bearing with Restriction...
Next, performed a direct sweep by (13) and (14) and were found sweep matrixes X2, X3, …, Xn and
vectors Y2, Y3, …, Yn.
For perform of reverse sweep requires vector Pn. To determine it used the flow rate balance
equation (5) through the gaps in narrow sectors [φj – τ/2, φj + τ/2]. The general formula for determining
flow rate through a cross section Z is given by [6]
Qk , j ( Z )
R
k
2S
Mj W
2
³W H
Mj
3
k
(M )
wPk ( Z , M )
dM ,
wZ
(15)
2
where Rk – inner or outer radius of the ring (inner radius Rk = 1, outer radius Rk = Rm).
Since the integration step τ is a small quantity, instead of (15) used simplified formula
Qk , j ( Z )
wPk ( Z , M j )
W Rk 3
H k (M j )
.
2S
wZ
(16)
Dependences (6) and (7) provided a formula for local flow rates at the inlet of t-area
Qt , j (0)
W
2S L1
(17)
H t3, j Pnj .
So far as function (7) does not depend on the longitudinal coordinate, in accordance with (16)
local flow rates at the entrance of cavity 5
(18)
Qm, j (0) 0.
Local flow rates on the boundary Z = L2 determined by means of the finite-difference formula for
the derivative at the right endpoint [8]
wP( L2 , M j )
3Pnj 4 Pnj1 Pnj 2
.
2Q
wZ
Substituting (19) in (16), found
W
Qc , j ( L2 )
H 3j 3Pnj 4 Pnj1 Pnj 2 .
4SQ
(19)
(20)
In the absence of external load the pressure in channel (8) and at entrance to all areas does not
depend of circumferential coordinate and is equal to
Pnj
F
Const ,
j
0,1, 2,..., m 1, m ,
So full flow rate through the throttling slit is determined by the formula [6]
Qn
An 1 F .
(21)
In the coaxial arrangement of the movable elements in c-area Qc = 0 and in t-area in view of (6)
and (15) total flow rate through one end of the bearing is equal to
Qt 0
1
2S
2S
§ H t30 F ·
³0 ¨© L1 ¸¹dM
F
L1
H t30 .
# 103 #
(22)
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir A. Kodnyanko. Static Characteristics of Active Hydrostatic Two-Row Radial Bearing with Restriction...
Substituting founded flow rate dependences in equation (5), obtained a formula for calculating the
parameter of throttling slit
F
An 1 F L1
H t30 ,
F H t30
.
L1 1 F An
(23)
Taking this into account the local flow rates through the throttling slit are determined by the
formula
W An
(24)
Qn , j
1 Pnj .
2S
After substituting (17), (18) (20), (24) into equation (15) had obtained a system of equations
An 1 Pnj H t3, j
Pnj L1
(j
H 3j
2Q
3P
j
n
4 Pnj1 Pnj 2 0.
(25)
0,1, 2,..., m)
In matrix-vector form (25) is equal to
BPn 4 Pn 1 Pn 2
(26)
ɋ,
where B diagonal matrix, C – vector
B
diag{b0 , b1 , b2 ,..., bm }, C
§
H t3, j
3 ɫ j ¨ An ¨
L1
©
2Q
, bj
H 3j
ɫj
An (c0 , c1 , c2 ,..., cm ),
·
¸¸ .
¹
Equations (9), (10) for i = n – 1 and (26) are given a system of matrix equations
­ Pn APn 1 Pn 2 0,
°
® X n P n Yn P n 1 ,
° BP 4 P P
ɋ
n 1
n2
¯ n
for determination of unknown vector
Pn
B E DX n C DYn ,
1
4 E A.
D
Using (10) performed reverse sweep and founded vectors Pn-1, Pn-2,…, P1, P0.
The resulting solution allowed define load-carrying capacity of lubricating area by the general
formula [6]
Wk
2 Rk
S
S
a
0
0
³ CosM ³ P (Z , M )dZ dM ,
k
where a – length of area.
For t-area according to (6)
Wt
2
S
L1
CosM ³ P (M )
S³
n
0
0
L1 Z
dZ dM
L1
# 104 #
L1
S
S
³ P (M )CosM dM
n
0
L1
S
J1 ,
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir A. Kodnyanko. Static Characteristics of Active Hydrostatic Two-Row Radial Bearing with Restriction...
where
S
J1
³ P (M )CosM dM.
(27)
n
0
For m-area according to (7)
2 Rm
Wm
S
S
L1
0
0
³ CosM ³ Pn (M )dZ dM
2 Rm L1
S
J1 ,
Integral (27) calculated by mains of quadrature Simpson's rule [9]
J1
W
m
¦ k P Cos( jW ),
3
j
j n
(k j
1, 4, 2, 4, 2,..., 2, 4,1).
j 0
For c-area in the calculation of the load capacity used cubature Simpson’s rule [10]
Wɫ
2
S
S
L2
³ CosM ³ P(Z , M )dZ dM
0
0
2QW
9S
m
ª
n
¦ «¬k Cos( jW )¦ k P
j
j 0
i i
i 0
j
º
».
¼
Full load capacity of bearing calculated by formula
W
2 Wt Wc .
To determine the full volumetric flow rate through the bearing used formula
Q
4Qt
2
S
S
3
³ Ht
0
wPt
dM
wZ
S
2
H 3 P dM
S L1 ³0 t n
2W m
¦ k j H t3, j Pnj .
3S L1 j 0
The force balance equation of the ring, commits to the membrane-hydrostatic suspension small
radial displacement, based on the use of Hooke's law [11]. As mentioned above, the direction of ring
movement and external force vector are opposite. Therefore, the eccentricity εm has non-positive value
(εm ≤ 0) and is associated with the acting forces on the ring by relation
Fm Wm Wt
(28)
0,
where
Fm
Km / İm ,
– elastic force of membrane, Km = Const – dimensionless compliance of the membrane material.
Calculation procedure
In calculations were used the input dimensionless parameters: pressure coefficient χ∈ [0, 1],
lengths L and L1, ring radius Rm, clearances Ht0 and Hm0, compliance coefficient of membranes Km, n
and m numbers of segments for finite-difference grid. Main performance characteristics of interest in
the present study are load capacity W, flow rate Q, eccentricities ε and εt.
In calculation of dependences were varied eccentricity ε m with small step throughout
the range [–Hm0, 0]. Successively changing this parameter and using obvious connection of
eccentricities
# 105 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir A. Kodnyanko. Static Characteristics of Active Hydrostatic Two-Row Radial Bearing with Restriction...
İ
İ m İt
for each fixed εm by bisection numerical method [12] found the solution of nonlinear equation (28) in
form
F m Wm Wt
G (İ t )
(29)
0
for one unknown variable εt∈[a, b], where a = 0, b = Ht0. Process was stopped when b – a < 10 -5.
Simulation results
In Fig. 2 shows plots of function G(εt) of equation (29) for different values of eccentricity εm.
It is seen that the curves on the interval εt∈[0, Ht0] are continuous, have a unique point of
intersection with the horizontal axis (the root of the equation) or on this segment does not intersect it
at all (no root). The latter refers to specific cases of contact of rings working surfaces and the shaft, in
which the bearing is losing efficiency.
The effect on the characteristics of the circumferential lubrication overflows in interrow area
can be seen in Fig. 3, which shows the comparative dependences ε(W), obtained in [4] on the basis of
one-dimensional model (1-D model) and calculated out in this work by using two-dimensional model
lubricant flow (2-D model).
It is seen that for bearings with smooth working surfaces using a one-dimensional flow model of
interrow area gives overestimation of load capacity of 1.5 – 2 times.
Distortion of characteristics is particularly noticeable at moderate and high loads and for low
capacity Km of membranes, when bearing still has positive compliance K = ∂ε/∂W. With the increase
Km error in the calculation of compliance K in small eccentricities area observably decreases and at
bearing work on modes of zero and negative susceptibility (K ≤ 0) one-dimensional model gives fairly
satisfactory results.
G(İt)
E
0,08
İm = -0.3
0,06
0,04
İm = -0.6
0,02
0
-0,02
-0,04
-0,06
İm = -0.9
-0,08
-0,10
0
0,1
0,2
0,3
0,4
0,5
Fig. 2. Plots of function G(εt)
# 106 #
0,6
0,7
0,8
0,9
İFt
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir A. Kodnyanko. Static Characteristics of Active Hydrostatic Two-Row Radial Bearing with Restriction...
İE
2D
1D
0,8
Kİm = -0
0,6
İm4 =
0,4
0,2
0
İ12
m=
-0,2
İm8 =
-0,4
0
0,2
0,4
0,6
0,8
1,0
1,2
F
W
Fig. 3. Comparative dependences ε(W) for a smooth (2-D) and profiled (1-D) bearings for various values of membrane compliance Km, χ = 0.52, L = 1.5, L1= 0.3, Rm = 1.2, Ht0 = 1
For moderate and large eccentricities the error is large for any values of the parameter Km. An
analysis of the data for the various teachings of elongation L, at high loads, even for short bearing (L <
0.6) on the charts there is a noticeable difference curves.
Analysis of the calculated dependences shows that the decisive influence on the static characteristics
has a pressure parameter χ, which defines setting for input throttling slit resistance. It is known that for
conventional bearing function of eccentricity ε(χ), the load capacity W(χ) and compliance K(χ) are of
an extreme character.
In Fig. 3 shows the curves of ε(χ) for being studied bearing.
It is seen that they have the same character. The calculations showed that, regardless of
displacement of mobile elements minimum eccentricity is provided at approximately the same value
of the parameter χopt.
So for the graphs presented in Fig. 4, the smallest value of the eccentricity and compliance occurs
at χ = 0.52. The only parameter that significantly affects the optimal value of χ, is a clearance Ht0. Thus,
when Ht0 = 1.5 optimal χopt = 0.24, at Ht0 = 1.25 χopt = 0.42, at Ht0 = 0.8 χopt = 0.68.
Attention is drawn to the fact that an increase in Km significantly expands the range of perceived
bearing loads (Fig. 3). In Fig. 5 shows the curves that show the relationship with the capacity of the
membranes Km percentage T of increase the maximum load capacity W under maximum load ordinary
bearing (Km = 0).
From the Fig. 5 graphs show that the dependences T(Km) are extreme, i. e. for a fixed set of
parameters it can specify a quite definite value of Km, in which the bearing will have the widest range
of carrying loads.
# 107 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
İE
0,35
İm = -0.1
0,30
-0.2
0,25
0,20
-0.3
0,15
-0.4
0,10
0,05
-0.5
0
-0,05
-0,10
-0,15
-0,20
0
0,1
0,2
0,3
0,4
0,5
0,6
0,7
0,8
0,9
Ȥ
F
Fig. 4. Dependences ε (χ) for different values of the eccentricity ε m of rings L = 1.5, L1 = 0.3, Rm = 1.2, Ht0 = 1,
Km = 12
E
T
0.2
40
L1 = 0.1
30
0.3
20
10
0.4
0
0
5
10
15
20
25
Fig. 5. Dependences T(Km) for different length L1 of the rings, L = 1.5, Rm = 1.2, Ht0 = 1
KF m
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir A. Kodnyanko. Static Characteristics of Active Hydrostatic Two-Row Radial Bearing with Restriction...
Decreasing width L1 of rings increases the maximum value of this relationship. Thus, when L1 =
0.4 the best result occurs when Km = 7, for L1 = 0.3 Km = 9, at L1 = 0.2 Km = 14, with L1 = 0.1 Km = 25.
For Fig. 5 data range of load capacity will be 28 – 43%. The increase of relative bearing length also
enhances load range. So, when L = 2 T = 48%, with L = 3 T = 52%, with L ≥ 4 T = 55%.
Conclusion
This study results allow us to conclude that the correct calculation of static characteristics of
hydrostatic radial bearings with smooth cylindrical surfaces and output flow rate restrictors over the
entire range of operating loads can be carried out on the basis of two-dimensional model of lubricant flow
only. However, for small eccentricities characteristics of bearing with zero or negative susceptibility
can be calculated with sufficient accuracy by the simplified method, based on one-dimensional motion
of lubricant flow [4]. Bearings with zero and negative compliances have capacity that is 20 – 50% more
than conventional bearings of the same overall dimensions. The hydraulic resistance setting of input
throttling slits decisive influence on the optimal static characteristics of the bearing. In this case the
optimal values of resistance for conventional and active bearings are practically identical.
Final conclusion on the performance of bearings can be obtained on the basis of dynamic quality
study and experimental investigation of its performance indicators.
References
[1] Kodnyanko V.A. Shatokhin S.N. // Engineering Science. 1978. № 6. P. 90–95
[2] Kodnyanko V.A., Shatokhin S.N. // Engineering Science. 1980. № 6. P 108–112.
[3] Pat. 2370680 (2008). Russian Federation. IPC. 2009. F16C 32/06. Bull. № 29.
[4] Kodnyanko V. A. / / Journal of Siberian Federal University. Machinery and technology.
Krasnoyarsk, 2009. T. 3. № 3. P. 444–453.
[5] Kodnyanko V. A. / / Moscow, Moscow State Technical University Stankin. № 16. 2001. http://
magazine.stankin.ru/arch/n_16/index.shtml
[6] Constantinescu V.N. Gas lubrication. Moscow, Mashinostroenie. 1968. P. 709.
[7] Samarsky A. A., Gulin A.V. Numerical methods. Moscow, Nauka. 1989. P. 432.
[8] Demidovich B. P., Maron I. A., Shuvalova E. Z. Numerical Methods of Analysis / / Moscow,
Fizmatgiz. 1963. 660 P.
[9] Dwight G. Tables of integrals and other mathematical formulas. Moscow, Nauka. 1973. 228 P.
[10] Kahaner D., Mouler C., Nash S. Numerical Methods and Software. Springer-Verlag. 2001.
575 P.
[11]. Sedov L. I. Continuum Mechanics. Moscow, Nauka. Vol 1. 1970. 452 P.
[12]. Krylov V. I., Babkov V. V. Monasyrsky P. I. Computational methods. Moscow, Nauka. 1976.
312 P.
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir A. Kodnyanko. Static Characteristics of Active Hydrostatic Two-Row Radial Bearing with Restriction...
Статические характеристики активного
гидростатического двухрядного
радиального подшипника
с ограничением выходного потока смазки
В.А. Коднянко
Сибирский федеральный университет,
Россия 660041, Красноярск, пр. Свободный, 79
Рассмотрена конструкция активного гидростатического радиального подшипника с гладкими
цилиндрическими рабочими поверхностями и ограничителями выходного смазочного потока
смазки в виде подвижных колец с мембранным подвесом. Устройство в несколько раз менее
энергоемко по сравнению c известными устройствами с регуляторами расхода. Подшипник
обладает отрицательной и нулевой податливостью (бесконечной жесткостью), поэтому
может быть использован в металлорежущих станках для подавления негативного влияния
деформации упругой системы на точность обработки.
На основе двухмерной модели смазочного потока разработана математическая модель,
метод и методика расчета несущей способности и расхода смазки подшипника. Установлено,
что расчет статических характеристик подшипника во всем диапазоне действующих
нагрузок может быть корректно выполнен только на основе двухмерной модели. При малых
эксцентриситетах характеристики нулевой и отрицательной податливости могут быть
с удовлетворительной точностью рассчитаны по упрощенной методике, базирующейся
на одномерном движении смазочного потока. Подшипники нулевой и отрицательной
податливости обладают грузоподъемностью, которая на 20 – 50 % больше, чем у обычных
подшипников тех же габаритных размеров. Настройка гидравлического сопротивления
входной питающей щели решающим образом влияет на оптимальные статические
характеристики подшипника. Оптимальные значения сопротивления щели для обычных и
активных подшипников практически совпадают.
Ключевые слова: энергосберегающая, гидростатическая опора, гидростатический подшипник,
нулевая податливость, бесконечная жесткость, отрицательная податливость, гладкие
цилиндрические поверхности.
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Journal of Siberian Federal University. Engineering & Technologies 1 (2012 5) 111-121
~~~
УДК 336.763
Extreme Value Theory
and Peaks Over Threshold Model
in the Russian Stock Market
Vladimir O. Andreeva*, Sergey E. Tinykovb,
Oksana P. Ovchinnikovaa and Gennady P. Parahinc
а
Oryol Regional Academy of State Service
5a Pobedy st., Oryol, 302028 Russia
b
Branch of Siberian Federal University, Zheleznogorsk
12a Kirov st., Zheleznogorsk, Russia
c
TelekomStroyService, Ltd.
61 Tsarev Brod st., Zarechensky, 302528 Russia 1
Received 6.02.2012, received in revised form 13.02.2012, accepted 20.02.2012
Traditional research methods adopts normal distributions as a pattern of the stock market behavior.
This paper utilized POT model of extreme value theory, and GPD distribution which can give more
accurate description on tail distribution of financial returns/losses. EVT and POT techniques are
applied to a series of daily losses of the RTS index (RTSI) over a 15-year period (1995-2009), RTSI
is total index of 50 largest Russian stocks. The focus is on the use of proposed methods to asses tail
related risk providing a modeling tool for modern risk management.
Keywords: Extreme Value Theory, General Pareto Distribution, Peaks over Threshold, Tail
Distribution, Value at Risk.
Introduction
The study of extreme events has attracted the special attention in connection with the global
crisis of 2008–2009. The Russian stock market has been dramatic volatile over 15-year period
(from 38 points on 05.10.1998 to around 2487 points on 19.05.2008 and back to about 498 points
on 23.01.2009). In irregular fi nancial market, it is necessary to set up models and systems to
evaluate and control risks. In this paper we focus on the extreme behavior of fi nancial series,
unraveling the volatilities in the fi nancial markets has always been an decipherable mystery. One
of the purposes of this chapter is to test the validity of a popular risk management instrument:
Value-at-Risk estimator in Russian equity market, which is a widely adopted technique in the
developed countries for quantifying market risk. We have to deal with extreme events when a risk
takes values from the tails of its probability distribution. In the field of market risk management
it is a great concern the day by day determination of the Value-at-Risk (VaR) [1]. VaR is a high
*
1
Corresponding author E-mail address: oragsob@mail.ru
© Siberian Federal University. All rights reserved
# 111 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir O. Andreev, Sergey E. Tinykov… Extreme Value Theory and Peaks Over Threshold Model in the Russian...
quantile of the distribution of losses (for example the 95th percentile): VaR p=F-1(p), where F is the
loss cumulative distribution function and p the selected probability level. Traditional procedure
calculating VaR based on normal distribution has limitations. VaR model reflects that it is to asses
the possible maximum loss under regular market environments. Risk managers have become more
concerned with events occurring under extreme market conditions [2,3]. This paper argues that
extreme value theory (EVT) and POT (Peaks Over Threshold) model provide tools for estimating
measures of tail risk under irregular volatility in market. We consider a fully parametric model,
based on the GPD (Generalized Pareto Distribution), which can be easily estimated by maximum
likelihood method [4,5].
I. Theoretical framework of the extreme value approach
Extreme value theory is a powerful and fairly robust framework to study the tail behavior of a
distribution. There have been a number of extreme value studies in the finance literature in recent
years: quantile estimation using the extreme value theory [6]; the estimation of the tails of loss severity
distributions and the estimation of the quantile risk measures for financial time series using extreme
value theory [7,8]; overview the extreme value theory as a risk management tool [9]; potentials and
limitations of the extreme value theory [10,11]; an extensive overview of the extreme value theory
for risk managers [12]; the estimation of tail-related risk measures for heteroskedastic financial time
series [13]; comprehensive source of the extreme value theory to the finance and insurance literature
[14,15].
POT model and Generalized Pareto distribution
We use of Extreme Value Theory to model the tail returns and then show how our EVT estimates
are incorporated into the risk measures. Two main approaches are proposed in the literature [16]: the
Block Maxima (BM) and the Peaks-over-Threshold models (POT). The group of models for threshold
exceedances are more modern and powerful than the BM models [16], we focus on this approach and
its application to the losses on the RTSI stock index. We apply the parametric POT method based on
the Generalized Pareto distribution (GPD) to describe tail behaviour. Begin by assuming that market
losses represent the realizations x of a random variable X over an enough high threshold u. More
particularly, if X has the cumulative distribution function F(x), we are interested in the cumulative
distribution function Fu(x) of exceedances of X over a high threshold u, i.e. the conditional excess
distribution function is defined as:
Fu ( x)
P( X u d x | X ! u )
F ( x u ) F ( u )
1 F ( u )
.
(1)
As to the sufficient large u, EVT provides us with a powerful key result, which states for a large
class of underlying distributions F(x) [2]:
Fu ( x) | G[ , E ( x), u o f ,
where Generalized Pareto Distribution is defined by:
1 (1 [ x / E ) 1/ [ , if [ z 0
G[ , E ( x) {
1 exp( x / E ), if [ 0
# 112 #
(2)
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir O. Andreev, Sergey E. Tinykov… Extreme Value Theory and Peaks Over Threshold Model in the Russian...
[0,f), ɟɫɥɢ [ t 0
{[0,E /[ ), ɟɫɥɢ [ 0
where x 
.
GPD subsumes three other distributions under its parameterization [2]. So, when tail index ξ=0,
we obtain a Type 1 (exponentially declining) distribution. If ξ<0, we have a Type 2 (power declining).
For ξ>0, we obtain a Type 3 (constant declining) distribution. Given these three types of distribution,
one of our tasks in this paper will be to uncover which type best describes the extremes of stock returns
on the emerging Russian market.
Identification of GPD parameters
Let (X1,X2,…,X k(u)) be a sequence of iid random variables from an unknown distribution F, Xi>u.
Shape ξ and scale β>0 parameters are then defined on the threshold u [16]. These GPD parameters can
be determined by maximum likelihood (ML) methods. The log likelihood function of the GPD for ξ≠0
is:
n
k (u )(log E ) (1 1/ [ )¦ log(1 [ xi / E ) ,
L([ , E )
(3)
i 1
where xi satisfies the constraints specified for xi. If ξ=0, the log likelihood function is:
n
k (u ) ln( E ) E 1 ¦ xi .
L( E )
(4)
i 1
ML estimates are then found by maximizing the log-likelihood function using numeral optimization
methods. We can get these ξ and β estimates through solving simultaneous equations (ξ≠0):
wL([ , E )
1
w[
[
wL([ , E )
wE
2
k (u )
xi
1
k (u )
xi
¦ log(1 [ E ) ( [ 1)¦ E [ x
i 1
k (u )
E
i 1
k (u )
1
( 1)¦
[
i 1
[ xi
E E[ xi
2
0
i
0
Tail evaluation formula
Assuming that u is sufficiently high, by combining expressions (2) and (3) the distribution function
F(x) for exceedances can be written as:
F ( x) (1 F (u )) Fu ( x) F (u )
1 k (nu ) (1 [E ( x u )
1[
(5)
) .
Fu(x) is a GPD and F(u) is given by [n-k(u)]/n; n is the total number of observations, k(u) the
number of observations above the threshold u, ξ and β are the parameters of the GPD.
Estimating VaR
For a given probability p>F(u) and threshold u, the value-at-risk (VaR) is calculated by inverting
the tail estimation formula (5):
VaR p
u [E {[ k (nu ) (1 p)][ 1}.
# 113 #
(6)
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir O. Andreev, Sergey E. Tinykov… Extreme Value Theory and Peaks Over Threshold Model in the Russian...
Choosing threshold value
Choice of the threshold u is the important issue to deal with: u too high results in too few
exceedances and consequently high variance estimators. On the other hand, u too small provides
biased estimators and the approximation to a GPD could not be feasible. It is possible to choose an
asymptotically optimal threshold by a quantification of a bias versus variance trade-off.
Mean excess function
One suggestion which is of immediate use in practice is based on the linearity of the mean excess
function e(u) for the GPD. From [2] we know that for a random value X with a GPD distribution
function Gξ,β the mean excess function is:
E( X u | X ! u)
e(u )
E [ u
1[
,
E u[ ! 0, [ 1.
(7)
It suggests a graphical approach for choosing u: choose u >0 such that e(x) is approximately
linear for x¸u. Using plots to compare resulting estimates across a variety of u-values, due to the usual
presence of multiple choice of the threshold, is recommended.
Hill plot
Let X1>X2>…>X n be the order statistics of positive random variables iid. The Hill estimator of the
tail index ξ using k+1 order statistics is defined by [17]:
[ˆ
1
k
¦
k
j 1
(ln X j , n X k , n ) .
(8)
Hill plot is a good instrument to find the optimal threshold [18]. Over a specific range of
exceedances, the Hill plot may be under the stationary series, and the turning point is a good choice of
optimal threshold. We use the following intuitive ideas:
(1) The sequence of the turning point is less than ~n/10 [19].
(2) The Hill estimator in the turning point has a relative large deviation from the fitted stationary
straight line.
(3) The turning point is the last sequence of point that satisfies the two conditions stated above.
Empirical results
We consider a extreme value approach, working on the series of daily log losses (negative returns)
of the Russian RTSI Index over a period of fifteen years (1995-2009). The Russian Trading System
Index (RTSI) comprises of 50 of the largest stocks capturing 85% of the total market capitalization of
the Russian Trading System exchange. The data used in this paper are obtained from RTS web site [21].
The empiric study uses the series of log daily losses of the RTSI Index, containing 3 447 trading days
(closing prices). Fig.1 shows the plot of daily dynamics of RTSI index values, and log daily losses.
Table 1 shows the summary statistics for the series of log daily changes. This table shows that
kurtosis value is 9.7024 and skewness value is 0.3752. Relative value of Normal distribution is 3 and
0, respectively. So we can see empirical distribution of log daily losses and normal distribution is not
compatible.
# 114 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir O. Andreev, Sergey E. Tinykov… Extreme Value Theory and Peaks Over Threshold Model in the Russian...
2500
2000
1500
1000
500
0
1995
2000
2005
2010
a) daily dynamics of index values
0.25
0.2
0.15
0.1
0.05
0
-0.05
-0.1
-0.15
-0.2
-0.25
0
500
1000
1500
2000
2500
3000
3500
b) log daily losses
Fig. 1. RTSI Index – sample period 01.09.1995 – 30.06.2009 (closing values)
In addition to this, Jarqua-Bera statistic shows that law of log daily losses is obviously different
from normal distribution. The JB test statistics is defined as [10]:
JB
n
6
( Kurtosis 3)
[ STD
].
6 24
2
2
The JB statistic has approximately a chi-squared distribution, with two degrees of freedom.
The Jarqua-Bera test depends on skewness and kurtosis statistics. If the JB test statistic equals
zero, it means that the distribution has zero skewness and kurtosis is about equal 3, and so it can be
concluded that the normality assumption holds. Skewness values far from zero and kurtosis values
far from 3 lead to an increase in JB values. The test returns the logical value h = 1 if it rejects the
null hypothesis at the p<0.05 significance level, and h =0 if it cannot. We have for data of Table 1:
JB value=6532.8, p~0, h=1. It means that we can reject the hypothesis that the distribution of daily
losses is normal.
# 115 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir O. Andreev, Sergey E. Tinykov… Extreme Value Theory and Peaks Over Threshold Model in the Russian...
Table 1. Summary statistics for daily losses in RTS
mean
min
max
-0.0007
-0.2020
0.2120
std
skewness
kurtosis
0.0289
0.3752
9.7024
variance
JB test
n
-0.0008
6532.8
h=1,p<0.001
3447
In Fig.2 we represent the Hill graph, which plots the Hill estimator of ξ, versus the k upper order
statistics (and threshold u, respectively). We select the last area to k~0.1*3447~350, where the Hill
estimator is more stable.
The mean excess function (7) allows to establish the behavior of the distribution tails [23]: we
choose threshold u looking at the linear shape (with positive slope) of the graph (Fig.3). Considering
Hill plot and the mean excess function, we choose u=0.0334 (the number of observation exceeding
threshold u is equal k=294).
The results of ML estimation of the GPD parameters (on chosen threshold u=0.0334) are ξ=0.1492
and β = 0.0206:
Maximum Likelihood (ML) estimates of ξ,β:
out =
par_ests: [0.1492 0.0206]
funval: -803.6979
par_ses: [0.0688 0.0018]
threshold: 0.0334
data: [1x294 double]
p_less_thresh: 0.9675
QQ-plot graph makes us able to evaluate the goodness of fit of the empirical series to a parametric
GPD model (Fig.4) [24]. Notice that a concave departure from the straight line in the QQ-plot (Fig.4a)
is an indication of heavy tailed distribution, whereas a convex departure is an indication of a thin tail.
After we get estimates ξ,β, use them in (5), get the formula for of tail evaluation:
p
F ( x)
1
k (u )
(1 0.1492 *
x 0.0334
n
k (u )
)
1
0.1492
0.0206
294, n=3447, k(u)/n=8.53%
Employ the result in (6), get the VaR formula on GPD model:
VaR p
0.0334 0.0206
0.1492
{[
3447
294
(1 p )]0.1492 1} .
In Table 2 we report 95%, 99%, 99.5%, 99,9% Value-at-Risk estimates of three different VaR
estimation methods. The performance of the different VaR estimation methods can be evaluated by
comparing the estimates with the actual losses observed, in particular by computing (and testing)
the number of VaR violations. VaR approaches based on the assumption of normal distribution are
# 116 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Hillplot
0.5
0.4
xi
0.3
0.2
0.1
0
0
50
100
150
200
250
Order Statistics
300
350
400
a) Hill estimator versus k upper order statistics (k≤350)
Hillplot
0.55
0.5
xi
0.45
0.4
0.35
0.3
200
250
300
Order Statistics
350
400
b) Hill estimator versus k upper order statistics (k≤350) (plot zooming)
Fig. 2. Hill estimator versus k upper order statistics (probability level p=0.95)
0.22
0.2
0.18
Mean Excess
0.16
0.14
0.12
0.1
0.08
0.06
0.04
0.02
-0.25
Fig. 3. Mean excess function
-0.2
-0.15
-0.1
-0.05
0
Threshold
0.05
0.1
0.15
0.2
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir O. Andreev, Sergey E. Tinykov… Extreme Value Theory and Peaks Over Threshold Model in the Russian...
9
8
Exponential Quantiles
7
6
5
4
3
2
1
0
-0.2
-0.15
-0.1
-0.05
0
0.05
Ordered Data
0.1
0.15
0.2
а) QQ-plot: empirical vs exponential distribution
16
14
GPD Quantiles; xi=0.1492
12
10
8
6
4
2
0
0
0.05
0.1
Ordered Data
0.15
0.2
b) QQ-plot: empirical vs GPD distribution ( ξ=0.1492, β=0.0206, u=0.0334)
Fig. 4. QQ-plot versus GPD distribution and exponential distribution
Table 2. VaR estimation for daily RTSI losses:one day horizon
VaR approach
p=0.950
p=0.975
p=0.990
p=0.995
p=0.999
Normal model
0.0394
0.0451
0.0520
0.0565
0.0663
Historical
simulation
GPD model
0.0452
0.0607
0.0849
0.1083
0.1771
0.0499
0.0611
0.0856
0.1062
0.1620
definitely to underestimate high percentiles, while estimates based on historical simulation face
with the problem of out of sample performance. The extreme value approach on GPD model seems
appropriate and easy to implement.
# 118 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir O. Andreev, Sergey E. Tinykov… Extreme Value Theory and Peaks Over Threshold Model in the Russian...
Conclusions
Since last century, volatility of international financial system is getting severe. A stable financial
system is so desirable. Therefore, risk management has aroused growing attention. As a measurement
of market risk, VaR has been widely used in risk management. However, derivation between VaR
estimation of normal hypothesis and abnormal distribution of practical benefit rate of financial
always cause the bigger error in estimation. Aiming at this problem, through GPD model which fits
tail distribution of financial products more accurately, this paper recalculates VaR by POT method.
Compared with traditional method of risk study, this paper has made some progress in research
approach and philosophy and more applicable in practice, which has been demonstrated by example of
the Russian market analysis.
We use software systems: EVIM [25,26], MATLAB [27] and LOGOS-EVT, developed by authors
of this paper.
References
[1] Jorison, P. (2001) Value at Risk: the New Benchmark for Managing Financial Risk, 2nd ed,
McGraw-Hill
[2] Embrechts P., Kluppelberg C. and Mikosch T. (1997) Modelling extremal events: for insurance
and finance, Springer, Berlin.
[3] Fernandez, V., 2005. Risk management under extreme events. International Review of Financial
Analysis 14, 113-148.
[4] Kluppelberg, C. (2004) Risk management with extreme value theory. In: Finkenstaedt, B. and
Rootzen, H. (Eds) Extreme Values in Finance, Telecommunication and the Environment, pp.101-168.
Chapman and Hall/CRC, Boca Raton.
[5] Gilli, M., Kellezi, E., 2006. An application of extreme value theory for measuring financial
risk. Computational Economics 27, 207-228.
[6] de Haan, L., Jansen, D. W., Koedijk, K. G., and de Vries, C. G. (1994). Safety first portfolio
selection, extreme value theory and long run asset risks. In J. Galambos, J. Lechner and E. Simiu, eds.,
Extreme value Theory and Applications, Kluwer, Dordrecht, pages 471–488.
[7] McNeil, A. J. (1997). Estimating the tails of loss severity distributions using extreme value
theory. ASTIN Bulletin, 27, 1117–137.
[8] McNeil, A. J. (1998). Calculating quantile risk measures for financial time series using extreme
value theory. Manuscript, Department of Mathematics, ETH, Swiss Federal Technical University.
[9] Embrechts, P., Resnick, S., and Samorodnitsky, G. (1998). Extreme value theory as a
risk management tool. Manuscript, Department of Mathematics, ETH, Swiss Federal Technical
University.
[10] Embrechts, P. (1999). Extreme value theory in finance and insurance. Manuscript, Department
of Mathematics, ETH, Swiss Federal Technical University.
[11] Embrechts, P. (2000a). Extreme value theory: Potentials and limitations as an integrated
risk management tool. Manuscript, Department of Mathematics, ETH, Swiss Federal Technical
University.
[12] McNeil, A. J. (1999). Extreme value theory for risk managers. Internal Modeling and CAD
II, Risk Books, pages 93–113.
# 119 #
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir O. Andreev, Sergey E. Tinykov… Extreme Value Theory and Peaks Over Threshold Model in the Russian...
[13] McNeil, A. J. and Frey, R. (2000). Estimation of tail-related risk measures for heteroscedastic
financial time series: An extreme value approach. Journal of Empirical Finance, 7, 271–300.
[14] Embrechts, P., Kluppelberg, C., and Mikosch, C. (1997). Modeling Extremal Events for
Insurance and Finance. Springer, Berlin.
[15] Embrechts, P. 2000. Extremes and Integrated Risk Management. Risk Books and UBS
Warburg, London.
[16] McNeil A.J., Frey R. and Embrechts P. Quantitative risk management: Concepts, techniques
and tools. 2005. Princeton University Press.
[17] B. M. Hill, A simple general approach to inference about the tail of a distribution, Annals
Statistics 3 (1975) 1163–1173.
[18] Reiss, R.-D.and Thomas, M.(1997). Statistical Analysis of Extreme Values. Birkhauser, Basel
nd
(2 ed. 2001).
[19] Zhou, Chunyang, Wu, Chongfeng, Liu, Hailong and Liu, Fubing,A New Method to Choose
the Threshold in the Pot Model. Available at SSRN.
[20] Loretan, M., Phillips, P., 1994, Testing the Covariance Stationarity of Heavy tailed Time
Series, Journal of Empirical Finance 1, 211-248.
[21] www.rts.ru
[22] Brys G, Hubert M, Struyf A. A Robustification of the Jarqua-Bera Test of Normality.
COMPSTAT 2004 Symposium, Section: Robustness, 2004.
[23] Erik Brodin, Claudia Klueppelberg Extreme Value Theory in Finance. In: Everitt, B. and
Meiinick, E. (Eds.) Encyclopedia of Quantitive Risk Assessment.
[24] Coles, S., 2001. An Introduction to Statistical Modeling of Extreme Values. Springer-Verlag,
London.
[25] Gencay, R., Selcuk, F., and Ulugaulyagci, A. (2003). EVIM: a software package for extreme
value analysis in Matlab. Studies in Nonlinear Dynamics and Econometrics, 5:213-239.
[26] www.sfu.ca/~rgencay/
[27] www.mathworks.com
Copyright ОАО «ЦКБ «БИБКОМ» & ООО «Aгентство Kнига-Cервис»
Vladimir O. Andreev, Sergey E. Tinykov… Extreme Value Theory and Peaks Over Threshold Model in the Russian...
Применение метода превышений порога
и теории экстремальных значений
для моделирования рисков
на российском фондовом рынке
В.О. Андреева, С.Е. Тиняковб,
О.П. Овчинниковаа, Г.П. Парахинв
а
Орловская региональная
академия государственной службы
Россия 302028, Орел, ул. Победы, 5а
б
Железногорский филиал СФУ
Россия, Железногорск, ул. Кирова, 12а
в
ТелекомСтройСервис
Россия 302528, Орловская обл.,
п. Зареченский, ул. Царев Брод, 61
При моделировании поведения фондового рынка обычно используется нормальное
распределение Гаусса. Однако, метод превышений порога (РОТ) теории экстремальных
значений (EVT), обобщенное распределение Парето (GPD) позволяют более точно описывать
финансовые доходности (потери) от операций с ценными бумагами, особенно на хвостах
вероятностных распределений. В данной работе описывается применение предложенных
методов к моделированию и анализу потерь на основе индекса РТС, представляющего собой
усредненную цену акций 50 крупнейших российских компаний, за 15-летний период (1995 - 2009
г.). Особое внимание уделено использованию предложенных методов для оценки экстремального
рыночного риска на хвостах распределений, что позволяет получить современный инструмент
моделирования для системы управления рисками.
Ключевые слова: теория экстремальных значений, обобщенное распределение Парето, метод
превышений порогового значения, вероятностное распределение на хвосте, стоимость под
риском.
Документ
Категория
Наука
Просмотров
501
Размер файла
3 178 Кб
Теги
170, журнал, университета, техника, технология, 2012, сер, сибирской, федеральное
1/--страниц
Пожаловаться на содержимое документа