Weakening of Compressive Strength of Granite by Piezoelectric Actuation of Quartz Using High-Frequency and High-Voltage Alternating Current: A 3D Numerical Study

Piezoelectric excitation of quartz mineral phase in granite using high-frequency and high-voltage alternating current (HF-HV-AC) is a potential new weakening pretreatment in comminution of rock. The present study addresses this topic numerically by quantifying the weakening effect on the compressive strength of granite. For this end, a numerical method based on a damage-viscoplasticity model for granite failure under piezoelectric actuation is developed. The rock material is modelled as heterogeneous and isotropic. However, the piezoelectric properties of quartz are anisotropic. The governing global piezoelectro-mechanical problem is solved in a staggered manner explicitly in time. Numerical simulations predict that the weakening effect on compressive strength of granite is 10% with the excitation frequency and voltage of 274.4 kHz and 150 kV of the pretreatment. As the weakening effect takes place at a natural frequency of the numerical rock sample, the quartz content has only a slight effect on the frequency at which maximum weakening occurs. Moreover, the weakening effect depends strongly on the orientation of the quartz crystals. In a more practical application of simulating low-rate compression of a sphere-shaped rock sample, a weakening effect of 8% after the HF-HV-AC pretreatment was predicted. Substantial damage can be induced on granite by piezoelectric excitation by high-voltage and high-frequency alternating current. This pretreatment by piezoelectricity of quartz can be used to weaken granite samples before mechanical comminution. The frequency of the excitation needs to match one of the natural frequencies of the rock sample. A weakening effect of 10% on the compressive strength can be achieved at frequency of 274 kHz and voltage of 150 kV. Experimental validation of the theoretical/numerical predictions is needed. Substantial damage can be induced on granite by piezoelectric excitation by high-voltage and high-frequency alternating current. This pretreatment by piezoelectricity of quartz can be used to weaken granite samples before mechanical comminution. The frequency of the excitation needs to match one of the natural frequencies of the rock sample. A weakening effect of 10% on the compressive strength can be achieved at frequency of 274 kHz and voltage of 150 kV. Experimental validation of the theoretical/numerical predictions is needed.


Introduction
Comminution and excavation of rocks and ores suffer from low energy efficiency and excessive tool wear.Indeed, comminution processes use more than half of the total energy used in mines (Klein et al. 2018;Napier-Munn 2015).These shortcomings of the traditional comminution methods, based on mechanical breakage, have spurred intensive search for new, more energy efficient comminution methods.The focus of research has lately been on the unconventional techniques, where a nonmechanical agent, such as an electric pulse of a thermal jet, is applied either alone or as a rock weakening pretreatment prior to mechanical comminution (Aditya et al. 2017;Sefiu et al. 2020).As a pretreatment, such methods can mitigate tool wear, thus enhancing the cost efficiency of comminution.
Electricity is exploited in a class of unconventional comminution methods.There are two varieties of this approach, the first being the electro-pulse drilling/excavation based on the electric breakdown phenomenon (Cho et al. 2016;Shi et al. 2014;Walsh and Vogler 2020;Li et al. 2020), and the second, potential new method exploits the piezoelectric properties of quartz mineral present in rocks such as quartzite and granite.This new method exploits the converse piezoelectric effect by electric actuation of the quartz phase.
The piezoelectricity in rocks, illustrated in Fig. 1, has been studied to some extent experimentally (Parkhomenko 1971;Bishop 1981;Tuck et al. 1977;Ghomshei and Templeton 1989).However, the question whether quartz bearing rocks truly exhibit a piezoelectric effect, through the presence of a quartz texture (i.e. a trend in the orientation of quartz crystal c-axes), remains controversial.According to Parkhomenko (1971), the quartz bearing rocks do display a texture.Notwithstanding, Tuck et al. (1977) reached the opposite conclusion and considered the weak piezoelectricity to be a statistical effect or caused by isolated large single/few quartz grains.Nevertheless, the measurements by Ghomshei and Templeton (1989) do seem to suggest the existence of piezoelectric textures in these rocks.In any case, rocks are aggregates of difference minerals and grains and, thus, the piezoelectric effect in quartz bearing rocks is orders of magnitude weaker than that of a single quartz crystal (Parkhomenko 1971;Bishop 1981).The magnitude of the electric field required to cause 10 MPa of stress in a quartz crystal by converse piezoelectric effect is ~ 500 kV/ cm.As this value exceeds the electric breakdown strength of granites, it seems that exploiting piezoelectric phenomenon of quartz is not a feasible pretreatment method.
However, Saksala (2021) demonstrated by numerical simulations that actuation using alternating current (AC) of high frequency (HF) and high voltage (HV) leads to cracking of granite rock samples.More precisely, tensile cracks were induced on cylindrical (numerical) rock samples by sinusoidal AC excitation at the frequency of ~ 100 kHz and the amplitude of ~ 10 kV.In order to induce cracks, the frequency of the excitation needs to match one of the natural frequencies of the sample.Therefore, the fracture mechanism is related to the resonance phenomenon of the sample.Saksala et al. (2022) continued this work, in a conference paper, demonstrating by numerical modelling that the tensile strength of rock can be weakened by 12% with the HF-HV-AC pretreatment.
The present paper continues the theoretical-numerical research on this method and tests, for the first time, how much the compressive strength of a granite can be weakened by the HF-HV-AC pretreatment.For this end, a numerical method including a staggered method to solve the governing piezoelectro-mechanical equations at global level, and a damage-viscoplasticity model to describe the rock material failure, is developed.A more practical comminutionoriented problem of crushing a spherical rock sample by diametral compression is also addressed.It should be stressed that this study, as a forerunner in the field, is theoretical speculative in nature and, therefore, no experimental studies on the topic seem to exist.Hopefully, thereby, this study motivates experimentalists to test the feasibility of the HF-HV-AC pretreatment for comminution.

