Advances in Civil Engineering
Volume 2021 (2021), Article ID 9855067, 21 pages
https://doi.org/10.1155/2021/9855067

Research Article

Numerical Simulation of Mining-Induced Damage in Adjacent Tunnels Based on FLAC3D

Guorui Wang,1,2 Peixian Li,¹ Qiang Wu,1,3 Ximin Cui,¹ and Zhixiang Tan⁴

¹College of Geoscience and Surveying Engineering, China University of Mining and Technology (Beijing), Beijing 100083, China
²Institute of Land and Resources Investigation and Monitoring of Ningxia Hui Autonomous Region, Yinchuan 750000, China
³National Coal Mine Water Hazard Prevention Engineering Technology Research Center, Beijing 10003, China
⁴School of Environment Science and Spatial Information, China University of Mining and Technology, Xuzhou 221116, China

Correspondence should be addressed to Peixian Li; pxlicumt@126.com

Received 16 April 2021; Accepted 4 June 2021; Published 25 June 2021

Academic Editor: Loke Foong

Copyright © 2021 Guorui Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The Bayueshan (BYS) tunnels are adjacent tunnels which comprise key engineering structures important for Yurong highway from Chongqing to Sichuan, China. The tunnels were established crossing over the old goafs of a coal mine. The tunnels have been damaged three times between 2014 and 2016. To evaluate the reasons for the damage, fast Lagrangian analysis of continua in three dimensions () was used to simulate the underground coal-mine excavation yearly from 2012 to 2015. The simulation results show that the subsidence of the ground surface increased annually. The maximum subsidence of the left tunnel was 505 mm and that of the right tunnel was 512 mm. The excavation led to increase in the tunnels’ stress concentration. The simulations further show that uneven horizontal displacement and subsidence were the major factors in inducing tunnel construction defects. Poor geological conditions, water, and vehicle dynamic loads were also important factors that induced damage. Several engineering suggestions for ways to maintain the tunnels are given in this paper, and the research provides a reference for safety evaluation for tunnels crossing over an old goaf.

1. Introduction

Underground mining can change rock’s original stress condition, and the rock mass will move to achieve a new balance [1, 2]. Mining-induced subsidence will cause deformation of the construction above a goaf or destruction [3–5]. In China, with rapid econiomic development in recent decades, numerous highways and railways have been established in which tunnels must pass through a coal-mining area or over old goafs [6–8]. There are mainly three major methods used to predict coal-mine subsidence of tunnels. The first is based on the model of stochastic medium theory [9, 10], the second on the theory of the elastic-plastic medium theory [11, 12], and the third and most used is numerical simulation [13–15]. The stochastic medium method is a semiempirical influence function, and it can be used to obtain reliable subsidence results with the parameters of field observation [1], but it lacks a reasonable explanation of deformation mechanisms. The elastic-plastic medium theory is often based on a hypothesis and simplified model, but it cannot reflect the actual geological and mining conditions, and its results are often significantly different from actual ones [16, 17]. Numerical simulation is currently the method that is used mostly in tunnel damage analysis [18–20]. Anasys, , UDEC, and 3DEC are mostly used numerical simulation software [13, 21, 22]. Chongqing Bayueshan (BYS) tunnel is an important length of highway engineering from Chongqing to Sichuan, China. The BYS tunnel was established over a coal mine, and some parts extend through the mine’s old goaf. The coal mine is still in production, and the underground mining and old goaf have caused the tunnels to crack, beginning in 2014. To evaluate the deformation influence on the tunnels and provide more information on tunnel maintenance, a fast Lagrangian analysis of continua in three-dimension () model was devised to simulate the tunnel deformation and damage. The research results are used as a basis for BYS tunnel maintenance design engineers.

2. Geological Conditions

According to the bore hole, the research area strata are mainly constituted by the Quaternary Xujiahe Formation of the Triassic upper system and the Triassic Jialingjiang Formation. Xujiahe Formation is the main coal-bearing stratum, as shown in Figure 1. Lithologies of the rock are gray black shale, sandy shale, brown yellow mud stone, and gray post. There are six lithological segments, the thicknesses of which are in the range of approximately 400600 m. This formation is parallel nonconformity with underlying strata. The investigated area contains three segments: the second, third, and fourth segments.

