Effect of Dynamic Stress Produced by Rock Blasting on the Optimal Dimensioning of Room and Pillars in Horizontal Layers

Abstract

An optimal dimensioning of the room and pillars of the square or rectangular section is finding adequate transversal section support for the pillar; in that sense, all forces acting on the pillar must be evaluated. The pillar has its ends subjected to normal (σ) and shear stresses (τ). In addition, the pillar can also be subjected to dynamic stresses from the seismic waves generated by rock blasting. Thus, when calculating the dimensions of a mine pillar, the stresses imposed by seismic waves must be considered. The combination of lithostatic and dynamic stresses generated by blasting can affect the pillar, generating regions, or areas that will exhibit the formation of cracks, microcracks, chipping, and other detrimental processes for its stability. In that regard, it is appropriate to assess the dynamic effects of the energy delivered in each wave pulse so that these dynamic quantities can be included in the dimensioning of the pillars. This paper presents the physical and mathematical formulation describing the dimensioning issues of room and pillars and its solution using numerical modeling. The numerical modeling is performed using nonlinear mathematical programming and simulation with FLAC 3D (Fast Lagrangian Analysis of Continua in 3D). In this study, it was proven that the stress generated by the seismic wave represents a dynamic stress of nearly 5.5 MPa. 'Dimensioned pillars considering wave stresses allows it to recover an average 1.6% less when it is compared with pillars only under overlying rock mass'.

1. Introduction

Mining is an activity whose purpose is to extract minerals that society needs and uses economically. Mining requires that the financial resources be effectively allocated [1]. Mining operations include the activities of drilling, rock blasting, loading, and transportation [2–4]. Rock excavation is the first activity in mineral production. It can be achieved using mechanical methods or with the use of explosives. Traditionally, at industrial-scale production, explosives are used for rock excavation. The explosives must be distributed in holes systematically in the rock mass to be blasted. This task is made to fragment the rock by releasing the explosive’s energy. It follows the principles of physic dynamics, cinematic, and the expansion of hypercompressed gases. Not all the energy released by the explosives is used in the fragmentation process, resulting in collateral effects such as fly rock, airblast, and seismic waves.

As a product of the detonations, the seismic waves will result in the loading and unloading of the rock mass in a repetitive manner, similar to the conditions of a body subjected to laboratory tests of simple cyclical compression [5, 6]. The repetitive application of stresses can be detrimental for the geotechnical parameters of the materials (Young’s modulus (E), shear modulus (G), constant K, and yield stress (σ_e)). Additionally, the magnitude of the tensile strength (σ_t), compression (σ_c) of the rock mass, and the factor of safety (FS) of mine structures are constantly subjected to variations as a consequence of the above-mentioned dynamic stresses [7].

Regarding a mine that uses the room and pillar mining method, the mine pillar is a structure subjected to lithostatic stresses (primary stresses), corresponding to the weight of a prismatic volume of an overlying rock that acts on its top [8]. The conventional sizing of a mine pillar does not include the effect of the wave stress produced by rock blasting on these structures. For this reason, the methodologies that allow integrating these dynamic loads in the sizing of mine pillars must be adopted.

The purpose of finding the pillar’s optimal dimensions is to obtain the safest geometry for the room and pillars system with the maximum recovery of the minerals in the deposit; hence, a detailed evaluation of the dynamic stresses acting on the pillar must be performed; accordingly, when dynamic stress is simplified under the premise “its value is negligible,” such assumption characterizes an error in the dimensioning of a pillar, since pillar stays overloaded (underdimensioned because of simplified estimates to find its dimensions). Thereby, there is a risk that pillar is working in transient regime close to its yield stregth limit, and ‘any increasing’ stress on the pillar would drive its strength parameters to break. Therefore, it is necessary to carry out rock mass monitoring with geophones/seismograph for measuring dynamic problems whose aim is controlling the energy released by rock blasting.

2. Methodology