Theory of the Numerical Method
The numerical method to simulate rock failure under piezoelectric excitation is presented in this section.First, the piezoelectric rock failure description based on a damage-viscoplasticity model is given.Then, the finite element discretized piezoelectro-mechanical problem and the time stepping method to solve it are presented.The basics of piezoelectricity related to rocks can be found e.g. in Parkhomenko (1971).

Rock Constitutive and Failure Description
The rock is described here as an elastic (up to failure) isotropic and heterogeneous material.This means that the rock constituent minerals are taken as isotropic, but each mineral phase has its specific material properties.However, the piezoelectric properties are taken, as they truly are (see Fig. 1), anisotropic.The failure model for rock consists of a viscoplastic part and a damage part.It should be noted that the 3D embedded discontinuity approach used by Saksala (2021) and Saksala et al. (2022) cannot be used here as it cannot adequately describe the compressive failure.In the present approach, the (perfectly) viscoplastic part indicates the stress states leading to rock failure and inelastic deformation while the damage part quantifies the failure by separate isotropic damage variables in compression and tension.
The stress states leading to failure and inelastic flow are indicated by a bi-surface yield criterion consisting of the Hokka power law criterion (Hokka et al. 2016) for compressive (shear) fracture, called here the MH criterion, and the modified Rankine (MR) criterion for the tensile fracture.The MH criterion matches the experimental data for Kuru granite under confined dynamic compression better than the Hoek-Brown criterion (Hokka et al. 2016).Perfect viscoplasticity is assumed, as the damage model accounts for stiffness and strength degradation.These criteria are written as where the symbols are as follows: σ i is the ith principal stress of the stress tensor ; σ c , σ t are the dynamic uniaxial compressive and tensile strengths of the material; λMH and λMR are the rates of the internal variables in compression and tension, respectively; and s MH and s MR are the constant viscosity moduli in compression and tension, respectively.Moreover, the Macaulay brackets, i.e. the positive part operator, are used in the MR criterion while B and n in the MH criterion are experimentally calibrated parameters.The rate sensitivity of the model is provided by adding the product of viscosity and the rate of the internal variable to the static values, σ c0 , σ t0 , of the compressive and tensile strength.This method is reliable, as demonstrated by the experimental results by Zhao (2000) showing that the dynamic failure envelope of Bukit Timah granite can be obtained from the quasi-static one by translation.Equation (4) presents the loading-unloading conditions of Kuhn-Tucker form, which means that the present formulation is based on the viscoplastic consistency model by Wang et al. (1997).
The damage part is formulated with separate scalar damage variables in tension and compression due to the asymmetric behavior of rocks in these stress regions.Following Grassl and Jirásek (2006), the damage variables are driven by viscoplastic strain.Consequently, there is no need for additional damage loading functions, which simplifies the implementation considerably.Moreover, if damage was driven by total strain, spurious stress effects would occur in cyclic loading.The components of the damage model resulting in exponential softening can be written as The symbol meanings here are: A t and A c control the final value of the damage variables; β t and β c , controlling the initial slope and the amount of damage dissipation, are defined by the fracture energies G Ic and G IIc ; h e is a characteristic length of a finite element; vp eqvt is the equivalent viscoplastic strain in tension defined, in the rate form, as the trace of the viscoplastic strain rate tensor, ̇ vp , using the Macauley brackets so that tensile damage evolution occurs only if the volumetric viscoplastic principal strain is positive; vp eqvc is the equivalent viscoplastic strain in compression, being defined with the deviatoric part, ̇ vp , of ̇ vp .Moreover, Eq. ( 8) is the Koiter's rule for bi-surface plasticity.The colon in (7) denotes the double contraction operator for tensors, i.e.A ∶ A = A ij A ij .Finally, the plastic potential g MH in (9) proposed by Saksala  et al. (2017) with ψ is the dilation angle of the rock.
At this point, it is informative to illustrate the bi-surface yield criterion in the principal stress space.As the MH criterion does not depend on the intermediate stress space, the criteria are plotted in σ 1 -σ 3 stress space in Fig. 2.
It can be observed in Fig. 2 that, in the present case, the MH criterion actually matches the tensile strength, 10 MPa, in uniaxial tension quite well but clearly overestimates the bi-axial tensile strength.Moreover, the flow direction of the MH criterion would be wrong in uniaxial tension, i.e. it would deviate from the Rankine flow direction (see Fig. 2b), which is the correct one.Still more, if the tensile strength was lower, say 7 MPa, the MH criterion would overpredict it by more than 40%.Therefore, the modified (rounded) Rankine criterion is needed as a tensile cutoff.Figure 2b also illustrates the corner plasticity direction at the intersection of the yield surfaces.In this case, the flow direction, (5) 1 , m p = sin upon stress return mapping, is calculated by the Koiter's rule in Eq. ( 8), while the return mapping is performed by the standard cutting plane algorithm (Wang et al. 1997;Simo and Hughes 1998).The return mapping scheme is further elaborated in Appendix 1.
The final model components are the nominal-effective stress relation, which specifies how damage variables operate on the stress, and the constitutive law (within the small deformation framework) between the effective stress and the elastic strain: where is the nominal stress, i.e. the one returned onto the failure surface when the trial stress violates the criteria.Moreover, e is the elasticity tensor, and tot is the total strain.Finally, d is the third-order piezoelectric constants tensor (its 3 × 6 matrix form can be inferred from Fig. 1), and = −∇ e is the electric field with e being the scalar electric potential.
The present formulation combines the damage and viscoplasticity parts of the model in the effective stress space (Grassl and Jirásek 2006) allowing to separate the damage and viscoplasticity computations so that the stress return mapping is first performed in the effective principal stress space.Then, the damage variables are updated and, finally, the nominal stress is calculated by (11).This formulation requires that the damage variables are driven by the (visco)plastic strain. (10)

