r. jecL et aL. DOUBLE DIFFUSIVE NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER WITH THE BOUNDARY DOMAIN INTEGRAL METHOD v. vukadin THE DEVELOPMENT AND ANALYSIS OF 3D GEOLOGICAL/ GEOMECHANICAL MODELS FOR THE KOZJAK PUMPED/STORAGE HYDROELECTRIC POWER PLANT j. Kortnik OPTIMIZATION OF THE HIGH SAFETY PILLARS FOR THE UNDERGROUND EXCAVATION OF NATURAL STONE BLOCKS g. vizintin et aL. THE DEVELOPMENT OF A "DRIVE-IN" FILTERS DEWATERING SYSTEM IN THE VELENJE COAL MINE USING FINITE-ELEMENT MODELLING a. strukeLj et aL. PREDICTION OF THE PILE BEHAVIOUR UNDER DYNAMIC LOADING USING EMBEDDED STRAIN SENSOR TECHNOLOGY ■H 1— ■H 9 LO OO ■H N S S I fiCTfl GGOTGCHNICñ SLOVGNICñ Founders ISSN: 1854-0171 ustanovLteLjL Univerza v Mariboru, Fakulteta za gradbeništvo University of Maribor, Faculty of Civil Engineering Univerza v Ljubljani, Fakulteta za gradbeništvo in geodezijo University of Ljubljana, Faculty of Civil and Geodetic Engineering ■ i :. Univerza v Ljubljani, Naravoslovnotehniška fakulteta University of Ljubljana, Faculty of Natural Sciences and Engineering Slovensko geotehniško društvo Slovenian Geotechnical Society Društvo za podzemne in geotehniške konstrukcije Society for Underground and Geotechnical Constructions urednLškL odbor edLtorLaL Board LzdajateLj IpubLLsher Univerza v Mariboru, Fakulteta za gradbeništvo University of Maribor, Faculty of Civil Engineering odgovornL urednLk edLtor-Ln-chLef Ludvik Trauner Univerza v Mariboru urednLka co-edLtors Stanislav Škrabl Univerza v Mariboru Bojan Žlender Univerza v Mariboru TehnLčna urednLka Desk edLtors Bojana Dolinar Univerza v Mariboru Borut Macuh Univerza v Mariboru Lektor proof-Reader Paul McGuiness NakLada cLrčuLatLon 500 izvodov - issues TLsk prLnt Tercia tisk d.o.o. Ptuj Revija redno izhaja dvakrat letno. Članki v reviji so recen-zirani s strani priznanih mednarodnih strokovnjakov. Baze podatkov v katerih je revija indeksirana: SCIE - Science Citation Index Expanded, JCR - Journal Citation Reports / Science Edition, ICONDA - The international Construction database, GeoRef Pri financiranju revije sodeluje Javna agencija za knjigo Republike Slovenije. Darinka Battelino Università degli Studi di Trieste József Farkas Budapesti Müszaki és Gazdaságtudományi Egyetem Theodoros Hatzigogos Aristotle University of Thessaloniki Rolf Katzenbach Technische Universität Darmstadt Zlatko Langof Univerzitet u Sarajevu Jakob Likar Univerza v Ljubljani Janko Logar Univerza v Ljubljani Bojan Majes Univerza v Ljubljani Milan Maksimovic Univerzitet u Beogradu Borut Petkovšek Zavod za gradbeništvo Slovenije Mihael Ribičič Univerza v Ljubljani César Sagaseta Universidad de Cantabria Stephan Semprich Technische Universität Graz Abdul-Hamid Soubra Université de Nantes Ivan Vaniček Ceské vysoké ucení technické v Praze Franjo Veric Sveučilište u Zagrebu Address NasLov urednLštva ACTA GEOTECHNICA SLOVENICA Univerza v Mariboru, Fakulteta za gradbeništvo Smetanova ulica 17 2000 Maribor Slovenija Telefon / Telephone: +386 (0)2 22 94 300 Faks / Fax: +386 (0)2 25 24 179 E-pošta / E-mail: ags@uni-mb.si spLetnL nasLov http:IIwww. fg.uni-mb.siIjournal-ags I web Address The journal is published twice a year. Papers are peer reviewed by renowned international experts. Indexation data bases of the journal: SCIE - Science Citation Index Expanded, JCR - Journal Citation Reports / Science Edition, ICONDA - The international Construction database, GeoRef Financially supported also by Slovenian Book Agency. VSEBINA Ludvik Trauner UVODNIK Renata jecL Ln drugI NARAVNA KONVeKCIJA V POROZNEM SLOJU ZARADI DVOJNE DIFUZIJE Z ROBNO OBMOČNO INTEGRALSKO METODO Vladimir vukadin RAZVOJ IN ANALIZA 3D GEOLOŠKO/GEOMEHANSKIH MODELOV ZA ČRPALNO ELEKTRARNO KOZJAK Jože KortnLk OPTIMIRANJE IN SPREMLJAVA VISOKIH VARNOSTNIH STEBROV PRI PODZEMNEM PRIDOBIVANJU BLOKOV NARAVNEGA KAMNA Goran VLžLntLn Ln drugi RAZVOJ VTISNIH FILTROV ODVODNJEVAL-NEGA SISTEMA PREMOGOVNIKA VELENJE S POMOČJO METODE KONČNIH ELEMENTOV 3- Andrej Štrukelj Ln drugL NAPOVED OBNAŠANJA PILOTA MED DINAMIČNIM OBREMENILNIM TESTOM Z UPORABO TEHNOLOGIJE VGRAJENIH SENZORJEV DEFORMACIJ gj 10. ŠUKLJETOVI DNEVI NAVODILA AVTORJEM CONTENTS Ludvik Trauner EDITORIAL Renata Jecl et al. DOUBLE DIFFUSIVE NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER WITH THE BOUNDARY DOMAIN INTEGRAL METHOD Vladimir vukadin THE DEVELOPMENT AND ANALYSIS OF 3D geological/geomechanical MODELS FOR THE KOZJAK PUMPED-STORAGE HYDROELECTRIC POWER PLANT E jože Kortnik OPTIMIZATION OF THE HIGH SAFETY PILLARS FOR THE UNDERGROUND EXCAVATION OF NATURAL STONE BLOCKS E Goran Vizintin et at. THE DEVELOPMENT OF A "DRIVE-IN" FILTERS DEWATERING SYSTEM IN THE VELENJE COAL MINE USING FINITE-ELEMENT MODELLING [65 nndrej Štrukelj et at. PREDICTION OF THE PILE BEHAVIOUR UNDER DYNAMIC LOADING USING EMBEDDED STRAIN SENSOR TECHNOLOGY 10th LUJO SUKLJE MEMORIAL DAY INSTRUCTIONS FOR AUTHORS L9 E UVODNIK Veseli smo, da se je število avtorjev, ki želijo objaviti svoje prispevke v reviji Acta Geotechnica Slovenica, v zadnjem letu zelo povečalo. V uredništvu smo se zato odločili, da povečamo tudi število člankov v posamezni številki revije. Tokrat smo izbrali pet zelo zanimivih prispevkov, ki obravnavajo predvsem inovativne pristope pri reševanju praktičnih problemov s širokega področja geotehnike. V prvem prispevku avtorjev R. Jecl, J. Kramer in L. Škerget je predstavljena robno območna integralska metoda, ki predstavlja eno izmed numeričnih metod za reševanje problemov povezanih z gibanjem tekočine v porozni snovi. Dinamika tekočin v porozni snovi je pomembno področje v znanosti in številnih inženirskih panogah, prenosni pojavi v porozni snovi so predmet številnih raziskav predvsem zaradi mnogih aplikacij na področjih gradbeništva, procesnega strojništva, kemije, ekologije. V članku je podana numerična rešitev vodilnih enačb s katerimi je opisano gibanje viskozne nestisljive tekočine v porozni snovi z uporabo ustrezno modificirane metode robnih elementov. V članku avtorja V. Vukadina so prikazani rezultati geoloških, geofizikalnih in geotehničnih preiskav na območju predvidene gradnje črpalne elektrarne ČHE Kozjak (Slovenija) in njihova uporaba za razvoj in izdelavo 3D geološko/ geomehanskih modelov. S pomočjo prostorske analize je mogoče lažje razumeti odnose med posameznimi enotami in elementi, ter hitro ugotoviti najbolj kritična območja s stališča geotehnične konstrukcije in jih tudi stabilnostno preveriti. Avtor tretjega prispevka J. Kortnik predstavlja postopke za načrtovanje, optimiranje in spremljavo visokih varnostnih stebrov pri podzemnem pridobivanju blokov naravnega kamna. v kamnolomih pisanega apnenca Hotavlje I in Lipica II v Sloveniji. Avtorji četrtega članka G. Vižintin, M. Veselič, A. Bombač, E. Dervarič, J. Likar in Ž. Vukelič obravnavajo razvoj vtisnih filtrov odvodnjevalnega sistema premogovnika Velenje s pomočjo metode končnih elementov. Podzemna dela v premogovniku Velenje ves čas ogroža prisotnost podzemne vode. Nevarnost vdorov vode se bo še posebej povečala po letih 2012 in 2017, ko bodo zaradi posedanja terena, ki nastaja zaradi odkopavanja in subsidence, opustili površinska črpališča. Zaradi tega bo morala biti na novo postavljena strategija odvodnjevanja. V članku, ki so ga prispevali avtorji A. Štrukelj, H. Vrecl-Kojc, M. Pšunder in L. Trauner so podani rezultati dinamičnega obremenilnega preizkusa na testnem pilotu v območju avtoceste Slivnica - Hajdina v bližini Maribora. V pilot je bilo vgrajenih devet posebej izdelanih senzorjev za meritve specifičnih deformacij, ki so bili razporejeni na enakomernih razdaljah dveh metrov vzdolž osi pilota. Omenjeni senzorji so omogočili spremljanje deformacijskega stanja v betonu pilota med samim preizkusom ter, na osnovi izmerjenih rezultatov, ocenitev normalne napetosti vzdolž osi pilota in trenje vzdolž plašča pilota ter normalne napetosti pod nogo pilota. V reviji je priloženo še povabilo za udeležbo na jubilejni 10. Šukljetov dan, ki bo 25. septembra na Brdu pri Kranju. Kot pobudnik Šukljetovih dnevov in organizator le teh vam z veseljem sporočam, da bo tokratno srečanje precej praznično obarvano, saj so na ta jubilej povabljeni tudi vsi dosedanji vabljeni predavatelji iz tujine (M. Jamiolkowski, G. Sanglerat, D.D. Potts, H. Brandl, R. Katzenbach, S. Semprich, S. Leroueil, D. Žnidarčic, A.S. Nossan). Posebno zanimivo pa bo tudi letošnje vabljeno predavanje strokovnjaka svetovnega slovesa Prof. Dr. Chandra S. Desai-a z naslovom »Constitutive Modelling and Computer Methods in Geotechnical Engineering«. Ludvik Trauner EDITORIAL We are very glad to see that the number of authors wishing to publish their contributions in our journal Acta Geotech-nica Slovenica has grown in the last year. Therefore, our editorial board has decided to enlarge the number of articles in each issue. This time we have chosen five, extremely interesting, contributions that provide us with new, innovative approaches to the solutions of practical problems from the broad field of geotechnics. In the first contribution, written by R. Jecl, J. Kramer and L. Škerget, the boundary-domain integral method, which provides one of the numerical methods for solving the transport phenomena in porous media, is presented. Fluid dynamics in porous media is an important topic in many branches of engineering and science; transferable phenomena in porous media are the object of numerous research activities, mostly a result of the variety of applications in the fields of civil engineering, process mechanical engineering, chemistry and ecology. In the article, a numerical solution for the governing equations is described, which is then used for a description of the flow of a viscous incompressible fluid in a porous medium, using an appropriate extension of the boundary-element method. In the second article, the author, V. Vukadin, provides us with the results of geological, geophysical and geotechnical research in the area of the future construction of the Kozjak pumped-storage hydroelectric power plant (Slovenia) as well as their usage in the development and production of 3D geological/geomechanical models. With the help of a space analysis, the relations between single units and elements are easier to understand and it is also possible to establish quickly the most critical areas from the point of view of geotechnical constructions and to test their stability. The author of the third article, J. Kortnik, presents procedures for the planning, optimization and monitoring of high safety pillars for the underground excavation of natural stone blocks at the Hotavlje I colourful limestone quarry and at the Lipica II quarry in Slovenia. ® . # The authors of the fourth article, G. Vižintin, M. Veselič, A. Bombač, E. Dervarič, J. Likar and Ž. Vukelič, deal with the development of a "drive-in" filter dewatering system with the help of finite elements in the Velenje coal mine. The underground works in the Velenje coal mine are constantly threatened by the presence of groundwater. In particular, after 2012 and 2017, the risk of water inrushes will increase, when a series of surface-drilled wells will be abandoned due to excavation works and mining-subsidence effects. Therefore, the dewatering schemes need to be completely reviewed. In the article by A. Štrukelj, H. Vrecl-Kojc, M. Pšunder and L. Trauner, the results of a dynamic loading test of a pile on the highway section Slivnica-Hajdina near Maribor are provided. In order to measure specific deformations, nine specially made strain sensors were built in and distributed evenly at a distance of two metres along the whole pile axis. These sensors enabled monitoring of the deformation behaviour in the concrete pile body during the performance of the test as well as, based on the measured results, an estimation of the shear stresses along the pile axis, the friction along the pile shaft and the normal stresses below the pile toe. Attached to the Journal please find the invitation to attend the 10th jubilee of the symposium Šukljetov dan (Šuklje's Day) which will be held on the 25th September at Brdo near Kranj. As initiator and organizer of the symposium I am glad to inform you that the September meeting will be quite a festive one as it will gather all past invited lecturers from abroad ( M. Jamiolkowski, G. Sanglerat, D.D. Potts, H. Brandl, R. Katzenbach, S. Semprich, S. Leroueil, D. Žnidarčic, A.S. Nossan). This year's invited lecture will be delivered by the world-renowned scientist, Professor Dr. Chandra S. Desai, who will speak about Constitutive Modeling and Computer Methods in Geotechnical Engineering. I hope that the lecture will be of great interest to you. Ludvik Trauner Editor-in-chief NARAVNA KONVEKCIJA V POROZNEM SLOJU ZARADI DVOJNE DIFUZIJE Z ROBNO OBMOČNO INTEGRALSKO METODO_ RENATA JECL, JANJA KRAMER Ln LEOPOLD ŠKERGET o avtorjih Renata Jecl Univerza v Mariboru, Fakulteta za gradbenišvo Smetanova ulica 17, 2000 Maribor, Slovenija E-pošta: renata.jecl@uni-mb.si Janja Kramer Leopold Škerget Univerza v Mariboru, Univerza v Mariboru, Fakulteta za gradbenišvo Fakulteta za strojništvo Smetanova ulica 17, 2000 Maribor, Slovenija Smetanova ulica 17, 2000 Maribor, Slovenija E-pošta: janja.kramer@uni-mb.si E-pošta: leo@uni-mb.si izvleček V prispevku je predstavljena robno območna integralska metoda, ki predstavlja eno izmed numeričnih metod za reševanje problemov povezanih z gibanjem tekočine v porozni snovi. Dinamika tekočin v porozni snovi je pomembno področje v znanosti in številnih inženirskih panogah, prenosni pojavi v porozni snovi so predmet številnih raziskav predvsem zaradi mnogih aplikacij na področjih gradbeništva, procesnega strojništva, kemije, ekologije. V članku je podana numerična rešitev vodilnih enačb s katerimi je opisano gibanje viskozne nestisljive tekočine v porozni snovi z uporabo ustrezno modificirane metode robnih elementov. Rešen je konkreten primer naravne konvekcije zaradi dvojne difuzije v horizontalnem poroznem sloju, ki je popolnoma zasičen s tekočino. Prenosni pojavi so procesi, ki opisujejo kako se določena veličina (masa, gibalna količina, temperature, snov) prenaša v porozni snovi. Naravna konvekcija je eden od najpogosteje obravnavanih prenosnih pojavov v porozni snovi, nastane pa kot posledica razlik v temperaturi in/ali koncentraciji snovi. Te razlike privedejo do sprememb v gostoti in posledično nastanka vzgonskih sil (termičnih in/ali snovskih), ki povzročijo gibanje. Zasledimo lahko veliko študij konvektivnega toka, ki ga povzročata tako termična kot snovska vzgonska sila - ta pojav imenujemo tudi naravna konvekcija zaradi dvojne difuzije. Praktične primere s tega področja srečamo pri problemih geotermalne energije, pri toplotni izolaciji in potovanju kontaminantov v podtalnici ali v deponijah. V pričujočem prispevku bo najprej predstavljen problem prenosa toplote in snovi v vertikalni smeri, kar pomeni, da obravnavamo horizontalni sloj v katerem nastopi konvektivni tok, ko je Rayleighovo število nad kritično vrednostjo, drug primer pa je problem prenosa toplote v horizontalni smeri in snovi v vertikalni smeri. Oba primera sta rešena na osnovi splošnega numeričnega algoritma, ki temelji na metodi robnih elementov. Osnova matematičnega modela so Navier-Stokesove enačbe za čisto tekočino, ki so zapisane na makroskopskem nivoju preko reprezentativnega elementarnega poroznega volumna (upoštevajoč dejstvo, da je samo del tega volumna na voljo za tok tekočine), kar vodi do modificiranih Navier-Stokesovih enačb za porozno snov. Omeniti velja, da v modelu uporabljamo kompletno, s postopkom povprečenja izpeljano, gibalno enačbo in ne le Darcy-jevo enačbo, kar je praksa pri večini študij, ki jih najdemo v literaturi. Ker vodilne enačbe predstavljajo vezan sistem močno nelinearnih difuzivno konvektivnih parcialnih diferencialnih enačb je potrebno za reševanje uporabiti robno območno integralsko metodo, ki je ena izmed izpeljank klasične metode robnih elementov. Enačbe se najprej transformirajo z uporabo hitrostno-vrtinčne formulacije, v nadaljevanju pa se z uporabo ustreznih Greenovih funkcij oziroma z metodo utežnih ostankov zapišejo v splošni integralski obliki, kasneje pa še v diskretni obliki za celotno računsko območje. Osnova numeričnega modela je programski paket BEEAS, Boundary Element Engineering Analysing System, ki je dopolnjen za reševanje prenosnih pojavov v porozni snovi. Izpeljan numerični algoritem je testiran na zgoraj omenjenih primerih, grafično in tabelarično so predstavljeni in komentirani rezultati za različne vodilne parametre (Rayleighovo število, Darcyjevo število, Lewisovo število, vzgonski koeficient) in podana je primerjava z nekaterimi objavljenimi numeričnimi študijami. Zaključki, kot so dobro ujemanje rezultatov izračuna z nekaterimi objavami v literaturi in tudi fizikalno smiseln potek tokovnih, temperaturnih kot tudi koncentracijskih polj, potrjujejo pravilnost zastavljenega matematičnega modela in uporabljene numerične sheme ter dokazujejo, da je robno območna integralska metoda učinkovita alternativa ostalim, v svetovnem prostoru že uveljavljenim numeričnim metodam za reševanje zapletenih prenosnih pojavov v porozni snovi. Ključne besede porozna snov, robno območna integralska metoda, naravna konvekcija zaradi dvojne difuzije, Darcy-Brinkmanova enačba DOUBLE DIFFUSIVE NATURAL CONVECTION IN A HORIZONTAL POROUS LAYER WITH THE BOUNDARY DOMAIN INTEGRAL METHOD RENATA JECL, JANJA KRAMER and LEOPOLD ŠIKERGET About the authors Renata Jecl University of Maribor, Faculty of Civil Engineering Smetanova ulica 17, 2000 Maribor, Slovenia E-mail: renata.jecl@uni-mb.si Janja Kramer University of Maribor, Faculty of Civil Engineering Smetanova ulica 17, 2000 Maribor, Slovenia E-mail: janja.kramer@uni-mb.si Leopold Škerget University of Maribor, Faculty of Mechanical Engineering Smetanova ulica 17, 2000 Maribor, Slovenia E-mail: leo@uni-mb.si Abstract We present the boundary-domain integral method, one of the numerical methods for solving the transport phenomena in porous media. The results for the case of double diffusive natural convection in a porous horizontal layer, which is fully saturated with an incompressible fluid, are obtained. Modified Navier-Stokes equations were used to describe the fluid motion in porous media in the form of conservation laws for mass, momentum, energy and species. Several results for different cases of double diffusive natural convection in a porous horizontal layer are presented and compared with some published studies in which calculations with other numerical methods were performed. Keywords porous media, boundary domain integral method, double diffusive natural convection, Darcy-Brinkman equation 1 INTRODUCTION The main purpose of this article is to present the Boundary Domain Integral Method as a numerical method for solving problems encountered with the flow through porous media. Fluid dynamics in porous media is an important topic in many branches of engineering and science, as well as in many fields of practical interest. It plays a fundamental role in ground-water hydrology and soil mechanics, as well as petroleum, sanitary and chemical engineering. New computational methods and techniques have allowed us to model and simulate various transport phenomena in porous media, which means our understanding of them is improving constantly. The aim of this work is to obtain a numerical solution for the governing equations describing the flow of a viscous incompressible fluid in a porous medium, using an appropriate extension of the boundary-element method. The numerical scheme was tested on a natural convection problem within a porous square cavity as well as in the porous layer, where different temperatures and concentration values are applied on the horizontal walls. The results for the different governing parameters (the Rayleigh number, the Darcy number, the buoyancy ratio and the Lewis number) are presented, commented on and compared with previously published work. 1.1 TRANSPORT PHENOMENA IN POROUS MEDIA Fluid transport phenomena in porous media refer to those processes related to the transport of fluid momentum, mass and heat, through the given porous media. These processes, which are encountered in many different branches of science and technology, are commonly the subject of theoretical treatments that are based on methods traditionally developed in classical fluid dynamics. Natural convection is the most commonly studied transport phenomena in porous media and also a process to which we pay full attention. There have been several reported studies dealing with natural convection due to thermal buoyancy forces, mainly because of its importance in industrial and technological applications, such as geothermal energy, fibrous insulation, etc. Less attention, however, has been dedicated to the so-called double diffusive problems, where density gradients occur due to the effects of combined temperature and compositional buoyancy. There are some important engineering applications for this phenomenon, such as the transport of moisture in fibrous insulations or grain-storage insulations and the dispersion of contaminants through water-saturated soil. Three main configurations are studied when dealing with double diffusive natural convection in an enclosure filled with porous media [1], [2]: - temperature and species concentration or their gradi ents are imposed horizontally along the enclosure, either aiding or opposing each other, - temperature and species concentration or their gradients are imposed vertically, either aiding or opposing each other, - temperature or its gradient is imposed vertically and species concentration or its gradient is imposed horizontally or vice versa. The present analysis focuses on the use of the Boundary Element Method for solving problems of fluid flows in porous media driven by coupled thermal and solutal buoyancy forces. The Darcy-Brinkman formulation is used for modeling fluid flow in porous media, thus enabling a satisfactory non-slip boundary condition on those surfaces that bound the porous media. 1.2 BOUNDARY ELEMENT METHOD FOR FLUID DYNAMICS IN POROUS MEDIA Fluid flows in porous media have been studied both experimentally and theoretically over recent decades. Different numerical methods have been introduced to obtain the solutions for some transport phenomena, e.g., the finite-difference method (FDM), the finite element method (FEM), the finite volume method (FVM), as well as the boundary element method (BEM). ). The main comparative advantage of the BEM, the application of which requires the given partial differential equation to be mathematically transformed into the equivalent integral equation representation, which is later to be discretized, over the discrete approximative methods, is demonstrated in cases where this procedure results in boundary integral equations only. This turns out to be possible only for potential problems, e.g., inviscid fluid flow, heat conduction, etc. In general, this procedure results in boundary-domain integral equations and, therefore, several techniques have been developed to extend the classical BEM. The dual reciprocity boundary element method (DRBEM) represents one of the possibilities for transforming domain integrals into a finite series of boundary integrals, see for example [3] and [4]. The key point of the DRBEM is an approximation of the field in the domain by a set of global approximation functions and the subsequent representation of the domain integrals of these global functions by boundary integrals. The discretization of the domain is only represented by grid points (poles of global approximation functions) in contrast to FDM meshes. However, the discretization of the geometry and the fields on the boundary is piecewise polygonal, which gives this method greater flexibility compared to the FDM methods in coping with boundary quantities. In the DRBEM all the calculations reduce to the evaluation of the boundary integrals only. Another more recent extension of the BEM is the so-called boundary domain integral method (BDIM), see [5], [6], [7] and [8]. Here, the integral equations are given in terms of variables on the integration boundaries, as well as within the integration domain. This situation arises when we are dealing with strongly nonlinear problems that have prevailing domain-based effects, for example, diffusion-convection problems. Navier-Stokes equations are commonly used as a framework for the solution of such problems, since they provide a mathematical model of the physical conservation laws of mass, momentum and energy. The velocity-vorticity formulation of these equations results in a computational decoupling of the kinematics and kinetics of the fluid motion from the pressure computation, see [9]. Since the pressure does not appear explicitly in the field functions' conservation equations, the difficulty associated with the computation of the boundary pressure values is avoided. The main advantage of the BDIM, compared to the classical domain-type numerical techniques, is that it offers an effective way of dealing with boundary conditions on the solid walls when solving the vorticity equation. Namely, the boundary vorticity in the BDIM is computed directly from the kinematic part of the computation and not through the use of some approximate formulae. One of the few drawbacks of the BDIM is the considerable CPU time and memory requirements, but they can be drastically reduced by the use of a subdomain technique, see [10]. Convection-dominated fluid flows suffer from numerical instabilities. In domain-type numerical techniques upwinding schemes of different orders are used to suppress such instabilities, while in the BDIM the problem can be avoided by the use of Green's functions of the appropriate linear differential operators, which results in a very stable and accurate numerical description of the coupled diffusion-convective problems. There are no oscillations in the numerical solutions, which would have to be eliminated by using some artificial techniques, e.g., upwinding, artificial viscosity, as is the case of other approximation methods. 2 GOVERNING EQUATIONS Due to the general complexity of the fluid transport process in a porous medium, our work is based on a simplified mathematical model in which it is assumed that: - the solid phase is homogeneous, non-deformable, and does not interact chemically with respect to the fluid, - the fluid is single phase and Newtonian; its density does not depend on pressure variation, but only on the variation of the temperature, - the two average temperatures, Ts for the solid phase and Tf for the fluid phase, are assumed to be identical - the porous media is in thermodynamic equilibrium, which means it is described by a single equation for the average temperature T = Ts = T} , - no heat sources or sinks exist in the fluid-saturated porous media; the thermal radiation and viscous dissipation are negligible, - the porosity and permeability are assumed to be constant throughout the whole cavity, - the density of the fluid depends on the temperature and concentration variations and can be described by p = p0[1 -[3t(T— T0) — ftC(C-C0) ,where the subscript 0 refers to a reference state, f3T is the volumetric thermal expansion coefficient ¡3T = —1/ p [dp/ 5T] , and ¡3C is the volumetric concentration expansion coefficient 3c =—V p [dp/ C . The transport phenomenon in the porous media is described using modified Navier-Stokes equations in the form of conservation laws for mass, momentum, energy and species. The equations are written at the macroscopic level, derived by averaging the microscopic equations for the pure fluid over the porous representative elementary volume and expressed by continuity, momentum, energy and species equations [11]: dvi dx, 0 dC dv,C d dt dx. dx. dC i- dx. D- (4) — = 0, (1) 1 dv 1 dVjVi --L + 0 dt 0 dx. 1 dp v d =---- + Fg,--v, +- p„ dx, K dx. v 2 d r i dvT d -0 c/ +(1 - 0), ] T + Cf = - dT dx. , (2) (3) The parameters used above are as follows: v, volume-averaged velocity, x, the i-th coordinate, 0 porosity, t time, p density, v kinematic viscosity, dp/dxi the pressure gradient, g, gravity, K permeability of the porous media and è,, the strain rate tensor è.. = l/2(dv,./dxj +dvj/dxl andfinally F is the normalized density difference function and is given as F = (P - Po)/Po =[[ (T-T) + I3C (C - C0) .Furthermore, cf = (pc) and cs = (pc)s are the isobaric specific heat per unit volume for the fluid and solid phases, respectively, T is the temperature, Xe is the effective thermal conductivity of the porous media given as Xe = 0\f + (l — 0)Xs, where Xf and As are the thermal conductivities for the fluid and solid phases, respectively. In the final equation C stands for the concentration, and D for the mass diffusivity. Equation (2) is known as the Darcy-Brinkman equation and consists of two viscous terms: the Darcy viscous term (third on the r.h.s.) and the Brinkman viscous term, which is analogous to the Laplacian term in the Navier-Stokes equations for a pure fluid (fourth on the r.h.s.). A non-slip boundary condition on a surface that bounds the porous media is satisfied by using the additional Brinkman term. There are several situations where it is convenient to use the Brinkman equation, e.g., when one wishes to compare the flows in porous media with those in a pure fluid or to investigate the convective flows in the context of solidification, where permeability and porosity are variables in space and time. If the included parameter K (the permeability) tends towards infinity ( K ) the Brinkman equation transforms into the classical Navier-Stokes equation for a pure fluid. On the other hand, if the permeability tends to zero ( K ^ 0 ), the Brinkman term becomes negligible and the equation reduces to a classical Darcy equation [1]. A convective flow in a horizontal porous layer is possible above the critical Rayleigh number. In the case of double-diffusive convection, where the density differences are the result of combined temperature and concentration gradients, the critical Rayleigh number is a function of the Darcy number Da, the Lewis number Le and the buoyancy coefficient N [1]. In vertical cavities maintained to horizontal temperature and concentration gradients the flow in the cavity is always unicellular. In the case of a horizontal porous layer, where the temperature and concentration differences are imposed on the horizontal walls, the flow structure becomes multi-cellular and is also called a Rayleigh-Benard flow structure [2]. 3 NUMERICAL METHOD The Boundary Domain Integral Method (BDIM), an extension of the classical Boundary Element Method (BEM), is used in order to solve the general set of equations. The discretization of the surface and the domain is required since the boundary and domain integrals are represented in the obtained set of integral equations. The above given equations (1), (2), (3), (4) should be modified in order to use the BDIM. Firstly, the kinematic viscosity in the momentum equation is partitioned into constant and variable parts, such as v = v + V, so the Brinkman term is also divided into two parts, and the equation becomes: dv ' dv ' ,v ' 1 dp vé -L H--J—L =---— + Fg¡--— v \ dt dx, p0 dx¡ K d2 dx,dx, —^ - (5) do¡ ' ~oT ' do' Vj ~oX~ do' '¡jkSk dF dx, dv V vé ' dx, K ° JL\v doÛ+dL dx, I dx, I dx, , (10) where w \ is the modified vorticity w \ = and f' is any contribution arising as a r es ult of non-linear material properties, given as f.. = û(Vx è ). The mathematical description of the transport phenomena in fluid flow is completed by providing suitable Dirichlet and Neumann boundary conditions [12], as well as some initial conditions for the energy transport and species transport equations (see Fig. 1): - dT _ T = T on the boundary r ; d— n. = q on the boundary r2 ; T = T0 in the domain Q, (11) where the term vj is now the modified velocity v' = v, /4 . In the same way as the kinematic viscosity, the thermal diffusivity ap , which is defined as ap = \/Cj , and the mass diffusivity D are divided into the constant and variable parts ap = ap + ap , D = D + D . By including the expression for the heat-capacity ratio a = 4 + (1 — 4)cs / cf in the energy equation, the formulations (3) and (4) can be rewritten as: a 8T_ é dt dC dt ~ dv1T dx,. dv 'C d 2T dx,. dxjdxj D d2c é dxjdxj _d_ dx,. dx, aP dT é dx, D dC dx, (6) (7) - dC _ C = C on the boundary r ; -n. = q on the _ dx. 1 boundary r2 ; C = C0 in the domain fi , (12) In the present analysis, a two-dimensional problem is considered, therefore all subsequent equations will be written for the case of planar geometry. The linear elliptical Laplace differential operator can be used for the velocity equation (9): L[] = d 2( dx. dx. (13) and the following relationship for the kinematics can be obtained: In the next step the above-stated governing equations are transformed by the use of the velocity-vorticity formulation (VVF), and consequently the computational scheme is partitioned into its kinematic and kinetic parts [6]. The vorticity vector o>t is introduced, which represents the curl of the velocity field: ,=* o where e jjk is the unit permutation tensor. The kinematic part is represented by the ellipticalvelocity vector equation: d2 v ' L [v ''i+b¡=äx-öx;+b¡=0 - (14) where b i stands for the pseudo-body source term. An integral representation of the velocity vector can be formulated by using the Green theorems for scalar functions or the weighting residuals technique, rendering the following vector integral formulation: c(C)v ',(£) + f v ', q'dr= f^-i-udT+ f bu dQ , (15) i i dn d2v \ dxjdxj i'k du\ dx, = 0 , (9) and the kinetics is governed by the vorticity, energy and species transport equation. The vorticity transport equation is obtained as a curl of the Brinkman momentum equation (5): where £ is the source point, u is the elliptical Laplace fundamental solution and q is its normal derivative, e.g., q = du*/dn . The fundamental solution u* is given by the expression: u = — ln 2^ r(£,s) where r is the vector from the source point £ to the reference field point s. Equating the pseudo-body force as: », = , if . (IT) we obtain the following integral formulation: c(Ov\(?)+Jv, qd r = ¡dtu'd r + e« Sd7u'd ° -(18) Q 9Xj The derivative of the vorticity in the pseudo-body source terms can be eliminated by using the Gauss divergence theorem. The integral representation of the kinematics is now [13]: c(£ rn) x v 9)+n«) xf v' qd r = , (24) ' = n(£)x f (q* x n)x v 'dr + n(£)x f u X q'd0 r o in order to obtain an appropriate non-singular implicit system of equations for the unknown boundary vorticity or tangential velocity component values to the boundary [13]. The formulations for the vorticity, temperature and concentration can generally be written as a non-homogeneous elliptical diffusion-convective equation [8]: d2u dv'u dxj dxj dxj At -JL + » = 0 , (25) C(.0V\(?) + J v\qd r = J-'-u'd r + et] J u' n]u'd r—et] J u' q* r r dn r a ,(19) d a By using these expressions for the vorticity definition, unit tangent, and normal vector, e.g., dv j/dn = dv j/dx^n^, io' = e v dv i /dxt = dvy jdx — dvxjdy, n = (nx, ny) and t = (—ny, nx) for i, j = 1,2 and applying the continuity equation (1), the following relationship can be derived: an+e,un=—e, <») The boundary integrals on the right-hand side of equation (19) can be rewritten as: c(0v\(0 + Jv\q'dT = —ev f^t-udT — e9 Ju'qjdQ. ,(21) r r dt n With further mathematical reformulations and by applying the Gauss theorem, the resulting integral representation can be written as follows: c(0v\(0 + fv\qdT = ejfv'jq*dr — efffw'qjdn ,(22) r r n where q* is the tangential derivative of the fundamental solution q* = du" /dt, or in the form of the elliptical integral vector formulation: c(? )V ) + J V' q'd r = J (q x n)xv' d r + J u X q *d Q .(23) The tangential form of equation (22) can be written as follows: where u is taken as the vorticity u ', the temperature T and the concentration C, respectively, p is defined by considering the conservation laws and constitutive hypothesis always divided as p = P + P , and bt stands for the pseudo-body source term. Since the fundamental solution exists only for steady diffusion-convective Partial Differential Equations (PDEs) with constant coefficients, the velocity field is decomposed into an average constant vector vj ' and a variable vector p ', such that v\ = v V + V \. Thus, the following integral formulation can be obtained: c(0u(0 + pJuq dr = J(pq — uv'n)u*dr + JbudQ .(26) where u* is now the fundamental solution of the steady diffusion-convective PDE with a first-order reaction term [5], in the form of: 1 u = - 2nP K o(^r )exp 2p (27) du * du Furthermore, q = — and q = — represents the dn dn normal derivative of the fundamental solution and a suitable field function. The parameter ^ is defined as ¡ = — \2 2 p +koL = P — \2 2 p + f . where v'2 = v'. 3 = 1/p At, K0 is the modified Bessel function of the second kind of order 0 , and r is the magnitude of the vector from the source to the reference point, i.e., r = \xi(0 — xi. The following integral representations for the vorticity, temperature and concentration kinetics are obtained according to equation (26): a a c(Z)w'(Z) + f w'Q*dr =1 f(q — w'v'n + ejgjFnj + f]n] ))dr + r v r + V f(w 'vj — e jgjF -y qj — fj) Qjd n+V f ^w U' d n + — f wF .u*dn 7At J F—1 vAt> , (28) c(H)T(0 + fTQ*dr = d tt^q — Tv'n J an J I m r I aaP —=~ f \~irqj—Tv a m ■ aaP r Qjd n u d r — p n j ' aapAt n f Tf—U d n , (29) c(Z )C (Z)+f c q ' d r = D f\D q—Cv'n D■ u dr— — m f| -qt — Cv' Qjdf CF U D j I d, j j ^j DA t j F—1 D- D At- d n , (30) where Q* = dU*/dn represents a normal and Q * = dU jdxj a space derivative of the modified elliptical diffusion-convective solution U* defined as U = vu in the momentum equation, U* = apjdu* in the energy equation and U* = D/du* in the species equation. Furthermore, q = du/dn represents a normal and q. = dujdxj a space derivative of the field function (vorticity, temperature, concentration). For an approximate numerical solution, the integral equations are written in a discretized manner, where the integrals over the boundary and the domain are approximated by the sum of the integrals over all the boundary elements and internal cells, respectively. The variation of the field functions within each boundary element or internal cell is approximated by the use of the appropriate interpolation polynomials [14]. After applying the discretized integral equations to the boundary and internal nodes, the implicit matrix systems can be obtained, firstly for the kinematic computational part as: [H]{v\} = e j [H,][v — e j [Dj]{w'}. (31) Furthermore, we obtain the implicit systems for the vorticity kinetics: -1 [Dj ] [w v — ej — y q; — f} +1 [B]j vdw'} + [B]K—.} ,(32) then for the heat energy kinetics [H][T } = A- [G]J ^q — Tv' | —— \DJ ]J ^ — Tv 1 Jl 1 aap Jj d j aaP 1 j J[ d m ,(33) + -J- [B]{T }p_1 and finally for the species kinetics: [H][C}=D [o]{ D q—Cv'j—D j D q; — cy}+Da IB][C}f—i .(34) In the above equations the matrices [H], [Ht], [G], [Dj] and [5] are the influence matrices composed of those integrals that took over the individual boundary elements and internal cells. This system of discretized equations is solved by coupling these kinetic and kinematic equations and considering the corresponding boundary and initial conditions. Since the implicit set of equations is written simultaneously for all the boundary and internal nodes, this procedure results in a very large and fully-populated system matrix, influenced by diffusion and convection. The consequence of this approach is a very stable and accurate numerical scheme with substantial computer time and memory demands. The subdomain technique is used to improve the economics of the computation, where the entire computational domain is partitioned into subdomains, to which the same described numerical procedure can be applied. The n final system of equations for the entire domain is then obtained by adding the sets of equations for each subdomain and considering the compatibility and equilibrium conditions between their interfaces, resulting in a sparser matrix system, suitable for solving using iterative techniques. In the present case, each subdomain consists of four discontinuous 3-node quadratic boundary elements and one 9-node corner continuous quadratic internal cell [6]. 3.1 SOLUTION PROCEDURE The obtained set of discretized equations (31)-(34) is solved by coupling the kinetic and kinematic equations. In this way the solenoidality of the velocity field for an arbitrary vorticity distribution is assured. The following numerical algorithm has to be performed in order to obtain a final solution: (1) Start with some initial values for the vorticity distribution. (2) Kinematic computational part: - solve the implicit set for the boundary vorticity values (equation (31)), - determine the new domain velocity values, (3) Energy kinetic computation part: - solve the implicit set for the boundary and domain values (equation (33)), (4) Species kinetic computation part: - solve the implicit set for the boundary and domain values (equation (34), (5) Vorticity kinetic computational part: - solve the implicit set for the unknown boundary vorticity flux and internal domain vorticity values (equation (32)), (6) Relax new values and check the convergence. If the convergence criterion is satisfied, then stop; otherwise go to step 2. 4 TEST EXAMPLES - RESULTS AND DISCUSSION The convective flow in a horizontal porous layer is possible above the critical Rayleigh number. In the case of double-diffusive convection, where the density differences are a result of combined temperature and concentration gradients, the critical Rayleigh number is a function of the Darcy number Da, the Lewis number Le and the buoyancy coefficient N [14]. In vertical cavities maintained at horizontal temperature and concentration gradients, the flow in the cavity is always unicellular. In the case of a horizontal porous layer, where the temperature and concentration differences are imposed on the horizontal walls, the flow structure becomes multi-cellular and is also called a Rayleigh-Benard flow structure [15]. Most of the studies regarding double-diffusive convection or thermohaline convection (the case where the constituent is salt) in a horizontal porous layer are focused on the problem of convective instability. There are many studies dealing with the onset of convection on the basis of the linear stability theory [16], [17] or nonlinear perturbation theory [18]. In these studies, the critical Rayleigh numbers for the onset of convective flows are predicted. The theoretical and numerical study of heat and mass transfer affected by a high Rayleigh number Benard convection in a porous layer heated from below is obtained in [19]. The numerical results and a scale analysis of the flow in a porous medium are presented, where the buoyancy effect is due entirely to temperature gradients. Some further numerical results for a double-diffusive convection in a horizontal porous layer with two opposing buoyancy sources can be found in [20]. The influence of the governing parameters (the Rayleigh number, the Lewis number, and the buoyancy ratio) on the overall heat and mass transfer is discussed for the case of a square cavity (aspect ratio equal to 1). Double-diffusive convection in a horizontal layer with some numerical results is also discussed in [21]. In this study the critical values of the Rayleigh numbers for the onset of convective motion are predicted on the basis of a nonlinear parallel flow approximation. All the above-mentioned numerical results are obtained on the basis of the Darcy flow model, which is more convenient for porous media with a low permeability. The Brinkman extended Darcy model, which accounts for friction due to macroscopic shear is more appropriate when describing fluid flows in the porous matrix, when the inertia effects are negligible. It was used in [14] to investigate the onset and development of double-diffusive convection in a horizontal porous layer with uniform heat and mass fluxes specified at the horizontal boundaries. The obtained analytical solution is compared to some numerical results obtained for different values of Ra, N, Le and Da. Only a few recently published studies deal with the cross relation between the temperature and the solutal gradient (the temperature or concentration are imposed horizontally and the concentration or temperature gradient imposed vertically) [22], [23], [24]. In those studies it is shown that the competition between the thermal and solutal force produces very complex flow patterns. In our case the system under consideration is a horizontal layer of width L and height H, filled with homogenous, nondeformable porous media, which is fully saturated with a Newtonian fluid. In the first case, the horizontal walls are subjected to different temperature and concentration values (TB, CB at the bottom boundary and Tn, Cn at the upper boundary), while the vertical walls are adiabatic and impermeable (Figure 1). In addition, the second case - the example of cross gradients of the temperature and concentration - is presented in Figure 4, where the horizontal walls are maintained at different concentrations (CB at the bottom boundary and Cn at the upper boundary) and the vertical walls at different temperatures (TL at the left boundary and TR at the right boundary). The fluid, saturating the porous media, is modeled as a Boussinesq incompressible fluid, where the density depends only on the temperature and concentration variations P = Po(l-Pt(T-T>)-(3C(C-C0)) . f and ¡3C are the volumetric thermal and concentration expansion coefficients. The numerical results for the heat and mass transfer induced by the double-diffusive natural convection in a horizontal porous layer subjected to vertical gradients of temperature and concentration are presented. The governing parameters of the problem are: - the porosity $ , - the Darcy number Da = K / H2, - the aspect ratio A = L / H, - the modified (porous) thermal Rayleigh number Ra = Kg PT AT H / av , where fT is the thermal expansion coefficient, a the thermal diffusivity and v viscosity, - the Lewis number Le = a / D , - the buoyancy ratio N = f3C AC / f3T AT , where K is the permeability, D and H are the width and height of the layer, respectively, g is the acceleration due to gravity, AT and AC are the temperature and concentration differences between the upper and lower boundaries, a is the thermal diffusivity and D is the mass diffusivity. For the aspect ratio A = 1, a non-uniform computational mesh 20x20 was used with a ratio between the longest and shortest elements of r = 6 , and for A = 2 and A = 4 20x10 subdomains were used. The computational mesh was denser near the solid walls, where large gradients in the temperature, concentration and velocity profiles were expected. Time-steps ranging from At = 1016 s (steady state) to At = 10-4 s (time dependent) were employed, and the convergence criteria (representing the difference between two progressive values in each time step) were determined as s = 5 x 10-6 for all cases. It should be noted that in the case of N = 0 the buoyancy effect is due entirely to temperature gradients. The mass transfer in this case is due to the temperature field and the concentration differences between the horizontal boundaries. In the case of positive values for the buoyancy ratio (N > 0), the thermal and solutal buoyancy forces aid each other (aiding convection) and for negative values of the buoyancy ratio (N < 0) the solutal and thermal effects have the opposite tendencies (opposing convection). The obtained numerical model was tested for different values of the governing parameters. The results for total heat and mass transfer through the horizontal layer are given by the values of the Nusselt (Nu) and Sherwood numbers (Sh). Furthermore, some simulations of the temperature and concentration fields are presented. Tu = 0 , Cu=0 ,v, = 0 H Figure 1. Geometry of the first problem and the boundary conditions. The validation of the code was accomplished by a comparison with some published numerical experiments. However, it should be noted that there are not many published studies giving numerical results for the described problem. Table 1 presents some results for A=1, Da = 10-5 and a different Rayleigh number, Lewis number and buoyancy ratio. The values of the overall heat and mass transfer are compared to the published results, where the numerical calculations based on the Darcy model are obtained [20]. The first result is for the case of Ra = 600, Le = 1 and N = 0 , which means that only the thermal buoyancy force is present. The overall heat and mass transfers, which are represented by Nu an Sh, are identical. The other two cases are for the Rayleigh number 100, the buoyancy ratio 0.2 and the Lewis numbers 10 and 30. In this case, both the thermal and solutal buoyancy forces are present and they aid each other. The values of the Sherwood numbers are now higher than those of the Nusselt numbers, which is a result of a higher Lewis number. The presented results are in agreement with the published ones. Table 2 presents the influence of the Darcy number, where the governing parameters are as follows: an aspect ratio A = 1 , the Rayleigh number Ra = 100, the Lewis number Le = 10 and the buoyancy coefficient N = 0.2. It is evident that with a decrease of the Darcy number the values of the Nusselt and Sherwood numbers increase. The Darcy number describes the influence of the additional Brinkman viscous term in the momentum equation. At higher values of the Darcy number, the effect of viscous forces is significant and generally slows down the convective motion. With a decrease of the Darcy number Da ^ 0 the viscous effect becomes smaller, so the convective motion becomes stronger, which can be observed from the presented results. In Table 3 the results for the case of the opposing fluxes of heat and mass (N < 0) and the parameters A = 1 , Ra = 100, Le = 1, Le = 3 and N = -0.1 are presented. In the case where 10-3 < Da < 10-1 the values of the Nusselt and Sherwood numbers equal 1.0, which means the heat and solute transfer are governed by diffusion. The Rayleigh number, which depends on the Darcy number, is not high enough for the onset of convective motion, which begins in this from values of Da < 10 4. In the case of Le = 1 the values of the Nusselt and Sherwood numbers are identical. The Sherwood number increases with the Lewis number; however, on the other hand, the Nusselt number is almost independent of the Lewis number. Table 1. Comparison of results with numerical experiments reported in the literature. Nu Sh This study Literature [20] This study Literature [20] Ra=600, Le=1, N=0 7.01 6.6 7.01 - Ra=100, Le=10, N=0.2 2.76 2.4 8.84 - Ra=100, Le=30, N=0.2_2.50_2.5_14.80_15 Table 2. Nu and Sh numbers for different values of Da and A = 1 , Ra = 100, Le = 10, N = 0.2 Da 10-1 10-2 10-3 10-4 10-5 10-6 Nu 1.00 1.53 2.26 2.63 2.76 2.79 Sh 1.00 3.93 6.13 7.93 8.84 9.10 Table 3. Nu and Sh numbers for different values of Da and A = 1 , Ra = 100, N = -0.1 Da 10-1 10-2 10-3 10-4 10-5 10-6 Nu 1.00 1.00 1.00 2.30 2.40 2.42 Le = 1 Sh 1.00 1.00 1.00 2.30 240 2.42 Le = 3 Nu 1.00 1.00 1.00 2.32 2.42 2.44 Sh 1.00 1.00 1.00 4.05 4.30 4.36 Table 4. Nu and Sh numbers for different values of Da and A = 4, Ra = 300, Le = 0.1, N = -2 Da 10-1 10-2 10-3 10-4 10-5 Nu 1.00 2.05 2.82 3.12 3.50 Sh 1.00 1.02 1.04 1.05 1.06 Table 4 presents the results for a shallow layer with the governing parameters as follows: the aspect ratio A = 4, the Rayleigh number Ra = 300, the Lewis number Le = 0.1 and the buoyancy ratio N = -2. In this case again the thermal and solutal buoyancy forces oppose each other. Because of the higher value of N, the combined buoyancy is dominated by the buoyancy due to concentration gradients. From the table it is evident that with any decrease in the Darcy number the value of the Nusselt number increases. In the cases where the Lewis number decreases Le ^ 0, the values of the Sherwood number tend to unity (Sh ^ 1), which implies that the mass transfer is dominated by diffusion. The same conclusions are also published in [14]. A direct comparison of the results is impossible, because of the different governing parameters (porosity, aspect ratio) in both studies. The graphical simulations of the velocity, temperature and concentration fields are presented in Figure 2 for the square geometry (A = 1) with the parameters Ra = 100, Da = 10-4, Le = 10 and N = 0.2 in Figure 3 for the case where A = 2 for the same parameters. The flow structure consists of a unicellular circulation filling up the entire layer. The temperature and concentration field present the classical stratified structures of the natural convective flows. The temperature and concentration boundaries are observed at the top and bottom walls. The governing parameters for the second case are A = 2, Ra =100, Da = 10-3, Le = 10, and N = 1. The flow in the horizontal layer where A = 2 becomes multi-cellular (with two cells) as it is obvious from Figure 3. The flow consists of rising hot fluid in the centre of the layer and colder fluid sinking along the vertical walls. In the centre of the domain, a higher solute concentration is found than along the adiabatic and impermeable side walls. Thin temperature and composition boundary layers are evident at the top and bottom walls. Figure 3. Streamlines, isotherms and isoconcentrations in a horizontal layer for A = 2, Ra = 100, Da = 10-3, Le = 10, N=1. The next example presents the convection under cross gradients of temperature and concentration. The geometry and boundary conditions are shown in Figure (4). Figure 2. Streamlines, isotherms and isoconcentrations for A = 1, Ra = 100, Da = 10-4, Le = 10, N=0.2 . CB=1 f=0 ' dy Vi = 0 H v, =0 Figure 4. Geometry of the second problem and the boundary conditions. In this case the temperature differences are imposed on vertical walls and the species concentration differences on horizontal walls. Because a differentially heated cavity has a tendency to mix the fluid, while the solutal effect has a tendency to prevent such a mixing, it is expected that the competition between the solutal and thermal forces will produce complex flow patterns. The graphical simulations for the streamlines, temperature and concentration fields are presented in Figures 5 and 6 for the parameters A = 2, Ra =100, Da = 10-3 and Le = 10 and the values of the buoyancy coefficient from N = 0 to N = -4 . In the case when N = 0 the flow remains unicellular, while with any decrease in the buoyancy coefficient the flow forms two slow cells. The temperature and concentration fields have typi- cally stratified structures. With any decrease in N, the isotherms and concentration lines become more linear, which means the convective mechanism is suppressing. For the case when N = -4, heat and mass transfer take place mainly by conduction, which is evident from the temperature and concentration fields in Figure 6, where the isotherms and isoconcentration lines are almost parallel to the vertical and horizontal walls, respectively. The values of the Nusselt and Sherwood numbers for the presented examples are given in Table 5. The results are compared to reference [24], where the same problem is solved for planar and spatial geometries with the use of the finite-difference method. Good agreement between the results can be obtained, which proves the accuracy of the mathematical model and also of the numerical scheme. Figure 6. Streamlines, isotherms and isoconcentrations for A = 2, Ra = 100, Da = 10-3, Le = 10 and N = -2 (left side) and N = -4 (right side). Table 5. Nu and Sh numbers for different values of N and A = 2, Ra = 100, Le = 10, Da = 10-3. N -4 -2 -1.5 -1 0 Nu This study 0.49 0.50 0.50 0.52 0.81 Literature [24] - 0.50 0.50 0.55 - Sh This study 1.04 1.17 1.30 1.71 3.39 Literature [24] - 1.20 1.30 1.90 - 4 CONCLUSION A numerical approach, based on the boundary domain integral method (BDIM), which is an extension of the boundary element method (BEM), was applied for the solutions of the transport equations in porous media. The modified Navier-Stokes equations, Brinkman-extended Darcy formulation with the inertial term included, were employed to describe the fluid motion in porous media. The porous media is saturated with a viscous, incompressible fluid. A velocity-vorticity formulation of the governing equations is adopted, resulting in the computational decoupling of the kinematics and kinetics of the fluid motion from the pressure computation. Since the pressure does not appear explicitly in the field function conservation equations, the difficulty connected with the computation of boundary pressure values is avoided. The proposed numerical procedure is applied to the case of natural convection in a porous layer where the horizontal walls are subjected to different values of temperature and concentration and in a porous layer, where the horizontal walls are subjected to different values of concentration and the vertical walls to different values of temperature (convection under the cross gradients of temperature and concentration). The examples for different values of the modified Rayleigh number, the Darcy number, the Lewis number and buoyancy coefficient are presented and compared with the published studies. Very good agreement between the results obtained with alternative numerical methods (the finite difference method, the finite volume method, and the finite element method) can be observed. Since the final set of equations results in a very large and fully populated system matrix, the computer time and memory demands are large. The need to improve the economics of the computation, when using the boundary domain integral method, is still one of the challenges for researchers. One of the efficient mathematical tools, developed especially for saving computational time and computer storage with the boundary element method, is the wavelet transform [25]. However, it can be stated that the boundary domain integral method, extended in a way that also enables an investigation of the fluid-transport phenomena in a porous media, appears to possess the potential to become a very powerful alternative to existing numerical methods, e.g., finite differences or finite elements, as a means to obtain solutions to the most complex systems of nonlinear partial differential equations, when attacking some unsolved problems in engineering practice. REFERENCES [1] Nield, D.A. and Bejan, A. (2006). Convection in Porous Media. 3rd.ed., Springer, Berlin. [2] Vafai, K. (2005). Handbook of Porous Media. Taylor&Francis, Boca Raton, London, New York, Singapore. [3] Pérez-Gavilán, J.J. and Aliabadi, M.H. (2000). A Galerkin Boundary Element Formulation with Dual Reciprocity for Elastodynamics. International Journal of Numerical Methods in Engineering, 48, 1331-1344. [4] Blobner, J., Hribersek, M. and Kuhn, G. (2000). Dual Reciprocity BEM-BDIM Technique for Conjugate Heat Transfer Computations. Computer Methods in Applied Mechanics and Engineering, 190, 1105-1116. [5] Skerget, L., Alujevic, A., Brebbia, C.A. and Kuhn, G. (1989). Natural and Forced Convection Simulation using the Velocity-vorticity Approach. Topics in Boundary Element Research, 5, 49-86. [6] Skerget, L., Hribersek, M. and Kuhn, G. (1999). Computational Fluid Dynamics by Boundary-Domain Integral Method. International Journal of Numerical Methods in Engineering, 46, 1291-1311. [7] Jecl, R., Skerget, L. and Petresin, E. (2001). Boundary Domain Integral Method for Transport Phenomena in Porous Media. International Journal of Numerical Methods in Engineering, 35, 39-54. [8] Jecl, R. and Skerget, L. (2003). Boundary element method for natural convection in non-Newtonian fluid saturated square porous cavity. Engineering Analysis in Boundary Elements, 23, 963-975. [9] Wu, J.C. (1982). Problem of General Viscous Flow. Chapter in Developments in BEM, Vol. 2., P.K. Banerjee, R.P. Shaw (eds.), Elsevier Applied Science Publication, London, UK. [10] Hribersek, M. and Skerget, L. (1996). Iterative Methods in Solving Navier-Stokes equations by the Boundary Element Method. International Journal of Numerical Methods in Engineering. 39, 115-139. [11] Bear, J. and Bachmat, Y. (1991). Introduction to Modeling of Transport Phenomena in Porous Media. Kluwer Academic Publishers, Dordrecht, Boston, London. [12] Skerget, L., and Jecl, R. (2002). Boundary element method for transport phenomena in porous medium. Chapter in Transport phenomena in porous media II, Derek B. Ingham, Ioan Pop (eds.), Elsevier Science. [13] Skerget, L. and Samec, N. (2005). BEM for the two-dimensional plane compressible fluid dynamics. Engineering Analysis in Boundary Elements, 29, 41-57. [14] Brebbia, C.A. and Dominquez, J. (1992). Boundary elements. An Introductory Course. McGraw-Hill, New York. [15] Amahmid, A., Hasnaoui, M., Mamou, M. and Vasseur, P. (1999). Double-diffusive parallel flow induced in a horizontal Brinkman porous layer subjected to constant heat and mass fluxes: analytical and numerical studies. Heat and Mass Transfer, 35, 409-421. [16] Kladias, N. and Prasad, V. (1989). Natural convection in horizontal porous layers: Effects of Prandtl and Darcy numbers. Journal of Heat Transfer, 111, 926-925. [17] Nield, D.A. (1968). Onset of thermohaline convection in a porous medium. Water Resources Research, 4, 553-560. [18] Nield, D.A., Manole, D.M. and Lage, J.L. (1993). Convetion induced by inclined thermal and solutal gradients in a shallow horizontal layer of a porous medium. Journal of Fluid mechanics, 257, 559-574. [19] Rudraih, N., Srimani, P.K. and Friedrich, R. (1982). Finite amplitude convection in a two component fluid saturated porous layer. International Journal of Heat and Mass Transfer, 25, 715-722. [20] Trevisan, O.V. and Bejan, A. (1987). Mass and heat transfer by high Rayleigh number convection in a porous medium heated from below. International Journal of Heat and Mass Transfer, 30, 2341-2356. [21] Rosenberg, N.D. and Spera, F. J. (1992). Thermohaline convection in a porous medium heated from below. International Journal of Heat and Mass Transfer, 35, 1261-1273. [22] Mamou, M., Vasseur, P., Bilgen, E. and Gobin, D. (1995). Double-diffusive convection in an inclined slot filled with porous medium. European Journal of Mechanics B/Fluids, 14, 629-652. [23] Mohamad, A. A. and Bennacer, R. (2001). Natural convection in a confined saturated porous medium with horizontal temperature and vertical solutal gradients. International Journal of Thermal Sciences, 40, 2001. [24] Kalla, L., Vasseur, P., Beji, H. and Duval, R. (2001). Double diffusive convection within a horizontal porous layer salted from the bottom and heated horizontally. International Communications in Heat and Mass Transfer, 28, 1-10. [25] Mohamad, A. A. and Bennacer, R. (2002). Double diffusion natural convection in an enclosure filled with saturated porous medium subjected to cross gradients; stably stratified fluid. International Journal of Heat and Mass Transfer, 45, 3725-3740. [26] Ravnik, J., Skerget, L. and Hribersek M. (2004). The wavelet transform for BEM computational fluid dynamics. Engineering Analysis in Boundary Elements, 28, 1303-1314. RAZVOJ IN ANALIZA 3D GEOLOŠKO/GEOMEHANSKIH MODELOV ZA ČRPALNO ELEKTRARNO KOZJAK_ VLADIMIR VUKADIN o avtorju Vladimir Vukadin IRGO - Inštitut za rudarstvo, geotehnologijo in okolje Slovenčeva 93, 1000 Ljubljana, Slovenija E-pošta: vlado.vukadin@irgo.si izvleček Za določitev geoloških in geotehničnih pogojev gradnje črpalne elektrarne ČHE Kozjak smo izvedli obsežne geološke, geofizikalne in geotehnične preiskave. Rezultati preiskav so bili uporabljeni za razvoj in izdelavo 3D geološko/geomehanskih modelov. Modeli so bili izdelani z združitvijo vseh razpoložljivih podatkov v značilne geomehanske enote, ki so bile v naslednjem koraku postavljene v prostor skupaj s predvideno geotehnično konstrukcijo. Tak pristop je bil izbran zaradi kompleksnosti geotehničnih objektov in velikosti preiskovanega območja. S pomočjo prostorske analize je mogoče lažje razumeti odnose med posameznimi enotami in elementi, ter hitro ugotoviti najbolj kritična območja s stališča geotehnične konstrukcije in jih tudi stabilnostno preveriti. V članku sta prikazana modela za akumulacijski bazen in strojnico. Predstavljen je tudi pregled opravljenih preiskav skupaj z rezultati, ter razpravo o obsegu in vrsti preiskav z vidika izdelave 3D geološko/geomehanskih modelov. Ključne besede geološke in geomehanske preiskave, 3D geološko-geomehanski modeli, ČHE Kozjak THE DEVELOPMENT AND ANALYSIS OF 3D GEOLO- gical/geomechanical MODELS FOR THE KOZJAK PUMPED-STORAGE HYDROELECTRIC POWER PLANT VLADIMIR VUKADIN About the author Vladimir Vukadin IRGO - Institute for Mining, geotechnology and environment Aškerčeva 12, 1000 Ljubljana, Slovenia E-mail: vlado.vukadin@irgo.si Abstract In order to specify the geological and geotechnical conditions for the construction of the Kozjak pumped-storage hydroelectric power plant (HPP), the site area has undergone considerable geological, geophysical and geomechanical investigations. Due to the complexity of the geotechnical constructions and the size of the investigated site a novel approach was proposed, where the results of the investigation were used to develop 3D geological/ geomechanical models. The models were created with the integration of all the available results and data into typical geomechanical units that were subsequently given spatial dimensions. Through an analysis of the models the relationships between the different units were better understood and the critical areas were quickly identified and verified. In this paper the models for the accumulation reservoir and the engine room are presented. The overview of the executed investigations and the results are given together with a discussion on the scope and the number of executed investigations from the standpoint of the 3D geological/geomechanical model's creation. Keywords Geological and geomechanical investigations, 3D geological and geomechanical models, Kozjak HPP 1 INTRODUCTION First, some basic technical information about Kozjak HPP is presented, so that the size and the complexity of the proposed construction can be better understood. This is followed by a description of the type and the number of investigations for the accumulation reservoir and the engine room. The information about the general geological setting for the wider area of the investigated site is followed by a description of the 3D-model development for the accumulation reservoir and the engine room. An example of the use of the 3D models for detecting the critical areas of the reservoir and the engine room is shown. A discussion on the scope and the number of executed investigations from the standpoint of a 3D geological/geomechanical model development is also presented. Finally, the conclusions are presented. 2 KOZJAK HPP- BASIC TECHNICAL DATA The future Kozjak pumped-storage hydroelectric power plant (HPP) is to be built on the left bank of the River Drava, near the Fala hydroelectric power plant, some 20 km away from Maribor in the northeast part of Slovenia. The general location of the proposed power plant is shown in Figure 1 (see next page). The pumped-storage power plant is composed of three key facilities: the accumulation reservoir, the pipeline and the engine room. The turbine of the power plant will be located at the bottom of a 85-m-deep engine-room shaft at approximately 200 metres above mean sea level. The shaft will be located in the field, in close proximity to the accumulation reservoir of Fala HPP, so that the water from there can be used to fill up the reservoir of Kozjak HPP. The accumulation the reservoir of Kozjak HPP will be located 995 meters above mean sea level, on top of the hill named Kolarjev vrh. The reservoir will have an area of approximately 300,000 m2 (880 m x 350 m) and will hold approximately 3 million cubic metres of water. It will be protected by a 25-m-high embankment. The engine room and the reservoir will be connected by a 2.4-kilometre-long pipeline that will run from 85 to 300 metres below the surface and will Figure 1. Ge neral loca tion of the proposed power plant. have a diameter from 3.6 to 4 meters. The characteristic longitudinal cross-section of Kozjak HPP is shown in Figure 2. 3 OVERVIEW OF THE EXECUTED INVESTIGATIONS 3.1 INTRODUCTION The idea of constructing Kozjak HPP is more than 30 years old, and preliminary investigations were carried out by [1] and revised by [2]. Since then the quality of the investigation techniques and, especially, the requirements in accordance with the standards and design approaches have changed. Therefore, new investigations were necessary. From November 2006 to November 2008, geological, geophysical and geomechanical inves- tigations were carried out at the location of Kozjak HPP. The purpose of these investigations was to identify the geotechnical conditions for the construction of all three facilities: the engine-room shaft, the accumulation reservoir and the underground pipeline. The investigations for the pipeline are not presented in this paper. Based on their different operational purposes, design requirements and sizes, the number of investigations and the geotechnical considerations were somewhat different for each facility. Some of the investigations, however, were carried out for Kozjak HPP as a whole, and these will be presented first, followed by the investigations carried out for each reservoir and engine room. 3.2 INVESTIGATIONS CARRIED OUT FOR THE GENERAL AREA OF KOZJAK HPP In order to obtain information on the tectonic and lithological, engineering geological and hydro-geological characteristics at the different facilities, general investigations of the wider area were carried out. These general investigations included the following. Morphotectonic analysis The morphotectonic analysis of the digital model of relief (DMR) was carried out in order to determine the most common azimuths of the lineaments present. An area of 10x10 km with a point density of 12.5x12.5 m was studied. The results showed the occurrence of four distinct families of possible structural lineaments in the region. The recognized directions and traces served as the input information, used and verified in structural-geological and engineering geological mapping as possible fault lines. Figure 2. Longitudinal cross-section of Kozjak HPP. Structural geological, engineering geological and hydro-geological mapping An area of 12 km2 around the Kozjak HPP was mapped in detail on a topographical scale of 1:5000. The reviews of results from a previous mapping were carried out and the results from Mioc et al. [3] and Masic et al. [1] were used. Due to the considerable topographic diversity and the thickness of the vegetation cover, the method of profiling and outcrop tracking was used. The mapping consisted of the lithological and engineering geological identifications of the rocks and soils. In addition, the structural characteristics of the rocks, such as the strike and dip of the foliation and lineation, the orientation of the joints, the orientation of the mesoscopic folds, and the orientation and presence of the tectonically disturbed fault zones, were logged. For the purpose of the faults' kinematic analysis, the slips on the fault zones were measured as well. During the engineering geological mapping, the thickness of the soil cover was measured where possible and the GSI at the outcrops was determined. The hydrogeological mapping consisted of field observations, the logging of springs and other sources of water discharges. 3.3 INVESTIGATIONS CARRIED OUT FOR THE RESERVOIR BASIN The accumulation reservoir is to be located at the top of the hill and protected by a 25-m-high earth-filled embankment. The embankment will be constructed with local material obtained from the excavation of the basin bowl, so detailed investigations of the local material quality were necessary. An earthquake sensibility and stability assessment was one of most important tasks, due to the position of the accumulation basin and the height of the embankment. The investigations were carried out according to those requirements and they included the following. Structural drilling and excavation pits In the area of the accumulation reservoir 16 structural boreholes with a total length of 353 m were drilled. Intact quality samples were required, so diamond-bit crowns with double core barrels and bentonite rinse were used. The lengths of the individual boreholes varied from 15 to 28 m, and they were all placed at the toe or beneath the crown of the proposed embankment. In addition to boreholes, four pits up to a depth of 5 m were excavated: two beneath the crown of the proposed embankment and two inside the basin area. Geological logging, in-situ measurements and geophysical investigations The obtained cores and excavation pits were photographed and logged in detail in terms of the lithology, the foliation dip, the discontinuities dip, the spacing, the filling and the roughness of the joints, the RQD, etc. At each metre of core several measurements of the uniaxial compressive strength were carried out using a Schmidt hammer. The rock-classifications indexes GSI and RMR were determined for each metre of the borehole. In the excavation pits a penetrometer was used to measure the undrained shear strength of the soils. From the boreholes and the excavation pits the soil and rock samples were taken for further geomechanical, petrologic and mineralogical laboratory investigations. A total of 42 pressiometric measurements were carried out in the 16 boreholes using an OYO rock elastometer. Each investigation contained three measurements with three load/unload loops in the selected 3.5-m-long section. The values of the pressiometric elasticity loading and unloading modules were obtained from those measurements. The hydraulic conductivity of the soils and rocks was obtained from 16 slug tests. They were carried out in 12 boreholes and on average each test lasted for 10 hours. The purpose of the geophysical investigations in the area of the accumulation was to obtain information about the tectonic setting, the dynamic elastic properties and to define the design of the earthquake-acceleration parameters. The following investigations were carried out: six 240-m-long seismic refraction profiles (measurements of both the p and s waves), a micro-tremor with 32 triaxial point setting and a full spectral analysis in all three directions as well as down-hole measurements in 11 boreholes (measurements of the p and s waves). Laboratory investigations From the boreholes and the excavation pits a total of 180 samples of soils and rocks were taken. In the geomechanical laboratory the following investigations were carried out: the granularity, the natural and dry volume weight natural moisture, the atterberg limits, the investigation of the strength of the rocks (point load tests, uniaxial compression test, triaxial shearing tests, direct shearing tests along discontinuities) and the soil (direct shear test), the investigation of the rocks' elastic properties (E-v) and the edometric tests on the soils. In addition, special investigations for the usage of local material for the embankments were carried out. This included water-adsorption tests, methyl-blue tests, suction tests, a mineralogy investigation, frost-sensibility tests, crushing sensibility tests, etc. 3.4 INVESTIGATIONS CARRIED OUT FOR THE ENGINE-ROOM SHAFT The engine room will be constructed in an 85-m-deep and 30-m-wide circular shaft, and at the bottom of the shaft a junction with an underground pipeline is designed. For an excavation pit of these dimensions, the stability of the excavation has to be a major concern. The engine room is located 80 metres below the water level of the Fala HPP accumulation basin, on one side, and at the foot of the 800-m-high slopes on the other side. Under such a setting the underground water inflows under pressure had to be taken into consideration and investigated. The following investigations were carried out. Structural and rotary-percussion drilling Around the outer perimeter of the engine-room shaft, three structural boreholes, each 110-m-deep, were drilled. One rotary-percussion borehole at the centre of the excavation with no core recovery to a depth of 90 m was also drilled. This borehole was used for geophysical and hydro-geological investigations. Geological logging, in-situ measurements and geophysical investigations For the geological logging the same principles were applied as for the reservoir area. In the boreholes, 33 pressiometric investigations and a few short pumping tests were carried out. In addition, the measurements with a multi-parametric mini probe were made: the temperature (T), the pH, and the electric conductivity (n) were measured along the borehole. Along with the pumping tests, Lugeon tests were also carried out in the boreholes. Due to the instability of the borehole walls, some tests could not be completed; however, in total, 25 tests were carried out in three boreholes. Based on the hydro-geological investigations performed in the boreholes, at least two different origins of the water sources were confirmed. The geophysical investigations were aimed at detecting the faults that could otherwise not be detected due to an alluvium cover. The following methods were used: electrometric profiling of the specific electrical resistance parameters (eight 100-m-long profiles), the electrical resistivity method or misse-a-la-masse (in four boreholes), seismic refraction profiling (two 240-m-long profiles) and seismic tomographic profiling in three directions with the hydrophones located in all four boreholes. Logging with a video camera was also carried out in all four boreholes. Laboratory investigations From the boreholes a total of 90 samples of rock were taken for the laboratory investigations, where similar investigations were carried out as for the accumulation reservoir. Petrologic analyses were carried out on 19 specimens prepared from the core samples of the boreholes. The optical effects in the scanning minerals under scanning light enable a qualitative identification of the minerals and their structure. On the basis of their structure and texture, the sequence of the minerals' separation can be determined, and by considering the stability fields of the individual minerals and the observed reactions between them, the environment of their formation and the geological history of the investigated rock can be identified. Microscopic analyses were required for the calibration of the macroscopic observations, the identifications and the descriptions, which were carried out for all the cores. 4 RESULTS 4.1 GENERAL GEOLOGICAL CONDITIONS The Kozjak HPP will be located in the Kobansko region, which is part of the Central Alpine metamorphic rocks complex. The region has passed the medium-grade of metamorphism, with the rocks having undergone subsequent retrograde metamorphism. The most frequently occurring rock in this territory is mica schist, with transitions into gneiss; in the lower part of the complex the rocks are strongly carbonized and marbles are often present. The accumulation reservoir will be completely located inside the phyllonite block, which is the highest-lying unit. It is built up of strongly altered slates and gneisses. The rocks show extremely strong signs of shear deformations: a mylonithic structure, dense foliation and penetrative shear belts, which are usually heavily weathered, with good cleavage along the foliation. This unit outcrops in the highest parts of the region and has sub-horizontal contact with the underlying unit. The unit represents a sub-horizontal shear zone, which is a maximum of 200 m thick and runs along the main boundary plane between the Pohorje extensional complex and Kozjak [4]. A geological map, which is the result of geological mapping, is shown in Fig. 3. In addition to lithology, the fault systems obtained with mapping and geophysical investigations are also shown in Fig. 3. It is clear that there is a major fault system running along the north-eastern section of the proposed reservoir. Figure 3. Geological map of the investigated area. The pipeline is going to be built in a characteristic sequence of metamorphic rocks, which appears through the entire region of Kozjak and the neighboring Pohorje massif. It is predominantly built up of gneisses and subordinately of mica schists. This unit outcrops in the major part of the mapped territory. Inside this unit, lenses and bodies of amphibolites, pegmatite with a mylonithic structure and quartzite appear in the mapped territory. Inlays and lenses can be quickly wedged out laterally. The usual dip of the metamorphic rock strata is in the direction 180-220°, with a strike angle of 20-40°, although the foliation inside these rock strata changes quickly. Only a few major faults are projected to cross the proposed pipeline direction. The designed engine room is located on the Drava river terrace, where no outcrops were present. Since the surface is covered with quaternary sediments, we obtained the data from the boreholes and geophysical investigations. The terrace is made of up to 10-m-thick gravel deposits that are covered with debris on the hill side of the terrace. Beneath the gravel, metamorphic rocks are found with prevailing marbles, schists and gneisses that are also frequently carbonized. No fault lines are shown in Figure 3 in the area of the engine room due to the scale of the map and the fact that the terrain is covered with alluvium deposits. The tectonic element could not be traced on the surface, but geophysical investigations and boreholes showed the presence of faults that cross the location of the engine room, as will be shown later. 4.2 DEVELOPMENT AND ANALYSIS OF A 3D MODEL FOR THE RESERVOIR AREA Before the 3D geological/geomechanical model was developed, an analysis of all the available data and results had to be carried out. As was shown in the previous section the reservoir will be built entirely inside the phyllonite complex. The first step in the model development requires that the geological units are interpreted from the standpoint of their engineering and geomechanical properties. Based on the results from in-situ measurements, laboratory results, the mappings of tectonic disturbance and the weathering of cores, the phyllonite was divided into four geomechanical units. In Fig. 4 an example of the core taken from the borehole is shown, where all four typical units are present. Figure 4. Typical geomechanical units of phyllonite shown in the borehole core. The first unit (1) or proluvium, which is not shown in Figure 4, represents a surface layer where the phyllonite pieces are still found in the clay matrix, and which is deposited under a layer of humus. The thickness of proluvium changes from 1 m to 4 m and has no engineering significance, since it will be removed. Under the proluvium there is usually a layer of heavily tectoni-cally disturbed and weathered phyllonite (unit 2) with a strong clay content and with only traces of original foliation still present. The thickness of the stratum is highly variable and changes from a few metres up to 20 m in the area of the fault zones. The next unit of phyllonite (unit 3) represents moderately tectonically disturbed phyllonite with distinctive foliation and the original structure still present. The thickness of this unit is variable and is often interchanged with unit 4, which represents a moderately to weakly tectonically disturbed phyllonite with a more distinctive structure. The last unit (unit 5) represents weakly disturbed to undisturbed, unweathered phyllonite, which was only found in the lower part of some boreholes and was identified in other areas by the correlations with geophysical investigations. The ranges of the proposed strength and deformability properties for each unit are summarized in Table 1. The strength parameters for unit 2 were obtained as a median value with a 95% level of confidence, based on 13 laboratory direct shear tests, and the deformability properties were similarly obtained from a statistical analysis of 17 pressiometer measurements. For units 3, 4 and 5 the Hoek & Brown strength criterion [5], [6] was applied using RockLab software for a determination of the strength and deformability parameters of the rock mass. The geological strength indexes (GSIs) [7] were taken conservatively for each unit at the lower margins. The uniaxial compressive strength was obtained from laboratory tests. The intact elastic modulus of the rock was obtained from the unloading pressiometeric measurements. They were two to three times higher than the modulus measured in the laboratory. The decision to take pressiometric measurements as Table 1. Geomechanical parameters for typical units in the area of the reservoir. Unit GSI UCS Proposed strength parameters Strength parameters measured in the laboratory Elastic modulus c P c P E MPa kPa O kPa O MPa 2 - 0.1-0.4 (0.23) 5.2 32.5 0-17.5 30.5-40.8 32 3 10-20 2.1-9.6 (5) 14-36 30-36 27-33 39-44.3 60-90 4 20-40 8.8-21.2 (14.6) 37-84 44-51 31-47 47-55.6 306-1071 5 30-50 18/.f"67.4 92-275 53-56 25-143 41-65.6 3541-13368 an input value was made because the laboratory samples were influenced by foliation and discontinuities, so the values of the intact modulus were affected. For comparison, the strength parameters obtained from the triaxial tests are presented in Table 1. It is clear that there is generally good agreement between the values determined using the Hoek & Brown strength criterion and the laboratory values for units 4 and 5. For unit 3 the values obtained using the Hoek & Brown strength criterion are lower than the ones obtained in the laboratory. The reason for that could be the conservative approach with the GSI assignment and also due to the tendency towards the better quality samples selection in the laboratory. The same discrepancy for unit 3 was also noted when the deformability parameters were compared. Based on all the information obtained with the investigations, a 3D geological/geomechanical model was developed using GoCad (Geological Object Computer Aided Design) software. GoCad is a powerful tool used for manipulation, geostatistical analyses, kinematic analyses and the visualization of georeferenced geological, geomechanical and geophysical data. The model is constructed using points, linear elements (wells), curves, planes and volume elements, which in addition to spatial information can also carry other relevant data. The development of the model for the reservoir was carried out in the following steps. 1. The digital model of the terrain was formed from the point elements. 2. The embankment was constructed according to the design requirements with curve elements and planes. 3. The boreholes were inserted into the model as line elements. Each borehole element contained the information (markers) about the position of the boundaries between the different geomechanical units. 4. The geophysical refraction profiles, together with the geophysically interpreted fault planes, were inserted into the model as raster and plane elements. 5. The information obtained from the geological mapping, including the trust line and the fault planes, was inserted into the model and reinterpreted, based on the results from the geophysical investigations. 6. The plane boundaries of the different geomechanical units were interpolated between the boreholes and the fault planes. An example of a model under development is shown in Fig. 5, with surface, boreholes (containing geomechani-cal unit markers) and fault planes already inserted. Figure 5. 3D geological/geomechanical model of reservoir under development. Figure 6. Final 3D geological/geomechanical model of the accumulation reservoir area. The final 3D geological/geomechanical model of the accumulation area is presented in Figure 6. It includes a thrust plane between the gneiss and phyllonite, as well as the expansions of the planes between the geomechanical units 1, 2, 3 and 5. The geomechanical units 3 and 4 have not been separated due to their irregular interchange in both the vertical and horizontal directions. Thinner segments of the geomechanical unit 2, inside units 3 and 4, have also not been separated. In the area of fault zones, the occurrence of the geomechanical unit 2 was considered up to a depth of 20 m, and the thickness of the fault zone was taken at 10, for smaller, and 20 m, for more prominent, faults. The width of these fault zones was determined on the basis of the geological mapping, the geophysical investigations and the results of drilling and excavations. The spatial representation of the geomechanical data enabled us to predict the spatial expansion of the geome-chanical units in the areas where the investigations were not carried out. It also gave us a better understanding of the relations between the different geomechanical units and their relationship with the proposed design solutions. With the use of the visualization tools in GoCad, we were able to investigate the different areas of the proposed embankment. In Figure 7, a cross-section of the model in the direction SW-NE is shown. It is clear that in this area we have a thick layer of unit 2 so that the embankment is founded completely inside this unit. With further analyses, a similar area was also found in the NW part of the reservoir, as was expected. This area of the reservoir is crossed by a major fault zone that can also be seen in the model and is shown in Figure 6. With the use of the spatial geological/geomechanical model 2D cross-sections were drawn in most critical areas and analyzed. The results of the analysis carried out showed that the embankment can be built from local material [8], but due to the earthquake stability the slope of the embankment cannot exceed 26°. 4.3 DEVELOPMENT AND ANALYSIS OF A 3D MODEL FOR THE ENGINE ROOM The geological conditions around the engine room are more diverse in terms of lithology in comparison to the accumulation reservoir. Schists, gneisses are quickly changing with their more carbonized varieties; marbles and pegmatites are sometimes also present. The difference between the different varieties is in the mineral composition that is a consequence of the metamorphoses level and the origin rock. Based on lithology and petrologic investigations, ten different types of meta-morphic rocks were identified. In Figure 8 photographs of the cores taken from the two different boreholes at the engine-room area are shown. As can be seen in the left-hand picture, the changes from gneiss to marble and back are seamless. Accumulation reservoir H Critical area Figure 7. Critical cross-section of the reservoir model. However, of all ten identified lithological types only a few are relevant due to their predominant occurrence and were according to the lithology, the laboratory results, the geophysical and pressiometric results, divided into four geotechnical units, with all of them shown in Figure 8. The gneisses and schists represent unit A, their carbonized varieties represent unit B, the marbles unit C, and in unit D all three types together with cataklazits are present. They are interchanging at a very fast rate so they cannot be divided into sensible geomechanical units. Each of these units is further divided into subunits in accordance with the GSI values. The range of the geome-chanical properties for each unit summarized in Table 2 was obtained with the Hoek & Brown strength criterion using RockLab. Like with the reservoir area, the GSI values for each unit were taken conservatively at the lower margins. The input parameter for the intact elastic modulus was obtained from the average value of the laboratory measurements and the unloading pressiomet-ric elastic modulus. The uniaxial compressive strength for each unit was determined as an average value of the laboratory measurements and the results obtained with the Schmidt hammer, corrected to the core diameter. For units A and D several different average values were given. For example, in unit A a value of 3 MPa was used for the fault zone and in unit D different values are assigned to rocks with a different degree of disturbance. The lowest values of the GSIs (values between 5 and 10) and the uniaxial compressive strength were used for a characterization of the fault zones running through the respective unit. The range of the proposed values of the strength and deformability properties for each unit is summarized in Table 2 (on next page). Figure 8. Most common lithological units in the area of the engine room. Table 2. Geomechanical parameters for typical units in the area of the engine room. Unit GSI UCS Proposed strength parameters Strength parameters measured in the laboratory Elastic modulus ac MPa c kPa P O c kPa P O E MPa A 5-10 10-20 20-30 30-50 50-70 70-80 3 51.7 24 62 340 280 814 1952 32 28 41 55 51 54 408 55 64 256-384 3517 9381 30556 47971 B 30-50 50-70 70-80 68 442 1074 2566 53 50 54 264 58 7042 22937 36009 C 50-70 70-80 57 1154 2198 49 53 4274 37.5 62316 97831 D 5-10 10-20 30-50 10 35 79 226 469 20 33 42 - - 64 306 14528 For comparison, the strength parameters obtained from the laboratory are also presented in Table 2. Even though not all the subunits were examined in the triaxial shear apparatus, where the data was available a generally good agreement was found between the values determined by the Hoek & Brown strength criterion and the laboratory values for the representative subunit. The 3D geological/geomechanical model of the engine-room area was created following the same procedure as for the development of the reservoir model. It was created using the results from boreholes' loggings, inputted as line elements (with markers), seismic tomography profiles (raster images and fault planes) and the results of electric scanning. The model under development is shown in Figure 9 and the final model that includes the fault planes and the boundaries between the different units is shown in Figure 10. Figure 9. 3D geological/geomechanical model of engine room under development. A final model was analyzed in order to access the most critical area and possible failure mechanisms that could endanger the excavation. The most critical cross-section, drawn from the model, is shown in Figure 11a, where a sub-vertical fault can be seen, crossing the engine-room shaft. A spatial analysis revealed that 85% of the entire rock mass around the engine room is composed of rock with a GSI higher than, or equal to, 30, with the unit B being the most dominant unit. The other 15% of the rock mass is found in the faults and fault zones or in the first several metres of weathered rock below the alluvium. Based on the analysis a further simplification of the spatial model was possible and the 3D model for the numerical analysis, shown in Figure 11b, was developed. W CROSSECTION 2 E 0 5 10 15 2D 25 35 4Û 45 50 55 60 65 70 75 30 35 90 95 LCD LB5 L1G Escxvatiou bottom ciliitilc .....: : __/W 7 'Oil V m \ ' "it H a) b) Figure 11. a) Most critical cross-section in the engine-room area, b) 3D numerical model. A 3D numerical analysis was carried out with the simplified model using the 3D Plaxis Foundation. The geomechanical characteristics used in calculations are shown in bold in Table 2. The excavation to a depth of 85 m in 2.5-m excavation steps, with a primary lining installation and dewatering, was modelled. The results of the analysis showed that if the appropriate support measures (anchors, installation of lining, dewatering and length of the excavation step) are exercised, the engine room can be safely executed. With the use of the 3D geological/geomechanical model an efficient spatial analysis was possible, resulting in simplifications to the numerical analysis. 4.4 DISCUSSION ON THE SCOPE OF THE EXECUTED INVESTIGATIONS The investigations at Kozjak HPP were carried out in order to obtain the necessary geological and geome-chanical characteristics for the safe and economical construction of the reservoir, the pipeline tunnel and the engine room. For the accumulation reservoir an area of 300,000 m2 was covered with 16 boreholes and 4 excavation pits, amounting to one investigation point per 15,000 m2. Because all the boreholes were located in pairs at the base and below the crown of the embankment, the distance between the borehole pairs is a more appropriate measure. A 2,500-m-long embankment was covered with a borehole pair at every 300 meters, leaving the interior of the reservoir with an area of 120,000 m2, covered by only 6 refraction profiles and 4 excavation pits. Due to the relatively uniform geological conditions and the fact that the embankment, which is the dominant geotechnical structure, was sufficiently covered by the investigations, we can conclude that the executed investigations were adequate. In order to obtain more detailed information for the excavation of the interior of the basin, a few more excavation pits, boreholes or refraction profiles, could be executed. But with the help of a 3D model a good estimate was also obtained in this area. The engine room is relatively small in comparison to the reservoir. An area of 3,600 m2 was investigated with four deep boreholes and geophysical investigations. In addition to the results from four 40-m boreholes, the refraction seismic profiles from a previous investigation phase were available [1]. From all those investigations, relevant information for the excavation of the engine-room shaft and for the development of the 3D model were obtained. Due to the different origins of the water, the possibility of water inflow through the fault systems has not been fully explained. The already planned investigations that include the drilling of several additional boreholes at different depths and the execution of pumping tests in the boreholes with parallel monitoring and the analysis of water, should give an appropriate answer also to those questions. 5 CONCLUSIONS In this paper the investigations carried out in the area of the future Kozjak pump-storage hydroelectric power plant were briefly described and discussed. In the area of the reservoir a few more excavation pits or boreholes inside the proposed basin could be carried out, but in general the investigations were adequate. The same can also be stated for the engine-room area, with the exception of the underground water investigations. In general it can be concluded that based on the executed investigations, sufficient data is available in order to access the important geotechnical aspects of the construction of the Kozjak pumped-storage hydroelectric power plant. The results of all the investigations were synthesized in the form of 3D geological/geomechanical models and 3D models for the reservoir and the engine room were presented. The use of the 3D geological/geomechanical models was very helpful in determining the most critical areas and understanding the mechanisms that can endanger the structure. Those critical sections can then be investigated and a 2D stability analysis or full 3D numerical analysis, as was demonstrated for the engine room, can be executed. Models are also important in areas where we have to make extrapolations of the geological and geomechanical characteristics due to a lack of investigations, as was the case in the interior of the basin. The use of the 3D geological/geomechanical model was shown to be very useful and is probably necessary when complex geotechnical problems are being addressed. From the geological and geotechnical aspects the main hazards connected with the reservoir are linked with the static and dynamic stability of the embankments in the areas where the phyllonite unit 2 is thick (the subsurface areas and the area of the fault zones). Those problems could be addressed with the reinforcement of the embankment toes with pile walls, on the one hand, or with lowering the height and slopes of the embankments, on the other. The main hazards connected with the engine room are linked with the possibility of a strong water inflow during the excavation, which has still not been fully investigated. There are also some possibilities of local block instabilities, which could be prevented by the appropriate excavation sequence and temporary support measures (shotcrete, anchors) ACKNOWLEDGMENTS This work was financed by DEM d.o.o. and the author would like to express his thanks to Mr. Igor Cus and Mrs. Alenka Prnaver, who are the project managers of the Kozjak HPP project. The author would also like to express thanks to colleges from IRGO, and especially to Dr. Marko Vrabec and Nina Jurecic, with whom he worked on the 3D geological/geomechanical models. REFERENCES [1] Masič, M., Ribičič, M., Mervič, I. (1980). Geološko-geotehnično poročilo za potrebe projekta ČE Kozjak. Geološki zavod Ljubljana. [2] Škrabl, S. et al. (2004). Študija geotehničnih pogojev za potrebe projektiranja ČHE Kozjak. Univerza v Mariboru, Fakulteta za gradbeništvo, št. 21-10/04. [3] Mioc, P. and Žnidarčič, M. (1972). Tolmač in list Slovenj Gradec. L33-44, Zvezni geološki Zavod Jugoslavije, Beograd. [4] Fodor L., Balogh K., Dunkl I., Horvath P., Koroknai B., Marton E., Pecskay Z., Trajanova M., Vrabec Ma., Vrabec Mi., and Zupančič, N. (2004). Deformation and exhumation of magmatic and metamorphic rocks of the Pohorje-Kozjak Mts. (Slovenia): constraints from structural geology, geochronology. Proceedings of the 9th Meeting of the Czech Tetconic Studies Group, 2nd Meeting of the Central European Tetconic Group, Lučenec, Slovakia, June 22-25, 2004. Geolines, 17, 31-32. [5] Hoek, E. and Brown, E.T. (1980). Empirical streng-hth criterion for rock Masses. J. Geotech. Eng. Div., ASCE, 106 (GT9), 1013-1035. [6] Hoek, E. and Brown, E.T. (1997). Practical estimates of rock mass strength. Int. J. Rock Mech& Geomechanics Abstracts, 34 (8), 1165-1186. [7] Marinos, P. and Hoek, E. (2000). GSI: A geologically friendly tool for rock mass strenght estimation. International Conference on Geotechnical & Geological Engineering (GeoENG2000), Melbourne, Australia, 1422-1442. [8] Vukadin, V. et al. (2007). Geološko-geomehansko poročilo za potrebe izdelave projektne dokumentacije ČHE Kozjak. Mapa 2: Akumulacijski bazen, Ljubljana. OPTIMIRANJE IN SPREMLJAVA VISOKIH VARNOSTNIH STEBROV PRI PODZEMNEM PRIDOBIVANJU BLOKOV NARAVNEGA KAMNA JOŽE KORTNIK o avtorju Jože Kortnik Univerza v Ljubljani, Naravoslovnotehniška fakulteta, Oddelek za geotehnologijo in rudarstvo Aškerčeva 12, 1000 Ljubljana, Slovenija E-pošta: joze.kortnik@ntf.uni-lj.si Izvleček Podzemno pridobivanje blokov naravnega kamna se je kot prvo v Sloveniji, pričelo poskusno uvajati leta 1993 v kamnolomu pisanega apnenca Hotavlje I in leta 2001 v kamnolomu Lipica IIpredvsem zaradi geološke zgradbe nahajališča, stanja kamnoloma, velikih količin odkrivke v primeru širjenja površinskega dela kamnoloma in zaradi vse večjih potrebah po surovini naravnem kamnu. Podzemno pridobivanje blokov naravnega kamna se izvaja z nahajališču prilagojeno oz. modificirano komorno-steberno odkopno metodo z pravilno/nepravilno razporejenimi visokimi varnostnimi stebri. Ker se podzemno pridobivanje blokov naravnega kamna izvaja relativno plitvo pod površino 10-40 m, je vrednost primarnega vertikalnega napetostnega stanja relativno nizka (<1.0 MPa) in s tem tudi znatno povečana možnost izpada klinov oz. blokov iz stopa odprtih podzemnih prostorov. F preteklih letih smo posebno pozornost namenili vgradnji napetostnih in deformacijskih sistemov za kontrolo načrtovanih dimenzij (širine in višine) velikih odprtih podzemnih prostorov (komor) in dimenzij visokih varnostnih stebrov ter sprotno spremljavo ter identifikacijo pojavov nestabilnosti v stropu in bokih velikih odprtih prostorov (komor). F članku so predstavljeni postopki za načrtovanje, optimiranje in spremljavo visokih varnostnih stebrov pri podzemnem pridobivanju blokov naravnega kamna. Ključne besede naravni kamen, visoki varnostni steber, komorno-steberna odkopna metoda, podzemno pridobivanje, kamnolom OPTIMIZATION OF THE FOR THE UNDERGROUND RAL STONE BLOCKS JOŽE KORTNIK HIGH SAFETY PILLARS EXCAVATION OF NATU- About the author Jože Kortnik University of Ljubljana, Faculty of Natural Sciences and Engineering, Department of Geotechnology and Mining Aškerčeva 12, 1000 Ljubljana, Slovenia E-mail: joze.kortnik@ntf.uni-lj.si Abstract For the first time in Slovenia, the underground excavation of natural stone blocks was introduced on a trial basis at the Hotavlje I colourful limestone quarry in 1993, and in 2002 also at the Lipica II quarry. This was primarily because of the geological structure of the site, the quarry's condition, the potentially large amounts of the overburden in the event of an expansion of the surface part of the quarry, and the increasing needs for this raw material, i.e., natural stone. Underground The underground excavation of natural stone blocks is done using a modified room-and-pillar excavation method that is adjusted to each site's characteristics, with regularly or irregularly distributed high safety pillars. Since the underground excavation of natural stone blocks is performed at a relatively shallow level under the surface, i.e., at a depth of only 10-40 m, the value of the primary vertical stress state is also relatively low (<1.0 MPa). This significantly increases the risk of wedge-shaped pieces or blocks falling out of the ceiling in open, underground spaces. In previous years, special attention was paid to the installation of stress-strain systems for controlling the planned dimensions (width and height) of large, open, underground spaces (rooms) and the dimensions of the high safety pillars, along with continual monitoring and identification of the instability phenomena in the ceiling and sides of the large open spaces (rooms). The paper presents the procedures for the planning, optimization and monitoring of high safety pillars for the underground excavation of natural stone blocks. Keywords natural stone, high safety pillars, room-and-pillar mining method, underground mining, quarry 1 INTRODUCTION The underground excavation of natural stone blocks is not an idea generated by the modern information society; it dates back to the times of the ancient Romans. Evidence found in South East Britain, in the county of Devon [11], i.e., in the now abandoned Beer quarry, which has been converted into a museum, shows us that the underground excavation of natural stone blocks was done here, even B.C., by the ancient Romans. In Europe, the underground excavation of natural stone blocks is nowadays mostly done in several quarries in Italy (Carrara, the Apuan Alps, Bolzano, etc.), Great Britain (Avon, Somerset, Dorset, etc. [2, 13]), Greece (Dionysos - Athens), Portugal (Solubema - Lisbon), Croatia (Kanfanar - Pazin), etc. In Slovenia, the underground excavation of natural stone blocks was first introduced on a trial basis at the Hotavlje I colourful limestone quarry in 1993, and in 2002 it was also implemented at the Lipica II limestone quarry. The Marmor Hotavlje (MH) company, as one of the leading Slovenian stone-cutting companies, began the organized excavation of natural stone at the Hotavlje I quarry in 1948, but the actual beginnings of the excavation of natural stone blocks at the Hotavlje I quarry date back to the 1800s. Natural The natural stone found here, the so-called "Hotaveljčan" limestone, is colourful (grey, red, pink, and sometimes almost black) and has white calcite veins as well as the remnants of individual corals and algae [14]. The MH management decided to introduce underground excavation, primarily due to the geological structure of the site, the condition of the quarry, the large amounts of the overburden in the event of an expansion of the surface part of the quarry, and the increasing needs for natural stone as a raw material. Marmor Sežana, which has been the leading stone-cutting company in the Karst region for over half a century, began its excavation of natural stone at the Lipica II quarry in 1986. The Lipica II quarry excavates two types of natural stone, which were named by the Karst stone-cutters as "Lipica Unito" (homogenous stone) and "Lipica Fiorito" (rose stone) [15]. In terms of size, the Lipica II quarry ranks among the largest Slovenian natural stone quarries. For similar reasons to the case of the Hotavlje I quarry, the Lipica II quarry also decided on a trial underground excavation of natural stone blocks in 2001 and introduced it in 2002. In both quarries, the underground excavation of natural stone blocks is done using the modified room-and-pillar excavation method, which is adjusted to the characteristics of the sites with irregularly spaced safety pillars. Since in both cases the underground excavation is done at a relatively shallow depth of about 34 m, the value of the primary vertical stress state is relatively low (<1.0 MPa). This significantly increases the risk that wedge-shaped pieces or blocks may fall out of the ceiling in open underground spaces. When planning the underground excavation, special attention therefore had to be paid to the engineering-geological mapping, which was initially done for the external surfaces of the future area of the underground spaces (i.e., galleries, transverse roads and niches, and, after deepening, also the rooms) and the structure of the productive layer. On the basis of these data, the predominant dike systems, which are important for the stability and consequentially also the safety of the underground spaces, were determined. However, this issue will not be addressed in greater detail in the paper. 2 PLANNING AND OPTIMIZATION OF THE ROOM-AND-PILLAR MINING METHOD In many underground quarries for the excavation of natural and technical stone, a modified room-and-pillar method is used, which is adjusted to the conditions of the site. It employs regularly (Figure 9) or irregularly spaced safety pillars that have as small as possible width-to-height ratios r [6]. This type of mining method enables the use of self-supporting rock as the support element, which is achieved by forming safety pillars of the appropriate dimensions during the excavation in order to ensure the stability of the hanging wall on the ceiling and the required span of open spaces between these pillars for safe access and operation of an underground excavation site. In planning the dimensions (the surface area and the height) of the safety pillars as well as the span of the open spaces (rooms) between them, mining engineers usually have to strike a balance between the requirements of the mining rights' owner, i.e., the client who has commissioned the excavation project, and the demands for the maximum possible yield of natural stone blocks, on the one hand, and the need to ensure a certain level of safety that is still adequate, on the other [9]. This is because, due to the limited amount of natural stone available at any site, the largest possible dimensions of open underground spaces and the smallest dimensions of the safety pillars are needed in order to obtain the maximum volume of natural stone blocks, although this also tends to reduce the safety of the underground spaces in an inverse proportional manner. By reducing the base dimensions of the safety pillars, their strength is reduced as well, while the loads and the risk of overloading the safety pillars increase in an inversely proportional manner [10]. Figure 9. Typical room-and-pillar layout with regularly spaced pillars for non-metal (stone) mines. The stability of the ceiling of large underground spaces and the safety pillars in them depend primarily on the quality of the geomechanical properties of the rock and the intensity of the tectonic activity at the site. Differences in stability occur due to the different mass mass-to to-volume ratios of the rock, the different loads created by the overburden and the different tectonic activities. This necessitates detailed planning of the excavation and the implementation of appropriate support measures in order to ensure safe working conditions. Therefore, the geomechanical properties of the rock, the primary stress state pz of the hanging wall and the tectonics of the site need to be known in detail in order to ensure the safe and stable excavation of natural stone blocks. In the case of the underground excavation of natural stone blocks at the Lipica II and Hotavlje I quarries, the use of a modified room-and-pillar method adapted to the two sites is planned, with irregularly distributed high safety pillars and a width-to-height of r < 1 (r= Wp/h). 2.1 DETERMINING THE DIMENSIONS OF HIGH SAFETY PILLARS WITH LOW WIDTH WIDTH-TO TO-HEIGHT RATIOS The strength of the safety pillars primarily depends on the strength of the rock and the value of a pillar's width-to-height ratio. Any reduction in the width-to-height ratio may cause a reduction in the total strength of a safety pillar. In existing underground mines of natural and technical stone, including the Lipica II and Hotavlje I quarries, the following common characteristics have been observed concerning safety pillars with a low width-to-height ratio (recommendations) [16]. The initial height of the safety pillars usually amounts to 4.5 m, but with deepening of the underground spaces it increases (3.0 m); in places it may reach values of up to 18.0 m or even more. Already, thin safety pillars are usually further weakened by each deepening of the basic level due to a reduction of their basic horizontal cross-section; they may collapse if the strength is reduced below the respective permissible limit due to a reduction in the value of the pillar's width-to-height ratio r. High safety pillars are those in which the value of the width-to-height ratio is r < 1. Initially, the stable safety pillars may become unstable with a gradual deepening of the basic panel. In many underground mines of natural and technical stone, deepening of the panel was stopped exactly because of the worsening of the condition of the rock in the deepened panels. The weakening or collapse of a single safety pillar may cause a chain reaction because it overloads the neighboringneighbouring safety pillars and leads to Figure 10. The dominant safety-pillar rupture mode involves; (a) spalling from the pillar surfaces, (b) shear fractures transecting the pillar, (c) lateral bulging or barrelingbarrelling of the pillar surface, (d) a pillar with set of natural transgressive fractures, (e) a pillar with a developed foliation or schistosity parallel to the principal axis of loading (formation of kink bands) [5]. settlement of the entire area of hanging rock above the pillars. This risk is especially high with thin safety pillars. Safety pillars whose width-to-height ratio is gradually reduced during excavation are especially sensitive to being weakened by transverse or vertical discontinuities (cracks or fractures). The geomechanical conditions in underground mines for natural and technical stone are normally very good, but they may quickly deteriorate because of an unexpected appearance of discontinuities, especially in the case of thin safety pillars (see figure 10). 2.2 DETERMINING THE LARGEST SPAN OF OPEN SPACES (ROOMS) BETWEEN SAFETY PILLARS Due to the increasing needs for raw materials and the desire to lower the excavation costs, mining companies tend to require an increase in the span of open, underground spaces. Because of the different geological conditions in individual quarries of natural and technical stone, the ceilings of large open rooms are most often composed of one or several layers of rock that may be parallel or inclined with respect to the surface. The stability of the ceilings in open spaces primarily depends on the geomechanical properties of the rock, the loading of the ceiling caused by the overburden (arch), and the tectonic conditions at the site. The usual spans of the open spaces in underground quarries of natural and technical stone amount to about 13.0 m, but they may be as high as 18.0 m. They also depend on the minimum dimensions of the excavation equipment (primarily the diamond diamond-belt sawing machine and the recommended dimensions of the natural stone blocks required for further processing [12]. With an increase in the span, the stability of the ceiling in an underground open space decreases because: - the bending stresses cause bending and twisting of the ceiling and lead to the appearance of shear cracks, - an increase in the bending and shear stresses within the ceiling may also cause critical damage to the intact rock or lead to the ceiling's collapse, - there is a higher probability that the ceiling of an underground space will be perforated by open dikes. As is the case with the planning of the dimensions of safety pillars, here as well one relies on practical experience acquired with the already-used ceiling spans in existing underground spaces, i.e., the trial-and-error approach. Methods to determine the appropriate spans of the ceiling of the underground spaces in quarries of natural and technical stone for various locally specific rock conditions and stress states are still in their initial stages of development. The methods most often used nowadays are, therefore, those which have been developed for coal mines and metal mines, while the actual conditions applying to the underground excavation of natural and technical stone are taken into account only indirectly. In the planning of the underground excavation of natural stone blocks using the room-and-pillar excavation method, it is also important to determine the Table 1. Calculation of the natural stone yield for the correct distribution of the predicted safety pillars for the Hotavlje I and the Lipica II quarry. W;1 W;2 2Dmodel Wr12 Pillar A area Tributary area Extraction/ Yield ratio Width-to-height ratio r h=4.5m/7.5m/10.5m [m] [m] [m] [m] [m2] [m2] [/] [/] 16.8 16.8 8.4 16.8 282.24 846.72 1 : 3.0 3.73/2.24/1.60 16.8 14.0 7.6 16.8 235.20 799.68 1 : 3.4 3.73/2.24/1.60 14.0 14.0 6.4 16.8 196.00 752.64 1 : 3.8 3.11/1.87/1.33 14.0 14.0 7.0 14.0 196.00 588.00 1 : 3.0 3.11/1.87/1.33 14.0 11.2 6.2 14.0 156.80 548.80 1 : 3.5 3.11/1.87/1.33 11.2 11.2 5.0 14.0 125.44 509.60 1 : 4.1 2.49/1.49/1.07 11.2 11.2 5.6 11.2 125.44 376.32 1 : 3.0 2.49/1.49/1.07 11.2 8.4 4.8 11.2 94.08 344.96 1 : 3.7 2.49/1.49/1.07 8.4 8.4 3.6 11.2 70.56 313.60 1 : 4.4 1.87/1.12/0.80 8.4 8.4 4.2 8.4 70.56 211.68 1 : 3.0 1.87/1.12/0.80 8.4 5.6 3.4 8.4 47.04 188.16 1 : 4.0 1.87/1.12/0.80 5.6 5.6 2.2 8.4 31.36 164.64 1 : 5.3 1.24/0.75/0,53 5.6 5.6 2.8 5.6 31.36 94.08 1 : 3.0 1.24/0.75/0,53 Figure 11. Regular distribution of safety pillars [5]. appropriate dimensions (width and height) of the large open spaces (galleries, transverse roads/crosscuts, niches and rooms) and the high safety pillars (see Figure 11) in order to achieve the optimal natural stone yield ratio (utility value). Table 1 presents the yields for different dimensions of pillars and open spaces. In practice, the yield ratios of 1 : 4.4 and 1 : 4.1 were found to be potentially the best for high safety pillars with respect to the mechanical properties of the rock. 3 THE PLANNING AND DESIGN OF HIGH SAFETY PILLARS The strength of safety pillars has been studied for decades by many researchers. The majority of studies in the past were focused on the research of safety pillars in coal mines, but some were also applied to rocks. As a result of these studies, it was found that the strength of a safety pillar is proportional to the strength of the rock in which it is located, and inversely proportional to its thickness: the thinner the pillar, the smaller is its load-carrying capacity. Among the methods used for planning safety pillars, the following two groups predominate: - Analytical methods, which are based on the mathematical principles of the mechanical behaviour of rocks and are computationally less difficult to execute. However, in spite of the possibility of fostering a better understanding of the mechanics of safety pillars, these methods have not been widely used in practice. Their primary disadvantage lies in the use of certain prescribed values (constants), which are difficult or almost impossible to determine in practical work. - Numerical methods, which use modern, numerical techniques and are computationally more demanding, are intended for modelling the loads on safety pillars and presenting changes in the rock's stress and strain states. Furthermore, they enable the modelling of special conditions by taking into account the faults and dikes, as well as the inclusion and assessment of the effect of weakened areas on the overall stability. Nowadays, numerical models play a very important role in the planning of safety pillars for special conditions. 3.1 ANALYTICAL METHOD FOR DETERMINING THE LARGEST SPAN OF OPEN SPACES (ROOMS) BETWEEN SAFETY PILLARS Analytical investigations are usually based on a determination of the static equilibrium of the rock. In such analyses, the average stress state is first determined within the support elements (i.e., the safety pillars) and then compared to the average value of the rock's strength (see Figure 12). Generally, stone pillars are less stable if the overburden is substantial because of the higher stress. Pillars are also less stable as the width-to-height ratio decreases, for example, in benching operations. The stress levels within pillars can be approximated by using the tributary area theory [5]. The average axial pillar stress level [19]; & p = P ■ g ■H ■ W + Wp, )2 + Wt P2 W ■W "pi "p 2 (1) w„ +w, Ww, +W. = pz pi p2 (Wpi Wp2 where Wp1 [m] Wp2 [m] Wrl [m] W2 [m] P [kg/m3] g [m/s2] H [m] pz [MPa] the pillar's width the pillar's length the room's width (gallery, crosscut, niche, room, etc.) the room's length (gallery, crosscut, niche, room, etc.) the density of overburden strata the acceleration due to gravity the thickness of the overburden the vertical normal component of the pre- mining stress field Figure 12. Cross-section over a horizontally positioned productive layer of uniform thickness, which is excavated by forming long rooms with a width of Wr and intermediary pillars with a width of Wp [5]. The pillar's stress levels are affected by the overburden and the relationship between the area supported by the pillar and the area of the pillar. The relationship is illustrated by comparing the post post-mining vertical stress levels as the overburden and the extraction ratio increase. The most generally accepted techniques for determining pillar strength, defined as the ultimate load per unit area of a pillar, use empirical equations based on survey data from actual mining conditions. The failings of the empirical method stem from an inability to extend these equations beyond the specific material properties, sizes, shapes and overburdens found in the survey data. Bien-iawski wrote that the strength of safety pillars depends upon three elements [4]: - the size or volume effect (strength reduction from a small laboratory specimen of rock to the full size safety pillars), - the effect of the pillar's geometry (shape effect), - the properties of the pillar's material. For non-coal pillars [8], empirical formulas have largely been derived from some form of the following power equation for the safety pillar's strength Sp , w; S; = O "¡7 (2) where a, b [MPa] the pillars strength [MPa] the pillar rock's uniaxial compressive strength [m] the pillar's height [/] the exponents determining the pillar's strength from its volume and shape strength rises rapidly. At higher width-to-height ratios (r > 1), strength increases occur at diminishing rates [9]. In other words, at some point the pillar would begin to display some plastic behaviour [3]. Pillar stability is most endangered at low width-to-height ratios. As typical stone safety pillars reach a width-to-height ratio of r > 1.5 [9], they begin to exhibit an almost indestructible character. Factor of safety Fs , Fs =-; >1,6 O p (3) The low factor of safety provided by this prospective layout indicates that a redesign is necessary to achieve the required factor of 1.6 [5]. The options are to reduce the room span, thereby reducing the pillar's stress level, to increase the pillar's width, or to reduce the pillar's (and mining) height. The selection of an appropriate safety factor can be based on a subjective assessment of the pillar's performance or a statistical analysis of the failed and stable cases. As Fs decreases, the probability of failure of the pillars can be expected to increase. In practical terms, if one or more pillars are observed to be failed in a layout, it is an indication that the pillar stress is approaching the average pillar strength, causing the weaker pillars to fail. The relationship between Fs and the failure probability, however, depends on the uncertainty and variability of the system under consideration. The value of the factor of safety Fs was calculated using Equation 3 for different values of the width-to-height ratio r and for different values of the uniaxial compressive strength of the pillar rock at a depth of 40 m below the surface. This is presented in Figure 14 (see next page). S p a c h This equation considers both the material's strength and the safety pillar's shape to calculate the pillar's strength. Table 2. Exponents determining the pillar's strength from its volume and shape (see Equation 2). Source a b Subject medium Hedley and Grant 0.5 0.75 Quartzite pillars; Ura- (1972) [7] nium mines near Elliot Lake, Canada; for w/h < 4.5 Stacey and Page (1986) [18] 0.5 0.70 In the planning of underground excavations of natural stone blocks using the room-and-pillar method, cautious use of the results of the empirical equation is required (2). At low width-to-height ratios (r < 1), the pillar's 3.2 NUMERICAL ANALYSIS OF THE STABILITY OF SAFETY PILLARS AND CEILING IN OPEN UNDERGROUND SPACES Nowadays, various program packages are available for the numerical analyses of the stability of safety pillars and ceilings of open underground spaces (FLAC 2D, FLAC 3D, PLAXIS, etc.). They are based on the finite finite-element method (FEM), difference method (FDM), the distinct distinct-element method (DEM), etc. The Fast Lagrangian Analysis of Continua (FLAC 2D) is a two-dimensional explicit finite finite-diferen-cedifference method (FDM). FLAC is well accepted by social mining and rock mechanics engineering. The main advantage of this method is the integration of the surrounding roof and floor conditions on the stone safety pillar strength [1]. Hedley/Grant h=4.5m Stacey/Page h=4.5m Hedfey/Grant h=7.5m Stacey/Page h=7.5m Hedley/Grant h=10.5m Stacey/Page h=10.5m Hedley/Grant h=13.5m Stacey/Page h=13.5m Hedley/Grant h=16.5m Stacey/Page h=16.5m 0,25 0,50 0,75 1,00 1,25 1,50 1,75 2,00 2,25 2,50 2,75 3,00 3,25 3,50 3,75 width-to-height ratio r Figure 13. Comparison between the pillar's width-to-height ratio and the average pillar strength for several different empirical equations (table 2. on previous page) based on a power function (Rock uniaxial compressive strength ffc=100 MPa). 13,0 12,0 11,0 10,0 £ 9,0 £ 8,0 « 7,0 0 L_ 6,0 0 0 b,0 Li. 4,0 3,0 2,0 1,0 0,0 0,00 —— ?c = 45 MPa — ?c = 40 MPa 1 1,99 7 c = 35 MPa c = 30 MPa c = 25 MPa 7 — ?c = 20 MPa 6 35 jr 4,5 3 y/ 3,99. --" 5 ,33 ; —** 2 B2 ,2,34 1 1 0,25 0,50 0,75 1,00 1,25 width-to-height ratio r 1,50 1,75 2,00 Figure 14. Factor of safety for shallow underground excavation spaces (thickness of overburden H=40 m, pillar width W=12m). For the numerical analysis the FLAC 3.3 software package was used. The purpose of the numerical analysis was to determine the stability of the planned dimensions of the underground rooms, to make a comparison between the deepening of the levels in monolithic rock without failure and in rock that has failed because it had cracks and dikes, and to provide the geotechnical foundations for the planning dimensions in an underground excavation, along with continued surface excavation at the Lipica II quarry. Figure 15. Stability analysis of the safety pillars using the FLAC 2D software package. According to data from the literature, the ratio of the horizontal to the vertical component of the primary stress state varies widely in the area close to the surface. For underground room depths of up to about 100 m, the value of the coefficient k (k = ah/av) is between 1.3 and 3.5. In the majority of cases close to the surface, the value of the horizontal stress ah is greater than that of the vertical stress av [8]. Underground rooms at the Lipica II and Hotavlje I quarries are located close to the surface, and the lowest height of the overburden is 39 m or 36 m. The value of the coefficient k = 1.3 was therefore used in the numerical analysis. On the basis of previously presented data, several models were made to present the deepening of a pair of galleries-rooms with widths of 13 m, 15 m, 16 m and 20 m in compact and failed limestone. Based on the modelling results, it was concluded that: - A gallery-room with a width of up to 13 m is stable up to a width-to-height ratio of r = 0.28 if the geomechanical properties of compact limestone are used; if those of failed limestone are taken into account, then up to a width-to-height ratio of r = 0.48, - A gallery-room with a width of 15 m is stable up to a width-to-height ratio of r = 0.36 if the geomechanical properties of compact limestone are used; if those of failed limestone are taken into account, then up to a width-to-height ratio of r = 0.76, - A gallery-room with a width of 20 m is stable up to a width-to-height ratio of r = 1.78 if the geomechanical properties of compact limestone are used; - The models showed that galleries-rooms with a width of 13 m, 15 m or 20 m can be deepened up to a width-to-height ratio of r = 0.28 by using support measures in the form of local anchoring of the lower half of gallery-room sides. Galleries-rooms with a flat ceiling remain stable even if a factor of safety Fs = 1.6 is used. The factor of safety is taken into account in the model so that the geome-chanical properties of limestone are reduced by the corresponding percentage of the safety factor. Table 4. The following geomechanical properties were used in the models: Parameters Compact limestone Model 3. Fractured limestone Model 1. in 2. Very fractured limestone E Module of Elasticity 14.7 GPa 9.8 GPa 5.6 GPa y Density 26.3 kN/m3 26.3 kN/m3 26.3 kN/m3 v Poisson's ratio 0,3 0,3 0,3 T Tensile strength 1 MPa 0,5 MPa 0,5 MPa f Angle of internal friction 52° 35° 29° c Cohesion 1.8 MPa 1.2 MPa 0.7 MPa 4 IN-SITU MEASUREMENTS In-situ control measurements for the room-and-pillar mining method include measurements of the stress state (2D stressmeter) as well as the strains (EL-beam sensors, multipoint extensometer, meter for determining the displacement of open dikes) within safety pillars and in ceilings of large, open, underground spaces. Only the results of in-situ control measurements done at the Lipica II quarry are shown below due to the more intense excavation of natural stone blocks and the longer monitoring of these measurements. To perform control measurements of the changes in the stress state of high safety pillars, a 2D stressmeter (VW (vibrating wire) biaxial stressmeter model 4350-1) manufactured by Geokon was used to monitor the main stresses in a single vertical plane perpendicular to the axis of the drill hole. Measurements of the main stresses are enabled by three VW sensors, which are oriented at 60° angles within the probe. The stressmeter also has a sensor for temperature measurements. The stressmeter body is made of a steel cylinder with a maximum external diameter of 57.1 mm [17]. Technical characteristics of 2D stressmeter (Model 4350 BX): Standard stressmeter range 70 MPa Resolution1 14 do 70 kPa Accuracy ±0.1% F.S. Temperature range (253 - 353 K) -20°C do +80°C Borehole diameter BX (60 mm) 1 Depends on rock modulus For transferring data from the stressmeter, a memory unit (datalogger CR10 module, AVW1, SC32B) is used for the data capture, along with the appropriate software (the PC200W software package). The data capture is done automatically, using the time interval set in the program (1 min, 60 min or 240 min). Figure 16., 17. VW Biaxial Stressmeter and its position in the borehole. The VW1 stressmeter used to monitor the changes in the primary stresses in the vertical plane perpendicular to the axis of the drill hole was installed in the SP02 safety pillar. The site of the stressmeter installation corresponds to the site of monitoring the primary stresses in the numerical model for the case of deepening of the gallery pairs. The results of the measurements of the stress state in the SP02 safety pillar are shown in Table 5. Table 5. Average measured values of the main stresses in the SP 02 safety pillar. VW1 stressmeter Temperature Sig_1 Sig_2 k width-to-height ratio r [°C] [MPa] [MPa] [/] 1.80 5.4/14.5 -2.60 -2.09 1.24 1.10 3.0/15.1 -3.03 -2.29 1.32 0.80 2.7/15.6 -3.62 -2.67 1.35 # Major principal stress Sig_1 in safety pillar SP02 —Flac Model 1. -•-Flac Model 2. —*— Flac Model 3. —•—VW1 stressm. (aw.) Figure 18. Comparison of measured and calculated (with the aid of numerical modelling) values of the primary vertical stresses Sig_1 in the SP02 safety pillar. Minor principal stress Sig_2 in safety pillar SP02 —♦— Flac Model 1. -•- Flac Model 2. Flac Model 3. —•—VW1 stressm. (avr.) Figure 19. Comparison of measured and calculated (with the aid of numerical modelling) values of the primary horizontal stresses Sig_2 in the SP02 safety pillar. With the reduction of the width-to-height ratio r, in-situ measurements of the stress state exhibit a trend of increase in the primary stresses. Even at the same width-to-height ratio r, the primary stresses in the SP02 safety pillar oscillate slightly due to temperature changes in the rock (summer/winter). The values measured in situ are within the range of the results obtained with analytical calculations and numerical modelling, while the results for the primary horizontal stresses deviate slightly. Up to a width-to-height ratio r = 1, the results of the in in-situ measurements of the stress state shown in Figures 18 and 19 remain within the range of the results obtained with numerical modelling. With the continued deepening, i.e., a reduction of the width-to-height ratio r < 1, in-situ measurements indicate a trend of greater increase of the primary stresses than were calculated during the numerical modelling. This necessitates additional analyses and prompt monitoring of the in-situ measurements during the continued deepening, i.e., a reduction of the width-to-height ratio r. With a reduction of the width-to-height ratio r, the ratio of the vertical to horizontal components of the primary stresses in the safety pillar is within the range of 1.24 to 1.35 (1.30). If open dikes appear in a safety pillar, relatively simple dike displacement meters are additionally installed in order to monitor the sliding surfaces within the dikes (Figure 21). 0 10 30 40 50 Figure 20. Plain view of the Lipica II underground quarry (VW1 - vibrating wire stressmeter, safety pillars SP01 to SP04). Figure 21. Manual measurement of the displacements of open dikes within safety pillars. A dike dike-displacement meter consists of three screws installed along a dike and arranged in the form of a equilateral triangle (see Figure 20). The measurements of dike displacements are done manually with an adjustable gauge, by measuring the changes in the distance between the screws. At the Lipica II quarry, three dike displacement meters are installed: one in the P2 P2-to to-SP02 line, one on the right-hand side of the G2 gallery and the third one on the right-hand side of the G1 gallery (see Figure 20). Because of difficult access to the meters after the first deepening, the measurements of dike displacements were done periodically. The results of the measurements for the relative displacements of open dikes in safety pillars SP02 to SP04 did not show any active displacements of the sliding surfaces in the pillars. The results of analytical calculations and numerical modelling showed that in the case when the geome-chanical properties of the compact limestone were taken, a gallery-room with a width of up to 13 m is stable up to a width-to-height ratio of r = 0.28, without having to employ any additional support measures. If the geome-chanical properties of failed limestone were used, then such a gallery-room is stable up to a width-to-height ratio of r = 0.48. On the basis of the modelling results, the width of the portal portion of the pillar was also estimated (along the cross-section that is perpendicular to the surface levels), which had to be greater than 13 m. The results of in-situ measurements of the stress state at the Lipica II quar ry in the SP02 safety pillar confirmed the results of the numerical modelling. The measurements of the dike displacements also do not indicate any displacements of the sliding surfaces in the area of open dikes in the safety pillars SP02 to SP04. For the time being, no methodology is available for dimensioning high safety pillars with a low width-to-height ratio for underground quarries of natural and technical stone. The experience and results of measurements currently obtained in both Slovenian quarries that employ the underground excavation of natural stone will be beneficial in the development of a new methodology for the implementation of this underground excavation method in other natural stone quarries that are suitable for its use. The pillar-design guidelines developed through the observational, analiticalanalytical and numerical simulations discussed above will require further field confirmation. This approach can help to form a part of the comprehensive pro-active mine safety ground-control plan for underground natural stone mines. 5 CONCLUSIONS AND FUTURE perspectives In the planning of an underground excavation of natural stone blocks using the room-and-pillar excavation method, special attention needs to be paid to the determination of the appropriate dimensions (width and height) of large, open, underground spaces (rooms) and high safety pillars, as well as the installation of appropriate systems for continual monitoring and identification of the instability phenomena in their ceilings. Due to the large heights (even in excess of 20 m) of such open, underground spaces, deepening of the plane renders access to the ceiling for any repair work or the installation of additional supports more difficult or even impossible. LITERATURE [1.] Anon, (1998). FLAC: Fast Lagrangian Analysis of Continua User's Guide. Itasca Consulting Group, Inc., Minneapolis. [2.] Bajželj, U., Kortnik, J., Petkovšek, B, Fifer, K., and Beguš, T. (1999). Okolju prijazno podzemno pridobivanje naravnega kamna. RMZ 46/2, 203-214. [3.] Barron, K. (1984). An Analytical Approach to the Design of Coal Pillars. CIM Bull., 77, 868, 37-44. [4.] Bieniawski, Z.T. (1984). Rock Mechanics Design in Mining and Tunneling. Balkema, Rotterdam, 193-209. [5.] Brady, B.H.G., Brown, E.T. (1985). Rock Mechanics for Underground Mining. George Allien&Unwin (Publisher) Ltd., 316-350. [6.] Esterhuizen, G.S., Dolinar, D.R., and Ellenberger, J.L. (2008). Pillar Strength and Design Methodology for Stone Mines. Proceedings of the 27th International Conference on Ground Control in Mining, Morgantown, West Virginia University, 241-253. [7.] Hedley, D.G.F. and Grant, F. (1972). Stope-and-Pillar Design for the Elliot Lake Uranium Mines. CAM Transaction, 75, 121-128. [8.] Hoek, E. and Brown, E.T. (1997). Practical Estimation of Rock Mass. Int. J. Rock Mech., 34, 8, 11651186. [9.] Iannacchione, A.T. (1999). Pillar Design Issues for Underground Stone Mines, Proceedings of the 18th International Conference on Ground Control in Mining, Morgantown, West Virginia University, 271-281. [10.] Iannacchione, A.T. and Prosser, L.J. (1997). Roof and Rib Hazard Assessment for Underground Stone Mines. SME Preprint 97-113, SME Annual Meeting, Denver, CO, 5. [11.] Old Beer quarry (2008), http://www.beerquar-rycaves.fsnet.co.uk/BEERQ~S/Beer Quarry Cx.html. [12.] Kortnik, J. (2007). Optimiranje in spremljava varnostnih stebrov pri podzemnem pridobivanju blokov naravnega kamna. Zbornik strokovnega posvetovanja rudarjev in geotehnologov ob 40. skoku čez kožo, Ljubljana, 46-54. [13.] Kortnik, J. and Bajželj, U. (2005). Underground mining of natural stone in Slovenia, 20th World Mining Congress 2005, Iran, Tehran, 277-286. [14.] Marmor Hotavlje (2008), web page: http://www. marmor-hotavlje.si/. [15.] Marmor Sežana (2008), web page: http://www. marmorsezana.com/. [16.] Mine Safety & Health Administration (2005), web page: http://www.msha.gov/stats/. [17.] Instruction manual Model 4350BX Biaxial stressmeter (2008), http://www.geokon.com/. [18.] Stacey, T.R., Page, C.H. (1986). Practical Handbook for Underground Rock Mechanics, Series on Rock and Soil Mechanics, Trans Tech Publications, 12, 53-63. [19.] Wagner, H. and Frommer, T. (2004). Changing the Mining Method in an Austrian Magnesite Mine. Proceedings of the 20th World Mining Congress, Tehran, Iran, 211-215. RAZVOJ VTISNIH FILTROV ODVODNJEVRLNEGR SISTEMA PREMOGOVNIKA VELENJE S POMOČJO METODE KONČNIH ELEMENTOV GORAN VIŽINTIN, MIRAN VESELIC, ANDREJ BOMBAŽ, EVGEN DERVARIČ JAKOB LIKAR Ln ŽELJKO VUKELIŽ o avtorjih Goran Vižintin Univerza v Ljubljani, Naravoslovnotehniška fakulteta Aškerčeva 12, 1000 Ljubljana, Slovenija E-pošta: goran.vizintin@guest.arnes.si Evgen Dervarič Univerza v Ljubljani, Naravoslovnotehniška fakulteta Aškerčeva 12, 1000 Ljubljana, Slovenija E-pošta: evgen.dervaric@ntf.uni-lj.si Miran Veselič Univerza v Ljubljani, Fakulteta za gradbeništvo in geodezijo Jamova 2, 1000 Ljubljana, Slovenija E-pošta: miran.veselic@guest.arnes.si Jakob Likar Univerza v Ljubljani, Naravoslovnotehniška fakulteta Aškerčeva 12, 1000 Ljubljana, Slovenija E-pošta: jakob.likar@ntf.uni-lj.si Andrej Bombač Univerza v Ljubljani, Fakulteta za strojništvo Aškerčeva 6, 1000 Ljubljana, Slovenija E-pošta: andrej.bombac@fs.uni-lj.si Željko Vukelič Univerza v Ljubljani, Naravoslovnotehniška fakulteta Aškerčeva 12, 1000 Ljubljana, Slovenija E-pošta: zeljko.vukelic@ntf.uni-lj.si Izvleček Podzemna dela v premogovniku Velenje ves čas ogroža prisotnost podzemne vode. Zaradi zapletene geološke zgradbe premogovnika je bilo izvrtanih veliko strukturnih vrtin in v njih izdelanih geofizikalnih karotažnih meritev. Nad debelim lignitnim slojem se nahaja večslojni pliocenski in pleistocenski vodonosnik, ki ga izmenično sestavljajo peščene in glinaste plasti. Leta 1981 je bil pliocenski vodonosnik razdeljen v tri vodonosne sisteme glede na hidravlične tlake v njih, črpalne teste, kemizem vode in geofizikalne lastnosti: a) prvi vodonosni paket (Pl1), b) vodonosnik 20 - 80 m nad premoško plastjo (Pl2) in c) zgornji pliocenski vodonosnik (Pl3). Za rudniška dela je najpomembnejši vodonosnik saturiranih peskov (Pl1). Hidravlični tlak v tem vodonosnem paketu direktno vpliva na varnost v spodaj ležečem rudniku. Ti vodonosniki so tudi najbolj pod vplivom odvodnjevanja. Ne glede na to pa so črpalne vrtine izdelane tako, da zajemajo cel paket Pl2 in ponekod tudi dele Pl3. Hidravlični tlaki v teh vodonosnikih lahko dosežejo vrednost do 35 barov, zato so pred leti pričeli z obsežnim programom zniževanja tlakov nad deli, kjer potekajo rudarske aktivnosti. Za opazovanje tlakov so bili na teh območjih izdelani več nivojski piezometri. Za predikcijo razvoja tlakov so bili uporabljeni številni 3D matematični modeli končnih diferenc (FDM). Praksa je pokazala, da so ti modeli za določanje regionalnega trenda toka podzemne vode dobri, glavna slabost pa se pokaže na točkah opazovanja tlakov podzemne vode, ki se nahajajo na mestih odvovodnjevanja z vtisnimi filtri, kjer prihaja do manjših napovedi znižanja od dejanskih. Ta slabost FDM je že dolgo poznana in se pojavlja zaradi večje odvisnosti napovedanega znižanja od velikosti celice kot pa od tokovnega polja. Nevarnost vdorov vode v premogovnik Velenej se bo še posebej povečala po letih 2012 in 2017, ko bodo zaradi posedanja terena, ki nastaja zaradi odkopavanja in subsidence, opustili površinska črpališča. Zaradi tega bo morala biti na novo postavljena strategija odvodnjevanja. Uničene odvodnjevalne objekte prvega reda bo potrebno nadomestiti z vtisnimi filtri, ki bodo izvrtani iz rudniških prog okoli odko-pnih čel. Za dimenzioniranje vtisnih filtrov uporaba FDM ni primerna glede na napake, ki nastanejo pri napovedih znižanj in črpanj. To je leta 2007privedlo do uporabe metode končnih elementov (FEM) za predikcijo znižanja podzemne vode in količin črpanja v območjih, kjer bodo rudniška dela izpostavljena nevarnostim vdorov vode iz površine. Na osnovi predikcije s FEM metodo smo lahko napovedali velikost in opremljenost vtisnih filtrov. Ključne besede vtisni filtri, matematično modeliranje podzemne vode, rudniške vode, rudniška hidrogeologija, geofizikalna karotaža & modeliranje geoloških struktur THe DeveLOPMeNT of r "drive-in" filters DEWATERING SYSTEM IN THE VELENJE COAL MINE USING FINITE-ELEMENT MODELLING GORAN VIŽINTIN, MIRAN VESELIC, ANDREJ BOMBAC, EVGEN DERVARIŽ JAKOB LIKAR and ŽELJKO VUKELIŽ About the authors Goran Vižintin University of Ljubljana, Faculty of Natural Science and Technology Aškerčeva 12, 1000 Ljubljana, Slovenia E-mail: goran.vizintin@guest.arnes.si Miran Veselič University of Ljubljana, Faculty of Civil and Geodetic Engineering Jamova 2, 1000 Ljubljana, Slovenia E-mail: miran.veselic@guest.arnes.si Andrej Bombač University of Ljubljana, Faculty of of Mechanical Engineering Aškerčeva 6, 1000 Ljubljana, Slovenia E-mail: andrej.bombac@fs.uni-lj.si Evgen Dervarič University of Ljubljana, Faculty of Natural Science and Technology Aškerčeva 12, 1000 Ljubljana, Slovenia E-mail: evgen.dervaric@ntf.uni-lj.si Jakob Likar University of Ljubljana, Faculty of Natural Science and Technology Aškerčeva 12, 1000 Ljubljana, Slovenia E-mail: jakob.likar@ntf.uni-lj.si Željko Vukelič University of Ljubljana, Faculty of Natural Science and Technology Aškerčeva 12, 1000 Ljubljana, Slovenia E-mail: zeljko.vukelic@ntf.uni-lj.si Abstract During the mining operations at the Velenje coal mine, groundwater has been presenting a constant threat to underground works. The hydrogeological setup is so complex that a lot of structural drilling and well-logging operations were needed in the past to clarify it. Above the lignite seam is a Pliocene and Pleisticene multilayer aquifer system, composed mainly of permeable sand layers and impermeable clay layers. In 1981 the Pliocene aquifers were divided into three packages. Based on the water-table data of each aquifer, pumping tests, chemical analyses of the groundwater and the geophysical properties the Pliocene aquifers directly above the seam, together with impermeable layers, were divided into: a) the first water-bearing sands (Pl1), b) the aquifers 2080 m above the coal seam (Pl2) and c) the upper Pliocene aquifers (Pl3). For the mining operations the most important aquifer of saturated sands is Pl1. The hydraulic pressure of the groundwater in these sands directly affects the safety of the mining. These aquifers are mostly affected by the dewatering activities, too. However, the dewatering wells are constructed in such a way as to capture the whole Pl2 and, somewhere, even a part of the Pl3 complex, too. The water pressure in this multilayer aquifer can reach over 35 bars, so a massive program of drawdown activities has been needed and is still in place to decrease the water table in the area related to the mining operations. Special, multilevel observation wells are used to monitor the water level. A number of 3D finite-difference models (FDMs) were used to estimate the regional groundwater drawdown. It was observed that the FDMs performed well when predicting the regional situation, but the model-predicted drawdown was lower than the observed values at observation points in the area where the dewatering operations using "drive-in" filters have taken place in the past. This is a well-known problem of the FDM: the drawdown being rather a function of the cell size than of the flow net. The risk of water inrushes will increase, especially after 2012 and 2017, when a series of surface-drilled wells, connected into the mine'spumping-line batteries, will be abandoned due to excavation works and mining-subsidence effects. Consequently, the dewatering schemes had to be completely reviewed. The destroyed, first-order dewatering structures will have to be replaced by a series of "drive-in" filters, drilled from the mine roadways in the area of the planned longwall face operations. For the drive-in filter-system design the FDM does not seem to be appropriate. This is especially so if the error in the drawdown andpumpingflowprediction is taken into account. That led, in 2007, to the selection of the finite-element method (FEM) for the prediction of the groundwater drawdown and the water pumping rates in the areas were the underground works will encounter the risk of a water inrush. Based on the FEM prediction the sizing and the layout pattern of the "drive-in" filters were made. Keywords drive-in filters, groundwater mathematical modelling, mining water, mining hydrology, geophysical well login & solid modelling 1 INTRODUCTION The Velenje coal mine is one of the largest and most modern collieries in Europe. It is situated in the northeastern part of central Slovenia in the Šaleška Valley tectonic depression, filled up with more than 1000 m of Pliocene and Pleistocene sediments (Fig. 1). The sediment material represents the complete sedimentation cycle of the terrestrial, swamp and lake phases in both directions. This sequence is sometimes interrupted with fluvial sediments transported from the northern and north-western side [7]. Coal extraction in the Velenje coal mine has been uninterrupted since the end of the 19th century, during which time several stopping methods for the excavation of wide coal seams have been tested. In the first half of the 20th century the room-and-pillar method as well as block caving were used. Since 1947, however, the long- wall mining method with improvements has been practiced; this is known as the Velenje longwall method[9]. The Velenje coal mine is producing lignite from a coal seam that has huge dimensions. The coal seam is a more or less continuous body: 8.3 km long, 2.5 km wide and more than 160 meters thick. The depth to the seam varies from 200 to 500 meters. In the central part, the exploitable seam is up to 168 m thick. In the central part the seam is the most deeply buried (approximately 450 m) and in the marginal parts it is closer to the surface (approximately 100 m) (Fig. 2). In the lower part of the coal seam the ash content gradually increases, while the upper margin is very sharp. The coal seam is covered by marl with fossil snail shells, followed by mudstone, sometimes laminated or massive (Fig. 2). Within this mudstone, intercalations of water-bearing sands and gravel of changeable thickness appear. This impermeable layer between the coal seam and the first sands on top of it is called the isolation or the protective layer. This isolation prevents water and mud inrushes into the mine openings [9]. However, large spaces are opened during the mining works, leading to a failure of the protecting layer, which allows the lowest sand and gravel layers of the hanging wall multi-aquifer system to come into contact with the lignite seam. Figure 1. Geological map of the surrounding area of the Velenje lignite mine underground works. The red lines representing the coal mine's active mine roadways. Because of this, the groundwater pressure in these aquifers, as the most relevant part of the hanging-wall multi-aquifer system, has to be constantly monitored [10] and groundwater modelling has to be used to predict the underground hydrological state. In the past, FDM numerical models were used. The most successful predictions until now were made with Visual ModFlow. Even though Visual ModFlow has been the best prediction tool and it is widely used in Velenje, its limitations were impeding its use in predicting the hydrological state on a short scale [8]. This problem has been very well known for a long time and different approaches need to be found. The main problem is that the FDM produces groundwater drawdown, which is more a function of the cell size than of the flow net [2] [3]. 2 HYDROGEOLOGICAL SETTING In general, in the Velenje coal mine three hydrogeologi-cal systems can be distinguished [12]. The upper one is a Quaternary system, which lies on the Pliocene aquifers. These Pliocene aquifers lying over the coal seam are the most important for mining safety. It was in 1981 that the roof aquifers were divided into the Pliocene and quaternary aquifers [12]. The division was based on water-table data in a single aquifer, pumping tests, chemical analyses of the water and the geophysical properties [15]. The Pliocene aquifers were then further divided into three packages. The aquifers directly above the seam and isola-tive layer, the first water-bearing sands (Pl1), the aquifers placed 20-80 m above the coal seam (Pl2) and the upper Pliocene aquifers (Pl3) [7]. From the point of view of safety criteria, the first and second sand layers (P11 & Pl2) are the most important. The high water pressure in the sands can directly affect the underground works below. The bottom elevation of the first sands (Pl1) in the years before the excavations is presented in Fig. 3. There are some difficulties in determining the aquifer thickness. The first problem is the heterogeneity of the whole aquifer complex, consisting of sand, gravel, mudstone, clay, and the various mixture layers [7]. Overall, the vertical and lateral permeability is, therefore, low. The second problem is to delineate the part of the complex affected by the dewatering operations, and the third one is the insufficient amount of data in the area of the missing coal seam [7]. The coefficient values were acquired through numerous pumping tests performed on the dewatering and observation wells. For the area where no wells were drilled, no values exist. Based on pumping test data, different zones of hydraulic conductivity could be established (Fig. 6). As can be seen, the hydraulic conductivity coefficient values are in range from k = 1.74 x 10-7 m/s up to k = 2.88 x 10-6 m/s (Fig. 6). Regarding the specific yield, only a few values are given. For modelling purposes a representative value of the specific storage S0 = 5 x 10-5 m-1 for the whole area was chosen. Because there were no data available for the effective porosity, we had to adopt an assumed value of nef = 5%. According to the references this is the value for dense, silty sands [6]. However, this value can be lower in the zones where more silt is present [7]. Figure 2. Geological cross-section (west-east) of the Velenje lignite mine. The multilayer aquifer system is composed of sand strata lying directly on the lignite seam [9]. 3 TIME EFFECTS OF MINING WORKS ON THE AQUIFERS Underground works will increase the open spaces, which are to be filled in with the material lying over the coal seam. The basic concept of excavating coal by the longwall method is that the affected area of exploitation extends above the supported roof coal of the face, thus enhancing the natural forces that break and crush the coal and/or assist the natural process. The excavation face is divided into the lower-excavation part and the upper-excavation part. The lower part is 3-4 m high and is protected by a hydraulic shield support, thus enabling mechanized coal production with shearers and haulage using chain conveyors. The upper-excavation part is 7-17m high and is exposed to dynamic stresses, which, in combination with blasting, cause the coal to disintegrate and crumble onto the conveyor. The direct roof crumbles into the cavity and consolidates over time, so the excavation of the lower panel is enabled. The winning of the upper part can be continuous or time delayed [9]. The mining operations will also affect the pumping rates. The pumping wells shown in the red rectangles in Fig. 3 & 4 will be disused in 2012 as a result of the coal-excavation operations affecting these wells (Fig. 3). As a consequence, the bottom-wall elevation of the sands Pl1 will change from that given in Fig. 3 to the state presented on the elevation contour map in Fig. 4. As has already been pointed out, the groundwater pressure in the sands Pl1 & Pl2 is the most important issue. From the measurement data using multilevel piezometers it was possible to reproduce the groundwater table (Fig 5). It is obvious, however, that the groundwater depression cone is the result of long-term dewatering in the pumping wells. To keep the effects of the pumping well battery abandoning at a minimum, the so-called "drive-in" filters, i.e., the pumping wells drilled vertically from underground roadways into the roof aquifers, will be applied (Fig. 5.). The emplacement and operation of "drive-in" filters will not start simultaneously, but their start and stop times will depend on the excavating processes with the longwall method. In Table 1 the start and stop times of the "drive-in" filters are shown. Table 1. The working scheme for "drive-in" filters is mainly dependent on the mining operations along the longwall. Name of "drive-in" filter Prescribed water head [m.s.l.] Start to work [day] Stop to work [day] VF-1 -54.70 0 510 VF-18 -64.00 14 477 VF-2 -57.60 28 477 VF-17 -61.40 42 442 VF-3 -50.30 56 442 VF-16 -50.10 70 409 VF-4 -41.90 84 409 VF-15 -43.90 98 372 VF-5 -31.20 112 372 VF-14 -35.00 126 337 VF-6 -20.80 140 337 VF-13 -23.60 154 302 VF-7 -8.60 168 302 VF-12 -13.40 182 267 VF-8 2.50 196 267 VF-11 -4.00 210 230 VF-9 16.30 224 230 VF-10 4.10 230 230 Figure 4. 3D bottom-wall elevation map for the Pliocene water-bearing sands Pl1 in 2012. As is clear from Graph 1, the central pumping line is now pumping more than 300 l/min. This stable pumping line is to be disused due to mining operations and has to be replaced with another one, drilled from the surface, or with a series of in-mine drilled "drive-in" filters. There are many objections to a stable pumping line, but the most important is that the pumping line will face the same end as the central pumping line with the mining works still going on. So, only the "drive-in" filters can be constructed for dewatering the part where longwall face mining operations are in progress. For this reason the model has to be very accurate in terms of groundwater-pressure predictions. Due to the very well-known problems of the FDM explained before, the FEM was selected. 5 MODELLING AND "DRIVE-IN" FILTER DEWATERING SYSTEM DESIGN To predict the effects of sand subsidence, of the shutdown of the existing pumping line due to mining operations and of a replacement of the shutdown wells with "drive-in" filters around the longwall face and underground roadways (Fig 5 & Fig. 7), a 3D FEM FeFlow 5.1 modelling program for flow, mass and heat transport was used [1]. For the "drive-in" filter design there are many hydraulic and geo-chemical parameters which have to be taken into consideration [13]. In the case of the Velenje coal mine, the geo-chemical parameters have less impact on the operation of the "drive-in" filters than the hydraulic ones. If we take into consideration that the "drive-in" filters will be in operation for less than 510 days the chemical or geo-chemical factors can be neglected. From this, the decision was taken that only the flow has to be modelled. The main objective was to model the maximum flow rates at the sites of the "drive-in" filters and to analyze the influence of the geometry of the water-bearing Pl1 sands on the flow pattern. The modelling in such a case has to be as much as possible realistic, so it was important to understand the processes that will take place in 2012 along the longwall. When the problem was analyzed more closely it was Figure 5. Groundwater-level contour map for the water-bearing sands Pl1 & Pl2 in 2007 (with north to the top). The underground vertical pumping wells are marked with the labels VF (black points). found that in fact two situations are the most important for the prediction: the first one, which has to show how the hydrological situation will evolve when the "drive-in" filters start to operate; and the second one, which has to show the groundwater-level rebound after the "drive in" filters stop. The design documents and papers state that the mining operations along the longwall face will take 510 days (Table 1). This includes the preparation works and the excavation works. Here, the time needed to move from one longwall face to another also has to be taken into consideration. From the mining documentation it seems that this time is mostly a function of the mining-operations management rather than a technical problem. However, the time needed to move from one face to another is also very important, especially if it is taken into account that no pumping via the "drive-in" filters will be possible during that period. The situation will be even more critical due to the fact that the central pumping line also has to be switched off because of the excavation works. To make the scenario more realistic it also had to be taken in account that the first sand water-bearing layer (Pl1) geometry will change with time. From a previous, detailed mining surveying [9] it is known that the subsidence caused by the groundwater pumping and excavation works amounts to nearly 80% of the excavating space height. Starting from this experience and feeding the current data to Surfer 8.0, the upper and lower borders of the water-bearing sands (Pl1) were first constructed for the current time situation (Fig. 3). Then, according to the planned future mining operations, defined in the mine design documents, a prediction for the Pl1 sands' geometry was made for 2012. The results of the gridding process are shown in Fig. 4. After the geometry data needed for the model were defined, in the next step the water-balance data were analysed. This is particularly important because the hydrological system of the Velenje mine is now in a quasi-steady state[11]. For this a series of measured flow data for the central pumping line battery and south line pumping battery were taken into consideration (Graph 1). The grids already presented in Fig. 3, Fig. 5, Fig. 6 and the data from Graph 1were used to calibrate the model for the current hydrological (Fig. 5) and mining situation (Fig. 3). After the calibration on the current hydrological and mining situation (Fig. 5) the calculated and calibrated heads were used as the initial condition for predicting the groundwater levels in 2012 (Fig. 4). For the flow boundary the Cauchy type was selected on the north, north-west and north-east parts of the model (Fig. 7). The prescribed "in-flow" in the model domain was equal to the amount of groundwater pumped out with the current pumping wells. The secondary constraint was also used in the form of minimum prescribed heads, which are in accordance with the Puping rates from dewattering wells 800.0 £ J 600X) a 400.0 □ Central pumping line battery ■ South pumping line battery Graph 1. Pumping rates for the central and south pumping line battery. In 2012 the central pumping line battery will cease to work. 3500 4000 4500 5000 5500 6000 0500 7000 7500 8000 GK y (meters) Figure 6. Map of the hydraulic conductivity zones of the modelled Pl1 sands deduced from a series of pumping tests. The zones were defined in one of the previous models in the ModFlow 4.2 FDM (North is orientated vertically). heads on the head map presented in Fig. 5. The water abstraction from the old remaining wells (aligned in the South pumping line) was selected as a single source/sink boundary (Fig. 7). As can be seen in Graph 1, the values of the pumping rates are more or less stable on both pumping lines batteries. For the pumping rates of each well their average rates from 2000 to 2006 were selected. For the estimation of the "drive-in" filters the pumping rates Dirichlet boundary conditions (BCs) were used. The heads of the Dirichlet BCs were prescribed to the maximum groundwater table drawdown by free-flow dewatering (Table 1), because of the possibility that the water level is dropping below the initial level and that instead of pumping, the recharge will begin at the Dirichlet BCs. The additional flow constraints, as the minimum outflows of greater than, or equal to, 0 m3/day were used. In the model, the starting time and the end of the operation time for the underground vertical wells were also taken into consideration. For the reason specified before, the development of the "in-flow" rates on the "drive-in" filters are the most important for the filter design. So, the model results were mainly used to predict the rates on spots of the "drive-in" filters. Using the flow budget analyser of the FeFlow the flow rates on the Dirichlet (BCs) were reconstructed, which were representing the "drive-in" filters. The model calculations are presented in Graph 2. The values in this graph are in range from a few decilitres per second up to a few litres per second. The results seem to be in line with the results of the pumping test during the 1980s [12], where in some cases pumping rates up to 15 l/s were measured. From the report of that pumping test [12] it can also be seen that if the pumping rates exceeded 5 l/s the wells were clogged by sand. The predicted "in-flow" rates are the largest at the beginning (Graph 2). This is typical of the start of the dewatering process. With time the "in-flow" rates start to decrease until the moment when the "drive-in" filters stop functioning one after another. After that, a strong increase in the predicted "in-flow" rates can be seen for the wells that remain operational. The amount of water inflow to the "drive-in" filters is also dependent on their position and the geometry of the strata. In the case that the geometry of waterbearing sands is very well known an inflow prediction can be made. If, for the sake of strata heterogeneity, the flow at a particular "drive-in" filter is different from that predicted it is usually compensated by another filter in the vicinity. Figure 7. FEM groundwater model. The red points represent the active pumping wells (south pumping line battery), the blue are the "drive in" filters (Dirichlet BCs) and the pink points represent the Cauchy BCs (North is vertically orientated). Graph 2. Predicted "in-flow" rates of the planned "drive-in" filters using Dirichlet BCs in the mathematical model. Figure 8. Predicted groundwater-table depression in 2012 after 230 days of mining operations on a longwall face. All the "drive-in" filters are in operation. (North is vertically orientated). For the designing of the "drive-in" filters a criterion of maximum flow velocity has to be taken into account [4]. The maximum water velocity through the open spaces of the screen (i.e., the screen opening) should not exceed 0.03 m/s [16]. Based on the model prediction it was possible to estimate the maximum expected flow rate to the "drive-in" filters so the calculation could be made after the equation for designing the necessary screen opening on the screened part of the "drive-in" filters [16]: A = F - Qf Vdop - l (1) where Af effective open area per metre of screen F well-clogging factor Qf flow through the open filters [m3/s] l length of the screens in sands vdop maximum admissible screen entrance velocity The percentage of screen opening can be calculated [16] from: P = - Af n - D„ (2) where P percentage of screen opening D diameter of well ("drive-in" filters") As can be seen, the maximum flow rates, which can go up to 13.5 l/s, are in the initial stages of dewatering around the longwall, but according to the Graph 2 the flow rates at the "drive-in" filters start to decrease with time. So, the "drive-in" filters have to be dimensioned for the period of maximum flow rates. Using equations 1 and 2, this calculation was performed for the clogging factors 1, 2 and 3. Graph 3. Predicted drawdown curves for the observation points between the "drive-in" filters positioned along the roadways. Table 2. Calculated percentage of screen opening. Dc [mm] F=1 F=2 F=3 Percentage of screen opening 48.3 32.95 65.90 98.85 60.3 26.39 52.79 79.18 76.1 20.91 41.83 62.74 88.9 17.90 35.81 53.71 108 14.74 29.47 44.21 The parameters used for the calculations are as follows: Qf 13.5 l/s l 10 m vdop °.°3 m/s A length of 10 m was selected for the screens in the water-bearing sands (Pl1). This value is based on the Velenje mine experience. The reasons are the diameter of the roadways and the time that would be necessary to increase the diameter. This also means that the drilling equipment has to be of limited dimensions and the drilling operations are restricted by difficult working conditions. According to the values in Table 2 a criterion of 0.03 m/s is very difficult to achieve. So, on one hand, we have a result that is showing that the dewatering can be very successful with the "drive-in" filters and, on the other hand, we can see that, physically, almost no small diameter screens are able to meet such requirements. According to Table 2, the required small-diameter screen-opening percentage, as calculated from data of the models, exceeds the technically possible 25 % if any additional screen clogging was taken into account. However, there is also another condition that we have to take into consideration, and this is related to the maximum velocity of the groundwater flow through the open spaces of the intergranular porous media[5] around the well. This is the well-known Sichardt criterion, which stipulates that the maximum withdrawal from the well also has to take into consideration the flow velocity at which a failure of the aquifer structure due to a displacement of fine fractions could immediately occur. So, the velocity should not exceed the velocity calculated using the Sichardt criterion. The calculation has to be performed with the following equations[5]: and ^ = ^ / 15 (3) Qsc = A * (4) Where A Open area, equal to A=2nrwl, with rw being the screen/well radius in [m2] k Coefficient of hydraulic conductivity QSic Sichardt criterion flow through the open filters [m3/s] VSic Sichardt criterion screen entrance velocity A simple calculation for a 10-m-long "drive-in" filter and a hydraulic conductivity value of 2.88 10-6 m/s, combining equations (3) and (4), produces the following result: QSjc = 0.83 l/s = 49.5 l/min For these reasons a special guideline has to be made, which addresses the initial and final phase of the dewa-tering operation along the longwall works. Following the construction of a "drive-in" filter, its activation has to start slowly, enabling a controlled removal of the fine aquifer fractions from its close surroundings, creating an area of enhanced hydraulic conductivity, augmenting its effective radius and open area and establishing the required filter productivity. So, the "drive-in" filters will need to be fitted with a special device that prevents water rushing in from the sands. The valves also have to be used in the final stages of dewatering. Because of the mining works some "drive-in" filters will stop working, resulting in a flow-rate increase for those remaining. This procedure was already successfully applied to some "drive-in" filters [11]. It is also very interesting that, instead of all the planned "drive-in" filters along roadways, a smaller number can produce the same results. By inspecting Graph 2 and Fig. 8 we get the impression that not all the "drive-in" filters are necessary for the required draw down. However, the Sich- ardt criterion has to be taken into account, so their number, as used for modelling, is already low. This is seen in Table 2 and during the calculation of the Sichardt criterion. In fact, by modelling with the FEM it was possible to predict the zones in which the "drive-in" filters will produce the maximum dewatering. As an example, in Fig. 8 the maximum depression for the 230th day after the "drive-in" filters start up is shown, as a result of mathematical modelling. As can be seen, the main depression is around the "drive in" filters. According to the development of the water-table depression it can be seen that the dewatering due to the "drive-in" filters is successful (Fig. 8). However, by inspecting the groundwater levels on Fig. 8 and Graph 3 it is clear that after the pumping on the "drive-in" filters is stopped, the water level starts to rise, and in less than 200 days the water pressure in Pl1 would increase dramatically. This is an effect that should not be forgotten when planning to move from one longwall face to another. CONCLUSIONS Using the FEM methods a more accurate prediction of groundwater table was made. The drawdown predicted by the FEM is locally much higher than when using the FDM. From our previous experience [11] and from the literature, [8] and [2], it can be seen that the calculated drawdown using the FDM is more a function of the cell size than of the true flow situation. This can be avoided with a large decrease of cell size, but it would result in a FDM numerical instability and slow numerical computing. So, the FEM proved to be a better tool on the local scale. With the FEM the prediction of the flow rates on the "drive-in" filters is much better. The modelling shows that there is only one problem when using the "drivein" filters. The velocity of the flow has to be carefully controlled, not to exceed those two criteria used for the calculation. Accomplishing this goal is difficult, because the space for the drilling machines is limited and the working conditions for the drilling-team members are bad. So, instead of drilling long "drive-in" filters a careful activation technique, with a special device for water-rates control has to be used. Based on our working experience from the Velenje coal mine, from a few tests of the "drive-in" filter installations [11] [12], it can be said that the model's results are in accordance with reality. So, it can be concluded that the planned dewatering program has sense only if a special measure of activation will be used, and if the time to switch from one longwall face to another is short enough. references [1] Anderson, M. P., (2005). Heat as a ground water tracer. Ground Water, 43(6), 951-968 [2] Bear, J., Verruijt, A. (1992). Modelling Groundwater Flow and Pollution. D. Reidel Publishing Company, Dordrecht, Holland. [3] Čenčur Curk, B. and Witthuser, K. (2000). Field study of flow and transport from soil to unsaturated fractured rock V: Vlahovic, I.(ur.), Biondic, R. (ur.). 2. hrvatski geološki kongres, Cavtat - Dubrovnik, 17-20.05.2000 - Second Croatian Geological Congress, Cavtat - Dubrovnik, 17-20.05.2000. Zbornik radova. Zagreb: Institut za geološka istraživanja - Institute of Geology. [4] Domenico, A. P. and Schwartz, F. (1990). Physical and Chemical Hydrogeology. John Wiley & Sons, Canada. [5] Custodio, E. and Llamas, M. R. (2005). Idrologia sotterranea. Dario Flaccovio Editore, Italia. [6] Kruseman G. P. and de Ridder, N. A. (1994). Analysis and Evaluation of Pumping Test Data (Second Edition, Completely Revised). International Institute for Land Reclamation and Improvement, Netherlands. [7] Lajlar, B. and Supovec, I. (1997). Prognosis of the dewatering effects in the Pliocene aquifers in Velenje colliery using mathematical model. V: Veselič, M. (ur.), Norton, P. J. (ur.). Mine water and the environment: proceedings. Ljubljana: IRGO; [Granada]: IMWA, Vol. 1 [8] Marsily, de, G. (1986). Quantitative Hydrogeology. Groundwater Hydrology for Engineers, Masson, Editeur, Paris. [9] Mavec, M. and Supovec, I. (1998). Pliocene aquifer dewatering in Velenje coal mine and its effects on land subsidence. V: Norton, P. J. (ur.), Veselič, M. (ur.). Mine water and the environment: proceedings. Johannesburg: IMWA, 1, 75-86. [10] Supovec, I. and Veselič, M. (1989). Poročilo o modeliranju pliocenskih vodonosnikov v RLV. Geološki zavod Slovenije, Ljubljana. [11] Supovec, I., Lenart, M., and Jamnikar, S. (2006). Napoved učinkovitosti odvodnjevanja pliocenskih peskov z baražno progo po kadunji. Raziskovalno razvojna naloga, končno poročilo. [12] Veselič, M. and Supovec, I. (1986). Analiza odvod-njevanja pliocenskih krovninskih vodonosnikov s poskusnimi vodnjaki ob jamski progi na koti - 72 v letih 1984-1986. Geološki zavod Slovenije, Ljubljana. [13] Veličkovic, B. (2005). Colmatation as one of the processes in interaction between the ground water and surface water. Architecture and Civil Engineering, 3, 2, 165 - 172. [14] Vulic, M. and Uranjek, G. (2007). A contribution to construction monitoring with simultaneous application of various types of observations - Prispevek k spremljanju objektov s simultanimi meritvami različnih tipov. RMZ-mater. Geoenviron., 54, 2, 247-263. [15] Vukelič, Z., Šporin, J., and Vižintin, G. (2004). Pore pressure. RMZ-mater. geoenviron., 51, 4, 21172125. [16] Walton, C.W. (1970). Groundwater resource evaluation. McGraw-Hill B.C. NAPOVED OBNAŠANJA PILOTA MED DINAMIČNIM OBREMENILNIM TESTOM Z UPORABO TEHNOLOGIJE VGRAJENIH SENZORJEV DEFORMACIJ_ ANDREJ ŠTRUKELJ, MIRKO PŠUNDER, HELENA VRECL-KOJC Ln LUDVIK TRAUNER o avtorjih Andrej Štrukelj Univerza v Mariboru, Fakulteta za gradbeništvo Smetanova 17, 2000 Maribor, Slovenija E-pošta: andrej.strukelj@uni-mb.si Mirko Pšunder Univerza v Mariboru, Fakulteta za gradbeništvo Smetanova 17, 2000 Maribor, Slovenija E-pošta: mirko.psunder@uni-mb.si vodilni avtor Ludvik Trauner Univerza v Mariboru, Fakulteta za gradbeništvo Smetanova 17, 2000 Maribor, Slovenija E-pošta: trauner@uni-mb.si Helena Vrecl-Kojc Univerza v Mariboru, Fakulteta za gradbeništvo Smetanova 17, 2000 Maribor, Slovenija E-pošta: helena.vrecl@uni-mb.si Izvleček Na avtocestni povezavi Slivnica - Hajdina v bližini Maribora je bil na testnem pilotu izveden standardni dinamični obremenilni preizkus. Poleg ustaljenega postopka izvedbe preiskave je bila dodatno testirana na novo razvita tehnologija monitoringa namenjena spremljanju deformacij znotraj težko dostopnih betonskih konstrukcijskih elementov, še posebej pilotov in drugih vrst globokih temeljev. V testni pilot je bilo zato dodatno vgrajenih devet posebej izdelanih senzorjev za meritve specifičnih deformacij, ki so bili razporejeni na enakomernih razdaljah dveh metrov vzdolž osi pilota. Z omenjenimi senzorji smo v testnem pilotu spremljali deformacijsko stanje v betonu pilota med samim preizkusom ter na osnovi izmerjenih rezultatov ocenili normalne napetosti vzdolž osi pilota in trenje vzdolž plašča pilota ter normalne napetosti pod nogo pilota. Med preizkusom je bilo izvedenih sedem udarcev s prostopadno utežjo določene mase, pri čemer je bila utež spuščena z različnih višin. Pri vsakem od sedmih obtežnih primerov so bile med drugim izmerjene tudi specifične deformacije vzdolž celotne osi pilota povzročene s primarnim udarcem uteži ob pilot ter po odbojih uteži še s sekundarnim in terciarnim udarcem. Za vsakega od naštetih dogodkov so bile poiskane maksimalne vrednosti za vsa merilna mesta. Prav tako pa so bile po umiritvi sistema izvrednotene preostale oziroma trajne vrednosti specifičnih deformacij, ki so bile po prvem obtežnem primeru največje, pri ostalih šestih pa so bile bistveno manjše. Krivulje, ki povezujejo omenjene skupine maksimalnih vrednosti kažejo pričakovani trend upadanja vrednosti v smeri od vrha pilota proti nogi. Iz izmerjenih vrednosti so bile izvrednotene tudi največje vrednosti nateznih specifičnih deformacij, ki so posledica prehoda valovne fronte mimo posameznega senzorja. Ob upoštevanju elasticitetnega modula betona so bile izvrednotene normalne napetosti v betonu pilota vzdolž njegove osi. Na osnovi sprememb rezultante normalnih napetosti na območju med posameznimi merskimi mesti vzdolž pilota pa je bil ocenjen tudi potek rezultante strižnih napetosti oziroma specifičnega trenja ob plašču pilota. Ocena je bila narejena ob predpostavki konstantne debeline pilota. Rezultate meritev smo primerjali tudi z računalniško simulacijo napetostnega in deformacijskega stanja v pilotu in tleh, zato je bila neposredno ob pilotu izvedena tudi vrtina, ki je omogočila natančno klasifikacijo posameznih slojev zemljine in določitev njihovih mehanskih lastnosti. Ključne besede piloti, globoko temeljenje, dinamični obremenilni preizkus, tehnologija meritev deformacij, uporaba elasto-plastičnega modela zemljine, metoda končnih elementov PREDICTION OF THE PILE BEHAVIOUR UNDER DYNAMIC LOADING USING EMBEDDED STRAIN SENSOR TECHNOLOGY_ ANDREJ STRUKELJ, MIRKO PSUNDER, HELENA VRECL-KOJC and LUDVIK TRAUNER About the authors Andrej Štrukelj University of Maribor Faculty for Civil Engineering, Smetanova ulica 17, 2000 Maribor, Slovenia e-mail: andrej.strukelj@uni-mb.si Mirko Pšunder University of Maribor Faculty for Civil Engineering, Smetanova ulica 17, 2000 Maribor, Slovenia e-mail: mirko.psunder@uni-mb.si Helena Vrecl-Kojc University of Maribor Faculty for Civil Engineering, Smetanova ulica 17, 2000 Maribor, Slovenia e-mail: helena.vrecl@uni-mb.si did not have the purpose of developing an alternative method of pile loading tests. The presented monitoring technology proved itself as a very accurate and consistent. It gave in the first place the possibility of a closer look at the strains and stresses of the most unapproachable parts of different types of concrete structure elements especially piles and other types of deep foundations. Keywords piles, deep foundations, dynamic loading test, strain measurement technologies, elasto-plastic modelling, finite-elements method corresponding author Ludvik Trauner University of Maribor Faculty for Civil Engineering, Smetanova ulica 17, 2000 Maribor, Slovenia e-mail: trauner@uni-mb.si Abstract A standard dynamic loading test of the pile was performed on the highway section Slivnica - Hajdina near Maribor, Slovenia. Parallel to standard testing procedures the new monitoring technology based on specially developed strain sensors installed inside the pile body along the pile axis was introduced. On the basis of the measured results the normal strains along the pile axis were measured. Taking into consideration the elastic modulus of the concrete the normal stresses in the axial direction of the pile were also calculated and afterwards the shear stresses along the pile shaft have been estimated as well as the normal stresses below the pile toe. The estimation was made by considering a constant value for the pile diameter. The measured results were also compared with the computer simulation of the pile and the soil behaviour during all the successive test phases. The strain measurements inside the pile body during the standard dynamic loading test in present case 1 INTRODUCTION The bearing capacity and settlement of vertical loaded pile can be estimated using different methods [2, 6]. In order to identify the pile behaviour and its interaction with the surrounding soil layers, knowledge of the state of the strains along the pile axis inside the pile body is of essential importance. To make this possible, a chain of measurement points has to be included outside the pile structure and it should be fixed to the appropriate position on the pile reinforcement before the reinforcement is placed in the pile pit. After that the standard procedure for concreting the pile must be executed. The first opportunity for testing this idea was an estimation of the testing pile shaft's resistance, as described in the paper Štrukelj et. al. [8]. The basis of this estimation was the measurement of the normal strains of the pile in its axial direction at measuring points distributed over equal distances along the pile axis. These strains are proportional to the axial forces in the pile. When the course of the axial force along the pile axis is known, the resistance of the pile shaft can be estimated. The measurement points were distributed over equal distances of 1.0 m, starting 0.5 m above the pile toe. Since the strain gauges, electrical contacts and communication cables are very sensitive, any moisture and mechanical loading could be very harmful to their performance. Therefore, the measuring points were protected with special care. To ensure the mechanical protection of the chain of strain gauges together with the communication cables the whole measuring system was built inside a steel tube made of two standard C-profiles. The interior surface of the steel tube also served as the grounding surface to which the strain gages could be glued. For transportation reasons the tube was made of three segments that could be easily put together. The mutual connection of the segments was ensured by a special system of joints. The tubes external surface was degreased and made rough to ensure the adhesion with the concrete. Each strain gauge inside the tube was protected against the moisture. When the measuring tube was finished, it was connected with the pile reinforcement and put into the pit that was prepared for the pile concreting. The measurement system fulfilled the expectations and the measured results were accurate and stable. The only disadvantage of this measurement system was its non-flexibility. It could only be used to build the measurement chains where the measurement points are placed along one line and oriented in the same direction. The loading test included seven loading cases. Each represented the drop of a steel weight from a different height to the top of the test pile and the second and third hits after the first and second repulsions [9]. The location of the testing pile was selected on a construction site of the highway section Slivnica -Hajdina in Slovenia. The geological conditions of the testing pile's location and the positions of the nine sensors are shown in Fig.1, and described in Section 2.1. In the following sections the measuring equipment, its installation, the performance of dynamic loading test, and finally the evaluation of the measurement results is presented. In the last part of the paper the comparison of the measured results with results of numerical axis-symmetry analyses using the finite-element method (FEM) is presented. 2 PRELIMINARY WORKS ON THE TESTING SITE 2.1 GEOLOGICAL FIELD AND LABORATORY INVESTIGATIONS Geological conditions of the wider location of the testing site were acquired from the geotechnical report for the design of a highway crossover foundation [7]. The strength and deformability parameters were defined on the basis of field investigations by a standard penetra- tion test and probe measurements, as well as by laboratory testing of the samples taken by sounds of depth up to 25 m from seven different locations. The region belongs to the eastern part of the River Drava field that is mainly a plain area with only small differences in height. The relief evolution is based on the accumulation river-denudation with river influences in the past, which have deposited a layer, more than 10 m thick, of gravel GP, sandy gravel GP-GM, and some lenses of sandy clay CL with boulders over the Miocene base of sandy marl. On the location of the testing pile one additional sounding to the depth of 15.0 m was performed. The strength and deformability parameters at some depths of this additional sound were defined on the basis of field investigations by standard dynamic penetration tests, and also on the basis of the laboratory testing of samples. The ground water level in this sounding was encountered at approximately 3.30 m below the surface. The cross-section of the ground space with the pile (Fig. 1.) is composed of a thin, 1.0 m layer of embankment that was built on an original space of eight characteristic layers of different thicknesses. Fig. 1 shows the geological conditions of the testing pile's location and the positions of the nine sensors. Table 1 presents the strength parameters and the classification (0 is the internal friction of the soil, c is the cohesion of the soil, and Eoed is the oedometer elasticity modulus of the soil) of the layers presented in Fig. 1, which were determined on the basis of laboratory and field testing on additional sound samples. 2.2 MEASURING EQUIPMENT Besides the standard equipment that is needed for an estimation of the bearing capacity of the pile on the basis of a dynamical loading test [1], the additional specially produced strain sensors were built into the pile body. They were placed at equal distances of 2.0 m along the pile axis, starting at 1.5 m from the pile toe. The patented sensor design used for this purpose is very efficient, easy to build in, and robust enough to stand the water pressure and all the possible mechanical burdens during the pile concreting. The final solution for the strain sensor was to build a complete strain sensor for a single measurement point for all the layers of protection coatings and wiring. Such a sensor can be placed in the desired position in a very short time. It is insensitive to moisture and dust and its vital parts are very well protected against mechanical damage. The basis for such a sensor is a standard reinforcement bar of length about 150 cm and diameter 16 mm. In the falling height 1.0-3.0m V 3.30 o csi I o csi loading equipment pile head _ embankment CL with : jSG-9 lenses ML ______SM_________ SG-8 CL with lenses CH and SM 2,50 3,« B,70 18,10 pile toe Figure 1. Section of the testing site with the disposition of the soil layers and the measuring points on the pile. Table 1. The strength parameters of the soil layers. depth [m] [°] c [kPa] Eoed [kPa] soil classification 0.0-1.0 5.4/14.5 0.0 40000 embankment 1.0-2.5 23.0 20.0 1800 CL with lenses ML 2.5-3.4 25.4 0.0 2700 SM 3.4-5.3 24.7 20.0 19000 CL with lenses CH and SM 5.3-8.7 24.7 21.8 14000 CL-CH 8.7-10.7 33.0 0.0 11700 GM 10.7-15.0 38.8 0.0 43200 GM-GP 15.0-18.0 34.7 0.0 20700 GP 18.0 30.0 300.0 450000 marl middle of the reinforcement bar its surface is grinded on one or both sides, depending on the number of strain gages that should be installed. These strain gauges can be connected in a full, half, double-quarter or quarter Wheatstone bridge [4], [5]. The connecting cable is protected by a polyethylene tube and fixed to the reinforcement bar. The protection coating consists of two layers of polyurethane varnish, a layer of special silicone putty and a layer of permanently plastic sealant putty coated by aluminium foil. This combination of protection was tested and remained waterproof even at 30 m under the surface of water. The final layer represents the physical protection and can be made of polyethylene tube (Fig. 2) when the dimensions of the measured concrete elements are not too small to be significantly weakened by the built-in sensor. Otherwise, the physical protection of the measurement area can be made of cement mortar with the addition of an acrylic emulsion. This type of strain sensor proved to be very reliable, easy to build and cost effective, and can also be used together with the appropriate equipment and software [11] for the purposes of monitoring the other parts of structures. The successive phases of the sensors' preparation and their placing in the planned positions of the pile reinforcement are shown in Figs. 2-5. 2.3 CONSTRUCTION OF THE TESTING PILE The testing pile was constructed as a bored reinforced concrete pile [2] of 80 cm in diameter. The final length of the testing pile was 21.60 m. The length of the underground part of the pile was 19.60 m. As the thickness of the sandy clay, sandy gravel and gravel layer was about 18.0 m, the testing pile was embedded about 2.0 m into the marl base. After the pile was cleaned up to the level of the working area a 2.0 m long reinforced lengthening piece was additionally concreted on the top of the pile head (Fig. 6). Figure 6. The top of the testing pile before the installation of the loading equipment. Figure 8. The configuration of the measurement equipment during the loading test. 2.4 INSTALLATION OF THE DYNAMIC LOADING EQUIPMENT The gravitational weight was applied for loading, which was mounted on a steel guidance with an inner diameter of 80 cm and guiders with a length of more than 5.0 m (Fig. 7). The weight could fall from a height of 3.0 m. The assembly and lifting of the weight were managed by a 40 tons portable derrick. The dropping of the weight was performed with full gravity falling of a steel ballast weighing 10.5 tons on the head of the pile. The steel ballast was built in a 3.5 tons heavy steel frame (Fig. 7). At the top of the pile between the concrete surface and the falling steel ballast two wooden boards were placed to prevent the concrete being crushed. Figure 7. The testing site during the dynamic load test. Figure 9. Measured results displayed online on the computer screen during the loading test. 3 THE PERFORMANCE OF THE DYNAMIC LOADING TEST During the test the steel ballast with a mass of 10,500 kg was dropped seven times to the head of the pile from different heights (from 1.0 m to 3.0 m). The exact data about the height of the gravity falling of the steel ballast during each loading phase are presented in Table 2. During the test all the sensors built inside the pile body were connected to the data-acquisition unit (Fig. 8). Because of the explicitly dynamic nature of the measured events the chosen loading rate was 2400 Hz. The additional certainty of the measured results was achieved by installing two separate strain gauges in each sensor. Each strain gauge on each measuring point was connected independently to the data acquisition and the measured result in the individual measuring point is believed to be accurate when both strain gauges from each sensor give the same result. The measured results were stored on a computer hard disc during the measurement and they were also displayed online on the computer screen (Fig. 9) during the measurement procedure [9]. 4 EVALUATION OF THE MEASUREMENT RESULTS The review of the measured strains versus time for the three characteristic loading cases (1st, 4th and 7th) is presented in Fig. 10. From the sequence of three loading cases it can be seen that the amplitude of the strains increases with the increased height of dropping the steel ballast. The designated area of the last of three diagrams in Fig. 10 represents the most interesting part of the strain distribution during the primary hit of the steel ballast in the seventh loading case. This part of the measured signal is shown in the detailed view in Fig. 11 and represents the peak of the measured values during the entire test. From the time histories of the strains recorded during each loading case it is clear that the steel ballast after first hit to the top of the pile repulsed. The amplitude of the measured strains after the second hit after the first repulsion is about 25 percent of the strain amplitude of the primary hit. During each loading case the strains were recorded until the system had completely come to rest. All the measured signals were precisely analyzed afterwards. The result of this analysis is shown in Fig. 12, which represents the review of the peak strain values at all the measurement points for the first, second and third hits of the dropping mass for all the loading cases. Considering the value of the elasticity modulus for the pile's concrete is 34 GPa; the peak values of the normal stresses in the pile's axial direction (Fig. 13) were evaluated from the strains given in Fig. 12. After the equilibrium condition was used, and it was considered that the pile diameter remained constant along its whole length, the shear stresses along the pile shaft were also evaluated for all the loading cases. The results of this evaluation are presented in Fig. 14. The measurement results gave the values of the remaining strains after each loading case as well. The integrals of these values along the whole pile axis gave the values for the pile contraction for each loading case. Since the settlements of the pile head were geodetically measured after performing each loading case, the pile-toe settlements could also be evaluated. These values are presented in Table 2. Table 2. The evaluation of the pile-toe settlements. number of loading case height of gravity falling [m] pile contraction [mm] cumulative contraction [mm] pile-head settlements [mm] cumulative pile-head settlements [mm] pile-toe settlements [mm] 1 1.0 0.65 0.65 4.00 4.00 3.35 2 1.0 0.04 0.70 2.00 6.00 5.30 3 1.5 0.11 0.81 1.50 7.50 6.69 4 2.0 0.17 0.97 1.50 9.00 8.03 5 2.5 0.14 1.11 4.50 13.50 12.39 6 2.5 0.23 1.35 3.00 16.50 15.15 7 3.0 0.15 1.49 1.50 18.00 16.51 Figure 10. Measured strain signal recorded after the first, fourth and seventh (the last) falls of the dropping mass. In all the signals the strain peaks after the first hit of the dropping mass as well as the second and third hits after both repulsions can be clearly seen. -1300 J-0.07 Figure 11. Detailed view of the first peak of the strains in the measured signal of the seventh fall of the dropping mass. 0.0S 0.09 0.1 0.11 0.12 0.13 0.14 time :s] Figure 12. The review of the peak values of the measured results where the peak strain values at all the measurement points for the first, second and third hits of the dropping mass for all the loading cases are shown. In the figure the maximum tensile strains occurring during the hit-wave propagation along the pile are also presented. Figure 13. A review of the peak values of the normal stresses in the pile's axial direction calculated from the measured peak strain values at all the measurement points for the first, second and third hits of the dropping mass for all the loading cases. Figure 14. Curves of the peak values of the shear stresses on the pile shaft for all seven loading cases. The continuous lines were calculated from the measured results, whereas the dashed lines are extrapolated. 5 FINITE-ELEMENT ANALYSIS The results of field investigations were compared with a set of numerical analyses using the finite-element method (FEM). Axis-symmetric analyses of the interaction between a bored reinforced pile loaded with vertical dynamic loading and the ground were performed for this purpose. The analyses considered a 20.0 m long pile that is placed in a length of 18.0 m in eight different soil layers and below this to a depth of 2.0 m in a marl base. The ground water level is 3.30 m under the ground level. The reinforced concrete of the pile was considered to be linear elastic with a Young's modulus E = 34.0 GPa, a Poisson's ratio v = 0.2 and a unit weight y = 25 kN/m3. The strength properties of the ground (see Table 1) were determined on the basis of laboratory and field testing results of additional sound samples. This soil half-space was designed in numerical analyses by Hardening-Soil material model with isotropic hardening. This model considers nonlinear elastic hyperbolic dependence between stresses and strains; it enables the consideration of increasing soil yielding as a function of ground stresses, dilatation and cap yield surface, and is not based on theory of plasticity [3]. The parameters in the elasto-plastic Hardening-Soil model Eoedref = E50ref and Eurref = 3 E50ref where Eoedref is the tangent stiffness for the primary oedometer loading at the reference pressure, and Eurref is the unloading/reloading stiffness [3]. The cross-section of the analyzed model and the corresponding finite-element mesh is presented in Fig.15. Fig. 16 presents the absolute displacement after each loading phase of seven blows to the point A (at the pile head) and the point B (at the pile toe). Each loading phase is combined with two separate dynamic stages, where the first stage performs the dynamic force or stroke, and the second stage performs the relaxation time of 0.2 sec. This time was the relaxation period of the pile displacement also at the actual state according to the measurement results (see Fig. 10). The loading for the simulated case was calculated from the measured peak stress value shown in Fig.13. The results of the FEM analysis; the vectors of the vertical displacements after the last blow where umax = 22.70 mm at the pile head and umax = 19.80 mm at the pile toe (see Fig. 15). 0.00 5.00 10.00 15.00 20.00 25.00 30.00 Figure 15. Cross-section of the model with finite-element mesh and soil layers (Fig. left) and vectors of vertical displacements after the last blow (Fig. right). Figure 16. Absolute displacement after each blow of the point A (at the pile head) and point B (at the pile toe). 6 COMPARISON OF THE MEASUREMENT RESULTS WITH THE FINITE-ELEMENT ANALYSIS RESULTS The finite-element analyses of the axis-symmetrical FE model of pile and surrounding soil gave the results that are comparable with the results of methods based on measurements. Fig. 17 shows the relation between the vertical loading and the movements of the pile obtained with the analysis on the axis-symmetrical model compared to the measured values, and to the results of CAPWAP (Case Pile Wave Analysis Program) which separates static and damping soil characteristics and also allows an estimation of the side shear distribution and the pile's end bearing. CAPWAP is based on the wave equation model, which analyses the pile as a series of elastic segments and the soil as a series of elasto-plastic elements with damping characteristics, where the stiffness represents the static soil resistance and the damping represents the dynamic soil resistance. Typically the pile top force and velocity measurements acquired under high strain hammer impacts can be analyzed utilizing the signal matching procedure yielding forces and velocities over time and along the pile length. [1], [10]. i 15 20 25 4000 force F [kN] 8000 12000 16000 20000 24000 -CAPWAP-pile head ■ FEM-pile head FEM-pile toe ■ measured-pile head measured-pile toe Figure 17. Comparison of analyzed (FEM) and measured values of the vertical load (F) versus the settlement (s) of the pile head and the pile toe. The cumulative contraction of the pile after each blow is comparable. The difference between the pile-head and pile-toe settlement during the first three blows is higher, while after the third blow the difference between the measured and the analyzed results is decreasing. It can be concluded that the results of the settlements of the field measurements and the FEM analysis are within a comparable range. the results obtained by both measurement methods and the computer simulation is acceptable. In addition, the further analysis of the results shows that the agreement of the calculated and measured values is even better after using some parameters for the soil model employed in a computer simulation and obtained from the analysis of the measured results. 7 CONCLUSIONS The present paper describes the development and application of a new measuring technology that proved itself to be very useful for investigating the behaviour of piles and also other embedded structures as well as the surrounding soil loaded by different kinds of static and dynamic loads. The presented results show that the state of the strains inside the pile body can be obtained in a very accurate and cost-efficient way even for events that show an extremely dynamic nature. In addition, using the described measurement method with the appropriate distribution of measurement points and measurement directions, the most complex states of strain can be followed. In our paper the presented example shows that the extreme values of the measured strains, and the calculated stresses based on these strains, are significantly higher than the values given in an official report of the standard dynamic pile test. However, it should be considered that the values obtained on the basis of the strain measurements are the peak values of the dynamic response and not the reduced values used for the simulation of the static behaviour of the pile on the basis of a dynamic test. Additionally, it should be stated that the purpose of the demonstrated measurement technology is in the first place to create the possibility of a closer look at the behaviour of the most unapproachable parts of essential construction elements, like piles and other types of deep foundations. The parallel measurements of the strains inside the pile body during the standard dynamic loading test therefore did not have the purpose of developing an alternative measuring procedure for the dynamic loading test. In this case the standard loading test served only as an opportunity for testing the applicability, accuracy and the consistency of the new measurement technology. Even though the conclusions were written before, the results of the field investigations were compared with a set of numerical analyses using the finite-element method (FEM). In the presented computer simulation the model with the axial symmetry was used for the interaction between a bored, reinforced pile loaded dynamically in the vertical direction. The comparison of ACKNOWLEDGMENT The present research work is a part of a project supported by EUREKA, a pan-European network founded to enhance European competitiveness through its support to businesses, research centres and universities, which carry out pan-European projects to develop innovative products, processes and services. references [1] ASTM (1994). Standard test for high-strain dynamic testing of piles, Designation: D4945-89. Philadelphia. [2] Bowles, J.E. (1996). Foundation analysis and design. McGraw -Hill, New York. [3] Brinkgrave, R.B.J., Vermeer, P.A. (1998). PLAXIS 2D users manual, version 8. Balkema, Rotterdam. [4] Hoffman, K. (1989). An Introduction to Measurements using Strain Gages. Hottinger Baldwin Messtechnik GmbH, Darmstadt. [5] Hoffman, K. (1996). Hinweise zum Applizieren von Dehnungsmeßstreifen (DMS), 4. erweiterte Fassung. Hottinger Baldwin Messtechnik GmbH, Darmstadt. [6] Škrabl, S. (2002). Bearing capacity and settlement of vertically-loaded piles. Proc. Interntl. Deep Foundations Congress, Geotechnical special technical publication, 116, Reston, VA: Geo-Institute, American Society of Civil Engineers, Orlando, Florida, 53-63. [7] Štern, K. (2007). Geološko-geotehnično poročilo, AC odsek Slivnica - Draženci - Gruškovje, arhiv.št. 50-3065/06. Geoinženiring d.o.o., Ljubljana. [8] Štrukelj, A., Škrabl, S., Štern, K., and Logar, J. (2005). The assesment of pile shaft resistance based on axial strain measurements during the loading test. Actageotechnica Slovenica, 2(2), 12-23. [9] Trauner, L. and Štrukelj, A. (2007). Poročilo o rezultatih meritev specifičnih deformacij v testnem pilotu na 0091 Slivnica-Hajdina sklop A od km 0+000 do km 2+000 med izvedbo dinamičnega obremenilnega preiskusa, opr.št. 006/07-C0TEG-AŠ, Faculty of Civil Engineering, University of Maribor, Maribor. [10] Strniša, G. (2007). Poročilo o meritvi nosilnosti uvrtanega pilota, opr.št.136-01-2007. SLP d.o.o., Ljubljana. [11] Rebolj, D., Čuš Babič, N., Magdič, A., Podbreznik, P., and Pšunder, M. (2008). Automated construction activity monitoring system. Advanced engineering informatics, 22, 4. We wish to invite you to our traditional, this year already the 10th LUJO ŠUKLJE MEMORIAL DAY This year's international lecturer is Prof. Desai, a worldwide renowned expert. Invited national lecturer, Dr. Logar, Slovenian expert and teacher, invited designers and contractors to present their topics about geotechnics in the development of the Port of Koper. It is our wish to get together at the 10th Lujo Suklje memorial day in the largest possible number and in this way pay tribute to the memory of a great man, Prof. Suklje, and show our loyalty to the field of science that gives us our daily bread. 10th Lujo Suklje memorial day will take place in Hotel Kokra, at the Protocol Centre Brdo pri Kranju, on Friday, 25 September 2009. Conference programme 13.00- 14.00 Participant registration 14.00- 14.30 Opening 14.30- 15.30 Constitutive Modelling and Computer Methods in Geotechnical Engineering, Prof. Dr. Chandra S. Desai, University of Arizona, Tucson, USA 15:30- 16:15 Break for coffee and snack 16.15- 17.45 Geotechnics in the development of the Port of Koper, Assist Prof. Dr. Janko Logar, University of Ljubljana, FGG and guests: Lilian Battelino, Institute for Water of RS, Gorazd Strniša, SLP and Martin Pregelj, Primorje 18.00- 18.45 Guided tour through the Protocol Centre Brdo 19.30 Festive dinner in the restaurant of Hotel Kokra Organising Committee Chairperson: Prof. Dr. Ludvik Trauner Members: Mr. Gorazd Humar Dr. Stanislav Lenart, Secretary of SloGeD Dr. Ana Petkovšek, President of SloGed Address ORGANISING COMMITTEE OF THE 10th DAY OF LUJO SUKLJE University of Ljubljana Faculty of Civil and Geodetic Engineering Jamova 2 1000 LJUBLJANA Tel.: +386 1 476 85 41, Fax.: +386 1 425 06 81 E-mail: ana.petkovsek@fgg.uni-lj.si stanislav.lenart@zag.si # NAVODILA AVTORJEM Članki so objavljeni v angleškem jeziku s prevodom izvlečka v slovenski jezik. VSEBINA ČLANKA Članek naj bo napisan v naslednji obliki: - Naslov, ki primerno opisuje vsebino članka in ne presega 80 znakov. - Izvleček, ki naj bo skrajšana oblika članka in naj ne presega 250 besed. Izvleček mora vsebovati osnove, jedro in cilje raziskave, uporabljeno metodologijo dela, povzetek izidov in osnovne sklepe. - Uvod, v katerem naj bo pregled novejšega stanja in zadostne informacije za razumevanje ter pregled izidov dela, predstavljenih v članku. - Teorija. - Eksperimentalni del, ki naj vsebuje podatke o postavitvi preiskusa in metode, uporabljene pri pridobitvi izidov. - Izidi, ki naj bodo jasno prikazani, po potrebi v obliki slik in preglednic. - Razprava, v kateri naj bodo prikazane povezave in posplošitve, uporabljene za pridobitev izidov. Prikazana naj bo tudi pomembnost izidov in primerjava s poprej objavljenimi deli. - Sklepi, v katerih naj bo prikazan en ali več sklepov, ki izhajajo iz izidov in razprave. - Literatura, ki mora biti v besedilu oštevilčena zaporedno in označena z oglatimi oklepaji [1] ter na koncu članka zbrana v seznamu literature. OBLIKA ČLANKA Besedilo naj bo pisano na listih formata A4, z dvojnim presledkom med vrstami in s 3.0 cm širokim robom, da je dovolj prostora za popravke lektorjev. Najbolje je, da pripravite besedilo v urejevalniku Microsoft Word. Hkrati dostavite odtis članka na papirju, vključno z vsemi slikami in preglednicami ter identično kopijo v elektronski obliki. Enačbe naj bodo v besedilu postavljene v ločene vrstice in na desnem robu označene s tekočo številko v okroglih oklepajih. ENOTE IN OKRAJŠAVE V besedilu, preglednicah in slikah uporabljajte le standardne označbe in okrajšave SI. Simbole fizikalnih veličin v besedilu pišite poševno (npr. v, T itn.). Simbole enot, ki sestojijo iz črk, pa pokončno (npr. Pa, m itn.). Vse okrajšave naj bodo, ko se prvič pojavijo, izpisane v celoti. SLIK6 Slike morajo biti zaporedno oštevilčene in označene, v besedilu in podnaslovu, kot sl. 1, sl. 2 itn. Posnete naj bodo v kateremkoli od razširjenih formatov, npr. BMP, JPG, GIF. Za pripravo diagramov in risb priporočamo CDR format (CorelDraw), saj so slike v njem vektorske in jih lahko pri končni obdelavi preprosto povečujemo ali pomanjšujemo. Pri označevanju osi v diagramih, kadar je le mogoče, uporabite označbe veličin (npr. v, T). V diagramih z več krivuljami mora biti vsaka krivulja označena. Pomen oznake mora biti razložen v podnapisu slike. Za vse slike po fotografskih posnetkih je treba priložiti izvirne fotografije ali kakovostno narejen posnetek. PREGLEDNICE Preglednice morajo biti zaporedno oštevilčene in označene, v besedilu in podnaslovu, kot preglednica 1, preglednica 2 itn. V preglednicah ne uporabljajte izpisanih imen veličin, ampak samo ustrezne simbole. K fizikalnim količinam, npr. t (pisano poševno), pripišite enote (pisano pokončno) v novo vrsto brez oklepajev. Vse opombe naj bodo označene z uporabo dvignjene številke1. SEZNAM LITERATURE Vsa literatura mora biti navedena v seznamu na koncu članka v prikazani obliki po vrsti za revije, zbornike in knjige: [1] Feng, T. W. (2000). Fall-cone penetration and water content ralationship of clays. Geotechnique 50, No. 2, 181-187. [2] Ortolan, Ž. and Mihalinec, Z. (1998). Plasticity index-Indicator of shear strength and a major axis of geotechnical modelling. Proceedings of the Eleventh Danube-European conference on soil mechanics and geotechnical engineering, Poreč, 25 -29 May 1998. [3] Toporišič, J. (1994). Slovenski pravopis. 2nd.ed., DZS, Ljubljana. PODATKI O AVTORJIH Članku priložite tudi podatke o avtorjih: imena, nazive, popolne poštne naslove, številke telefona in faksa, naslove elektronske pošte. Navedite kontaktno osebo. SPREJEM ČLANKOV IN AVTORSKE PRAVICE Uredništvo si pridržuje pravico do odločanja o sprejemu članka za objavo, strokovno oceno mednarodnih recenzentov in morebitnem predlogu za krajšanje ali izpopolnitev ter terminološke in jezikovne korekture. Avtor mora predložiti pisno izjavo, da je besedilo njegovo izvirno delo in ni bilo v dani obliki še nikjer objavljeno. Z objavo preidejo avtorske pravice na revijo ACTA GEOTECHNICA SLOVENICA. Pri morebitnih kasnejših objavah mora biti AGS navedena kot vir. Rokopisi člankov ostanejo v arhivu AGS. Vsa nadaljnja pojasnila daje: Uredništvo ACTA GEOTECHNICA SLOVENICA Univerza v Mariboru Fakulteta za gradbeništvo Smetanova ulica 17 2000 Maribor Slovenija E-pošta: ags@uni-mb.si INSTRUCTIONS FOR AUTHORS The papers are published in English with a translation of the abstract into Slovene. FORMAT OF THE PAPER The paper should have the following structure: - A Title that adequately describes the content of the paper and should not exceed 80 characters; - An Abstract, which should be viewed as a mini version of the paper and should not exceed 250 words. The Abstract should state the principal objectives and the scope of the investigation and the methodology employed, it should also summarise the results and state the principal conclusions; - An Introduction, which should provide a review of recent literature and sufficient background information to allow the results of the paper to be understood and evaluated; - A Theoretical section; - An Experimental section, which should provide details of the experimental set-up and the methods used for obtaining the results; - A Results section, which should clearly and concisely present the data using figures and tables where appropriate; - A Discussion section, which should describe the relationships shown and the generalisations made possible by the results and discuss the significance of the results, making comparisons with previously published work; - Conclusions, which should present one or more conclusions that have been drawn from the results and subsequent discussion; - References, which must be numbered consecutively in the text using square brackets [1] and collected together in a reference list at the end of the paper. LAYOUT OF THE TEXT The text should be written in A4 format, with double spacing and margins of 3 cm, to provide editors with space to write in their corrections. Microsoft Word for Windows is the preferred format for submission. One hard copy, including all figures, tables and illustrations and an identical electronic version of the manuscript must be submitted simultaneously. Equations should be on a separate line in the main body of the text and marked on the right-hand side of the page with numbers in round brackets. UNITS AND ABBREVIATIONS Only standard SI symbols and abbreviations should be used in the text, tables and figures. Symbols for physical quantities in the text should be written in Italics (e.g. v, T, etc.). Symbols for units that consist of letters should be in plain text (e.g. Pa, m, etc.). All abbreviations should be spelt out in full on first appearance. FIGURES Figures must be cited in consecutive numerical order in the text and referred to in both the text and the caption as Fig. 1, Fig. 2, etc. Figures may be saved in any common format, e.g. BMP, JPG, GIF. However, the use of CDR format (CorelDraw) is recommended for graphs and line drawings, since vector images can be easily reduced or enlarged during final processing of the paper. When labelling axes, physical quantities (e.g. v, T) should be used whenever possible. Multi-curve graphs should have individual curves marked with a symbol; the meaning of the symbol should be explained in the figure caption. Good quality black-and-white photographs or scanned images should be supplied for illustrations. [3] Toporišič, J. (1994). Slovenski pravopis. 2nd.ed., DZS, Ljubljana. AUTHOR INFORMATION The following information about the authors should be enclosed with the paper: names, complete postal addresses, telephone and fax numbers and E-mail addresses. Indicate the corresponding person. ACCEPTANCE OF PAPERS AND COPYRIGHT The Editorial Committee of the Slovenian Geotechnical Review reserves the right to decide whether a paper is acceptable for publication, to obtain peer reviews for submitted papers, and if necessary, to require changes in the content, length or language. Authors must also enclose a written statement that the paper is original unpublished work, and not under consideration for publication elsewhere. On publication, copyright for the paper shall pass to the ACTA GEOTECHNICA SLOVENICA. The AGS must be stated as a source in all later publication. Papers will be kept in the archives of the AGS. TABLES Tables must be cited in consecutive numerical order in the text and referred to in both the text and the caption as Table 1, Table 2, etc. The use of names for quantities in tables should be avoided if possible: corresponding symbols are preferred. In addition to the physical quantity, e.g. t (in Italics), units (normal text), should be added on a new line without brackets. Any footnotes should be indicated by the use of the superscript1. For further information contact: Editorial Board ACTA GEOTECHNICA SLOVENICA University of Maribor Faculty of Civil Engineering Smetanova ulica 17 2000 Maribor Slovenia E-mail: ags@uni-mb.si LIST OF REFERENCES References should be collected at the end of the paper in the following styles for journals, proceedings and books, respectively: [1] Feng, T. W. (2000). Fall-cone penetration and water content ralationship of clays. Geotechnique 50, No. 2, 181-187. [2] Ortolan, Ž. and Mihalinec, Z. (1998). Plasticity index-Indicator of shear strength and a major axis of geotechnical modelling. Proceedings of the Eleventh Danube-European conference on soil mechanics and geotechnical engineering, Poreč, 25 -29 May 1998. NAMEN REVIJE Namen revije ACTA GEOTECHNICA SLOVENICA je objavljanje kakovostnih teoretičnih člankov z novih pomembnih področij geomehanike in geotehnike, ki bodo dolgoročno vplivali na temeljne in praktične vidike teh področij. AIMS AND SCOPE ACTA GEOTECHNICA SLOVENICA aims to play an important role in publishing high-quality, theoretical papers from important and emerging areas that will have a lasting impact on fundamental and practical aspects of geomechanics and geotechnical engineering. ACTA GEOTECHNICA SLOVENICA objavlja članke s področij: mehanika zemljin in kamnin, inženirska geologija, okoljska geotehnika, geosintetika, geotehnične konstrukcije, numerične in analitične metode, računalniško modeliranje, optimizacija geotehničnih konstrukcij, terenske in laboratorijske preiskave. Revija redno izhaja dvakrat letno. ACTA GEOTECHNICA SLOVENICA publishes papers from the following areas: soil and rock mechanics, engineering geology, environmental geotechnics, geosynthetic, geotechnical structures, numerical and analytical methods, computer modelling, optimization of geotechnical structures, field and laboratory testing. The journal is published twice a year. AVTORSKE PRAVICE COPYRIGHT Ko uredništvo prejme članek v objavo, prosi avtorja(je), da prenese(jo) avtorske pravice za članek na izdajatelja, da bi zagotovili kar se da obsežno razširjanje informacij. Naša revija in posamezni prispevki so zaščiteni z avtorskimi pravicami izdajatelja in zanje veljajo naslednji pogoji: Fotokopiranje V skladu z našimi zakoni o zaščiti avtorskih pravic je dovoljeno narediti eno kopijo posameznega članka za osebno uporabo. Za naslednje fotokopije, vključno z večkratnim fotokopiranjem, sistematičnim fotokopiranjem, kopiranjem za reklamne ali predstavitvene namene, nadaljnjo prodajo in vsemi oblikami nedobič-konosne uporabe je treba pridobiti dovoljenje izdajatelja in plačati določen znesek. Naročniki revije smejo kopirati kazalo z vsebino revije ali pripraviti seznam člankov z izvlečki za rabo v svojih ustanovah. Elektronsko shranjevanje Za elektronsko shranjevanje vsakršnega gradiva iz revije, vključno z vsemi članki ali deli članka, je potrebno dovoljenje izdajatelja. Upon acceptance of an article by the Editorial Board, the author(s) will be asked to transfer copyright for the article to the publisher. This transfer will ensure the widest possible dissemination of information. This review and the individual contributions contained in it are protected by publisher's copyright, and the following terms and conditions apply to their use: Photocopying Single photocopies of single articles may be made for personal use, as allowed by national copyright laws. Permission of the publisher and payment of a fee are required for all other photocopying, including multiple or systematic copying, copying for advertising or promotional purposes, resale, and all forms of document delivery. Subscribers may reproduce tables of contents or prepare lists of papers, including abstracts for internal circulation, within their institutions. Electronlc Storage Permission of the publisher is required to store electronically any material contained in this review, including any paper or part of the paper. ODGOVORNOST Revija ne prevzame nobene odgovornosti za poškodbe in/ali škodo na osebah in na lastnini na podlagi odgovornosti za izdelke, zaradi malomarnosti ali drugače, ali zaradi uporabe kakršnekoli metode, izdelka, navodil ali zamisli, ki so opisani v njej. RESPONSIBILITY No responsibility is assumed by the publisher for any injury and/or damage to persons or property as a matter of product liability, negligence or otherwise, or from any use or operation of any methods, products, instructions or ideas contained in the material herein.