THE SYSTEMATIZATION OF ANALYTICAL AND NUMERICAL METHODS OF LANDSLIDE STABILITY CALCULATION

S obzirom na kompleksnost geometrije reljefa zemljine površi, kosine su među problematičnijim geološkim formama u geotehnici. Kosine karakteriše nagla promena geometrije terena (denivelacija), s predispozicijom promene ove geometrije usled dejstva različitih faktora. Najčešći i najsloženiji vid narušavanja tla i geometrije kosine odnosi se na stabilnost terena – bilo prirodnih padina ili veštačkih kosina. Svako narušavanje postojeće ravnoteže na padinama ili kosinama izaziva pomeranja pod uticajem gravitacije: klizanje, odronjavanje ili tečenje površinskog dela tla, ali i dubljih delova stenske mase. Za ovako uspostavljeno klizanje, u geološkoj i geotehničkoj terminologiji i nomenklaturi, ustaljen je termin – klizište [11]. Uslovi za nastanak i razvoj klizišta jesu: geotehnički, geološki, geomorfološki, hidrogeološki, meteorološki, vegetacioni, antropogeni, dejstvo zemljotresa, dejstvo akumulacija, vibracije usled saobraćaja i drugi. U poslednjih sto godina, zabeležen je znatan broj katastrofalnih klizišta, nastalih kao posledica dejstva zemljotresa, erupcije vulkana, nagomilavanja snega, višednevnih i intenzivnih kiša i uragana [16]. Zbog formiranja ovih klizišta, poginulo je nekoliko stotina hiljada ljudi koji su u najvećem broju slučajeva imali sagrađene


INTRODUCTION
Given the complexity of geometry of the relief of the earth's surface, slopes represent one of the problematic geological forms in geotechnics.The slopes are characterized by a sudden change in geometry of the terrain (denivelation) with a predisposition to the change of this geometry due to the effects of various factors.The most common and most complex type of soil disturbance and slope geometry is the stability of the terrain, whether natural slopes or artificial slopes.Any disturbance of the existing balance on the slopes causes displacement under the influence of gravity: sliding, erosion or flowing the surface of the soil, but also the deeper parts of the rock mass.For the established sliding, in the geological and geotechnical terminology and nomenclature, the term "landslide" is established [11].Conditions for the formation and development of landslides are: geotechnical, geological, geo-morphological, hydro-geological, meteorological, vegetation, anthropogenic, earthquake effects, accumulation effects, traffic vibrations, etc.
Cilj istraživanja prikazanog u ovom radu jeste da se detaljnije sistematizuju metode proračuna klizišta i algoritmi modeliranja, s posebnim osvrtom na numeričke analize stabilnosti.hundred thousand people were killed, who in most cases had their own homes built on or near the landslide.According to [29], the most catastrophic landslides recorded over the last hundred years, were: Haiyuan landslides in China 1920, Vargas tragedy in Venezuela 1999, Nevado del Ruizdebris flows in Colombia 1985, Nevados Huascaran debris avalanche in Peru 1970, North India flood mudslides in India 2013, Khait rock slide in the USSR in 1949 and the like.In the cadastral register of landslide records in Serbia, more than 2200 active, currently calm and reactive landslides have been recorded [15].There is a significantly lower number of currently calm, compared to active and reactive landslides.Very often, in practice, problems with the stability of the landslide are encountered, when after the conducted field investigations, laboratory-geomechanical tests, defining the causes and conditions of landslide formation, it is necessary to define repair measures.However, in order to successfully manage all project situations of the landslide stability analysis and landslide repair, it is necessary to have high quality mathematical models and landslide methods.Previous experience shows that there is a need for the implementation of more complex (more realistic) mathematical models for practical purposes and further improvement of existing landslide stability methods.
One of the first papers in which the aspects of several analytical and also numerical methods of landslide stability are adequately theoretically considered is the paper [9], where the classification was carried out taking into account the limit formulation and the displacement formulation.There are two options for the limit formulation: upper bound solution and lower bound solution, where the method of characteristics for displacement belongs to the group of upper bound solutions, and the method of characteristics for stress belongs to the group of lower bound solutions.The methods defined by the limit formulation are in fact the Limit Equilibrium Method (LEM), of which the upper bound solutions are most commonly applied.In relation to them, the methods defined by the displacement formulation are in fact Displacement Formulation Methods (DFM), or numerical methods.The paper [20] gives an overview of the numerical methods of landslide stability, using the Finite Difference Method (FDM) formulation.A study of the performances of several different landslide stability methods is presented in [27], while in [21] the methods of landslide stability are presented taking into account the deterministic approach, the theory of reliability and optimization.The use of numerical methods for analyzing the stability of slopes in the alternating serpentine rock mass was shown in [25], where, among other things, these parameters were used: the geological strength index and deformability of the rock mass, and for the criterion of failure Hoek Brown's criterion was applied.Consideration of the complex problem of landslide stability, but also from the aspect of hazard analysis, vulnerability analysis, risk assessment and risk management, are presented in [6], where, based on the aforementioned theories, a GIS integral model for landslide analysis is proposed.
The aim of the research presented in this paper is to further systematize the methods of landslide calculations and modelling algorithms with a special emphasis on numerical stability analyses.