Finite Element Form of Piezoelectro-Mechanical Problem and Its Solution
The equations governing the HF-HV-AC pretreatment are the standard piezoelectro-mechanical equations, i.e. the balance of linear momentum and the piezoelectro-static balance, written in the finite element discretized form as (Allik and Hughes 1970) with In these equations, symbols are as follows: t is time; A is the finite element assembly operator; e ϕ ( = e ϕ e ) is the gradient of the electric potential interpolation matrix e ϕ ; e u is the gradient of displacement interpolation matrix e u ; M is the (lumped) mass matrix; C is the damping matrix with a constant α; is the electric charge; int t and ( 12) where I is the unit matrix and 0 and r are, respectively, the permittivity of vacuum and the relative permittivity of the mineral phase in the rock.By these relations, the piezoelectric coupling term approaches zero as quartz (being the only piezoelectric phase) phase undergoes damage while, by ( 19), the isotropic permittivity of all three minerals approaches to that of vacuum (Eqs.( 18) and ( 19) are mineral phase specific).It is emphasized that all mineral phases in the numerical rock description have permittivity, i.e. nonzero matrix while only quartz is piezoelectric, having nonzero matrix d.The physical meaning of these matrices can be readily understood from the constitutive equation of linear elastic piezoelectric material: D i = e ikl kl − ij E j , where D i is the electric displacement (C/m 2 ), kl is the mechanical strain tensor, e ikl is the piezoelectric coupling tensor (C/m 2 ), E k is the electric field (V/m), and ij is the dielectric constants tensor (F/m = C/Vm).
The solution of this problem with an electric potential boundary condition (the HF-HV-AC actuation) is here performed explicitly in time using a (globally) non-iterative staggered (partitioned) solution method (Ramegowda et al.

2019; Saksala 2021):
where the Euler scheme is used here for time stepping.An implicit, globally iterative method would not be practical, as the high frequency of the AC voltage to be used in the simulations necessitates an extremely short time step.