This paper uses a novel methodology to include the stresses generated by blasting around the pillars for its dimensioning. The methodology can be summarized as follows:(1)Analysis and evaluation of recorded seismograph information, lithological information, and geotechnical data of the area under study are achieved(2)Characterization of the dynamic properties of the rock mass using the Fourier Transform to understand the mechanism of energy dissipation in the rock mass domain includes the application of conservation theories of energy, mass, as well as determination of the particle velocity responsible for displacements and deformation, and contrasting this information via mathematical models and simulation with a commercial software (I-Blast V.8)(3)Formulation and development of equations integrate the sizing mine pillar model(4)Design one algorithm that provides the solution for dimensioning room and pillars through nonlinear mathematical programming and simulation with FLAC 3D(5)Finally, there is analysis of results

2.1. Seismic Waves Generated by Rock Blasting

The study of oscillations and vibrations is a fundamental part of dynamic studies in physics [9]. In general, the vibrational wave is attenuated and can be characterized as a damped vibration where a portion of the energy delivered by the wave is ‘absorbed’ by the rock material [6]. The deformation produced in a material is proportional to the applied tension and can be represented by the postulates of Hooke’s law (σ = Eε, with σ being the stress, and E and ε are Young’s modulus and unit strain, respectively).

The detonation of an explosive charge in the rock mass in a blast hole imposes a dynamic loading in the vicinity of the structures, and, depending on the magnitude of these dynamic loads, the rock mass will begin to fail [10].

2.2. Conservation of Movement, Conservation of Mass, and Energy

The work developed by force ‘F’ is equivalent to the variation of energy (dW = −dU); such a statement is supported by the reduced equation of Bernoulli:where V, P, ρ, and K are the velocity, pressure, density, and constant, respectively.

The gases responsible for the fracture are dissipated by exerting an internal pressure in the middle; from the Bernoulli equation, it is verified that the pressure at any point in the rock mass can be found once the pressure (inside the hole) and the velocity at an initial point are known.

