<^2 c. s. Desai constituive modeling and computer methods in geotechnicrl engineering g. jeromel et al. an analysis op the geomechrnical processes in coal mining using the velenje mining method M.jesmanl et al. effect of plasticity and normal staess on the undrained shear modulus op clayey soils j. floser 6 n. Gosar determinrtion of vs3o for seismic ground clrssificrtion in the ljubljana area, slovenia ISSN: 1854-0171 flCTR GGOTGCHNICfl SLOVGNICn ustanovitelji Founders uredniški odbor editorial Board 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 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 Univerza v Mariboru, Fakulteta za gradbeništvo University of Maribor, Faculty of Civil Engineering odgovorni urednik editor-in-chief Ludvik Trauner University of Maribor uredniki co-editors Bojana Dolinar University of Maribor Borut Macuh University of Maribor Stanislav Škrabl University of Maribor Helena Vrecl Kojc University of Maribor Darinka Battelino UUniversity of Trieste Heinz Brandl Vienna University of Technology József Farkas Budapest University of Technology and Economics Theodoros Hatzigogos Aristotle University of Thessaloniki Rolf Katzenbach Technical University Darmstadt Zlatko Langof University of Sarajevo Jakob Likar University of Ljubljana Janko Logar University of Ljubljana Bojan Majes University of Ljubljana Milan Maksimovic University of Belgrade Borut Petkovšek Slovenian National Building and Civil Engineering Institute Mihael Ribičič University of Ljubljana César Sagaseta University of Cantabria Patrick Selvadurai McGill University Stephan Semprich University of Technology Graz Abdul-Hamid Soubra University of Nantes Küchi Suzuki Saitama University Antun Szavits-Nossan University of Zagreb Ivan Vaniček Czech Technical University in Prague Chandrakant. S. Desai University of Arizona Pedro Seco e Pinto National Laboratory of Civil Engineering Address Lektor 1 proof-Reader Paul McGuiness Naklada circulation Naslov uredniš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 500 izvodov - issues spletni naslov ueb Address print Uni Založba d.o.o. Tisk izdajatelj publisher Bojan Zlender University of Maribor posvetovalni uredniki Advisory editors http://www. fg.uni-mb.si/journal-ags 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. 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 CONTENTS fr !V Ludvik Trauner UVODNIK Chandrakant s. Desai KONSTITUTIVNO MODELIRANJE IN RAČUNALNIŠKE METODE V GEOTEHNIKI Gregor jeromeL in drugi ANALIZA GEOMEHANSKIH PROCESOV ODKOPAVANJA PREMOGA Z VELENJSKO ODKOPNO METODO Mehrab jesmani in drugi VPLIV PLASTIČNOSTI IN NORMALNE NAPETOSTI NA NEDRENIRANE STRIŽNE MODULE GLINENIH ZEMLJIN [óoj janez Rošer in nndrej Gosar DOLOČITEV VS30 ZA SEIZMIČNO KLASIFIKACIJO TAL NA OBMOČJU LJUBLJANE 3 NAVODILA AVTORJEM Ludvik Trauner EDITORIAL chandrakant s. Desal CONSTITUTIVE MODELING AND COMPUTER METHODS IN GEOTECHNICAL ENGINEERING A Gregor jeromeL et aL. AN ANALYSIS OF THE GEOMECHANICAL PROCESSES IN COAL MINING USING THE VELENjE MINING METHOD E Mehrab jesmanl et aL. EFFECT OF PLASTICITY AND NORMAL STRESS ON THE UNDRAINED SHEAR MODULUS OF CLAYEY SOILS E janez Roser and Andrej Gosar DETERMINATION OF VS30 FOR SEISMIC GROUND CLASSIFICATION IN THE LJUBLJANA AREA, SLOVENIA INSTRUCTIONS FOR AUTHORS E L9 UVODNIK Acta Geotechnica Slovenica je po šestih letih izhajanja postala odmevna mednarodna revija, ki jo indeksirajo priznane baze podatkov: SCIE, ICONDA in GeoRef. Z namenom, da bi še izboljšali kvaliteto njenih vsebin, smo v letošnjem letu razširili člane uredniškega odbora s svetovno priznanimi strokovnjaki. Posebno veseli in ponosni smo, da so se nam pridružili posvetovalni uredniki Darinka Battelino, Heinz Brandl, Chandrakant S. Desai in Pedro Seco e Pinto in člani uredniškega odbora Patrick Salvadurai, Kuchi Suzuki in Antun Szavits-Nossan. Vsem se iskreno zahvaljujemo za sodelovanje. Regent's profesor dr. Chandrakant S. Desai z Oddelka za gradbeništvo in mehaniko Univerze v Arizoni, Tucson, objavlja v pričujoči izdaji zanimiv pregledni članek z naslovom »Konstitutivno modeliranje in računalniške metode v geotehniki«, ki povzema njegovo bogato znanje in izjemne prispevke v bazičnih in aplikativnih raziskavah na področju konstitutivnega modeliranja materialov, laboratorijskega preizkušnja in računalniških metod v gradbeništvu, povezanih z geomehaniko, geotehniko in konstrukcijsko mehaniko/konstrukcijskim inženiringom. Vsebino članka je posvetil prof. Luju Šukljetu, slovenskemu pionirju za mehaniko tal. Temo tega prispevka je dr. Desai predstavil na jubilejnem 10. Šukljetovem dnevu (2009), kjer je prejel tudi najvišje priznanje Slovenskega geotehniškega društva. Avtorji Gregor Jeromel, Milan Medved in Jakob Likar v članku predstavljajo numerični model, ki omogoča poglobljene analize geomehanskih procesov v krovnini, talnini ter premoškem sloju pri časovnem razvoju podetažnega odkopavanja premoga, ima širši okvir uporabnosti in je izjemnega pomena za ugotavljanje intenzivnosti ter obsega rušnih procesov pri podetažnem odkopavanju premoga, za realna načrtovanja odkopavanja premoga v premogovniku v naslednjih desetletjih ob povečani varnosti rudarjev. Mehrab Jesmani, Hamed Faghihi Kashani in Mehrad Kamalzare v članku prikazujejo rezultate direktnih nedreniranih strižnih preiskav porušenih glin, ki so bile izvedene z namenom proučevanja vpliva indeksa plastičnosti in normalne napetosti na njihove strižne lastnosti in strižni modul. Avtorja Janez Rošer in Andrej Gosar obravnavata raziskave povprečne hitrosti strižnega valovanja v vrhnjih 30 m sedimentov (Vs30) na območju Ljubljane za seizmično klasifikacijo tal. To območje sodi med seizmično najbolj aktivne v Sloveniji in ga je v preteklosti prizadelo več večjih potresov. Obstoječe potresne mikrorajonizacije v Ljubljanskem bazenu, zaradi pomanjkanja geofizikalnih podatkov in podatkov iz vrtin, niso zanesljive. Ugotovljena prostorska porazdelitev Vs30 zato predstavlja pomembno dopolnitev že obstoječe mikrorajonizacije. Ludvik Trauner Glavni urednik EDITORIAL After six years of publication, Acta Geotechnica Slovenica has become an established international journal that is indexed by important databases: SCIE, ICONDA and GeoRef. With the aim of improving the quality of its contents, we increased the number of editorial members with worldwide recognized experts during this year. We are especially glad and proud that Advisory Editors Darinka Battelino, Heinz Brandl, Chandrakant S. Desai and Pedro Seco e Pinto as well as Members of the Editorial Board Patrick Salvadurai, Kuchi Suzuki and Antun Szavits-Nossan have joined us. We want to thank them all for their cooperation. In this issue, regent's professor dr. Chandrakant S. Desai from the Department of Civil Engineering and Engineering Mechanics, the University of Arizona, Tucson, Arizona, USA, publishes an interesting review article with the title Constitutive Modeling and Computer Methods in Geotechnical Engineering, which summarizes his rich knowledge and exceptional contributions with regard to basic and applicative research in the field of the constitutive modelling of materials, laboratory testing and computer methods in construction in connection with geomechanics, geotechnics and constructional mechanics/constructional engineering. The content of his article is dedicated to professor Lujo Suklje, the Slovenian pioneer of soil mechanics. Dr. Desai presented the theme of his contribution at the jubilee 10. Suklje Day (2009), where he was awarded with the highest honour of the Slovenian Geotechnical Society. The authors Gregor Jeromel, Milan Medved and Jakob Likar introduce in their article a numerical model that allows for in-depth analyses of the geomechanical processes occurring in the hanging wall, the footwall and in the coal seam during sub-level coal excavation The model is broadly applicable and highly relevant for analysing the intensity and the level of the caving processes in sub-level coal mining as well as for making realistic plans for coal excavation with workers' safety in mind. In their contribution, Mehrab Jesmani, Hamed Faghihi Kashani and Mehrad Kamalzare introduce the results of direct shear tests performed in order to study the effect of the plasticity index and the normal stress on the shear behaviour and the shear modulus of remoulded clays. The authors Janez Roser and Andrej Gosar deal with research on the average shear-wave velocity in the upper 30-m sediments (Vs30) for a seismic classification of the ground in the Ljubljana region. This area belongs to one of the most seismically active in Slovenia, and as a result it was subject to some major earthquakes in the past. The existing micro-zonation studies in the Ljubljana basin are inadequate, since there is a lack of borehole and geophysical data. So, the determined space distribution Vs30 represents an important supplement to the already-existing microzonation. Ludvik Trauner KONSTITUTIVNO MODELIRANJE IN RAČUNALNIŠKE METODE V GEOTEHNIKI CHANDRAKANT S. DESAI o avtorju Chandrakant S. Desai The University of Arizona, Department of Civil Engineering and Engineering Mechanics Tucson, Arizona. ZDA E-posta: csdesai@email.arizona.edu Izvleček Pri postopkih analize in projektiranja geotehničnih problemov so v ospredju računalniške metode. Konstitutivni modeli, ki podajajo obnašanje geoloških materialov in stikov mejnih ploskev, imajo pomembno vlogo pri rešitvah, pridobljenih z uporabo računalniških metod ali drugih postopkov pridobivanja rešitev. O konstitutivnih kot tudi računalniških modelih obstaja obširna literatura, ta članek pa obravnava koncept porušnega stanja ("disturbed state concept" - DSC) za konstitutivno modeliranje in metodo končnih elementov za računalniške rešitve. Koncept porušnega stanja (DSC), poenoten in hierarhičen pristop, določa poenoten model za podajanje obnašanja geoloških materialov in stikov mejnih ploskev. Pomembni dejavniki, kot so elastično in plastično območje ter območje lezenja, obremenilna stanja sprememba volumna (krčenje in raztezanje), porušitev (mehčanje in poškodbe ali utrjevanje), termični učinki, delno nasičenje in likvifakcija, so lahko vključeni v isti model DSC. Zaradi njegove hierarhične narave, lahko iz DSC izvedemo poenostavljene modele za določene aplikacije. Uspešno je bil uporabljen za definiranje obnašanja mnogih geoloških materialov in stikov mejnih ploskev. Razviti so bili postopki za določitev parametrov za DSC modele, ki temeljijo na laboratorijskih preizkusih. Na nivoju vzorcev so bili glede na podatke laboratorijskih preizkusov potrjeni različni modeli iz DSC. Potrjeni so bili na ravni problema praktične mejne vrednosti, s primerjavo med opaženim obnašanjem na terenu in/ali simuliranimi problemi v laboratoriju ter napovedmi ob uporabi računalniških postopkov (končnih elementov), ki so bili dopolnjeni z DSC, predstavljeni pa so bili v različnih publikacijah s strani Dr. Desaia in sodelavcev, kot je navedeno v referencah. V članek so vključeni trije tipični primeri tovrstnih potrditev na nivoju praktičnih problemov. DSC lahko poda poenotene in učinkovite modele za širok spekter geomehanskih in drugih inženirskih materialov ter stikov mejnih ploskev. Ključne besede konstitutivno modeliranje, koncept porušnega stanja - DSC, geološki materiali, stiki mejnih ploskev posvetilo Prof. Lujo Šuklje je bil pionir mehanike tal v Sloveniji, z mednarodno priznanimi raziskovalnimi deli. Pomembno je prispeval k razumevanju reoloških lastnosti tal, vključno z lezenjem in anizotropijo, visokoplastičnimi modeli, konsolidacijo in temeljenjem. Njegova knjiga Reološki vidiki mehanike tal, objavljena leta 1969, je pustila neizbrisen pečat v inženirski literaturi. Ta članek posvečam, z vsem strokovnim spoštovanjem, spominu na prof. Šukljeta. CONSTITUTIVE MODELING AND COMPUTER METHODS IN GEOTECHNICAL ENGINEERING CHANDRAKANT S. DESAI About the author Chandrakant S. Desai The University of Arizona, Department of Civil Engineering and Engineering Mechanics Tucson, Arizona. USA E-mail: csdesai@email.arizona.edu and powerful models for a wide range ofgeomechanical and other engineering materials, and interfaces/joints. Abstract Computer methods are in the forefront of the procedures for analysis and design for geotechnical problems. Constitutive models that characterize the behavior of geologic materials and interfaces/joints play a vital role in the solutions obtained by using computer methods or any other solution procedure. The literature on both constitutive and computer models is wide; attention in this paper is devoted to the disturbed state concept (DSC) for constitutive modeling and the finite element method for computer solutions. The disturbed state concept, a unified and hierarchical approach, provides a unified framework for characterization of the behavior of geologic materials and interfaces/joints. Important factors such as elastic, plastic and creep responses, stress paths, volume change (contraction and dilation), disturbance (softening and damage or stiffening), thermal effects, partial saturation and liquefaction can be included in the same DSC framework. Because of its hierarchical nature, simplified models for specific applications can be derived from the DSC. It has been applied successfully for defining behavior of many geologic materials and interfaces/joints. Procedures for the determination of parameters for the DSC models based on laboratory tests have been developed. Various models from the DSC have been validated at the specimen level with respect to laboratory test data. They have been validated at the practical boundary value problem level by comparing observed behavior in the field and/or simulated problems in the laboratory with predictions using computer (finite element) procedures in which DSC has been implemented; these have been presented in various publications by Desai and coworkers, and are listed in the References. Three typical examples of such validations at the practical problem level are included in this paper. It is believed that the DSC can provide unified Keywords constitutive modeling, disturbed state concept - DSC, geologic materials, interfaces/joints Dedication Professor Lujo Suklje was the pioneer of Slovenian Soil Mechanics, whose contributions were recognized internationally. He has made outstanding contributions in rheological properties of soils including creep and anisotropy, viscoplastic models, consolidation and foundation engineering. His book, Rheological Aspects of Soil Mechanics published in 1969, has left an indelible mark on the literature in engineering. I dedicate, with profound respects, this paper to the memory of Professor Suklje. 1 INTRODUCTION In the continuing pursuits for realistic models for geologic materials and interfaces or joints, a great number of simple, and simplified to complex constitutive models have been proposed. It is not the intention here to present a detailed review of such models, which includes elasticity, plasticity, creep, damage, microme-chanics, micro-cracking leading to fracture and failure, and liquefaction. Such reviews are available in many publications, including Desai (2001b). 1.1 SCOPE The modeling approach based on the disturbed state concept (DSC) is the main objective in this paper. The DSC provides a hierarchical and unified framework to include other models such as elasticity, conventional plasticity, continuous hardening plasticity, creep, and micro-cracking leading to softening as special cases. Also, the number of parameters required in various versions of the DSC model is smaller or equal to the number in other available models of similar capabilities. The DSC model has been implemented in two- and three-dimensional problems for dry and saturated geologic materials under static and cyclic (dynamic) loadings. Since the DSC model results in a special form of nonsymmetrical matrix equations, it is necessary to develop special schemes for the implementation of the DSC models in computer (finite element) procedures. The descriptions of the DSC model and the finite element procedure are divided in the following components: - Basic concept and hierarchical framework. - Parameter determination and calibration from laboratory and/or field tests. - Brief description of a special scheme for implementation of the DSC in the finite element procedure. - Validations with respect to tests used for calibration and independent tests, and applications of the finite element method for solutions of a number of practical problems in geotechnical engineering. 2 BASIC CONCEPT AND HIERARCHICAL FRAMEWORK Details of the DSC are given in various publications (e.g. Desai 1974, Desai 1987, Desai and Ma 1996, Desai and Toth 199, Katti and Desai 1995, Park and Desai 2000), which are referenced in Desai (2001b). A brief description of the framework is given below with attention to the calibration for DSC parameters, numerical implementation and validations. 2.1 EQUATIONS The DSC is based on the idea that the observed behavior of a material can be expressed through the behavior of reference components in the deforming material. The reference components are usually called relative intact (RI) and fully adjusted (FA). The behavior of the RI part, which may not be affected by microstructural discontinuities, can be defined by using a continuum model, e.g. elasticity, plasticity or viscoplasticity. The RI part is transformed continuously to the FA part, distributed at random locations, Fig. 1. The RI and FA parts are coupled and interact to yield the observed response, which is expressed in terms of the RI and FA Ri n ri ""FA FA -, D=0 ÊS fc- — FA / © Ê ) 0 ËL. '-•RI © D„ Df D D=1 D=0 (or D0) D=0 D>0 (a) Clusters of RI and FA parts i ^ relative intact a ^ observed c ^ fully adjusted (b) Symbolic representaion of DSC (c) Shematic of stress-strain response Figure 1. Representations of DSC. D^D^l O R c behavior, connected through the disturbance function, which provides for the coupling. Based on the above concept, the incremental constitutive equations for DSC are derived as da" =(1-D) C'de + DCcdec + dn(ac-a') (1a) ~ ~ ~ (1b) where a = stress vector, e = strain vector, a, i and c denoteobserved, RI and FA behavior, respectively, C and Cc = constitutive matrices for the RI and FA parts, respectively, D = disturbance and d denotes increment. The RI behavior can be expressed by using a suitable continuum model, e.g., linear elastic, nonlinear elastic, elastic-plastic or viscoplasticity. Here, the hierarchical single surface (HISS) plasticity (associative) model (80), derived from the general basis (Desai 1980), is used in which the yield function, F (Desai et al. 1986), is given by F = Jin-(-a (1-b sr; = 0 (2) where J2D = second invariant of the deviatoric stress tensor, J1 = J + 3R, J1 = first invariant of the stress tensor, R = parameter related to the cohesive strength, c, Fig. 2, y = parameter related to the ultimate yield surface, ft = parameter related to the shape of F in a1 - a2 - a3 stress space, Sr = stress ratio = ((27 / 2) (J3D / J2/2), J3D= third invariant of the deviatoric stress tensor, and a is the growth or hardening function. In Eq. (2), the stress invariants are nondimensionalized with respect to pa , the atmospheric pressure constant. In a simple form, the growth or hardening function is given by a=" (3a) X1 where a1 and n1 = hardening parameters and £ is the trajectory of or accumulated plastic strains: x= /(dee; ■ deep)1/2 (3b) in which eipj = plastic strain tensor, which is the sum of the deviatoric and volumetric plastic strains. The FA part can be defined by assuming that it has no strength like in classical damage model (Kachanov 1986), or it has hydrostatic strength like in classical plasticity, or it has strength corresponding to the critical state model. It can be also defined by using the behavior of other states such as zero suction for partially saturated soils. For soils, it is found appropriate to use the critical state equations to define the FA response (Roscoe, et al. 1958, Desai 2001b): ec = ec„-\ln (4a) (4b) where c denotes the critical state, e = void ratio, eco initial void ratio and m, and X are the critical state parameters. or Ultimate Envelope (a) (Jjp)« - J, (b) Octahedral plane; (0 < 0.756 for convexity) Figure 2. Plots of F in different stress spaces. 2.2 DISTURBANCE FUNCTION The disturbance function, D, can be defined on the basis of observed stresses, volumetric strains, pore water pressures or nondestructive properties such as S- and P-wave velocities. In terms of stress, D is expressed as D = (5a) S — s The disturbance, D, is also expressed as (Fig. 3) D = Du (1 — e-AxD) (5b) where Du = ultimate disturbance, A and Z are disturbance parameters, and £D is the trajectory of deviatoric plastic strains. Figure 4. Hierarchical versions in DSC. Figure 3. Disturbance for softening, stiffening behavior and critical disturbance. 2.3 HIERARCHICAL APPROACH Various continuum models can be derived from the DSC equations, Eq. (1). If D = 0, i.e., the material does not experience any microstructural changes leading to discontinuities, Eq. (1) reduce to that for standard continuum models: d S = C d ei (6) where C = constitutive matrix for continuum model, e.g., elasticity, plasticity, etc. Various models such as von Mises, Drucker-Prager, critical state, Matsuoka and Nakai (1974), Lade and Kim (1988), etc. can be obtained as special cases of the HISS model, Eq. (2). Fig. 4 shows various versions that can be derived from the DSC. 2.4 interface/joints Since both the solid material (soil) and the interfaces are coupled, it is appropriate to use the same mathematical framework for both. One of the major advantages of the DSC is that it provides the same framework for the solid material and the interfaces. It has been formulated for interface/joint behavior with suitable yield function for RI behavior and critical state models for FA behavior (Desai et al. 1984, Desai and Ma 1992, Desai et al. 1995, Desai 2001b, Desai et al. 2005). This paper essentially addresses the soil behavior. 2.5 ADVANTAGES OF DSC/HISS MODELS The DSC allows for discontinuities experienced by a deforming material, that can often result in degradation or softening; it can also allow stiffening and healing in deforming materials. It is essential to include the existence and effect of discontinuous material in the model; then only can it provide consistent and improved characterization compared to other models based on continuum plasticity proposed and used to account for such behavior (Mroz et al. 1978), Pestana and Whittle 1999, Elgamal et al 2002). Because the DSC allows for the coupling between the RI and FA responses, it can avoid such difficulties as spurious mesh dependence that occurs in classical damage models (Kachanov 1986). Also, if such a coupling is considered in the classical damage model by introducing additional enrichments (e.g. Bazant 1994, Muhlhaus 1995), the resulting models may be complex. In the fracture mechanics approach, usually it is required to introduce cracks in advance of the loading. In contrast, the DSC does not need a priori introduction of cracks, and initiation and growth of micro-cracking, fracture and failure can be traced at appropriate locations depending upon the geometry, loading and boundary conditions, on the basis of critical disturbance, Dc , obtained from test results (Desai 2001b,Park and Desai 2000, Pradhan and Desai 2006). The DSC allows identification of and initiation and growth of microstructural instability or liquefaction by using the critical disturbance (Desai et al. 1998, Desai 2000, Park and Desai 2000, Pradhan and Desai 2006). The DSC/HISS models also possess the following advantages: 1. The yielding in the HISS model is assumed to be dependent on total plastic strains or plastic work. Hence, in contrast to other models such as critical state and cap in which yielding is based on the volumetric strains, the HISS model includes the effect of plastic shear strains also. 2. The yield surface allows for different strengths along different stress paths. 3. Because of the continuous and special shape of the yield surface, it allows for dilatational strains before the peak stress, which can be common in many geomaterials, Fig. 2(a), and later Fig. 7(b). 4. The DSC allows for degradation and softening, and it can allow also for stiffening or healing. 5. Introduction of disturbance can account for (a part) the nonassociative behavior, i.e., deviation of plastic strain increment from normality. 6. The DSC allows intrinsically the coupling between RI and FA parts and the nonlocal effects. Hence, it is not necessary to add extra or special enrichments, e.g., microcrack interaction, gradient and Cosserat schemes, etc. (Muhlhaus 1995). 7. The DSC is general and can be used for a wide range of materials: geologic, concrete, asphalt, ceramics, metal, alloys and silicon (Desai 2001b), if appropriate test data is available. 8. It can also be used for repetitive loading which may involve a large number of loading cycles; a procedure is described below. 2.6 REPETITIVE LOADING: ACCELERATED PROCEDURE Computer analysis for 2- and 3-D idealizations for time dependent problems can be time consuming and expensive, especially when significantly greater number of cycles of loading need to be considered. Therefore, approximate and accelerated analysis procedures have been developed for a wide range of problems in civil (pavements) (Huang 1993), mechanical engineering and electronic packaging (Desai and Whitenack 2001). Here, the computer analysis is performed for only a selected initial cycles (say, 10, 20), and then the growth of plastic strains is estimated on the basis of empirical relation between plastic strains and number of cycles, obtained from laboratory test data. A general procedure with some new factors has been developed and applied for analysis of pavements and chip-substrate systems in electronic packaging (Desai 2007; Desai and Whitenack 2001). From experimental cyclic tests on engineering materials, the relation between plastic strains (in the case of DSC, the deviatoric plastic strain trajectory, £D), and the number of loading cycles, N, can be expressed as, Fig. 5: Xd ( N ) = x (Nr N Nr (7) where Nr = reference cycle, and b is a parameter, depicted in Fig. 5. The disturbance equation (5b) can be written as D = D 1 _ exp (-A {Xd (N)Z } (8) ln(UN)) In(iV) Figure 5. Plastic strain trajectory vs. number of cycles for approximate accelerated analysis. Substitution of (£D(N)) from Eq. (7) in Eq. (8) leads to N = N D,, D _ D (9) Now, Eq. (9) can be used to find the cycles to failure, Nf, for chosen critical value of disturbance = i (say, 0.50, 0.75, 0.80). The accelerated approximate procedure for repetitive load is based on the assumption that during the repeated load applications, there is no inertia due to "dynamic" effects in loading. The inertia and time dependence can be analyzed by using the 3-D and 2-D procedures; however, for million cycles, it can be highly time consuming. Applications of repeated load in the approximate procedure involve the following steps: 1. Perform full 2-D or 3-D FE analysis for cycles up to Nr , and evaluate the values of £D(Nr) in all elements (or at Gauss points). 2. Using Eq. (7) compute £D(N) at selected cycles in all elements. 3. Compute disturbance in all elements using Eq. (8). 4. Compute cycles to failure Nfby using Eq. (9), for the chosen value of Dc . 5. The above value of disturbance allows plot of contours of D in the finite element mesh, and based on the adopted value of Dc , it is possible to evaluate extent of fracture and cycles to failure, Nf Loading-Unloading-Reloading: Special procedures and integrated in the codes to allow for loading, unloading and reloading during the repetitive loads; details are given in (Shao and Desai 2000, Desai 2001b). PARAMETER DETERMINATION AND CALIBRATION FROM LABORATORY AND /OR FIELD TESTS 3.1 PARAMETERS Based on the RI (Eq. 2), FA (Eq. 4) and disturbance, D (Eq. 5), the parameters required are: Elastic E and v (or K and G) RI: Plasticity y, ft, n, a1, n1 FA: Critical state e0, m, \ Disturbance Du , A and Z where E = elastic modulus, v = Poisson's ratio, K = bulk modulus, and G = shear modulus. The above parameters are averaged from the parameters obtained from test data under different confining pressures, and other facors. If the variation is severe, then a parameter needs to be expressed as a function of the factor, which increases the number of parameters. The number of parameters can increase if additional factors such as temperature and rate dependence are considered (Desai 2001b). The following Table 1 gives some of the specialized versions with related number parameters. It may be noted that the number of parameters in the DSC model is less than or equal to those in any other available model of comparable capabilities. Table 1. Specialized Versions and Parameters from DSC/HISS Model. Specialization Number of parameters Elastic 2 Classical Elastic-plastic: - von Mises 3 - Mohr Coulomb 4 - Drucker-Prager 4 Continuous Hardening: - Critical State 5 HISS 50 Model (associative plasticity) 8* HISS S1 Model (nonassociative plasticity) 9 HISS S0+vp Viscoplastic 11 HISS 52 Anisotropic 12 HISS 50 Disturbance 12 * If the cohesive strength is not included, the number of parameters reduces by 1. Note: The number of parameters increases by about 2 if pore water pressure is considered and by about 2 if temperature is included. 3.2 PARAMETER DETERMINATION The laboratory tests used for the calibration include consolidation (hydrostatic), triaxial and multiaxial and/or shear tests. The field test may include nondestructive, e.g., S-wave measurements. Direct shear and simple shear (e.g., using the CYMDOF device, Desai and Rigby, 1997) tests are used for DSC parameters for interfaces/joints (Desai 2001b). The values of elastic moduli are found from the slopes of unloading responses. For instance, E is found from unloading slopes of s — s3 vs. e1 plots from triaxial tests. Similarly, v can be found from plot of £v vs. e1 (ev = volumetric strain), G can be found from shear stress vs. shear strain plot, and K can be found from mean pressure (p) vs. ev plot. Determinations of various elastic parameters are depicted in Fig. 6. a) c) Ej (Ej-EJ) b) d) Figure 6. Elastic constants from test data; (a) stress-strain behavior, (b) volumetric-strain behavior, (c) shear stress-strain behavior, and (d) pressure -volumetric stain behavior. The parameters y and ft are related to the ultimate response, which is often adopted from the asymptotic value (about 5 to 10% higher than the peak stress) for a given stress-strain curve, Fig. 7(a). They are found from the specialized expression for F (where a = 0), Eq. (2), at the ultimate stress conditions, based on tests under different confining pressures: F = J 2D -7 Ji2 (1-b Sr r5° = 0 (10) By substituting different values of ultimate stresses, Fig. 2, y and ft are found by using procedures such as the least-square method. The phase change parameter is found from the point (b) in the stress-strain behavior, Fig. 7(b), where the response transits from compressive to dilative. One of the expressions for n is " C3 (V£d) 2 1- J2D J_ Ji Fs7 (11) —c (FA) i (ep) 0 «) (a) ^-ec (d) • Loose e \ (v) Dense \ c) rw ec In (j;/3 Pa) (e) Figure 7. Behavior of loose and dense geo-materials: (a) stress-strain behavior, (b) volumetric-strain behavior, (c) void ratio-pressure behavior, (d) Ji - ^jJ2D - critical state behavior, and (e) void ratio-critical state behavior. n= The hardening parameters are determined based on the following expression from Eq. (3). £na + T}linX= in a1 (12) Here, the stress-strain curve, Fig. 8(a) is divided into increments and £ is found as the accumulated plastic strain for a given increment. Then a for a given £ is found from Eq. (1), i.e., F = 0. The plot of (Fig. 8b) Ina vs. ln£ yield the (average) values of a1 and qv The critical state parameters, Eq. (4), are obtained from the plots of „JJ2D — J1 for ; and A, Figs. 7(d) and 7(e). and M/i / Pa for a) Cna b) ■■■tx-. Figure 8. Determination of hardening or growth parameters, a1 and nl. The disturbance parameter, Du, is found from Eq. (5b) by substituting the value of the residual (FA) stress, Fig. 7(a). Then the values of A and Z are obtained by plotting in(-in(Du -D)/Du) vs. ln£D , Fig. 9, based on the following expression from Eq. (5): Zin(XD - InA = In —In D,, — D D (13) Figure 9. Determination of disturbance parameters, A and Z. The values of D can be found from static or cyclic response, Fig. 1(c) or Fig. 10. 3-3 PARAMETERS FOR TYPICAL MATERIALS AND INTERFACES/JOINTS The DSC models have been used to model a wide range of materials, e.g., clays, sands (dry and saturated), glacial tills, partially saturated soils, rocks, concrete, asphalt concrete, metals, alloys (solder), silicon., and various interfaces and joints. Application for interfaces includes clay-steel, sand-steel, sand-concrete and chip-substrate. Laboratory and/or field tests have been used to calibrate the DSC parameters by following the foregoing proce- Figure 10. Cyclic tests for disturbance and liquefaction. c dures. Details of the parameters have been presented in a number of publications, and also given in Desai (2001b). Parameters for the relevant geologic materials and interfaces are given in the following Examples 1 to 3. 3.4 HIERARCHICAL VERSIONS Various plasticity models such as nonassociative (S1) and anisotropic hardening (S2) can be obtained by modifying or adding "corrections" in the HISS model for the factors that influence the specific behavior, to the basic associative (S0) model (Desai et al. 1986, Somasundaram and Desai 1988, Wathugala and Desai 1993). Since the introduction of disturbance can account for (some) nonassociative response, the DSC model with the basis associative (S0) version is found to account for behavior of most geologic materials. 3.5 OTHER CAPABILITIES IN DSC 3.5.1 Creep Models A schematic of the creep behavior is shown in Fig. 11. Various creep models, e.g., Maxwell, viscoelastic (ve), elastoviscoplastic (evp), and viscoelastic- viscoplastic (vevp), can be obtained as specialization of the Multi-component DSC (MDSC) (Desai, 2001b). For example, the evp model can be defined based on the Perzyna's theory (1966), in which viscoplastic strain increment, e"p, is expressed as e = r^)°F (14) oa where r = fluidity parameter, and < > has the meaning of a switch-on-switch-off operator as F if — > 0 —o F if — < 0 F (15) Fo = reference value of F (e.g., yield stress or pa) and 0 = flow function. One of the expressions of 0 is f = (16) In above equations, r and N are the viscous parameters. The values of r and N are found from a creep test, Fig. 11; details are given in (Desai et al. 1995, Desai 2001b). Various creep models such as viscoelastic (ve), viscoplastic (vp), viscoelastic -viscoplastic (vevp) can be obtained from the MDSC; they are depicted in Fig. 12. The Overlay or mechanical sub layer models (Duwez 1935, Zienkiewicz el al. 1972, Pande et al. 1977) can be considered to be special case of the MDSC approach. 3.5.2 Liquefaction Liquefaction in saturated soils (Fig. 10) subjected to dynamic (or static) loading can be identified in the DSC by locating the critical value of the disturbance, Dc, Eq. (5), Fig. 3, that represents the microstructural instability. There is no additional parameter needed for liquefaction. Details of the development and use of the DSC for liquefaction are given in various publications (Desai et al. 1998b, Desai 2000, Park and Desai 2000, Pradhan and Desai 2006). Figure 11. Schematic of creep behavior. Unit 1 Unit 2 ...................................... Unit k Figure 12. (a) MDSC creep models, (b) viscoelastic (ve), (c) viscoplastic (evp), and (d) viscoelastic-viscoplastic (vevp) models (Desai 2001b). 3.5.3 partially saturated soils The DSC model can be applied to model partially saturated soil by including factors such as suction and saturation (Geiser et al. 1997, Desai 2001b). 3.5.4 stiffening and Healing The disturbance, D, can simulate softening or degradation. It can also define stiffening or healing (Fig. 3) in the microstructure; for example, chemical and thermal (causing particle fusing) effects can increase in the bonding in the microstructure (Desai et al 1998a, Desai, 2001b). 3.5.5 Thermal effects The thermal effects can be introduced in the DSC model. Very often, the parameters are expressed as a function of temperature, T, as (Desai 2001b) P = Pr (17) where p = any parameter, and pr is the value of the parameters at a reference temperature, Tr (e.g., the room temperature), and c = parameter. 3.5.6 Rate effects Behavior of many materials is affected by the rate of loading. For a plasticity model, the rate effect can be introduced by modifying the static yield function determined from tests conducted under very slow or static rate ( Baladi and Rohani 1984, Sane et al. 2009): Fd = F - fr) / Fo (18) where Fd represents the yield function defined for 'dynamic ' rate, Fs represents the yield function at static or very small rate, Fo is the used to nondimensionalize the terms, and FR is the rate dependent function. Details for the determination FR and Fs are given in (Sane et al. 2009). 3.6 COMMENTS ON VRRIOUS ISSUES 3.6.1 physical Meanings and curve fitting A constitutive model for defining behavior of a material requires a certain level of curve fitting for calibration of parameters based on the test data. It is useful and desirable that such curve fitting is reduced as much as possible, particularly if the behavior is influenced by many factors. One of the ways is to develop parameters that have physical meanings, e.g., the parameters are associated with physical states during the response, such as ultimate condition, and change in volumetric strain. Most of the parameters in the DSC have such physical meanings; hence, the need for curve fitting is reduced and the procedures for the determination of parameters are simplified. 3.6.2 Simple and Simplified Models It is often said that a constitutive model should be simple, particularly for practical applications. A functional representation (e.g., parabola, hyperbola, spline, etc.) of a given (single) stress-strain response as a nonlinear elastic model, is considered to be a "simple" model. Such a function like the hyperbola extended to include the effect of the confining pressure may still be called "simple". However, since the soil behavior is affected by other practical factors such as irreversible (plastic) deformations, stress path, volumetric strains, creep and type of loading, the resulting model may not be simple! It is apparent that models by using only mathematical functions (e.g. hyperbola, spline, etc) to represent stress-strain curves are not able to account for such practical factors. Then, it is better to develop and use "simplified" models. A simplified model should be derived on the basis of basic principles of mechanics, and is intended to characterize the behavior affected by significant factors important for practical problems. It is often obtained by deleting less significant factors from a general model developed by considering the physical nature of the problem, loading and other relevant conditions. In other words, a "simplified" model could include the effect of the significant factors, in which the number of parameters could increase with the inclusion of each additional factor. An important issue is how many parameters the model involves, and what type of tests are required to determine the parameters. The other important issue is how many significant factors are vital for practical problems and must be included in the model? It is natural that the number of parameters would increase with the number of factors. Hence, for a practical application, one must adopt a model that contains the influence of the significant factors for a particular application. Then the notion of adopting just a simple model that is not capable of accounting for practically important factors would have no basis! Indeed, for practical applications, it is desirable to develop as simplified model as possible with smaller number of parameters. The DSC/HISS model has been developed such that its simplified versions can be adopted depending upon the need of the practical application. Moreover, the DSC/HISS model is hierarchical in nature and allows selection of a model for required factors by introducing additional factors and related parameters in the same basic model. The hierarchical approach in the DSCC/HISS model can reduce the number of parameters, and simplify the procedure for their determination based on standard (e.g., triaxial and shear) tests. 3.6.3 dsc and Damage Models The classical damage approach (Kachanov 1986) is based on the physical cracks (damage) in a deforming material, and the resulting effect on the observed behavior. The DSC is different from the damage approach in that it considers a deforming material as a mixture of two or more reference components (e.g., RI and FA states); then the observed behavior is expressed in terms of the behavior of the components, which interact to yield the observed behavior. The disturbance represents the deviation of the observed behavior from the reference state(s). Damage model can be considered as a special case of DSC if the interaction is removed, e.g., the FA state is assumed to possess no strength. However, such a model, without interaction between RI and FA states, may not be realistic and possesses certain (computational) difficulties (Desai 2001b). Moreover, the DSC approach is general and can include both softening and stiffening effects. 4 FINITE ELEMENT FORMULATION INCLUDING SPECIAL SCHEMES The procedure for implementation of the DSC model in computer (finite element) methods is given below. The finite element equations for two- or three-dimensional problems (Fig. 13) can be derived as (Desai, 2001b): fBT a" dV = Q (19a) V n + 1 ~ n + 1 where B = strain-displacement transformation matrix, 5 = observed stress vector at n + 1 step and the load ~ n+1 vector: Q =f Nt X dV +f Nt T d S/ (19b) ~ n+1 J ~ ~ J ~ ~ where N = matrix of interpolation functions, X = body force vector, T = surface traction vector, V denotes volume and S~ denotes a part of the surface. Figure 13. Finite Element Discretization: 2-D and 3-D elements. Now with S+1 =S + dS+ (20) where n denotes increment, the incremental equations f BT do" dV = Q -f Bt O dV J ~ ~ n+1 ~ n+i J ~ n (21) = Q - Qb = d Qb where Q represents the internal balanced load at step n, Fig. 14, and d Qb denotes the out-of-balance or residual load vector.1 Now, by substituting Eq. (20) in Eq. (19a), the incremental equations with the DSC model are: J BTCDSC B dV • dqn+i = dQb+i (22) where CDSC = (l - Dn ) C + Dn ( 1 + a) C^ + RT S; ~ n ~ n ~ S =sc -a' , dDn = Rt de (Desai, 2001b), and ~ n ~ n ~ n ~ ~ n dS = ( 1 + a) Cc de' and a is the factor for relation between RI and FA strains. Now, we can write Eq. (22) as kDSC dqn+1 = dQ6 (23) ~ n ~ n+1 where kDSC = f BT CDSC ■ B ■ d V . ~ rj ~ ~ n ~ y A number of procedures are possible for the solution of Eq. (23) (Desai, 2001b). Since kDSC is nonsymmetrical, it may cause computational problems. Simplified and approximate procedure used commonly is briefly described below. This procedure is based on the solution of the RI behavior in which the stiffness matrix is symmetric: k dq' = dQn+1 = Qn+1 - Q b(i) (24) are n Figure 14. Incremental iterative analysis. Figure 15. Iterative procedure in solution based on RI equations. Fig. 15 shows the schematic of the procedure where the solution, Eq. (24) pertains to the continuously hardening (RI) behavior. Once the incremental displacements, d q' , for the RI behavior are determined, the strains d s' ) and stresses (d S ) are determined. Then the ~ n+lj \ ~ n + 1/ observed stresses are found from the following equation: ds+) =(1-Dn+1) S + D'n+1 S+1 (25) ~ n+1 v ' ~ n+1 where j denotes iterations depicted in Fig. (15) (Desai, 2001b). Then the observed (applied) load is found as: d Qa = f Br sa(m) dV (26) ~ n+1 ~ ~ "+1 where m denotes the last (converged) iteration. 5 VALIDATIONS The DSC equations (1) are integrated to predict the observed laboratory behavior for a wide range of materials (Desai, 2001b). For most of the materials, the DSC model provides very good correlations with the test data for the specimen level tests. Also, the DSC model has provided very good correlations for independent test data, which were not used for calibration of the parameters. Details of validations at the practical boundary problem level are presented in a number of publications, e.g. Desai (2001b); three typical examples of applications are given below. 5.1 COMPUTER CODES A computer code have been developed for (a) the determination of DSC parameters based on laboratory data. Here, the user can input points on various stress-strain-volumetric responses, and the code yields the parameters. This code also includes the integration of Eq. (1) by using the determined parameters so as to back predict stress-strain-volumetric responses toward validations at the specimen level. Computer codes have been developed based on the above finite element method for static and dynamic analysis including liquefaction of practical boundary value problems, e.g. ( i ) DSC-SST2D (Desai 1999), (ii) DSC-DYN2D (Desai 2001c), and (iii) DSC-DYN3D (Desai 2001a). These and other codes have been used to solve a wide range of problems in civil, mechanical and electronic engineering. The predictions compare very well with laboratory and/or field tests for practical problems. Three typical example problems for application of these codes are given below. 5.2 EXAMPLE i: DYNAMIC BEHAVIOR OF PILE TESTED IN CENTRIFUGE Fig. 16 shows the pile tested in the centrifuge test facility (Wilson, et al., 2000). The pile-Nevada sand system was subjected to base motion shown in Fig. 17. Various items such as pore water pressures, accelerations, displacements and bending moments were recorded with time by using instruments installed as marked in Fig. 16. Finite element analysis was performed using a procedure based on generalized Biot's theory in which the DSC model was implemented (Desai, 2001c; Pradhan and Desai, 2006). The DSC model parameters for the Nevada sand were determined from laboratory triaxial tests reported by The Earth Technology Corporation for VELACS (Arulmoli, et al. 1992). The test data included three undrained monotonic triaxial tests and three undrained cyclic triaxial tests with continuing pressures, a3 = 40, 80 and 160 kPa at relative density, Dr = 60%. Details of the evaluation of parameters for the Nevada sand are given in Pradhan and Desai (2002); the values of the parameters are presented in Table 2. The interface test data for the aluminum pile-Nevada sand were not available. Hence, a neural network approach was used for estimation of the parameter (Pradhan and Desai, 2002). Here, two sets of data for training of the neural network were used. Set 1 consisted of Ottawa sand parameters using triaxial tests (Park and Desai, 2000), and the Ottawa sand-steel interface using the cyclic multi-degree-of-freedom (CYMDOF) device (Alanazy, 1996 and Desai and Rigby, 1997). The second set consisted of parameters for a marine clay from triaxial and multiaxial testing (Katti and Desai, 1995) and marine clay-steel interface using the CYMDOF device (Shao and Desai, 2000). It was assumed that the steel and aluminum behave approximately similar in interfaces with the soil. Also, the Nevada sand gradation curve falls between those of the Ottawa sand and the marine clay. Finally, the parameters for Nevada sand-aluminum interface estimated from the neural network were found and are given in Table 2. The sand and interface (aluminum-sand) were modeled using the DSC with the parameters given in Table 2. The aluminum pile itself was assumed to be linear elastic with elastic modulus = 70.0 GPa and Poisson's ratio = 0.33. The single pile, Fig. 16, was made of aluminum pipe with diameter of 0.67 m with wall thickness of 72 mm; the total depth of sand was 20.7 m and the embedment length of the pile was 16.8 m. The pile-system was idealized as plane strain (Pradhan and Desai, 2006; Anandara-jah, 1992). The finite element mesh for the total domain in the centrifuge test for the pile is shown in Fig. 18(a); details of mesh around the pile are given in Fig. 18(b). ■ pore pressure 11 bending moment r=i- displacement » acceleration Figure 16. Schematic of pile in centrifuge test and instrumentation (from Wilson et al. 2000). 0.25 0.20 -0.20 - -0.25 4—-1-—1-1-—---------------- 0 5 10 15 20 25 Time (sec) Figure 17. Base motion acceleration (from Wilson et al. 2000). Table 2. Parameters for Nevada Sand and Sand-Aluminum Interface. Material Parameter Symbol Nevada Sand Nevada Sand - Aluminium Interface RI: Elastic E V 41.0 MPa 0.316 14.6 MPa 0.384 Y 0.0675 0.246 ß 0.00 0.00 RI: Plasticity 3R n 0.00 4.10 0.00 3.360 ai 0.1245 0.620 ni 0.0725 0.570 m 0.22 0.304 FA: Critical State X 0.02 0.0278 eO 0.712 0.791 Du 0.99 0.990 Disturbance Z 0.411 1.195 A 5.02 0.595 The FE mesh contained 302, 4-noded elements with 342 nodes. Interface elements were provided for the vicinity of pile by assuming small thickness of the sand. Boundary BC (Fig. 18a) was restrained in the y-direction, and AB and CD were free to move in both x- and y-directions. The base of mesh was subjected to the acceleration-time history shown in Fig. 17 (Wilson, et al., 2000). The initial (in situ) and pore-water pressures at the centers of all elements were introduced by using: S = gh ; S = Ko s , Ko = n/(1- n , and P =gwh where and Sh = effective vertical and horizontal stresses at depth, h; ys = submerged unit weight of soil; yw = unit weight of water; p = (initial) pore water pressure, and Ko = co-efficient of earth pressure at rest. (27) 5.2 RESULTS To identify the influence of interfaces between structure (pile) and soil, analyses were performed (a) without interfaces and (b) with interfaces. Only typical results are presented herein. Fig. 19 shows comparisons between measured pore water pressure in element 139 in the vicinity of the pile, and predictions with and without interface. The observed pore water pressure approaches the initial effective vertical pressure (S), showing potential liquefaction, after about 9 seconds. The analysis without interface shows that the computed pore water pressure remains far from Sv, i.e., it does not show liquefaction. However, the analysis with the interface shows that after about 9 seconds the pore water pressure is almost equal to s/, indicating liquefaction. Thus, provision of interface including related relative motions is considered to be essential for realistic prediction pore water pressures, liquefaction and behavior of the pile-soil system. Interface elements - 133 to 143 and 162 to 172 Pile elements - 146 to 159 12 - 22,667 m a) 3.7 m 20.7 m 159 130 143 156 172 1S5 128 141 170 183 xx, 139 169 181 — F.mb edmenl E cptb = 16.8m [24 ¡37 166 179 131 146 162 IIS 131 144 160 173 1.07 m b) 3.7 ra Figure 18. Finite element mesh for pile test. 0 9 1» V 0 9 1! 27 Tine, sec Time, see £ t; S. 0 9 IB 27 Time, sec Figure 19. Comparisons of pore water pressures for element 139 in vicinity of pile. Table 3 shows time to liquefaction in the conventional procedure when the pore water pressure, Uw, reaches the initial effective pressure s time for reaching the critical disturbance, Dc = 0.86 (determined from the plot of disturbance vs. time), is considered to designate microstructural instability. Although the values of the two times are not too different, the times to reach Dc are lower than those to reach S. This can imply that the microstructural instability may occur earlier than the time to liquefaction obtained from the conventional procedure. The conventional procedure suggests that the soil does not possess any strength when liquefaction occurs. In fact, it is observed that the soil can retain a small strength when the liquefaction is reached; this is also shown in observations during the Port Island, Kobe, Japan earthquake (Desai, 2000). Hence, it could be stated that the DSC based on the critical disturbance, Dc , could provide more accurate time to liquefaction (microstructural instability). Table 3. Times to Liquefaction from Conventional and Disturbance Methods. Element Conventional Uw =S(s) Disturbance D = Dc (s) 143 1.74 1.23 130 1.83 1.29 104 2.55 1.92 78 3.36 2.67 52 3.81 2.94 26 9.03 8.22 5.3 EXAMPLE 2: REINFORCED EARTH RETAINING WALL The DSC parameters for backfill (soil) and interface for a Tensar reinforced wall, Fig. 20 (USDOT 1989), were defined based on comprehensive laboratory tests, triaxial for soil, and interface shear for interfaces between soil and (Tensar) reinforcement (Desai and El Hoseiny, 2003). The DSC parameters for the backfill soil and interfaces are shown in Table 4. The finite element procedure with the DSC model was used to analyze the earth retaining wall reinforced by using Tensar geogrid. Coarse and fine meshes were used for the analyses; it was found that the fine mesh analysis yields improved predictions compared to those from the coarse mesh. Fig. 21 shows a part of the fine mesh (with 11 layers); it contained 1188 nodes, and 1370 elements including 480 interface, 35 wall facings and 250 bar (reinforcement) elements. The analysis included simulation of construction sequences, in which the backfill was constructed in 11 layers, as was done in the field. The reinforcement was placed on a layer after it was compacted. The completion of the sequences of construction is referred to as "end of construction." Then the surcharge load due to traffic of 20 kPa was added uniformly on the top of the mesh, which was referred to as "after opening to traffic." The in situ stress was introduced by using the co-efficient K„ = 0.40. Table 4. Parameters to Analysis of Reinforced Wall. Material Parameter Symbol Soil Interface RI: Elastic E or kn v or ks fi(^s)ai 0.30 f2(*n)b fs^y Y 0.12 2.30 P 0.45 0.00 RI: Plasticity n ai 2.56 3.0E-5 2.80 0.03 ni 0.98 1.00 kd 0.20 0.40 Du 0.93 - Disturbance Z A 0.37 1.60 - Angle of friciton and adhesion f/S/c« f = 40° S = 34° ca = 66 kPa Unit weight Y 18.00 kN/m3 - Co-efficient of earth pressure at rest K0 0.400 - « E = 62 x 103 of8 b kn (normal stiffness) = 18 x 103 S'29 c ks (shear stiffness) = 30 x 103 s'28 d k = nonassociative parameter Figure 20. Tensar Reinforced wall (USDOT 1989). 5.4 RESULTS Only typical results are presented herein. Figure 22(a) shows comparison between predictions and field measurements for the vertical stress at elevation, 1.53 m, at the end of construction; those after opening to traffic are shown in Fig. 22(b). This figure shows also the overburden stress, and the trapezoidal vertical stress used in design calculations. The FE predictions agree well with the measured values. Fig. 23 shows comparison between predictions and measurements for the horizontal stress carried by the geogrid near the wall face, at different elevations. The measurements are obtained from the strain gages installed on the geogrid. The predictions for the geogrid stresses compare well with the measurements. s _ * r~.......—^ ~ ... - 4 Field -•-Pralkdicn -Ov«ftaKtf*fl ---- TrtpookU --1-■- O 0,5 1 1.6 2 25 3 3.5 Distance from Wall Face (m) (a) End of construction Figure 22. Comparisons between field measurements and predictions for vertical stress at elevation 1.53 m. r y \ ]—Geogrid -FEM p N. i hp " 0 5 10 15 2C Horizontal Stress (kPa) Figure 23. Comparison of field measurements and prediction for horizontal stress carried by geogrid near wall face. 5.5 EXAMPLE 3. MODELING FOR GLACIAL TILL FOR MOTION OF GLACIERS Constitutive modeling of glacial till and till-ice interfaces play a significant role in predicting the movements of ice sheets or glaciers. The glacial till occurs between the ice and the base rock, and its deformations can contribute significantly to the motion of glaciers. A comprehensive laboratory triaxial testing for shear and creep behavior of two glacial tills, Tiskilwa and Sky Pilot, were performed and used to define the DSC model (Desai et al., 2008, Sane, et al., 2008). The triaxial shear tests were performed under confining pressures, S = 20, 50, 100, 200, 300, 500 kPa. Tests were also performed for S = 150 and 400 kPa for independent validations of the model. Creep tests with constant (vertical) stress of 20, 40, 60 and 80% of the peak stress observed in the shear tests under S3 = 100, 200 and 400 kPa; tests under s/ = 300 kPa were performed for independent validation. These tests were used to find the parameters for the DSC model; Table 5 gives the parameters for two tills. Note: Creep parameters for elastic-viscoplastic (evp) and viscoelastic-viscoplastic (vevp) models are given in Sane, et al. (2008). The conventional plasticity models such as the MohrCoulomb have been used for behavior of glacial and till, and for identifications of the motion of glaciers. However, based on laboratory tests and finite element computations, it was found that the Mohr-Coulomb criterion that assumes failure and motion at maximum stress reached at very small (elastic) strains may not provide realistic behavior of till and initiation of motion. Table 5. Parameters for Glacial Tills. Material Parameter Symbol Tiskilwa Till Sky Plot Till RI: Elastic E 57.75 MPa 50.50 MPa 0.45 0.45 V Y 0.0175 0.0092 ß 0.36 0.52 RI: Plasticity n 3.60 6.85 a¡ 2.0 E-4 1.35 E-7 n¡ 0.27 0.14 m 0.10 0.09 FA: Critical State X 0.21 0.16 < 0.48 0.52 Du 1.00 1.00 Disturbance Z 2.68 1.00 A 20.20 5.50 Figure 24. Photographs of triaxial test samples (a) before the test, and (b) after test for Sky Pilot till, and (c) and (d) after test,for Tiskilwa till. It was observed that the glacial tills experience the disturbance at distributed locations in the test specimen, and fails when the critical disturbance, Dc , occurs in critical fraction of volume of the specimen. Fig. 24 shows the deformed configurations of till specimens which do not show specific failure planes, as is assumed in conventionally. In fact, it is considered that the "failure" in terms of the initiation of motion of the ice sheet may occur long after the peak stress, and at approximately the onset of the residual condition at Dc. Figure 25. (a) Rate of change of displacement with applied stress predicted by finite element analysis, and (b) Rate of change of £ with applied strain for triaxial test Tiskilwa till under confining stress = 100 kPa. Note: After the induced peak shear stress = 60 kPa, stress reduces to about 10 kPa in the post peak region, which is indicated by the dotted line in (a) above. Figure 26. Simulated section of ice sheet on till and finite element mesh with loading. (Not to scale.) Fig. 25(b) shows the plot of plastic strain trajectory (£) vs. applied axial strain for the behavior of the triaxial specimen under typical confining stress, S = 100 kPa. The instability causing failure or initiation of motion appears to occur when Dc = 0.85 is reached, which indicates the transition from linear to nonlinear behavior. The above is also seen, Fig. 25(a), when the computed displacement vs. applied shear stress reaches Dc = 0.85 at the transition from linear to nonlinear behavior occurs. The results in Fig. 25(a) were obtained from the finite element analysis of a simulated problem involving ice sheets, till layer and bed rock, Fig. 26. The details of the finite element analysis are given in Sane, et al., 2008 and Desai, et al., 2008. Thus, the DSC provides a realistic model for the behavior of glacial till and for the prediction of movements of ice sheets or glaciers, a subject which could have significant influence on global climate change. 6 CONCLUSIONS Accurate modeling of the behavior of geologic materials and interfaces/joints are essential for realistic prediction of the behavior of geotechnical problems, e.g., by using computer (finite element) procedures. The disturbed state concept (DSC) provides unified and powerful constitutive models for solution of geotechnical and other engineering problems. This paper presents a description of the DSC including theoretical details, determination of parameters, and validations of the model at specimen level, and practical boundary value problem level. Typical examples of application of finite element procedures with the DSC model are presented and include comparisons between predictions with observations from field and laboratory simulated problems. The DSC model is considered to be unique and provides for successful modeling of a wide range of engineering materials, and interfaces and joints. ACKNOWLEDGMENT The results presented herein have been partly supported by a number of grants from the National Science Foundation, Washington, DC. Participation and contributions of a number of students and colleagues are acknowledged; some are included in the references. REFERENCES Alanazy, A.S. (1996). Testing and modeling of sand-steel interfaces under static and cyclic loading. Ph.D. Dissertation, Dept. of Civil Eng. and Eng. Mechanics, The University of Arizona, Tucson, AZ. Anandarajah, A. (1992). Fully coupled analysis of a single pile foundation in liquefiable sands. Report, Dept. of Civil Eng., The Johns Hopkins University, Baltimore, Md. Arulmoli, K., Muraleetharan, K.K., Hossain, M.M., and Fruth, L.S. (1992). VELACS: Verification of liquefaction analysis by centrifuge studies. Soil Data Report, The Earth Technology Corporation, Irvine, CA. Baladi, G.Y., and Rohani, B. (1984). Development of an Elastic Viscoplastic constitutive relationship for earth materials. Chapter 2 in "Mechanics of Materials", Edited by Desai, C. S., and Gallagher, R.H., John Wiley and Sons Ltd, United Kingdom, 23-43. Bazant, Z.P. (1994). Nonlocal damage theory based on micromechanics of crack interactions. J. Eng. Mech., ASCE, 120, 593-617. Desai, C. S. (1974). A consistent finite element technique for work-softening behavior. Proc., Int. Conf on Comp. Meth. in Nonlinear Mech., J.T. Oden (editor), Univ. of Texas, Austin, TX, USA. Desai, C.S. (1980). A general basis for yield, failure and potential functions in plasticity. Int. J. Numer. Analy. Methods Geomech., 4, 361-375. Desai, C. S. (1987). Further on unified hierarchical models based on alternative correction or disturbance approach. Report, Dept. of Civil Eng. & Eng. Mech., The Univ. of Arizona, Tucson, AZ, USA. Desai, C.S. (1999). DSC-SST2D-computer code for static, dynamic, creep and thermal analysis: Solid, structure and soil-structure problems. Reports and Manuals, Tucson, AZ, USA. Desai, C.S. (2000). Evaluation of liquefaction using disturbed state and energy approaches. J. of Geotech. and Geoenv. Eng., ASCE, 126, 7, 618-631. Desai, C.S. (2001a). DSC-SST3D-computer code for static and coupled dynamic analysis: Solid, (porous) structure and soil-structure problems. Report and Manuals, Tucson, AZ, USA. Desai, C.S. (2001b). Mechanics of Materials and Interfaces: The Disturbed State Concept, CRC Press, Boca Raton, FL, USA.. Desai, C.S. (2001c). DSC-DYN2D-dynamic and static analysis, dry and coupled porous saturated materials. Report and Manuals, Tucson, AZ, USA. Desai, C.S. (2007). Unified DSC constitutive model for pavement materials with numerical implementation. Int. J. of Geomech., 7, 2, 83-101. Desai, C.S., Dishongh, T., and Deneke, P. (1998a). Disturbed state constitutive model for thermome-chanical behavior of dislocated silicon with impurities. J. of Applied Physics, 84, 11, 5977-5984. Desai, C.S., and El-Hoseiny, K.E. (2005). Prediction of field behavior of reinforced soil wall using advanced constitutive model. J. of Geotech. and Geoenviron. Eng., ASCE, 131, 6, 729-739. Desai, C.S., and Ma. Y. (1992). Modelling of joints and interfaces using the disturbed state concept. Int. J. Num. and Analyt. Methods in Geomech., 16, 9, 623-653. Desai, C.S., Park, I.J., and Shao, C. (1998b). Fundamental yet simplified model for liquefaction instability. Int. J. Num. Analyt. Meth. in Geomech., 22, 7, 721-748. Desai, C.S., Pradhan, S.K., and Cohen, David (2005). Cyclic testing and constitutive modeling of saturated sand-concrete interfaces using the Disturbed State Concept. Int. J. of Geomechanics, 5, 4, 286-294. Desai, C.S., and Rigby, D.B. (1997). Cyclic interface and joint shear device including pore pressure effects. J. of Geotech. and Geoenviron. Eng., 123, 6, 568-579. Desai, C.S., Sane, S.M., and Jenson, J.W. (2008). Constitutive modeling and testing of glacial tills for prediction of motion of glaciers. Proc., Int. conf. on Conf. Methods and Advances in Geomech., Goa, India, 4494-4508. Desai, C.S., Samtani, N.C., and Vulliet, L. (1995). Constitutive modeling and analysis of creeping slopes. J. Geotech. Eng., ASCE, 121, 43-56. Desai, C.S., Somasundaram, S., and Frantziskonis, G. (1986). A hierarchical approach for constitutive modeling of geologic materials. Int. J. Numer. Anal. Methods in Geomech., 10, 3, 225-257. Desai, C.S., and Toth, J. (1996). Disturbed state constitutive modeling based on stress-strain and nondestructive behavior. Int. J. Solids & Struct., 33, 11, 1619-1650. Desai, C.S., and Whitenack, R. (2001). Review of models and the disturbed state concept for thermo -mechanical analysis in electronic packaging. J. Electron. Packaging, 123, 1-15. Desai, C.S., Zaman, M.M., Lightner, J.G., and Siriwar-dane, H.J.( 1984). Thin-layer element for interfaces and joints. Int. J. Num. and Analyt. Methods in Geomech., 8,1, 19-43. Duwez, P. (1935). On plasticity of crystals. Physical Reviews, 47, 6, 494-501. Elgamal, A., Yang, Z., and Parra, E. (2002). Computational modeling of cyclic mobilityand post-liquefaction site response. Soil Dynamics and Earthquake Eng., 22, 259-271. Geiser, F., Laloui, L., Vulliet, L., and Desai, C.S. (1997). Disturbed state concept for partially saturated soils. Proc., 6th Intl. Symp. on Num. Models in Geome-chanics, Montreal, Pietrusczka, S. and Pande, G.N. (Eds.), Canada, Balkema, Netherlands. Huang, Y.H. (1993). Pavement Analysis and Design, Prentice-Hall, Englewood Cliffs, NJ. Kachanov, L.M. (1986). Introduction to Continuum Damage Mechanics, Martinus Nijhoff Publishers, Dordrecht, The Netherlands. Katti, D.R., and Desai, C.S. (1995). Modeling and testing of cohesive soil using the disturbed state concept. J. of Eng. Mech., ASCE, 121, 5, 648-658. Lade, P.V., and Kim, M.K. (1988). Single hardening constitutive model for frictional materials. III. Comparisons with experimental data. Computers and Geotechnics, 6, 31-47. Matsuoka, H., and Nakai, T. (1974). Stress-deformation and strength characteristics of soil under three different principal stresses. Proc., Jap. Soc. Civil Engrs., 232, 59-70. Mroz, Z., Norris, V.A., and Zienkiewicz, O.C. (1978). An anisotropic hardening model for soils and its application to cyclic loading. Int. J. Num. & Analyt. Methods in Geomech., 2, 208-221. Muhlhaus, H.B. (Ed.) (1995). Continuum Models for Materials with Microstructure, John Wiley, UK. Pande, G.N., Owen , D.R.J., and Zienkiewicz, O.C. Overlay models in time dependent nonlinear material analysis," Computers and Structures, 7, 435-443. Park, I.J., and Desai, C.S. (2000). Cyclic behavior and liquefaction of sand using disturbed state concept. J. Geotech. Geoenv. Eng., ASCE, 126, 9. Perzyna, P. (1966). Fundamental problems in viscoplas-ticity. Adv. in Appl. Mech., 9, 243-277. Pestana, J.M., and Whittle, A.J. (1999). Formulation of a unified constitutive model for clays and sands. Int. J. Num. & Analyt. Methods in Geomech., 23 (12), 1215-1243. Pradhan, S.K., and Desai, C.S. (2002). Dynamic soil- structure interaction using disturbed state concept and artificial neural networks for parameter evaluation. NSF Report, Dept. of Civil Eng. and Eng. Mechanics, The University of Arizona, Tucson, AZ. Pradhan, S.K., and Desai, C.S. (2006). DSC model for soil and interface including liquefaction and prediction of centrifuge test. J. of Geotech. and Geoenviron. Eng., ASCE, 132, 2, 214-222. Roscoe, A.N., Schofield, A., and Wroth, P.C. (1958). On yielding of soils. Geotechnique, 8, 22-53. Sane, S.M., Desai, C.S., and Rassaian, M. (2009). Rate dependent constitutive modeling of lead free solders in electronic packaging. Journal of Electronic Packaging, ASME, submitted. Sane, S.M., Desai, C.S., Jenson, J.W., Contractor, D.N., Carlson, A.E., and Clark, P.U. (2008). Disturbed state constitutive modeling of two Pleistocene tills. Quaternary Science Reviews, 27, 267-283. Shao, C., and Desai, C.S. (2000). Implementation of DSC model and application for analysis of field pile tests under cyclic loading. Int. J. Num. Analyt. Meth. Geomech., 24, 6, 601-624. Somasundaram, S., and Desai, C.S. (1988). Modelling and testing for anisotropic behavior of soils. J. Eng. Mech., ASCE, 114, 1473-1496. Suklje, L. (1969). Rheological aspects of soil mechanics. Wiley-Interscience, London. United States Department of Transportation (USDOT) (1989). Tensar geogrid reinforced soil wall. Experimental Project 1, Ground modification techniques", FHWA-EP-90-001-005, Washington, D.C., Wathugala, G.W., and Desai, C.S. (1993). Constitutive model for cyclic behavior of clays-theory, Part I. J. of Geotech. Eng., 119, 4, 714-729. Wilson, D.N., Boulanger, R.W., and Kutter, B.L. (2000). Observed seismic lateral resistance of liquefying sand. J. Geotech. Geoenviron. Eng., ASCE, 126, 10, 898-906. Zienkiewicz, O.C., Nayak, G.C., and Owen, D.R.J. (1972). Composites and overlay models in numerical analysis of elasto-plastic continua. Proc. Int. Symp. On Foundations of plasticity, Warsaw, Poland. ANALIZA GeOMeHANSKIH PROCESOV ODKOPAVANJA PREMOGA Z VELENJSKO ODKOPNO METODO GREGOR JEROMEL, MILAN MEDVED Ln JAKOB LIKAR o avtorjih Gregor Jeromel Premogovnik VELENJE d.d. Partizanska 78, 3320 Velenje, Slovenija E-pošta: gregor.jeromel@rlv.si Milan Medved Premogovnik VELENJE d.d. Partizanska 78, 3320 Velenje, Slovenija E-pošta: milan.medved@rlv.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 izvleček S poglobljeno geomehansko analizo podzemnega odkopavanja premoga s podetažno širokočelno odkopno metodo lahko podamo medsebojne povezave med fizikalnimi in mehanskimi parametri nastopajočih geoloških materialov v odvisnosti od intenzivnosti odkopavanja premoga. Obseg in intenzivnost odkopavanja premoga ima vpliv na napetostne in deformacijske spremembe v hribinah in premogovem sloju na širšem območju odkopnega prostora. Način podetažnega odkopavanja premoga zahteva večkratne porušitve krovninskih plasti, ki so večkrat ponovno komprimirane in vsaki posebej predstavljajo krovnino pri odkopavanju od zgoraj navzdol. Ponavljajoči procesi rušenja in komprimacije gledani z gledišča teorije plastičnosti predstavljajo relativno manj raziskano področje, saj so pri vsakem tovrstnem procesu prisotne strukturne spremembe naravnih večkrat porušenih in ponovno skomprimiranih materialov v krovnini. Z intenzivnostjo odkopavanja so v neposredni odvisnosti tudi vplivi na druge jamske prostore, ki so locirani v neposredni bližini odkopa kakor tudi v večji oddaljenosti. Obsežne napetostne in deformacijske spremembe v okolici in v jamskih objektih posredno vplivajo tudi na varnost zaposlenih, saj obstaja realna nevarnost porušitve podpornega sistema. Zato je kontrolirano odkopavanje premoga in dobro poznavanje geomehanskih lastnosti vseh nastopajočih materialov ter procesov izjemnega pomena za načrtovanje in vodenje rentabilne proizvodnje ob za delo varnih razmer pri delu zaposlenih pri kompleksnem procesu odkopavanja premoga. Izdelava numeričnega modela, ki omogoča poglobljene analize geomehanskih procesov v krovnini, talnini ter premoškem sloju pri časovnem razvoju podetažnega odkopavanja premoga, ima širši okvir uporabnosti in je izjemnega pomena za ugotavljanje intenzivnosti ter obsega rušnih procesov pri podetažnem odkopavanju premoga, za realna načrtovanja odkopavanja premoga v premogovniku v naslednjih desetletjih ob povečani varnosti rudarjev. Ključne besede odkopavanje premoga, podetažna odkopna metoda, širokočelna odkopna metoda, rušni procesi, metoda končnih diferenc, FLAC, matematični model AN ANALYSIS OF THE GEOMECHANICAL PROCESSES IN COAL MINING USING THE VELENJE MINING METHOD_ GREGOR JEROMEL, MILAN MEDVED and JAKOB LIKAR About the authors Gregor Jeromel Premogovnik VELENJE d.d. Partizanska 78, 3320 Velenje, Slovenia E-mail: gregor.jeromel@rlv.si Milan Medved Premogovnik VELENJE d.d. Partizanska 78, 3320 Velenje, Slovenia E-mail: milan.medved@rlv.si Jakob Likar University of Ljubljana, Faculty of Natural Sciences and Technology Aškerčeva 12, 1000 Ljubljana, Slovenia E-mail: jakob.likar@ntf.uni-lj.si Abstract With in-depth geomechanical analyses of sub-level mining using the longwall mining method we can identify the relationships between the physical and mechanical parameters of geological materials, depending on the intensity of the coal extraction. The extent and the intensity of the mining operations impose impacts on the stresses and cause deformation changes in the rocks and in the coal seams on a broader area of excavations. The method of sub-level coal extraction requires multi-caving of the hanging-wall layers, which are recompressed, and in sub-level stoping each represents a hanging wall. The repeating processes of caving-in and compression, from the aspect of the theory of plasticity, have been relatively little researched because every such process brings about structural changes in natural, multi-caved and recom-pressed materials in the hanging wall. The intensity of the coal extraction has direct impacts on the surrounding and distant mining areas. Extensive stress and deformation changes in the surrounding area, and in the mine, represent a safety hazard for the employees, since the supporting system in the mine roadway could collapse. Therefore, a controlled excavation of the coal, and a good understanding of the geomechanical properties of all the materials and processes involved, is extremely important for planning and managing economic production, while also ensuring safe mining operations. A numerical model that allows for in-depth analyses of the geomechanical processes that occur in the hanging wall, the footwall and in the coal seam during sub-level coal excavation, is broadly applicable and highly relevant for analysing the intensity and the level of the caving processes in sub-level coal mining, and for making realistic plans for coal excavation with workers' safety in mind. Keywords coal mining, sub-level mining method, longwall coal mining method, caving processes, finite-difference method, FLAC, mathematical model 1 INTRODUCTION Caving processes, which are the result of sub-level coal mining, significantly encroach into the natural rock structure and consequently change the deformation and stress fields, and to a certain extent influence the existing objects in the mine. This largely depends on the primary stress situation and the geomechanical properties of the rocks, as well as on the operational time of the mine. All these factors are important, since coal excavation is a complex technological process that needs to function harmoniously in all segments and is vital for normal mine operations. The support systems that provide the stability of mine objects and reduce the impacts of caving processes need to be adapted to the actual geotechnical conditions, such as deflections and changes to constantly changing stress fields, as well as the position and time functions. Understanding the geotechnical properties of geological structures with a different composition of rocks and soils in the hanging wall, footwall and the properties of the coal seams are, in addition to the mining-impact factors, crucial for rational and safe mining. 2 GEOLOGICAL STRUCTURE OF THE VELENJE MINE REGION The Velenje coal region lies in a synclinal valley, which is a depression between the Smrekovec and Šoštanj fault. It is a lignite deposit, formed during the Pliocene era. The depression is meshed with local faults from different ages and runs in different directions. The present valley has been formed by subsiding and simultaneous backfilling. The coal deposit is lens-shaped, 8.3 km deep and 2.5 km wide, with a thickness of more than 160 m and a depth from 200 to 500 m. The deposit is closest to the surface on the edges and deepest in the middle [7]. From the geomechanical point of view, the strength and rigidity of the lignite layer is much higher than the materials in the hanging wall and in the footwall of the deposit. This fact is important, both for planning excavations as well as for designing the most suitable excavation method, which needs to be adapted to the natural geological and geomechanical conditions. Also, the caving-in of the materials in the adjacent hanging wall is different from the subsiding hanging-wall layers lying above, where clay and silt layers are first plastically reshaped, while the layers situated higher up are bent and subsided to the surface, where large depression lakes have been formed as a consequence of the mining activities. 3 RESULTS OF PREVIOUS RESEARCH The basic goal in designing sub-level coal extraction is to obtain the maximum quantities of coal per unit length of the excavation face. In the process of advancement of the excavation in the lower and upper sections, the layers above and behind gradually collapse to a certain height - i.e., caving height - and subside to the bottom of the longwall floor. Depending on the natural and mining technological factors, the layers eventually become compressed, until they reach a sufficient density and hence strength. In geotechnical terms this new rock-soil has specific physical and mechanical properties. The rate of backfilling of the space that has been excavated, and subsiding of the materials from the layers above, depends on the geological and geotechnical features of the rocks, the duration of compression of the caved-in rock layers, until the balance of the primary strength has been reached. There were some studies done in the past in which we tried to determine the distance from the excavation face to the place where the primary stress situation becomes re-established, i.e., to the point where the deformation due to compression develops no more. In-situ measurements of the stresses and deformations are practically impossible due to the inaccessibility of the area. Several researchers have made assumptions with indirect prediction methods: e.g., Wilson [13] assumes that in the part which has been caved in behind the excavation site, vertical stresses linearly increase from the 0 point, reaching the final value, which is at a distance where no more impacts of instability of the mine roadways due to previous mining operations can be detected. Later on, he found that the arch of the rock mass pressure changes in intervals between 0.2 and 0.3 times the depth of the excavation under the surface. Some other authors, e.g., King and Whittaker [6], Choi and McCain [2], and Mark [9], suggested similar interval ratios, i.e., 0.12. It should be noted that some studies based on numerical methods showed completely different results. For example, Trueman and Thin [12] proved that the impact area of the coal extraction decreases with the increase of the depth of the excavation NCB [11], and Figure 1. Geological structure of the coal ACTA GeOTeCHNICA SLOVENICA, 2010/l deposit in Velenje - west-east direction. after extensive measurements and research in situ, they found that the ratio between the size of the compression pressures and their geometrical distribution in relation to the primary pressure stress depends to a great extent on the geotechnical properties of the caved-in material, which is reasonable since the subsidence of the hangingwall layers in the old part depends on these properties. 4 CHARACTERISTICS OF THE VELENJE SUB-LEVEL MINING METHOD (VELENJE MINING METHOD - VMM) The Velenje mining method (VMM), which is characterised by a continuous caving-in of the hanging-wall layers, started in 1952. Based on continuous research, this method has undergone various technological improvements. With a better understanding of the geomechanical processes in coal mining, several improvements have been made in order to increase the production. The basic characteristic of this method is that the area of the coal excavation is extended to the area above by means of a hydraulic shield support. With this method we can make use of the primary and secondary stress conditions and their transformations, which cause breaking and crushing of the coal seam, plastic transformations of the protective clay layer and bending and subsiding of soil layers lying above. Figure 2. Longwall operation in the Velenje coal mine. Due to its geological and geomechanical characteristics, as well as other technological mining features, the Velenje mining method is considered as a more complex method in terms of ensuring safety according to European and world standards. Modern mining equipment has been developed to operate safely and to contain explosions. It is composed of sections with hydraulic support, a face conveyor, with standard or side charging, an electronically operated cutting machine and conveyance with a crusher. The main operations involve cutting the excavation wall on the longwall, shifting the chain conveyor, shifting the hydraulic support, and coal extraction from the upper excavation section [7]. After several cuts and shifts have been made in the lower excavation section the coal extraction from the upper excavation section runs in subsequent steps. With the existing equipment, this means 2 to 3 cuts through hard and tough coal, and 3 to 6 cuts in the caved-in and crushed coal. The width of the excavation round depends on the basic dimensions of the excavating equipment and the geomechanical properties of the coal. A controlled extraction of the crushed coal from the upper excavation section can be reached by sliding the coal down via bent roof beams to the front part of the lower excavation section [7]. VMM is a completely mechanised method, but with some automated phases. The method can offer high productivity if natural factors of caving-in in the upper excavation section are provided. In terms of organisation and technological solutions, the Velenje longwall mining method has been constantly upgraded since its beginnings and has become well known worldwide and considered as the most important method of excavating thick coal seams. In technological and organisational terms this method is constantly being optimised, particularly in terms of increasing the production, increasing the yield of a seam, the safety of the employees and the humanisation of the work. 5 LABORATORY ANALYSES OF THE GEOMECHANICAL PROPERTIES OF ROCKS Numerous studies and measurements of various geomechanical and hydrogeological parameters have been made so far, both on the laboratory scale and in situ [10]. It was found that the mechanism of compres- mjs* Csi jl shearer hydraulic shield support Figure 3. Excavation scheme of the VMM [5]. sion of the caved-in coal and clay, which is the result of the direct caving in of the hanging-wall material above the longwall mining, is directly related to the time variable. e=e(t) + e(p) (1) where e - is the vertical strains of the caved hanging-wall ground layers behind the excavation (/), t - is the time (days), p - is the vertical hanging-wall pressure (MPa). Fig. 4 presents the caving-in and recompression processes in the hanging wall behind the advancing longwall face for a twice-caved-in hanging wall. The same situation would occur in thrice-, or multi-caved-in hanging walls. A laboratory analysis [10] was made to show the changes in the volume and the shear strength of the crushed rock. During the compression of the materials from the hanging wall, using the instrument in Fig. 4, the main impact was due to vertical pressures, which were identical to those which we found in situ. The analysed values of the shear strength showed that the clay layers are well homogenised and with a relatively high cohesion and a low angle of inner friction, while crushed coal, with practically zero cohesion, has a relatively high angle of inner friction. © 0 (minj (T ¡-O (min) G i=0 (min) 0 (min) (D —> 0" "., a*" -> a""* O 1 —> 0""» —> ct"". Figure 4. Changes in the volume properties of the multi-caved-in hanging wall behind the advancing longwall face. Figure 5. A device for measuring the shear strength of the caved ground mass at high vertical hanging-wall pressures [10]. A - a sample of caved ground mass, B - piston, C - shear box filled with sample , D - vertical cylinder for forming vertical load, E - cylinder for forming shear strain in the sample, F - displacement transducer for measuring shear strains, G - hydraulic pump, H - manometer. Measurements of the compression properties of crushed clay 1 50 100 150 200 t (min) - time 250 300 Vertical stress 0 = 1.38 M Pa -m- 0 = 2.77 MPa -A- 0 = 4.15 MPa -*- 0 = 5.55 MPa a = 6.93 MPa 350 Figure 6. Changes of compressive strains versus time at different normal loads of caved ground of clay [10]. 0 Fig. 6 presents the compression of caved clay from the hanging wall with compressive strains after exerting a load with normal stress in relation to the compression time. As can be seen from the graph, 98% strain under pressure stress conditions occurs instantly, depending on the vertical load, and only 2 % during the remaining time. From these results we can infer that the vertical ground mass load plays the principal role in the compression of caved-in clay and that the rheological impact is small. A more detailed analysis shows that there is a logical relationship between the deformations and the intensity of the load; the vertical load increases rapidly as the excavation advances from the observation point. In addition to this, our analyses and observations to date have shown that multi-caving in the hanging-wall layers (in some cases also five-times caved the same hangingwall ground) does not depend on the basic structure of the caved ground mass but on the composition of finegrained clay particles. Measurements of the compression properties of crushed coal, mixed with the hanging-wall overburden layers, have shown that compression processes are similar. The process also depends on the intensity of the vertical load of the hanging-wall layers, which continuously subside onto the caved hanging-wall layers. In this case too, the intensity of the vertical pressures, which continuously increase behind the advancing longwall, have the most significant impact on the compression process. 6 GEOMECHANICAL MEASUREMENTS In order to obtain a real picture of the impact of excavations on the surrounding rocks and adjacent mine objects, which are the constituent parts of the system of longwall mining, the secondary stress conditions were monitored by measuring the compressive stresses. We measured the time-dependent secondary stresses using Gloetzl pressure cells, which were installed in front and behind the excavation face to monitor the stress changes in relation to the distance to the longwall face. The measuring points and the methods of collecting the data had to be carefully selected due to possible risks of methane, which require special safety precautions. Since in geomechanical terms the caving processes in single-or multiple-caved hanging-wall ground layers are quite complex events where natural materials are completely cave. Measurement data taken in the mine are of very important because they provide the basic information for a realistic and good-quality analysis of the compression processes. With in-situ measurements we are faced with secondary impacts, as well as measurement errors, which among other factors, occur due to the so-called "scale effect". The sizes of the particles in the hanging wall, which have been caved-in behind the advancing excavations, are to some extent unpredictable. Regardless of the fact that the caving process is continuous, this means that behind the longwall in the goaf there are no empty spaces left. The measurements in the mine are carried out based on the previous, long-term experiences of professionals. Previously, many types of measuring cells were in use in order to overcome the "scale effect"; however, the most reliable are the Gloetzl pressure cells. Fig. 8 shows the measuring principle of the pressure stresses in the region behind the advancing excavation, where the compression processes are the most intensive. Measurements of the compression properties of crushed lignite 100 150 200 t (min) - time Vertical stress ■ a = 1.38 MPa a = 2.77 MPa a = 4.15 MPa a = 5.55 MPa a = 6.93 MPa 350 Figure 7. Changes of the vertical strains with time under different normal loads of caved lignite [10]. Figure 8. Principles of measuring the increase of the compressive stresses within the region of the advancing mining. The decisive factors in analysing the compression processes are the results obtained from measurements of the increasing vertical stresses in the caved hanging wall behind the excavation. Of particular importance are the parameters related to changes in the vertical stresses, with regard to the distance from the measuring point at the longwall. It is clear that an understanding of the caving process depends on the results of measurement deformations in the coal in the hanging wall related to the distance from the caving and geomechanical properties' changes after the caving-in and re-compression phase, and are the most important for a determination of the mine-exploitation parameters for the next sub-level panel. The conditions under which measurements for analysing the compression process of the caved-in hanging wall are performed are complex and demanding. The reasons are as follows: the inaccessibility of the measuring points, the distribution of the measuring equipment around the excavation site, the complex installation of the equipment, the methane conditions, the extreme deformations and the large shear displacements that hinder the drilling process and the installation of the probes, which all consequently shorten the time of the observations. 7 NUMERICAL MODEL We made a special 3D analysis of the stress conditions in the surroundings and a 3D numerical analysis of the excavation of the level. The size of the finite-difference mesh was 300 m x 300 m x 300 m. The distribution of the elements reached a depth of 2 m. Thus, the model consisted of 355,500 elements. To draw the finite-difference mesh we used the Autodesk AutoCAD software application. It allowed us to make a transformation of the basic geometry of the model into a recognisable form for the FLAC 3D code [3]. The size of the primary stress conditions that we used in the 3D analyses was determined on the basis of in-situ measurements carried out in previous years. The vertical component of the primary stress was: Pv (2) where: yi - is the unit weight of each hanging-wall ground layer (kN/m3), H' - is the thickness of each hanging-wall ground layer (m). Measurement results of vertical stresses in the caved hanging wall in the goaf Distance into the goaf [m] Figure 9. Measurement results of the vertical pressure increases in the caved hanging-wall ground layers in the goaf behind the longwall [4]. The ratio of the vertical stress is expressed as a sek / apri, where a sek is the secondary stress and apri is the primary stress. Figure 10. The finite-difference model with mesh using the FLAC 3D code. Figure 11. Geometrical parameters in modelling the area of the coal exploitation. 38. ACTA GEOTECHNICA SLOVENICA, 2010/1 The ratio between the vertical and horizontal components of the primary stresses was 0.85, determined by more measurements in previous years. The geometrical design of the excavation on the level consists of the segments of excavation on the lower excavation section, with lengths of 4 m, and the upper section, with lengths of 11 m. To simulate the development of the coal exploitation we used, for every two cuttings into the longwall face, the amount of one-third of the upper excavation section to reach the required 15 m excavation height. The scheme in Fig. 11 shows the method of using the geometrical parameters for constructing the model and for modelling the excavation. 7-1 RESULTS OF THE COMPUTER SIMULATION OF THE LONGWALL ADVANCEMENT The numerical simulation of the advancement of the sub-level coal exploitation in the upper and lower excavation sections in the ideal hanging wall, coal seam and footwall ground-layer conditions, showed a relatively good agreement with the values measured before, particularly when we compare the stress changes within narrower areas of the longwall advancement. From the results of the calculations we can conclude that in front of the longwall excavation face, which was moved ahead, as well as on both sides of the longwall panel, the vertical compressive stresses increase significantly, while the intensity of the stresses drops quickly with the distance from the longwall excavation face. Also, the vertical stresses in front of the excavation start increasing at x = 4h to 5h (h = excavation height ) in front of the excavation face, while in the region with the distance of xmin = 0h to 2h excavation height in front of the excavation face, the vertical stresses reach the highest values, as shown in Fig. 13 (see next page). A marked increase of the vertical stresses in front of the excavation face represents approximately 180 % of the vertical stress component in the primary stress conditions. This value does not change in computer simulations of the advancement of mining, which means that it is constant at a rate of advancement of 6 m/day. Figure 12. Distribution of the vertical stresses in the area of coal exploitation. FLAC3D 3.00 Step 131178 Model Perspective 14:4;: 18 Wed Jul 08 2009 Center: <: 1.5006+002 : 1.5003+002 Z: 1,500e+002 Dlst 3.735e+002 Rotation: ■-:: 90.000 Y: 270.000 Z: 0.000 Mag.: 1 Aria.: 22.500 Contour of SYY M agfac = 0.000e+000 Gradient Calculation H-1,8201 e+007 to -1.7850S+007 ■ -l.7000e+007 to-1.6150e+007 □ ■I .58006+007 to -1.44506+007 H-1.38006+007 to-1.27506+007 Q-l19006+007 to-1.10506+007 R-1,0200e+007 to-9.3500e+008 -8.50006+005 to -7.55006+008 -8.80006+006 to-5.95006+008 -5,10006+005 to-4.25006+008 H -3,40006+005 to -2.55006+008 ■ -1.70006+005 to-8.50006+005 U O.OOOOe+OOO to 1.81776+005 Interval = 8,5e+005 Itasca Consulting Group. Inc. Minneapolis. MN US- u Figure 13. Distribution of the vertical ground stresses in the wider region of the longwall coal exploitation. Figure 14. Distribution of the vertical stresses in the region of longwall coal production. ACTA GCOTeCHNICA SLOVENICA, 2010/1 Later on, the results of the present analyses show that in the region behind the advancement of the mining there is a large distress field. Based on a computational analysis, zero additional stresses at the bottom of the longwall and close to the beginning of goaf were found. A gradual increase in the vertical stresses occurs during the process of compression until it reaches the values of the vertical components of the primary stress state. In the pillars on both sides of the longwall the calculations show a marked impact of the distress zone, which reaches a length of 15 m as a result of the coal production. In the inward direction the values jump up, reaching a peak value of the stresses, i.e., 170 % of the primary stress state, while afterwards, the impact of the stresses decreases, eventually reaching the values of the primary stress field. 7.2 COMPARISON OF THE RESULTS OF CALCULATIONS USING THE 3D NUMERICAL MODEL WITH THE MEASURED PARAMETERS The analysis of the stress and strain changes in the thick coal seam simulating the sub-level coal exploitation in a large area shows that the steep increase of stresses, i.e., 2-3 for the excavation height in front of the excavation face and a steep decrease of stresses 0-2 for the excavation height in front of the excavation face, reflect well the natural conditions and the measured values. It is clear that the stresses increase by approximately 80 % of the primary stress state, which in fact is the largest vertical secondary stress, i.e., 1.8 x y x H. From the results of the present analyses the conclusions indicate that behind the longwall excavation face, the stresses increase gradually and steadily and reach the value of the primary stress state at a distance of approximately 17.5 times the excavation heights (if the excavation height is 10 m, the distance is 175 m). The calculated parameters, e.g., the additional vertical pressure in front and behind the excavation face and the vertical strains in the hanging wall layers and in the coal seams, are in relatively good agreement with the measurement results. For a clear understanding the measured values of additional ground pressures in front and behind the longwall face are carried out carefully and with proper confidence. In the middle of the pillar the stresses increase and sometimes reach as high as 220 % of the primary vertical stress, which is also in agreement with the in-situ measurements. However, the distribution of vertical pressures, which was found in the intermediate pillar based on the results of the numerical simulation, has an Figure 15. Distribution of the vertical ground stresses and the vectors of ground movements in the wider area of coal exploitation. Changes of vertical stress conditions in front and behind the longwall face re u 01 > o o CO o: A 1 5 5 1 U.O -150 -50 50 150 Distance from the longwall face [m] Measured values of vertical stress: h = 10 m, v= 4 m/day - Computer simulation: h = 15 m, v= 6 m/day Figure 16. Vertical stress-change measurements versus the distance from the longwall face. Ratio of the vertical stress is expressed as a sek / apri, where a sek is the secondary stress and apri is the primary stress. informative character because the model was simplified. The distribution of additional vertical pressures in the intermediate pillar, which was measured partially by pressure cells, differs from the calculated values obtained by the 3D numerical simulation. The visual presentation of the results of the measurements and the numerical simulation of the excavation in the Fig. 16 gives an actual picture of the events in front of the longwall face and in the goaf behind it. The shape of the curve of the vertical stress changes obtained with the computer simulations takes the form of a curve obtained by the measurements; however, we cannot compare it with the relation to the distance from the caving, since the computer simulation can only give computational steps. 7.3 MATHEMATICAL MODEL Numerous analyses of the stress conditions in-situ and further confirmations by numerical simulations, allowed for the application of a constitutive model. This constitutive model describes the changes of the stress and deformation conditions in relation to the advancing speed of the coal production and the distance from the longwall face. A modified form [14] also includes the impact of the speed of the advancement of the coal extraction. The distribution of the stresses in the goaf is presented using the following equation: where: is the vertical stress (MPa) at a distance of lx(m), is the unconfined compressive strength (MPa), is the initial bulking factor (/), is the maximum surface subsidence (m), is the mining height (m), , are the coefficients that depend on the strata lithology composition of the hanging-wall ground layers, is the distance from the longwall excavation face (m), is the speed of the longwall advance (m/day). Table 1. Values of the material parameters of the marl for primary and secondary creep after loading. Strata Unconfined Coefficient lithology compressive strength C3 C4 hard > 40 MPa 1.2 2 medium 20-40 MPa 1.6 3.6 soft < 20 MPa 3.1 5 The maximum allowable deflection of the hanging wall is expressed by the equation: gHb?J(hb - 0.05h1'2) 0.4h — " 10.39sL042(k-1) + gHb& -0.05/Î12 (4) 10.39s!.042 [b-1] S. + h+h - 0.05h1'2 c3h + c4 1 - e X X [hb87 + 0.05hL2b7'7 ]- b8'7 S. + h+ - »05h"" c3h + c4 c X where H is the depth of mine exploitation working (m), Y is the unit weight of the hanging-wall ground layer (kN/m3), b is the initial bulking factor (/), ffc is the unconfined compressive strength of the ground layer (MPa), h is the mining height (m), c3, c4 are the coefficients that depend on strata lithology in the hanging wall. Fig. 17 shows how the vertical stresses would increase in the goaf with the distance from the observation point, using the proposed mathematical model. The distribution of the stresses in front of the face is expressed with the following equation: sX1 = —0.0000064hx4 - 0.00032hx3 - (5) - 0.022hx + (1 + 0.045h)h01 and with the equation: sx2 = —0.0003x4 — 0.0128x3 — 0.2101x2 — 1,2564x — 0.0019 (6) where ox1 is the vertical stress in front of the longwall face at the distance x1 (MPa), ox2 is the vertical stress in front of the longwall face at the distance x2 (MPa), h is the mining height (m), x is the distance from the longwall face (m). Changes in the curve occur at the point where: sx1 = sx 2 ^ x ^ h = 15m ^ x = —5h (7) h = 10m ^ x = —2.75h This curve shows the location of the prediction points where the ground pressure starts to increase in the surrounding areas of the longwall face. The graph in Increase of vertical stresses in the goaf 1.2 1 0.8 0.6 0.4 0.2 0 100 200 300 400 Distance into the goaf [m] 500 600 -1 m / day -2 m / day c %), the Shear strain rate, the Soil structure, and the Loading history, etc. In this paper, undrained, direct shear tests were conducted to study the effect of the plasticity index (PI) and the normal stress (a) on the shear behavior and the shear modulus of remolded clays. The results show that the normalized shear modulus at a constant strain will generally increase as the a and PI increase, and the common empirical equations for undisturbed soils at y = 0~0.1 might be applicable for the disturbed soils too. G = 625- OCRk Keywords plasticity index, normal stress, shear modulus, disturbed clayey soils 0.3 + 0.7e' S axs (1) G = 6908(2.17 e)2 S2 (2) 1 + e G = 3230(2.97 e)2 S1 (3) 1 + e Where pa is the atmospheric pressure and K is a constant coefficient depending on PI. The laboratory reports presented by Humphries & Wahls (1968) [2] showed that equation (2) is usable for normally consolidated clay with an ordinary sensitivity. The effects of primary and secondary consolidation on the shear modulus of soil were studied by Hardin & Black (1968) [1] and Humphries & Wahls (1968) [2], and it was revealed that the shear modulus increases during the primary and secondary consolidation. The relationship between the shearing modulus and the other parameters was suggested by Hardin & Drenvich (1972) [3] as equation (4) for undrained, over consolidation clay with an ordinary and relatively low sensitivity. G = 3230(2'97 6)2 (OCR)K 1 + e i t' 2 (4) where K is the PI-dependent coefficient. These researchers presented a relationship for the maximum shearing module and the actual stress and strain of the soils, according to equation (5), and amended it for the dynamic state as a general equation (Eq. 6) for undisturbed soils. OCRk (5) 1 G max+ 7 G = G 7h = (1 +7 ) 1 + ae (6a) (6b) where a and b are constant values depending on the type of soil and the number of cycles of loading and yr is the base strain, which is defined as 7r = G max Wahls & Marcuson (1972) [4] studied the effects of confining the effective stress on the variation of the maximum shearing modulus and suggested equation (7), based on the void ratio (in the range between 1.5 and 2.5). G = 445(——e)l. i 1 + e c (7) where cC shows the effective confining stress. Kokusho et al. (1982) [5] performed some triaxial dynamic tests on the undisturbed samples and presented equation (8) for the normally consolidated clays in the restrictive limits e = 1.5~4 and y <10-4. G = 90 (732 - e)2 . 6 1 + e c (8) Richard & Athanaspoulos (1983) [6] studied the effects of the confining pressure of cohesive soils in undrained conditions on the maximum shearing modulus. Based on their observations, it was demonstrated that the lowering of this pressure causes a sudden fall of Gmax. The percentage of moisture is another important factor that affects Gmax variations, and researchers such as We et al. (1984) [7], Qian et al. (1991) [8] and Marinho et al. (1995) [9] studied the effect of this parameter on the Gmax variation and presented the moisturizing limit when the soil reaches its maximum shearing modulus. Studying the effects of PI and OCR on the dynamic shearing modulus of saturated soils was performed by Dobry & Vucetic (1991) [10], and it was shown that by increasing PI in a constant strain, the value of G/Gmax increases. An increase in the shearing strain speed in normally consolidated clay soils causes a decrease in the dynamic, undrained G/Gmax value. Lanzo et al. (1997) [11] performed some dynamic tests on remolded clay and sand samples and showed that the diminishing diagram of Gsecant/Gmax vs. Y shifts upward by increasing the normal stress and the ratio of over consolidation. In other words, in a constant shearing strain, by increasing the normal stress and the ratio of over consolidation, the Gsecant/Gmax increases. In the clay and sand samples, the increase in the shearing strain speed increased the secant shearing modulus in the undrained condition. Zhou Yan Guo (2004) [12] made some tests on saturated clay samples and showed that irregularity in the soil structure leads to a decrease in the hardness and the undrained resistance of the soil. A series of dynamic tests on large-grain soils was conducted by Kalinski & Hardin (2005) [13] and showed that the dimensions of the granules are effective in the maximum shearing module. Zhou Yan Gou et al (2005) [14] showed that the effect of the history of dynamic loading is not less than the effect of the history of non-dynamic loading in saturated sands in undrained conditions. Romo & Ovando (2006) [15] presented a relation that showed that the frame of the apparatus subject to the test could be effective in calculating Gmax. Bergado et al. (2006) [16] presented the case history of a laboratory evaluation of the interface shear-strength properties of various interfaces encountered in a modern-day landfill with the emphasis on a proper simulation of the field conditions and the subsequent use of these results in a stability analyses of a linear system. Li (2007) [17] implemented the nonlinear shear strength criteria of a power-law type in a finite-element slope-stability analysis. Guetif et al. (2007) [18] proposed a method for evaluating an improvement of the Young's modulus of the soft clay in which a vibrocompacted stone column is installed. From the numerical results the degree of improvement of the Young's modulus of the soft clay was estimated. Also, the zone of influence of the improved soft clay was predicted. Basudhar et al. (2007) [19] made an experimental study on the circular footings resting on a semi-infinite layer of sand that was reinforced with geotextiles. The study highlighted the effect of the footing size, the number of reinforcing layers, the reinforcement placement pattern and the bond length, and the relative density of the soil on the load-settlement characteristics of the footings. Sabermahani et al. (2009) [20] conducted a series of 1-g shaking-table tests on 1-m-high reinforced-soil wall models. The distribution of the shear-stiffness modulus (G) and the damping ratio (D) of the reinforced soil along the wall height were assessed. Most studies performed by other researchers discussed the effects of other parameters on the shearing modulus or cover the effects of PI and a on undisturbed samples and in dynamic conditions; therefore, for studying the effects of PI and the normal stress on the static shearing behavior of undisturbed clay and studying the reliability of the existing equations on undisturbed clay soils, in this paper, the results of laboratory tests performed on samples with different values of PI and various normal stresses have been presented. An admissible calibration has been established between the dynamic equations and the results of the static test. In fact, the main purpose of this research is to make a calibration between the dynamic equations and the dynamic tests results of previous investigations, and the results of static tests that have been performed in this study. However, the results of this investigation can be used in different scopes of geotechnical engineering and in various soil structures, such as embankment dams that contain solid structures. For example, a concrete culvert, where plastic clay has been used at the interface of the concrete and the soil; or in water-isolating trenches with plastic clay. 2 PROPERTIES OF SAMPLES AND INITIAL TESTS A soil that indicates the specifications of clay soils used in the practical project was used as the base soil. For making samples with different values of PI, the base clay soil was mixed with different percentages of bentonite and after passing it through a sieve number 40, some tests were performed on samples in order to determine the relative density (Gs), the hydrometer granulation, the Atterberg limits (LL, liquid limit, PL, plastic limit and PI, plasticity index) and the modified compaction of soil in accordance with ASTM standard. The weight percentages of bentonite along with the identifying tests are listed in Table 1 and a diagram of the hydrometer granulation is shown in Fig. 1. It should be noted that the percentage of the bentonite inserted in Table 1 shows the weight percentage of the bentonite used relative to the final weight of the reconstructed samples. 3 PROCESS OF LABORATORY RESEARCHES TO STUDY THE SHEARING PROPERTIES In this section, the method of making samples and the stages of direct shear tests are introduced. Each one of the remolded samples with different values of PI was compacted in the direct shear apparatus with 95% of dry density and with the optimum moisture in different layers according to the ASTM D-3080 standard and was immersed in water to become saturated. After saturation, the samples were loaded to remove any swelling and the samples returned to their initial dimensions (the amounts of stresses that remove the swelling are shown in Table 2). After the obviation of the swelling, normal loading was made on the three frames made of different values of PI and were immediately sheared in the direct shear apparatus at high speed. The validity of the direct shear test on saturated clay has been presented in different texts (Bowles, 1968 [21] and 1984 [22], Murthy, 1977 [23], Punmia et al. 2005 [24]). In order to determine 0 01 0 1 Particle size (mm) Figure 1. Hydrometer grain-size distribution test results of reconstructed samples. the rate of load growth, the Ra coefficient is defined in equation 9. Normal stress in direct shear test Table 3. Normal stress for different R„. Stress that removes swelling (9) Table 1. Result of identifying test. Result Sample 1 Sample 2 Sample 3 Bentonite Percentage 0 10 20 Soil Type CL CL CL Liquid Limit (LL) 28 36 43 Plastic Limit (PL) 17 19 20 Plastic Index (PI) 11 17 23 Gs 2.65 2.63 2.60 "opt (%) 15.0 13.7 12.6 Yw (kWm3) 18.42 19.13 20.24 The ratio of the increase in load for three frames of each sample is considered to have constant values equal to 1.44, 2.88 and 4.32. The saturation moisture and the amount of stress for removing the swelling is shown in Table 2 and the normal stress for different Ra (different samples) is specified in Table 3. Table 2. The saturation moisture and the stress for removing the swelling in the direct shear tests. Sample Saturation moisture Stress for removing (%) swelling (kN/m2) Sample 1 19.9 136.3 Sample 2 22.4 272.6 Sample 3 24.8 354.0 Sample Normal Stress (kN/m2) (R = 1.44) Normal Stress (kN/ m2) (R = 2.88) Normal Stress (kN/ m2) (R = 4.32) Sample 1 196.13 392.26 588.40 Sample 2 392.26 784.53 1176.80 Sample 3 509.94 1019.89 1529.84 4 RESULTS AND DISCUSSION The results obtained from the samples' identification index tests and the engineering tests analysis as explained before will be discussed in this section. 4.1 RESULTS OF THE IDENTIFICATION INDEX TESTS In this section, in order to obtain more suitable and accurate results, the values obtained from the testing on samples with 30% and 40% bentonite are presented too. 4.1.1 Effects of the Percentage of Ben-tonite on the Atterberg Limits of the samples As observed in Fig. 2, by increasing the percentage of bentonite in the samples, the Atterberg limits (LL, PL, and PI) are increased. Based on the obtained results, the relationship of the Atterberg limit and the percentage of bentonite of the samples is a relatively linear relation and the rate of increase in the liquid limit with an increase in 60 50 40 o, J Cm 30 20 10 ♦ Liquid Limit LL = 0.66(B%)+29 R: = 0 988 ' ■ Plastic Limit Plastic Index ♦ PI = 0.48(B%)+ 11 6 R- = 0 981 PL= 0 18(B%)+ 16 6 R: = 0 98" 0 ? 10 15 20 25 30 35 40 45 Bentoiiite Percentage (Bi o) Figure 2. Atterberg limits versus percentage of bentonite. the percentage of the bentonite of samples exceeds the rate of increase in the plastic index and the plastic limits. The relations obtained based on a linear interpolation among the values of the Atterberg limit and the percentage of bentonite are shown in Fig. 2. It should be noted that a 10% increase in the amount of bentonite of the samples will cause an increase of 6.45, 1.69 and 4.76 in the liquid limit, the plastic limit and the plasticity index, respectively. 4.1.2 Effect of pi on the Optimum and saturation Moisture Based on the obtained results, the increase in the PI of the samples causes an increase in the percentage of the optimum and saturation moisture. It can be seen in Fig. 3 that the process of the increase in the optimum and saturation moistures vs. increase in the PI value follows linear relations and for each 10 percent increase in the PI, the optimum moisture increases 2.51% and saturation moisture increases 4%, on average. The linear formula between the optimum and saturation moisture, and the percentage of bentonite is shown in Fig. 3. 4.1.3 Effects of pi on the Maximum Dry Density Based on the results of the modified compaction tests, an increase in the values of PI in the samples caused a decrease in the maximum dry density. As can be seen in Fig. 4, the variation of the parameters is almost linear, and for each 10% increase of PI, the maximum dry density decreases by 1.12 kN/m2. 10 15 20 ^ Plastic Limit. (%) Figure 3. Optimum and saturated moisture versus different percentage of PI. 20 30 Plastic Limit ("») Figure 4. Dry density versus different percentage of PI. 4.2 OBTAINED RESULTS FROM DIRECT SHEAR TESTS 4.2.1 Effects of pi on Shearing Resistance Parameters Fig. 5 shows the effect of PI on the shear strength parameter. As is clear, an increase in the values of PI leads to an increase of the cohesion and a decrease of the internal friction angle. The relationship between the percentage of bentonite and the cohesion and internal friction angle reveals that for each 10% increase in the amount of bentonite of samples, the cohesion increases 40 kN/m2 and the internal friction angle decreases by 9.30 degree. The increase in the cohesion and the decrease in the internal friction angle might be attributed to the increase in the clay minerals and the placement of bentonite particles among the large granules of soil. 4.2.2 Effects of Ra on the variation of the G/Gmax versus y The variation of the normalized shear module (G/Gmax) vs. the values of the shear strain (y) for different Ra is shown in Fig. 6. It can be seen that with an increase in the amount of shear strain (y), the normalized shearing module (G/Gmax) decreases. On the other hand, in a constant shear strain as Ra increases, there is an increase in the amount of G/Gmax . In Fig. 6(a) the difference between the variations of G/Gmax in a constant shear strain for different Ra is shown. A comparison between these two values shows that as the normal stress increases, the change in the amount of G/Gmax in a constant strain decreases, and this process continues until the failure point is reached. 4.2.3 Effects of pi on variation of g/ Gmax versus v max 1 The variation of the normalized shearing module (G/Gmax) vs. the shear strain values (y) for different values of PI is shown in Fig. 7. As is clear, in a constant Ra and a constant shear strain, the samples with bentonite showed less G/Gmax than the samples without bentonite (which has been referred to as base soil). The reason for this could be because of the imbalance made in the soil structure and placement of bentonite among its particles. In addition, a comparison between samples 2, 3 (mixed with bentonite), shows an increasing process of G/ Gmax as a result of the increase in the value of PI at a constant shearing strain. 4.2.4 variation of Gmax versus pi and Ra Fig. 8(a) and (b) show the relationship between Gmax and PI and Ra , respectively. As is clear, with an increase in PI, the Gmax increases and the effects of Ra on the process of this increase in such a way that in constant PI the increase in R„ increases (the hardness of soil). The variation a max x ' of vs. PI and R„ is shown in a 3D contour (Fig. 9). max a 4.2.5 verification of the Test Result In order to determine the efficiency of the existing reliable relations on disturbed soils and comparing the variation of the G/ Gmax versus y resulting from these relations, and the laboratory results, a new parameter (A) is defined as below: Figure 5. Effect of PI on the shear-strength parameter. a o a) o o a) 0.00001 0.0001 0.001 0.01 Shfiai ill uji ■ o.ooi o.oi S]:--;ii stümi (-y) b) b) -PI -11,0 = 332 21) k».'ju2 -PT = 17,0 = 784.53 kN,'in: •Pl-23,0=1019.89 kN 1112 0.001 0 01 Sa^.n -Irion I' I 0 001 0.01 .Sj;U s'l.lll. {*() c) -Rcr—I 32,a-1529 84kN/ilil2 J.........1 0.00001 0 0001 0.001 0.01 Sheai fliMli ■ Figure 6. Variation of normalized shear module versus shear strain for different R , (a). PI=11, (b). PI=17, (c). PI=23. <5 ö c) -PI= 11, G- 588.40 kN/ml -PI- 17,0- 1176.80 LN/in2 PI= 23, cr= 1529.84 kN/m2 0 00001 0.0001 o.ooi 0.01 S.i'..11 .■■■...Ii Figure 7. Variation of normalized shear module versus shear strain for different PI, (a). R=1.44, (b). R=2.88, (c). R=4.32. %D = G G G G max (Eq) max (lab) G G max( Eq) x100 (10) It should be noted that the Hardin and Drnevich (1972) [3] equations (Eq. 5, 6), which are known as reliable equations, were used in this comparison. In Fig. 10 to 13, the variation of A versus y for different values of Ra and PI are shown. With respect to these figures, for strains in the range 0~0.1, the A has relatively small values that show the reliability of the presented methods of this study for disturbed soils in small shear strains (less than 0.1); however, in higher shear strains (more than 0.1) the value of A increases rapidly. 4.2.6 effect of Ra on the variation of a versus y The increase of Ra (which leads to an increase in the normal stress) for constant PI, leads to a rapid growth in the value of A, which begins from greater strains. This shows the greater compatibility of the Hardin and Drnevich (1972) [3] relations with the lab tests for an increase in the normal stress (Fig. 10 and 11). 4.2.7 effect of pi on the variation of a versus y As was shown in Fig. 12 and 13, in sample 1 with no bentonite addition, the value of A begins to grow at larger shearing strains than the two other samples (bentonite containing). The process of these changes could be attributed to the lack of bentonite particles among the base soil particles, the maintenance of the initial structure of the soil and subsequently, greater compatibility between its behavior with the Hardin and Drnevich (1972) [3] relations. A comparison between samples 2 and 3 in Fig. 12 and 13 shows that with an increase in the values of PI, the size of A starts the rapid increase from the larger strains, which is a proof of the good compatibility between the Hardin and Drnevich (1972) [3] equations and the laboratory tests results. 5 MATHEMATICAL FORMULATIONS From the previous sections it can be seen that the plasticity index and the normal stress are the most effective parameters for Gmax . In this section, one- and two-variables functions are introduced to formulate the Gmax . Equations 11 and 12 are single-variable functions that show the relationship between Gmax, PI and a. These functions are valid for the following values of a (kN/m2) and PI: 196 < a < 1530 and 11 < PI < 23. Gmax = —0.674PI2 +47.84PI + 908.5 (11) Gmax = 0.322 c( / ) n = 1,2,3..., N (5) In this equation the only unknown is the phase velocity c(f), which can be obtained from the inversion of the autocorrelation coefficients. Two assumptions are made for those methods. Firstly, the Rayleigh wave microtremors are dominated by the fundamental mode and, secondly, the subsoil structure under the array of observation is parallel. These methods cannot discriminate between the different modes (fundamental and higher) but will only show the predominant one (Okada, 2003). The SPAC method requires a circular array for data collection, but in urban areas it is difficult to deploy a circular array. The benefit of the ESAC method is that it does not require a circular array and is applicable to arbitrary arrays (L-shape, T-shape, cross-shape). A significant advantage of the SPAC and ESAC methods over the ReMi method is the extraction of the scalar wave velocity irrespective of the direction to the source and an omni-directional wave field (that is typical for urban areas) leads to better estimates of the scalar velocity (Asten, 2001). The identification of the dispersion curve in the SPAC and ESAC methods is not performed by picking the phase velocity and thus the results are not exposed to a subjective choice. 5 DATA ACQUISITION AND ANALYSES Microtremor array surveys pose practical challenges in urban areas as very little open space is available for setting up the arrays. Thus the proper locations within the city were selected on the basis of aerial orthophoto images and field observations. The study area is limited by the Ljubljana highway ring and has an area of about 55 km2. The selected locations correspond to city parks, school playgrounds and agricultural areas. The microtremor records were acquired at thirty such sites (Figure 12). The microtremor array surveys were conducted by using a SoilSpy Rosina (Micromed) multichannel seismic digital acquisition system that offers a good signal-to-noise ratio (Micromed, 2008a). Twenty-four 4.5-Hz vertical geophones were planted to form a 2D array. For the ESAC analysis we used L-shape arrays and for the ReMi analysis the longer linear legs of the same arrays (Figure 4). Our task was to image the shear-wave velocity of the subsurface layers down to at least 30 m. Therefore, the frequency content of the records had to be low enough to obtain the phase velocities at large wavelengths as well. The geophone geometry was chosen according to the local situation. The corresponding dimensions of the arrays were between 42 and 75 m in length, with regular and irregular geophone inter-station distances (varying from 1 to 5 m), depending upon the available space. This information is summarized in Table 2 for all the 30 sites. The ambient noise was sampled during the day time for approximately 15 min at a 512-Hz sampling rate. We avoided performing measurements on artificial soil (e.g., asphalt, concrete and pavement), while stiff artificial soils may severely affect the HVSR curves (Castellaro and Mulargia, 2009b). The ReMi analysis was then performed using Soil Spy Rosina and Grilla software (Micromed, 2006; Micromed, 2008b). For both the ReMi and ESAC method the record was divided into 10-second, non-overlapping time windows and dispersion curves for each were computed. In the case of the ReMi non-informative time windows were removed and the final dispersion curve is the average of each time window dispersion curves. No time windows were removed in the case of the ESAC method and the final dispersion curve is also Table 2. Data-acquisition parameters, fundamental soil frequency f0 and Vs30 results Point Coordinate X Coordinate Y Array L-shape geometry * Geophone inter-station distance [m] HVSR f0 [Hz] s Vs30 ** [m/s] SS 01 458976 99345 16 + 8 varying from 1 to 5 2.7 0.08 285 SS 02 460625 98439 16 + 8 3 1.2 0.01 225 SS 03 463246 97829 16 + 8 3 1.1 0.13 137 SS 04 461607 100062 16 + 8 3 2.3 0.08 266 SS 05 461943 101264 16 + 8 3 3.5 0.07 396 SS 06 463986 102671 16 + 8 3 2.6 0.11 530 SS 07 457571 101120 16 + 8 varying from 1 to 5 3.7 0.01 347 SS 08 458979 102780 16 + 8 varying from 1 to 5 2.9 0.05 357 SS 09 459227 103263 16 + 8 3 3.2 0.08 350 SS 10 459413 102716 16 + 8 varying from 1 to 5 3.7 0.02 339 SS 11 460528 99546 16 + 8 3 2.3 0.01 256 SS 12 464198 103737 16 + 8 5 2.3 0.06 557 SS 13 462971 103350 16 + 8 3 2.2 0.05 561 SS 14 465082 97149 16 + 8 3 2.1 0.02 235 SS 15 464495 98035 16 + 8 3 3.7 0.03 224 SS 16 462369 98766 16 + 8 3 2.2 0.10 149 SS 17 459481 100711 16 + 8 3 3.1 0.04 302 SS 18 460723 100853 16 + 8 3 3.8 0.07 232 SS 19 463018 100758 17 + 7 varying from 2 to 3 3.6 0.01 332 SS 20 463185 101590 14 + 5 varying from 2 to 3 2.5 0.26 585 SS 21 462576 102243 18 + 6 varying from 2 to 3 2.6 0.09 595 SS 22 460860 103347 17 + 7 3 2.1 0.16 515 SS 23 461579 104405 18 + 6 3 2.1 0.03 605 SS 24 465818 102771 18 + 6 3 3.1 0.06 545 SS 25 467418 102529 18 + 6 3 3.8 0.13 549 SS 26 466896 101553 18 + 6 3 3.3 0.05 465 SS 27 465641 101549 18 + 6 3 2.9 0.15 432 SS 28 465969 99974 18 + 6 3 4.4 0.16 393 SS 29 467257 100051 18 + 6 3 3.3 0.11 456 SS 30 458712 101060 16 + 8 3 2.7 0.01 336 * The number of geophones in the longer and in the shorter leg of the L-shaped array ** Typical uncertainty in Vs30 determination is of about 10-20% (Mulargia and Castellaro, 2009) obtained as an average of all the time-window dispersion curves. An iterative grid-search procedure was used in both methods in order to find the values of the phase velocity that gives the best fit to the data. During the deployment of the seismic array, singlestation ambient microtremor noise was recorded as well at three or four locations along the array for 15 minutes at a 128-Hz sampling rate with a three-component, high-sensitivity Tromino tromograph (Micromed, 2005) developed for high-resolution digital measurements of the seismic noise. The data analysis to obtain the HVSR was performed using Grilla software (Micromed, 2006). The recorded noise was divided into 45 non-overlapping time windows, each being 20 seconds long. The signal of each window was corrected for the base line, padded with zeroes and tapered with a Bartlett window. The relevant spectrogram was smoothed through a triangular window with a frequency-dependent half-width of 5% of the central frequency. For each time window both horizontal spectra (NS and EW direction) were geometrically averaged and then divided by the vertical spectra (Equation 2) to compute the HVSR curve. When transient signals were observed, the relevant portions of the recordings were rejected. This procedure influences only the reliability of the confidence interval and does not bias the average HVSR (D'Amico et al., 2004). The final HVSR curve with the relevant 95% confidence interval is obtained by averaging all the time-window HVSR curves. The underground shear-wave velocity structure was obtained through joint modelling of the dispersion curve and the HVSR curve by using the Phase Velocity Spectra Module in the Grilla software (Micromed, 2008b and Micromed, 2007) that assumes a vertically heterogeneous, one-dimensional elastic model. This assumption was verified by taking several HVSR recordings along the array and comparing them. Some examples are shown in Figure 5, i.e., measurements along SS-02 (panel a), SS-16 (panel b), SS-19 (panel c) and SS-24 (panel d). All four cases show essentially identical HVSR curves and thus the assumption of 1D appears to be very well satisfied. The only difference is observed in the case of SS-24 in the frequency range 15 to 40 Hz (Figure 5 (d)), which corresponds to a depth of 3 to 5 m. We estimate that this difference is of minor importance. The visual comparison between the ReMi and ESAC final dispersion curves was made to check their coher- ency. Finally, an interactive modelling of the final dispersion curve and the HVSR curve with a trial-and-error procedure to obtain a shear-wave velocity versus depth profile was performed. Each layer in a one-dimensional ground model is defined by its thickness (h), shear-wave velocity (Vs), longitudinal wave velocity (Vp) and mass densities (y). It should be noted that Vp and the density p are of minor importance in determining the shape of the dispersion curve. The density p values have relatively small variations with depth, ranging from 1.6 to 2.M03 kg/m3. Moreover, the Vp value depends on Vs through the Poisson modulus v (Equation 6), which was fixed to 0.35 during the process of modelling. In the case of this study the constraint is provided by the Vs of the first layer, inferred from the ReMi and ESAC methods applied at the same site. V = I^ (6) Vp V 1-n () In the starting models we introduced, at first, a solution of a maximum four layers. Then we estimate the thicknesses and velocities for the top two layers according to the HVSR and dispersion curves, fixing those values and then modelling the deeper layers. Once we had a satisfactory fit with the HVSR curve and the dispersion curves, we fixed the geometry of the layers and varied the layer velocities to estimate the best fit. During the modelling we fitted a synthetic model, not only to the main resonance frequency of the HVSR curve, which is in general the largest peak, but also to other stable humps and troughs in the HVSR curve. KVSR companion ilonq !h* f'lf SS-02 HVSR i dmpflmim Vonq tht itray SS-19 HVSR companion *tong1h* arjaj SS-T6 »buwwicy (US MvSfl iwnpBitHjn afong the aita* SS-21 Figure 5. Comparison of the HVSR curve recorded at different places along the array; thicker lines are the average HVSR and the thinner lines represent the 95% confidence interval Finally, the average shear wave velocity values in the upper 30 m (Vs30) were calculated in accordance with the following expression (CEN, 2004): V= 30 Eh V i=1,N vi (7) where hi and Vi denote the thickness (in metres) and the shear-wave velocity (at a shear strain level of 10-5 or less) of the i-th layer, in a total of N, existing in the top 30 m. In the end the Vs30 data were used for the preparation of the shear-wave velocity distribution map of the study area in the topmost 30 m of sediments. 6 RESULTS AND DISCUSSION The results of the microtremor measurements carried out in Ljubljana city include the HVSR curves, the Rayleigh wave-dispersion curves and the corresponding shear-wave velocity profiles. In Table 2 are the fundamental resonant soil frequency and Vs30 values computed for each site. Some characteristic examples are shown in Figures 6 to 11; panel (a) shows the experimental HVSR curves (in a red colour, the black lines present the 95% confidence interval) and the fitted synthetic HVSR curves in the blue coloured solid line. The comparison between the experimental dispersion curves from the ReMi and ESAC methods with the theoretical one is shown in panel (b). The dispersion curves derived from the ESAC are shown with yellow circles and the ReMi dispersion curves are shown as contour plots. The light blue circles indicate the theoretical Rayleigh wave phase dispersion curves obtained by joint modelling and corresponding to the shear-wave velocity profile shown in panel (c). The main purpose of the HVSR measurements was not to estimate the soil's fundamental frequency of resonance but to operate a joint fit of the HVSR and dispersion curves and to verify the 1D soil assumption underlying array techniques. However, the values are in good agreement with results from Gosar et al. (2010). Clear peaks were obtained in the whole southern part of the city (i.e., Figures 3 and 9) and also at some sites in the northern part (i.e., Figure 8). Generally, the northern part of the city is characterized by HVSR curves showing two or more peaks with modest amplitude (Figures 7, 10, and 11), denoting subsoil locally characterized by two or more significant impedance contrasts at different depths. A second peak occurs at frequencies above 10 Hz, suggesting a very shallow seismic interface. This ascertainment corresponds to known information about the complex geology of stiffer cemented gravel horizons inside the sedimentary cover in this area. For most sites, the velocity of the Rayleigh dispersion increases continuously as the frequency decreases, indicating a continuously increasing shear-wave velocity with depth. The phase velocity estimates from the ReMi and ESAC methods are very consistent in the high-frequency range, i.e., from 15 to 35 Hz. In the frequency range below 5 Hz both methods were less effective due to the limited geometry of field lay-out and the limitations of the 4.5-Hz geophones. Furthermore, the depth penetration of an array is more often limited by the presence of strong impedance contrast at some depth rather than by the array geometry. In the case of the presence of a strong impedance contrast the most part of waves would be reflected and keep travelling in the shallow layer only, and thus not provide any information about the structure below that shallow layer. Usually, in the absence of strong impedance contrast it is feasible to inspect depths of 4 times the array length. The results from the ReMi method were severely limited below 10 Hz where the dispersion curves become less coherent, but generally allowed good coverage of the phase velocity in the high-frequency range up to 50 Hz (evident from Figures 6 to 10), except in the case SS-27 (Figure 11). The low-frequency limit of the dispersion curve from the ESAC is in agreement with that from the HVSR curve, as at the resonant frequency the vertical component of the Rayleigh wave vanishes and therefore the phase-velocity values are reliable only from that frequency onwards (Scherbaum et al., 2003). On the other hand, in some cases the ESAC method was unable to provide clear information above 20 Hz (e.g., SS-11, SS-12 and SS-27). From the ReMi contour plot for site SS-03 the fundamental mode of the dispersion curve, as well as two higher modes, is observed (Figure 6). Only the fundamental mode is visible in the case of SS-06, SS-10 and SS-11 (Figures 7 to 9). The fundamental modes from the ReMi method and the resulting curves from the ESAC method, which only shows the predominant one, agree well. Hence it follows that in most cases the higher mode does not severely mask the fundamental mode in the ESAC curves. There were only a few exceptions, where the higher mode in a certain frequency range is dominant over the fundamental one. The Rayleigh dispersion result of the ESAC method at site SS-12 shows that the dispersion in the frequency range 20-40 Hz corresponds to a higher mode, which is evident from the comparison with the dispersion curves of the ReMi contour plot (Figure 10). Joint modelling with the HVSR data has proved very useful to constrain the model parameters and reduce the ambiguity, which cannot be avoided if only the dispersion curves are used to derive a velocity profile. SS-03 Synthetic vs experimental HVSR Vs profile Figure 6. Site SS-03 (a) The experimental HVSR curve in red (black lines present the 95% confidence interval) and the fitted synthetic HVSR curve in blue (b) The comparison between experimental dispersion curves from the ReMi (contour plot) and ESAC (yellow circles) methods with the theoretical one (blue circles) (c) Corresponding shear-wave velocity profile SS-06 Synthetic vs experimental HVSR Mix HVatiib to 11 M: (ullifl Jjnqr 0 0 iiOKSl Dispersion curves ■50 ■60 10 IS 20 35 30 Frequency [Hz] JO 45 50 300 Vs profile 600 900 J200 1500 (c) Vs30 = 530 m/s Vs [m/s] SS-10 Synthetic vs experimental HVSR Vs profile Figure 8. As in Figure 6, but for site SS-10 SS-11 Synthetic vs experimental HVSR Vs profile »¡.aa-stom 100 2M 300 400 500 600 700 SS-1 2 Synthetic vs experimental HVSR Vs profile Figure 10. As in Figure 6, but for site SS-12 SS-27 Synthetic vs experimental HVSR The calculated shear-wave velocity distribution down to a depth of 30 m in Ljubljana city is presented in Table 2 and in Figure 12. Black circles on the figure mark the location of the 30 measuring points. The Vs30 distribution map is overlaid over the existing microzonation based on the EC8 (Zupančič et al., 2004). On the basis of Vs30 values the study area can be classified as ground type B (360-800 m/s), C (180-360 m/s), D and E (< 180 m/s). The northern part of the study area, located in the Ljubljana Field, has higher shear-wave velocity values, within the range 360-600 m/s, and thus falls into the ground type B. The low-velocity sites are in general related to Ljubljana Moor and correspond to the ground types C, D and partly E. According to our study, the sites have been generally classified into higher soil classes, i.e., having better seismological and geotechnical properties than indicated by previous study. A comparison of the geological map (Figure 2) and the Vs30 map shows very good agreement, reflecting the fact that low velocities in the southern part correspond to soft lacustrine and fluvial Quaternary sediments. The Vs30 values are then increasing towards the north, where stiffer material (Pleistocene and Holocene glaciofluvial deposits) is dominant. Also, a comparison with the sediment resonant frequency map (Figure 3) clearly highlights the relations between the fundamental soil frequency, the shear-wave velocity and the thickness (Equation 3) known from sparse geotechnical data. 7 CONCLUSIONS Although no damaging earthquakes have occurred since 1963 in the vicinity of Ljubljana, this area is one of the most seismically active regions in Slovenia. For this reason, microzonation studies are very important for earthquake engineering design. The Eurocode 8, which constitutes the standard for seismic-resistant design in Figure 12. Vs30 distribution map (coloured contour plot) and previous soil classification for Ljubljana (after Zupančič et al., 2004) according to EC8 ground types (b & w hatching plot) Europe and also in Slovenia, adopts as a key quantitative parameter, the average shear-wave velocity in the first 30 m of subsoil, commonly termed as Vs30. In this paper the results of the shear-wave velocity determination in the top 30 m (Vs30) based on joint modelling of the HVSR, ReMi and ESAC microtremor data are presented. The ReMi and ESAC methods provided very similar phase-velocity dispersion curves at all 30 measuring sites in Ljubljana city. The combination of both methods made it possible for the dispersion curves of the fundamental mode Rayleigh wave to be determined for the frequencies from approximately 2 Hz and up to 50 Hz. The dispersion curves were further used to develop one-dimensional, shear-wave velocity profiles from which the Vs30 values were calculated. Simultaneous HVSR modelling has proved to be very useful to constrain model parameters and reduce ambiguity, which cannot be avoided if only dispersion curves are used to derive a velocity profile. Based on the Vs30 values the study area can be classified as ground types B, C, D and E. In general, the study area could be divided into two parts with different average velocities. In the northern part, higher velocities were obtained, ranging from 360-600 m/s, and thus correspond to ground type B. Lower velocities are characteristic for the whole southern part of Ljubljana. The values there are below 300 m/s with the lowest shear-wave velocity of 137 m/s and the area is therefore classified as ground types C, D and E. The results suggest an interesting difference in comparison with the previous microzonation based on EC8 (Zupančič et al., 2004). The velocities obtained in this study indicate that soils have better seismogeological conditions and can be classified in general to one ground type class better according to EC8 that in the previous study. It should be stressed that the previous microzonation did not involve any direct shear-wave velocity measurements. Although Vs30 is perhaps not the most suitable parameter to define seismic site response, it is presently requested by many national seismic regulations. Nevertheless, the spatial distribution of Vs30 has provided valuable information to integrate and supplement existing seismic microzonations of Ljubljana. The true progress of this study is the acquisition of several Vs profiles that are much more informative than Vs30 alone and can be used for direct, one-dimensional modelling of the seismic ground motion. ACKNOWLEDGMENTS The authors are indebted to Tjaša Mlinar for her help with the field measurements. Thanks are due to Polona Zupančič for the compilation of the seismicity and geological map in Figures 1 and 2. The study was realized with the support of research program P1-0011 and grant 1000-05-310228 financed by the Slovenian Research Agency. REFERENCES Aki, K. (1957). Space and time spectra of stationary stochastic waves, with special reference to micro-tremors. Bulletin of the Earthquake Research Institute 35: 415-456. Aki, K. and Richards, P.G. (2002). Quantitative Seismology, 2nd Edition. University Science Books, 700 pp. Anderson, J.G. (2007). Physical Processes That Control Strong Ground Motion. In Hiroo Kanamori, Volume Editor; Gerald Schubert, Editor in Chief (Ed.), Treatise on Geophysics, Vol. 4, Earthquake Seismology (Vol. 4, 513-565). Amsterdam: Elsevier. Asten, M.W. (2001). The spatial auto-correlation method for phase velocity of microseisms - another method for characterisation of sedimentary overburden: in Earthquake Codes in the Real World, Australian Earthquake Engineering Society, Proceedings of the 2001 Conference, Canberra. Bard, P.Y. (1999). Microtremor measurements: a tool for site effect estimation? In: Irikura, K., Kudo, K., Okada, H., Sasatami, T. (Editors) The effects of surface geology on seismic motion. Balkema, Rotterdam, pp 1251-1279. Beaty, K.S. and Schmitt, D.R. (2003). Repeatability of multimode Rayleigh-wave dispersion studies. Geophysics 68: 782-90. Castellaro, S., Mulargia, F., Rossi, P.L. (2008). Vs30: Proxy for Seismic Amplification? Seismological Research Letters, Vol. 79, No. 4: 540-543. Castellaro, S., Mulargia, F. (2009a). Vs30 estimates using constrained H/V measurements. Bulletin of the Seismological Society of America 99: 761-773. Castellaro, S., Mulargia, F. (2009b). The effect of velocity inversions on H/V. Pure and Applied Geophysics 166: 567-592. CEN (2004). Eurocode 8—design of structures for earthquake resistance. Part 1: general rules, seismic actions and rules for buildings. European standard EN 1998-1, December 2004, European Committee for Standardization, Brussels D'Amico, V., Picozzi, M., Albarello, D., Naso, G., Tropenscovino, S. (2004). Quick estimates of soft sediment thicknesses from ambient noise horizontal to vertical spectral ratios: a case study in southern Italy. Journal of Earthquake Engineering 8: 895-908. Gosar, A., Rošer, J., Šket-Motnikar, B., Zupančič, P. (2010). Microtremor study of site effects and soil-structure resonance in the city of Ljubljana (central Slovenia). Bulletin of Earthquake Engineering 8: 571-592. Grad, K. and Ferjančič, L. (1974). Basic geological map of Yugoslavia 1:100.000 - sheet Kranj. Geological Survey of Slovenia, Ljubljana. Lapajne, J. (1970). Seismic microzonation of Ljubljana - geophysical investigations. Geological Survey of Slovenia, Ljubljana, 16 pp. Lapajne, J., Šket-Motnikar, B., Zupančič, P. (2001). Design ground acceleration map of Slovenia. Potresi v letu 1999: 40-49. Louie, J. (2001). Faster, better: shear-wave velocity to 100 meters depth from refraction microtremor arrays. Bulletin of the Seismological Society of America 91: 347-364. Medvedev, S.V. (1965). Inženjerska seizmologija. Grad-jevinska knjiga, Belgrade, 271 pp. Mencej, Z. (1989). The gravel fill beneath the lacustrine sediments of the Ljubljansko barje. Geologija 31-32: 517-553. Micromed (2005). Tromino, portable seismic noise acquisition system, user's manual. Micromed, Treviso, 101 pp. Micromed (2006). Grilla ver. 2.2, spectral and HVSR analysis - user's manual. Micromed, Treviso, 47 pp. Micromed (2007). An introduction to the H/V inversion for stratigraphic purposes. Micromed, Treviso, 31 pp. Micromed (2008a). SoilSpy Rosina ver. 2.0 - user's manual. Micromed, Treviso, 81 pp. Micromed (2008b). An introduction to the phase velocity spectra module in Grilla. Micromed, Treviso, 16 pp. Mucciarelli, M. and Gallipoli, M.R. (2001). A critical review of 10 years of microtremor HVSR technique. Bolletino di Geofisica teorica ed applicata 42: 255-266. Mulargia, F. and Castellaro, S. (2009). Experimental Uncertainty on the Vs(z) Profile and Seismic Soil Classification. Seismological Research Letters, Vol. 80, No. 6: 985-988. Nakamura, Y. (1989). A method for dynamic characteristics estimation of subsurface using microtremor on the ground surface. Quarterly Report of Railway Technical Research Institute (RTRI) 30: 25-33. Nakamura, Y. (2000). Clear identification of fundamental idea of Nakamura's technique and its applications. 12WCEE, Auckland. Nogoshi, M. and Igarashi, T. (1970). On the propagation characteristics estimations of subsurface using microtremors on the ground surface. Jour. Seismol. Soc. Japan 23, 264-280. Nogoshi, M. and Igarashi, T. (1971). On the amplitude characteristics of microtremor (Part 2) (in Japanese with English abstract). Jour. Seism. Japan 24, 26-40. Ohori, M., Nobata, A., Wakamatsu, K. (2002). A comparison of ESAC and FK methods of estimating phase velocity using arbitrarily shaped microtremor arrays. Bulletin of the Seismological Society of America 92, 2323-2332. Okada, H. (2003). The microtremor survey method. American Geophysical Monograph 12, Society of Exploration Geophysicist, 135 pp. Premru, U. (1982). Basic geological map of Yugoslavia 1:100.000 - sheet Ljubljana. Geological Survey of Slovenia, Ljubljana. Scherbaum, F., Hinzen, K.G., Ohrnberger, M. (2003). Determination of shallow shear wave velocity profiles in Cologne Germany area using ambient vibrations. Geophysical Journal International 152: 597-612. SIST EN 1998-1:2005/oA101 (2005). Eurocode 8, design of structures for earthquake resistance - part 1: general rules, seismic actions and rules for buildings, national Annex. Slovenian institute for standardization, Ljubljana. Thorson, J.R. and Claerbout, J.F. (1985). Velocity-stack and slant-stack stochastic inversion. Geophysics 50: 2727-41. Žlebnik, L. (1971). Pleistocen Kranjskega, Sorškega in Ljubljanskega polja (Pleistocene of Kranj, Sora and Ljubljana Field). Geologija 14, 5-52. Zupančič, P., Šket-Motnikar, B., Gosar, A., Prosen, T. (2004). Seismic microzonation map of the municipality of Ljubljana. Potresi v letu 2002: 32-54.