Numerical Examples
The simulations testing the weakening effect of the HF-HV-AC pretreatment are carried out in this section.First, the material and model properties are given.Second, the weakening effect on the compressive and tensile strength of a cylindrical specimen is tested.Third, weakening effect in a compression of a rock ball (sphere) is investigated.Moreover, a verification simulation is presented in Appendix 2, where a cantilever beam made of quartz is exposed to a tip point load.Comparison to the analytical solution by Gao and Li (2009) shows that the present model correctly predicts the surface charge density due to polarization via the direct piezoelectric effect.

Material Properties and Model Parameters
The material properties, representing a generic granite like rock, and the model parameters are given in  (Hokka et al. 2016) for Kuru granite.The viscosity moduli values are chosen so as, on one hand, not to cause significant strain rate hardening and, on the other hand, to stabilize the local stress integration.For more details, see Saksala (2021) and Saksala et al. (2017).
The piezoelectric constants and the form of matrix d are given in Fig. 1.These are values in the local crystal coordinate system shown in Fig. 1, which in the present finite element context is the local element coordinate system.Therefore, when other orientations are tested, matrix d is rotated to the desired orientation using e.g.Euler angles (see Saksala (2021) for more details).In this case, the rotation matrix = 3 (ψ) 2 (θ) 1 (φ) where the three rotation matrices i performs the Euler rotation indicated by their arguments as follows: The first rotation is about the z-axis; the second rotation is about the new x′-axis; and the third rotation is about the new z″-axis.

HF-HV-AC Treatment
Numerical piezoelectric treatments using sinusoidal AC excitation of the rock sample are carried out here.Figure 3 shows the boundary conditions for the HF-HV-AC pretreatment, the finite element mesh, and the numerical rock mineral textures.The medium surrounding the rock is ignored, i.e. perfect electric insulation is assumed at the boundaries of the model not covered with the electrodes.
The mineral textures are generated by randomly clustering finite elements with an assigned mineral number (coding in Fig. 3c).This can be achieved as follows: An array with a length of the number of finite elements in the mesh is first filled with integers from 1 to 3 coding the three mineral phases.Moreover, the number of each integer corresponds to the percentage of the constituent mineral in rock.Then, random permutation is performed to this array.Now, as each entry in the array implicitly codes the global element number in the mesh, a spatially mesoscopic description of heterogeneity results, when each mineral phase is assigned with the material properties in Table 1.
The element size in the mesh is ~ 1 mm, which should be small enough to capture the high frequency excitations here (see Saksala 2021).First, the local orientation of quartz crystals is set as in Fig. 3a, i.e. the c-axis and a-axis are, respectively, parallel to z-axis and x-axis (indicated by c↑↑z).Moreover, the loading amplitude is ϕ 0 = 50 kV, and the frequency is set to f = 274.4kHz, which is a natural frequency of the sample NumRock1 in Fig. 3c.This frequency was found altering, by trial and error, the resonance frequency 305 kHz of a similar sample with the same mineral contents, albeit with anisotropic elasticity, used by Saksala (2021).
The simulation results for 500 cycles of the excitation are shown in Fig. 4. Figure 4a and b shows the normalized potential field and the magnitude of the electric field when ϕ = ϕ 0 , at the crest of the sinusoidal wave.The electric potential and the electric field vary between the two extremes cases (the crest and the through of a wave), the first of which is shown here.The second, with ϕ = − ϕ 0 , is similar, albeit with negated colors since the voltage is negative at the through.
Figure 4c and d shows snapshots of the first principal stress distribution and the deformed mesh at t = 3.167E-4 s demonstrating one of the natural modes of the numerical specimen.These high principal stresses inflict damage, both tensile and compressive, on the specimen, as can be observed in Fig. 4e and f.The damage areas show some regularity reflecting the natural mode of deformation (Fig. 4d) related to this natural frequency.The maximum values of the damage variables are ω t,max = 0.76 and ω c,max = 0.64.However, the damaged areas are restricted close to the specimen surface, as demonstrated in Fig. 4g, which shows the elements on the mesh with ω t > 0.1.It should be mentioned that continuing simulation up to 5000 cycles did not result in more damage.Therefore, only 500 cycles of excitation are carried out henceforth.
Next, a simulation is carried out with the same setup except that the local quartz crystal coordinated is aligned so that c-axis is parallel to y-axis while a-axis is as before (denoted c↑↑y).In addition, a higher amplitude of 150 kV is also tested with this orientation.The simulation results for damage variable fields are shown in Fig. 5.With this orientation, the magnitude of damage is much smaller at ϕ 0 = 50 kV (note the range of the colorbar), but the damage pattern orientation displays horizontal bands (Fig. 5a).Moreover, when the voltage is increased to 150 kV, the damage levels increase drastically resulting in what can be interpreted as fragments separating from the specimen.
Next, more natural orientations of quartz crystals (finite elements representing quartz in the mesh actually) are tested.The cases are as follows: (Case 1) random rotation with the Euler angle θ; (Case 2) random rotation with the Euler angles θ and φ; and (Case 3) random rotation with all Euler angles ψ, θ and φ.In each case, the rotation angles vary With these random, and to a large extent more realistic, quartz crystal orientations, the damage levels are quite modest even with such a high amplitude of voltage as 150 kV.Indeed, the magnitude of tensile damage variable barely reaches 0.05 in Case 3. The compressive damage was negligible in this case so that it is not shown in Fig. 6.Interestingly however, the damage patterns exhibit characters similar to those when the crystal axes have clear orientation.This reflects the fact that, despite the different orientations of the quartz crystal axes, which have an effect only on the piezoelectric constants in the present isotropic material assumption, it is the natural mode of the numerical sample that is excited in each simulation here (since the frequency was the same).