Figure 1: Geological profile of left tunnel (reproduced from [10]).

The second segment is a light gray, thin to medium thick layer, constituted by coarse feldspar quartz sandstone. Its thickness is 160 m200 m.

The third segment is gray black, thin to medium thick layer, a mud stone and silt stone interbed, partly with carbon mud stone and thin coal seam. The coal seam is 0.15 m0.35 m. It is the main coal seam in this area. The layer thickness is approximately 40 m70 m. This segment belongs to soft rock. The rock has poor interlayer bonding and good joint development, and the rock mass has integrity.

The fourth segment is a light gray, thin to medium thick layer with coarse feldspar quartz sandstone and partly boulder clay. This segment thickness is 5070 m. This segment belongs to hard rock with simple interlayer bonding and good joint development, and the rock mass has integrity.

3. Simulation Model

3.1. Introduction of and Mohr–Coulomb Criterion

is a finite-difference program developed by Itasca Co. The software can simulate different materials with mechanical behavior of failure and plastic flow. It is good at simulating problems of failure and instability, especially numerical stability with large displacements, including collapse [23]. uses the explicit solution method to solve the nonlinear problem quickly with small memory overhead. It provides 12 material constitutive models: a null model, three elastic models, and eight plasticity models. The Mohr–Coulomb strength criterion is a nonlinear constitutive model often used in rock and soil engineering numerical simulation [15, 19] and is expressed by the following equation:

where and are the maximum and minimum principal stress, respectively, and and are the cohesion and internal friction angle, respectively. When , the material will fail by shearing. Under a normal stress state, the rock can fail by tensile strength. In the work described in this paper, the Mohr–Coulomb model was used for the rock and the null model for the excavation area.

3.2. Boundary Condition and Model Geometry

According to the geological conditions of two coal mines, the model was built with 1 : 1 geometric scale. The rock was simplified into layered rock mass with single thickness. The simulation models are free displacement only with a boundary displacement constraint. As shown in Figure 2, the Z displacement of the model bottom, the X displacement of the left and right sides, and the Y displacement of both the front and back sides are fixed. To obtain the subsidence, the top of the model is set to free. The load applied to the studied model is only under gravity. No tectonic stress is considered in this research.

(a)

(b)

(c)

Figure 2: Boundary conditions of the model.

Tunnel crack is a plane strain problem, and, to reduce the boundary effect, the geometry is decided by the deformation edge, which should be large enough for investigating area mining subsidence. The model geometry used in this paper is shown in Figure 3(a). Coal and rock layers in the model were generated with a wedge mesh. The twin-tunnels model is shown in Figure 3(b). The tunnels were generated with a radial cylinder and cylindrical shell. The wedge, radial cylinder, and cylindrical shell were all built in mesh elements in the software. The built model has 106,586 zones and 414,359 grid points. The excavation area used the null model.

(a)

(b)

(c)

Figure 3: Geometry of the model: (a) three-dimensional simulation model; (b) tunnel geometry; (c) relative position of coal seams and tunnel.

3.3. Rock Properties

The actual rock properties were not available to the authors, so several average values [24–26] were used in the model, as shown in Table 1. The parameters in the table are laboratory test-piece values, but naturally rock is a large rock mass and has numerous joints and cracks. The results of a significant amount of research show that strength of rock mass is approximately 1/101/30 of that of the test pieces [9, 27, 28]. In the present paper, the bulk and shear moduli applied to the model were both 1/15 of the laboratory values.

Table 1: Laboratory rock-mass properties.

Rock

()
()
()
()
()
()

Coal
1,420
6.1
1.9
2.11
2.6
29.5

Mudstone
2,437
15.6
8.9
2.2
1.8
38

Sandstone
2,690
10.7
6.15
1.55
1.38
41

Sprayed concrete
2,300
14.0
10.2
2.0
1.5
50

: density; : bulk modulus; : shear modulus; : cohesion; : tensile strength; : internal friction angle.

4. Results

4.1. Subsidence Characteristics of Mining

To obtain the damage mechanism, four simulation steps were designed according to the yearly excavation engineering conducted after the tunnel was completed. The subsidence of the surface is shown in Figure 4, illustrating that upon excavation, both the surface-subsidence value and size of the influence area are increased yearly. The maximum surface-subsidence values were 290 mm (2012), 335 mm (2013), 457 mm (2014), and 471 mm (2015).