GENERAL ON LANDSLIDE STABILITY CONSIDERATIONS AND FACTORS RELEVANT TO THE CALCULATION
Generally speaking, the slopes can be found in a stable equilibrium state, unstable imbalance state and indifferent semi-equilibrium state.A stable equilibrium state is characterized by the established relation of destabilizing and stabilizing forces, so if the effect of stabilizing forces is greater, the safety factor is greater.The unstable imbalance state is characterized by a disturbed relation between destabilizing and stabilizing forces, so the influence of destabilizing forces is dominant.The indifferent (indefinite) half-balance state represents a transition category between a stable equilibrium and an unstable imbalance state.The ratio of destabilizing and stabilizing forces of the indifferent semi-equilibrium state is significantly more disturbed than the stable equilibrium state, and it is sufficient that the small increment of destabilizing forces transforms this state into an unstable imbalance state.Given the complexity of the indifferent state and the possibility of a partial change in the slope geometry or the semi-forming of the landslide, the indifferent state is characterized by a set of several different half-equilibrium states.This is especially characteristic in situations where the disturbance of the stable equilibrium state occurs, without the complete movement of the sliding mass of the soil, but a new half-balance state can be established.A more detailed classification of the stable and instable states of the slopes, taking into account the transition categories, is shown in [24]: stable slope, potentially unstable slope, early demolition phase, medium demolition phase, partial or complete demolition and complete demolition, while the mechanisms of formation and development of landslides are: rotational model, translational model, model formed from different geometric shapes of blocks, model with sliding, rolling and falling of stone of different dimensions, model with significant sliding of the landslide soil, landslide formed due to heavy rainfall, landslide formed as voluminous torrential flow, landslide formed by soil flow and landslide formed by soil crawling with the appearance of cracks and clefts in the soil.One of the key parameters in landslide classification is the velocity of movement of the sliding mass, as well as the level of underground water in the soil.In general, it can be concluded that the landslides, which have a higher slope of the outer contour of the soil, have a higher velocity of movement of the sliding mass [18].This is due to the effect of gravitational forces.However, the consideration of the influence of the rate of movement of the sliding mass of the soil and the incremental increase in water in the soil requires the application of methods for estimating the stability of landslides in the time domain, which is considerably more complex than the usual methods of calculation.
The methodology of the analysis of the potential landslide consists of several segments: geodetic survey of terrain and data collection, geotechnical in-situ testing, the analysis of physico-mechanical parameters of soil in the laboratory and calculation of slope using mathematical methods in geo-technics.The methodology of the analysis of the formed landslide is based on a landslide repair project consisting of several nih rešenja, višekriterijumska optimizacija varijantnih rešenja i detaljna analiza tehnologije sanacije klizišta za optimalno izabrano rešenje.
Prilikom formiranja proračunskog modela klizišta, potrebno je razmotriti sve relevantne parametre i odrediti njihove vrednosti, budući da je konačno rešenje u segments: geodetic surveying of the terrain and data collection, geotechnical in-situ testing, analysis of physico-mechanical parameters of soil in the laboratory, reconstructive analysis of the previous landslide state, analysis of the factors that led to the formation of landslides, the consideration of variant solutions for landslide repair, calculations of variant landslide solutions using mathematical methods in geotechnics, economic analysis of variant solutions, multicriteria optimization of variant solutions, and detailed analysis of landslide repair technology for the chosen optimal solution.
During the landslide analysis, preliminary geotechnical in situ testing are carried out, by which, in the first place, an engineering-geological profile of the terrain is formed.Key geotechnical investigations carried out for the formation of the engineering-geological profile of the terrain are exploratory boreholes.Exploratory boreholes are made using core drilling technique, where soil samples are carefully classified for soil type identification according to the depth and the analysis of physico-mechanical characteristics of the soil.The execution of the exploratory borehole must be carried out deep enough to determine the sliding surface on the engineering-geological profile of the landslide.The number of required exploratory boreholes is in correlation with: landslide geometry, landslide dimensions, depths of sliding surfaces, geological variations, and the like.
When forming the calculated landslide model, it is necessary to consider all relevant parameters and determine their values, since the final solution is in a direct correlation with the selection and variation of parameter values.The relevant parameters of the calculated landslide model can be classified into five groups: parameters of the geometry of the calculated model, parameters of physical-mechanical characteristics of the soil, parameters of actions, special types of parameters and calculation parameters.The complexity of the slope geometry and the multi-layered soil depth are modelled by parameters of the geometry of the calculated landslide model.Since the layers can be complex geometries, and not only horizontal or approximately horizontal, the numerical methods of calculating the landslide are used in the calculation of complex multilayer soil geometry.In addition, these parameters include the parameters of the geometry of the structures located on or near the landslide, or they are constructions for which the landslide is being repaired, so that for them, the effect of the soil-structure interaction (SSI) is considered.The correct input of these parameters depends on the level of quality of the formed engineering-geological profile of the terrain.The parameters of the physical-mechanical characteristics of the soil are obtained from laboratory testing of samples, of which the following are distinguished: one-axial strength, direct shear strength, triaxial test and edometric compressibility test.For stability analyzes of slopes, parameters are important: soil weight, soil weight in saturated state, cohesion, internal friction angle, Young's elastic modulus, edometric modulus of compressibility, deformation module, Poisson's coefficient, reference shear modulus, dilatation, coefficient porosity and other.Depending on the type of constitutive model of soil behaviour, the relevant parameters are defined, whereas for the constitutive models describing the threedimensional stress state the number of parameters is considerably increased.Most often, Mohr-Coulomb's soil model is used as a constitutive soil model for analyzing landslide, while depending on soil type specificity, soft soil model or hardening soil model can be used, Cam-Clay model, Drucker-Prager model and others.There are also additional parameters that enhance the constitutive model of soil behaviour, such as, for example, parameters that additionally affect the change in strength and cohesion along the depth of the soil, the introduction of tensile strength of the soil and the definition of consolidation parameters.In addition, knowing the soil consistency (very soft, soft, medium, rigid, very rigid) is of significant interest.When defining the parameters according to EN 1997-1:2004 code, it is necessary to know the partial factors for: angle of internal friction, effective cohesion and undrained shear strength of soil [7].The parameters of the actions are defined: types of loads (concentrated, linear, surface, elementi veze kojima se mogu, između ostalog, modelirati i specifični konstruktivni elementi (model ponašanja može biti linearno-elastičan ili nelinearan), podzemna voda (direktno modeliranje horizontalnog ili promenljivog nivoa podzemne vode NPV, modeliranje pornih pritisaka, modeliranje sile uzgona, modeliranje prslina na površini tla ispunjenih vodom i nastalih usled zatezanja) i fazna gradnja/sanacija (modeliranje promene geometrije kosine po fazama gradnje/sanacije, modeliranje promene tla po fazama gradnje/sanacije, modeliranje promene nivoa podzemne vode po fazama gradnje/sanacije, modeliranje promene opterećenja po fazama gradnje/sanacije, modeliranje promene dejstva zemljotresa po fazama gradnje/sanacije, modeliranje promene projektne situacije po fazama gradnje/sanacije).Parametri proračuna umnogome definišu aspekte numeričkih analiza stabilnosti kosina: broj inkremenata kod inkrementalno-iterativne analize, broj iteracija kod inkrementalno-iterativne analize, broj korekcija matrice krutosti sistema, vrednosti tolerancija (za pomeranje, neizbalansirane/rezidualne sile i energiju) i faktor optimizacije (algoritam pretraživanja minimalnog faktora sigurnosti za veći broj kliznih površi).

Podela metoda proračuna stabilnosti klizišta
Metode proračuna stabilnosti klizišta generalno se mogu podeliti u četiri grupe: analitičke, numeričke, eksperimentalne i hibridne.U zavisnosti od toga koja će metoda biti primenjena, dobijaju se rešenja s manjim ili veći stepenom pouzdanosti, s tim što prednost treba dati numeričkim metodama.S obzirom na to što se analitičke i numeričke metode proračuna stabilnosti klizišta najviše primenjuju pri projektovanju i sanaciji klizišta, ali i za potrebe naučnih istraživanja, pregled istraživanjaprikazan u daljem tekstu rada -odnosi se samo na ove metode.U zavisnosti od načina dobijanja konačnog rešenja ispitivanja stabilnosti klizišta, moguće je sprovesti podelu na metode kojima se rešenje dobija putem jednog koraka ili jednokoračne analize (one step), putem više koraka ili višekoračne analize (step by step) i inkrementalno-iterativne nelinearne analize.Shodno prethodno definisanom, uvedena je podela na metode proračuna klizišta: spatial), types of load action (permanent, occasional, incidental), seismic effect (through seismic coefficients, where the response of the landslide is considered in the capacity domain or through the accelerogram, where the response of the landslide is considered in the time domain) and the project situation (permanent, occasional, incidental, seismic).Specific types of parameters are modelled: contour conditions (only stiffness or stiffness and damping components), interface zone, rigid bodies (they do not take into account the effects of their deformations, but only displacements in overall system displacements), specific friction on construction-ground relation (for structures located on or near the landslide or structures that are used for landslide repair), link elements which can be used to model specific structural elements (the behaviour model can be linear-elastic or non-linear), groundwater (direct modelling of horizontal or variable level of groundwater NPV, modelling of the pore stress, modelling of the lifting force, modelling of cracks on the surface of the soil filled with water and caused by tensioning) and phase construction/repair (modelling the slope geometry change by construction/repair phases, modelling the soil change by construction/repair phases, modelling the groundwater level change by construction/repair phases, modelling the load change by construction/repair phases, modelling the change of the earthquake effects by construction/repair phases, modelling the change in the project situation by construction/repair phases).The calculated parameters define, as far as possible, the numerical analysis of the slope stability: number of increments in the incrementaliterative analysis, number of iterations in the incremental-iterative analysis, number of corrections of the system stiffness matrix, tolerance values (for displacement, unbalanced/residual forces and energy) and optimization factor (algorithm for search of the minimal safety factor, for a greater number of sliding surfaces).