Uniaxial Compression and Tension Tests
The standard uniaxial compression and tension tests are now performed first on the intact numerical rock in order to demonstrate the model predictions and to provide the reference cases.Then, the same tests are carried out on the HF-HV-AC pretreated samples to test the weakening effect.A constant velocity boundary condition (BC) is applied at the top of the finite element mesh while the bottom is simply supported.For compression and tension, a velocity of 0.05 m/s downwards and 0.005 m/s upwards are, respectively, applied.The simulation results are shown in Fig. 7.The average stress curve (Fig. 7d) is calculated by dividing the reaction force by the cross-sectional area, while the average strain curve (Fig. 7e) is obtained after dividing the displacement at the top edge by the height of the sample.
The failure modes realized in the uniaxial compression simulations are, for both specimens, the experimentally attested multiple (shear) fracturing mode (Fig. 7a, b) classified by Basu et al. (2013).The experimental transverse splitting mode with a single failure plane is predicted in tension (Fig. 7c).The corresponding compressive and tensile strengths, as attested in the average stress-strain curves in Fig. 7d and e, respectively, are 151 MPa and 8.4 MPa.However, while the details of the failure modes are different, the stress-strain curves are virtually identical in compression for the numerical rock samples in Fig. 3c.The situation did not change upon testing with additional mineral textures.The same finding applies for tension.Thereby, only NumRock1 is considered in the compression and tension test simulations on pretreated rock.
Next, the tension and compression tests are repeated on HF-HV-AC pretreated samples.The simulation results in Fig. 8 for compression show that the weakening effect is quite mild except in the case c↑↑y of the quartz crystal orientation, when ϕ 0 = 150 kV, where the compressive strength (68 MPa) is 55% smaller than that of the intact rock.However, in this case, the numerical sample was already broken due to the pretreatment.In Case 1 of the random crystal orientations, the compressive strength is 136 MPa, i.e. 10% smaller than that of the intact rock.The failure modes for the pretreated samples differ both from each other and from the intact rock, which is a realistic result due to the different initial condition of the rock material in each case.
The results for tension test simulations in Fig. 9 reveal that the weakening effect of the pretreatments is weaker in tension that in compression.Even in the case with c↑↑y and ϕ 0 = 150 kV, i.e. when sample failed already during the pretreatment, the tensile strength is 5.1 MPa, which is 39% smaller than the tensile strength of the intact rock.In other cases, and in the cases with the results not shown, the weakening effect was insignificant.

