This work presents a numerical model to simulate intergranular damage in a granitic rock by oscillating piezoelectric excitation of quartz dispersed in the structure. The damage evolution at grain boundaries was assumed to be related to fatigue and it was modelled using cohesive elements and a damage evolution model formulated in terms of discipation of mechanical work. An explicit representation of the granular mesostructure was built, and it was subjected to high-voltage alternating-current exitation. The effect of the fatigue damage on the mechanical properties was quantified by simulated tension and compression tests. The numerical results show that the electrical treatment can potentially cause rock weakening due to fatigue, but the model needs to be calibrated with experimental data for a quantitative analysis.Peer reviewe

Hokka, Mikko

Kane, Pascal-Alexandre

Rubio Ruiz, Arturo

Saksala, Timo

Trepo - Institutional Repository of Tampere University

ABSTRACT: This work presents a numerical model to simulate intergranular damage in a granitic rock by oscillating piezoelectric excitation of quartz dispersed in the structure. The damage evolution at grain boundaries was assumed to be related to fatigue and it was modelled using cohesive elements and a damage evolution model formulated in terms of discipation of mechanical work. An explicit representation of the granular mesostructure was built, and it was subjected to high-voltage alternating-current exitation. The effect of the fatigue damage on the mechanical properties was quantified by simulated tension and compression tests. The numerical results show that the electrical treatment can potentially cause rock weakening due to fatigue, but the model needs to be calibrated with experimental data for a quantitative analysis. Keywords: Piezoelectric actuation, High voltage, High frequency, Alternating current, Rock fracture, Damage cyclic evolution. 1 INTRODUCTION Deep drilling operations in hard rock mines are expensive mainly because of the wear of drilling tools and low rate of penetration. There is extensive work in the search for new methods to induce damage in rocks using non-mechanical-contact techniques. A potential solution to weaken granitic formations is to induce damage by excitation of dispersed piezoelectric phases (quartz) using high voltage alternating current (HV-AC) (Saksala, 2021). A possible reason for such rock weakening is the nucleation of fatigue cracks at grain boundaries. This work aims at improving the understanding of this newly observed multiphysical phenomena. A 2D finite element (FE) model with piezoelectric coupling was developed to reproduce the phenomenon. The virtual mesostructure of the rock was built by importing Electron Back Scatter Diffraction (EBSD) data into the model. The intragranular mechanical behavior was isotropic linear elastic, while the load transfer between grains and the damage evolution at grain boundaries was described using cohesive elements (CE). A cohesive zone model with continuous time evolution of damage was implemented to describe the fatigue propagation at grain boundaries. The effects of the fatigue damage due to the HV-AC treatment on 15th ISRM Congress 2023 & 72nd Geomechanics Colloquium. Schubert & Kluckner (eds.) © ÖGG   A numerical analysis of weakening of a granitic rock by piezoelectric excitation of quartz Rafael Arturo Rubio Ruiz Faculty of Engineering and Natural Sciences, Tampere University, Tampere, Finland Timo Saksala Faculty of Built Environment, Tampere University, Tampere, Finland Alexandre Kane Sintef, Trondheim, Norway Mikko Hokka Faculty of Engineering and Natural Sciences, Tampere University, Tampere, Finland -1852-the mechanical behavior of the material was estimated using virtual uniaxial tension and compression tests.  2 THEORY AND MODELLING FRAMEWORK The studied phenomena consist of the cyclic excitation of the piezoelectric properties of quartz grains dispersed in the polycrystalline structure of a granitic rock. The hypothesis is that such cyclic excitation induces damage at grain boundaries, while the interior of the grains remains undamaged. A synthetic microstructure was generated to explicitly represent the polycrystalline structure in a FE model with CE to describe the load transfer and damage evolution at grain boundaries. Hereafter are presented the global governing equations of the piezoelectric-mechanical problem, the CE formulation with the constitutive traction separation relationship, and the numerical approach. 2.1 Global piezoelectric-mechanical problem The governing Equations of the piezoelectric-mechanical problem in matrix form are the following  𝐌𝐌?̈?𝒖 + 𝐂𝐂?̇?𝒖 + 𝐊𝐊𝒖𝒖𝒖𝒖𝒖𝒖 + 𝐊𝐊𝒖𝒖𝒖𝒖𝒖𝒖+ 𝐟𝐟ch,𝒖𝒖 = 𝐟𝐟𝒖𝒖,ext (1)  𝐊𝐊𝒖𝒖𝒖𝒖𝒖𝒖 + 𝐊𝐊𝒖𝒖𝒖𝒖𝒖𝒖 + 𝐟𝐟ch,𝒖𝒖 = 𝐟𝐟𝒖𝒖,ext, with (2)  𝐌𝐌 = 𝑨𝑨e=1nbe  �𝑵𝑵𝒖𝒖e,Tρ𝑵𝑵𝒖𝒖ed𝛀𝛀𝐞𝐞 ,𝑲𝑲𝒖𝒖𝒖𝒖 =  𝑨𝑨e=1nbe  �𝑩𝑩𝒖𝒖e,T𝐂𝐂𝐞𝐞𝑩𝑩𝒖𝒖ed𝛀𝛀𝐞𝐞,   𝐂𝐂 = 𝜉𝜉𝐌𝐌, (3)   𝑲𝑲𝒖𝒖𝒖𝒖 =  𝑨𝑨e=1nbe  �𝑩𝑩𝒖𝒖e,T𝛆𝛆𝑩𝑩𝒖𝒖e d𝛀𝛀𝐞𝐞 ,𝐊𝐊𝒖𝒖𝒖𝒖 = 𝑨𝑨e=1nbe  �𝑩𝑩𝒖𝒖e,T𝐞𝐞𝑩𝑩𝒖𝒖e d𝛀𝛀𝐞𝐞 ,𝐊𝐊𝒖𝒖𝒖𝒖 =  𝐊𝐊𝒖𝒖𝒖𝒖T  (4) were 𝑨𝑨e=1nbe  is the assembly operator over the bulk elements of the finite element mesh. 𝒖𝒖, ?̇?𝒖 and  ?̈?𝒖 are the displacement, velocity and acceleration vectors, while 𝒖𝒖 is the electric potential vector. 𝑩𝑩𝒖𝒖e  and 𝑩𝑩𝒖𝒖e  are the gradient of the interpolation matrices of displacement and electric potential, 𝑵𝑵𝒖𝒖e  and 𝑵𝑵𝒖𝒖e , respectively. 𝝔𝝔 and 𝒕𝒕� are respectively the electric field and the traction vector over the boundary of the computational domain, 𝚪𝚪. The constant 𝜉𝜉 is the proportionality constant between the lumped mass matrix, 𝐌𝐌, and the damping matrix, 𝐂𝐂. 𝐟𝐟ch,𝒖𝒖 and 𝐟𝐟𝒖𝒖,ext represent the internal cohesive forces and the flux of electric potential through cohesive elements, both described in Section 2.2. ρ is the mass density, while 𝐂𝐂𝐞𝐞, 𝛆𝛆 and 𝐞𝐞 are the elastic, dielectric and piezoelectric tensors. The mechanical behavior of the interior of the constituent minerals was isotropic linear elastic and plain stress was assumed in the 2D model. The constituent minerals were modelled as dielectrically isotropic, so the form of the dielectric tensor reduced to ε𝐈𝐈, were ε and 𝐈𝐈 stand for the dielectric constant and the identity matrix. The tensor 𝐞𝐞 is the result of the product between the elastic tensor 𝐂𝐂𝐞𝐞 and the piezoelectric tensor d, which has nonzero entries only in the quartz grains as,   𝐝𝐝 = �d11 0 −d110 −2d11 0� (5) -1853-2.2 Cohesive elements formulation 2.2.1 Cohesive elements mechanical formulation The finite element model presented in this work is a 2D simplification of the mesostructure. Thus, the interphases can be described using four-noded cohesive elements. The relation between nodal displacements and the cohesive forces are defined in terms of local normal and tangential coordinates of the CE. The local rotation matrix R is defined by the direction cosines of the normal of each CE and it transforms the nodal displacements from global to local coordinates.  The rotated vector of nodal displacements is related to the continuous interphase separation vector by the matrix 𝐇𝐇 which computes the relative nodal displacements and interpolates them along the cohesive element using the shape functions as,  𝜟𝜟(𝐱𝐱) = 𝐇𝐇𝒖𝒖� = 𝐇𝐇𝐇𝐇𝒖𝒖 = 𝐁𝐁ch𝒖𝒖 (6) Then the internal vector of cohesive forces in Equation (1) is calculated as,  𝐟𝐟ch,𝒖𝒖 = 𝑨𝑨che=1nche �𝐁𝐁ch,T𝒕𝒕ch(𝚫𝚫)dΓche (7) were 𝑨𝑨che=1nche  is the assembly operator over the cohesive elements of the finite element mesh and 𝒕𝒕ch is the traction vector at the interphase, which is determined by the traction separation relationship explained in Section 2.2.2. The presence of CE causes artificial compliance, which can be diminished by increasing the initial stiffness of the CE, which leads to a reduction of the critical time step in explicit integration schemes. Thus, the bi-penalty formulation presented by Hetherington and Askes (2014) was used avoid the reduction of the critical time step due to the presence of interphases. The mass bi-penalty method consists of adding a mass penalization matrix (Mp) to the global mass matrix in Equation (1). In this work the matrix Mp was defined with a proportionality of the tangent stiffness matrix of the cohesive element (Kch) obtained by linearization of the constitutive traction separation relationship explained in Section 2.2.2,  𝐌𝐌𝐩𝐩 = 𝑨𝑨che=1nche1R �𝑩𝑩ch,T𝑪𝑪ch𝑩𝑩ch𝑑𝑑Γche =1R𝐊𝐊ch, (8) with R being a proportionality constant that needs to meet the condition R < 4 Δ𝑡𝑡2⁄ , with Δt the critical time step of the problem without interphases. 𝑪𝑪ch is the linearized rigidity of the cohesive element. One important feature of Mp is that its entries add to zero, thus the penalization does not add any effective mass to the system. 2.2.2 Constitutive equations and fatigue dame evolution  In this work the constitutive equation relating interphase openings (𝚫𝚫) and the tractions at the interphases (𝒕𝒕ch) was the bilinear cohesive law, which relies on a scalar damage variable (D) and a scalar effective opening (Δeff), computed as,  𝐷𝐷 = 𝛥𝛥eff −  𝛥𝛥0eff𝛥𝛥𝑐𝑐eff −  𝛥𝛥0eff, with   𝛥𝛥eff =  �𝛥𝛥n2 + 𝛽𝛽2𝛥𝛥τ2 (9) were Δn and Δτ are, respectively, the normal and tangential components of 𝚫𝚫. The parameters 𝛥𝛥0eff and 𝛥𝛥𝑐𝑐eff are the elastic separation limit and the maximum separation before complete decohesion. The cyclic damage accumulation in the model is controlled by a damage evolution model formulated in terms of mechanical work dissipation. Then, the damage can evolve even though the effective separation is smaller than 𝛥𝛥0eff or the maximum recorded effective separation (𝛥𝛥max), -1854- ?̇?𝐷 =  κΔ̇eff(〈𝑡𝑡𝑛𝑛〉2 + 𝑡𝑡𝜏𝜏2)12Gch  (10) were κ is a material parameter determining the fraction of mechanical work that is dissipated by damage evolution during loading. The notation 〈∗〉 indicates the Macaulay brackets, which are used in the normal component of the traction vector to ensure that compressive loads do not contribute to the damage evolution. Gch is the cohesive fracture energy, which can be approximated as the area under the inelastic part of the traction-separation curve under monotonic tension. For a bilinear traction law, the parameter Gch equals 𝑡𝑡𝑐𝑐𝛥𝛥𝑐𝑐eff 2⁄ , with 𝑡𝑡𝑐𝑐 the yield stress of the cohesive element. Notice that the evolution of D in Equation (10) also causes the evolution of 𝛥𝛥max and its value is obtained by equating  𝛥𝛥efffrom Equation (9) using the updated value of D. The components of 𝒕𝒕ch depend on 𝛥𝛥eff, 𝛥𝛥max, the sign of Δn and the accumulated damage D according to,  𝒕𝒕ch = (1 − 𝐷𝐷) �𝛼𝛼 00 𝛼𝛼�𝚫𝚫, if   𝛥𝛥max ≤ 𝛥𝛥0eff  (11)  𝒕𝒕ch = 𝑡𝑡𝑐𝑐〈1 − 𝐷𝐷〉𝛥𝛥eff𝚫𝚫, if   𝛥𝛥max > 𝛥𝛥0eff and 𝛥𝛥eff = 𝛥𝛥max (12)  𝒕𝒕ch = 𝑡𝑡𝑐𝑐〈1 − 𝐷𝐷〉𝛥𝛥max𝚫𝚫, if   𝛥𝛥max > 𝛥𝛥0eff and 𝛥𝛥eff < 𝛥𝛥max  (13) were 𝛼𝛼 is the initial stiffness of the cohesive element. If in any of the previous loading cases 𝛥𝛥n  is negative, then the normal component of 𝒕𝒕ch (tnch) equals 𝛼𝛼𝛥𝛥n. 2.2.3 Cohesive elements electromechanical coupling The flux of electrical potential through the CE (Dch) was defined according to the work of Kozinov (2014), who described Dch as,  Dch =  −𝜀𝜀⟦𝜑𝜑⟧(𝚫𝚫 ∙ 𝐧𝐧 + δnat) (14) were n and ⟦𝜑𝜑⟧ are the normal vector and the jump of electrical potential across the CE, while δnat is the finite minimum thickness of the grain boundary to avoid singularities when 𝚫𝚫 ∙ 𝐧𝐧 approaches zero. The calculation of ⟦𝜑𝜑⟧ and 𝐟𝐟ch,𝒖𝒖 follow a similar procedure as the calculation of 𝚫𝚫 and 𝐟𝐟ch,𝒖𝒖 described in Section 2.2.1. 2.2.4 Computational approach and simulation set up The FE model presented in this work contains an explicit representation of mesostructure of a rock with 321 grains based on the example EBSD data on Mylonite from the MTEX MATLAB code (Bachmann et al., 2011), which is shown in Figure 1a. The constituents of Mylonite are Quartz, Biotite, Andesine, and Orthoclase. Both Andesine and Orthoclase are Feldspars, and therefore this rock is similar to Granite. The EBSD data for each grain included the position of the grain centroids, the phase identification and the Euler angles. A polygonal approximation of the mesostructure was generated using the position of the grain centroids as seeds for Voronoi Tessellations generated using the MATLAB Code Polymesher (Talischi et al., 2012), as shown in Figure 2b. The Euler angles of the quartz grains -1855-were used to rotate the piezoelectric tensor. The 2D polycrystalline structure was meshed using linear triangular elements. The elastic properties and densities of the constituent minerals and the material parameters governing the CE are shown in Table 1.  Different numerical rock samples were subjected to HC-AC treatments varying the value of the parameter 𝜅𝜅 to assess its impact on the rock weakening. During the HV-AC treatments one of the edges of the rock was grounded while the opposite edge was subjected to a harmonically varying voltage with a frequency of 300 kHz and a voltage amplitude of 50 kV. For all samples, 3000 cycles were simulated in 0.01 s. After the HV-AC treatments the rock samples were mechanically tested in tension and in compression to determine the effects of the fatigue damage on the mechanical behavior of the rock. The integration scheme used for the time integration of Equation (1) was the well-known forward explicit, while the Equation (2) was solved using a quasistatic scheme. Table 1. Mineral properties and CE parameters.  Mineral/Property E [GPa] 𝜈𝜈 𝜌𝜌 [kg/m3] 𝜀𝜀 [F/m] d11[pC/N] Quartz 100.9 0.06 2650 4.5 𝜀𝜀0 2.27 Biotite 102.8 0.21 3250 6.3𝜀𝜀0 0 Feldspar 74.4 0.31 2580 7.75𝜀𝜀0 0 𝜀𝜀0 [F/m] 𝑡𝑡𝑐𝑐 MPa] Gch[J/m2] 𝛽𝛽   8.854 × 10−12 8.0 40 0.16    3 NUMERICAL RESULTS The HV-AC treatment simulation set-up shown in Figure 1c was applied to rocks with different values of the parameter κ, ranging within 0–0.2. A sample with κ=0 was included in the analysis as an undamaged reference material. The obtained electrical potential field and the corresponding stress field at the maximum voltage of the electrical excitation are shown in Figure 2a and Figure 2b. Notice that the chosen excitation frequency is close to a natural frequency of the rock sample, causing larger deformations within the material, which is a phenomenon previously observed by Saksala (2021). In fact, several simulations were attempted at lower AC frequencies and the stress levels in the rock sample were 50% lower than those at a frequency close to a natural frequency of the sample. The extent of damage in the rock samples changed significantly when with different values of parameter 𝜅𝜅, which had different effects on the consequent tension and compression behavior. The tensile mechanical behavior of the treated rock changed notably with the HV-AC treatment and the effect of the parameter 𝜅𝜅 was very strong, as shown in Figure 2c. The tensile stiffness of the rock changed from 64 GPa for the undamaged rock to 5.5 GPa for the treated rock with 𝜅𝜅 equal to 0.2. Moreover, the tensile strength of the rock was reduced by 27% due to the HV-AC treatment. The compressive behavior was less sensitive to the HV-AC treatment than the tensile behavior, but the effects were still noticeable. The samples treated with 𝜅𝜅 equal to 0.01 and 0.02 had a very similar Figure 1. Mylonite EBSD data (a), polygonal approximation of the rock mesostructure (b), HV-AC treatment simulation set-up(c). 𝜙𝜙0 = 50 kV and 𝑓𝑓 = 300 kHz. -1856-behavior as the untreated rock. However, the samples with 𝜅𝜅 equal to 0.1 and 0.2 were weaker, with maximum reduction of 36% in stiffness and 21% reduction in strength, as shown in Figure 2d. 4 CONCLUSIONS A method to simulate the damage evolution in a granitic rock due to cyclic piezoelectric excitation of quartz was presented. The model included an explicit representation of the mesostructure of the rock, which was generated from EBSD data. The load transfer between grains and the intergranular fatigue damage was modeled using a damage evolution law based on dissipation of mechanical work. The HV-AC treatment of the rock samples showed that resonance might play an important role if the excitation frequency is close to a natural frequency of the treated sample.  The tensile behavior of the rock is more sensitive to the HV-AC damage than the compressive behavior, but in both cases the influence of the HV-AC treatment is noticeable. The extent of the damage and consequent detriment of the mechanical performance of the rock due to the HV-AC treatment is highly dependent on the value of the parameter 𝜅𝜅. Thus, the model needs to be calibrated with experimental data to draw quantitative conclusions.  REFERENCES Bachmann, F., Hielscher, R., & Schaeben, H. (2011). Grain detection from 2d and 3d EBSD data—Specification of the MTEX algorithm. Ultramicroscopy, 111(12), 1720–1733. https://doi.org/10.1016/j.ultramic.2011.08.002 Hetherington, J., & Askes, H. (2014). A mass matrix formulation for cohesive surface elements. Theoretical and Applied Fracture Mechanics, 69, 110–117. https://doi.org/10.1016/j.tafmec.2013.11.011 Kozinov, S., Kuna, M., & Roth, S. (2014). A cohesive zone model for the electromechanical damage of piezoelectric/ferroelectric materials. Smart Materials and Structures, 23(5), 055024. https://doi.org/10.1088/0964-1726/23/5/055024 Saksala, T. (2021). Cracking of granitic rock by high frequency-high voltage-alternating current actuation of piezoelectric properties of quartz mineral: 3D numerical study. International Journal of Rock Mechanics and Mining Sciences, 147, 104891. https://doi.org/10.1016/j.ijrmms.2021.104891 Talischi, C., Paulino, G. H., Pereira, A., & Menezes, I. F. M. (2012). PolyMesher: A general-purpose mesh generator for polygonal elements written in Matlab. Structural and Multidisciplinary Optimization, 45(3), 309–328. https://doi.org/10.1007/s00158-011-0706-z Figure 2. Electrical potential field (a) and stress field (b) both at maximum voltage of electrical excitation. Stress-strain curves from tension tests (c). Stress-strain curves from compression tests (d). -1857-

A numerical analysis of weakening of a granitic rock by piezoelectric excitation of quartz

https://trepo.tuni.fi/bitstream/10024/152614/1/ISRM2023_2.pdf

A numerical analysis of weakening of a granitic rock by piezoelectric excitation of quartz

Abstract

Similar works

Full text

Available Versions

Trepo - Institutional Repository of Tampere University