Methods of landslide stability calculation division
Methods of landslide stability calculation can, generally, be divided into four groups: analytical, numerical, experimental and hybrid.Depending on the method applied, solutions with a lower or a higher degree of reliability are obtained, while the priority should be given to numerical methods.Since the analytical and numerical methods of landslide stability calculation are mostly applied in the design and repair of landslides, but also for the needs of scientific researches, this exactly is why the overview of the researches, presented in the following text, applies only to these methods.Depending on the method of obtaining the final solution of the landslide stability test, it is possible to divide the methods according to whether the solution is obtained through one step or one-step analyses, through several steps or step-by-step analyses, and incrementally-iterative nonlinear analyses.According to the previously defined, a division of landslide calculation methods was introduced:  analitičke jednokoračne;  analitičke višekoračne (iteracije kliznih površi);  numeričke višekoračne (iteracije kliznih površi);  numeričke inkrementalno-iterativne (nelinearne) analize;  numeričke inkrementalno-iterativne (nelinearne) analize, s primenjivanjem numeričke integracije u vremenskom domenu.
 analytical one-step,  analytical step-by-step (iterations of sliding surfaces),  numerical step-by-step (iterations of sliding surfaces),  numerical incremental-iterative (nonlinear) analysis,  numerical incremental-iterative (nonlinear) analysis, by applying numerical integration in the time domain.