Compression of Rock Sphere
This test was originally developed for understanding the dynamic crushing of a rock ball (sphere) using the Hopkinson bar apparatus (Huang et al. 2014;Huang 2016).As the application in mind in the present paper is comminution by crushing, a lower impact velocity is applied here.The setup modelled here consists of two rigid platens between which the rock sample is compressed by moving the upper platen downwards with a constant velocity.The modelling principle is illustrated in Fig. 10.
The platen geometry is a virtual and its contact with the rock finite element mesh is modeled by imposing kinematic contact constraints between the rock nodes and the nodes representing the platens.These constraints are of form: u i,z − u RBn,z = b n , where u i,z and u RBn,z are, respectively, the platen node (i = 1,2) and rock ball contact node (n = 1 to number of contacts) displacements in z-direction, and b n is the initial distance them.Moreover, contact forces P 1 and P 2 are solved along with the global solution with a method described in Saksala (2010).Dashpots, i.e. viscous dampers, are required by the method to absorb possible stress waves generated when the velocity v 0 is significant.The equations of motion of the platen nodes are where m 1 and m 2 are computational, not real, masses of the platens and c b b A b is the impedance of the platen defined by the bar velocity of a stress wave in slender rod, density, and the computational area of the platen of the Hopkinson bar device.In the following simulations, values of steel are used for density (7800 kg/m 3 ) and the bar velocity c b = 5150 m/s, while m 1 = m 2 = 0.1 kg.Other material parameters and model parameters are as above.The interesting outcome of the simulations here are the contact (reaction) forces since they are the loading carried out by the consumable tools in comminution.Figure 11 presents simulation results for crushing the numerical rock balls in Fig. 10c when v 0 = 0.5 m/s.
The failure modes for both numerical rock balls exhibit the experimentally observed (Huang 2016) split of the sample in three or four slices.Moreover, some secondary (meaning mild value of the damage variable) crack-like tensile damage patterns can be inferred from both samples.Substantial compressive damaging can be observed (Fig. 11a) at the poles of the sphere (i.e. in the contact areas).The contact reaction forces, reaching a magnitude of 18 kN, show some fluctuations (Fig. 11d) due to the Finally, Fig. 11e shows the average of the contact forces for both rock samples, and as can be observed, they are almost identical.

HF-HV-AC Treatment
The HF-HV-AC pretreatment was carried out on the numerical rock ball (NumRock1) in Fig. 10c.The electrode configuration is shown in Fig. 12 along with the simulation results for cases c↑↑z and c↑↑y.Again, 500 cycles were simulated with the loading amplitude of ϕ 0 = 100 kV, and the frequency was set to f = 279 kHz.As can be observed in Fig. 12c and d, higher values (note the different ranges of colorbars) of damage variables are, again, generated when the c-axis of crystals are aligned with global z-axis.In this case, there is an interesting crack-like tensile damage formation on the equator of the rock sphere.When the c-axis is aligned with the global y-axis, the damaging is strongest at the pole area of the sphere where the electric field strength (Fig. 12b) has its maximum due to the edge effect in the voltage BC, i.e. the ground electrode.
Next, the same cases with random rotation of the quartz crystal c-axis as above are tested here as well.Moreover, one additional simulation for Case 2 is carried out with the electrodes aligned, instead of z-axis, with y-axis (designated as Case 4).The results are shown in Fig. 13.
In each case shown in Fig. 13, except Case 3, the tensile damage variable reaches almost 1 at some locations.However, Case 3 (Fig. 13c) is the most unfavourable to damaging due to the total randomness of the quartz crystal orientations.

Compression of the HF-HV-AC Pretreated Rock Balls
The compression test detailed above is now performed on the pretreated rock samples.The results are shown in Fig. 14 for failure modes, as represented by the tensile damage field, and for the corresponding average contact forces.The viewpoint, with respect to the global coordinates, is the same for each tensile damage plot as the one shown in Fig. 14a.
The failure modes for the HF-HV-AC pretreated samples display axial splitting of the sample into two or three parts.The average contact forces for the pretreated cases vary from 16.4 to 17.6 kN while the intact samples withstand 17.7 kN of average compressive force.Therefore, the weakening effect in these cases is quite low.

Effect of Quartz Content
The simulations above were carried out with the same mineral composition of the numerical rock.It is instructive to test the effect of quartz content, as it varies from granite to granite.Moreover, it was already demonstrated by Saksala (2021) that the quartz content influences the frequency at which the maximum damage occurs.For this reason, a set of simulations are performed here with the mineral composition 45% of quartz, 47% of feldspars and 7% of biotite.Changing the mineral composition, changes also the mineral texture, denoted as NumRock3.In this case, the frequency at which (local) maximum damage was caused is 285 kHz.Changing the quartz content has an insignificant effect, at least within the present modelling approach, as can be observed in Fig. 15e.More specifically, the maximum contact force is 17.8 kN for the intact rock while the strongest weakening occurs in Case1with 16.4 kN of maximum contact force.These figures are almost the same as the ones above when the quartz content was 33%.