Tesarik and Hustrulid [11], used hydrodynamics concepts to assess the peak particle velocity (theoretical velocity) at a distance from the origin due to an explosive charge. This methodology is based on the continuity equation (, where is the density, is the time, and is the velocity vector.

The explosive energy fragments the rock in the process, resulting in a material with the behavior of a flow (by an area element ). The flow of the material can be described as work with a domain (D), and the flow problem can be analyzed using the Stokes theorem (the first following equation) and Gren’s Potential theorem (the second following equation). By this way, the equilibrium relationships for the system in motion are expressed by the movement conservation equations (the third following equation) [3, 12–14].where the variable represents the stress acting on the elementary unit, is the volumetric mass (density), are the velocity vectors, is the gravitational component, is the gradient operator, and x, y, and z are the planar coordinates.

The conservation of mass is applied to fluids (flows) in motion, and its generalized equation is expressed through the integral surface of a vector function (, is an elementary area). This analysis will consider some basic criteria accepted by Tesarik and Hustrulid (2009): irrotational flow and rock incompressibility ( or ) have a behavior similar to an inviscid liquid.

Particles that carry the seismic wave are subjected to stages of excitement, and the wave passage is responsible for the potential energy of deformation and the kinetic energy [9].

The quantity of energy released in each wave pulse is related proportionally to the amplitude. Its magnitude depends on the parameters selected in the drilling operation, rock mass properties, and blasting sequencing [15]. PANEIRO states that the liberated energy is suddenly in a short interval of time in the form of high pressures and temperatures, which are sufficiently big enough to produce the fragmentation of the material, and the residual stresses dissipated in the vicinity produce the vibrations in the ground [16]. The internal energy (U) is generated from the sudden disintegration of the explosive, characterized as an “adiabatic process,” which is transformed into energy work and is manifested by rock fragmentation and subproducts.

The blasted rocks under this kinetic process follow the principles of the Newtonian flow. Considering the criterion of continuity and conservation of mass and energy, the equations derived from Navier-Stokes allow finding an implicit way through which the energy dissipates by the wave in the rock volume. This law describes the movement of the flow in the middle, and, for that, the following premises for the flow in question must be assumed: density and viscosity (μ) are constant.

2.3. Applied Methodology for Determining Dynamic Stress

2.3.1. Vibrations and Stress on the Ground

After detonating the explosive material in the blast hole, a set of stress waves are generated, traveling through the medium as vibrations. Many researches proposed equations (Duvall and Fongelson, 1962, Ambraseys and Hendrom, 1968, and Indian Standard Predictor, 1973) to estimate those particle velocities. The velocity can be measured using seismographs [17], and the signal information obtained from the seismograph records provides relevant data about the behavior of the ground.

The ‘maximum standby charge method or scaled distance’ is amply accepted to predict a vibration level at a certain point. One of the methodologies to determine the PPV, as well as the characteristics of the ground, is expressed by Devine’s equation [3].

With PPV being the peak particle velocity, D is the distance between the location of the explosive charge and the point of interest, Q is the amount of standby charge (explosive charge), and α and K are the attenuation and the ground coefficients, respectively. In addition, for determining the PPV in a short field (close to the explosive charge), the Holmberg and Persson equations can be used, and the PPV may also be determined by the hydrodynamic theory as mentioned in Section 2.2; this criterion has a physical-mathematical basis in fluid dynamics [18].

Zhang admits to an empirical relation to find the deformations made by a wave from a distance ‘r’ in the detonation point (nf(t) to f(t) = B or f(t) = B, where K is the rock mass constant, B takes the value of 3 (however, this value has a mathematical basis and can be found from the rock mass characterization and the function decoded via FFT), the constant r_h is the radius of the hole, and a constant β is the attenuation coefficient [19].

2.3.2. Explosive Energy

The shock of explosive energy relates to the detonation velocity of the explosive through the following expression:

With P_d, being the detonation pressure (GPa), is explosive density; VOD is explosive detonation velocity (m/s); and γ indicates the specific heat ratio of the gases produced by the detonation [20]. On the other hand, the energy delivered by the explosive to fragment the rock is a function of the pressure mentioned above; the explosive density is corrected in the Chapman–Jouguet plane () and Gurney’s factor () for explosive mixtures with ANFO (Ammonium Nitrate-Fuel Oil) and for mixtures with emulsion [19].

2.3.3. Frequency Analysis of Wave Signals

The frequency analysis of the seismic waves can be done using the Fourier Transform. The Fourier Transform allows converting the time-domain waveform in the frequency-domain [21]. This is done using where x(m) is the magnitude data (frequency information); x(n) is the waveform signal in a discrete form; , m, n, t, and N are the angular frequency; and m is the component associated with each sample spectrum, n is the component associated with the parcels that characterize the series of coefficients, t is the time, and N is the amount of data.

2.3.4. Assessment of the Stresses Generated by Energy Dissipation through Waves on Support Structures: Mine Pillar

A general equation for the propagation of a plane wave (on the x axis) in a rock medium, continuous and isotropic:

This equation is a simplification of a 3D wave, whose expression is defined by

Deformations are proportional to the wave velocity in the rock mass, which is with ‘’ being the wave velocity at point of the analysis; additionally, the seismic energy can be expressed as a function of the velocity. Therefore, it is possible to write an equation that represents the behavior of the disturbance based on the following mathematical model: , where U_o is the fundamental amplitude, k is the constant attenuation, c is the peak particle velocity, and x is the distance [11].

In the same order of ideas, the paper in [9] establishes a relation of energy associated with the deformation of the material and the acting stress . This differential expression is associated with variable amplitude; the author considers a solution based on a complex conjugate model to represent the decay of the wave with time and formulates the following function: and its implicit solution:  = 0; the constants A and ∝ allow adjusting the solution to the values of displacement and velocity, while p is a value associated with variable angular frequency, γ is a constant related to the viscosity of the material, and is the natural frequency.

According to the expressions proposed by the authors mentioned above, it is evident that the variable amplitude can be written as follows: These different relations and the formulation of Devine can be equated to define the magnitude of the stress of the acting wave that produces the deformation.

On the other hand, Goodman sustains that knowing Young’s modulus values and Poisson relation (μ), the deformations are defined for a state of three-dimensional stresses in an elastic solid: isotropic and linear [5]. From the concepts of stress-deformation and the theories by the above-mentioned researchers, it is possible to infer the following equation, which describes the state of stress produced by the passage of the wave in the rock mass. Although the design of the following equation is to consider Hooke’s law, it has complex solution elements (attenuation coefficient); its solucion depends on the elastic parameters of the rock mass, expression representing the displacemente component , time, and the position of the pulse .

Note that, in a more restricted concept in (9), the variable x seems to be the position of the pulse responsible for the deformation (different from the general concept of ‘wave path’), and, in this understanding, x is replaced by the affine x’, not to create ambiguity in the understanding; thus, x’ becomes an ‘independent’ and differentiable variable. At the same time, it is a function linked to the dimensions of the spans and the pillars. Analogously, the variable t can also be conveniently rearranged as an expression of velocity (speed of the wave in the process of attenuation and the dimensions of the pillar and spans).

2.4. Stresses Acting on the Pillar

2.4.1. Dynamic Stress on the Mine Pillar

The function defines the behavior of the ‘displacements’ of the particles when subjected to dynamic stresses, and the variables x, x´, t, and express the pulse position, time, and speed.

There is a function that expresses the deformations associated with the variable tension ; such relations are linked implicitly to the velocity of the wave dissipated in the rock mass. The generalized expression of Hooke’s law allows us to relate the ‘stresses and strains’ since the size must meet the stability criteria:

On the other hand, (10) can be written as a function of the characteristics of the traveling wave propagating in the rock mass. The component represents the dynamic stress, and its equation is a complex derivate expression, whose solution is composed of a real and an imaginary part

In an attempt to facilitate the presentation, the math presentation of (11) is summarized throughout symbolic operators (a, A, B, C, D, E, and F).

The geometric variables (W_p, L_p, W_o, and L_o) are conjugated of the tributary area dimensions; E_m, is the rock mass elasticity modulus; is rock mass density; is density of explosive; E_U is the useful energy provided after the detonation of the explosive, is a factor associated with the amount of charge (ANFO) per hole, r is the radius of the hole, and is called Poisson’s ratio.

Figure 1 is a schematization of the process of mineral extraction by the room and pillar method, in which the origin of the coordinates (O) is located in the center of the room section (W_oH_p), x, u(x = 0, t)). The dynamic tension (σ_o) active in the pillar (Figure 1) is defined by the quantity of charge delivered in each pulse of energy when the wave runs through the distance ‘x-vt’ in Figure 2 on the horizontal axis OX. Therefore, the value of σ_o, must be added in the lithostatic stress acting on the pillar (σ_p).

Figure 1 

Wave stress in the mineral extraction process with explosives in the room and pillars method.

Figure 2 

Dynamic stress dissipated in the pillar after the hole detonation started (adapted from [19]).

3. Case Study

The sizing problem of the pillars presented in this document refers to a real situation that happened in a mining company situated in Brazil. The company performs mechanized mining of mineral resources. In this study, an emphasis has been put on the room and pillar sizing, considering the stresses delivered by the dissipated waves of explosive detonations; after these produced energies are quantified, they will soon be incorporated in the formulation that integrates the algorithm to find optimal dimensions of the section of pillar support. It is opportune to quote the resistance criterion of Hoek and Brown’s 2002 version, which considers the empirical factor ‘D’; such a factor ‘D’ correlates the damage produced by the excavation of rocks with explosives and the resistance of the rock mass in an empirical way. However, with the support of physical tools—mathematics by what is cited in the methodology—the tools of the wave are studied and included in the sizing process of the room and mine pillars. Thus, the sizing problem is patterned and the results obtained represent information with bigger accuracy in contrast with simplified models or empirical foundations, in which the wave voltages are disregarded.

Regarding the characteristics of the deposit, the ore layer is at a depth between 300 and 350 meters and has a specific weight of 0.027 MN/m³ to a barren rock (overburden), whose mineral occurrence is given in layers of 2 m of thickness in a horizontal position. The roof that comprises the spans shows a sedimentary structure (slab) of 1 m of thickness, and the mineralized rock is a quartzite; in Table 1, the lithological materials as well as the geotechnical data for these materials are broken down. The mining company in question uses the ANFO (total energy of the explosive: 912 kcal/kg) explosive, which has a specific mass of 900 kg/m³. The approximate loading ratio of 3.3 kg/m³ (4.5 kg/hole) for a stope whose average section is 4.3 m × 2 m, with holes of 45 mm in diameter and 45 holes, is distributed in the blasting stope.

Table 1 

Data and geotechnical information.

3.1. Analysis of Ground Characteristics

For this case study, the seismograph data was recorded in a stope, with a standard section of 4 m × 2 m; the sensors were installed at distances between 30 and 75 m from the stope. Table 2 includes the information corresponding to (1) ‘peak velocity (V)’ corresponding to the channel ‘L.’ The L channel represents seismographic data with higher intensity information and, therefore, the wave energy is increased and (2) scaled distance is illustrated in the same table (SD (), with D being the distance between the stope and the sensor = 70 m, Q = 70 kg). The parameters for the rock mass are ‘K = 153.04 and α = −1.463’. The information for K and α is used to analyze the wave velocity in a short field in the software I-Blast V. 8.

Table 2 

Seismograph records and scaled distance information.

3.2. Analysis of Seismographic Data through Fourier Transform

The information in Table 3 includes the result of the analysis of the seismograph records after being analyzed using (7). The table includes (a) phase angle , (b) amplitude (Uo) ≈ 0.45 m, and the fundamental frequency (f) = 48.7 Hz with the dominant angular frequency ( = 306 rad/s.

Table 3 

Data corresponding to the decoded wave spectrum.

3.3. Analysis of the Velocities Responsible for the Stresses in the Pillar

The software I-Blast V.8 simulates the wave propagation, derived from information contained in Section 3.1 and Table 2. In this way, the stress found in the centroid of the pillar is an implicit function of liberated energies, which are expressed in scalar magnitudes of velocity. It must be highlighted that the simulation provides the peak velocity of a particle; therefore, these wave magnitudes represent the maximum value of the wave pulse velocity.

3.4. Lithostatic Stresses and Effective Stresses on the Pillar

These are stresses related to the lithological load exerted by the rock mass on the pillar that works as a support at the applied depth. Table 4 presents the average lithostatic stress, as well as the pillar effective stress (transversal section of pillar support ≈ 10.3 m² and tributary area ≈ 72 m²). Such information corresponds to the historical data provided by the geomechanical department.

Table 4 

Lithostatic stress and pillar-induced stresses.

3.5. Approach to the Problem

The problem of pillar dimensioning, including the dynamic stresses, was solved using numerical modeling through nonlinear mathematical programming (PNL). It is possible to establish those recovery functions of deposits whose variables are the geometric unknowns, which must be respected in order to guarantee the mining arrangement. The combination of the recovery functions and the geometric restrictions must be written in a function of structural interdependence that characterizes the algorithm as a whole, which will allow obtaining the sought answers. In addition, the information found in the mathematical model, in which the dynamic forces are incorporated, will be evaluated by simulation.

3.5.1. Math Formulation

Betts shows a solution set, in which the structure establishes the objective function and the restrictions to be respected [22]; it is as follows: evaluate the n–vector: x^T = [x₁, x₂, …, x_n], to maximize (or minimize) the function F = F(x_i).

Subject to ‘m’ restrictions, C_L ≤ C(x) ≤ C_U.

The limits of coherence: x_L ≤ x ≤ x_U ∀x 𝜖 ℝⁿ/ℝ ⟶ C^r, where x^T is the n-dimensional vector of the project variables; F(x) is the objective function to be optimized; C(x) corresponds to the restrictions that condition the solution; the values of x represent the permissible magnitudes associated with the coherence of the model and decide the space of solutions in the domain Rⁿ, whose values for the variables are of the real type within the closed-convergent bound that encompasses the solution, and the expressions are continuous and differentiable functions (C^r).

3.5.2. Objective Function

The math expression of the objective to be reached must involve the corresponding variables to the geometric data of the mining arrangement (see Figure 3); in this way, the recovery of the function (Rec) is expressed by Rec = Rec (W_p, W_o, L_p, and L_o), with W_p, W_o, L_o, and L_p being the pillar width, span width, pillar length, and span length, respectively.

Figure 3 

Arrangement of mining via room and pillars method: ‘plan view.’

Considering a regular arrangement and a uniform distribution of the pillars in the area corresponding to the mineral block, as well as assuming that the entire thickness (H_p) of the ore layer is mined, the recovery function (objective function) is defined as follows: ; when L_p = W_pe W_o = L_o, the expression is reduced by [23–25].

3.5.3. Tributary Area Theory (TAT)

This method is based on the proposition of Newton’s third law [6, 24–26]. For this, the strength of pillar equals the weight of the volume of overlying prismatic-shaped rock, overlying the vertical reaction exerted over an influential area (At) or “tributary” (such dimensions include the pillar itself and its surroundings going up to half of the adjacent spans) (Figure 3).

The average stress that acts over the pillar, as given by , is compared with the pillar strength increased by a safety factor; the average stress must be less than or equal to the strength of the pillar [27].

In the formulation above, it is assumed that the pillar section (A_p) is subjected to an uniform stress; this assumption takes the following premises: (1) there is flat and extensive topography in relation to depth and (2) there is no redistribution of stresses in the mine confines, and it is assumed that all pillars are under the same axial load [12].

3.5.4. Pillar Strength

The magnitude of pillar strength is bound to the shape effect [28]; in this understanding, the expressions that give the pillar strength must integrate the geometric variables of the pillar and soon be equated based on the TAT, characterizing a limit restriction.

On the other hand, the requested resistance in the pillar must consider the safety margin (FS_pilar ≈ FS) so that the structure does not have to face the ultimate stress regime (yield stress). Given the above, the equation that satisfies this requirement is .

A lot of research proposed mathematical relations deduced from retroanalysis and by laboratory tests. The two expressions that are broadly accepted and used for pillar design projects are the equation of Obert and Duvall [29] as in the first following equation and the equation of Salomon and Munro (1967) as the second following equation (14), well known for its application:

With being the strength of a rock volume of unit dimensions, W_p/H_p is the relation of pillar slenderness (W_p, H_p being pillar width and height, respectively). It is worth highlighting that (11) has been countersigned through algebraic proof by Pariseau [24]; such a demonstration must be sustained according to Mohr–Coulomb’s rupture criterion.

Regarding pillars with a rectangular support cross section, Mark and Chasse (1997) proposed an expression that involves the pillar edge dimensions, shown in the following equation. This expression is a derivative obtained from the Bieniawski equation:

3.5.5. Strength of Foundations

The acting load (lithostatic weight) at the top of the pillar is ‘passed on’ (equilibrium principle) to the floor; consequently, there is a demand for strength from the loaded pillar. It is verified that the floor is charged by the ‘retransmitted charge’ and the charge is due to the mass of the pillar itself. In order to calculate the resistance capacity of the foundation to hold loads due to 'overlying weight', it is advised to use the classic expression applicable to shallow foundations, as proposed by Terzaghi [6]: Q_b is the load capacity of the foundation in units of force/area  is the angle of friction of the foundation rock c is the cohesion of the foundation rock γ is the specific weight of the foundation rock W_p is the width of the pillar Sγ = 1.0 − 0.4(W_P/L_P) and Sq = 1.0 + sin, with L_p being the length of the pillar

Therefore, the equation that characterizes this restriction is written as.

3.5.6. Stability of the Layers That Compose the Roof of the Room Spans

In the estimated width of the spans, two considerations are made [2, 26]: (1) the access widths must guarantee the operability of equipment and production systems as well as the supply of services (e.g., air for ventilation) and (2) the structure and stability of the roof must be the same as the other parts.

In another view, the roof can be formed by a slab or a set of stratified and sedimentary slabs, arranged in layers, and can be determined as a structure flexion through concepts of an elastic beam [30]. Goodman uses the mechanic concept of the materials, specifically the elastic beam theory, to analyze behavior [5]; Pariseau formulates that the slab behavior submitted to pure flexion considering the plane strain theory [24]. For that, the restriction of stability cited by the aforementioned authors is written in the following manner [30], in which the analysis of the ‘n’ slabs willing to flex under an equivalent charge ‘Q’ is as follows:where Q is a load value (stress) distributed by unit of length; E_i and E_j are the modulus of elasticity of the first slab and the upper slabs, respectively; and h_i and h_j are their respective thickness. The Q value is a function of the ‘n’ “stable” slabs:  If Q(1) is bigger than Q(2), it follows that the second layer (superior) is loading at first (bottom) If Q(1) its smaller than Q(2), it follows that the second layer is supporting part of the weight of the first

This analysis must be realized for the entire structure of the ‘n’ layers that make up the ceiling (1, 2, 3, …, n), and, in the absence of tension between the layers Q(n) and Q(n + 1), there is decoupling between the superior and bottom slabs.

The equation for the restriction of width (or length) of the span is defined as follows:

Here, h is the size of the layer; is the tensile strength; Q is the stress exerted by the set of slabs that will inflect and is the security factor for the slabs. By what was already mentioned, the minimum span coherence constraint is written like this: ; L_min is the minimum width, usually defined as a function of the largest equipment that will travel inside the mine, with a gap on each side of the road.

3.5.7. Numerical Modeling and Simulation

The mining area under study is in a domain presenting 3 types of lithology. The information and input data (geometric arrangement) are provided by the optimization through the algorithm (as described in Section 3.1); therefore, the formulations corresponding to the hydrodynamic theory are placed in the model script. The information for the modeling is the following:(1)For square section:(a)Pillar section: 2.66 × 2.66 m² and room of 4.0 × 4.0 m²(b)Pillar section: 3.21 × 3.21 m² and room of 5.24 × 5.24 m²(2)For rectangular section:(a)Pillar section: 2.63 × 2.66 m² and room of 4.10 × 4.14 m²(b)Pillar section: 3.08 × 3.08 m² and room of 5.24 × 5.24 m²

Items ‘a’ refer to pillars affected by wave energy, whereas pillars in item ‘b’ do not consider dynamic energy (for both sections).

3.6. Formulation of the Optimal Sizing Problem vs. Conventional Sizing Problem

Maximize

Subjected to

Pillar strength constraint on the ‘x’ axis is

Pillar constraint on the ‘y’ axis is

Foundation stability constraint is

Stability constraint of the slab on the ‘x’ axis is

Stability constraint of the slab on the ‘y’ axis is

Coherence constraint is

3.7. Simulation with FLAC 3D (Fast Lagrangian Analysis of Continua in 3D)

It should be highlighted that the Mohr–Coulomb failure criterion was adopted for assessing this model (see Figure 4). The values related to the geometric data are found in the algorithm, as previously described in Section 3.5.7; therefore, the geotechnical information and the wave equations are inserted in the script (Sections 3 and 3.2). The FISH language of FLAC 3D allows creating routines to generate the model interactively; this simplicity facilitates scenario analysis, as well as providing data and obtaining answers quickly. The routine that characterizes the numerical model is established as follows:(1)Input data obtained from the Mathcad algorithm (mathematical modeling)(2)Geometry construction, constitutive model, properties, and constraints(3)Mechanical response (output) of the rock material to the total load stresses(4)Evaluation(5)Analysis with established safety factor (1.3 for the pillar and 1.25 for the slabs)(6)Validation

Figure 4 

Simulation model in process.

4. Results

The described model and the details in Section 3.5.1 are placed in the standard format for mathematical nonlinear programming, as given in Section 3.6. The objective function structure with the restrictive equations represents the sentences that compose the algorithm (model), whose equations depict physical behavior for the mine pillar under static load and dynamic stress. The analytical model described was analyzed in the Mathcad V. 14.0 software and the subsequent evaluation with the FLAC 3D software, for which the variables mentioned above were considered (Sections 3.5.2 and 3.5.7).

On the other hand, the responsible velocity for the dynamic stresses corresponding to the absolute values was analyzed in the centroid of the pillar support section, and these found velocities were analyzed using the numeric methods given by the Newton theory of fluxions and contrasted with the I-Blast V. 8 software, as detailed in the Methodology Section. In Table 5 are the optimal values for the dimensions of the geometric arrangement of the pillar and a discriminated analysis for a pillar with square and rectangular cross section support.

Table 5 

Optimal sizing of rooms and pillars support when the wave affect is or is not considered.

As shown in Figure 5, the dynamic stresses due to the seismic wave represent an additional maximum load of 5.99 MPa; therefore, the pillars designed considering these dynamic stresses are more robust; this represents on average 3.7 tons more ore mass contained in each pillar compared to the pillars conventionally designed.

Figure 5 

Dynamic stresses in the square and rectangular pillar section.

5. Conclusion

As presented in Table 5, the results corresponding to the recovery factors for the mining arrangement, when considering the wave stress, are relatively inferior in relation to the conventional dimensioning, in which the dynamic stress delivered by wave pulses is neglected. This difference is due to the fact that the conventional sizing follows simplified criteria; hence, its attractive recovery (recovery factor) is detrimental to the pillar quality to respond to stress request. The dimensioned pillars, in a conventional manner, will offer an average recovery of 1.6% more compared to pillar dimensioning, when wave stresses are considered. The percentage of extra recovery represents a value of 3.7 tons on average per pillar; this mineral mass is an apparent gain in the dimensioned pillars.

As shown in Section 4, the dynamic stresses are around 4–5.99 MPa, a value under the mathematical model; although these stresses are much lower than the strength of the pillar, when added to the lithostatic stresses, they could create plasticization processes in several portions in the weakened structure of the pillar. Thus, compressive stress values were also obtained (Figures 6(a) and 6(b)) under the values obtained from the algorithm using Mathcad (see Table 5) in both sections (square x rectangular).

(a)

(b)

Figure 6 

Pillar under combined stresses (lithostatic and dynamic stresses). (a) Square section dimensions obtained with the Obert and Duvall equation (1967), as in Section 3.5.4 (2.66 m × 2.66 m). (b) Rectangular section dimensions obtained with the Mark and Chase equation (1997), as in Section 3.5.4 (2.66 m × 2.63 m).

All velocities were calculated through I-Blast software and hydrodynamic theory application using Mathcad algorithm. Both methodologies demonstrated coherence. P-wave velocity: 424.5 mm/s and 444.64 mm/s (I-Blast and hydrodynamic theory, respectively). Such velocities refer to information on the pillar centroid. Furthermore, velocities answers were contrasted and verified with heat map (Figure 7, energy released). Both cases also provided similar information 1 Value of 2 for the ANFO explosive agent.

Figure 7 

Energy released after sequential detonation (4.01 m × 2 m section).

Pillars with rectangular support cross sections represent the best support option as compared to square section pillars, since the load is better distributed in the former. The wave stress (P-wave) and the dimensions of the widths of the room spans are inversely proportional; these results are in agreement with the velocitythat is contrasted by the simulation and the numerical analysis of the potential velocity. In both cases (support), there is a difference of 0.73% and 0.66% when the wave effect is or is not considered severally (see Table 5).

Data Availability

The information and data that made this study possible will be made available by the authors upon reasonable request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors are grateful to DNA Blast group who provided information and I-Blast software for this paper.