Analytical methods of landslide stability calculation
The key factor in landslide analysis is landslide stability calculation, so as to identify whether the landslide is in the state of balance, whether there is a risk of its losing the balance, or if it is not in the state of balance.In general, with analytical landslide stability methods, the ground is divided into vertical blocks, and for each block corresponding forces are determined, whereby the sliding surface can be circular or polygonal.Depending on the mathematical model of the calculation of the forces acting between the blocks and the shapes of the blocks, there are many analytical methods developed, and the methods of landslide stability to which the practice and the science became accustomed with are those according to: Sarma, Spencer, Janbu, Morgenstern-Price, Shahunyants, Bishop, Fellenius/Petterson and the like.Figure 1 gives a schematic presentation of the ground division into blocks for general analysis of slope stability with a polygonal and a circular slide planes.The corresponding forces for all the blocks are: n normal forces Niacting on each individual block, ns hear forces Ti which act on the edge of the slide plane of each individual block, n-1 normal forces Ei acting between the blocks, n-1 shear forces Xi which act between the blocks, n-1 geometric places zi acted on by Ei forces and n the geometric places li where forces Ni act.In total, 6n-2 unknowns which should be determined from 4n equations (equilibrium conditions).It is obvious that 2n-2 unknowns have to be either approximated or predetermined.Sarma-ina metoda zasniva se na podeli tla na blokove koji nisu strogo vertikalni, već imaju određeni ugao zakošenja, pri čemu su Ei i Xi normalne i smičuće sile između blokova, Ni i Ti -normalne i smičuće sile koje deluju po ivici klizne površi svakog pojedinačnog bloka, Wi -sopstvena težina bloka, KhWi -horizontalna sila kojom se obezbeđuje postizanje graničnog stanja [28].Kh faktor predstavlja odnos horizontalnih i gravitacionih ubrzanja.Na slici 2 prikazana je podela tla na blokove za analizu stabilnosti kosine prema Sarmainoj metodi.
The Sarma method is based on the division of the ground into blocks that are not strictly vertical, but rather have a certain inclination angle, where Ei and Xi are normal and shear forces between the blocks, Ni and Ti normal and shear forces acting on the edge of the sliding surface of each individual block, Wi the block's self weight, KhWi horizontal force which ensures reaching the limit state [28].The Kh factor represents the ratio of horizontal and gravitational accelerations.Figure 2 shows the division of the ground into blocks for the analysis of slope stability according to the Sarma method.
where Fx,i and Fy,i are components of the horizontal and vertical forces projections, rx,i and ry,i arms of the forces Fx,i and Fy,i, respectively, of the PWi resultant of the force of the pore pressure to the divided blocks, The Spencer method is based on the limit equilibrium of the slope, by reaching the balance of forces and moments acting on individual blocks [30].Figure 3 shows the division of the ground into blocks for slope stability analysis according to the Spencer method.Slika 3. Podela tla na blokove za analizu stabilnosti kosine prema Spencer-ovoj metodi [30] Figure 3. Division of the ground into blocks for slope stability analysis according to the Spencer method [30] S ciljem postizanja rešenja problema granične ravnoteže kosine, koja je podeljena na blokove, uvedene su određene pretpostavke: ravni -kojima su podeljeni blokovi -ostaju vertikalne i tokom proračuna, linija dejstva sopstvene težine bloka Wi prolazi kroz centar itog segmenta klizne površi i predstavlja se tačkom M, normalna sila Ni deluje u centru i-tog segmenta klizne površi u tački M i ugao dejstva sile Ei, koja deluje između blokova, jeste konstantan za sve blokove i jednak je δ.Algoritam proračuna stabilnosti kosine prema Spencerovoj metodi zasniva se na izrazima: In order to achieve a solution to the problem of the limit equilibrium of the slope, which is divided into blocks, certain assumptions have been made: the planes, which divide the blocks, remain vertical during the calculations as well, the line of action of the block's self weight Wi passes through the centre of the i-th segment of the sliding surface and it's represented as the point M, the normal force Ni acts in the centre of the i-th segment of the slide plane at the point M and the angle of action of the force Ei, which acts between the blocks, is constant for all the blocks and equals δ.The algorithm of the slope stability calculation according to the Spencer method is based on the following expressions: gde je Ui rezultanta pornog pritiska na za i-ti segment klizne površi, M1i -momenat sila Fx i Fy oko tačke M. Izraz (5) predstavlja relaciju između efektivne i totalne vrednosti normalnih sila koje deluju duž klizne površi.Izraz (6) predstavlja relaciju između normalnih i smičućih sila segmenta klizne površi (Mohr-Coulomb-ovi uslovi).Preformulacijom izraza (7) i ( 8) dobija se: where Ui is the resultant of the pore pressure for the i-th segment of the slide plane, M1i is the moment of forces Fx and Fy around the point M. The expression (5) represents the relation between the effective and the total value of the normal forces acting along the sliding surface.The expression (6) represents the relation between the normal and shear forces of the sliding surface segment (Mohr-Coulomb conditions).By reformulating the expressions ( 7) and (8), we get: Primenom izraza (10) mogu se odrediti sve sile Ei koje deluju između blokova za date vrednosti δi i Fs.Preformulacijom izraza (9), dobija se: By applying the expression (10), all the forces Ei acting between the blocks for the given values δi and Fs can be determined.By reformulating expression (9) we get: Primenom izraza (11), mogu se odrediti svi kraci sile z za date vrednosti ugla δi.Faktor sigurnosti Fs određuje se primenom iterativnog algoritma: inicijalna vrednost za ugao δ jeste δ=0, faktor sigurnosti Fs, za datu vrednost ugla δ, određuje se prema izrazu (10), imajući u vidu to što je En+1=0 na kraju klizne površi, ugao δ se određuje iz izraza (11), koristeći vrednosti za silu E -koja je određena iz prethodnog koraka analize, pri čemu je vrednost zn+1=0 i prethodna dva koraka analize iterativno se ponavljaju sve dok vrednost ugla δ, u dve uzastopne iteracije, ne postane jednaka.Da bi algoritam iteracija bio dovoljno stabilan, potrebno je intervenisati kako bi se otklonila nestabilna rešenja.Ove nestabilnosti javljaju se kada se u izrazima (10) i ( 11) pojave situacije deljenja s nulom.U izrazu (11) ovakva situacija može se pojaviti za vrednosti ugla δ=π/2 ili δ=-π/2, pa se rešenje mora tražiti za interval ugla δ=[-π/2;π/2].Deljenje s nulom u izrazu (10) pojavljuje se u slučaju: By applying the expression (11) all the moment arms of the force z for the given values of the angle δi can be determined.The safety factor Fsis determined using an iterative algorithm: the initial value for the angle δ is δ=0, the safety factor Fs for the given value of the angle δ is determined according to the expression (10), taking into account that En+1=0 at the end of the sliding surface, the angle δ is determined from the expression (11) using the values for the force E, which is determined from the previous step of the analysis, where the value zn+1=0 and the previous two steps of the analysis are repeated iteratively until the value of the angle δ, during two consecutive iterations, becomes equal.In order for the iteration algorithm to be stable enough, it is necessary to intervene with the aim of eliminating any unstable solutions.These instabilities occur when expressions (10) and (11) show the situation of the division by zero.In expression (11) such a situation can occur for the values of the angle δ=π/2 or δ=-π/2, so the solution should be sought for the interval of the angle δ=[π/2;π/2].Division by zero in expression (10) appears in the case of: Radi sprečavanja nestabilnosti rešenja, potrebno je sprovesti proveru parametra mα prema izrazu: In order to prevent the solution instabilities, it is necessary to perform a parameter check mα according to the expression: Pre nego što se započne sa iterativnom analizom, potrebno je pronaći najveću kritičnu vrednost Fs,min koja zadovoljava prethodne uslove.Vrednosti faktora sigurnosti Fs koje su ispod ove kritične vrednosti Fs,min pripadaju oblasti nestabilnog rešenja.Prva iteracija započinje s vrednošću faktora sigurnosti Fs koja je tek nešto veća od Fs,min, tako da su i preostale vrednosti faktora sigurnosti Fs -koje se određuju proračunomuvek veće od Fs,min.
Janbu-ova metoda jeste procedura verifikacije Before beginning iterative analysis, it is necessary to find the highest critical value of Fs,min that satisfies the previous conditions.The values of the safety factors Fs below this critical value Fs,min belong to the area of unstable solution.The first iteration starts with the value of the safety factor Fs, which is just slightly higher than Fs,min, so the remaining values of the safety factors Fs, which are determined by the calculation, are always higher than Fs,min.
The Janbu's method is a procedure of verifying the stabilnosti granične ravnoteže kosina, a zasniva se na uspostavljanju ravnoteže sila i momenata koji deluju na pojedine blokove [19].Na slici 4 prikazana je podela tla na blokove za analizu stabilnosti kosine prema Janbuovoj metodi.stability of the slopes' limit equilibrium, and it is based on establishing the balance between forces and moments acting on individual blocks [19].Figure 4 shows the division of the ground into blocks for slope stability analysis according to Janbu's method.
Slika 4. Podela tla na blokove za analizu stabilnosti kosine prema Janbu-ovoj metodi [19] Figure 4. Division of the ground into blocks for slope stability analysis according to the Janbu's method [19] Radi postizanja rešenja problema granične ravnoteže kosine koja je podeljena na blokove, uvedene su određene pretpostavke: ravni -kojima su podeljeni blokovi -ostaju vertikalne i tokom proračuna, linija dejstva sopstvene težine bloka Wi prolazi kroz centar itog segmenta klizne površi i predstavlja se tačkom M, normalna sila Ni deluje u centru i-tog segmenta klizne površi u tački M i vertikalna pozicija zi dejstva sile Ei, koja deluje između blokova, jednaka je nuli za krajnje tačke klizne površi.Izbor vertikalne pozicije zi dejstva sile Ei ima značajan uticaj na dobijanje konvergentnog rešenja.Ukoliko se loše pretpostave vertikalne pozicije zi, može nastupiti divergencija rešenja, uz prethodno znatno povećanje vremena proračuna.Vertikalne pozicije zi dejstava sila Ei usvajaju se da su jednaki trećini visine blokova na koje je podeljena kosina.Ukoliko nastupi divergencija rešenja, potrebno je korigovati vrednosti zi, tako što se one blago povećavaju kod blokova pasivne zone (kod nožice kosine) i blago smanjuju kod blokova aktivne zone (kod vrha kosine).Algoritam proračuna stabilnosti kosine prema Janbu-ovoj metodi zasniva se na izrazima: In order to reach a solution to the problem of the limit equilibrium of the slope, which is divided into blocks, certain assumptions have been made: the planes which divide the blocks, remain vertical during the calculation as well, the line of action of the block's self weight Wi passes through the centre of the i-th segment of the sliding surface at the point M, the normal force Ni acts in the centre of the i-th segment of the slide plane at the point M and the vertical position zi of the action of the force Ei, which acts between the blocks, is equal to zero for the end points of the sliding surface.The choice of the vertical position zi of the effect of the force Ei has a significant influence on obtaining a convergent solution.If the vertical positions of zi are inaccurately assumed, divergence of the solution can occur, with a significant increase in the calculation time.Vertical positions of zi action of the forces Ei are assumed to be equal 1/3 of the blocks height, to which the slopes are divided.If there a divergence of the solution occurs, it is necessary to correct the zi values, by slightly increasing them with the passive zone blocks (at the foot of the slope) and slightly decreasing them with the blocks of the active zone (at the top of the slope).The algorithm of the slope stability calculation according to the Janbu's method is based on the expressions: Preformulacijom izraza (16) i (17) dobija se: Reformulation of the expressions ( 16) and (17) gives: Preformulacijom izraza (18) dobija se: Reformulation of the expression (18) gives: Faktor sigurnosti Fs određuje se primenom iterativnog algoritma: inicijalne vrednosti svih uglova su δi=0 i pozicije zi su usvojene da su jednake trećini visine blokova, faktor sigurnosti Fs, za datu vrednost ugla δ, određuje se prema izrazu (19), uzimajući u obzir da je En+1=0 na kraju klizne površi, ugao δ se određuje iz izraza (20) koristeći vrednosti za silu E, koja je određena iz prethodnog koraka analize i prethodna dva koraka analize iterativno se ponavljaju, sve dok vrednost ugla δ u dve uzastopne iteracije ne postane jednaka.Otklanjanje nestabilnih rešenja sprovodi se isto kao i u slučaju Spencer-ove metode.
The safety factor Fsis determined using an iterative algorithm: the initial values of all angles are δi=0 and the positions zi are assumed to be equal to 1/3 of the blocks' height, the safety factor Fs for the given angle δ value, is determined according to the expression (19), taking into account that En+1=0 at the end of the sliding surface, the angle δ is determined from the expression (20) using the values for the force E, which is determined from the previous step of the analysis, and the previous two steps of the analysis are iteratively repeated until the value of the angle δ in two consecutive iterations is equal.Removing any unstable solutions is conducted in the same way as with Spencer's method.
In order to reach a solution to the problem of the limit equilibrium of a slope, which is divided into blocks, certain assumptions (similar to the Spencer's method) have been made: the planes, which divide the blocks, remain vertical during the calculations as well, the line of action of the block's self weight Wi passes through the centre of the i-th segment of the sliding surface and it's represented as the point M, the normal force Ni acts in the centre of the i-th segment of the sliding surface at the point M, and the angle of action of the force Ei (acting between the blocks) is different for all the blocks and equals δ=0 for the end points.The assumption of the value of the angle δi is established by using the halfsine function.Figure 6 shows a spectrum of half-sine functions.The choice of the form of the function has less influence on the quality of the final solution, but with the choice of an appropriate form of the function, contributes to the convergence of the solution.The angle δi is determined by multiplying the value of the half-sine function f(xi) and the parameter λ.Slika 6. Polusinusna funkcija za pretpostavke o vrednosti ugla δi [26] Figure 6.A half-sine function for assumptions about the value of the angle δi [26] Algoritam proračuna stabilnosti kosine, prema Morgenstern-Price-ovoj metodi, zasniva se na izrazima koji su identični izrazima (5÷11) kod Spencer-ove metode.Faktor sigurnosti Fs određuje se primenom iterativnog algoritma: inicijalna vrednost uglova δi je δi=λf(xi), faktor sigurnosti Fs, za datu vrednost ugla δ, određuje se prema izrazu (10), uzimajući u obzir da je En+1=0 na kraju klizne površi,ugao δ se određuje iz izraza (11) koristeći vrednosti za silu E, koja je određena iz prethodnog koraka analize (zn+1=0), pri čemu se vrednost polusinusne funkcije f(xi) zadržava kao konstantna kroz iteracije, a iterira se parametar λ i prethodna dva koraka analize iterativno se ponavljaju sve dok vrednost ugla δ u dve uzastopne iteracije ne postane jednaka.Kako bi se sprečila numerička nestabilnost rešenja, sprovode se kontrole prema izrazima (12) i ( 13).
The algorithm of the slope stability calculation according to the Morgenstern-Price's method is based on the expressions that are identical to expressions (5÷11) in the Spencer's method.The safety factor Fs is determined by using an iterative algorithm: the initial value of the angles δi is δi=λf(xi), the safety factor Fs for the given value of the angle δ is determined according to the expression (10), taking into account that En+1=0 is at the end of the sliding surface, the angle δ is determined from the expression (11) using the values for the force E, which is determined from the previous step of the analysis (zn+1=0), while the value of the half-sine function f(xi) is kept constant through iterations, and the parameter λ is iterated and the previous two steps of the analysis are iteratively repeated until the value of the angle δ is equal in two consecutive iterations.In order to prevent the numerical instability of the solution, controls are conducted according to the expressions ( 12) and (13).
The Shahunyants's method for verifying the stability of the limit equilibrium of slopes is based on a similar principle as the previous methods [31].Figure 7 shows the division of the ground into blocks for slope stability analysis according to the Shahunyants's method.In order to reach a solution to the problem of the limit equilibrium of the slope, which is divided into blocks, certain assumptions have been made: the planes, which divide the blocks, remain vertical during the calculation, and the angle of action of the force Ei, acting between the blocks, equals zero (the forces act horizontally).Slika 7. Podela tla na blokove za analizu stabilnosti kosine prema Shahunyants-ovoj metodi [31] Figure 7. Division of the ground into blocks for slope stability analysis according to the Shahunyants's method [31] Algoritam proračuna stabilnosti kosine prema Shahunyants-ovoj metodi započinje transformacijom sila Px,i i Py,i u pravcu normale (N) i tangente (T) klizne površi: The algorithm of the slope stability calculation according to the Shahunyants's method begins with the transformation of the forces Px,i and Py,i in the direction of the normal (N) and the tangent (T) of the sliding surface: Sile koje deluju duž segmenata klizne površi proračunavaju se prema: The forces acting along the sliding surface segments are calculated according to: Jednačina ravnoteže upravno na ravan segmenta klizne površi glasi: The equation of equilibrium perpendicular to the plane of the sliding surface segment is: dok jednačina ravnoteže u ravni segmenta klizne površi glasi: while the equation of equilibrium in the plane of the sliding surface segment is: Uvođenjem izraza ( 23) u (25) dobija se: By introducing the expression ( 23) into (25), we get: dok se uvođenjem izraza (24) u (26) dobija: whereas, by introducing the expression ( 24) into (26), we get: Nakon sređivanja izraza (27), dobija se: After arranging the expression (27), we get: i.e.: S obzirom na to što je: Taking into the account the following: dobija se da je izraz (29): we get that the expression ( 29) is: a dodatnom modifikacijom izraza (31) dobija se: and the additional modification of the expression (31) gives: Primenom izraza (32) sile koje deluju između blokova Ei određuju se prema: By applying the expression (32), the forces action between the blocks Ei are determined according to: Sada se u proračun stabilnosti kosine uvodi faktor sigurnosti Fs, dok se PQ,i sile razlažu na sile koje doprinose klizanju PQ,i,sd (aktivne sile) i sile koje ne doprinose klizanju PQ,i,ud (stabilizujuće sile): Now, the safety factor Fs is introduced into the slope stability calculation, while the PQ,i forces are broken down into the forces contributing to the sliding PQ,i,sd (active forces) and the forces that do not contribute to sliding PQ,i,ud (stabilizing forces): PQ,i je pozitivno kada doprinosi klizanju kosine, a negativno kada ne doprinosi klizanju kosine, tako da se izraz (34) može pisati u formi: PQ,i is positive when it contributes to the sliding of the slope, and negative when it does not contribute to the sliding of the slope, hence, the expression (34) can be written in the form: Na kliznoj površi vrednost sile E0 jednaka je nula, dok za E1 važi: On the sliding surface, the value of the force E0 equals zero, whereas the following applies to E1: a za E2: and to E2: Slično se mogu prikazati i izrazi za sve sile koje deluju između blokova, pri čemu je En=0: The expressions for all the forces acting between the blocks can be presented in a similar way, where En=0: tako da se iz ovog izraza može direktno prikazati faktor sigurnosti Fs u formi: so that, from this expression, the safety factor Fs can be directly presented in the following form: Faktor sigurnosti Fs prema Fellenius/Petterson-ovoj metodi određuje se na osnovu izraza: The safety factor Fs, according to the Fellenius/Petterson's method, is determined on the basis of the expression: dok se prema Bishop-ovoj metodi određuje na osnovu izraza: whereas, according to Bishop's method, it is determined on the basis of the expression: (41)

Numerical methods of landslide stability calculations
Landslide stability calculation using numerical methods is based on methods of domain discretization, such as: In these methods, the soil is considered as a linearelastic, elasto-plastic and non-linear material.The Finite Element Method (FEM) is mostly used for solving the problem of numerical slope stability analysis, so a large number of software has implemented algorithms based on this method.Figure 8 shows the mesh of finite elements of the discrete numerical model of the slope and the set of points obtained by optimizing the slope safety factor according to the Finite Element Method (FEM).The slope is modelled by using surface finite elements with an integrated mathematical formulation for the analysis of the plane strain.When modelling and analyzing slope stability, two important aspects need to be taken into account: discretization and approximation.Discretization refers to the problem of the ground domain division into finite elements of sufficiently small dimensions for which the criteria of the relation between the diagonal and the angles of the quadrangle finite element or the relations of the sides of the triangle finite element must be respected.In the area of the contact Diskretizacija se odnosi na problem podele domena tla na konačne elemente dovoljno malih dimenzija za koje se moraju poštovati kriterijumi odnosa dijagonala i uglova četvorougaonog konačnog elementa ili odnosi stranica traouganog konačnog elementa.U oblasti kontakta tla sa elementima za plitko ili duboko fundiranje, koji se koriste prilikom sanacije klizišta, potrebno je izvršiti progušćenje mreže konačnih elemenata.Takođe, progušćenje se sprovodi i u zoni klizne površi, na mestima diskontinuiteta i otvora u tlu i slično.
Establishing a connection between the basic finite elements, which form the domain of the ground, with the increased density mesh of finite elements is carried out by using transition elements.As transition elements, the most commonly used are triangular finite elements.A very important aspect is also establishing the compatibility of finite elements nodes through conformity/nonconformity analysis, especially with transition finite elements, whereby it should not be allowed for certain nodes, in combination of basic and transition finite elements, to be left unconnected or partially connected.Figure 9 shows the 2D numerical slope models with generated finite element mesh and increased density over selected domains.
Compared to the 2D model of the slope, which is the one mostly used for practical purposes, by using the 3D model of the slope, more complex geometric models can be modelled, with a spatially more complex and significantly more variable geology in a smaller area.Figure 10 shows the 2D and 3D numerical models of the slope with a separate representation of the sliding mass of the soil and the spatial model of the sliding surface.For the modelling of the 3D model of slopes, solid or tetrahedral finite elements are used, whereby modelling of the ground domain by spatial finite elements requires significantly higher hardware capacities.With prismatic finite elements, a minimum of 2x2x2 numerical integration is applied over Gaussian quadratures [8].
In certain cases, when the ground domain is of considerable dimensions and a slightly complex geometry, the finite elements mesh of the 3D model of the slope can even have a several million finite elements, so in these cases, the most commonly used is parallel processing technique.With this type of problems, the mesh of finite elements and the numbering of the nodes of the elements are additionally optimized, since optimizing the numbering of finite element nodes reduces the bandwidth of the system stiffness matrix and concentrates the members of the system stiffness matrix around the diagonal.Figure 13 shows 3D numerical slopes models of a slightly complex geometry with the separate sliding mass of soil.Modelling the sliding surface, when analyzing the stability of 3D slopes models, can be carried out, as it has already been presented, by using 3D spatial finite elements or even 2D surface finite elements.
Figure 14 shows the 3D numerical models of the slopes of a slightly complex geometry with separately shown sliding mass of the soil and the locations of the calculated points of the safety factors, obtained through optimization for the concave and convex sliding surface.The concave sliding surface is formed from 3D spatial finite elements, while the convex sliding surface is formed by combining 3D spatial and 2D surface finite elements.
Proširena metoda konačnih elemenata (XFEM), za razliku od metode konačnih elemenata (FEM), ima mogućnost primene poboljšane nelinearne analize i proračuna postnelinearnog ponašanja sistema.Takođe, kod ove metode, prilikom formiranja klizišta, može se modelirati razvoj: prslina, pukotina i raseda u tlu.Prsline u tlu, u opštem slučaju, modeliraju se kao razmazane, dok se kod visokozahtevnih problema formiranja klizišta primenjuju algoritmi modeliranja diskretnih prslina.Model diskretnih prslina u tlu zahteva implementaciju algoritama mehanike kontakta, dok se model razmazanih prslina u tlu rešava nelinearnom analizom trajektorija ekstremnih vrednosti glavnih napona u tlu.Metoda graničnih elemenata (BEM) ima značajnu primenu u geotehnici, budući da se primenom ove metode brže dobijaju rešenja, u odnosu na metodu konačnih elemenata (FEM), pri čemu je i nivo kvaliteta konačnog rešenja zadovoljavajući.S obzirom na to što postoji nekoliko algoritama u okviru metode graničnih elemenata (BEM), oni se -u najvećem broju slučajeva -Modelling of the softening and discontinuity in the soil is carried out by correcting the parameters of the constitutive model of soil behaviour and eliminating the connection of certain finite elements or even the reducing of a number of finite elements, which are located in a special zone of refined finite element mesh.The approximation refers to the choice of the optimal type of the finite element through which the field of soil displacement in the slope model is effectively modelled.In this case, there is a number of developed finite elements types in which unknowns are determined by: force, displacement or combined (mixed).For interpolation functions, an isoparametric formulation is used, while the nodes for the numerical integration calculation are mapped: in the angles, in the interior and/or on the contour of the final element.Also, the aspect of approximation refers to: numerical modelling of contour and transition conditions, modelling of material behaviour and modelling of effects -loads.
The eXtended Finite Element Method (XFEM), compared to the Finite Element Method (FEM), offers the possibility of applying an improved nonlinear analysis and the post-non-linear system behaviour calculation.Also, with this method, during the formation of the landslide, it is possible to modelled the development of: cracks, gaps and splits in the soil.In general, cracks in the soil are modelled as smeared, while with the highly demanding problems of landslide formation, the modelling algorithms for discrete cracks are applied.The model of discrete cracks in the ground requires the implementation of algorithms of contact mechanics, while the model of smeared cracks in the soil is solved by nonlinear analysis of the main stress in the soil for extreme values trajectory.The Boundary Elements Method (BEM) has a significant application in geotechnics, since the application of this method gives solutions faster than the Finite Elements Method (FEM), while the quality of the final solution is also satisfactory.zasnivaju na diskretizaciji granične oblasti (kontura) graničnim elementima.Unutrašnjost oblasti najčešće se ne diskretizuje, pa ovakve metode pripadaju grupi bezmrežnih metoda.Metoda diskretnih elemenata (DEM) zasniva se na razmatranju ravnotežnog stanja pojedinačno za svaki konačni element.U poređenju s metodom konačnih elemenata (FEM), gde se ravnotežno stanje razmatra na globalnom nivou preko kompletnog numeričkog modela, kod metode diskretnih elemenata (DEM) jednačine kretanja definišu se posebno za svaki konačni element, tako da se mogu pratiti međusobno nezavisna polja pomeranja konačnih elemenata.Na slici 15 prikazan je 2D numerički model kosine prema metodi diskretnih elemenata (DEM) sa identifikovanom zonom iniciranja klizišta.
Since there are several algorithms within the Boundary Elements Method (BEM), they are mostly based on the discretization of the boundary area (contours) by the boundary elements.In most cases, the intrinsic domain is not discretized, so such methods belong to the group of mesh free methods.The Discrete Element Method (DEM) is based on the analysis of the equilibrium state for each finite element individually.In comparison to the Finite Element Method (FEM), where the equilibrium state is considered globally, through a complete numerical model, with the Discrete Element Method (DEM), the motion equations are defined for each finite element individually, so that the independent fields of finite elements movement can be traced.Primenom ove metode, može se pratiti inkrementalni razvoj klizišta, tako da se kao konačna vrednost proračuna dobija spektar faktora sigurnosti.Takođe, ova metoda primenjuje se i za 3D modeliranje složenih formi kosina, pri čemu je razvijen niz algoritama za topologiju i kompaktnost elementa kojima se formira 3D model kosine.Na slici 16 prikazan je postupak formiranja 3D numeričkog modela kosine prema metodi diskretnih elemenata (DEM) i odgovarajuće inkrementalne proračunske faze.
For this algorithm to be effectively applied in practice, the connections between the finite elements are modelled by the contact elements with the possibility of implementing different nonlinear behaviours.The contact elements define the stiffness components under the pressure, and the tensile stresses can also be defined or even eliminated.When modelling the contact between two points of the model, two states occur: active (the contact is established with the involvement of certain stiffness) and inactive (the contact is not established with the involvement of little stiffness or without the introduction of stiffness effects).In order to efficiently model the effects of contact elements interaction, it is necessary to apply the geometric nonlinear incrementaliterative analysis.Due to the non-linear behaviour of the contact element, where the change of the state can be followed by a major change in stiffness, serious difficulties can arise in ensuring the convergence of the kontrole inkrementalnog priraštaja silama.
nonlinear solution.In that sense, it may be more beneficial to use the procedure for controlling the incremental increase of displacements, rather than the procedure for controlling the incremental increase of forces.
Numeričke inkrementalno-iterativne (nelinearne) analize stabilnosti klizišta zasnivaju se na formulaciji nelinearnog problema sistemom nelinearnih algebarskih jednačina oblika [2], [5]: The introduction of the contact mechanics in the analysis of the development of large plastic deformations and the displacement of the landslide soil mass is also carried out with the eXtended Finite Element Method (XFEM), similar to the Discrete Element Method (DEM).In the formulation of the problem itself, it is considered that during incremental calculation phases occurs such a change in the geometry of the contact zone, that the initial generated mesh of finite elements is corresponding to the configuration of the mesh of finite elements for any incremental situation.This eliminates the use of additional algorithms for search for a favourable configuration in connecting the mesh nodes in i-th incremental analysis, or even the use of an adaptive method for correcting the mesh of finite elements of the system [34].
Numerical incremental-iterative (nonlinear) landslide stability analyses are based on the formulation of a nonlinear problem through a system of non-linear algebraic equations of the form [2], [5]: i.e.: gde su {u} nepoznati parametri pomeranja, {F} generalisani spoljašnji uticaji (opterećenja) u čvorovima sistema.Jednačine problema (42) umesto za ukupno opterećenje, rešavaju se za niz posebnih inkrementalnih opterećenja.U okviru svakog inkrementa, pretpostavlja se da je sistem jednačina linearan.Na taj način, rešenje where {u} is the unknown displacement parameters, {F} generalized external effects (loads) in the system nodes.The equations of the problem (42) instead of for the total load, are solved for a series of specific incremental loads.Within each increment, it is assumed that the equation system is linear.In that way, the solution of a nelinearnog problema dobija se kao zbir niza linearnih (inkrementalnih) rešenja.Nelinearan problem može da se prikaže izrazom: nonlinear problem is obtained as the sum of a series of linear (incremental) solutions.A non-linear problem can be represented by the expression: i.e.: gde je {P} vektor unutrašnjih generalisanih sila modela koje su funkcija vektora generalisanih pomeranja {u}, λ parametar inkrementalnog opterećenja (odnos inkrementalnog i kompletnog opterećenja).U skladu s konceptom inkrementalnog rešenja jeste: where {P} is the vector of the internal generalized model forces, which are the function of the generalized displacement vector {u}, {λ} the incremental loading parameter (the ratio of incremental and total load).In accordance with the concept of incremental solution, we have: Iz izraza (46) određuju se inkrementi vektora pomeranja za inkremente opterećenja i tangentnu matricu krutosti modela klizišta, koja se formuliše za referentno stanje na početku inkrementa.Referentnom stanju na početku prvog inkrementa odgovara linearna matrica krutosti klizišta (inicijalna matrica krutosti).Opšti i-ti korak inkrementalnog postupka obuhvata: formiranje tangentne matrice krutosti [Kt]i numeričkog modela klizišta, određivanje inkremenata vektora opterećenja {ΔF}i numeričkog modela, određivanje inkremenata vektora generalisanih pomeranja {Δu}i rešavanjem sistema linearnih algebarskih jednačina za tangentnu matricu krutosti, određivanje inkremenata uticaja u konačnim elementima (deformacije, naponi), i određivanje ukupne vrednosti generalisanih pomeranja inkrementalnim (kumulativnim) sabiranjem.Pomeranja posle m-tog inkrementa određena su izrazom: From the expression (46), the increments of the displacement vector for loading increments of the load and the tangent stiffness matrix of the landslide model stiffness are determined, which is formulated for the reference state at the beginning of the increment.The reference state at the beginning of the first increment corresponds to the linear matrix of the landslide stiffness (initial stiffness matrix).The general i-th step of the incremental procedure includes: the formation of a tangent stiffness matrix [Kt]i of the numerical landslide model, determining the load vector increment {ΔF}i of the numerical model, determining the vector of generalized displacements increments {Δu}i by solving the system of linear algebraic equations for the tangent stiffness matrix, determining the increments of the impact in the finite elements (deformations, tensions), and determining the total value of generalized displacements by incremental (cumulative) addition.Displacements after the m-th increment are defined by the expression: Razlog za pojavu greške inkrementalnog rešenja jeste sprovedena linearizacija u okviru inkrementa.Veličina greške može da se odredi iz uslova ravnoteže na kraju inkrementa.Kao posledica linearizacije, javljaju se neuravnotežena (rezidualna) opterećenja koja su mera odstupanja inkrementalnog rešenja od tačnog.Vektor rezidualnog opterećenja može se prikazati kao odstupanje od ravnoteže: The reason behind the occurrence of the incremental solution error is the linearization conducted within the framework of the increment.The error dimensions can be determined from the balance conditions at the end of the increment.As the linearization consequence, unbalanced (residual) loads occur, that are the measure of deviation of the incremental solution from the exact one.The residual load vector can be represented as a deviation from balance: Korekcija greške postiže se dodavanjem rezidualnog opterećenja na spoljašnje opterećenje u sledećem inkrementu: Error correction is achieved by adding the residual load to the external load in the following increment: Najbolji rezultati postižu se ako se kombinuje inkrementalni i iterativni postupak.U prvoj iteraciji, pojavljuju se rezidualna opterećenja zbog neispunjavanja uslova ravnoteže.Ako se naredne iteracije realizuju samo s rezidualnim opterećenjima, uz korekciju tangentne matrice krutosti, postupak može da konvergira uz minimiziranje rezidualnog opterećenja.Pri formulisanju iterativne metode, polazi se od izraza za razvoj u Taylor-ov red vektora rezidualnih sila u okolini pomeranja {u}j: The best results are achieved if the incremental and iterative processes are combined.In the first iteration, residual loads appear due to unfulfilled balance conditions.If the following iterations are realized only with residual loads, with the correction of the tangent stiffness matrix, the process can converge, with the minimization of the residual load.When formulating the iterative method, it is started with the expression for development in the Taylor series of the residual forces vector in the vicinity of the displacement {u}j: Iz uslova da rezidualno opterećenje ispunjava uslove ravnoteže {R}j+1=0, važi: From the condition that the residual load meets the balance conditions {R}j+1=0, follows: Poslednja dva izraza predstavljaju osnovu iterativne metode.Kombinacijom inkrementalne i iterativne metode dobija se Newton-Raphson-ova inkrementalno-iterativna metoda (slika 17).

Kompleksno 3D geometrijsko modeliranje i numeričke metode proračuna stabilnosti klizišta
Standardni pristup u modeliranju terena i klizištainkorporiranog u terenu -zasniva se na korišćenju tehnike 2D prezentacije primenom situacionog plana i vertikalnih poprečnih preseka.Na osnovu definisanih tipova slojeva tla po dubini i njihovih fizičko-mehaničkih The correction of the system stiffness matrix is carried out after each applied time step, and according to the previously presented Newton-Raphson's method.Using the NDA analysis with the HHT method and the NR method for calculating the 2D and 3D landslide models, the most reliable solutions for estimating the nonlinear system response are obtained.Using this method allows us to consider the influence of the level of underground and surface water increase dynamics, as well as the effect of the earthquake, incrementally scaling the accelerometer.System (landslide) response is represented as the function of change of the safety factor Fs in time, and not only as a unique (discrete) value.
Generalno razmatrajući modeliranje površi u prostoru može se sprovesti primenom matematičkih funkcija, mapiranja i diskretnih vrednosti.Najviše se koristi tehnika mapiranja terena s rasterskom mrežom (ortogonalna, poluortogonalna, radijalna i zakrivljena) za formiranje mape terena, ali je primena diskretnih vrednosti i formiranje polilinija, površi i solida u prednosti, pa se za ovakvu grafiku koristi termin vektorska grafika.Izohipse terena i klizne površi, u opštem slučaju, predstavljaju se primenom polilinija i splajnova.Da bi se geometrijski i matematički modelirao skup tačaka koji formira jednu kliznu površ u 2D koordinatnom sistemu, potrebno je physico-mechanical characteristics, an analytical and/or numerical calculation of the slope stability is carried out.In the case of a spatially slightly complex terrain model and slightly complex landslide geometry, the question of 2D modelling and the reliability of the corresponding analyses can be debatable.However, even in situations where it can be reliably asserted that the 2D presentation of the terrain and landslide by using the situational plan and vertical cross-sections is reliable, the following questions remain open: can the presentation of the terrain and the landslide be further improved in accordance with modern information technology and whether the volume of the soil forming the landslide can be reliably determined?The solution to these issues can be found in 3D visualization of terrain and landslide, whereby the most sophisticated solution, by using 4D visualization (3D + dynamic simulation), can present the problem of landslide sanation, from the initial state, through phase solutions to the final solution.3D modelling of the terrain and landslide is used for geometric presentation and numerical analysis through using surfaces and solids.Geometric 3D presentation, in most cases, has a greater degree of visualization of the final solution, while numerical 3D analysis aims to use the surfaces and solids to model the terrain and landslide, so that each geometric-numerical element also has in itself an integrated mathematical formulation of the problem.It is presumed that the realistic effect of the geometric presentation can be additionally achieved in numerical modelling of the terrain and landslide, however, in these situations the time of the calculation is further increased, so that, in very complex models and with a very large number of finite elements, the calculation is reduced to the application of the parallel processing technique.However, geometric 3D modelling for the presentation of terrain and landslides is much more useful for soil volume calculations, since modelling the landslide as a solid can quickly determine the appropriate volume, even in situations of very complex solid models.The process of complex 3D modelling the terrain and landslide is based on: the previous identification of a larger number of sliding planes for the corresponding number of engineering-geological profiles, integration of these with the 2D situational plan of the landslide and the construction of the landslide 3D model in the geometric presentation software (CAD).For the integrated sliding planes, a sliding surface is formed in space, while for the modelled landslide in space a solid landslide model is formed.The modelling of the sliding surfaces in space is based on the application of complex curved surface that forms a grid of quadrangles, while the solid model of the landslide is generated using primitives and techniques for editing primitives: extrusion, cutting, expanding, unifying, extraction, intersection, and the like.Figure 18 shows the generated complex geometric-numerical 3D terrain models for landslide stability analysis.
graphics.The terrain isohypse and sliding surfaces, in general, are represented using polylines and splines.In order to geometrically and mathematically model the set of points that forms a single sliding surface in the 2D coordinate system, interpolation is required.Those simpler interpolation models are based on the application of mathematical functions in closed form.However, interpolation of a large number of points, that form a single sliding surface in 2D coordinate system, is based on the application of parametric functions, where the solution is not defined in a closed form, but in a set of functions.The connection of these functions is established by the conditions of the tangent equivalence for curves on the left and right at each point of interpolation.This way, a smooth interpolated curve is obtained, so the application of the spline has turned out to be among the best parametric functions.
In the case of 3D terrain model and sliding plane, more precisely the sliding surface, they are modelled in the space using NURBS curves (Non-Uniform Rational Basis Spline).NURBS curves are defined by the control nodes and the node vector.In general, NURBS curves and the corresponding surfaces are the generalization of B-splines and Bezier's curves and surfaces.The control nodes define the shape of the surface, in particular, the sliding surface, while the node vector determines where and how the surface touches the control nodes.However, even with the application of NURBS surfaces, a problem may arise in interpolation, if the adequate interpolation parameters are not selected for certain control nodes, and which are discrete values of a set of sliding planes.Excessive interpolation deviations can occur so that the 3D terrain and landslide model can be approximated in a similar manner as the principle in regression analysis applies, whether they are linear or nonlinear.Minimizing the previous problem is achieved by increase in the density of the mesh of finite elements through the introduction of new inter elements.In general, the most reliable, but also visually rougher solutions are achieved by applying quadrangles whose nodes directly connect discrete nodes (linear interpolation) of the terrain and landslide.The mesh refinement is achieved by interpolation by triangles.As previously mentioned, the presentation of the terrain is carried out, in fact, by using a wireframe plane model with the addition of 3D planes, while the landslide modelling is carried out by using a solid (3D geometric body).The differentiation of the landslide in relation to other parts of the terrain can be carried out by allocation and display of landslide only, irrespective of the terrain, with the possibility of 4D continuous translation and rotation in space, and rendering, so that the terrain transparency is increased in relation to the landslide.Slika 18. Generisani kompleksni geometrijsko-numerički 3D modeli terena za analizu stabilnosti klizišta [17] Figure 18.Generated complex geometric-numerical 3D terrain model for landslide stability analysis [17] 4 ZAVRŠNE NAPOMENE Primenom sprovedene sistematizacije analitičkih i numeričkih metoda proračuna stabilnosti klizišta, može se efikasno razmotriti koji tip metode se može primeniti u fazama preliminarnih i finalnih analiza za naučna istraživanja i stručne projekte.Autori su napravili sopstvenu sistematizaciju metoda proračuna stabilnosti klizišta, s tim što pojedine metode mogu pripadati i prelaznim kategorijama.Posebno je to slučaj kod onih metoda koje se zasnivaju na direktnoj analizi stabilnosti za odgovarajuću kliznu površ i kod metoda koje koriste iteracije kliznih površi primenom optimizacionih algoritama.

FINAL REMARKS
By applying the conducted systematization of analytical and numerical methods of landslide stability calculation it can effectively be considered which type of method can be applied in the phases of preliminary and final analyzes for scientific research and professional projects.The authors have made their own systematization of the methods of landslide stability calculation, but some methods can also belong to transition categories.This is especially the case with those methods that are based on a direct stability analysis for the corresponding sliding surface and for methods using sliding surface iterations by applying optimization algorithms.
Key problems in modelling and numerical analysis of nowadays landslides could be presented through several aspects:  generalization of insufficient number of sampling and obtaining appropriate quality laboratory tests of physical-mechanical characteristics of soil and constitutive models of soil behaviour for a complete landslide,  the application of the geometric-numerical presentation of the landslide through 3D models (in certain situations, higher safety factors can be obtained due to the wedging of the landslide during the soil sliding),  it is necessary to further improve the methodology of landslide stability verification based on mathematical models and analysis of incremental displacement of the landslide, by monitoring the deformations, but not the forces and moments,  implementing parallel processing techniques for practical purposes (increasing: hardware capacities through multi-core processing and resources through storage of the memory).

Kristina BOZIC-TOMIC Nenad SUSIC Mato ULJAREVIC
According to the analysis of a large number of scientific papers, the authors of the paper presented their own original systematization of the analytical and numerical methods of landslide stability calculation, with a large part of them still to be further improved, implemented and tested on complex 3D landslide models.Methods for calculating the stability of the landslide are classified into five groups: analytical singlestep, analytical multi-step (iterations of sliding surfaces), numerical multi-step (iterations of sliding surfaces), numerical incremental-iterative (nonlinear) analysis and numerical incremental-iterative (nonlinear) analysis, applying numerical integration in the time domain.By using the systematization method of calculating the stability of the landslide it can be very effective to consider which type of method is optimal for landslide analysis and which type of method should be considered in the phase of preliminary and final analysis for scientific research and expert projects.

Slika 1 .Figure 1 .
Figure 1.Division of the ground into blocks for general analysis of slope stability: a) polygonal sliding surface, b) circularsliding surface[10] of the internal friction along the sliding surface of the individual blocks, i c the average cohesion value along the sliding surface of individual blocks.The slope safety factor Fs is determined by iteratively reducing the parameters c and tgφ, so as to reach the factor Kh value (zero or greater than zero).

Slika 8 .Figure 8 .
Figure 8. 2D numerical model of a slope: a) a finite elements mesh of the discrete numerical model of a slope accordingto the Finite Element Method (FEM); b) a set of points obtained by optimizing the slope safety factor according to the Finite Element Method (FEM)[12] Figure 15.2D numerical model of the slope: a) numerical model of the slope according to the Discrete Elements Method (DEM), b) identification of the landslide initiation zone according to the Discrete Elements Method (DEM)[22]