Discussion
First, the weakening predicted in uniaxial compression and tension tests are addressed.A 10% weakening is predicted for compressive strength in Case 1 with 150 kV of voltage.
It should be noted that this case is a more realistic one (compared to the case with c↑↑z, for example), as piezoelectric effect in rocks may arise when the optical axes of quartz grains (c-axis) are randomly distributed in a common plane, while the electrical axis (a-axis) coincides with the normal of this plane (Parkhomenko 1971).The highest weakening in these simulations was predicted in the case with c↑↑y at 150 kV of voltage resulting in a partial fragmentation of the sample.However, this case, as having perfect alignment of quartz crystals, is not expected to be found in natural rocks.
Second, the predicted weakening effects in the more practically oriented problem of crushing the rock ball by diametral compression are quantified in Table 2 to facilitate discussion.Despite the high voltage, the predicted weakening effects are somewhat mild, i.e. barely 8% (NR3) in the most optimistic case tested.In addition to Case 1, notable weakening (6.2% in Case 2 and 7.3% in Case 4 for NR2) was predicted also in Case 2, which further randomize the local quartz crystal axes.However, if there is no texture or preferred orientation of the quartz crystal axes at all (Case 3), the weakening effect is insignificant, ~ 1%.
As to the mineral composition, it has an effect, as was already shown in Saksala (2021).More precisely, changing the quartz content, the resonance frequency changes.Indeed, as demonstrated in the numerical simulations above, when the quartz content was increased, the resonance frequency increased from 279 to 285 kHz.Moreover, it was found that lowering the quartz content to 20%, decreased the resonance frequency to 272 kHz.However, the results for this case were not presented because the weakening effect was similar to those in other cases of quartz content.
The most notable weakness of the present modelling approach should also be addressed.Namely, the present continuum approach based on a damage-viscoplasticity model and standard finite elements describes the fractures and cracking in the smeared sense, which may not always be realistic for brittle materials such as granite rock.This drawback can be highlighted by comparing the study by Saksala et al (2022) to the present one with respect to uniaxial tension.Saksala et al (2022) predicted 12% weakening effect in tension for a similar cylindrical sample as here, while only 2.4% weakening is predicted above (Fig. 9) in the same case of the quartz crystal orientation and at the same voltage of 50 kV.However, Saksala et al. modelled the rock as anisotropic, and described the cracking rock with the embedded discontinuity finite element method, which is superior to the present approach in fracture description.

Concluding Remarks
This paper presented a numerical study to quantify the weakening effects by piezoelectric actuation of quartz mineral phase in granite rock using high-voltage and high-frequency alternating current.The developed numerical method consists of a damage-viscoplasticity model and a globally noniterative staggered scheme to solve the coupled piezoelectro-mechanical problem explicitly in time.This method was applied in standard uniaxial tension and compression test simulations as well as in simulation of diametral compression of rock balls.Considering these applications, following closing remarks are delivered.
In uniaxial tension and compression tests: Actuating a cylindrical numerical granite specimen at one of its natural frequencies, 274.4 kHz in the present case, with the voltage of 50-150 kN, considerable damage can be induced.However, the weakening effect of the pretreatment-induced damage fields on the tensile and compressive strengths of the sample were quite modest.Moreover, the amount of damage inflicted, and hence, the weakening effect depends strongly on the orientation of the local quartz crystal coordinates so that more randomness in orientation results in less damage and weakening.In a realistic case, where the optical axes of quartz grains (c-axis) are randomly distributed in a common plane, while the electrical axis (a-axis) coincides with the normal of this plane, a weakening effect of 10% in compressive strength was predicted after treating the sample 500 cycles with 150 kV of amplitude.
In diametral compression of a rock ball: In this important problem with respect to comminution practice, the predicted weakening effect due to the piezoelectric pretreatment was similar to that in the uniaxial compression.Namely, 7.9% weakening in the average contact (reaction) forces was predicted in a realistic case of the quartz crystal orientations.Moreover, the weakening effect was insensitive to different quartz contents of the numerical rock.Altering quartz content also changes the natural frequency of the specimen at which the (local) maximum damaging occurs so that the higher the quartz content, the higher the frequency.
It should finally be mentioned that the present continuum damage-viscoplasticity model describes cracking in brittle materials in a smeared sense ignoring, to a large extent, the crack tip effects.Therefore, the present simulations provide a lower limit or conservative estimates for the weakening effect of the HF-HV-AC pretreatment.Advanced rock fracture models, taking more properly into account the crack tip geometry effects, could probably make these predictions more optimistic.In any case, experiments are sorely needed to ultimately decide whether this kind of pretreatment is a feasible one for comminution.