(a)

(b)

(c)

(d)

Figure 4: Subsidence contour of coal mine (units: mm): (a) 2012; (b) 2013; (c) 2014; (d) 2015.

Figure 5 shows the subsidence contour of the tunnel profile. With tunnel underground excavations of , , and , the tunnels’ subsidence values and influence area sizes all increased yearly. The subsidence of the sectional elevation contours is shown in Figure 5. The slice contour shows the change process of the rock after coal-mine excavation. As shown in the figure, the maximum subsidences are directly over the top of the working panel and decrease from the coal roof to the surface. The subsidences of the tunnel increased yearly, and there were not uniformly distributed. The uneven settlement often causes stress concentration and structural failure.

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

Figure 5: Subsidence of tunnel slice: (a) Left 2012; (b) Left 2013; (c) Left 2014; (d) Left 2015; (e) Right 2012; (f) Right 2013; (g) Right 2014; (h) Right 2015.

4.2. Stress Fields of Tunnels

The maximum principal stresses of the tunnels are shown in Figure 6. The stress concentrations are located at the roof of the both tunnels, and the influence areas are LK35LK35 + 160 and RK35 + 13RK35 + 156. The maximums of the maximum principal stresses are shown in Table 2. The XX stresses (stresses along the tunnel) are shown in Figure 7. The XX-stress concentrations are located at LK34 + 994LK35 + 153 for the left tunnel and at RK34 + 987RK35 + 123 for the right tunnel. The YY stresses (stresses perpendicular to tunnel) are shown in Figure 8. The YY-stress concentrations are located at LK35 + 53  LK35 + 154 for the left tunnel and at RK35 + 31RK35 + 134 for the right tunnel. The ZZ stresses (stresses along the direction of plumb) are shown in Figure 9. The ZZ-stress concentrations are located at LK35 + 212LK35 + 333 for the left tunnel and at RK35 + 189RK35 + 370 for the right tunnel. As shown in the stress-contour figures, both the roof and floor of the tunnel were bearing stress concentration. Both tunnels have a stress-concentrated area and exhibit structural stress imbalance.

(a)

(b)

(c)

(d)

Figure 6: Maximum principal stress contours of tunnels (unit: MPa): (a) 2012; (b) 2013; (c) 2014; (d) 2015.

Table 2: Maximums deformation values of tunnels.

Item
Tunnel
2012
2013
2014
2015

Value
Position
Value
Position
Value
Position
Value
Position

(MPa)
Left
14.02
124
13.87
125
13.13
113
13.6
114

Right
16.96
101
17.07
100
19.18
100
19.61
101

(MPa)
Left
9.05
133
9.02
133
10.21
133
10.35
133

Right
11.45
100
12.07
101
13.09
100
13.22
100

(MPa)
Left
6.57
130
6.72
130
6.72
131
6.84
131

Right
11.2
100
11.74
101
12.74
100
12.91
100

(MPa)
Left
6.25
124
8.20
–7
8.20
–9
8.43
–7

Right
6.96
100
6.38
–24
8.07
–20
7.00
–21

(mm)
Left
83
294
102
330
149
335
161
335

Right
84
252
102
365
156
368
152
368

(mm)
Left
–55
123
–67
122
–96
123
–99
122

Right
–48
182
–56
182
–75
180
–75
174

(mm)
Left
–328
193
–372
193
–493
193
–505
193

Right
–332
180
–375
181
–498
181
–512
181

(a)

(b)

(c)

(d)

Figure 7: SXX contours of tunnels (unit: MPa): (a) 2012; (b) 2013; (c) 2014; (d) 2015.

(a)

(b)

(c)

(d)

Figure 8: SYY contours of tunnels (unit: MPa): (a) 2012; (b) 2013; (c) 2014; (d) 2015.

(a)

(b)

(c)

(d)

Figure 9: SZZ contour of tunnels (unit: MPa): (a) 2012; (b) 2013; (c) 2014; (d) 2015.

4.3. Displacement Characteristics of Tunnels

