This paper deals with analysis of rock massif in connection with a planned
deposit of nuclear waste. This problem has emerged with renessaince of nuclear power
plants which are now considered as a stable and relatively clean source of electrical energy.
The deposits are located deep in rock massif and analysis of stress redistribution caused by
cave drifting and temperature changes is required. Difficulties with solution of systems of
non-linear algebraic equations were discovered and two variants of the arc-length method
were tested

Koudelka, Tomás

Kruis, Jaroslav

English

UPCommons. Portal del coneixement obert de la UPC

521IV International Conference on Computational Methods for Coupled Problems in Science and EngineeringCOUPLED PROBLEMS 2011M. Papadrakakis, E. Oñate and B. Schrefler (Eds)ANALYSIS OF ROCK MASSIF BASED ON THE THEORYOF DAMAGEJAROSLAV KRUIS∗ AND TOMÁŠ KOUDELKA∗∗Department of Mechanics, Faculty of Civil EngineeringCzech Technical University in PragueThákurova 7, 166 29 Prague, Czech Republice-mail: jk@cml.fsv.cvut.cz, koudelka@cml.fsv.cvut.czKey words: Damage, Heat Transfer, Coupled Problems, Non-linear Systems of Equa-tions, Arc-length MethodsAbstract. This paper deals with analysis of rock massif in connection with a planneddeposit of nuclear waste. This problem has emerged with renessaince of nuclear powerplants which are now considered as a stable and relatively clean source of electrical energy.The deposits are located deep in rock massif and analysis of stress redistribution caused bycave drifting and temperature changes is required. Difficulties with solution of systems ofnon-linear algebraic equations were discovered and two variants of the arc-length methodwere tested.1 IntroductionRenaissance of nuclear power plants in many countries is accompanied with solutionof problem with nuclear waste. One method considered is the excavation of large under-ground caves where the waste can be stored with the help of very ingenious technology.The artificially made caves are located deep under the surface in rock massif. This topicis in the center of attention of many research group [1, 2, 3].The final analysis of deposit of nuclear waste has to contain not only mechanical partbut also transport processes have to be studied in connection with possible leakage. There-fore, hydro-thermo-mechanical analysis will be required. The main emphasis will be laidon the coupling of damage parameters with material diffusivity and permeability. In thispreliminary study, only one way coupling, where the heat transfer influences the mechan-ical analysis, is taken into account.Special attention has to be devoted to the solver of non-linear algebraic equations whichare obtained by discretisation of the problem solved by the finite element method. Thearc-length method is used in this paper because of the softening branch of the stress–straindiagrams. Some difficulties were observed and the classical spherical arc-length methodwas compared with the linearized method.1522Jaroslav Kruis and Tomáš KoudelkaThreedimensional model of the whole rock massif with excavated caves is usually verylarge, i.e. it contains at least hundreds of thousands degrees of freedom. Such large prob-lems are hardly solvable on single processor computers. In fact, there are two possibilitieshow to solve these problems efficiently. One is based on domain decomposition meth-ods which can be performed on parallel computers [4, 5, 6]. Second possibility is basedon adaptive methods where the mesh and the degree of approximation polynomials arechanged with respect to error estimates [7].The paper is organized as follows. Section 2 describes very briefly material models used.Section 3 summarizes variants of the arc-length method. Section 4 deals with numericalresults.2 Material ModelsAs was mentioned before, this contribution deals with preliminary analysis of rockwhere the mechanical load cases are accompanied by thermal load. The mechanical anal-ysis is based on the damage theory. The simplest isotropic damage model was not usedwith respect to different strength in tension and compression. The used orthotropic dam-age model is described in the next subsection.2.1 Orthotropic Damage ModelThe main drawback of the scalar isotropic damage model is that it uses only onedamage parameter for all principle directions regardless of tension or compression. Oncethe damage parameter caused by exceeding limit strain in one principle direction evolves,it reduces stiffness in all remaining principle directions even though they should not beinfluenced. This drawback is not significant in the case of the one-dimensional stress statesuch as pure bending but it becomes more important especially for the three-dimensionalstress state.That led to development of the more advanced damage model which can describebetter the 3D problems. In reference [8], the authors proposed general anisotropic modelfor concrete which contains nine material parameters. The laboratory measurements ofthe required material parameters has to be performed but it caused difficulties for certaincases. Additionally, the model required a significant number of internal variables thathave to be stored. These difficulties led to development of a simplified version of themodel which is based on six material parameters - three for tension and another threeparameters for compression.The model is based on the following stress-strain relationσα = (1 − H(εα)Dtα − H(−εα)Dcα)[(3K − 2G)εv + 2Gεα], (1)where the index α stands for the index of principle components of the given quantity. Themodel defines two sets of damage parameters Dtα and Dcα for tension and compression,respectively. In the equation (1), the symbol H() denotes the Heaviside function, K isthe bulk modulus, G is the shear modulus and εv stands for volumetric strain.2523Jaroslav Kruis and Tomáš KoudelkaThere are many evolution laws that can be used for Dt and Dc description. In ourproblems, the two evolution laws for the damage parameters are used similarly to thelaws used in the scalar isotropic damage model. The first law gives better results forcompression but the determination of the material paremeters is more complicated. Itcan be written in the formDβα =Aβ(|εβα| − εβ0 )Bβ1 + Aβ(|εβα| − εβ0 )Bβ, (2)where β represents indices t or c which are used for tension and compression. Aβ, Bβare material parameters controlling the peak value and slope of the softening branch andεβ0 is the strain threshold. The second law involves correction of the dissipated energywith respect to the size of elements and it describes tension better. It is defined by thenon-linear equation (3) which can be solved using the Newton method(1 − Dβα)E|εβα| = fβ exp(−Dβαh|εβα|wβ0). (3)In the above equation, fβ represents the tensile or compressive strength and wβ0 controlsthe initial slope of the softening branches. More details about the implemented modelscan be found in [9, 10, 11].2.2 Models for Transport ProcessesWith respect to limited space, detailed description of material models for transportprocesses are not included in this paper. They can be found e.g. in references [12, 13, 14].Efficient computer implementation of transport processes can be found in [15].3 Arc-length MethodWith respect to the softening part of stress–strain diagram, the systems of non-linearalgebraic equations are solved by one of arc-length methods. Detailed description of themethod can be found in references [16, 17].The equilibrium condition of a structure after discretization by the finite elementmethod has the formf int(d) = f c + λf p (4)where d denotes the vector of nodal displacements, f int denotes the vector of internalforces, f c denotes the vector of constant prescribed forces, λf p denotes the vector ofproportionally changing forces and λ denotes the scalar load-level multiplier. The vectorof unbalanced forces has the formr(d, λ) = f c + λf p − f int(d) (5)3524Jaroslav Kruis and Tomáš Koudelkaand it is the residual. The relationship between the forces f and the displacements dhas to be obtained by an incremental process where increments of length along the curvef—d are used. Within each increment, iteration process is needed in order to attain theequilibrium state. Let the i-th increment be known, i.e. the vector di and the parameterλi are known and r(di, λi) = 0. Expansion of the residual has the formr(di+1, λi+1) = r(di, λi) +∂r(di, λi)∂dδdi +∂r(di, λi)∂λδλi = −Ki,0δdi,1 + f pδλi,1 = 0 (6)Let the vector δdi,1 be in the formδdi,1 = δλi,1vi,1 (7)The length of arc in the first iteration within the increment can be written(δdi,1)T δdi,1 + ψ2(δλi,1)2fTp f p = (δλi,1)2vTi,1vi,1 + ψ2(δλi,1)2fTp f p = (∆l)2 (8)where the scaling parameter ψ was defined. The increment of the scalar load multiplierhas the formδλi,1 = ± ∆l√vTi,1vi,1 + ψ2fTp f p(9)New system of equations has to be solved in second iterationf c + (λi + δλi,1)f p − f int(di + δdi,1) − Ki,1δdi,2 + f pδλi,2 = 0 (10)Cumulative quantities are defined∆di,j = ∆di,j−1 + δdi,j (∆di,1 = δdi,1) (11)∆λi,j = ∆λi,j−1 + δλi,j (∆λi,1 = δλi,1) (12)and they are schematically depicted in Figure 1. Equation (10) can be split into twosystemsKi,1ui,2 = f c + (λi + ∆λi,1)f p − f int(di + ∆di,1) (13)Ki,1vi,2 = f p (14)where the previously defined notation is used. The length of arc has now the form‖∆di,1 + ui,2 + δλi,2vi,2‖2 + ψ2‖∆λi,1f p + δλi,2f p‖2 = (∆l)2 (15)which is the quadratic equationa1(δλi,2)2 + a2(δλi,2) + a3 = 0 (16)4525Jaroslav Kruis and Tomáš Koudelkafdf(d)dif c + λifp∆lδdi,1∆di,1δλi,1fpδdi,2∆di,2Figure 1: Non-linear force–displacement relationship.with coefficientsa1 = vTi,2vi,2 + ψ2fTp f p (17)a2 = 2vTi,2(∆di,1 + ui,2) + 2∆λi,1ψ2fTp f p (18)a3 = (∆di,1 + ui,2)T (∆di,1 + ui,2) + (∆λi,1)2ψ2fTp f p − (∆l)2 (19)The increment δλi,2 is obtained from the quadratic equation (16). New values are againsubstituted to the residual and equality to the zero vector is checked. The algorithm issummarized in Table 1 and it is called the spherical arc-length method. If the scalingparameter ψ is equal to zero, the method is called the cylindrical arc-length method.Solution of the quadratic equation (16) is straightforward but only one root has to beused for next computation. One of the criteria used has the formcos θ =∆dTi,j+1∆di,j(∆l)2→ max (20)Substitution of (11) and (12) leads to the formcos θ =1(∆l)2∆dTi,j(∆di,j + ui,j+1 + δλi,j+1vi,j+1) (21)Both roots of the equation (16) are substituted to the expression (21) and the root leadingto the larger value is selected.Linearized form of the arc-length leads to the expressionδλi,j+1 =−12li,j − ∆dTi,jui,j+1∆dTi,jvi,j+1 + ψ2∆λi,jfTp f p(22)and no root selection procedure is needed.5526Jaroslav Kruis and Tomáš KoudelkaTable 1: Algorithm of the Arc-length Methodλ0 = 0, d0 = 0For i = 0, 1, 2, . . .∆λi,0 = 0, ∆di,0 = 0, ri,0 = 0For j = 0, 1, 2, . . .ui,j+1 = K−1i,j ri,jvi,j+1 = K−1i,j f pa1 = vTi,j+1vi,j+1 + ψ2fTp f pa2 = 2vTi,j+1(∆di,j + ui,j+1) + 2∆λi,jψ2fTp f pa3 = ‖∆di,j + ui,j+1‖2 + (∆λi,j)2ψ2fTp f p − (∆l)2a1(δλi,j+1)2 + a2(δλi,j+1) + a3 = 0 ⇒ δλi,j+1δdi,j+1 = ui,j+1 + δλi,j+1vi,j+1∆di,j+1 = ∆di,j + δdi,j+1∆λi,j+1 = ∆λi,j + δλi,j+1ri,j+1 = f c + (λi + ∆λi,j)f p − f int(di + ∆di,j)if ‖ri,j+1‖ < ε, stopλi+1 = λi + ∆λidi+1 = di + ∆di6527Jaroslav Kruis and Tomáš KoudelkaFigure 2: Difficulties with non-linear solvers in softening branch.4 Numerical resultsTwo problems were numerically modelled in this preliminary study. First problem wasrepresented by a compression of cylinder which was used for material parameter identi-fication. Cylindrical rock samples were tested in laboratory [1] and the Young modulus,stress-strain relationships and the compressive strength were obtained. Difficulties wereobserved during the numerical analysis in softening branch of stress-strain relationship.The spherical arc-length method converges very slowly while the linearized arc-lengthmethod exhibits incorrect unloading–loading cycles. The situation is depicted in Figure2 where the red curve is obtained by the spherical method and the black curve by thelinearized method.Second problem was represented by a rock massif with excavated caves. Two differentcaves were taken into account. Distribution of the damage parameter is depicted for bothcases in Figure 3.5 ConclusionsAnalysis of rock massif based on the orthotropic damage material model was performed.Load cases were represented by the dead load and thermal load. Identification of materialparameters was done with the help of laboratory tests. Difficulties with solution of non-linear systems of equations were observed.7528Jaroslav Kruis and Tomáš KoudelkaFigure 3: Distribution of the damage parameter in the vicinity of caves.Acknowledgement:Financial support for this work was provided by project number P105/10/1682 of CzechScience Foundation. The financial support is gratefully acknowledged.REFERENCES[1] Staub, I. and Andersson, J.C. Äspö Pillar Stability Experiment. Swedish NuclearFuel and Waste Management Co., Sweden, (2004).[2] Blaheta, R. and Byczanski, P. and Kohut, R. and Kolcun, A. and Šňupárek, R. Large-scale modelling of T-M phenomena from underground deposition of the spent nuclearfuel. In EUROCK 2005: Impact of Human Activity on the Geological Environment,pages 49-55, (2005).[3] Blaheta, R. and Byczanski, P. and Jakl, O. and Kohut, R. and Kolcun, A. andKrečmer, K. and Starý, J. Large scale parallel FEM computations of far/near stressfield changes in rocks. Future Generation Computer Systems-The International Jour-nal of Grid Computing Theory Methods and Applications, (2006) 22:449–459.[4] Kruis, J. Domain Decomposition Methods for Distributed Computing. Saxe-CoburgPublications, Kippen, Stirling, Scotland, (2006).[5] Kruis, J. Domain Decomposition Methods on Parallel Computers. In Topping, B. H.V. and Mota Soares, C. A., editors, Progress in Engineering Computational Technol-ogy, pages 299–322, Saxe-Coburg Publications, Stirling, Scotland, UK, (2004).[6] Kruis, J. The FETI Method and its Applications: A Review. In Topping, B.H.V.and Iványi, P., editors, Parallel, Distributed and Grid Computing for Engineering,pages 199-216, Saxe-Coburg Publications, Stirling, Scotland, UK, (2009).[7] Šoĺın, P. and Segeth, K. and Doležel, I. Higher-Order Finite Element Methods. Chap-man & Hall/CRC, USA, (2004).8529Jaroslav Kruis and Tomáš Koudelka[8] Papa, E. and Taliercio, A. Anisotropic Damage Model for the Multiaxial Static andFatigue Behaviour of Plain Concrete. Eng. Fract. Mech. (1996) V. 55, N. 2 :163–179.[9] Koudelka, T. and Krejč́ı, T. An Anisotropic Damage Model for Concrete in Cou-pled Problems. In Topping, B.H.V. and Papadrakakis, M., editors, Proceedings ofthe Ninth International Conference on Computational Structures Technology, Civil-Comp Press, Stirlingshire, UK, (2008), paper 157.[10] Krejč́ı, T. and Koudelka, T. and Šejnoha, J. and Zeman, J. Computer Simulation ofConcrete Structures subject to Cyclic Temperature Loading. In Topping, B.H.V. andCosta Neves, L.F. and Barros, R.C., editors, Proceedings of the Twelfth InternationalConference on Civil, Structural and Environmental Engineering Computing, Civil-Comp Press, Stirlingshire, UK, (2009), paper 131.[11] Koudelka, T. and Krejč́ı, T. and Šejnoha, J. Analysis of a Nuclear Power PlantContainment. In Topping, B.H.V. and Costa Neves, L.F. and Barros, R.C., editors,Proceedings of the Twelfth International Conference on Civil, Structural and Environ-mental Engineering Computing, Civil-Comp Press, Stirlingshire, UK, (2009), paper132.[12] Schrefler, B.A. and Lewis, R.W. The Finite Element Method in the Static and Dy-namic Deformation and Consolidation of Porous Media. 2nd Edition, John Wiley &Sons, (1998).[13] Šejnoha, J. and Bittnar, Z. and Krejč́ı, T. and Kruis, J. Computational Aspects inThermo-Hydro-Mechanical Analysis of Porous Media, Part I: Transport Processes.In Topping, B. H. V. and Bittnar, Z., editors, Computational Structures Technology,pages 153–182, Civil-Comp Press Ltd, Stirling, Scotland, UK, (2002).[14] Černý, R. and Maděra, J. and Koč́ı, J. and Vejmelková, E. Heat and moisture trans-port in porous materials involving cyclic wetting and drying. In Fourteenth Interna-tional Conference on Computational Methods and Experimental Measurements, pages3–12, WIT Press, Wessex, UK, (2009).[15] Kruis, J. and Koudelka, T. and Krejč́ı, T. Efficient computer implementation of cou-pled hydro-thermo-mechanical analysis. Mathematics and Computers in Simulation(2010) 80:1578–1588.[16] Crisfield, M.A. Non-linear Finite Element Analysis of Solids and Structures. JohnWiley & Sons Ltd, Chichester, UK, (1991).[17] Bittnar, Z. and Šejnoha, J. Numerical Methods in Structural Mechanics. ASCE Press,New York, USA, (1996).9

Analysis of rock massif based on the theory of damage

https://upcommons.upc.edu/bitstream/2117/327284/1/Coupled-2011-45_Analysis%20of%20rock%20massif.pdf

Analysis of rock massif based on the theory of damage

Abstract

Similar works

Full text

Available Versions

UPCommons. Portal del coneixement obert de la UPC