fi. s. Nossan ADVANCES AND UNCERTAINTIES IN THE DESIGN OF ANCHOAED RETAINING MALLS USING NUMEAICAL MODELLING s. škrabl. THE LIMIT VALUES AND THE DISTAIBUTION OF THAEE-DIMENSIONAL PASSIVE EAATH PRESSURES s. Lenart THE RESPONSE OF SATUAATED SOILS TO A DYNAMIC LOAD t. pllberšek & n. umek GAEEN'S FUNCTION FOA TANGENTIALY LOADED HOAIZONTALY LAYEAED HALF-SPACE r— H O I fiCTfl GCOTCCHNICñ SLOVGNICñ ISSN: 1854-0171 ustanovLteLjL Founders urednLškL odbor edLtorLaL Board Univerza v Mariboru, Fakulteta za gradbeništvo University of Maribor, Faculty of Civil Engineering š S / 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 LzdajateLj IpubLLsher Univerza v Mariboru, Fakulteta za gradbeništvo University of Maribor, Faculty of Civil Engineering odgovornL urednLk edLtor-Ln-chLef Ludvik Trauner Univerza v Mariboru urednLka co-edLtors Stanislav Škrabl Univerza v Mariboru Bojan Žlender Univerza v Mariboru TehnLčna urednLka Desk edLtors Bojana Dolinar Univerza v Mariboru Borut Macuh Univerza v Mariboru Lektor proof-Reader Paul McGuiness NakLada cLrčuLatLon 500 izvodov - issues TLsk prLnt Tercia tisk d.o.o. Ptuj Revija redno izhaja dvakrat letno. Članki v reviji so recen-zirani s strani priznanih mednarodnih strokovnjakov. Baze podatkov v katerih je revija indeksirana: SCIE - Science Citation Index Expanded, JCR - Journal Citation Reports / Science Edition, ICONDA - The international Construction database, GeoRef Pri financiranju revije sodeluje Javna agencija za raziskovalno dejavnost republike Slovenije. Darinka Battelino Università degli Studi di Trieste József Farkas Budapesti Müszaki és Gazdaságtudományi Egyetem Theodoros Hatzigogos Aristotle University of Thessaloniki Rolf Katzenbach Technische Universität Darmstadt Zlatko Langof Univerzitet u Sarajevu Jakob Likar Univerza v Ljubljani Janko Logar Univerza v Ljubljani Bojan Majes Univerza v Ljubljani Milan Maksimovic Univerzitet u Beogradu Borut Petkovšek Zavod za gradbeništvo Slovenije Mihael Ribičič Univerza v Ljubljani César Sagaseta Universidad de Cantabria Stephan Semprich Technische Universität Graz Abdul-Hamid Soubra Université de Nantes Ivan Vaniček Ceské vysoké učeni technické v Praze Franjo Veric Sveučilište u Zagrebu Address NasLov urednLštva ACTA GEOTECHNICA SLOVENICA Univerza v Mariboru, Fakulteta za gradbeništvo Smetanova ulica 17 2000 Maribor Slovenija Telefon / Telephone: +386 (0)2 22 94 300 Faks / Fax: +386 (0)2 25 24 179 E-pošta / E-mail: ags@uni-mb.si spLetnL nasLov http://www. fg.uni-mb.si/journal-ags I web Address The journal is published twice a year. Papers are peer reviewed by renowned international experts. Indexation data bases of the journal: SCIE - Science Citation Index Expanded, JCR - Journal Citation Reports / Science Edition, ICONDA - The international Construction database, GeoRef Financially supported also by Slovenian Research Agency. vsebina fr !v Ludvik Trauner uvodnik nntun szavlts Nossan napredek in nezanesljivosti pri numeričnem modeliranju sidrnih konstrukcij Stanislav Škrabl mejne vrednosti in razporeditev 3d pasivnih zemeljskih pritiskov s Stanislav Lenart odziv zasičenih zemljin na dinamično obremenitev [50j Tomaž Pliberšek Ln Andrej umek greenova funkcija tangencialno obremenjenega horizontalno slojevitega polprostora s3 navodila avtorjem contents Ludvik Trauner editorial a Rntun szavits Nossan advances and uncertainties in the design of anchored retaining walls using numerical modelling e stanisLav Škrabl the limit values and the distribution of three-dimensional passive earth pressure StanisLav Lenart the response of saturated soils to a dynamic load Tomaž Pliberšek in nndrej umek green's function for tangentialy loaded horizontaly layered half-space e instructions for authors e uvodnik Člani uredniškega odbora revije Acta Geotechnica Slovenica z veseljem sporočamo, da smo v petem letu izhajanja revije dosegli pomemben cilj, to je njeno vključitev v Thomsonovi bazi Science Citation Index Expanded in Journal Citation Reports/Science Edition, kjer sedaj poteka preverjanje člankov za določitev faktorja vpliva. Ob tej priliki se zahvaljujemo vsem, ki ste nas podpirali in vseskozi trdno verjeli v uspeh revije. Druga novost je odločitev uredniškega odbora revije in vodstva Slovenskega geotehniškega društva (SloGeD), da bodo odslej teme vabljenih predavanj na Šukljetovih dnevih predstavljene v obliki člankov v reviji Acta Geotechnica Slovenica. Šukljetove dneve vsako leto organizira SloGeD v spomin na pionirja slovenske geotehnike, akademika prof. dr. Luja Šukljeta. V duhu tega velikega raziskovalca so izbrane tudi vsebine vabljenih predavanj. Doslej so bili na Šukljetove dneve povabljeni naslednji predavatelji iz tujine: M. Jamiolkowski (2000), G. Sanglerat (2001), D.D.Potts (2002), H. Brandl (2003), R. Katzenbach (2004), S. Semprich (2005), S. Leroueil (2006), D. Žnidarčič (2007) in A. Szavits Nossan (2008). Pričujoča številka revije prinaša štiri zanimive prispevke avtorjev Antuna Szavits Nossana, Stanislava Škrabla, Stanislava Lenarta ter Tomaža Pliberška in Andreja Umeka. V prispevku A. Savitz Nossana je predstavljena možnost napovedovanja horizontalnih premikov in notranjih statičnih količin v sidranih podpornih konstrukcijah za zaščito izkopov z uporabo standardnih terenskih in laboratorijskih preiskav ter komercialnega programa s končnimi elementi. Slednji vsebuje konstitutivni model zemljine, ki simulira osnovne aspekte obnašanja tal na lokaciji gradbene jame. Članek S. Škrabla obravnava izvirni pristop določanja kritične razporeditve in mejnih vrednosti pasivnih zemeljskih tlakov za tri-dimenzionalne primere po metodi mejne analize in teorema zgornje vrednosti. Za določanje kritične razporeditve pasivnih tlakov vzdolž višine podporne konstrukcije je uporabljena metoda mejne analize z množico tri-dimen-zionalnih kinematično dopustnih hiperboličnih rotacijskih porušnih mehanizmov po metodi postopnega določanja intenzitete pasivnih tlakov od zgoraj navzdol. Avtor tretjega prispevka S. Lenart predstavlja dva najbolj izrazita načina deformiranja dinamično obremenjenih zasičenih zemljin in sicer likvifakcijo s tečenjem in ciklično mobilnost. Oba pojava sta bila preiskana na meljnih peskih in prodno peščenih meljih, ki izhajajo z območja potopljenega nasipa železniške proge zaradi novozgrajenega akumulacijskega bazena na reki Savi v Boštanju in velikega plazu, ki se je sprožil na področju Stože v Julijskih Alpah. Članek T. Pliberška in A. Umeka obravnava nov pristop k evaluaciji integralne predstavitve Greenove funkcije za slojevit pol-prostor, ki je na površini obremenjen s harmonično tangencialno točkovno silo. Ludvik Trauner I .I o i m i 11 rc>r\ in IA editorial I am very pleased to be able to report that Acta Geotechnica Slovenica has been selected for coverage in Thomson Reuters products and custom information services. Beginning with Vol.4 (1) 2007, Acta Geotechnica Slovenica will be indexed and abstracted in the Science Citation Index Expanded (also known as SciSearch®) and the Journal Citation Reports/Science Edition. I would like to take this opportunity to thank all of you who have supported us and contributed to the continued success of our journal. Other news is the decision of journal's editorial board and the Slovenian Geotechnical Society (SloGeD) that the invited lectures of Suklje's Days, which are organized every year by SloGeD in memory of the pioneer of Slovene geotechnics, academician professor Lujo Šuklje, will be published in Acta Geotechnica Slovenica. In this way the journal will have the opportunity to be the first to publish the results of these important scientific studies. The publication of the contents of the invited lectures in the form of papers is a great advantage for the authors also, who will now have the opportunity to acquaint a broader circle of interested readers with their achievements, since these articles can be found in the Thomson databases. In recent years, the invited lecturers on Suklje's Days were M. Jamiolkowski (2000), G. Sanglerat (2001), D.D. Potts (2002), H. Brandl (2003), R. Katzenbach (2004), S. Semprich (2005), S. Leroueil (2006), D. Žnidarčic (2007) and A. Szavits Nossan (2008). The first issue of Year 5 contains four interesting articles authored by Antun Szavits Nossan, Stanislav Škrabl, Stanislav Lenart, and Tomaž Pliberšek and A. Umek. In his paper, A. Szavits Nossan presents a prediction of the horizontal displacements and the internal forces in an anchored wall for the protection of an excavation using standard field and laboratory tests and a finite-element program. The last of these includes a constitutive soil model that can simulate the key aspects of the soil's behaviour at a construc-# tion site. The paper of S. Škrabl deals with a novel approach to the determination of the critical distribution and limit values of three-dimensional passive soil pressures acting on flexible walls, following the upper-bound method within the framework of the limit-analysis theory. The method of limit analysis with a set of three-dimensional cinematically admissible hyperbolic translational failure mechanisms is used to determine the critical distribution of the passive pressures along the retaining structure's height. The author of the third paper, S. Lenart, treats the two most marked types of deformation behaviour for dynamically loaded saturated soil, i.e., flow liquefaction and cyclic mobility. Both phenomena were researched in silty sand and lacustrine carbonate silt, which are found in the area of a submerged railway line, due to the newly built Sava-river accumulation reservoir in Boštanj, and where the large landslide occurred in the Stože area of the Julian Alps. The topic of the paper by T. Pliberšek and A. Umek is a novel evaluation of the integral representation of a surface Green's function for a layered half-space, loaded on its surface by a harmonic tangential point force. Ludvik Trauner napredek in nezanesljivosti pri nume-ričnem modeliranju sidrrnh podpornih konstrukcij_ ANTUN SZAVITS NOSSAN o avtorju Antun Szavits Nossan Univerza v Zagrebu, Fakulteta za gradbeništvo Kačičeva 26, 10 000 Zagreb, Hrvaška E-pošta: szavits@grad.hr Izvleček V članku je predstavljena možnost napovedovanja horizontalnih premikov in notranjih statičnih količin v sidranih podpornih konstrukcijah za zaščito izkopov z uporabo standardnih terenskih in laboratorijskih preiskav ter komercialnega programa s končnimi elementi. Slednji vsebuje konstitutivni model zemljine, ki simulira osnovne aspekte obnašanja tal na lokaciji gradbene jame. V prispevku je prikazano, da se mora uporabnik dobro seznaniti s konstitutivnim modelom vključenim v program ter da predstavlja odločilen del modeliranja izbira primernih parametrov zemljin za numerične analize. Za izbiro primernih parametrov za simulacijo dejanskih pogojev prisotnih med gradnjo je koristna izvedba numeričnih simulacij standardnih laboratorijskih preizkusov, ki jih je potrebno primerjati s poznanim obnašanjem zemljine. V članku je prikazano, da izmerjene hitrosti strižnih valov, iz katerih lahko določimo strižno togost tal pri majhnih deformacijah, lahko uporabimo tudi za določitev statične togosti tal za velikosti deformacij obravnavane geotehnične konstrukcije, tako v koherentnih kot tudi nekoherentnih tleh. Raziskovalno delo je bilo izvedeno za primer iz geotehnične prakse z detajlno analizo zaščite izkopa s sidrano armirano betonsko steno v relativno togi zemljini. Deformacije stene so bile merjene z vgrajenim inklinometrom. Pretežni del članka predstavlja izbiro parametrov konstitutivnega modela, še posebej parametrov togosti tal. Za potrditev ocene zmanjšanja sekantnega deformacijskega modula zaradi povečanja mobilizirane strižne trdnosti za trde gline z objavljenimi empiričnimi odnosi iz literature je bila uporabljena simulacija triosnega konsolidacijskega nedreniranega preizkusa. Prikazano je, da je s takšno izbiro parametrov togosti v konstitutivnem modelu tal možno dobiti sprejemljivo napoved deformacij sidrane stene. Čeprav je predstavljen samo en primer uspešne analize, le ta daje vzpodbudo, saj prikazuje možnost relativno zanesljive napovedi deformacij samo na osnovi terenskih in laboratorijskih preizkusov in z uporabo razpoložljivih računalniških programov z realnim modelom zemljine. Ključne besede sidrana stena, model tal, strižna togost, numerično modeliranja, merjene deformacije advances and uncertainties in the design of anchored retaining walls using numerical modelling ANTUN SZAVITS NOSSAN About the author Antun Szavits Nossan University of Zagreb, Faculty of Civil Engineering Kaciceva 26, 10 000 Zagreb, Croatia E-mail: szavits@grad.hr Abstract This paper describes research on the prediction of horizontal displacements and internal forces in an anchored wall for the protection of an excavation, using standard field and laboratory tests and a finite-element programme with a soil model that can simulate the key aspects of soil behaviour at a construction site. It is important to be acquainted with the constitutive model incorporated in the programme, and the selection of the appropriate soil parameters for the numerical analysis is a crucial part of the modelling. As a result, it is useful to carry out numerical simulations of standard laboratory tests with well-known soil behaviour in order to select the relevant parameters for the simulation of the actual construction process. It is shown in this paper that the measurements of the shear-wave velocities, which can provide the soil's stiffness at very small strains, can also be useful for determining the static stiffness at a magnitude of the strains relevant for the geotechnical structure under consideration, for both cohesive and noncohesive soils. The research was carried out by a detailed analysis of a case history involving an anchored, reinforced concrete wall supporting the walls of an excavation in a relatively stiff soil. The wall displacements were monitored using an installed inclinometer. The major part of the paper is devoted to an analysis of the selection of parameters, especially the stiffness parameters. The simulation of the triaxial, consolidated, undrained tests was used in order to assess the reduction of the secant stiffness modulus with an increase of the relative mobilized shear strength for the hard clay layer according to the published empirical evidence. It is shown that by selecting the appropriate stiffness parameters for the soil model used in the numerical analysis, it is possible to get an acceptable prediction of the anchored-wall displacements. This is just one example of a successful analysis, but it is encouraging in the way that it shows how it is possible to make reliable predictions based on standard field and laboratory tests and with the use of an available computer programme with a realistic soil model. Keywords anchored wall, soil model, shear stiffness, numerical modelling, measured displacements 1 INTRODUCTION Anchored, retained structures are often used as temporary protection for deep excavations in urban areas. Their role is to ensure the stability of the soil around the excavation and to prevent any damage to surrounding buildings that might be caused by the excavation. The successful design of such structures depends a great deal on a realistic solution to the interaction between the structure, the anchors and the soil, taking into consideration the mechanical characteristics of the surrounding soil as well as the manner and the sequence of the construction. Gaba et al. [1] gave an overview of the available numerical methods, together with an assessment of their advantages and drawbacks. A detailed solution to the interaction problems is becoming increasingly more accessible with the use of commercial, numerical tools based primarily on the finite-element method, which allows for the use of complex, constitutive soil models [2], [3]. There are, however, serious problems with the practical use of these tools. Schweiger [4] describes a detailed benchmarking experiment in which several experts were invited to numerically model the behaviour of an anchored diaphragm wall. The results were scattered over an alarmingly wide range, which is not acceptable in practice, due to the selection of different constitutive models and soil parameters. De Vos and Whenham [5] have shown the results of a survey among a large number of users of geotechnical finite-element programmes that show the problems they were encountering. The first item on the list of problems is the determination of soil parameters (23% of answers), followed by the determination of the initial conditions in the soil, the selection of the constitutive soil model, the interpretation of the results, the numerical discretization, the boundary conditions and the selection of the type of analysis. The first three items represent the core of the geotechnical design, supported by numerical modelling, and so they appeared at the top of the list in more than 50% of the answers. Gaba et al. [1] state, among others, the following reasons for these problems: the inadequate constitutive models, where the simple ones are not realistic; the data on soil strength; the stiffness and initial stresses that are not of sufficient quality; the insufficient user experience with the particular programme; and the inadequate modelling of the undrained conditions in cohesive soils. They claim that "Ground movements cannot be predicted accurately. It is essential that optimum use is made of precedent in comparable conditions through the use of good-quality case-history data. Case-history-based empirical methods of prediction are to be preferred to the use of complex analyses, unless such analyses are first calibrated against reliable measurements of well-monitored comparable excavations and wall systems." In any case, finite-element analyses should be used with caution, but they remain the only tool in cases of unusual structures for which there is no comparable experience. Studies in which complex numerical models are calibrated against the monitoring data of a case history can be helpful in resolving the above-mentioned problems related to the use of commercial finite-element programmes for geotechnical structures. This paper describes such a case history and the subsequent numerical modelling. The case history comprises an excavation protected by an anchored, retaining structure, of which there are several examples constructed recently in Zagreb, Croatia. Standard geotechnical investigations of average quality were carried out along with measurements of the shear-wave velocities with respect to depth. It was intended to use these measurements for the prediction of the anchored-wall displacements, based on the significance of this aspect of soil behaviour, which has recently gained attention [6], related to the soil-structure interaction [7] and particularly to the interaction of the soil with the anchored walls [8]. Shear-wave velocities provide a direct in-situ measure of the soil stiffness without the necessity to retrieve undisturbed soil samples or use problematic correlations. The anchored-wall displacements were monitored during construction, and the excavation was successfully completed. Subsequent numerical analyses were carried out using the finite-element programme Plaxis V8 [9], which is widely used in Croatia. Its option of small strain was used in order to take advantage of the shear-wave velocity measurements and the resulting soil stiffness at very small strains. It was decided, for practical reasons, to use Plaxis V8, even though sophisticated analyses of anchored walls at small strains have been reported [10], but using a commercially unavailable programme. Designers in Croatia are familiar with the use of Plaxis for modelling anchored structures. Their predictions of displacements based on the standard recommendations for the selection of soil parameters usually turn out as a significant overestimation in comparison with the measured wall displacements. As a result they use a higher soil stiffness, based on the argument of available data on similar structures in similar soils. This type of reasoning, which is not based on serious studies, makes the use of complex finite-element calculations questionable, because they do not seem to have a significant advantage over, for example, the method of a beam resting on elasto-plastic springs, where the springs' characteristics are determined empirically from displacement measurements on similar anchored walls. 2 THE CASE HISTORY The excavation, 14.5 m deep, is located in a rapidly expanding commercial area in Zagreb, and is intended for the construction of underground storeys of a commercial building. An existing, old, brick house, sensitive to soil displacements, is located near to the excavation. A 17.5-m-high and 0.6-m-thick wall of reinforced concrete, embedded in the soil 4 m below the bottom of the excavation, provided protection. The wall was cast in place prior to the excavation works. Three rows of BBR 1860/1660 pre-stressed ground anchors were installed at a horizontal distance of 2.5 m in each row. The upper, first row anchors consist of 4 strands of high-strength steel, 0.6" in diameter. The second and third row anchors consist of 5 strands. Each anchor in the first two rows was pre-stressed to 500 kN, whereas the anchors in the third row were pre-stressed to 650 kN. Inclinometer measurements of the relative horizontal wall displacements were taken during the excavation works. The inclinometer tube was installed in the wall concrete along its whole height at the location of the brick house. The vertical excavation section with the wall and the neighbouring house is shown in Fig. 1. The ground surface at the location is horizontal and the underlying ground is horizontally layered. The surface layer is around 2 m thick and it consists of medium dense fill and clay underlain by a layer of poorly graduated medium dense gravel down to a depth of 14 m. Below this depth is a thick layer of hard, overcon-solidated clay. The geotechnical field investigation was carried out in several 30-m-deep boreholes. Disturbed and undisturbed samples were retrieved and SPT measurements were taken. The shear-wave velocities were measured in two boreholes using the down-hole method. The underground water level was determined in the gravel layer at 7 m below the ground surface. Standard classification tests were carried out in the laboratory on disturbed samples, and undisturbed clay samples were used for the triaxial, consolidated, undrained (CIU) and unconsolidated, undrained (UU) tests. The CIU tests were performed with pore-water pressure measurements in order to determine the effective shear-strength parameters. The undrained shear strength was determined in the UU tests and with the use of a pocket penetrometer. The undisturbed clay samples were also used for oedometer tests. The results of the field and laboratory tests are presented in Fig. 2. The SPT blow count N was corrected by the standard hammer impact energy of 60% and the normalized vertical effective stress, where pref = 100 kPa, according to Skempton [11] The full line in Fig. 2 represents the selected characteristic value of the design parameter (N1)60 according to Eurocode 7 [12]. The same characteristic value of Figure 2. Soil profile with the fines and sand content in the gravel, the water content (w0), the liquid limit (wL) and the plastic limit (wp), the undrained shear strength (cu), the corrected SPT blow count (Nl and (Nl)60) and the shear-wave velocity (vs). this parameter was selected for the gravel for reasons of simplicity, even though a larger value could have been selected for the gravel above the water level. It also seemed reasonable to select a unique value of this parameter for the entire clay layer. The characteristic value of the shear modulus for very small strains, G0, was determined from the shear-wave velocity through G0 = pvs2 , where p is the soil density. The distribution of this modulus with depth was assumed according to the following expression where Gr0ef is the reference shear modulus at a vertical effective stress of 100 kPa. Two distinct values of G0ef were allocated to the entire layers of gravel and clay. No such parameter was allocated to the thin surface layer because it was assumed that its influence on the behaviour of the anchored wall was negligible. The full line in Fig. 2 shows the design characteristic shear-wave velocities, which result from the above assumptions. The characteristic value of the effective angle of internal friction ' for the gravel layer was determined through the correlation with (N1)60 proposed by Hatanaka and Uchida [13] ) = 200 15.4(Ni)60 (3) which gives a characteristic value of ) (3) rd = r ■ ^* r. = r ■ e (4) (5). The radius and the centre of the arbitrary half cone in the r-z plane are: R* = r* ■ sinh[($ — $»)tan*| (6) where R* , r* and 9* denote the cone diameter in the cross-section of the plane 9-z, and the polar coordinates of the apex of the hyperbolic half-cone. Figure 1. Cross-section of the failure mechanism. Figure 2. The scheme of the spatial failure mechanism. All hyperbolic half-cones whose infinite set represents the lateral surface of the failure mechanism are also kinematically admissible when the additional geometry condition is satisfied: Vr >- -Ar < r ; §0 < [t/2 — (7) which ensures that there exists no half-cone with its apex on the vertical wall (1-0, see Fig. 2) that could intersect the vertical wall under point (r0,). Since all the hyperbolic half-cones are kinematically admissible, then using the additional condition (7) the lateral surface, which is the envelope of the infinite set of all half-cones defined by expressions (8), (9) and (10), is also kinematically admissible. r^ = r,cosh[($ — $„)tan 0] — r,sinh[($ — $„)tan 0]sin(e,) (8), zs = r,sinh[(§ — §,)tan 4 ]cos(e„) (9), tanh[(§ — §,)tan0] + tan0 tan s, = arcsin(dR, /dr) = arcsin 1 + tanh[(§ — §,)tan0]tan0 tan (10). Considering r* = r1 and = §1 , expressions (8), (9) and (10) define the coordinates of the envelope on the leading half-cone. The coordinate zf of the lateral failure surface can be expressed: 1 sin § zf = zs, = r,sinh[(§ — §,)tan 4 ]cos(e») (11) Vr > rs A r < r1e (§—§1)tan 4 (12) Zf = z(r, §) = -^2rr cosh[(§ — §1)tan0 ] — r2 — r12 5 WORK EQUATION The considered failure mechanism on the width b is limited on the left by a vertical wall, on the right by a curved surface in the shape of a log spiral, and above by an even surface on which the surcharge can act. Both lateral surfaces are defined by the curved surfaces of the leading half-cone and the envelope of all the other hyperbolic half-cones (see Fig. 2). At each point on the so-formed failure surface the normal vector of the surface encloses with the plane r-z shear angle 0 and also defines the direction of the normal stress to the surface (see Fig. 3). dN = a dA, dT0 = dN tan0, dQ0 = ^dN2 + dT0 (13) where a and A denote the normal stress and the area of the lateral surface, and N and T^ denote the resultant x 0 Figure 3. The forces on the failure surface. values of the normal and the shear-stress components on the spatially formed failure surface. The shear cones on the failure surface define all the real or admissible directions of the forces dT^ and dQ^ (see Fig. 3). According to the upper-bound theorem the analyses should consider those directions of shear-strength activation that are kinematically admissible and ensure the highest possible value of the passive pressure for the chosen failure mechanism. The considered spatially formed failure body is certainly symmetrical in the symmetry plane r-9 that runs through the centre of the rectangular wall surface, and should be in equilibrium, considering all the forces that act on it. Certainly, all the forces dQ^ act in the plane r-z, and so they do not cause any momentum around the z axis, f epi ( )(cos£cos § sin£ tan § tan §„ sin3 § sin2 § 2 '1" f f (1 + 2zf / b)sin § r2drd§ = 0 §0 y0/cos§ which runs through the coordinate system's origin. Like in the 2D analyses, the equilibrium condition of all the momentums around the z axis is chosen for the work equation. From Fig. 3 it is evident that the maximum possible passive pressures arise when the shear force dT^ acts at each point of the failure surface in a direction that is defined by the cross-section of the normal plane through the centre of the hyperbolic half-cone and the tangent plane to the failure plane through the considered point. The coefficients of the individual parts of the passive pressure epy and e^ (let us call them the coefficients of passive pressure distribution) in the 3D problem are not constant along the wall height h. Certainly, they increase non-linearly with increased ratios of b/h. If y ^ 0, 0 # 0, S # 0 and q = c = 0 the work equation can be given in the following integral form: 0 r1e(J-J1)»n* )d§ — f f (1 + 2zf / b)sin § r2drd§- (14) § § x / sin § 10 And when y = 0, 0 # 0, S # 0, q # 0, and c = 0 it can be given in the following integral form: J v2(- cos* cos $ sin* )d$ - j(1 + 2zf /b)y2-sin$d$ = 0 (15) J 1 cos3 $ sin3 $ sin2 $ At point y = y0 and when b/h = the factor of the passive pressure distribution is epy = 0, and the appurtenant value of the factor of the passive pressure distribution epq is determined with a 2D model considering the geometry condition (7). The unknown functions efy = e (0, S, b/h) and epq = epq (0, S, b/h), which are the minimal possible solutions of the integral expressions (14) and (15) for all real ratios b/h, define the distribution of the passive pressures along the wall height. The minimum values of e^ and e can be determined numerically for an individual in advance for known ratios of b/h. The geometry model (height h = 1, unit weight y = 1 and ratio b/h) and the soil characteristics (shear angle 0 and the friction between the soil and the wall S) were used in our analysis. 6 NUMERICAL ANALYSIS AND RESULTS The values of the passive pressure distribution factors efy and eq are determined gradually from the top of the wall downwards for different ratios of b/h (b/h = 100, 75, 50, 25, 20, 16 down to 0.25), as can be seen in Fig. 4. It is assumed in the analysis that the passive pressures increase linearly between the individual calculation points upwards of the wall height. For each calculating point along the wall height there is an exactly determined spatial failure surface, which ensures the smallest possible value of the factors of the passive pressure distribution, e„„ and e_, for the chosen ratio b/h. pY pq In step m of the passive pressure determination, the minimum values of the factors e° to e™-1 and ePq to e"71 p ! p ! pq pq are known from the preceding steps. The appurtenant known momentums can be determined with the expres- The numerical resolving of the integral equations (14) and (15) is performed by dividing the analysed region in the x-y plane into an arbitrary number of triangular and rectangular finite elements. These are suitable for Gauss's numerical integration as well as for the calculation of the integral over the area of the plane y = y0 , where one-dimensional Gauss's numerical integration elements (see Fig. 4) are used. fm-i m- , , ,(yi+1 -y-1) Ay,-1 + y, + y+1) . * fp-, = £ep-,(y-y^——cos*--—+— x,sin<5 (16), fT = e^ (^) [cos * ^^ - x0 sin * 1 + ' pq m—1 £ • ¡=1 pq (y,+1- y,-1) pq Ay,-1 + yt + y,+1) . , cos*--x sin * (17), $ $ $ Figure 4. Passive pressure distribution and the scheme of the numerical integration. /■m _ PY — /•m _ pq — ( ym y m—1) Pi 2 ( ym y m—1) 2 cos ë(ym—1 + 2ym ) — r sinë cosë( ym—1 + 2 ym ) — v sin ë (18), (19), where f' 1 and f' 1 define the momentums of the J pj J pq already known values of the passive pressures, and f™ and fpq , the momentums of the passive pressures for e'^ = 1 and e'q = 1, according to the origin of the coordinate system x-y-z. The appurtenant momentum of the unit weight of the ground (y = 1) and the surcharge (q = 1), above the failure surface are determined using expressions (20) and (21). gp, —-EE Kwik (1 + 2*; / bK sin 0 >k, (20) j—1 k—1 where A'yy denotes the area of the triangular or rectangular element j in the plane x-y (see Fig. 4), wjk is the weight coefficient for Gauss's integration point k, zf is the coordinate z on the envelope of the hyperbolic half-cones, rjk is the radius of the integration point k on element j in the plane x-y, and is the appurtenant angle of the radius rjk. In the numerical integration of the considered problem in plane x-y, 514 rectangular and 42 triangular elements with 9 and 6 Gauss's integration points were used (see Fig. 4). — -EE LW (1- -2z"k / b)rlk sin 0 lk (21), where Lxy denotes the length of a one-dimensional integration element on the ground surface y = y0 in the plane x-y (see Fig. 4), wlk is the weight coefficient for Gauss's integration point k, zf is the coordinate z of the integration point on the envelope of the hyperbolic half-cones in the plane y = y0, rlk is the radius of the integration point k on element I in the plane x-y, and Figure 5. Set of spatial failure surfaces. l—1 k—1 is the appurtenant angle of the radius rlk. In the numerical integration of the considered problem in the plane, 42 one-dimensional integration elements with 3 Gauss's integration points were used. The unknown values of the passive pressure distribution factors are determined using: em — PY g - f &PY J p P fp em — pq g - f &pq J p pq PY fp (22), pq In the numerical procedure determining the minimal value of the passive pressure distribution factors e'mi and e"mq , the starting failure surface in the optimization procedure is determined with the initial values of the parameters ^ and $2, which should satisfy the following boundary conditions: *0 >0, y0 >0, > (n/2)-< (23). Mathematical optimization was used to determine the unknown parameters §1 and $2 of the critical failure surface, which defines, in the considered calculation step, the minimal value of the unknown factor of the passive pressure distribution, emY the wall. and em , at the toe of pq The Solver Optimization Tool (Microsoft Excel) with the generalized-reduced-gradient method was used in the minimization process. The result of the gradual determination of the passive pressure distribution factors from the top of the wall downwards are the numerical values of the factors epY and epq, and a set of spatial failure surfaces that are presented in Fig. 5 for the case when 0 = 40° and S/0 = 1. The values of the factors of the passive pressure distribution, epy and epq , for different values of 0 , 8/0 and b/h are presented in Figs. 6 and 7. The values of the comparative passive pressure coefficients, K'pY and K'pq, and the distances of the handling points of the resultants, ay , and a.q , from the surface of the backfill soil are presented in Tables 1 and 2. The values of the handling points are given for individual shear angles and given ratios b/h, where the numerically obtained results for different shear ratios S/0 do not deviate by more than 0.5% from their average value. The appurtenant values of the substitutive coefficient and the distances of the resultants from the surface of the backfill soil are determined with the expressions (24) to (27). m m 1000 q 100 : 10 S/0 =1 1000 0.25 0.5 1 2 b/h 5 10 25 1000 1000 100 10 : S/0=o -------^ 0.25 0.5 1 2 b/h 5 10 25 Figure 6. The factors of passive pressure distribution epy for different values of 0 , S/0 and b/h. m—1 Ky = Y,(e'pi (y — yo)( y+1— y—1)+em (ym— y>)(ym— ym—1) (24) i=1 m—1 aY = ^(eP7( y — yo)( y.+1— y.—1)( y.—1 + y. + y.+1—3 y,)/3) + '=1 (25) emr,( ym—yo)( ym—ym—1)( ym—1 +2 ym— 3 yo)/3 m—1 K M = e°pq ( y1 — yo)/2 + ^ (e„ ( yi+1 — yi—1 )/2 + empq ( ym — ym—1 ) /2 (26) .=1 m—1 = e0pq ( y — yo)( y1 — yo)/6 + ^(ePq ( y.+1 — y._1)( y.—1 + y{ + y.+1 — 3 yo)/6) + i=1 C (ym — ym—1)(ym—1 + 2ym — 3yo)/6 (27) In the analyses of the spatial stability problems the theorem of corresponding states (Caquot and Kérisel 1948, Soubra and Regenass 2000) is still valid. The comparative coefficient of the passive earth pressure due to cohesion ( K'pc ) can be determined by using the comparative coefficient of passive earth pressure due to the surcharge ( K'pq ). K Pc =■ KM — 1/cos(<^) tan(^) (28) The values of K*pc for the purely cohesive soil (c > 0 and 0 = 0) with different ratios of ca/c and with a centre of gravity of epc the pressures measured from the top of the wall are given in Table 3. Table 1. Values of K*py for different values of the parameters f, 8 and b/h with the centre of gravity of the epy pressures measured from the top of the wall b/h f 8/f center of (deg) 0 1/3 1/2 1/3 1 gravity 15 3.6279 4.1527 4.3864 4.6365 5.1458 0.712 20 5.3933 6.4350 7.0004 7.5984 8.8877 0.712 25 7.9261 10.0117 11.2456 12.6106 15.7108 0.725 0.25 30 11.8711 15.9410 18.5946 21.7237 29.5912 0.730 35 20.0486 27.6728 33.3372 40.5337 60.4985 0.734 40 43.0671 57.3693 69.2222 85.9233 139.1175 0.738 45 116.4677 149.3839 177.6761 220.0038 375.5334 0.741 15 2.6711 3.0260 3.2012 3.3775 3.7313 0.699 20 3.7238 4.4311 4.8126 5.2130 6.0592 0.706 25 5.2089 6.5721 7.3711 8.2474 10.2079 0.712 0.5 30 7.4363 9.9335 11.6437 13.5710 18.3323 0.718 35 11.9863 16.5711 19.9663 24.2567 36.0095 0.724 40 24.6495 32.8951 39.7426 49.3782 79.7705 0.729 45 64.2513 82.5392 98.8063 121.9540 208.4111 0.734 15 2.1892 2.4647 2.6014 2.7383 3.0092 0.687 20 2.8862 3.4211 3.7071 4.0047 4.6270 0.693 25 3.8439 4.8396 5.4156 6.0425 7.4308 0.698 1 30 5.8439 7.0201 8.1663 9.4902 12.6793 0.704 35 7.9191 11.0232 13.2817 16.1123 23.6998 0.710 40 15.4367 20.6573 25.0023 31.1047 50.0679 0.716 45 38.1335 49.1298 58.6535 72.9298 124.8263 0.723 15 1.8479 2.1801 2.2961 2.4114 2.6380 0.678 20 2.4651 2.9099 3.1455 3.3890 3.8937 0.683 25 3.1579 3.9637 4.4240 4.9214 6.0131 0.687 2 30 4.1095 5.5370 6.4305 7.4462 9.8331 0.691 35 5.9395 8.2498 9.9397 12.0356 17.5319 0.696 40 10.8269 14.5378 17.6323 21.9686 35.1861 0.702 45 25.0634 32.4150 38.8084 48.4179 83.0027 0.707 15 1.7980 2.0064 2.1091 2.2105 2.4084 0.672 20 2.2106 2.5986 2.8021 3.0112 3.4411 0.674 25 2.7423 3.4307 3.8182 4.2341 5.1395 0.676 5 30 3.4441 4.6399 5.3717 6.1970 8.1113 0.679 35 4.7302 6.5872 7.9347 9.5873 13.8072 0.682 40 8.0584 10.8652 13.2127 16.4877 26.2246 0.686 45 17.2046 22.3721 26.8925 33.7108 57.8700 0.691 15 1.7483 1.9478 2.0456 2.1422 2.3290 0.670 20 2.1253 2.4935 2.6857 2.8827 3.2865 0.671 25 2.6034 3.2508 3.6129 4.0005 4.8411 0.672 10 30 3.2223 4.3370 5.0128 5.7721 7.5262 0.673 35 4.3270 6.0335 7.2664 8.7689 12.5571 0.675 40 7.1344 9.6407 11.7401 14.6610 23.2244 0.677 45 14.5765 19.0171 22.9160 28.8085 49.4740 0.681 15 1.6984 1.8886 1.9817 2.0736 2.2518 0.667 20 2.0396 2.3876 2.5686 2.7541 3.1334 0.667 25 2.4644 3.0696 3.4067 3.7670 4.5479 0.667 2D 30 3.0000 4.0321 4.6525 5.3492 6.9591 0.667 35 3.6901 5.4448 6.5993 7.9724 11.3870 0.667 40 4.5989 7.6224 9.8346 12.6613 20.3076 0.667 45 5.8284 11.1974 15.6822 21.9144 40.6109 0.667 Table 2. Values of K*pq for different values of the parameters 0, S and b/h with the centre of gravity of the epq pressures measured from the top of the wall b/h 0 S/0 center of (deg) 0 1/3 1/2 2/3 1 gravity 15 4.5768 5.2900 5.6081 5.9280 6.5309 0.609 20 7.0283 8.4840 9.2341 9.9907 11.5050 0.621 25 10.7112 13.5556 15.1308 16.8518 20.5128 0.630 0.25 30 16.2767 21.6076 25.0415 28.8780 37.6802 0.637 35 26.9519 36.2461 42.7186 51.2493 72.8142 0.642 40 52.6835 68.3932 81.3285 98.5627 151.3097 0.646 45 116.5809 149.1632 176.9264 217.2506 349.3207 0.650 15 3.1488 3.5825 3.8054 4.0025 4.4002 0.582 20 4.5447 5.4631 5.9108 6.3815 7.3299 0.596 25 6.5956 8.3312 9.2904 10.3297 12.5242 0.608 0.5 30 9.6385 12.8412 14.8706 17.1286 22.2471 0.618 35 15.4244 20.8022 24.6109 29.5968 41.8059 0.625 40 29.1399 37.9850 45.2969 55.1065 84.9392 0.632 45 62.8494 88.6238 95.8063 118.0501 191.9431 0.638 15 2.4316 2.7450 2.9041 3.0424 3.3246 0.555 20 3.3004 3.9297 4.2434 4.5715 5.2181 0.569 25 4.5344 5.7128 6.3552 7.0497 8.4886 0.581 1 30 6.3194 8.4512 9.7630 11.2191 14.4578 0.592 35 9.6528 13.0803 15.5514 18.7062 26.1839 0.601 40 17.3562 22.7883 27.2809 33.3902 51.4902 0.610 45 35.9646 46.3732 55.2288 68.4366 112.5270 0.619 15 2.0703 2.3184 2.4440 2.5563 2.7758 0.534 20 2.6753 3.1633 3.4047 3.6552 4.1404 0.545 25 3.5038 4.3975 4.8749 5.3812 6.4281 0.555 2 30 4.6598 6.2408 7.1807 8.2196 10.4893 0.565 35 6.7730 9.2194 11.0219 13.1920 18.2517 0.575 40 11.4902 15.1828 18.2730 22.5430 34.4973 0.585 45 22.5171 29.1956 34.9820 43.6475 72.5151 0.598 15 1.8513 2.0600 2.1639 2.2564 2.4323 0.516 20 2.2974 2.6965 2.8917 3.0896 3.4692 0.522 25 2.8854 3.5970 3.9688 4.3579 5.1504 0.529 5 30 3.6641 4.9014 5.6070 6.3809 8.0320 0.536 35 5.0419 6.9032 8.2543 9.7985 13.3404 0.543 40 7.9825 10.6136 12.8683 16.0172 23.9640 0.551 45 14.4186 18.8594 22.8019 28.7468 47.8074 0.560 15 1.7775 1.9726 2.0678 2.1531 2.3132 0.509 20 2.1705 2.5382 2.7169 2.8958 3.2362 0.512 25 2.6793 3.3254 3.6605 4.0075 4.7070 0.516 10 30 3.3322 4.4441 5.0746 5.7474 7.1781 0.521 35 4.4630 3.1341 7.3062 8.6307 11.6245 0.526 40 6.7878 9.0886 11.0668 13.7561 20.3108 0.531 45 11.7060 15.4077 18.7082 23.7798 39.2151 0.538 15 1.6984 1.8836 1.9685 2.0050 2.1969 0.500 20 2.0369 2.3770 2.5400 2.7022 3.0107 0.500 25 2.4644 3.0468 3.3495 3.6573 4.2786 0.500 2D 30 3.0000 3.9871 4.5357 5.1180 6.3569 0.500 35 3.6903 5.3540 6.3516 7.4707 9.9784 0.500 40 4.5990 7.4305 9.3077 11.5115 16.7775 0.500 45 5.8284 10.7914 14.4498 19.0443 30.7851 0.500 Table 3. Values of K*pc for f = 0° and different values b/h and cjc with center of gravity of epc pressures measured from the top of the wall. b/h K* ^ pc center of c/c = 0 c/c = 1/3 c/c = 1/2 c/c = 1/3 cjc = 1 gravity 0.25 6.9282 7.4720 7.7231 7.9631 8.4156 0.6051 0.50 4.5691 5.0257 5.2356 5.4287 5.7541 0.5854 1.00 3.3302 3.7248 3.8942 4.0439 4.2737 0.5611 2.00 2.6822 3.0314 3.1760 3.3024 3.4925 0.5391 5.00 2.2783 2.5938 2.7217 2.8321 2.9997 0.5192 10.00 2.1402 2.4427 2.5646 2.6693 2.8249 0.5104 Table 4. Comparison of K*py and K*pq with KpY and Kpq for different values f, S/f and b/h. ■ pq ■py pq $ O 8/$ KpY (Soubra and Regenass 2000) Kpy (Skrabl and Macuh 2005) K*py (proposed) b/h=0.5 b/h=1.0 b/h=10.0 b/h=0.5 b/h=1.0 b/h=10.0 b/h=0.5 b/h=1.0 b/h=10.0 20 0.5 5.04 3.85 2.75 4.92 3.76 2.69 4.81 3.71 2.69 1.0 6.99 5.14 3.35 6.35 4.77 3.30 6.06 4.63 3.29 40 0.5 53.74 31.22 14.75 41.55 25.92 11.85 39.74 25.00 11.74 1.0 131.75 77.02 26.42 90.36 55.48 23.93 79.77 50.07 23.22 $ O 8/$ Kpq (Soubra and Regenass 2000) Kpq (Skrabl and Macuh 2005) K*pq (proposed) b/h=0.5 b/h=1.0 b/h=10.0 b/h=0.5 b/h=1.0 b/h=10.0 b/h=0.5 b/h=1.0 b/h=10.0 20 0.5 6.22 4.45 2.75 6.10 4.35 2.73 5.91 4.24 2.72 1.0 8.06 5.54 3.17 7.79 5.44 3.27 7.33 5.22 3.24 40 0.5 74.26 43.48 12.82 49.68 29.50 11.33 45.30 27.28 11.07 1.0 130.19 73.35 21.22 104.80 61.07 21.36 84.94 51.49 20.31 7 COMPARISON WITH EXISTING SOLUTIONS In the literature only 2D analyses of the soil-pressure-limit values using different approaches are presented, while the research results for 3D cases are very limited. The research results of 3D passive pressure analyses according to the theorem of the upper-bound value have been presented in Soubra and Regenass (2000), and Skrabl and Macuh (2005). A comparison of the results for the coefficients K*py and K*pq for 8/$ = 0.5 and 1, $ = 20° and 40°, b/h = 0.5,1, 10 is presented in Table 4. A comparison of the results indicates that, particularly at greater shear angles and greater ratios of 8/$, the differences between the values of passive-earth-pressure coefficients for the compared failure mechanisms are the greatest. The coefficient Kpy for the proposed trans-lational failure mechanism is up to 11.72% smaller than the same coefficient for the failure mechanism (Skrabl and Macuh, 2005) when b/h = 0.5, while the coefficient Kpq is up to 18.95% smaller for the same b/h = 0.5. For higher ratios of b/h the difference gradually decreases, and when b/h > 20 the solutions are almost equal. 8 CONCLUSIONS This paper presents a procedure for determining 3D passive earth pressures according to the kinematic method of limit analysis. The set of three-dimensional kinematically admissible hyperbolic translational failure mechanisms with lateral surfaces bounded by the envelope of the hyperbolic half-cones is used to determine the critical distribution of passive pressure along a flexible retaining structure's height. The intensity of the passive pressures is gradually determined with the previously mentioned translational failure mechanisms from above, downwards. Thus, the critical distribution, the trust point and the resultant of the passive pressures that can be activated at the limit state for the chosen kinematic model are obtained. Using the diagrams presented in Figs. 6 and 7 it is possible to determine the actual critical distribution of the passive pressure limit values for any arbitrary practical case (within the frame of given assumptions) that is applicable in geotechnical design. The results of the numerical analyses indicate that, when considering the upper-bound theorem and the set of three-dimensional kinematically admissible hyperbolic translational failure mechanisms, the passive-earth-pressure coefficients are lower than in the case of the hyperbolic translational failure mechanism and the translational mechanisms published in the literature for b/h < 10. The upper-bound values of the comparative passive-earth-pressure coefficients with a calculated pressure distribution are lower than the existing solutions with an assumed pressure distribution obtained using the upper-bound method within the framework of the limit analysis. This means that the classically presumed passive-earth-pressure distribution in 3D analyses is not acceptable, because it can actually not be activated. Furthermore, the trust point of the passive pressures resultant is independent of the friction between the retaining structures and the soil. Therefore, the presented results are applicable in geotechnical practice. N % R, r length of one dimensional integration element l on the ground surface; resultant value of normal stress component on spatial formed failure surface; resultant value of stress on spatial formed failure surface; cone diameter; polar co-ordinate; polar co-ordinate of the apex of the curved cone; co-ordinate appurtenant to gradient angle of the envelope; resultant value of shear stress component on spatial formed failure surface; wjk wik weight coefficients for Gauss's integration point k; co-ordinate appurtenant to gradient angle of the E* envelope; co-ordinate of the section of the envelope and the cl leading cone shaft in plane r-9; Y S el $ 9 9, unit weight of the soil; friction angle at the soil-structure interface; gradient of the envelope in point rEl which is defined in an arbitrary plane r-z; angle of internal friction of the soil; polar co-ordinate; polar co-ordinate of the apex of the curved cone. LIST OF SYMBOLS REFERENCES b c Ca epc epY epq /'m PI cm gq area of triangular or rectangular element j in plane x-y; width of the retaining wall; cohesion; adhesion along the soil-structure interface; factor of passive earth pressure distribution of the cohesion influence; factor of passive earth pressure distribution of the soil weight influence; factor of passive earth pressure distribution of the surcharge influence; fPY momentums of passive pressures for e^ = 1 f™q momentums of passive pressures for e^ = 1 gY momentums due to unit weight of the ground; momentums due to surcharge loading on the backfill surface; h height of the retaining structure; K * comparative coefficient of passive earth pressure of pc the cohesion influence; K * comparative coefficients of passive earth pressure of the soil weight influence; K * comparative coefficient of passive earth pressure of pq the surcharge influence; [1] Blum, H. (1932). Wirtschaftliche Dalbenformen und deren Berechnung. Bautechnik, 10(5): 122-135 (in German). [2] Brinch Hansen, J. (1953). Earth Pressure Calculation, Danish Technical Press, Copenhagen. [3] Caquot, A., and Kérisel, J. 1948. Tables for the calculation of passive pressure, active pressure and bearing capacity of foundations, Gauthier-Villars, Paris. [4] Chen, W. F. (1975). Limit analysis and soil plasticity. Elsevier Scientific Publishing Company, Amsterdam, The Netherlands. [5] Coulomb, C. A. (1776). Essai sur une application des règles de maximis et minimis à quelques problèmes de statique relatifs à l'architecture. Mémoire présenté à l'académie Royale des Sciences, Paris, Vol. 7, 343-382 (in French). [6] Drescher, A., and Detournay, E. (1993). Limit load in translational failure mechanisms for associative and non-associative materials. Géotechnique, London, 43(3): 443-456. [7] Duncan, J. M. and Mokwa, R. L. (2001). Passive earth pressures: Theories and tests. Journal of L xy r re* Geotechnical and Geoenvironmental Engineering ronmental Engineering Division, ASCE, 126(11): Division, ASCE, 127(3): 248-257. 969-978. [8] Janbu, N. (1957). Earth pressure and bearing [22] Terzaghi, K. (1943). Theoretical soil mechanics. capacity calculations by generalised procedure of Wiley, New York. slices. In Proceedings of the 4th International Confer- [23] Vrecl-Kojc, H., and Skrabl, S. (2007). Determi-ence, International Society of Soil Mechanics and nation of passive earth pressure using three- Foundation Engineering: 207-213. dimensional failure mechanism. Acta Geotechnica [9] Lee, I. K., and Herington, J. R. (1972). A theoretical Slovenica, 4(1): 10-23. study of the pressures acting on a rigid wall by a sloping earth on rockfill. Géotechnique, London, 22(1): 1-26. [10] Kérisel, J., and Absi, E. (1990). Tables for the calculation of passive pressure, active pressure and bearing capacity of foundations. Gauthier-Villard, Paris, France. [11] Kumar, J., and Subba Rao, K. S. (1997). Passive pressure coefficients, critical failure surface and its kinematic admissibility. Géotechnique, London, England, 47(1): 185-192. [12] Michalowski, R. L. (1989). Three-dimensional analysis of locally loaded slopes. Géotechnique, The Institution of Civil Engineering, London, England, 39(1): 27-38. [13] Michalowski, R. L. (2001). Upper-bound load estimates on square and rectangular footings. Géotechnique, The Institution of Civil Engineering, London, England, 51(9): 787-798. ^ [14] Mroz, Z., and Drescher, A. (1969). Limit plastic- ^ ity approach to some cases of flow of bulk solids. Journal of Engineering for Industry, Transactions of the ASME, 91: 357-364. [15] Ovesen, N. K. (1964). Anchor slabs, calculation methods, and model tests. Bull. No. 16, Danish Geotechnical Institute, Copenhagen: 5-39. [16] Salençon, J. (1990). An introduction to the yield design theory and its applications to soil mechanics. European Journal of Mechanics - A/Solids, Paris, 9(5): 477-500. [17] Shields, D. H., and Tolunay, A. Z. (1973). Passive pressure coefficients by method of slices. Journal of the Soil Mechanics and Foundation Division, ASCE, 99(12): 1043-1053. [18] Skrabl, S., and Macuh, B. (2005). Upper-bound solutions of three-dimensional passive earth pressures. Canadian Geotechnical Journal, Ottawa, 42: 1449-1460. [19] Sokolovski, V. V. (1965). Static of granular media. Pergamon Press, New York. [20] Soubra, A. H. (2000). Static and seismic earth pressure coefficients on rigid retaining structures. Canadian Geotechnical Journal, Ottawa, 37: 463-478. [21] Soubra, A. H., and Regenass, P. (2000). Three-dimensional passive earth pressure by kinematical approach. Journal of Geotechnical and Geoenvi- odziv zasičenih zemljin na dinamično oga6m6nit6v STANISLAV LENART o avtorju Stanislav Lenart Zavod za gradbeništvo Slovenije, Odsek za mehaniko zemljin in geotehnična opazovanja Dimičeva 12, 1000 Ljubljana, Slovenija E-pošta: stanislav.lenart@zag.si Izvleček V članku sta predstavljena dva najbolj izrazita načina deformiranja dinamično obremenjenih zasičenih zemljin. Likvifakcija s tečenjem in ciklična mobilnost sta pojava, ki pritegneta posebno pozornost zaradi velikih deformacij, ki ju spremljajo. Potopitev nasipa železniške proge zaradi novozgrajenega akumulacijskega bazena na reki Savi v Boštanju in velik plaz, ki se je sprožil na področju Stože v Julijskih Alpah predstavljata študiji primerov v Sloveniji, kjer sta bili analizirani tudi likvifakcija s tečenjem in ciklična mobilnost. Upoštevana je bila dinamična obremenitev zaradi železniškega prometa in morebitne potresne obtežbe. Materiali iz obeh lokacij: meljni pesek in prodnato-peščeni melj, so bili uporabljeni v obsežnih laboratorijskih preiskavah, katerih namen je bil določiti postopek za modeliranje generiranja dodatnega pornega pritiska v dinamično obremenjenih z vodo zasičenih zemljinah. Novejše ugotovitve kažejo, da je sprememba pornega tlaka povezana s količino disipirane energije, ki jo določajo histerezne zanke dobljene z dinamičnim obremenjevanjem. Na osnovi eksperimentalnih rezultatov je bila predlagana enačba za empirično zvezo, ki definira generiranje pornega tlaka med dinamičnim obremenjevanjem. Enačba je sestavljena iz dveh delov, prvi opisuje generiranje nepovratnih sprememb pornega tlaka, drugi del pa opisuje prirastke in upade tlaka porne vode znotraj enega obremenjevalnega cikla, ti. povratne spremembe pornega tlaka. Z uporabo predlaganega energijskega nume-ričnega modela je mogoče določiti dejanske efektivne napetosti in s tem napetostno pot dinamično obremenjene zemljine. Predlagani model za porne tlake se lahko uporabi tudi za modeliranje deformacijskega obnašanja. Z eksperimenti je bilo ugotovljeno, da ima dinamično obremenjena zemljina po nekaj ciklih obtežbe na začetku cikla zelo nizko togost, ki pa se kasneje poveča. Deformacija, ki se razvije za časa trajanja te podajne faze zemljine, predstavlja glavnino deleža skupne deformacije. Pojav in trajanje te faze sta prav tako povezana z energijo, ki se disipira med cikličnim obremenjevanjem. Odnos med disipirano energijo in dodatnim pornim tlakom ter kratkotrajnim tečenjem med pojavom ciklične mobilnosti, omogoča preprosto modeliranje odziva dinamično obremenjenih zasičenih zemljin. Ključne besede likvifakcija s tečenjem, ciklična mobilnost, dodatni porni tlak, disipirana energija the response of saturated soils to a dynamic load STANISLAV LENART About the author Stanislav Lenart Slovenian National Building and Civil Engineering Institute, Section for Soil Mechanics and Geotechnical Monitoring Dimiceva 12, 1000 Ljubljana, Slovenia E-mail: stanislav.lenart@zag.si Abstract This paper presents the two most significant types of deformation behavior for dynamically loaded, saturated soil. Flow liquefaction and cyclic mobility deserve special attention because of the large deformations that accompany these two phenomena. The submergence of a railway-line embankment due to the newly built Sava-river accumulation reservoir in Bostanj and the large landslide that occurred in the Stoze area in the Julian Alps are case histories in Slovenia where flow liquefaction and cyclic mobility were analyzed. The dynamic loading caused by railway traffic and possible seismic activity were taken into account. Material from these two sites, silty sand and lacustrine carbonate silt, were used in extensive laboratory research, with the objective to define a procedure for excess pore-water pressure-generation modeling in dynamically loaded saturated soil. It has been found recently that the change of the pore-water pressure is related to the dissipated energy density calculated from the hysteresis loops caused by dynamic loading. Based on the experimental results an empirical equation defining the generation of pore pressure during dynamic loading has been proposed. The equation is divided into two parts: the first part describing the residual pore-water pressure generation, and the second part describing the increment and decrement of pore-water pressure within the load cycle, the so-called temporary pore-water pressure change. The proper effective stresses and thus the stress path of the dynamically loaded soil can be defined by using the proposed energy-based numerical model. The proposed pore-pressure model can also be used in deformation-behavior modeling. It was observed from the experimental results that after a few cycles of dynamic loading the saturated soil starts to exhibit a very low stiff- ness at the beginning of a load cycle, after which it begins to strengthen. The strain developed during this softening phase represents the main share of the total strain. The occurrence and duration of this phase are related to the energy dissipated during the cyclic loading as well, and the relation between the dissipated energy, the excess pore pressure and the short-term flow during cyclic mobility, give us an opportunity for a simple response modeling of the dynamically loaded saturated soils. Keywords flow liquefaction, cyclic mobility, excess pore pressure, dissipated energy 1 INTRODUCTION Due to the existing hydro-power potential, a chain of five new hydro-power plants is planned for the lower part of the Sava river. The construction of the first power plant started at the Bostanj site in November 2002 and it was completed in May 2006 [1]. The railway connection between Ljubljana and Zagreb runs along the Sava river. The construction of the accumulation reservoir caused the submergence of the railway-line embankment and raised questions about possible changes in the response of newly saturated soils resulting from the dynamic load caused by the railway. Therefore, the stability analyses considering a new ground-water level were required by the owner of the railway [2, 3]. Field and laboratory tests were carried out before the upheaval of the water in order to get the input parameters for the analyses. Two possible types of saturated soil behavior were in question: liquefaction and cyclic mobility. The following extensive research gives an opportunity for a detailed study of the liquefaction potential of silty sand from the lower Sava river at Bostanj. The other case of material susceptible to liquefaction was found in the northern part of Slovenia - in the Julian Alps. Lacustrine material of glacial origin has been recognized as very sensitive to various factors, such as water content and loading conditions, including seismic effects. A strong earthquake in April 1998, with its epicenter in the Krn mountains and a magnitude of MWA=6.0, caused the collapse of an approximately 100-m-long section of the shore of the 20-km-distant Lake Bohinj. Saturated lacustrine soils were thought to be present. The same earthquake, in amplified form, caused serious damage to several buildings in the village of Mala vas, near Bovec, which were founded on saturated lacustrine soil [4]. Quite soon after the above-mentioned occurrences, in November 2000, a very severe landslide occurred in the area called Stoze. It resulted in debris flow in its lower part. Layers of lacustrine carbonate silt of a relatively small thickness were observed in the material displaced during the landslide, and the question arose as to whether the presence of these layers was responsible for the landslide [5]. An investigation of the static and dynamic liquefaction potential of the lacustrine carbonate silt was initiated as a result [6, 7]. Several findings from the research mentioned above are presented in this paper. Saturated silty sand from the Sava river at Bostanj and the lacustrine carbonate silt from the Julian Alps are materials that forced the author of this paper as well as the geotechnical society in Slovenia to accept the danger of liquefaction as well as the occurrence of cyclic mobility in Slovenia as real possibilities. The findings from these two recent case histories have helped to introduce a more cautious treatment of dynamically loaded saturated soil response in daily practice. They are used in this paper as examples. 2 LIQUEFACTION OF SOIL When dynamically loaded saturated soils are being considered, the term liquefaction is very important. Liquefaction is defined as the transformation of soil from the solid to the liquid state. It happens as a consequence of increased pore pressure and a reduced effective stress, mostly in saturated cohesionless soils. When such soil is subjected to rapid loading, e.g., earthquake loading or another kind of dynamic loading, the pore water is unable to drain in a very short time period. The loading conditions might be understood, therefore, as the undrained loading conditions. If it is not too dense, a cohesionless soil subjected to cyclic loading, especially cyclic loading in the shear mode, has a tendency to densify. As the pores between the soil grains are filled with water, which cannot drain sufficiently, the generation of excess pore pressure occurs. Figure 1 shows the changes in the soil skeleton caused by cyclic loading, which results in excess pore-water pressure being generated. The term static liquefaction (flow failure) refers to the rapid increment of pore-water pressure followed by a sudden loss of strength after the peak value of the deviator stress is reached, until a residual/steady-state strength is reached. Flow liquefaction appears when the residual strength of the soil is smaller than the static shear stress required for the equilibrium of a soil mass. The liquefied stress state, in that case, is represented by the initial effective confining pressure, decreased by the excess pore pressure. 2.1 simplified procedures for an evaluation of soil's liquefaction potential Based on the findings of many previous case histories, with and without the occurrence of flow liquefaction, some simplified procedures were developed for evaluating the liquefaction potential in a specific case [8, 9, 10]. The procedures are based on different field measurements. The most widely used among them are SPT, CPT and shear-wave velocity measurements. The liquefaction potential for silty sand from the Sava river at Bostanj and lacustrine carbonate silt from the Julian Alps was evaluated in the manner of SPT and shear-wave velocity measurements. An earthquake with a magnitude of M=7.0 was proposed as the strongest possible type of dynamic load, and it can be seen from the results (Figure 2) that the estimated danger of the occurrence of liquefaction depends strongly on the procedure, and that the results differ. Thus, more detailed research was needed. 2.2 the state-criteria approach As described above, the occurrence of flow liquefaction depends upon the ability of a material to contract itself. This contractive behavior is the reason for the pore-pressure increase. Therefore, it is of interest to know where is the boundary, called the critical void ratio (CVR) line, between the contractive and dilative soil state. The states of the tested soil - silty sand from the Sava river and lacustrine carbonate silt from Julian Alps - were defined in terms of the void ratio. Tests at different effective confining pressures resulted in the CVR lines shown in Figure 3. Figure 2. The danger of liquefaction occurring for silty sand from the Sava river and lacustrine carbonate silt from the Julian Alps. The CRR are based on (a) SPT measurements, (b) shear-wave velocity measurements. Log of vertical effective stress, loga'^ [kPa] Figure 3. CVR line as a boundary between the loose contractive states and the dense dilative states of the two tested materials. The generation of an excess pore-water pressure causes a decrease in the effective stresses. The effective-stress conditions leading to the occurrence of flow liquefaction can be most easily described in stress-path space. The response of five saturated cohesionless specimens isotropically consolidated to the same void ratio and different effective confining pressures, in undrained stress-controlled triaxial compression is shown in Figure 4 [11]. Regarding the initial states of the specimens according to the steady-state line, specimens 1 and 2 exhibit dilative behavior when the shearing starts, while specimens 3, 4 and 5 exhibit contractive behaviors, which is necessary for flow liquefaction. FLS Steady state point > p- Figure 4. Initial conditions susceptible to either flow liquefaction or cyclic mobility (adapted from [11]). The flow liquefaction is initiated at the peak of the stress path in the case of the latter three specimens. The locus of the points describing the effective stress conditions at the initiation of flow liquefaction is a straight line [11, 12] called the flow liquefaction surface (FLS). All the specimens reach the same steady-state point as they have the same void ratio. The FLS is truncated at the level of the steady-state point as flow liquefaction cannot occur if the stresses are below this point. The FLS therefore marks the boundary between the soil states at which either flow liquefaction or cyclic mobility can occur. While cyclic mobility is described in more detail in the next section, the flow-liquefaction potential of two investigated materials, silty sand from the Sava river at Bostanj and lacustrine carbonate silt from Julian Alps, is estimated on the basis of triaxial test results. Samples of lacustrine carbonate silt (Figure 5) were reconstituted at different initial states, while the samples of silty sand (Figure 6) were intact. The tested state of the silty sand seems unsuitable for contractive behavior and thus flow liquefaction is also not expected. Loose samples of lacustrine carbonate silt contracts remarkably during shearing. To trigger the flow liquefaction the static shear stress should exceed the shear strength of a soil in the liquefied state. nj 0. 50 40 p l b 30 w £> 1« 20 o & 5 10 i 45 kPa, 0.38 49 kPa, D.47 {¡j c'n=4B kPa / e=0.54 J 10 20 30 40 50 60 Effective stress, p l=(tJ1+2a3)/3 [kPa] Figure 5. Undrained triaxial test of isotropic, consolidated, reconstituted samples of lacustrine carbonate silt from the Julian Alps. 200 n 0. t/t ¡/I E TS 0 L_ B s 1 210 = 140 70 0 5 10 Axial strain, e, [%] 0 15 \ e=C .654 N, cr'0=196 kPa e =0.702 / a-'0=52 kPa / \ e =0.689 | 1 1 \.cr'o=Q7 KPa 1 1 ■ s i 250 50 100 150 200 Effective stress, pl={a1+2a3)/3 [kPa] Figure 6. Undrained triaxial test of isotropic, consolidated, intact samples of silty sand from the Sava river at Bostanj. 3 CYCLIC MOBILITY Cyclic mobility deformations are not of the flowing type and thus the damage would normally be smaller, but still severe, than in the case of flow liquefaction. Deformations due to cyclic mobility are developed incrementally during cyclic loading. The main reason for the dramatic increase of cyclic mobility deformations is the loss of stiffness caused by a decrease of the effective stresses. An excess pore pressure generated during dynamic loading moves the stress path from its initial position in the direction of the failure envelope (Figure 7). If the cyclic stress is large enough, the steady-state strength might be exceeded during the cyclic loading. If this happens near the FLS, the effective stress path can touch the FLS. Momentary instability can occur therefore, leading to significant strain development. If the static shear stress is smaller than the steady-state strength, the strain generally ceases when the shear stress returns to the values below the steady-state strength. If steady-state strength is not exceeded during the cyclic loading, the effective stress path approaches the so-called phase transformation surface (PTS) [13]. The PTS represents a kind of boundary between the dilative and the contractive behavior of loaded soil. Above the PTS a dilative tendency increases the effective confinement (and consequently the shear strength), while below the Figure 7. Generation of excess pore pressure due to the cyclic loading causes movement of the stress path from its initial position in the direction of the failure envelope. PTS the soil exhibits a contractive behavior and thus a tendency to generate excess pore pressure. Youd [14] clearly described the rearrangements of soil grains that happen in cyclically loaded soil when the stress path approaches and crosses the PTS from one side or the other. A significant shear strain may develop without an appreciable shear stress at the moment when the PTS is crossed (Figure 8). This, almost flowing, behavior of the soil when the stress path meets the PTS causes serious problems in the numerical modeling of cyclic mobility phenomena [15]. When the cyclic stresses are larger than the static shear stresses, stress reversal occurs. Thus, each load cycle includes compression and extension loading. Any excess pore pressure generated during the cyclic loading causes the movement of the stress path in the direction of a zero effective stress (origin of q-p graph). This state is called the initial liquefaction [16]. When the stress path reaches it, only further oscillations along the compression and extension portions of the drained failure surface are possible due to the continuation of the cyclic or monotonic loading [11]. Figure 8. A typical stress path for cyclically loaded soil and the shear-strain relationship when it crosses the phase transformation surface (adapted from [15]). It can be seen from the results of the two tested materials (Figure 5, Figure 6) that flow liquefaction is hardly likely to occur. There is a much higher risk of cyclic mobility, especially in the case of silty sand from the Sava river at Bostanj, which was found in the railway embankment. However, the dynamic load of a passing train might cause an increase in the pore pressure, leading to the occurrence of cyclic mobility. Two down-hole arrays with accelerometers and pore water-pressure sensors were established to enable monitoring of the dynamically loaded investigated soil response during and after the upheaval of the water and the saturation of the soil (Figure 9). The aim of the research was to define an effective procedure for an evaluation of the excess pore pressure. This would help to define a stress path during the dynamic loading and thus help with the prediction of any deformations. Figure 9. Railway embankment with a down-hole array and an accelerometer before the upheaval of the water. from contractive to dilative behavior and vice versa happen when the stress path crosses the PTS. Cyclic loading causes the irrecoverable contraction of the soil skeleton (Figure 1), which in the case of an undrained loading condition is accompanied by the permanent generation of excess pore pressure. It is obvious that if a proper model for pore-pressure changes during dynamic loading were to exist, it would enable simple access to the modeling of the response of dynamically loaded saturated soils. The idea for the solution was taken from metal-fatigue analyses, for which purpose a cumulative damage hypothesis was developed [17]. Using this hypothesis an irregular dynamic loading can be converted into an equivalent damaging quantity, which makes it possible to evaluate a stress path's movement from its initial position, approaching the phase-transformation surface and the failure envelope. The energy dissipated during dynamic loading is chosen as the equivalent damaging quantity. Quantification by seismologists of the energy released during earthquakes and the use of the energy dissipation in performance-based design in structural earthquake engineering argues for the use of energy for an excess pore-pressure evaluation procedure. The first steps toward the energy approach to excess pore-pressure evaluation were made using the relationships between the energy released during earthquakes and the sites where liquefaction occurred [18, 19]. Nemat-Nasser and Shokooh [20] presented governing differential equations relating energy dissipation to the densification of dry samples and to the generation of excess pore-water pressure in saturated samples. The dissipated energy density W was defined generally at time t as 4 AN ENERGY APPROACH TO EVALUATING THE EXCESS PORE PRESSURE The stress path moves from its initial position due to changes in the shear stress, which are caused by loading, and due to the effective pressure decreasing during the cyclic loading (Figure 7). The effective pressure usually decreases due to the excess generation of pore pressure. We can be sure that an excess generation of pore pressure during the cyclic loading actually leads to stress-path movements and thus to changes in the soil strength and the stiffness. Soil-grain rearrangements connected with soil-structure changes cause the pore pressure to increase when the soil contracts and decrease when it dilates. Changes 1 r W(t) — — / ajdej , (1) a0J0 where aj and e^ denote the stress and incremental strain tensors, respectively, t is the time in which the total dissipated energy is in question, and a'0 is an initial effective confining stress. In the case of the laboratory cyclic-loaded test results the dissipated energy density is defined as the area bounded by the hysteresis loops of the stress-strain curve (Figure 10). Complementing the theoretical framework, several laboratory tests [21, 22, 23] as well as field measurements [24] have been performed to prove the relation. A typical general form, derived from the proposed expressions, could be written as ru — a-Wb, (2) Figure 10. The dissipated energy per unit volume for a soil sample in case of cyclic triaxial test results. where ru denotes the pore-water pressure ratio (=Au/a'0), Au is the pore-water pressure change, while a and b are functions of the soil type, the relative density of the soil, the stress conditions, the initial soil-state parameters, etc. It should be mentioned again that the pore-water pressure change from Eq. 2 is caused only by the soil particles being rearranged and that this is a permanent change. Lenart [25] divided the excess pore-water pressure generated during the dynamic loading into two parts: the temporary pore-pressure change and the residual pore-pressure change. Temporary pore-pressure changes can be observed as oscillations of the pore pressure in a normal pore-water pressure curve. They are caused by the transmission of compressive stresses onto the pore water. The origins of the pore-pressure oscillations and their effect upon the deformation behavior of the soil during dynamic loading are described in another paper [26]. Knowing the time function F(t) of a normal component of dynamic loading to which the soil is subjected, it is possible to write the equation for evaluating the pore-pressure ratio, ru (Eq. 3). The parameters kr and kt are the residual and temporary pore-pressure parameters, which depend upon the type of soil and its state. Their evaluation procedure is described in more detail in [25]. Using the least-square method, the best agreement between the proposed relationship and the empirical results was found if a dissipated energy density, W, in Eq. 3 is raised to the power of e/10, where e means the base of the natural logarithm. As it is based on the dissipated energy density, the pore-pressure ratio calculated with Eq. 3 is independent of the loading frequency or the rate impacts in the case of the strain range typical for cyclic mobility or liquefaction [26]. ru =(W-[kr + ktF(t)] (3) The proposed equation was tested in the case of two investigated materials. Figure 11 shows the approximated pore-water pressure changes compared to the experimental results in the case of the undrained cyclic triaxial test of the lacustrine carbonate silt sample. The loading took the form of a sine wave with the frequency of the loading being 1 Hz. To prove the proposed equations' independence from the frequency and the rate of loading, an irregular loading test was performed on a sample of silty sand from the Sava river. The loading simulated a seismic load recorded during the Petrovac earthquake in Montenegro in 1979. Good agreement was obtained between the results of the numerical analysis and the results of the laboratory tests (Figure 12). Similar pore water pressure response due to dynamic loading was observed also, when other kind of materials were tested. An interesting increase in share of temporary pore water pressure changes was noticed in the case of highly porous snail soil [27]. Lacustrine carbonate silt. Stoze frequency of loading = 1 Hz 0,0E+D0 1.0E-02 2,0E-Q2 3:0E-Q2 4,CE-02 5,01-02 Dissipated energy ratio, w [-] Figure 11. Pore-water pressure curve during an undrained cyclic triaxial test of lacustrine carbonate silt. Figure 12. Results of an experiment compared to the modeling of an excess pore-water pressure generation. short-term flow in compression Figure 13. Evaluation of the energy dissipated during soil softening. 5 STRESS-STRAIN RELATION MODELING If the generated excess pore-water pressure during the dynamic loading of saturated soil can be evaluated in the proper way, most of the work is done already. This makes it possible to evaluate the exact position of the stress path during the dynamic loading and thus to take into account the progressive degradation of the stiffness and the strength of the soil in an effective stress analysis. The remaining problem is how to treat the phases of sudden soil softening and the significant regain in soil stiffness during cyclic loading at large deformations. Large deformations are limited by the soil hardening when the stress path crosses the PTS. A recent study [25] showed that this highly yielded segment can be limited by the amount of dissipated energy. Figure 13 presents an evaluation of the energy dissipated during the softening phase. It has been found through research [25] that the dissipated energy during the softening phase in a single load cycle is linearly related to the residual pore-pressure ratio. Using this finding one can define the residual pore-pressure ratio at which soil softening due to the cyclic mobility effect starts. The pore-pressure model presented in the previous section and the relation between the pore pressure and the energy dissipated during short-term flows when the cyclic mobility occurs were used [28] for simple stressstrain relation modeling. The results in case of the tested lacustrine carbonate samples are presented below. Figure 14. Typical simulation of the pore pressure (a), stress states (b), displacements (c) and stress-strain relation (d) for a cyclic triaxial test of a reconstituted sample of lacustrine carbonate silt. 6 CONCLUSION Due to the large accompanying deformations and the possibility of severe damage, the response of dynamically loaded saturated soils has long attracted attention. Flow liquefaction and cyclic mobility are two phenomena that are often confused with each other. Their characteristics are described using two case histories from Slovenia: lacustrine carbonate silt from the Julian Alps and silty sand from the Sava river in Bostanj. The effective stress decrease and the occurrence of large strain without any noticeable increase of the stress are common characteristics of flow liquefaction and cyclic mobility. An undrained loading condition, which is needed in both cases, is assumed to be present in saturated soil subjected to a rapid, dynamic load. Flow liquefaction appears suddenly when the residual strength of a soil is smaller than the static shear stress required for the equilibrium of a soil mass. On the other hand, cyclic mobility deformation develops incrementally during cyclic loading, mostly due to a decrease of the stiffness caused by a decrease of the effective stresses. It is important, therefore, to know the excess pore-water pressure generated during the cyclic loading of saturated soil, which impacts most upon the effective stresses in the soil. An energy approach to saturated-soils response modeling during a dynamic load is presented in this paper. The energy concept is based on the idea that part of the energy of a dynamic load is dissipated into the soil. The density of the dissipated energy is represented by the area of the hysteretic strain-stress loop. The dissipated energy density is related to the generated excess pore-water pressure. The latter was divided into the temporary and residual generated excess pore-water pressure. Such a formulation helps to model very precisely the pore-water pressure oscillations during irregular dynamic loading. The dissipated energy was evaluated during the soil-softening phase during cyclic mobility as well. Based on the observed linear relation between it and the residual excess pore-water pressure a promising attempt at modeling the response of dynamically loaded saturated soils was made. acknowledgments The author wishes to thank the Laboratory of Soil Mechanics and Prof. B. ¿lender of the Faculty of Civil Engineering, University of Maribor, Slovenia, for supporting the study with experimental results. The valuable help of G. Vilhar of the National Building and Civil Engineering Institute (ZAG) and the advice of Prof. J. Logar of the University of Ljubljana, were much appreciated. The financial support of the Ministry of Higher Education, Science and Technology of the Republic of Slovenia is gratefully acknowledged. RCFCRENCeS [1] Official web page of Holding Slovenske elektrarne HSE, http://www.hse.si, 2006 [2] Lenart S, Petkovsek B (2006). Assessment of liquefaction potential of a loose soil in the base of a railroad, Proceedings of the Xlllth Danube European Conference on Geotechnical Engineering, Ljubljana, Vol. 2, 965-970 [3] Lenart S, Fifer Bizjak K, Petkovsek B (2005). Report on the dynamic analysis of the critical cross section of the railway line at Bostanj accumulation basin. ZAG, P 810/05-750-2 (in Slovene) [4] Petkovšek, A. (2001). Geological-geotechnical Investigations of the Stože Landslide, Ujma, Ministrstvo za obrambo, Uprava RS za zaščito in reševanje, No. 14-15, 109-117 [5] Lenart S. (2006). Deformation characteristics of lacustrine carbonate silt in the Julian Alps. Soil Dynamics and Earthquake Engineering, 26, 131-142 [6] Lenart S., Žlender B. and Vilhar G. Liquefaction potential of lacustrine carbonate silt in the Julian Alps, Soil Dynamics and Earthquake Engineering, (to be published) [7] Žlender B. and Lenart S. (2005). Cyclic liquefaction potential of lacustrine carbonate silt from Julian Alps, Acta Geotechnica Slovenica, 2005/1, 22-31 [8] Seed, H.B. and Idriss, I.M. (1971). Simplified procedure for evaluating soil liquefaction potential, Jour. of the Soil Mechanics and Foundation Div., ASCE, Vol. 97, 1099-1119 [9] Andrus, R.D., Stokoe, K.H. (2000). Liquefaction Resistance of Soils from Shear-Wave Velocity, Journal of Geotechnical and Geoenvironmental Engineering, Vol. 126, No. 11, 1015-1025 [10] NCEER (1997). Proceeding of the NCEER Workshop on Evaluation of Liquefaction Resistance of Soils, Youd T.L. and Idriss I.M. (ed.), Technical Report NCEER-97-0022, National Center for Earthquake Engineering Research, State University of New York, Buffalo [11] Kramer, S.L. (1996). Geotechnical Earthquake Engineering, Publ. Prentice Hall, 348-422 [12] Vaid, Y.P. and Chern, J.C. (1983) Effect of static shear on resistance of liquefaction, Soils and Foundations, Vol. 23, No. 1, 47-60 [13] Ishihara, K. (1985). Stability of natural deposits during earthquakes, Proceedings, 11th International Conference on Soil Mechanics and Foundation Engineering, San Francisco, Vol. 1, 321-376 [14] Youd, T.L. (1977). Discussion, Soils and Foundations, 17, No. 1, 82 [15] Elgamal, A., Yang, Z., Parra, E., and Ragheb, A. (2003). Modeling of Cyclic Mobility in Saturated Cohesionless Soils, Int. J. Plasticity, 19, 883-905 [16] Seed, H.B. and Lee, K.L. (1966). Liquefaction of saturated sands during cyclic loading, J. Soil Mech. Found. Div, ASCE 92, SM6, 105-134 [17] Miner, M.A. (1945). Cumulative damage in fatigue, Trans. ASME 67, A159-A164 [18] Berrill, J.B. and Davis, R.O. (1985). Energy dissipation and seismic liquefaction of sands: Revised model, Soils and Foundations, Vol.25, No.2, 106-118 [19] Trifunac, M.D. (1995). Empirical criteria for liquefaction in sands via standard penetration tests and seismic wave energy, Soil Dynamics and Earthquake Engineering, Vol. 14, 419-426 [20] Nemat-Nasser, S. and Shokooh, A. (1979). A unified approach to densification and liquefaction of cohesionless sand in cyclic shearing, Can. Geotech. J., 16(4), 659-678 [21] Law, K.T., Cao, Y.L. and He, G.N. (1990). An energy approach for assessing seismic liquefaction # potential, Can. Geotech. J., 27, 320-329 [22] Liang., L. (1995). Development of an Energy Method for Evaluating the Liquefaction Potential of a Soil Deposit, PhD Dissertation, Case Western Reserve University [23] Figueroa, J.L., Saada, A.S., Liang, L., and Dahisaria, N.M. (1994). Evaluation of Soil Liquefaction by Energy Principles, Journal of Geotechnical Engineering, ASCE, Vol. 120, No. 9, 1554-1569 [24] Davis, R.O. and Berrill, J.B. (2001). Pore Pressure and Dissipated Energy in Earthquakes - Field Verification, Journal of Geotechnical and Geoenvi-ronmental Engineering, Vol. 127, No. 3, 269-274 [25] Lenart, S. (2006). Numerical model for pore pressure change evaluation in soils with high liquefaction potential, PhD Thesis, University of Ljubljana [26] Lenart, S. (2008). The use of dissipated energy at modeling of cyclic loaded saturated soils, Earthquake Engineering: New Research, Nova Science Publishers, Chapter 8 [27] Zlender, B. and Trauner, L. (2007). The dynamic properties of the snail soil from the Ljubljana marsh, Acta Geotechnica Slovenica, 2007/2, Vol.4, 48-61 [28] Lenart, S. and Logar, J. (2006). Numerical modelling of vertically loaded saturated soil, Proc. of XIII. Danube-European Conference on Geotechnical Engineering, 83-88 gr66nova funkcija tangencialno obrcmcnjc-nega horizontalno slojevitcgr polprostora TOMAŽ PLIBERŠEK Ln ANDREJ UMEK o avtorjih Tomaž Pliberšek Univerza v Mariboru, Fakulteta za gradbeništvo Smetanova 17, 2000 Maribor, Slovenija E-pošta: tomaz.plibersek@uni-mb.si Andrej Umek Univerza v Mariboru, Fakulteta za gradbeništvo Smetanova 17, 2000 Maribor, Slovenija E-pošta: umek@uni-mb.si izvleček V članku je obravnavan nov pristop k evaluaciji integralne predstavitve Greenove funkcije za slojevit pol-prostor, ki je na površini obremenjen s harmonično tangencialno točkovno silo. Enačbe gibanja razvežemo z uvedbo valovnih potencialov. Nevezane enačbe rešujemo tako, da jih s pomočjo predpostavljene odvisnosti od kotne koordinate in uvedbo Hanklove transformacije, pretvorimo v navadne diferencialne enačbe. Njihova rešitev in uvedba inverzne Hanklove transformacije vodita do integralne predstavitve za pomike na površini pol-prostora. Te nato izvrednotimo s pomočjo predlaganega tristopenjskega postopka. Najprej integrande razcepimo na dva dela, od katerih prvi po izvedeni integraciji vodi do singularnih pomikov, drugi del pa do regularnega dela pomikov. Opazimo, da je ta drugi del, potem ko uvedemo ustrezno izbrane razvejiščne reze in veljavnost integrandov s pomočjo analitičnega nadaljevanja razširimo s pozitivne realne osi Hanklovega parametra na njegovo celotno kompleksno ravnino, možno izvrednotiti s pomočjo konturne integracije za poljubno število slojev. Tako je ta del pomikov podan z rezidui integranda izvrednotenih v njegovih polih in integrali s končno razdaljo integracije vzdolž razvejiščih rezov. Te slednje z lahkoto izvrednotimo z zaželeno natančnostjo s pomočjo numerične integracije Izbrani nume-rični rezultati dopolnjujejo matematična izvajanja. Ključne besede elastodinamika, valovanje, Greenova funkcija, horizontalno slojevit polprostor, tangencialna koncentrirana sila grccn's function for trngcntirly loaded horizontrly lrycrcd hrlf-sprce TOMAŽ PLIBERŠEK and ANDREJ UMEK About the authors Tomež Pliberšek University of Maribor, Faculty of Civil Engineering, Smetanova 17, 2000 Maribor, Slovenia E-mail: tomaz.plibersek@uni-mb.si Andrej Umek University of Maribor, Faculty of Civil Engineering, Smetanova 17, 2000 Maribor, Slovenia E-mail: umek@uni-mb.si Abstract The topic of this paper is a novel evaluation of the integral representation of the surface Green's function for a layered half-space, loaded on its surface by a harmonic tangential point force. The equations of motions are reduced to wave equations by the introduction of wave potentials. The Hankel transform is applied to them and they are consecutively solved leading to the integral representation of surface displacements. They are consecutively evaluated by the proposed three step procedure. First the singularity is extracted. It is further noted that so obtained integrals, after suitably chosen branch cuts and analytic continuation of integrands are introduced, can be evaluated by contour integration for an arbitrary number of layers. They are, therefore, expressed by number of residues at the poles of integrand and the integrals along finite portions of the branch cuts. The latter ones can easily be evaluated to any desired accuracy leading to a closed form solution. Some numerical results corroborating the presented approach are given. Keywords elastodynamics, elastic wave propagation, Green's function, horizontaly layered half-space, horizontal point load 1 INTRODUCTION Modeling the elastodynamic characteristics of soil is required in a number of engineering problems, e.g. dynamic soil-structure interaction, dynamically loaded foundations etc. The soil is geometrically in prevalent number of cases modeled as a half space, which is, to be more realistic, endowed with some structure. As a starting point to determine the elastodynamic characteristics of soil one can use the fundamental solution or the Green's function. The use of the fundamental solution, which is known from the literature, results in the integrals over the whole surface of the half-space, interface between the soil and the superstructure and over the interfaces between the materials with different elastodynamic properties. Some of them are of infinite or semi-infinite extend. Their evaluations, which can be performed only numerically, represent the major difficulty of this approach. For practical evaluation the integration area has to be made finite, what results into introduction of the fictitious boundary, where the radiation conditions should be satisfied. The latter represent a demanding and not completely resolved problem e.g. Premrov [1]. The Green's function approach leads us to the integrals extending across the interface surface between the soil and the superstructure only. Their evaluation is straightforward and relatively easy to perform, ones the Green's function is given. The difficulty of the problem is transferred to the determining and evaluating the Green's function itself. The problematic of the homogeneous as well as the layered half-space has drawn the attention of considerable number of authors, not all of them can be mentioned here. The first elasticity solutions for whole-and half-space problems, static as well as dynamic ones have been obtained by Kelvin [2], Boussinesq [3], Cerruti [4], Lamb [5], Mindlin [6] and the others. The basic results on layered media were presented by Ewing et al. [7]. In more recent times the authors sought the solutions by two basically different approaches. On one side there are the approximate methods e.g. Luco's ray method [8, 9], Kausel's [10] thin layer method as well as BEM and FEM methods, their combinations e.g. Gaitanaros et al. [11], Triantafyllidis [12], Auersch [13] and refinements e.g. Aubry et al. [14]. On the other hand we have methods leading to exact solutions in the form of integrals of semi- or infinite extend e.g. Vostroukhov [15] and Jin [16]. Their evaluation in the concept of FFT [17] concludes the latter approach. It, however, appears that the use of FFT like integration is a successful method to evaluate the Hankel transforms inversion integrals in the cases, where there are no singularities in the resulting displacement field. It must be however noted that in the case of the Green's function the vital information about its singularity comes from the portion of integration path, where the integration variable is very large or infinite. This fact makes in the case of the Green's function the FFT like evaluation of integrals less efficient. Kobayashi et al. [18] considered the homogeneous, elastic half-space and succeeded to reduce the semi-infinite integrals representing the Green's function to a part containing the singularity, the residue at the Rayleigh pole and the integrals of finite extend along the properly chosen branch cut, which can be easily evaluated with any desired accuracy. Štrukelj et al. [19] succeeded to modify the Kobayashi's approach, so that it could be applied to the problem of vertically loaded layered half-space. The authors have first derived the Green's function for a single layer [20] and have demonstrated that under the assumption of infinite thickness of the layer their solutions lead to the Kobayshi's ones. Our decision to focus on the layer has also been motivated by the investigations of the mechanical properties of soils e.g. Žlender et al. [21] and [22], which are given at a point and can be easily generalized to a layer. This paper continues and completes the method to determine the surface Green's function for a layered half-space loaded with the harmonic force acting in any direction. It is motivated by the previous works by Kobayashi et al. [18, 23], Štrukelj et al. [19] and Pliberšek et al. [20]. The load acting in a general direction is decomposed into normal and tangential components with respect to the surface of the half-space, so that the problem of the latter one has to be dealt with only. First the general equations of motion for a single layer in Hankel transform domain are derived by adapting the Vostroukhov's [15] approach and then transformed back into the geometrical domain. These single layer solutions are then combined into the solution for a layered half-space making use of the continuity conditions on the interfaces between the layers and boundary conditions on the surface of the half-space. It is proven that the basic mathematical properties of these solutions do not depend on the number of layers. They are the same for a homogeneous half-space as well as for the half-space with any number of layers. The inverse Hankel transform integrals appearing in the Green's function can be therefore for any number of layers expressed with part proportional to r— , integrals of finite extend along the properly chosen branch cut and some residues at the poles of the integrand. The closed form solutions, obtained by our approach, for the Green's function of the tangentially loaded layered halfspace are consecutively presented graphically for some selected number of cases. 2 METHOD OF ANALYSIS 2.1 governing equations for a layer We consider a horizontally layered half-space, which consist of n layers on a half-space, as shown on the Fig. 1. The material properties of each layer and of the underlying half-space are assumed to be isotropic and homogeneous. Shear modulus b, Poisson's ratio vi, mass density p{ and the dumping coefficient b are the material constants of i-th layer. The homogeneous halfspace is labeled as the layer number H. The global, cylindrical co-ordinate system and the local cylindrical co-ordinate systems having their origins at the top of each layer are introduced. The model is subjected on its free surface to a concentrated tangential point load, which varies harmonically in time. Without loss of generality, it is assumed that it acts in the direction of positive x-axis. The governing equation for each homogeneous, elastic layer is the well known [24] reduced Navier's equation of motion in frequency domain: bv2ut + h±h-VV• u, = -w2ut ; i e[ 1,2,...,n,H]. (1) Pi Pi The boundary conditions on the free surface of the halfspace can be in the cylindrical coordinates written as: 4 (r, w) = a„ (r,0, w) • cos(^) = - ^MlM • cos(^) (2) l2j (r, f, z, w) = aêz (r, z, w) • sin (f ) = 2 • n • r Q(")• S (r) . 2 • n • r azz>1 (r, f ,0, w) = 0 . (4) sin(f) (3) where 8 (r ) is the Dirac delta function. The continuity equation along the interfaces of consecutive layers, where perfect bonding is assumed, and radiation conditions for r complete the definition of boundary value problem under consideration. The geometry of the above defined boundary-value problem is axi-symmetric. Therefore the 9-dependence of the problem is governed by the 9-dependence of the loading only. Due to very simple 9-dependence as a first step in the solution procedure reduced stresses: azj (r,z,w) = azj (r,z,w)- cos(^) ; i g[1,2,...,«,H] (5) anzi (r,z,w) = ai,zj (r,z,w)- sin($) ; i g[1,2,...,«,H] (6) azzj (r, z, w) = azzi (r, z, w)-cos($) ; i g[1,2,...,«,H] (7) and reduced displacements: ur t (r,§,z,w) = ur t (r,z,w)- cos('ff ) ; i g[ 1,2,...,«,H] (8) u,ffJ (r, •&, z, w) = (r, z, w)- sin($) ; i g[ 1,2,...,«, H ] (9) uz. (r,z,w) = uz. (r,z,w) - cos(^) ; i g[ 1,2,...,«,H] (10) are introduced. In the second step we make use of the Helmholtz wave potentials [25] to decouple the equations of motion (1): ui = V-( +Vx-4>i ; i G [ 1,2,...,n, H ] . (11) The vector potential ^ should in addition satisfy the constraint condition: V^ =0 . (12) As in the case of displacements and stresses we introduce the reduced wave potentials: ( (r,ê,z,w) = ( (r,z,w) -cos(ê) ; i G [ 1,2,...,n,H] (13) j (r,ê,z,w) = j (r,z,w) - sin(ê) ; i G[ 1,2,...,n,H] (14) ^ (r,ê,z,w) = (r,z,w)-cos(ê) ; i G [ 1,2,...,n,H] (15) (r, ê, z, w) = (r, z, w) -sin(ê) ; i G[1,2,...,n, H ]. (16) The substitution of Eqs. (13) to (16) into equations of motion (1) yields: d2^ dr2 1 r dr 1 - r2 • & d2& dz2 ) + K i ■ ^ = 0 ; i g[ 1,2,. .,«, H (17) d2 dr2 1 r d dr 2 - + r d2 dz2 • Vr 2 - +— • Vv,, + r -kl. • Vr . = 0 ; i G[ 1,2,... n, H ] (18) d2 dr2 1 + -r d dr 2 -- + r d2 dz2 • ' 2 - + -• Vr r • Ve. = 0 ; i G[ 1,2,... «, H ] (19) d2 dr2 1 + -• r d dr 1 - + r d2 dz2 • V, ,i + k 2 T ,i = 0 ; i G[ 1,2,.. ,«, H (20) Vr ,i~ -Âi + r • d'Vr.i dr 1 dVz,i dz = 0 ; i G[ 1,2,...,«,H] (21) We note that the above equation system contains two decoupled equations (17) and (20) and two still coupled equations (18) and (19). To decouple the latter ones we add them and subtract them and introduce two new reduced wave potentials Xi and Ki as: x, = Vr,, + Vv,; ; k = Vi- i>v (22) where i e [ 1,2,...,«, H] and the equations (18) and (19) are replaced by: d2 +1 Br1 r dr Dz2 d2 + 1 ± Br2 r dr r2 dz' x, + k2u • X, = 0 ; i G[ 1,2,...,«,H] (23) K + k2T;i • K = 0 ; i G [1,2,...,n,H] (24) and the compatibility condition becomes: 1 d(x, + K ) d'z,, K + r dr dz = 0 ; i G[ 1,2,...,«,H]. (25) The most efficient way to solve the boundary value problem under consideration is by the introduction of Hankel integral transform. 2.2 hankel transform and solutions in hankel transform domain To solve the equations of motion (17), (20), (23) and (24) for each layer with appropriate boundary conditions (2) to (4) and continuity conditions Hankel integral transform r ^ Ç : m fH" (Ç) = Hn [f (r )] = f f (r )• Jn (Ç r )• r • dr (26) and its inverse £ — r : to f (r) = H- [fHn (£)] = f fH«(£)■ J« (£r)• £ • d£ , (27) 0 are introduced [26]. n is the order of the transform and J«(£ r j is Bessel function of the first kind and order «. To transform the equations of motion to their canonical form the integral transforms of different orders are employed. Equations (17) and (20) are transformed through Hankel transform of order 1, equation (23) of order 0 and equation (24) of order 2. This yields the following system of equations: d (Ç ) dz2 d 2XH (ç) dz2 d2KH2 (Ç) -(( - kl),H1 (Ç) = 0 ; i G[ 1,2,...,«,H] (28) 2 - k2TJ yx?' (Ç) = 0 ; i G[ 1,2,...,«,H] (29) dz -((2 - k^)-iKiH2 (Ç) = 0 ; i g[ 1,2,...,«,H] (30) d2Vz/' (Ç) dz2 -(( - kh H/1 (Ç) = 0 ; i G[ 1,2,...,«, H ] . (31) The solutions of the above equations can be written as: = cw • ea' z! + C2,i • e~arz> ; ; i G[ 1,2,.. .,«, H ] (32) XiH0 = C3,i • e t rz' + C4,i • e^' ; i G [ 1,2,.. .,«, H ] (33) KH2 = C5,; • e •z' + C6,i • e-* z' ; ; i g[ 1,2,.. .,«, H ] (34) VH = C7,i • e « •z'+ C8,i • e^* z ; , i G [ 1,2,.. .,«, H] , (35) where: a =>/Ç2 - K, ; t ^ Ç2 - k t2,; ; i G[ 1,2,...,«,H] . (36) 2 0 The constraint condition the transformed wave potentials must satisfy is obtained from equation (25). This yields: dfe (£, z,w £ • Kf2 (£,z,w) — £ • xH (£,z,w) + 2 • i e[ 1,2,...,«, H ] dz = o ; (37) The yet unknown integration constants Cj ; i g[ 1,2,...,8] and j e[ 1,2,.,«, H ] will be determined from boundary, continuity and radiation conditions. For this purpose the reduced displacement components are expressed through the transformed wave potentials as: 1 to ■r,i ( ^ ^ W") = - • J £ • 0 ;>2„t h, £•TH1 (£, z, w) + £ • AHH (£, z, w)- G>KiH2 (£,z,w) dz 2 d2Az/1 (£,z,w) + £ dz2 dKH2 (£,z,w) dz •Jo(£r)- —£TH1 (£,z,w)+ £• AH (£,z,w)+ (38) • J2 (£ r )[• d£ ; i e[ 1,2,...,«, H ] w( r, z w)=2 •j£ o U2 T H, £ • tH (£,z,w) — (£,z,w) dz + 2 d2 A H (£, z, w) £ ' dz2 dKH (£,z,w) dz +£ • A/1 (£, z, w) J 0 (£ r ) + |— £ •TH1 (£, z, w) + +£•A,,H (£,z,w; J2 (£r)}• d£ ; i e[ 1,2,...,«,H] (39) ; ,i( r, z w) = J£ d^iHl (£,z,w) d^/1 (£,z,o dz dz (40) —£ • (£,z,w)h J, (£r)• d£ ; i e[ 1,2,...,«,H]. From the above expressions strains are derived and introduced into the constitutive equation for linear, homogeneous and isotropic solid. The pertinent reduced stresses are then given as: I TO &rz,i (r, Z, w) = y •] j £ 2 £ d^Hi (£ ) 2 d3aAz,iHl (£,z,w) 2 • £--(£, z, w)------ dz ' £ dz3 d%H2 (£,z, w - TH dz2 — £2 • KH2 (£,z,w; Jo(£r)■ +1—2 • £ • dTH1 (£, z, w) + 2 • £ • d A ,iHl (£, z, w) dz dz (41) + + £2 • KiH2 (£,z,w)- J2 (£r) [• d£| ; i G [ 1,2,...,«,H] o (r,j) = V, J*-JU.*+ *2. fH (*,z,j) + 2 l Jo ' H ' dz 2 d3(*,z, j) 3%h2 (*,z,w dz3 d - H, + -2 - * ïH (*, z, w) + , dz + *2 . KH2 (*,z, j)l- J2 (*/ dz2 d2KH (*,z,w J o (* r ) + (42) ■ + 2 - * dt.H1 (*,z,J dz2 ' dz ; i e[ 1,2,...,«,H] to ^zz,i (r, z, j) = Vi J * ■ \ + 2 - y, V, d20,Hl (*,z,j) dz2 - *2-r1 (*, z, J + 2 - *2 - ïH1 (*,z, w) - 2 - * -dKH2 (*,z, j) + dz (43) - 2 d2t,,Hl (*,z,w) dz2 - J, (*r)- d* ; i e[ 1,2,...,«, H]. o The equations (38) to (43) permit us to construct a liner, algebraic equation system in integration constants Cij as unknowns. It is worth-while to note that in these equations due to constraint condition (37) only three wave potentials are appearing, what results in six integration constants per layer or underlying half-space. If a branch cut in the complex £-plane is introduced, which makes Re(ai )> 0 and Re(^)> 0 on the real positive £-axis, then the radiation conditions demand: C, H = C5, H = C7,H = 0 . (44) The boundary conditions on the surface of the halfspace are given by equations (2) to (4) and for the continuity condition perfect bonding between the layers and the underlying half-space i.e. the continuity of displacements and stresses across the interfaces, is assumed. Thus the boundary and continuity conditions result in 6 - « + 3 equations for the same number of unknown integration constants Cij. For the reason of in the forthcoming paragraph introduced solution procedure we must study some properties of the matrix of this equation system and its submatrices. To simplify the further mathematical derivation some dimensionless variables and constants are introduced: * = n - kn ; = ^ ; 7i = ; va = V kT1 kT1 V1 (45) a = kT1 -a ; fii = kT1 -¡3i ; = kT1 -h ; i e[1,2,...,n,H]. We begin the derivation of the equation system with the boundary conditions defined on the upper surface of the first layer shown on the Fig. 2. Figure 2. Material and geometrical properties of the first layer. We introduce the expressions for stresses (41) to (43) into the transformed boundary conditions (2) to (4) and take into account the equation (45). This yields: m, • kTi ■ {-2 • n • a • [ci.i- ci.2 ] + (2 + Ji ) • [ci,3 + ci.4 ] + +-• Ai3 •[Ci, - Ci, ][ = n I n (46) Mi • kTi "I-2 n • ai •[Ci.i - Ci,2 ] + (2 + (•[Ci,3 + Ci,4 ] + +2 • n • Ai •[Ci.5 - Ci.6]} = 0 (47) mi • kti 1 (Ai + 2 • Mi1 2 2\ , 2 ----(ai -n ) + n 2 Mi [Ci.i + Ci.2] + -n • Ai • [Ci.3 - Ci.4 ]- Ai2 • Ci.5 + Ci.< (48) = 0 The other equations follow from continuity conditions on the interfaces between layers and the nth layer and the underlying half-space respectively. The interface between the ith and the (i+1)st layer is depicted in Fig. 3. Figure 3. Material and geometrical properties of ith and (i+1)st layer. Introduction of equations for the displacements and the stresses (38) to (43) into continuity equations leads to: f i (a, + 2 • mi ) •(a,2 - n2 )+n2 12 Mi °"zz.i , - °~zz,i. „ \Zi =hi IZ:+i =0 C 2 • e-a>4> ] - n • Jt • IC 3 • eA4 - Ci4 • e-J'4 i-J2 I Ci, • ea'4' + Ci5 • eA 4 +Ci6 • 4 I i ((+i + 2• M,+i) ¡a2 ^ 2 -Mi+i.i •] 2---(a,+i- n )+n I 2 M,+i n • Ji +i •[Ci +i.3 Ci + i.4 ] Ji+i Ci + i.5 + Ci + i.6 [Ci + i.i + Ci + i.2] = 0 (49) = -Mi.i • {-2 • n • a, • [C,.i • ea 4 - C,.2 • e a 4 ]+(( + J [Ci.3 • e^ + Ci.4 • e-M' ] + + n J • C .5 • e^ - Ci.6 • e^ | [ + M,+i.i • {-2 • n • ai+i • (50) [Ci + i.i - C + i.2 ]+(( + Ji+i ) • [Ci + i.3 + Ct +i.4 ] + - • Ji+i • [Ci + i.5 - Ci+i.6 H = 0 - (J,, = ^ • {-2 • n • a • [cy • ea - c,,2 • e-"j + ( + 3, -3, •t, + +2 • n 3 • IC,5 • e3rt - C,,6 • e-3■tl I }- M,+1,1 • { -2 • n • a -3, -t' fc,3 • e34> + C, C+1,1 - C,+1,2 j + +(n2 + 3i+1) [C.+1.3 + C,+1,41 + 2 • n • 3,+1 • [c,+1,5 - c,+wJ}=0 (51) -n -fC,3 • e3 'h + C, — \ s^ a, -t, s-* -a, -1, = { a • IC ,1 • e , , - C,,2 • e , , - 3 C • e3 - C,,6 • e^ U + -3, t -{ a+1 •[C.+1,1-C.+1,2j-n•[C,+1,3 + C,+1,4J-3+1 •[C,+1,5-C,+1,6J } =0 (52) = { n •[C,,1 • e a ^ + C,,2 • e~a ^ J- 3 fC,3 • e3 ^ - C,4 • e■t' 1 + n - - 32 n, fC,5 • e3 h + C. -3, t | n • [Ci+1,1 + C,+1,2 j 3+1 • [C;+1,3 C,+1,4J + ,5 1 ¡,6 2 -T n - - • 3+1 n [C, + 1,5 + C, + 1,6 J =0 (53) - w z ,, + 1 ' lz, = I z,+! = 0 {-n • [C,,1 • ea^ + C,,2 • e^^ J + 3, •[c,,3 • e3^ - C,,4 • e~3^ J + n fc„ • e34> + + C,,6 • e -3, t - {-n • [c,+1,1 + c,+1,2 J + 3,+1 • [c,+1,3- c,+1,4 J+n C, +1,5 +C, + 1,6 =0 (54) where , e [ 1,2,...,«,« +1 = H] and the equation (44) has to be accounted for. The above created system of 6 • n + 3 equations can be written in the matrix form: A • C = b (55) The matrix A in the above equation is a band matrix with the bandwidth of maximum 9 terms. C is a properly ordered vector of integration constants C{ j. The right side of the system is a vector, where each term except the first one equals to zero. As we have limited our interest to motion at the surface, only the first six integration constants are needed. They are obtained by Cramer's rule. The value of the determinant of the matrix A is determined by its development along the first row. According to the fact, that only first six terms of the first row are different from zero the expression of the determinant |A| has the following form: j=6 A = £(-1)j+1 • «1,j • \Aj\ , (56) where a1,j ; j e [ 1,2,.. .,6] represent the first six non-zero terms offirst rowofmatrix [A] and |Aj|; j e [ 1,2,.,6} represent corresponding sub matrices. Due to the fact that the vector on the right side of the system equations (55) has the following form: b= Q( n • • kT ; b, = 0 ; j e[2,3,... ,6« + 2J, (57) the six constants, which determine the surface motion of the half-space, can be written as: cu = (-1) j+1 • A • Q( \A\ n • M1 • kT1 j e[ 1,2,.,66. ; k e [1,2,5,6,7,8] (58) j=1 Introducing the equations (32), (33), (35) and (58) into equations (38) to (40) and evaluating the latter at z = 0 lead to the displacements on the surface of the half-space as: *(a,0) = Of- /¡77' n ((|-A,|l-^-1 (+A4I) + 2n' Ml cT,1 K lAl n - 2 '(n2-1)'(A5 |-| A / 0 (a n)' dn+/ A • [-n • (l A1-| a,|) + (59) ^vn2-! •(( l+lA 41)+n •(( |-| A6 |)j. /2 (a n). dn} (a,0)= Q(-) - fTO(-n) 2n • M1 cT,1 K lAl n - 2 '(n2-1)'(A5 H A6 - 'i/^-n (1-1A2I j-Vn2^ (+| A4 + / 0 (a n). dn+/A '[-n (1 hi A21) + '(a. I+IA41)+n ( - A6 |)j. /2 (a n). dn} (60) u (a,0) = ^ -./A .ff-Y! (l+| A21)-n .(A. - A41) + n. M1 cT,1 J0 lAl 1 1 1 ' (l+| A6|)l /1 (a n). dn, (61) where following change of variable has been introduced: a = kT,1. r = . (62) cT,1 For further analysis it is worth vile to note that in the above equations only three distinct integrals appear. They are denoted as I1 , I2 and I3 and are given as: I1 = /A r(HA2l)-^ (+|A41) + n—(n2 -1 n (|-| A6/0 (a n). dn = / B1 (n). / (a n). dn (63) I2 = /A .[-n (H a2 O+V^2-1 .((A l+l a4 |)+n (5 l-l a6 to . /2 (a n). dn = / B2 (n). / (a n). dn (64) I= to - / A (+1 A21)-n (1-1A41)+ - VnT-1. ((A |-| A61)1. /1 (a n). dn = / B3 (n). /1 (a n). dn. (65) The newly introduced functions Bi (n), i = 1,2,3 can be identified from the above equations. They also allow us to write the reduced surface displacements in a more compact form. They are given as: to to to to 0 to 0 0 0 ur (a,0)= Q( w) 2n - oj ( + A Q ( ) to to -- fb. (n) -J0(an) -dn + fb2(n)-j2(an) -dn 2n - J J ut ,1 (66) u, (a,0) = QM-•w-(-I, +12 2n - ¡¡1 cT Q(w) w 2n - cT to TO f B1 (n)- J 0 (a n) -dn+f B2 (n)- J2 (a n)- dn (67) Q( w) w T = Q( w) w t ,1 uz (a,0) = - —-I3 = n - ¡¡1 cT j n - ¡¡1 c. w f B3 (n) - J1 (an) - dn. (68) cT,1 0 3 eVfiLUfiTION OF THE INVeRSION INTEGRALS As it is known from the literature [27] and it was already demonstrated in the previous paper by the authors [19] the leading term of singularity depends on the local conditions only. Due to this fact the singularity at the surface of the layered half-space can be determined by considering the homogeneous half-space with the material properties of the uppermost layer, what considerably simplifies the analysis. The integrals f, i = 1,2,3 equations (63) to (65) are first reduced to the case of a homogeneous half-space and consecutively their limits are evaluated. This yield: w -Ic, = w -lim f a—0 J n n 2 ct 1 2 CT 1 ^ff—l \Fh (n) J0 (an)- dn oj -f11 J n- lim n n vn2—1' t ,1 0 w 3 - 2 7 vn2—1 |Fh (n) J (an) -dn (69) n 2 to 2 cT1 2 — 2 7 fJ 0 (a n) -dn = 2 — v. 2r 0 0 0 and in an analogous way: 2 c, -IS,2 = t,1 2c, t,1 -lim f a- 0 •vn2—1 jnF—1 \Fh (n) J0 (an)-dn = t1 2r (70) where: |Fh(n)\ = (2-n2—i)2 — 4-n2-'1 -Vn1—T7 . (71) The limit of the integral I3 as a — 0 is zero and therefore it remains regular for all values of a. The singulari- ties given by equations (69) and (70) are now subtracted from integrals I1 and I2 respectively and taken under the integral sign. This yields: I1 = I1 - IS,1 = JA {n 0Ax|-I A2 \)-4nf-1 ' (|A3 | +1A4 |) + + n-'(n2-1) OA5I-|A61)-(2-)J0(an)'dn to = J (n)'J 0 (a n)' dn (72) I2 = I2 - Is,2 = J A '{-n (-| a2 I)+Vnr-T 'OA3I+|a4 |) + +n' (A51-| A61)-vx}' j 2 (a n)' dn to = J B2 (n)'J (a n)' dn. (73) The horizontal displacements components are now given by: - / „, Q(— ) 1 W /-r -r ur (a,0) = ' - +--(I1 +12 2n' r cT a _ Q(w) f1 w 2n' ¡il I r cT1 1+—' JB1 (n)'J0 (an)' dn+JB2 (n)' J2 (an)' dn (74) (a,0)= Q(w) 1-v w 2n' ^ Q(—) 2n' ^ +--( +12 1-v W T.1 0 to to -JB1 (n)' J0 (an)' dn+JB2 (n)' J2 (an)' dn [. (75) TO 0 TO 0 0 0 0 0 We note that the singular terms, which appear in the horizontal displacement components, are given explicitly. The values of the integrals I1 and I2 are bounded for all values of the parameter a or r respectively. And what is also important for a numerical evaluation and of course for our later considerations B1 (n), B2 (n) and B3 (n) tend to zero as n ^to . 3.1 EXTENDING THE RANGE OF INTEGRATION To transform the integrals f, i = 1,2 and I3 into a form permitting their evaluation by contour integration in a complex n -plane we must make their integrands single valued and extend the range of integration from -to to to . Thus we note that the functions B1 (n), B2 (n) and B3 (n) are not single-valued due to the terms ai and 3 appearing in them. They are made single valued by introducing the branch cuts in the complex n -plane. In the selection of a suitable branch cut we are however limited by the following requirements: imposed radiation conditions, which require that the real parts of a{ (n) and 3 (n) are positive on the positive real n -axis; that it does not intersect the big semi-circle in the upper n -half-plane and by the demand that a{ (n) and 33i (n) are odd functions of n on the real n -axis. The latter is needed to extend the range of integration from semi-infinite to infinite. The branch-cut fulfilling the above stated requirements is shown in Fig. 4. Figure 4. Branch points of expressions and ßßi with introduced branch cut and the corresponding Rayleigh pole. For greater clarity of the figure some material damping has been assumed and expressed by complex shear module j = j • evf . We realize that the chosen branch-cut indeed makes ai (n) and fj. (n) odd functions of n on the real n -axis and does not intersect the big semi-circle in the upper n -half-plane. We further note that all terms of matrix A, equation (55), except the exponential functions: pi (n) =e =e ±äjtj _ ±iln -Y2 t q± (n) = ± ß>ti ±Jn2 11 e ' ' = e v ' ' and (76) which are neither odd nor even. To make the terms in matrix A determined with respect to evenness or oddness respectively we replace them by their analytic continuations, which do not change their values on positive n -axis and are even functions on the real n -axis. As functions satisfying the stated requirements we have chosen: ±nx3—F t, p. (n) = e 1 " = e ' ' and ±nJ7-W t: - q±(n) = e n = e±-t' . (77) The replacement of the functions p± (n) and q± (n) h equations (49)-(54) through p± (n) and q± (n) leads to the matrix A , which has on the integration path of inverse Hankel transform exactly the same values as the original matrix A. The terms of matrix A have interesting and for further analysis very useful properties. It can be seen from equations (46) to (48) that all the terms in the first three rows of matrix A are even functions of n on the real n -axis. The further rows of matrix A are coming from continuity conditions on interfaces between layers, six rows for each interface. From equations (49) to (54) it is easy to recognize that the first tree of these six rows contain only terms, which are even functions of n , and the next three only terms, which are odd functions of n on the real n -axis. Matrix A therefore has 3 • (n +1) rows, where all the terms are even, and 3 • n rows, where all the terms are odd and n is the number of layers. This implies that all the determinants |A| and |A;|, i s[ 1,2,...,6] are in the case of odd number of layers n odd functions and on the other hand in the case of even number of layers n are even functions on the real n -axis. We now derive functions B{ (n) exactly the same way as we have derived functions B{ (n) only that we make use of determinants |A| and |A.|, i e[1,2,...,6] , instead ofthe determinants |A| and A|, i e [1,2,.,6] . This yields: B, (n) = pA-[n-((-| A? Ij-Vn7-!-((+A |) + ? n—(n2 -(( - N) n j B? (n) = jA-[-n-(( - |A? D+Vn2—!-((+| a4 |) + (78) (79) and +n (I -| A |)J-v Bs(n) = n-k/n2-7? (Av,!+|A|)-n-|Aa|)-Vn7-!-((+|A6 A (80) where all the Bi (n), i €[ 1,2,3] are even functions of n on the real n -axis and Bt (n) = Bt (n), i € [1,2] and B3 (n) = B{ (n) on the real, positive n -axis. Taking into account these equalities, the equations (72), (73) and (65) can be written as: m m I = JB, (n) -Jo (an) - dn = JB, (n) J (an) - dn = J (81) 0 0 I? = /B? (n) -J? (an) - dn = / B? (n) -J? (an) - dn = (82) 00 and m m /3 = J B3 (n)- J, (an)- dn = / B3 (n)- J, (an)- dn = 6 3 .(83) 00 At this point we make use of integral representations of Bessel functions known from the literature e.g. Gradshteyn et al. [28]: Ji (an) = hi (an) + h (-an); i = ,,2,3 where: h0 (an) = — Jeiansin* dtf (85) (84) h, (a n) = , Hi (a n) (86) 1 1 h2 (a r,) = j- J ei(a n sin ^dti (87) and HJ (an'j is the Hankel's function of the first order and the first kind. Making use of the relationship (84) equation (8,) yields: I = I = Where we have made use of the change of variables n ^ -n and the symmetry of the function B, (n). In a very analogous fashion we obtain for the other two pertinent integrals given by equations (82) and (83) following expressions: I, = J, = m m Jb? (n) -J? («n) - dn = J (n)-h? (an) - dn = I? (89) 13 = J B5 (n)-h (an)-dn = J -65 (n)-h, (an)-dn = I3 . (90) The equations (88) to (90) clearly show that we have successfully replaced the original inverse Hankel transform integrals with the range of integration from 0 to +m with newly defined integrals having the range of integration from -m to +m . The integrals 1., i = ,,2,3 can be evaluated by contour integration in the complex n -plane as it will be shown in the next paragraph. 3.2 EVALUATION OF INTEGRALS BY CONTOUR INTEGRATION Integrals I,, I? and I3, are finally in the form permitting their evaluation by the contour integration in the complex n -plane. The most suitable contour is shown in Fig. 5. By the residue theorem it can be written: 1+ + I.r +1- + 4 + 4 + 4? = 2 n i Y res, ; i = ,,2, 3. (9,) It can be easily shown that the value of the integral along the big semi-circle in the upper n -half plane is identi- = J = J B, (n)-J 0 (an)- dn = J B, (n)- h0 (a n)- dn+J B, (n)- h0 (-a n)- dn 0 0 0 m -m m = J B, (n)- h0 (a n)- dn + J B, (-n)- h0 (a n)- (-dn) = J B, (n)- h0 (a n)- dn = J1,. (88) 0 0 m 0 0 cally zero. This is due to the behaviour of integrands, which are dominated by hi (an) functions. On the big semi-circle n can bi given as n = R eiV , where p takes the values from 0 to n . Equations (85) and (87) yield: Ml - Im(i7)>0 1 ^ lim h0 (a 7)= lim — | exp [iR(cos^ + i sin^)sin §1 d§ = ^ ' r—to n j 0 1 n = — lim I exp(-Rsin^> + iRcos^>)sin§d§ — 0 n r—to O 0 1 n lim h2(z) = lim — I expfiR(cos^ + isin^>)sin§ — 2i§1d§ = 'i—to r—to n j 0 1 n — lim | exp [(—Rsin^ + iR cos ^>)sin§ — 2i§l d§ — 0 . n R—TOO ^ I'll Im(7)>0 (92) (93) Figure 5. Integration path for the evaluation of integrals I, I2 and I . Some material damping is assumed to make the picture clearer. And the same behavior for h( (an) = H((() (a n'j, namely that lim H™ (R • e'V) 0 for 0 < vp < n is well docu- mented in the literature [29]. Therefore it can be concluded: ItR = 0; i = 1,2,3. (94) Taking in account the above equation the equation (91) can be rewritten as: 1 = I+ +1- = 2 n i £ res,. - - 4 - ^ ; i = 1,2,3. (95) It is clear from Fig. 5 that all three integrals appearing in the right hand term of the above equation have finite integration path. Therefore by the equation (95) our fundamental goal has been achieved. It can be further noted that if the integrals Iih( and Im are led along one and the other side of the branch cut and the value of the integral Ijr is equal zero. For the numerical calculation it is advantageous to express the integrals iih( and Iih2 through a sum of integrals of even shorter integration range stretching from one singularity on the branch cut to the other. These singularities are either branch points of functions ai and 3 defined through the equation (45) or the poles defined by the zeros of the determinant |A|, which lie on the branch cut. Introducing the equation (95) and considering the equations (88)-( 90) into equations (74), (75) and (68) yields the surface displacements as: 0 reS2 ) I1b1 I2b1 ^1b2 ^2b2 - , ^ Q(— ) I 1 — f„ . ur (a,0) = v y ••!- +---12ni • ^(res1 + 2n • ¿tj I r cr1 L (a,0) = Q(• I + — • f2n1 • ^(-reS1 + reS2 ) + V - ^M + Ab2 - ^2 1 [ 2n • ¿1 | r cT 1 L ^ JJ (96) (97) UZ (a,0) = ^^ — (2^reS3 - 13bl - 4 2 n • cT ,1 v (98) The above three equations are the final result of our analysis. In them the surface displacements of a layered half-space due to a tangential point load are given by a singular term, in the components, which become singular as r ^ 0 , a sum of the residues and a sum of integrals with the finite integration path. The latter ones can be evaluated numerically with any desired accuracy. We can further note that the residues represent the surface waves and the integrals are due to body waves. 4 NUMERICAL EXAMPLE As an illustrative example a one-layer half-space with the same geometrical and material characteristics as the one considered by Štrukelj et al. [19] has been chosen. It is shown in Fig. 1 with n = 1. On its surface a horizontal, harmonic point-load is applied. The material properties of this half-space are as follows: ratio of material densities in the underlying half-space and the layer is PhIP1 = 1.5, the ratio of shear modules ¡i1 = 2.0 , the Poisson's ratio for the layer is v1 = 1/3 and the Poisson's ratio for underlying half-space is vH = 1/4. The materials in the layer and in the underlying half-space are assumed to have no material damping. The ratio of the layer thickness h1 and the wave length of shear waves in the layer A1 is taken to be fy/ A1 = 2.0. As it can be seen from equations (96) to (98), the displacement components are expressed as sums of several terms. It is worthwhile to note that in the both horizontal displacements components the same terms appear with exception of the singular term. In the vertical displacement, however, completely different terms are forthcoming. The numerical effort can be, therefore, considerably reduced by computing first each of these terms separately and later combine them into displacement components as given by equations (96) to (98). In the Fig. 6 and 7 as an example three such characteristic terms are given. Our choice of the geometrical and material properties of the half-space is based on the fact that it is nearly impossible to obtain the data in the pertinent literature, with which our results could be compared to prove their valid- 0.15 0.05 -5> a Figure 6. Real parts of terms 2n i^^ res3 and 1, ity and accuracy. Therefore we have decided to make use of the principle of reciprocity in elastodynamics, which is well described in the literature e.g. by Achenbach [30]. It can be concluded from this principle, that the vertical displacements along the ray § = 0 in a half-space loaded with a horizontal, harmonic, unit point-load are equal to the radial displacements in an identical half-space loaded with a vertical, harmonic and unit point-load. Therefore we can state that for our above described choice of the layered half-space the displacement component uz (a,0), as given by equation (98), should equal the radial component of the surface displacements presented by Strukelj et al. [19]. It is however worth vile to note that this equivalence can not be seen from the two expressions before their numerical evaluation. The integrals in equation (98) have integrands, which are based on the determinant and sub-determinants of a 9x9 matrix A. The corresponding integrals are in the case of the vertically loaded half-space based on a 6x6 matrix. The real and imaginary parts of both displacement functions are shown in Fig. 8 and 9. It can be seen from both figures that they are in an excellent agreement with the results presented by Štrukelj et al [19]. The results of the evaluation of the radial displacement component ur (a,0) are shown in the Fig. 10 and 11. In exactly the same way the circumferential displacement component u 6 (a,0) can be evaluated through a different combination of terms appearing in ur (a,0). The results of this evaluation will not be presented in this paper. Figure 8. The real part of the displacement function (a,0) given by the dotted line and the real part of the function ur (a,0) presented through the solid line. Figure 9. The imaginary part of the displacement function uz (a,0) given by the dotted line and the real part of the function ur (a,0) presented through the solid line. Re[u,(a,0)) lm (u,(a,0H 5 CONCLUSION REFERENCES The displacement components on the surface of a horizontally layered half-space due to a tangential point-load were expressed through a combination of tree distinct Hankel's inverse integrals and trigonometric functions with the circumferential coordinate as argument. For the evaluation of these Hankel's integrals a novel three step procedure is employed. In the first step the singularity at the generic point from the integrals, where it exists is extracted and the resulting new integrals are made regular. In the second step we replaced the new integrands functions with their suitable analytic continuations, by which we were able to extend the integration range of Hankel's integrals to -to to +to . By this extension of the integration range we were in the last step permitted to evaluate them by contour integration. Through these three steps we were able to transform the Hankel's integrals into sum of three terms. The first one contains the singularity in the form C/r , the second one is given by a sum of the residues of the integrand and finally the third term consists of finite number of integrals along the suitable chosen branch-cut. The latter ones regular and finite in their integration range can be easily evaluated numerically. The results presented in this paper together with the our previous results, Strukelj et al. [19], constitute, what we believe, a robust and numerically efficient method to evaluate the displacements on the surface of the horizontally layered half-space due to a point force of any direction. The method of evaluation presented in this paper provides us with exact and closed form expressions for the singularities of displacement field, what makes our results very suitable to be used in soil-structure interaction problems. We are convinced that an even more efficient integration line around the branch cut from the one used in this paper can be developed. Before we could come up with a definite recommendation concerning the integration path more numerical research is needed. It is however believed that this problem is beyond the scope of this paper, where we succeeded to demonstrate that the Green's function for a layered half-space can be expressed as a combination of terms, which can be easily, especially in comparison with original Hankel's inversion integrals, and accurately evaluated. [1] Premrov M, Spacapan I. (2004). Solving exterior problems of wave propagation based on an iterative variation of local DtN operators. Appl. Math. Model., 28, 3, 291-305. [2] Kelvin (Lord) (1843). Displacements due to a point load in an indefinitely extended solid. Citation no. 66 in Love's Treatise on the Mathematical Theory of Elasticity, 4th ed. Dover Press, New York,. [3] Boussinesq J. (1878-1885). Citation no. 67 in Love's Treatise on the Mathematical Theory of Elasticity, 4th ed. Dover Press, New York. [4] Cerruti V. (1882). Citation no. 68 in Love's Treatise on the Mathematical Theory of Elasticity, 4th ed. Dover Press, New York. [5] Lamb H. (1904). On the propagation of tremors over the surface of an elastic solid. Phil. Trans. Royal Soc. Lond.; A 203, 1-42. [6] Midlin RD. (1936). Force at a point in the interior of a semi-infinite solid. J. Appl. Phys., 7, 195-202. [7] Ewing MW, Jardetzky WS, Press F. (1957). Elastic waves in layered media. McGraw-Hill Book Company New York. [8] Apsel J, Luco J.E. (1983). On the Green's functions for a layered half-space, Part I. Bull. Seismol. Soc. Am., 73(4), 909-929. [9] Apsel J, Luco J.E. (1983). On the Green's functions for a layered half-space, Part II. Bull. Seismol. Soc. Am., 73(4), 931-951. [10] Kausel E, Peek R. (1984). Dynamic loads in the interior of a layered stratum, An explicit solution. Bull. Seismol. Soc. Am., 74, 1459-1481. [11] Gaitanaros AP, Karabalis DL. (1988). Dynamic analysis of 3D flexible embedded foundations by a frequency domain BEM-FEM. Earthquake Engng. Struct. Dyn., 16, 653-674. [12] Triantafyllidis T. (1991). 3D time domain BEM using half-space Green's functions. Engng. Anal. Bound. Elem., 8 (3), 115-124. [13] Auersch L. (1994). Wave propagation in layered soils, Theoretical solution in wave number domain and experimental results of hammer and railway traffic excitation. J. Sound Vibration, 173(2), 233264. [14] Aubry D., Clouteau D. (1991). A regularized boundary element method for stratified media. In: Cohen G., Halpen L., Joly P., editors: Proceedings of the First International Conference on Mathematical and Numerical Aspects of Wave Propagation Phenomena, 660-668. [15] Vostroukhov A.V., Verichev S.N. (2004). Kok A.W.M., Esveld C., Steady-state response of a stratified half-space subjected to a horizontal arbitrary buried uniform load applied at a circular area. Soil Dyna. Earth. Engng., 24, 449-459. [16] Jin B., Lia H. (1999). Exact solution for horizontal displacement at center of surface of an elastic halfspace under horizontal impulsive punch loading. Soil Dyna. Earth. Engng., 18, 495-498. [17] Bringham E.O. (1974). The Fast Fourier Transform. Prentice - Hall, New York. [18] Kobayashi T, Sasaki F. (1991). Evaluation of Green's function on semi-infinite elastic medium. KICT Report No. 86, Kajima Technical Research Institute, Kajima Corporation. [19] Štrukelj A., Pliberšek T., Umek A. (2006). Evaluation of Green's function for vertical point load excitation applied to the surface of a layered semiinfinite elastic medium. Arch. Appl. Mech., 76, 465-479. [20] Pliberšek T., Štrukelj A., Umek A. (2005). Green's function for an elastic layer loaded harmonically on its surface. Acta Geotechnica Slovenica, 2, 4-21. [21] Žlender B., Lenart S. (2005). Cyclic liquefaction potential of lacustrine carbonate silt from Julian alps. Acta Geotechnica Slovenica, 1, 22-31. [22] Žlender B., Trauner L. (2007). The dynamic properties of the snail soil from the Ljubljana marsh. Acta Geotechnica Slovenica, 2, 48-61. [23] Kobayashi T. (1981). Evaluation of response to point load excitation on semi-infinite elastic medium (in Japanese). Tran. Archi. Institute., Japan, 302, 29-35. [24] Achenbach J.D. (1973). Wave propagation in elastic solids. North-Holland, Amsterdam, London. [25] Miklowitz J. (1980). The theory of elastic waves and waveguides. North-Holland, Amsterdam, New York, Oxford. [26] Sneddon I.H. (1972). The use of integral transforms. McGraw-Hill Book Company, New York, Toronto. [27] Pao Y.H., Mow C.C. (1971). Diffraction of elastic waves and dynamic stress concentrations. Crane Russak, New York. [28] Gradshteyn I.S., Ryzhik, I.M. (1994). Table of Integrals, Series, and Products, Academic Press Limited, London. [29] Abramowitz M., Stegun I.A. (1972). Handbook of mathematical functions. New York, Dover publications. [30] Achenbach J.D. (2003). Reciprocity in Elastody-namics, Cambridge University Press, Cambridge. NAVODILA AVTORJEM Članki so objavljeni v angleškem jeziku s prevodom izvlečka v slovenski jezik. VSEBINA ČLANKA Članek naj bo napisan v naslednji obliki: - Naslov, ki primerno opisuje vsebino članka in ne presega 80 znakov. - Izvleček, ki naj bo skrajšana oblika članka in naj ne presega 250 besed. Izvleček mora vsebovati osnove, jedro in cilje raziskave, uporabljeno metodologijo dela, povzetek izidov in osnovne sklepe. - Uvod, v katerem naj bo pregled novejšega stanja in zadostne informacije za razumevanje ter pregled izidov dela, predstavljenih v članku. - Teorija. - Eksperimentalni del, ki naj vsebuje podatke o postavitvi preiskusa in metode, uporabljene pri pridobitvi izidov. - Izidi, ki naj bodo jasno prikazani, po potrebi v obliki slik in preglednic. - Razprava, v kateri naj bodo prikazane povezave in posplošitve, uporabljene za pridobitev izidov. Prikazana naj bo tudi pomembnost izidov in primerjava s poprej objavljenimi deli. - Sklepi, v katerih naj bo prikazan en ali več sklepov, ki izhajajo iz izidov in razprave. - Literatura, ki mora biti v besedilu oštevilčena zaporedno in označena z oglatimi oklepaji [1] ter na koncu članka zbrana v seznamu literature. OBLIKA ČLANKA Besedilo naj bo pisano na listih formata A4, z dvojnim presledkom med vrstami in s 3.0 cm širokim robom, da je dovolj prostora za popravke lektorjev. Najbolje je, da pripravite besedilo v urejevalniku Microsoft Word. Hkrati dostavite odtis članka na papirju, vključno z vsemi slikami in preglednicami ter identično kopijo v elektronski obliki. Enačbe naj bodo v besedilu postavljene v ločene vrstice in na desnem robu označene s tekočo številko v okroglih oklepajih. ENOTE IN OKRAJŠAVE V besedilu, preglednicah in slikah uporabljajte le standardne označbe in okrajšave SI. Simbole fizikalnih veličin v besedilu pišite poševno (npr. v, T itn.). Simbole enot, ki sestojijo iz črk, pa pokončno (npr. Pa, m itn.). Vse okrajšave naj bodo, ko se prvič pojavijo, izpisane v celoti. SLIKe Slike morajo biti zaporedno oštevilčene in označene, v besedilu in podnaslovu, kot sl. 1, sl. 2 itn. Posnete naj bodo v kateremkoli od razširjenih formatov, npr. BMP, JPG, GIF. Za pripravo diagramov in risb priporočamo CDR format (CorelDraw), saj so slike v njem vektorske in jih lahko pri končni obdelavi preprosto povečujemo ali pomanjšujemo. Pri označevanju osi v diagramih, kadar je le mogoče, uporabite označbe veličin (npr. v, T). V diagramih z več krivuljami mora biti vsaka krivulja označena. Pomen oznake mora biti razložen v podnapisu slike. Za vse slike po fotografskih posnetkih je treba priložiti izvirne fotografije ali kakovostno narejen posnetek. PREGLEDNICE Preglednice morajo biti zaporedno oštevilčene in označene, v besedilu in podnaslovu, kot preglednica 1, preglednica 2 itn. V preglednicah ne uporabljajte izpisanih imen veličin, ampak samo ustrezne simbole. K fizikalnim količinam, npr. t (pisano poševno), pripišite enote (pisano pokončno) v novo vrsto brez oklepajev. Vse opombe naj bodo označene z uporabo dvignjene številke1. SEZNAM LITERATURE Vsa literatura mora biti navedena v seznamu na koncu članka v prikazani obliki po vrsti za revije, zbornike in knjige: [1] Feng, T. W. (2000). Fall-cone penetration and water content ralationship of clays. Geotechnique 50, No. 2, 181-187. [2] Ortolan, Ž. and Mihalinec, Z. (1998). Plasticity index-Indicator of shear strength and a major axis of geotechnical modelling. Proceedings of the Elev- enth Danube-European conference on soil mechanics and geotechnical engineering, Poreč, 25 -29 May 1998. [3] Toporišič, J. (1994). Slovenski pravopis. 2nd.ed., DZS, Ljubljana. PODATKI O AVTORJIH Članku priložite tudi podatke o avtorjih: imena, nazive, popolne poštne naslove, številke telefona in faksa, naslove elektronske pošte. Navedite kontaktno osebo. SPRCJCM ČLANKOV IN AVTORSKE PRAVICE Uredništvo si pridržuje pravico do odločanja o sprejemu članka za objavo, strokovno oceno mednarodnih recenzentov in morebitnem predlogu za krajšanje ali izpopolnitev ter terminološke in jezikovne korekture. Avtor mora predložiti pisno izjavo, da je besedilo njegovo izvirno delo in ni bilo v dani obliki še nikjer objavljeno. Z objavo preidejo avtorske pravice na revijo ACTA GEOTECHNICA SLOVENICA. Pri morebitnih kasnejših objavah mora biti AGS navedena kot vir. Rokopisi člankov ostanejo v arhivu AGS. Vsa nadaljnja pojasnila daje: Uredništvo ACTA GEOTECHNICA SLOVENICA Univerza v Mariboru Fakulteta za gradbeništvo Smetanova ulica 17 2000 Maribor Slovenija E-pošta: ags@uni-mb.si INSTRUCTIONS FOR AUTHORS The papers are published in English with a translation of the abstract into Slovene. FORMAT OF THE PAPER The paper should have the following structure: - A Title that adequately describes the content of the paper and should not exceed 80 characters; - An Abstract, which should be viewed as a mini version of the paper and should not exceed 250 words. The Abstract should state the principal objectives and the scope of the investigation and the methodology employed, it should also summarise the results and state the principal conclusions; - An Introduction, which should provide a review of recent literature and sufficient background information to allow the results of the paper to be understood and evaluated; - A Theoretical section; - An Experimental section, which should provide details of the experimental set-up and the methods used for obtaining the results; - A Results section, which should clearly and concisely present the data using figures and tables where appropriate; - A Discussion section, which should describe the relationships shown and the generalisations made possible by the results and discuss the significance of the results, making comparisons with previously published work; - Conclusions, which should present one or more conclusions that have been drawn from the results and subsequent discussion; - References, which must be numbered consecutively in the text using square brackets [1] and collected together in a reference list at the end of the paper. LAYOUT OF THE TEXT The text should be written in A4 format, with double spacing and margins of 3 cm, to provide editors with space to write in their corrections. Microsoft Word for Windows is the preferred format for submission. One hard copy, including all figures, tables and illustrations and an identical electronic version of the manuscript must be submitted simultaneously. Equations should be on a separate line in the main body of the text and marked on the right-hand side of the page with numbers in round brackets. UNITS AND ABBREVIATIONS Only standard SI symbols and abbreviations should be used in the text, tables and figures. Symbols for physical quantities in the text should be written in Italics (e.g. v, T, etc.). Symbols for units that consist of letters should be in plain text (e.g. Pa, m, etc.). All abbreviations should be spelt out in full on first appearance. FIGURES Figures must be cited in consecutive numerical order in the text and referred to in both the text and the caption as Fig. 1, Fig. 2, etc. Figures may be saved in any common format, e.g. BMP, JPG, GIF. However, the use of CDR format (CorelDraw) is recommended for graphs and line drawings, since vector images can be easily reduced or enlarged during final processing of the paper. When labelling axes, physical quantities (e.g. v, T) should be used whenever possible. Multi-curve graphs should have individual curves marked with a symbol; the meaning of the symbol should be explained in the figure caption. Good quality black-and-white photographs or scanned images should be supplied for illustrations. [3] Toporišič, J. (1994). Slovenski pravopis. 2nd.ed., DZS, Ljubljana. AUTHOR INFORMATION The following information about the authors should be enclosed with the paper: names, complete postal addresses, telephone and fax numbers and E-mail addresses. Indicate the corresponding person. ACCEPTANCE OF PAPERS AND COPYRIGHT The Editorial Committee of the Slovenian Geotechnical Review reserves the right to decide whether a paper is acceptable for publication, to obtain peer reviews for submitted papers, and if necessary, to require changes in the content, length or language. Authors must also enclose a written statement that the paper is original unpublished work, and not under consideration for publication elsewhere. On publication, copyright for the paper shall pass to the ACTA GEOTECHNICA SLOVENICA. The AGS must be stated as a source in all later publication. Papers will be kept in the archives of the AGS. TABLES Tables must be cited in consecutive numerical order in the text and referred to in both the text and the caption as Table 1, Table 2, etc. The use of names for quantities in tables should be avoided if possible: corresponding symbols are preferred. In addition to the physical quantity, e.g. t (in Italics), units (normal text), should be added on a new line without brackets. Any footnotes should be indicated by the use of the superscript1. For further information contact: Editorial Board ACTA GEOTECHNICA SLOVENICA University of Maribor Faculty of Civil Engineering Smetanova ulica 17 2000 Maribor Slovenia E-mail: ags@uni-mb.si LIST OF REFERENCES References should be collected at the end of the paper in the following styles for journals, proceedings and books, respectively: [1] Feng, T. W. (2000). Fall-cone penetration and water content ralationship of clays. Geotechnique 50, No. 2, 181-187. [2] Ortolan, Ž. and Mihalinec, Z. (1998). Plasticity index-Indicator of shear strength and a major axis of geotechnical modelling. Proceedings of the Eleventh Danube-European conference on soil mechanics and geotechnical engineering, Poreč, 25 -29 May 1998. NAMEN Revije Namen revije ACTA GEOTECHNICA SLOVENICA je objavljanje kakovostnih teoretičnih člankov z novih pomembnih področij geomehanike in geotehnike, ki bodo dolgoročno vplivali na temeljne in praktične vidike teh področij. ACTA GEOTECHNICA SLOVENICA objavlja članke s področij: mehanika zemljin in kamnin, inženirska geologija, okoljska geotehnika, geosintetika, geotehnične konstrukcije, numerične in analitične metode, računalniško modeliranje, optimizacija geotehničnih konstrukcij, terenske in laboratorijske preiskave. Revija redno izhaja dvakrat letno. AIMS AND SCOPE ACTA GEOTECHNICA SLOVENICA aims to play an important role in publishing high-quality, theoretical papers from important and emerging areas that will have a lasting impact on fundamental and practical aspects of geomechanics and geotechnical engineering. ACTA GEOTECHNICA SLOVENICA publishes papers from the following areas: soil and rock mechanics, engineering geology, environmental geotechnics, geosynthetic, geotechnical structures, numerical and analytical methods, computer modelling, optimization of geotechnical structures, field and laboratory testing. The journal is published twice a year. AVTGRSKE PRAVICE CGPYRIGHT Ko uredništvo prejme članek v objavo, prosi avtorja(je), da prenese(jo) avtorske pravice za članek na izdajatelja, da bi zagotovili kar se da obsežno razširjanje informacij. Naša revija in posamezni prispevki so zaščiteni z avtorskimi pravicami izdajatelja in zanje veljajo naslednji pogoji: Fotokopiranje V skladu z našimi zakoni o zaščiti avtorskih pravic je dovoljeno narediti eno kopijo posameznega članka za osebno uporabo. Za naslednje fotokopije, vključno z večkratnim fotokopiranjem, sistematičnim fotokopiranjem, kopiranjem za reklamne ali predstavitvene namene, nadaljnjo prodajo in vsemi oblikami nedobič-konosne uporabe je treba pridobiti dovoljenje izdajatelja in plačati določen znesek. Naročniki revije smejo kopirati kazalo z vsebino revije ali pripraviti seznam člankov z izvlečki za rabo v svojih ustanovah. elektronsko shranjevanje Za elektronsko shranjevanje vsakršnega gradiva iz revije, vključno z vsemi članki ali deli članka, je potrebno dovoljenje izdajatelja. Upon acceptance of an article by the Editorial Board, the author(s) will be asked to transfer copyright for the article to the publisher. This transfer will ensure the widest possible dissemination of information. This review and the individual contributions contained in it are protected by publisher's copyright, and the following terms and conditions apply to their use: Photocopying Single photocopies of single articles may be made for personal use, as allowed by national copyright laws. Permission of the publisher and payment of a fee are required for all other photocopying, including multiple or systematic copying, copying for advertising or promotional purposes, resale, and all forms of document delivery. Subscribers may reproduce tables of contents or prepare lists of papers, including abstracts for internal circulation, within their institutions. electronlc storage Permission of the publisher is required to store electronically any material contained in this review, including any paper or part of the paper. ODGOVORNOST Revija ne prevzame nobene odgovornosti za poškodbe in/ali škodo na osebah in na lastnini na podlagi odgovornosti za izdelke, zaradi malomarnosti ali drugače, ali zaradi uporabe kakršnekoli metode, izdelka, navodil ali zamisli, ki so opisani v njej. RESPONSIBILITY No responsibility is assumed by the publisher for any injury and/or damage to persons or property as a matter of product liability, negligence or otherwise, or from any use or operation of any methods, products, instructions or ideas contained in the material herein.