Appendix 1
The return mapping scheme related to the corner (visco) plasticity situation is elaborated in detail here.The generalized cutting plane algorithm is chosen (Simo and Hughes 1998).It is based on the standard elastic predictor-plastic corrector split.That is, trial state is first calculated by: where t+Δt = e t+Δt is the new total strain and ̄ pr trial is the trial principal stress.With the rate independent plasticity models, the consistency condition, ̇f = 0 , is usually exploited for solving the viscoplastic multiplier λ .However, with the viscoplastic consistency approach, this is not feasible, since the consistency condition has the following form in the general case (Wang et al. 1997) The consistency condition thus becomes a first order differential equation for λ .However, an alternative, algo- rithmic avenue opens upon considering that, at the end of the time step, condition f i t+Δt , λt+Δt,i = 0 (i = MH, MR) must be satisfied.Now, the rates of the internal variables can be replaced by λt+Δt,i = Δ t+Δt,i ∕Δt (Winnicki et al. 2001), yielding the vector form of the criteria as ( , Δ ) = f MH , Δ MH f MR , Δ MR T .Then, expanding this with the first term of the vector valued Taylor series gives: with where signifies gradient with respect to .Now, it is assumed that a potential corner (visco)plasticity case is realized, i.e. the trial state violates both of the criteria.The return mapping proceeds now as follows.
It is implicitly understood that the steps above are performed in the principal stress space so that, as the final step not shown, the corrected stress and strain are transformed to the xyz-coordinates.In addition, the algorithm outlined includes the procedure for the case in which a genuine corner plasticity situation does not realize, i.e. the return is not to the intersection of the yield surfaces but to either f MH or f MR .A need to return to a single surface is signalled by a negative value of the cumulative viscoplastic increment, i.e. either Δ (i+1)  MH < 0 or Δ (i+1)  MR < 0 .When such a situation occurs, the algorithm is restarted with a single surface (Step 3).This trial-and-error procedure is required since, according to Simo and Hughes (1998), the set of active criteria cannot be known in advance.This is because the positivity of failure criteria in the trial step does not guarantee that a genuine corner (visco)plasticity case has been realized.

Appendix 2
A problem of a cantilever beam made of quartz under a point load is solved here in order to verify the piezoelectro-mechanical part of the numerical method presented in Sect. 2. The beam with the finite element mesh consisting of 6930 tetrahedrons is shown in Fig. 16a.In the simulation, the point load is applied slowly by increasing it to a unity during 1000 time steps.The analytical solution of the problem is by Gao and Li (2009), which for the surface charge density reads: where P is the magnitude of the external force, while d 31 and d 36 are the piezoelectric coefficients (see Fig. 1).
Figure 16b-c shows the simulation results for the bending stress, potential (voltage), and the surface charge density.The distribution of voltage in Fig. 16c is similar to that reported by Gao and Li (2009).Unfortunately, they don't show the values of the potential so that quantitative comparison is impossible.However, the surface charge density in Fig. 16d is in a fair agreement with the theory (Eq.( 31)).

Fig. 2 a
Fig. 2 a Yield criteria in the principal stress space (σ t = 10 MPa, σ c = 165 MPa); b a detail illustrating the corner plasticity case at the intersection of the criteria

Fig. 7
Fig. 7 Simulation results for uniaxial compression and tension test: a tensile and compressive damage fields for NumRock1 in compression; b tensile and compressive damage fields for NumRock2 in com-

Fig. 8 Fig. 9
Fig. 8 Simulation results for uniaxial compression test on HF-HV-AC pretreated samples (NumRock1): a tensile and compressive damage fields for c↑↑z with ϕ 0 = 50 kV; b tensile and compressive damage fields for c↑↑y with ϕ 0 = 50 kV; c tensile and compressive damage

Fig. 14
Fig. 14 Simulation results for compression of the HF-HV-AC pretreated rock balls (f = 279 kHz): a tensile damage field for NumRock1 with c↑↑z and ϕ 0 = 100 kV; b tensile damage field for NumRock1 with c↑↑y and ϕ 0 = 100 kV; c tensile damage field for NumRock2 with Case 1 and ϕ 0 = 150 kV; d tensile damage field for NumRock2

Table 1 .
The elasticity constants (E, ν) are generated with the MTEX software (Mainprice et al. 2014), and the rest of the material properties are by Hokka et al. 2016, Newnham 2005 and Mahabadi 2012.The dielectric constants in Table 1 are given relative to the permittivity of vacuum 0 .Moreover, the MH criterion parameter values, B and n, match the quasistatic experiments

Table 2
Predicted contact forces in diametral compression of rock ball If f (i+1) MH < TOL & f (i+1)MR < TOL then go to Update.Else set i ← i + 1 and go to 1.