The X displacements (displacements along the tunnel direction) are shown in Figure 10. The X displacement increased yearly, and the mostly uneven X-direction horizontal movements were located at LK35 + 84LK35 + 107 and RK35 + 71RK35 + 127. The Y displacements (displacements perpendicular to the tunnel direction) are shown in Figure 11. The much larger uneven influence areas of the Y displacements were located at LK35 + 39LK35 + 154 and RK35 + 29RK35 + 137. The Z displacements (subsidences) are shown in Figure 12. The maximum subsidence of the left tunnel is 505 mm at LK35 + 193 and that of the right tunnel is 512 mm at RK35 + 181. The subsidence increased yearly, and the areas of subsidence that exceeded 450 mm were LK35 + 157LK35 + 250 and RK35 + 154RK35 + 259.

(a)

(b)

(c)

(d)

Figure 10: X-displacement contours of tunnels (unit: mm): (a) 2012; (b) 2013; (c) 2014; (d) 2015.

(a)

(b)

(c)

(d)

Figure 11: Y-displacement contours of tunnels (unit: mm): (a) 2012; (b) 2013; (c) 2014; (d) 2015.

(a)

(b)

(c)

(d)

Figure 12: Z-displacement contours of tunnels (unit: mm): (a) 2012; (b) 2013; (c) 2014; (d) 2015.

4.4. Statistics of Maximum Values and Locations

The maximum deformation values of principal stress, XX stress, YY stress, ZZ stress, and their displacements are shown in Table 2.

5. Tunnel Damage Analysis

The first crack in the tunnel occurred on September 4, 2014, at LK35 + 162LK35 + 171 after a heavy rain. The concrete surface was cracked (the crack was 1 cm wide and 4 m high) from the road surface to the arch, as shown in Figure 13.

Figure 13: Crack of LK35 + 168 (September 4, 2014) (source from [10]).

On November 15, 2014, the right tunnel was cracked at RK35 + 268 (the crack was 3 cm wide and 6 m high), as shown in Figure 14(a). On July 7, 2016, the tunnel was cracked after a heavy rain at LK35 + 220; the lining of the wall was cracked, and the fracture penetration is shown in Figure 14(b).

(a)

(b)

Figure 14: Images of tunnel cracks (source from [10]): (a) RK35 + 268; (b) LK35 + 220.

The numerical simulation results show that the maximum subsidence was laid up on the downhill side. The numerical simulation surface deformation law and subsidence value, range, and size basically conform to the conventional mining subsidence law. The subsidence, stress, and strain can be used as a reliable simulation results in the tunnel damage analysis.

Construction of the tunnel in which the cracked areas are located was completed on April 26, 2011. After being in service for approximately three years, the tunnel deformation and cracking mainly resulted from supplementary loading. Several noteworthy characteristics are as follows.

First, the tunnel damage mainly occurred after a heavy rain or after the rainy season. The research area was mainly formed by a fractured rigid rock mass. In fact, some fissures are connected with the goaf or old goaf. Underground mining is the major factor for the tunnel cracking, and the rain is an important cause of it.

Second, the cracked area did not locate above the goaf directly but located outside the goaf boundary.

Third, the tunnel support structure is also influenced the stability of the tunnel. Some cracks occurred near the support change area.

Fourth, all of the cracks occurred at the tunnel movement joint position. The joints moved together, or there was relative movement between joints which caused the lining to crack.

6. Discussion

The BYS tunnel cracks have damaged some structure of the tunnel, and the dropping cement has led a mining car to drop out. As the above simulation results and analysis show, the reason for the tunnels’ defects was mainly caused by the underground excavation of the ChangRong coal mine. It led to the tunnel subsidence and uneven horizontal movement. The rock movement caused the tunnel and the road at K35 + 335K35 + 500 to move, as well the tunnel lining to break or deformation seams to join together. In addition, the tunnel was always damaged after a new panel excavation. A delay of 35 months resulted, which was the time required for rebalancing the rock.

The tunnel was established in a coal-mine area and crosses over the mine’s old goaf, with the support type mostly IV or V. Poor supporting geological conditions comprise another reason for the tunnel damage.

Third, damage occurred three times during the rainy season or after a heavy rain, showing that the rain changed the original balance of the rock and induced the tunnel cracking. Rain can reduce tunnel strength. In addition, after the rainy season, the rock will shrink owing to water loss, which also induces tunnel deformation. Water is a factor that induces tunnel damage.

