# 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 г.). Особое внимание уделено использованию предложенных методов для оценки экстремального рыночного риска на хвостах распределений, что позволяет получить современный инструмент моделирования для системы управления рисками. Ключевые слова: теория экстремальных значений, обобщенное распределение Парето, метод превышений порогового значения, вероятностное распределение на хвосте, стоимость под риском.

1/--страниц