Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Simulation of Rock Mass Horizontal Displacements with Usage of Cellular Automata Theory.

Simulation of Rock Mass Horizontal Displacements with Usage of Cellular Automata Theory. coal panel. Obtained results were consistent with results obtained by geometric-integral theory. Keywords: horizontal displacements, rock mass deformation, mining area protection, cellular automata W artykule przedstawiono plaski model górotworu jako deterministyczny, skoczony automat komórkowy. Wykorzystana do opisu rozkladu deformacji wewntrz górotworu i na jego powierzchni teoria automatów pozwala w relatywnie prosty sposób uzyska profil niecki obnieniowej zgodny z profilem obserwowanym pomiarami geodezyjnymi na powierzchni terenu. Przedstawiony w pracy najprostszy model górotworu przedstawia jego plaski przekrój w postaci regularnej siatki komórek, które cile do siebie przylegaj i maj jednakowy ksztalt (Rys. 1). Dla modelu zdefiniowano podstawowe parametry automatu komórkowego takie jak: warunki brzegowe, decydujce o pocztku i kocu symulacji, ssiedztwo komórkowe, okrelajce przestrze w siatce komórek, w obrbie której dochodzi do bezporedniej wymiany informacji zapisanych w poszczególnych komórkach oraz funkcj przejcia, która decyduje o ostatecznej charakterystyce rozkladu symulowanego zjawiska w siatce automatu. W artykule zastosowano deterministyczn funkcj rozkladu. * SILESIAN UNIVERSITY OF TECHNOLOGY, FACULTY OF MINING AND GEOLOGY, INSTITUTE OF MINE SURVEYING AND MINING AREA PROTECTION W wyniku licznych prób modelowych stwierdzono, e stosowana funkcja rozkladu dla symulacji rozkladu obnie (Rys. 3) nie pozwala na symulacj ruchów poziomych jakociowo i ilociowo zgodnych z przemieszczeniami poziomymi obserwowanymi w rzeczywistoci. Jako rozwinicie dotychczasowej koncepcji budowy górotworu jako skoczony automat komórkowy, w pracy opisano funkcj rozkladu (Rys. 5), która pozwala, równolegle do symulacji obnie, symulowa przemieszczenia poziome zgodne z wynikami obserwacji geodezyjnych. Na podstawie wyników licznych symulacji komputerowych opisano podstawow matematyczn zaleno (wzór 11) okrelajc stosunek maksymalnych przemieszcze poziomych do maksymalnych obnie w niecce nadpelnej w odniesieniu do przyjtych parametrów opisanego modelu górotworu, tj.: odwzorowywanych w rzeczywistoci wymiarów komórki, wartoci tzw. przejcia glównego, glbokoci eksploatacji oraz parametru maksymalnego nachylenia (bdcego odpowiednikiem np. parametru tg w teorii Budryka-Knothego). Dla pokazania moliwoci zaproponowanego automatu komórkowego wykonano symulacj rozkladu deformacji wewntrz modelu górotworu i wyznaczono profil linii obnie, nachyle, przemieszcze poziomych i odksztalce poziomych powierzchni modelu dla przykladu abstrakcyjnego, wyeksploatowanego pokladu o okrelonych parametrach górniczo-eksploatacyjnych. W wyniku symulacji otrzymano nadpeln nieck obnieniow o ksztalcie opisywanym przez calk z funkcji Gaussa, w której rozklad przemieszcze poziomych byl zbieny z modelem wzorcowym. Slowa kluczowe: przemieszczenia poziome, deformacje górotworu, ochrona terenów górniczych, automat komórkowy 1. Introduction The Cellular Automata Theory, begun in the 1940s of the last century, found widespread application in many fields of science thanks to its simplicity and exceptional complexity of received results (Wolfram, 2002). The Cellular Automaton (CA) was originally discovered by Stanislav Ulam (Pickover and Clifford, 2009) and John von Neumann (Schiff & Joel, 2011) at the Los Alamos National Laboratory. Cellular automata are commonly used e.g. in computer science, to test the algorithms, in physics to simulate the spread of fire, gas, heat, etc., or to simulate the traffic (Chopard at al., 1996). It turns out that despite the huge potential in the field of rock mechanics and mining area protection, so far, the theory wasn't practically applied. Today, in Poland, predictions of the impact of underground mining on the surface mainly rely on methods which uses geometric-integral theory to describe the distribution of deformation and methods basing on models of continuous centre where the state of stresses and displacements defines a system of differential equations and the equation of state dependent on the adopted model. In previous publications there was proved (Bialek & Sikora, 2012; Sikora, 2010, 2011, 2013, 2014) that using cellular automata to describe the rock mass deformation, results can be very similar, and in some cases even identical with the results obtained with application of different methods, e.g. Knothe's method. So far, it has been shown that for certain conditions the theory of cellular automata allows in a relatively easy way to get subsidence trough profile consistent with the profile observed by geodetic measurements. Moreover it is possible to simulate directly the impact of heterogeneous build of rock mass and nonlinear properties of summing of the mining influences on the subsidence distribution process. Especially there was presented the method of simulation of the impact of faults on the distribution of subsidence. Normal faults, in a form of splits, were implemented on the edges of corresponding cells. Modified distribution function allows for limitation of subsidence propagation trough the split. Nonlinear properties of subsidence summation were implemented conditioning the value of subsidence propagation to the sustained tilt values. In a consequence it is possible to simulate e.g. the impact of multiple extraction on the subsidence distribution. Finally there was presented build concept of the spatial rock mass model as a cellular automaton. In this paper, which is a development of the prior solutions, there will be presented a method of horizontal displacements simulation. In consequence the solution will allow for the characteristics of the horizontal displacements and horizontal deformations within the rock mass and on its surface. The assumptions and the method of implementation were shown on the example of the simplest two-dimensional rock mass model (Sikora, 2011) built in relation to the deterministic finite automaton definition (Moore, 1956). 2. The model of the rock mass as a deterministic finite automaton By cellular automata as a rock mass model can be specified a number of commonly used algorithms, e.g. rock mass model described by T. Niemiec (Niemiec, 1985; Niemiec et al., 2008). However in practise this term is being used in some specific cases and applications. The simplest definition describes cellular automaton as a mathematical model that its structures are determined by the net of cells, theirs states, transitions and rules of transitions (Moore, 1956). In relation to the previous definition there can be build the simplest rock mass model as its flat cross-section presented by the net of discrete cells. Detailed built concept of that model and its variants were described in previous publications (Sikora, 2010, 2011). However there is a group of fundamental parameters describing a model of the rock mass as a deterministic finite automaton, including: the type of cell, the net of cells, so called cell's neighbouring, boundary conditions and a transition function that decide in most about the final characteristics of the simulated phenomena. Subsequently, these parameters are shortly discussed. 2.1. Cell type and the grid of cells Described case of the rock mass cross-section presume its segmentation into regular net of cells of determined number of rows m and columns n. All cells has the same shape of rectangle and the same size and also closely adhere to each other (Fig. 1). Cell's dimensions are corresponding with real dimensions of width Sk [m] and height Wk [m]. Thus formed table (grid of cells) represents rock mass in a flat system, where vertical axis is converging with columns (depth). Individual rows relate to the respective levels of rock mass, however the last row m (the highest) relates to the land surface. 2.2. Initial state In case of presented type of cellular automaton it is assumed that every cell starts in the same state, except for a finite number of cells in other states. It is called configuration (Schiff and Joel, 2011) or initial state. Fig. 1. Cellular automata grid representing horizontal section trough rock mass and extracted longwall (Sikora, 2013) In case of considered rock mass model, as a initial state, can be assumed two-dimensional grid of cells. For every single cell can be assigned a value representing the void. At the beginning values assigned to proper cells (initial subsidence) shows performed mining operations (equation 1). pi, j xm,n = agp where: pi, j xm,n wm,n gp a -- -- -- -- -- (1) elementary particle of the seam with coordinates i, j, a single cell in the grid of automaton with coordinates m,n, initial subsidence in the current cell with coordinates m,n, the average thickness of the seam, subsidence factor ­ the theory's parameter that shows the difference between the volume of extracted seam and the subsidence trough on the surface (subsidence at the surface depends on the strength and nature of the rocks overlying the seam). The above equation means that for every particular element of the seam pi, jj (i, j are its real coordinates) corresponds at least one cell xm,n in the CA's grid. For that cell it's being ascribed the value of extraction height (of the seam) gp multiplied by the value of subsidence factor a. The final state is the moment when the sum of void's volume in the row representing the surface of the mining area equals the sum of void's volume assigned to specified cells at the initial state. It should be noted that the initial conditions does not affect the cell shape. Theoretically cells are considered to be shapeless and it's not required to visualize them. There will be analysed only the effect of their work. 2.3. Cell's neighborhood As a neighborhood of the individual cell meant all cells that are directly related by a dependence function (Schiff and Joel, 2011). NW W SW NE E SE Fig. 2. Moore type neighborhood defined for the dark cell in the regular net of Cellular Automaton (own source) In the rock mass as a closest neighbourhood of a post-operational void can be understood the immediate surrounding area (cave zone). In the CA theory every single cell has the same neighbourhood. One of the most popular type of neighbourhood is the Moore's type (Fig. 2) which include all adjacent cells to the current cell. Formally this is the best type of neighbourhood that could be considered in case of mapped void in the rock mass model. The force of gravity causes the formed empty space in the rock mass (void) is being filled with overhanging rocks. Taking this assumption and the lack of existence of any additional forces affecting tightening of the void, there can be finally determined the cell's neighborhood. Therefore it's limited to three adjacent cells in the row above the main cell (cells NW, N and NE in the Fig. 2). In other words every cell will interact only with neighbor cells. 2.4. Propagation function in order to the rock mass subsidence simulation The manner of operation, compatibility features of the model with the real object and the state of individual elementary cell in a discrete time decide first and foremost the evolution rule of cellular automaton. The evolution of the model base on the boundary conditions defined for the automaton and the propagation function that determine intercellular relations in the space defined by the cell's neighborhood. In case of presented deterministic model (there is also stochastic model (Sikora, 2011)) fulfilment of the post-operational void rely on determination of partial distribution of the void's volume outflow in the range of cell's neighborhood. Presented algorithm in the Fig. 3 assumes that the partial value of a void is filled by the rocks volume from a cell placed directly above the void. It's so called main transition P. Other transitions take place symmetrically with respect to the main transition to cells from the neighborhood (on the left and right). For values of P = 0,5 presented characteristics of the propagation function leads to the subsidence distribution (on the surface and inside the rock mass) described by J. Litwiniszyn (1954) for the loose medium model (Sikora, 2011). Fig. 3. Sample characteristic of distribution function in deterministic model of rock mass in the scope of defined neighborhood for the individual cell in case of simulation of subsidence (Sikora, 2010) The simulation is performed separately for each cell (in the copy of main grid). The following cells, taken to simulation, are selected accordingly to the direction of mining operations. Fig. 4. Symmetrical distribution of subsidence from single cell (Bialek and Sikora, 2012) The results are being summed in a common table. This method reproduces real progress of mining operations. Also it is possible to depend the progressive deformation in the model on already passed value of subsidence, which in turn allow for simulation of nonlinear properties of subsidence summing (Bialek & Sikora, 2012). 2.5. Practical implementation of the Cellular Automaton Simulation of rock mass deformation within presented build concept (algorithm) needs to use computer with a dedicated software. For the purpose of the current rock mass model there was written a new program with usage of Visual Studio Integrated Development Studio (IDE). The algorithm was written in Visual Basic.NET programming language. The realization of the main presumption of CA, that during every iteration all cells in the grid must evolve, was done by simple loop "for". Entire grid of cells will evolve so many times until final state conditions will be met (section 2.2). In case of 2-dimensional model there are no special restrictions about the computer resources. Practical tests shows that modern PC-type computers can perform simulation, even for large grids, in a very short time (comparing e.g. with Knothe's method). 3. Fundamentals of simulation of vertical displacements The basis for any considerations related to the practical use of cellular automata to simulate subsidence of rock mass is the relationship (equation 2) binding cell's dimensions Sk and Wk, the depth of mining operation H, the value of maximal subsidence wmax = ag and the maximum slope Tmax of simulated full subsidence trough (the type of trough where can be observed theoretical maximum subsidence wmax) (Sikora, 2011; Bialek & Sikora, 2012). Tmax AP ag Wk Sk H (2) where: AP -- matching factor due to the value of the parameter P. Based on the equation 2, knowing the value of theoretical maximum subsidence (in the full subsidence trough) and the value of maximum slope Tmax can be determined the relationship between the dimensions of the cell and the number of rows in the grid that allows to build a cellular automaton as a model of the rock mass by means of which can be approximated subsidence profile to preset mining conditions. As the result of simulation there's being obtained subsidence distribution in the whole section of the cell's grid. It allows for the analysis of the subsidence trough profile, of the shape described by the integral of the Gaussian function (Sikora, 2010), at any level of rock mass. 4. The problem of simulating the distribution of horizontal displacement Parallel to the simulation of subsidence distribution, can be simulated horizontal displacements. Registration of results takes place in a separate table (grid of cells), which is a copy of a subsidence grid. The values of horizontal displacements in particular passes will depend, like in case of subsidence, on the value of main transition P, in a manner consistent with the algorithm shown in the Fig. 5. Fig. 5. Characteristic of the distribution function for simulation of horizontal displacements in the deterministic rock mass model in the scope of defined cell's neighborhood (own source) The above characteristics of the propagation function for horizontal displacements simulation shows that for assumed system of cells there are being registered only values of subsidence transferred to cells placed on the right and left in a row above the source cell. Transfer to the left cell is being registered as a negative value. Vertical transfer in case of horizontal displacement simulation is omitted (is not registered). As the effect of presented algorithm are obtained, directly and practically independently to the subsidence simulation, results qualitatively consistent with those observed in a reality. Unfortunately results are multiple smaller then observed in reality or being calculated with usage of different methods, e.g. Knothe's theory. However the maximal values are in the same places that in reference results. Although different combinations of the model's parameters it's impossible to obtain results that are also quantitative compatible (in case of present concept). For presented set of model assumptions can be simulated directly only vertical deformations accordingly with the observed in reality. Similar problem of the modeled rock mass displacements has been described, among others, by F. Dymek (1979), who considered analogous solutions for expanding and non-expanding type of loose medium. Described in the literature case of the non-expanding centre corresponds concept of rock mass model as a finite cellular automaton, in which does not appear the "pressure side" between cells. In mentioned article (Dymek, 1979) was presented the general state of stress and displacement of non-expanding loose center in case of plane strain state in the following form (equation 3): z x2 ; E 1 w ; z xx 2 2 x x2 ; 2 x (3) xz ; u w x z where: z -- is the depth value, x, E -- material constants in case of non-expanding centre (Dymek, 1979), xz -- shear stress. The author (Dymek, 1979) by assuming displacement boundary conditions approximately describes the displacement of the ceiling over extracted seam (equation 4), w x, 0 w0 const for x 0, a ; 0 for x 0, a ; 0 w x, (4) has received solution (equation 5), describing rock mass vertical deformation (subsidence) above extracted seam. w x, z wmax 2 x a (5) It should be noted that the presented equation after certain transformation, corresponds to a similar formula (equation 6), which was received by J. Litwiniszyn (1953) in the course of his research on stochastic center. w x, z wmax x z (6) where: -- is a variable (Litwiniszyn, 1953). Finally, substituting the calculated value of subsidence w (equation 5) to the system of equations (equation 3) can be determined the value of horizontal displacement u. From obtained results the author (Dymek, 1979) conclude that non-expanding centre model can be used to describe the state of stress and displacement of rock mass when the extracted seam lies horizontally on a sufficiently large depth. However, the solution leads to very flat horizontal displacements distributions in comparison to those observed by geodetic observations. Finally, in practice they're determined with usage of the parameter B (equation 7) matching vertical and horizontal displacements (Awierszyn, 1947). u x w x BT x (7) The value of the parameter B for the land surface is variously accepted in practice and literature (Knothe, 1984; Drzla, 1978; Popiolek, 1976; Popiolek & Ostrowski, 1981;Tajdu & Tajdu, 2015). Originally, according to W. Budryk (1953), the value of parameter B was determined basing on the following formula (equation 8): B 0, 4 H tan (8) where: -- the angle of main influences range. Finally, taking the value of the parameter B in accordance with the above equation (equation 8), the theoretical value of the maximum horizontal displacements umax, in case of full subsidence trough (and in the case of a half-plane), can be determined from the following relationship (equation 9): umax 0, 4 wmax mm (9) In principle, one could stop at the presented solution. However, due to the potential studies of the impact of the heterogeneity of the rock mass and nonlinear properties of mining deformations summation on the distribution of horizontal displacements (with usage of CA) there's a need to develop a new distribution function (Tajdu, 2014). In the case of cellular automaton in which the distribution of deformation takes place in the entire space, on the basis of numerous tests and studies of F. Dymek solution (1979), it was found the way to simulate horizontal movements both quantitatively and qualitatively consistent with a reference model. The solution may be additional dependence of the propagation function for the horizontal displacements to the ratio of mapped (assumed) cell dimensions Sk, Wk and the distance of the cell, considered in current iteration of simulation, to the cell representing the source void, for which there's being done simulation of elementary distribution of deformations. The schema of properly modified non-linear distribution function in case of simulation of horizontal displacements in a discrete time for the elementary cell (void) was shown in the figure 6. Fig. 6. Distribution function for simulation of horizontal displacements that make the propagation value conditional on the ratio of cell's projected dimensions and the distance L between the concerned cell wm,n and the elementary cell representing extracted panel pi for which is being performed partial simulation of horizontal displacement (own source) In the figure above there was presented part of the grid and the elementary void pi. In the distance of Lpi-m,n there's placed cell of index m,n. This cell can be called main cell. For that cell is being performed transition function in the range of defined cell's neighborhood (in the discrete time). Some part of the void's volume from cell is being ascribed to two cells on the right and left in the row above. The algorithm of transition, defined in the Fig. 6, is never changing. The algorithm is being repeated for every single cell due the fact that all cells must evolve (section 2.5). For a suitably modified distribution function there were performed a series of simulations of horizontal and vertical displacements in which were created full subsidence troughs. Every time there was determined ratio (equation 10) of the value of maximum horizontal displacements umax to the value of maximal subsidence wmax on the level representing the land surface. BU umax wmax (10) Individual simulations differed by adopted depth of mining operations, the value of main transition P, the value of maximum tilt aT (Bialek & Sikora, 2012) and relating values of adopted cell dimensions Sk and Wk. Sample results were shown in the Table 1. Maximal values of subsidence wmax and horizontal displacements umax at given calculation level were determined as the extreme values from corresponding to the current level's row of automaton. TABLE 1 Exemplary simulation results of subsidence and horizontal displacements for different parameters H [m] aT P Sk [m] Wk [m] BU 2,0 2,0 1,6 4,0 2,0 0,8 0,5 0,5 0,3 0,8 7,756 8,908 19,919 3,754 21,937 0,397 0,999 0,799 2,809 0,398 On the basis of numerous results of simulations there was determined formula describing the ratio of maximal horizontal displacements to the maximal subsidence (in the full subsidence trough) in relation to the cell's dimensions, mining operation depth, the value of parameter aT and the value of main transition P (equation 11). BU W H AU k Sk Wk 0,5 1 P aT (11) where: AU -- adjustment parameter depending on the depth of mining operations. The accuracy of the simulation is determined mostly by assumed cell's dimensions (Sikora, 2011). In turn, the accuracy of designation of the AU parameter mainly affects the accuracy of determining the value of the parameter BU. The average value of parameter AU for performed simulations was equal ­0,196 while the standard deviation 0,046. The parameter AU can be determined more accurately by a polynomial dependent on the depth of simulated mining operations H [m]. In case of adoption of the 3rd degree polynomial (equation 12) maximum absolute error of BU, calculated as the difference between the theoretical and practical value, does not exceed 0,01. AU 0, 000000001H 3 0, 000003H 2 0, 0024 H 0,31 (12) On the basis of directly obtained table with distribution of horizontal displacements can be determined table of horizontal strain , which is calculated in individual cell from the following formula (equation 13): um , n m,n 1 um , n Sk (13) Using equation 2 and 11 can be determined simulation parameters that allows directly generate subsidence trough of assumed slope and ratio of maximum subsidence to maximum horizontal displacement. 5. Example of simulation of subsidence and horizontal displacements In order to show the possibilities of presented model there was done simulation of subsidence and horizontal displacements of rock mass caused by extraction of abstract coal seam lying horizontally on the depth of 300 m and length of 800 m. It was assumed that maximum subsidence in subsidence trough will be wmax = ag = 1000 mm. There were also assumed values of maximum tilt factor aT = 2,0 (eg. equivalent of parameter tan in the theory of S. Knothe (1984)) and the value of parameter BU = 0,4. Applying equations 2 and 11, with usage of the last squares method, there were determined values of parameters Sk, Wk and P to met earlier assumptions. -300 0,5 0,4 0,3 0,2 0,1 0 -0,1 -0,2 -0,3 -0,4 -0,5 -0,6 -0,7 -0,8 -0,9 -1 -1,1 -1,2 -1,3 -1,4 -200 -100 0 100 200 300 400 500 600 700 800 900 1 000 1 100 [mm] [m] seam W [m] U [m] Fig. 7. The chart of subsidence W [m] and horizontal displacements U [m] of subsidence trough formed on the model's surface (own source) Interpolated results of simulation in the form of graphs of subsidence, slopes, horizontal displacements and horizontal strains on the surface of the rock mass model were shown in the following pictures (Figure 7 and 8). [mm/m] 0 -2 -4 -6 -8 -10 -300 -200 -100 0 100 200 300 400 500 600 700 800 900 1 000 1 100 [m] seam T [mm/m] [mm/m] Fig. 8. The chart of slopes T [mm/m] and horizontal strains [mm/m] of subsidence trough formed on the model's surface (own source) As it results from the illustrated example obtained values of rock mass subsidence are consistent with results that can be obtained e.g. by Knothe's method for corresponding values of parameters. Obtained subsidence trough as a result of simulation represents its profile that can be described by the integral of the Gaussian function. Maximum value of slope Tmax obtained by rock mass subsidence simulation, for assumed automaton parameters, is equal to the theoretical value resulting from S. Knothe theory (1984), that can be calculated from the following formula (equation 14): Tmax a g tan H 6,67 (14) The maximum value of horizontal displacements, determinated directly from simulation, umax = 400 mm is also consistent with the theoretical value, wherein horizontal displacements are defined by assuming proportionate to the slope accordingly to the postulate of S.G. Awierszyn (1947) (equation 7,8 and 9). Horizontal strains have been determinated on the basis of the horizontal displacements values. Maksimum value max is also consistent with the theoretical value resulting from the S. Knothe formula calculated from the following equation (equation 15): max 0, 6 ag r (15) 6. Summary In the article there were reminded fundamentals of the simplest rock mass model construction as a deterministic finite cellular automaton. The model presents flat cross-section trough the rock mass by its fragmentation into a regular grid of cells. Cells have the same shape and closely adhere to each other. Every cell has attributed dimension of width and length. Subsequently there were described the basic features of the model in relation to the theory of cellular automata. Characterized been cell's neighborhood and the initial conditions of the automaton. In relation to previous publications there were summarized existing solutions in terms of subsidence simulation including basic characteristics of the distribution function and developed formula binding cell's dimensions, the value of main transition P and the depth of mining operations with the value of maximum subsidence and slope in the resulting full subsidence trough on the model's surface. Based on current knowledge, by adding to the model second grid of cells, being a copy of a subsidence grid, there was presented concept of propagation function which allows to simulate parallel subsidence and horizontal displacements quantitatively and qualitatively consistent with a reference model basing on the geometric-integral theory. On the basis of numerous simulations, for assumed distribution function, there was presented dependence determining the ratio of maximal horizontal displacements to the maximal subsidence in the full subsidence trough in relation to specified cell's dimensions, the value of main transition P, the depth of underground mining operations and value of the parameter of maximal tilt aT. Finally the properties of presented model were verified on chosen example of simulation of subsidence and horizontal displacements propagation. Obtained results were quantitatively and qualitatively consistent with results that can be obtained, for corresponding parameters values, eg. with usage of geometric-integral theory. Presented solution characterize by a simplicity in relation to the complexity of obtained results. It should be noted that the algorithm doesn't need any special mathematical knowledge to apply it. It has been shown that for certain parameters values the results of calculations with usage of cellular automaton are consistent with those observed in reality both in case of horizontal and vertical displacements. These advantages and great potential for further development cause that the model has a chance to be practically applied in the future. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Archives of Mining Sciences de Gruyter

Simulation of Rock Mass Horizontal Displacements with Usage of Cellular Automata Theory.

Archives of Mining Sciences , Volume 61 (4) – Dec 1, 2016

Loading next page...
 
/lp/de-gruyter/simulation-of-rock-mass-horizontal-displacements-with-usage-of-6tvFftWk09
Publisher
de Gruyter
Copyright
Copyright © 2016 by the
ISSN
1689-0469
eISSN
1689-0469
DOI
10.1515/amsc-2016-0051
Publisher site
See Article on Publisher Site

Abstract

coal panel. Obtained results were consistent with results obtained by geometric-integral theory. Keywords: horizontal displacements, rock mass deformation, mining area protection, cellular automata W artykule przedstawiono plaski model górotworu jako deterministyczny, skoczony automat komórkowy. Wykorzystana do opisu rozkladu deformacji wewntrz górotworu i na jego powierzchni teoria automatów pozwala w relatywnie prosty sposób uzyska profil niecki obnieniowej zgodny z profilem obserwowanym pomiarami geodezyjnymi na powierzchni terenu. Przedstawiony w pracy najprostszy model górotworu przedstawia jego plaski przekrój w postaci regularnej siatki komórek, które cile do siebie przylegaj i maj jednakowy ksztalt (Rys. 1). Dla modelu zdefiniowano podstawowe parametry automatu komórkowego takie jak: warunki brzegowe, decydujce o pocztku i kocu symulacji, ssiedztwo komórkowe, okrelajce przestrze w siatce komórek, w obrbie której dochodzi do bezporedniej wymiany informacji zapisanych w poszczególnych komórkach oraz funkcj przejcia, która decyduje o ostatecznej charakterystyce rozkladu symulowanego zjawiska w siatce automatu. W artykule zastosowano deterministyczn funkcj rozkladu. * SILESIAN UNIVERSITY OF TECHNOLOGY, FACULTY OF MINING AND GEOLOGY, INSTITUTE OF MINE SURVEYING AND MINING AREA PROTECTION W wyniku licznych prób modelowych stwierdzono, e stosowana funkcja rozkladu dla symulacji rozkladu obnie (Rys. 3) nie pozwala na symulacj ruchów poziomych jakociowo i ilociowo zgodnych z przemieszczeniami poziomymi obserwowanymi w rzeczywistoci. Jako rozwinicie dotychczasowej koncepcji budowy górotworu jako skoczony automat komórkowy, w pracy opisano funkcj rozkladu (Rys. 5), która pozwala, równolegle do symulacji obnie, symulowa przemieszczenia poziome zgodne z wynikami obserwacji geodezyjnych. Na podstawie wyników licznych symulacji komputerowych opisano podstawow matematyczn zaleno (wzór 11) okrelajc stosunek maksymalnych przemieszcze poziomych do maksymalnych obnie w niecce nadpelnej w odniesieniu do przyjtych parametrów opisanego modelu górotworu, tj.: odwzorowywanych w rzeczywistoci wymiarów komórki, wartoci tzw. przejcia glównego, glbokoci eksploatacji oraz parametru maksymalnego nachylenia (bdcego odpowiednikiem np. parametru tg w teorii Budryka-Knothego). Dla pokazania moliwoci zaproponowanego automatu komórkowego wykonano symulacj rozkladu deformacji wewntrz modelu górotworu i wyznaczono profil linii obnie, nachyle, przemieszcze poziomych i odksztalce poziomych powierzchni modelu dla przykladu abstrakcyjnego, wyeksploatowanego pokladu o okrelonych parametrach górniczo-eksploatacyjnych. W wyniku symulacji otrzymano nadpeln nieck obnieniow o ksztalcie opisywanym przez calk z funkcji Gaussa, w której rozklad przemieszcze poziomych byl zbieny z modelem wzorcowym. Slowa kluczowe: przemieszczenia poziome, deformacje górotworu, ochrona terenów górniczych, automat komórkowy 1. Introduction The Cellular Automata Theory, begun in the 1940s of the last century, found widespread application in many fields of science thanks to its simplicity and exceptional complexity of received results (Wolfram, 2002). The Cellular Automaton (CA) was originally discovered by Stanislav Ulam (Pickover and Clifford, 2009) and John von Neumann (Schiff & Joel, 2011) at the Los Alamos National Laboratory. Cellular automata are commonly used e.g. in computer science, to test the algorithms, in physics to simulate the spread of fire, gas, heat, etc., or to simulate the traffic (Chopard at al., 1996). It turns out that despite the huge potential in the field of rock mechanics and mining area protection, so far, the theory wasn't practically applied. Today, in Poland, predictions of the impact of underground mining on the surface mainly rely on methods which uses geometric-integral theory to describe the distribution of deformation and methods basing on models of continuous centre where the state of stresses and displacements defines a system of differential equations and the equation of state dependent on the adopted model. In previous publications there was proved (Bialek & Sikora, 2012; Sikora, 2010, 2011, 2013, 2014) that using cellular automata to describe the rock mass deformation, results can be very similar, and in some cases even identical with the results obtained with application of different methods, e.g. Knothe's method. So far, it has been shown that for certain conditions the theory of cellular automata allows in a relatively easy way to get subsidence trough profile consistent with the profile observed by geodetic measurements. Moreover it is possible to simulate directly the impact of heterogeneous build of rock mass and nonlinear properties of summing of the mining influences on the subsidence distribution process. Especially there was presented the method of simulation of the impact of faults on the distribution of subsidence. Normal faults, in a form of splits, were implemented on the edges of corresponding cells. Modified distribution function allows for limitation of subsidence propagation trough the split. Nonlinear properties of subsidence summation were implemented conditioning the value of subsidence propagation to the sustained tilt values. In a consequence it is possible to simulate e.g. the impact of multiple extraction on the subsidence distribution. Finally there was presented build concept of the spatial rock mass model as a cellular automaton. In this paper, which is a development of the prior solutions, there will be presented a method of horizontal displacements simulation. In consequence the solution will allow for the characteristics of the horizontal displacements and horizontal deformations within the rock mass and on its surface. The assumptions and the method of implementation were shown on the example of the simplest two-dimensional rock mass model (Sikora, 2011) built in relation to the deterministic finite automaton definition (Moore, 1956). 2. The model of the rock mass as a deterministic finite automaton By cellular automata as a rock mass model can be specified a number of commonly used algorithms, e.g. rock mass model described by T. Niemiec (Niemiec, 1985; Niemiec et al., 2008). However in practise this term is being used in some specific cases and applications. The simplest definition describes cellular automaton as a mathematical model that its structures are determined by the net of cells, theirs states, transitions and rules of transitions (Moore, 1956). In relation to the previous definition there can be build the simplest rock mass model as its flat cross-section presented by the net of discrete cells. Detailed built concept of that model and its variants were described in previous publications (Sikora, 2010, 2011). However there is a group of fundamental parameters describing a model of the rock mass as a deterministic finite automaton, including: the type of cell, the net of cells, so called cell's neighbouring, boundary conditions and a transition function that decide in most about the final characteristics of the simulated phenomena. Subsequently, these parameters are shortly discussed. 2.1. Cell type and the grid of cells Described case of the rock mass cross-section presume its segmentation into regular net of cells of determined number of rows m and columns n. All cells has the same shape of rectangle and the same size and also closely adhere to each other (Fig. 1). Cell's dimensions are corresponding with real dimensions of width Sk [m] and height Wk [m]. Thus formed table (grid of cells) represents rock mass in a flat system, where vertical axis is converging with columns (depth). Individual rows relate to the respective levels of rock mass, however the last row m (the highest) relates to the land surface. 2.2. Initial state In case of presented type of cellular automaton it is assumed that every cell starts in the same state, except for a finite number of cells in other states. It is called configuration (Schiff and Joel, 2011) or initial state. Fig. 1. Cellular automata grid representing horizontal section trough rock mass and extracted longwall (Sikora, 2013) In case of considered rock mass model, as a initial state, can be assumed two-dimensional grid of cells. For every single cell can be assigned a value representing the void. At the beginning values assigned to proper cells (initial subsidence) shows performed mining operations (equation 1). pi, j xm,n = agp where: pi, j xm,n wm,n gp a -- -- -- -- -- (1) elementary particle of the seam with coordinates i, j, a single cell in the grid of automaton with coordinates m,n, initial subsidence in the current cell with coordinates m,n, the average thickness of the seam, subsidence factor ­ the theory's parameter that shows the difference between the volume of extracted seam and the subsidence trough on the surface (subsidence at the surface depends on the strength and nature of the rocks overlying the seam). The above equation means that for every particular element of the seam pi, jj (i, j are its real coordinates) corresponds at least one cell xm,n in the CA's grid. For that cell it's being ascribed the value of extraction height (of the seam) gp multiplied by the value of subsidence factor a. The final state is the moment when the sum of void's volume in the row representing the surface of the mining area equals the sum of void's volume assigned to specified cells at the initial state. It should be noted that the initial conditions does not affect the cell shape. Theoretically cells are considered to be shapeless and it's not required to visualize them. There will be analysed only the effect of their work. 2.3. Cell's neighborhood As a neighborhood of the individual cell meant all cells that are directly related by a dependence function (Schiff and Joel, 2011). NW W SW NE E SE Fig. 2. Moore type neighborhood defined for the dark cell in the regular net of Cellular Automaton (own source) In the rock mass as a closest neighbourhood of a post-operational void can be understood the immediate surrounding area (cave zone). In the CA theory every single cell has the same neighbourhood. One of the most popular type of neighbourhood is the Moore's type (Fig. 2) which include all adjacent cells to the current cell. Formally this is the best type of neighbourhood that could be considered in case of mapped void in the rock mass model. The force of gravity causes the formed empty space in the rock mass (void) is being filled with overhanging rocks. Taking this assumption and the lack of existence of any additional forces affecting tightening of the void, there can be finally determined the cell's neighborhood. Therefore it's limited to three adjacent cells in the row above the main cell (cells NW, N and NE in the Fig. 2). In other words every cell will interact only with neighbor cells. 2.4. Propagation function in order to the rock mass subsidence simulation The manner of operation, compatibility features of the model with the real object and the state of individual elementary cell in a discrete time decide first and foremost the evolution rule of cellular automaton. The evolution of the model base on the boundary conditions defined for the automaton and the propagation function that determine intercellular relations in the space defined by the cell's neighborhood. In case of presented deterministic model (there is also stochastic model (Sikora, 2011)) fulfilment of the post-operational void rely on determination of partial distribution of the void's volume outflow in the range of cell's neighborhood. Presented algorithm in the Fig. 3 assumes that the partial value of a void is filled by the rocks volume from a cell placed directly above the void. It's so called main transition P. Other transitions take place symmetrically with respect to the main transition to cells from the neighborhood (on the left and right). For values of P = 0,5 presented characteristics of the propagation function leads to the subsidence distribution (on the surface and inside the rock mass) described by J. Litwiniszyn (1954) for the loose medium model (Sikora, 2011). Fig. 3. Sample characteristic of distribution function in deterministic model of rock mass in the scope of defined neighborhood for the individual cell in case of simulation of subsidence (Sikora, 2010) The simulation is performed separately for each cell (in the copy of main grid). The following cells, taken to simulation, are selected accordingly to the direction of mining operations. Fig. 4. Symmetrical distribution of subsidence from single cell (Bialek and Sikora, 2012) The results are being summed in a common table. This method reproduces real progress of mining operations. Also it is possible to depend the progressive deformation in the model on already passed value of subsidence, which in turn allow for simulation of nonlinear properties of subsidence summing (Bialek & Sikora, 2012). 2.5. Practical implementation of the Cellular Automaton Simulation of rock mass deformation within presented build concept (algorithm) needs to use computer with a dedicated software. For the purpose of the current rock mass model there was written a new program with usage of Visual Studio Integrated Development Studio (IDE). The algorithm was written in Visual Basic.NET programming language. The realization of the main presumption of CA, that during every iteration all cells in the grid must evolve, was done by simple loop "for". Entire grid of cells will evolve so many times until final state conditions will be met (section 2.2). In case of 2-dimensional model there are no special restrictions about the computer resources. Practical tests shows that modern PC-type computers can perform simulation, even for large grids, in a very short time (comparing e.g. with Knothe's method). 3. Fundamentals of simulation of vertical displacements The basis for any considerations related to the practical use of cellular automata to simulate subsidence of rock mass is the relationship (equation 2) binding cell's dimensions Sk and Wk, the depth of mining operation H, the value of maximal subsidence wmax = ag and the maximum slope Tmax of simulated full subsidence trough (the type of trough where can be observed theoretical maximum subsidence wmax) (Sikora, 2011; Bialek & Sikora, 2012). Tmax AP ag Wk Sk H (2) where: AP -- matching factor due to the value of the parameter P. Based on the equation 2, knowing the value of theoretical maximum subsidence (in the full subsidence trough) and the value of maximum slope Tmax can be determined the relationship between the dimensions of the cell and the number of rows in the grid that allows to build a cellular automaton as a model of the rock mass by means of which can be approximated subsidence profile to preset mining conditions. As the result of simulation there's being obtained subsidence distribution in the whole section of the cell's grid. It allows for the analysis of the subsidence trough profile, of the shape described by the integral of the Gaussian function (Sikora, 2010), at any level of rock mass. 4. The problem of simulating the distribution of horizontal displacement Parallel to the simulation of subsidence distribution, can be simulated horizontal displacements. Registration of results takes place in a separate table (grid of cells), which is a copy of a subsidence grid. The values of horizontal displacements in particular passes will depend, like in case of subsidence, on the value of main transition P, in a manner consistent with the algorithm shown in the Fig. 5. Fig. 5. Characteristic of the distribution function for simulation of horizontal displacements in the deterministic rock mass model in the scope of defined cell's neighborhood (own source) The above characteristics of the propagation function for horizontal displacements simulation shows that for assumed system of cells there are being registered only values of subsidence transferred to cells placed on the right and left in a row above the source cell. Transfer to the left cell is being registered as a negative value. Vertical transfer in case of horizontal displacement simulation is omitted (is not registered). As the effect of presented algorithm are obtained, directly and practically independently to the subsidence simulation, results qualitatively consistent with those observed in a reality. Unfortunately results are multiple smaller then observed in reality or being calculated with usage of different methods, e.g. Knothe's theory. However the maximal values are in the same places that in reference results. Although different combinations of the model's parameters it's impossible to obtain results that are also quantitative compatible (in case of present concept). For presented set of model assumptions can be simulated directly only vertical deformations accordingly with the observed in reality. Similar problem of the modeled rock mass displacements has been described, among others, by F. Dymek (1979), who considered analogous solutions for expanding and non-expanding type of loose medium. Described in the literature case of the non-expanding centre corresponds concept of rock mass model as a finite cellular automaton, in which does not appear the "pressure side" between cells. In mentioned article (Dymek, 1979) was presented the general state of stress and displacement of non-expanding loose center in case of plane strain state in the following form (equation 3): z x2 ; E 1 w ; z xx 2 2 x x2 ; 2 x (3) xz ; u w x z where: z -- is the depth value, x, E -- material constants in case of non-expanding centre (Dymek, 1979), xz -- shear stress. The author (Dymek, 1979) by assuming displacement boundary conditions approximately describes the displacement of the ceiling over extracted seam (equation 4), w x, 0 w0 const for x 0, a ; 0 for x 0, a ; 0 w x, (4) has received solution (equation 5), describing rock mass vertical deformation (subsidence) above extracted seam. w x, z wmax 2 x a (5) It should be noted that the presented equation after certain transformation, corresponds to a similar formula (equation 6), which was received by J. Litwiniszyn (1953) in the course of his research on stochastic center. w x, z wmax x z (6) where: -- is a variable (Litwiniszyn, 1953). Finally, substituting the calculated value of subsidence w (equation 5) to the system of equations (equation 3) can be determined the value of horizontal displacement u. From obtained results the author (Dymek, 1979) conclude that non-expanding centre model can be used to describe the state of stress and displacement of rock mass when the extracted seam lies horizontally on a sufficiently large depth. However, the solution leads to very flat horizontal displacements distributions in comparison to those observed by geodetic observations. Finally, in practice they're determined with usage of the parameter B (equation 7) matching vertical and horizontal displacements (Awierszyn, 1947). u x w x BT x (7) The value of the parameter B for the land surface is variously accepted in practice and literature (Knothe, 1984; Drzla, 1978; Popiolek, 1976; Popiolek & Ostrowski, 1981;Tajdu & Tajdu, 2015). Originally, according to W. Budryk (1953), the value of parameter B was determined basing on the following formula (equation 8): B 0, 4 H tan (8) where: -- the angle of main influences range. Finally, taking the value of the parameter B in accordance with the above equation (equation 8), the theoretical value of the maximum horizontal displacements umax, in case of full subsidence trough (and in the case of a half-plane), can be determined from the following relationship (equation 9): umax 0, 4 wmax mm (9) In principle, one could stop at the presented solution. However, due to the potential studies of the impact of the heterogeneity of the rock mass and nonlinear properties of mining deformations summation on the distribution of horizontal displacements (with usage of CA) there's a need to develop a new distribution function (Tajdu, 2014). In the case of cellular automaton in which the distribution of deformation takes place in the entire space, on the basis of numerous tests and studies of F. Dymek solution (1979), it was found the way to simulate horizontal movements both quantitatively and qualitatively consistent with a reference model. The solution may be additional dependence of the propagation function for the horizontal displacements to the ratio of mapped (assumed) cell dimensions Sk, Wk and the distance of the cell, considered in current iteration of simulation, to the cell representing the source void, for which there's being done simulation of elementary distribution of deformations. The schema of properly modified non-linear distribution function in case of simulation of horizontal displacements in a discrete time for the elementary cell (void) was shown in the figure 6. Fig. 6. Distribution function for simulation of horizontal displacements that make the propagation value conditional on the ratio of cell's projected dimensions and the distance L between the concerned cell wm,n and the elementary cell representing extracted panel pi for which is being performed partial simulation of horizontal displacement (own source) In the figure above there was presented part of the grid and the elementary void pi. In the distance of Lpi-m,n there's placed cell of index m,n. This cell can be called main cell. For that cell is being performed transition function in the range of defined cell's neighborhood (in the discrete time). Some part of the void's volume from cell is being ascribed to two cells on the right and left in the row above. The algorithm of transition, defined in the Fig. 6, is never changing. The algorithm is being repeated for every single cell due the fact that all cells must evolve (section 2.5). For a suitably modified distribution function there were performed a series of simulations of horizontal and vertical displacements in which were created full subsidence troughs. Every time there was determined ratio (equation 10) of the value of maximum horizontal displacements umax to the value of maximal subsidence wmax on the level representing the land surface. BU umax wmax (10) Individual simulations differed by adopted depth of mining operations, the value of main transition P, the value of maximum tilt aT (Bialek & Sikora, 2012) and relating values of adopted cell dimensions Sk and Wk. Sample results were shown in the Table 1. Maximal values of subsidence wmax and horizontal displacements umax at given calculation level were determined as the extreme values from corresponding to the current level's row of automaton. TABLE 1 Exemplary simulation results of subsidence and horizontal displacements for different parameters H [m] aT P Sk [m] Wk [m] BU 2,0 2,0 1,6 4,0 2,0 0,8 0,5 0,5 0,3 0,8 7,756 8,908 19,919 3,754 21,937 0,397 0,999 0,799 2,809 0,398 On the basis of numerous results of simulations there was determined formula describing the ratio of maximal horizontal displacements to the maximal subsidence (in the full subsidence trough) in relation to the cell's dimensions, mining operation depth, the value of parameter aT and the value of main transition P (equation 11). BU W H AU k Sk Wk 0,5 1 P aT (11) where: AU -- adjustment parameter depending on the depth of mining operations. The accuracy of the simulation is determined mostly by assumed cell's dimensions (Sikora, 2011). In turn, the accuracy of designation of the AU parameter mainly affects the accuracy of determining the value of the parameter BU. The average value of parameter AU for performed simulations was equal ­0,196 while the standard deviation 0,046. The parameter AU can be determined more accurately by a polynomial dependent on the depth of simulated mining operations H [m]. In case of adoption of the 3rd degree polynomial (equation 12) maximum absolute error of BU, calculated as the difference between the theoretical and practical value, does not exceed 0,01. AU 0, 000000001H 3 0, 000003H 2 0, 0024 H 0,31 (12) On the basis of directly obtained table with distribution of horizontal displacements can be determined table of horizontal strain , which is calculated in individual cell from the following formula (equation 13): um , n m,n 1 um , n Sk (13) Using equation 2 and 11 can be determined simulation parameters that allows directly generate subsidence trough of assumed slope and ratio of maximum subsidence to maximum horizontal displacement. 5. Example of simulation of subsidence and horizontal displacements In order to show the possibilities of presented model there was done simulation of subsidence and horizontal displacements of rock mass caused by extraction of abstract coal seam lying horizontally on the depth of 300 m and length of 800 m. It was assumed that maximum subsidence in subsidence trough will be wmax = ag = 1000 mm. There were also assumed values of maximum tilt factor aT = 2,0 (eg. equivalent of parameter tan in the theory of S. Knothe (1984)) and the value of parameter BU = 0,4. Applying equations 2 and 11, with usage of the last squares method, there were determined values of parameters Sk, Wk and P to met earlier assumptions. -300 0,5 0,4 0,3 0,2 0,1 0 -0,1 -0,2 -0,3 -0,4 -0,5 -0,6 -0,7 -0,8 -0,9 -1 -1,1 -1,2 -1,3 -1,4 -200 -100 0 100 200 300 400 500 600 700 800 900 1 000 1 100 [mm] [m] seam W [m] U [m] Fig. 7. The chart of subsidence W [m] and horizontal displacements U [m] of subsidence trough formed on the model's surface (own source) Interpolated results of simulation in the form of graphs of subsidence, slopes, horizontal displacements and horizontal strains on the surface of the rock mass model were shown in the following pictures (Figure 7 and 8). [mm/m] 0 -2 -4 -6 -8 -10 -300 -200 -100 0 100 200 300 400 500 600 700 800 900 1 000 1 100 [m] seam T [mm/m] [mm/m] Fig. 8. The chart of slopes T [mm/m] and horizontal strains [mm/m] of subsidence trough formed on the model's surface (own source) As it results from the illustrated example obtained values of rock mass subsidence are consistent with results that can be obtained e.g. by Knothe's method for corresponding values of parameters. Obtained subsidence trough as a result of simulation represents its profile that can be described by the integral of the Gaussian function. Maximum value of slope Tmax obtained by rock mass subsidence simulation, for assumed automaton parameters, is equal to the theoretical value resulting from S. Knothe theory (1984), that can be calculated from the following formula (equation 14): Tmax a g tan H 6,67 (14) The maximum value of horizontal displacements, determinated directly from simulation, umax = 400 mm is also consistent with the theoretical value, wherein horizontal displacements are defined by assuming proportionate to the slope accordingly to the postulate of S.G. Awierszyn (1947) (equation 7,8 and 9). Horizontal strains have been determinated on the basis of the horizontal displacements values. Maksimum value max is also consistent with the theoretical value resulting from the S. Knothe formula calculated from the following equation (equation 15): max 0, 6 ag r (15) 6. Summary In the article there were reminded fundamentals of the simplest rock mass model construction as a deterministic finite cellular automaton. The model presents flat cross-section trough the rock mass by its fragmentation into a regular grid of cells. Cells have the same shape and closely adhere to each other. Every cell has attributed dimension of width and length. Subsequently there were described the basic features of the model in relation to the theory of cellular automata. Characterized been cell's neighborhood and the initial conditions of the automaton. In relation to previous publications there were summarized existing solutions in terms of subsidence simulation including basic characteristics of the distribution function and developed formula binding cell's dimensions, the value of main transition P and the depth of mining operations with the value of maximum subsidence and slope in the resulting full subsidence trough on the model's surface. Based on current knowledge, by adding to the model second grid of cells, being a copy of a subsidence grid, there was presented concept of propagation function which allows to simulate parallel subsidence and horizontal displacements quantitatively and qualitatively consistent with a reference model basing on the geometric-integral theory. On the basis of numerous simulations, for assumed distribution function, there was presented dependence determining the ratio of maximal horizontal displacements to the maximal subsidence in the full subsidence trough in relation to specified cell's dimensions, the value of main transition P, the depth of underground mining operations and value of the parameter of maximal tilt aT. Finally the properties of presented model were verified on chosen example of simulation of subsidence and horizontal displacements propagation. Obtained results were quantitatively and qualitatively consistent with results that can be obtained, for corresponding parameters values, eg. with usage of geometric-integral theory. Presented solution characterize by a simplicity in relation to the complexity of obtained results. It should be noted that the algorithm doesn't need any special mathematical knowledge to apply it. It has been shown that for certain parameters values the results of calculations with usage of cellular automaton are consistent with those observed in reality both in case of horizontal and vertical displacements. These advantages and great potential for further development cause that the model has a chance to be practically applied in the future.

Journal

Archives of Mining Sciencesde Gruyter

Published: Dec 1, 2016

References