Since the tunnels are also a highway road, vehicle dynamic loads have a long-term influence on the structure. Since the tunnel crosses over a goaf, vehicle dynamic load can change the mechanical equilibrium of the broken rock mass. This, in turn, will cause tunnel deformation and additional stress. Vehicle dynamic load also induces tunnel damage.

To prevent deformation and harm from reoccurring, the coal mine has been closed. More protective measures, such as reinforcing the goaf, are needed.

The numerical simulation results were used to analyse the crack caused by coal mine. It shows that the main reason to cause the tunnel crack is the long-term, small production excavation. The coal mine was forced to stop production after the tunnel cracked. And the tunnel was maintained. And all the goaf was reinforced with the grouting filling method. The simulation results provide a technical reference for the implementation of those measures.

7. Conclusions

The following conclusions are drawn from the results of the present study:

(1)To study the law of tunnel damage caused by mine excavation, a model was built based on the geological and topographic information. The surface and tunnel movement were consistent with the general mine subsidence law of surface and rock-mass deformation.(2)The tunnel stresses , , , , and were analysed, and results show that the stress-concentration area coincides with the damage position of the tunnel. Rock movement-induced subsidiary stress concentration is the major reason for the tunnel damage.(3)Old goaf deformation, poor geological conditions, water, and vehicle dynamic load also influenced the tunnel deformation and compromised its deformation resistance. Therefore, additional measures are needed to prevent the reoccurence of traffic hazards.

Data Availability

The data used to support this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors would like to thank the China Railway Fifth Survey and Design Institute Group Co., Ltd. for providing the basic geological data and pictures of BYS tunnel damage. This work was supported by the Fundamental Research Funds for the Central Universities (grant nos. 2020YQDC05 and 2021YQDC01), Ecological-Smart Mines Joint Research Fund of Natural Science Foundation of Hebei Province (grant no. E2020402086), and funds from the Natural Science Foundation of Ningxia Hui Autonomous Region (grant nos. 2021AAC03426 and 2021AAC03429).

References

G. He, L. Yang, G. Ling, F. Jia, and D. Hong, Mine Subsidence, China University of Mining and Technology Press, Xuzhou, China, 1991, https://books.google.ca/books?id=nqrYAAAACAAJ.

J. Cheng, F. Liu, and S. Li, “Model for the prediction of subsurface strata movement due to underground mining,” Journal of Geophysics and Engineering, vol. 14, no. 6, pp. 1608–1623, 2017.

A. M. S. Iwanec, J. P. Carter, and J. P. Hambleton, “Geomechanics of subsidence above single and multi-seam coal mining,” Journal of Rock Mechanics and Geotechnical Engineering, vol. 8, no. 3, pp. 304–313, 2016.

E. F. Salmi, M. Nazem, and M. Karakus, “Numerical analysis of a large landslide induced by coal mining subsidence,” Engineering Geology, vol. 217, pp. 141–152, 2017.

B. Ghabraie, G. Ren, X. Zhang, and J. Smith, “Physical modelling of subsidence from sequential extraction of partially overlapping longwall panels and study of substrata movement characteristics,” International Journal of Coal Geology, vol. 140, pp. 71–83, 2015.

L. Tong, L. Liu, Y. Qiu, and S. Liu, “Tunneling in abandoned coal mine areas: problems, impacts and protection measures,” Tunnelling and Underground Space Technology, vol. 38, pp. 409–422, 2013.

J. Zheng, G. Zhang, X. Li, L. Yang, and Y. Dong, “The technology of cutting through coal seam and preventing collapse in the right line of huayingshan highway tunnel,” Chinese Journal of Geotechnical Engineering, vol. 22, no. 5, pp. 559–561, 2000.

D. Jiang, “Study on stability of tunnel surrounding rock crossing coal gob strata,” Chongqing Jiaotong University, Chongqing, China, 2017, M.A Thesis.

L. Yong, “Research on the strength characteristics and macro mechanical parameters of jointed rock mass,” Beijing Jiaotong University, Beijing, China, 2013, Ph.D Thesis.

L. Peixian, Y. Lili, and Y. Dehua, “Study of tunnel damage caused by underground mining deformation: calculation, analysis, and reinforcement,” Advances in Civil Engineering, vol. 2019, Article ID 4865161, 18 pages, 2019.

K. M. Lee, X. Y. Hou, X. W. Ge, and Y. Tang, “An analytical solution for a jointed shield-driven tunnel lining,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 25, no. 4, pp. 365–390, 2001.

L. Chen, “Stability analysis of rock pillar in the uncovering coal process of highway tunnel though coal seam,” Chongqing University, Chongqing, China, 2009, Ph.D Thesis.

A. Mitelman and D. Elmo, “Modelling of blast-induced damage in tunnels using a hybrid finite-discrete numerical approach,” Journal of Rock Mechanics and Geotechnical Engineering, vol. 6, no. 6, pp. 565–573, 2014.

X. J. Han, X. Meng, X. Zhang, and Y. Zhang, “The deformation stability analysis of the tunnel in mined out areas based on creator and flac3d,” Journal of Water Resources and Architectural Engineering, vol. 12, no. 5, pp. 93–97, 2014.

S. Anubhav and H. Wu, “Modelling of non-linear shear displacement behaviour of soil-geotextile interface,” International Journal of Geosynthetics and Ground Engineering, vol. 1, no. 2, 2015.

G. Zheng, T. Cui, and X. Cheng, “Study of the collapse mechanism of shield tunnels due to the failure of segments in sandy ground,” Engineering Failure Analysis, vol. 79, pp. 464–490, 2017.

H. Zhang, A. C. Adoko, Z. Meng, H. Wang, and Y. Jiao, “Mechanism of the mudstone tunnel failures induced by expansive clay minerals,” Geotechnical and Geological Engineering, vol. 35, no. 1, pp. 263–275, 2017.

J. Yang, W. Lu, Y. Hu, M. Chen, and P. Yan, “Numerical simulation of rock mass damage evolution during deep-buried tunnel excavation by drill and blast,” Rock Mechanics and Rock Engineering, vol. 48, no. 5, pp. 2045–2059, 2015.

J. F. Labuz and A. Zang, “Mohr-coulomb failure criterion,” Rock Mechanics and Rock Engineering, vol. 45, no. 6, pp. 975–979, 2012.

J. Chen, X. Liu, and Q. Xu, “Numerical simulation analysis of damage mode of concrete gravity dam under close-in explosion,” KSCE Journal of Civil Engineering, vol. 21, no. 1, pp. 397–407, 2017.

S. Kuili and V. R. Sastry, “A numerical modelling approach to assess the behaviour of underground cavern subjected to blast loads,” International Journal of Mining Science and Technology, vol. 28, no. 6, pp. 975–983, 2018.

P. H. S. W. Kulatilake, Q. Wu, Z. Yu, and F. Jiang, “Investigation of stability of a tunnel in a deep coal mine in China,” International Journal of Mining Science and Technology, vol. 23, no. 4, pp. 579–589, 2013.

H. Xie, H. Zhou, J. Wang, L. Li, and M. A. Kwasniewski, “Application of FLAC to predict ground surface displacements due to coal extraction and its comparative analysis,” Chinese Journal of Rock Mechanics and Engineering, vol. 18, no. 4, pp. 397–401, 1999.

L. Huang, Y. Lu, D. Su, and D. Zhang, “Treatment technology of highway tunnel through steep inclined goaf,” Journal of Highway and Transportation Research and Development, vol. 29, no. 11, pp. 80–85, 2012.

C. Meifeng, H. Manchao, and L. Dongyan, Rock Mechanical and Engineer, China Science Press, Beijing, China, 2017.

L. Shiping, A Brief Notebook of Rock Mechanics, China University of Mining and Technology Press, Xuzhou, China, 1986.

E. Hoek and T. Brown, Underground Excavation in Rock, Austin and Sons Ltd, Hereford, UK, 2nd edition, 1988.

M. Gergi and A. tacheva E, “On the valuation of strength and resistance condition of the rock in natural rock mass,” in Proceedings of the 2nd Congress of the International Society for Rock Mechanics, pp. 365–374, Yugoslavian Science Press, Belgrade, Serbia, September 1970.