Stochastic mathematical model for rock cutting force generation

Introduction

For simulation-based modeling and optimization to be used more widely in mining-machine research and design, a sufficiently general yet relatively simple mathematical model is needed to describe the dynamic load components acting on working units during the cutting of a wide range of rocks. Such a model is an essential part of the integrated mathematical model, or digital twin, of a mining shearer and can be used for engineering analysis, design calculations, process simulation, and machine-parameter optimization [1].

Mean load levels acting on cutting tools and the energy demands of operating processes in mining shearers and drilling machines have already been studied in detail [2–4].

Considerable research attention has been devoted to mathematical modeling of rock cutting with allowance for dynamic load components. The cutting process is cyclic in nature: rock is separated from the rock mass in discrete portions, or chips, and this separation is accompanied by fracture of the newly formed core and disintegration of the detached volume [4, 5]. Studies [6, 7] analyzing cutting-force oscillograms showed that, over the duration of a single chip formation event, the cutting force varies linearly and can therefore be represented as a triangular pulse. This conclusion is supported by numerous experimental studies of rock cutting reported by different authors [8–10].

Reference [11] describes a model in which the pressure distribution in the contact zone between the pick and the rock is governed by elastic deformation beneath the pick and by irreversible deformation associated with rock crushing in this region due to the formation of a pre-fractured surface zone. Under this model, the occurrence of rock chipping at discrete moments in time is one of the causes of the oscillatory variation of the cutting-force components.

Numerous studies have modeled the rock cutting process using different numerical methods [12–14]. These include the finite element method (FEM) [15], the boundary element method (BEM) [14], the discrete element method (DEM) [16–18], and combined approaches such as FDEM [14]. These methods make it possible to simulate rock cutting processes with varying degrees of accuracy. At the same time, their application to mining-machine calculations remains problematic. Methods of this type involve a large number of parameters, typically about 10 to 20, many of which are difficult to determine either analytically or experimentally. Some can only be estimated through parameter identification, that is, model calibration. This is a labor-intensive procedure that requires a substantial amount of experimental data. In addition, implementation of such models usually relies on expensive commercial software packages such as ANSYS, LS-DYNA, and ABAQUS, most of which are developed by foreign vendors.

The rock cutting process is inherently random because the structure and mechanical properties of rock vary randomly along the cutting path [4–6]. A mathematical model of rock cutting must therefore be stochastic. Determining the parameters of such models requires extensive experimental data, which makes it essential to minimize the number of parameters that must be measured or identified. A reasonable compromise is therefore needed between the number of factors included in the model, the feasibility of estimating their parameters reliably, and the accuracy required to reproduce the random dynamic load components.

As is well established [3, 4], the load acting on a cutting tool during coal and rock cutting is resolved into three components: cutting force, feed force, and side force. For modeling dynamic processes in mining shearer power systems, the dynamic components of the cutting force are of primary importance. These components generate dynamic loads in the machine drive, and their spectra contain significant frequency components up to 40–50 Hz [4, 19, 20]. For most practical purposes, therefore, there is no need to reproduce the spectral composition of the cutting force acting on the cutting tool with high accuracy above 40–50 Hz.

It is also important to note that, in the formation of random dynamic loads in mining shearer power systems, the exact form of the cutting-force distribution for an individual pick is of little practical importance. According to the central limit theorem, the distribution of the load generated by several picks, reaching several dozen in modern shearers, approaches a normal distribution regardless of the distribution law governing the load acting on each individual pick.

These considerations provide a basis for developing sufficiently simple algorithmic models capable of reproducing, with the required accuracy, the random loads arising during coal and rock cutting.

The objective of this study is to improve the design quality and operating efficiency of mining shearers by developing an easy-to-use stochastic mathematical model of rock cutting as a component of the integrated mathematical model, or digital twin, of a mining shearer for engineering analysis, design calculations, process simulation, and machine-parameter optimization.

To achieve this objective, the following tasks were undertaken:

• A stochastic mathematical model of the cutting process was developed in the form of a flow of random events, namely elementary loading events and the corresponding chip formation events associated with rock fracture, requiring identification of only two or three parameters instead of the 10–20 parameters typically required by models based on different variants of the finite element method;
• Experimental studies of the cutting process were carried out using a shearer cutting a full-scale coal-cement block, and realizations of cutting force for a radial pick were obtained under different operating conditions;
• Statistical processing of the experimental data was performed, the model parameters were identified, and the adequacy of the model with respect to the real process was established.

Theory. Mathematical model

In this paper, rock cutting by a mining-machine tool is modeled as a flow of random events, namely elementary loading events involving individual volumes of rock mass. Each event ends with an element separating from the rock mass and subsequently disintegrating. The interval between successive elementary loading events in time or space is treated as a random variable.

The forces Z_i, generated by the interaction between the cutting tool and the rock during the i-th elementary loading event are represented, in the general case, by triangular pulses with random parameters:

where Δx_i(t) is the elastoplastic deformation of the rock volume involved in the given elementary loading event; C_ps,i is the linearized pseudo-stiffness¹ coefficient of the rock, characterizing the resistance of the material being fractured to penetration by the cutting tool; and P_chip,i is the force at which brittle fracture occurs and the rock volume under consideration separates from the rock mass.

The idealized graphical form of relationship (1) is shown in Fig. 1.

Fig. 1. Idealized dependence of force Z_i on deformation Δx_i in an elementary loading event 

In the general case, the parameters C_ps,i and P_chip,i depend on the properties of the rock being cut, the geometry of the cutting tool, and the cutting conditions. The parameters should, in principle, be treated as random variables that vary from one elementary loading event to the next. The statistical characteristics of C_ps,i are established either during parametric identification of the proposed mathematical models or through dedicated cutting experiments. The parameter P_chip,i, however, is handled differently because reliable and widely accepted procedures are available for calculating the mean load acting on a pick. Accordingly, in each elementary loading event, P_chip,i is selected so that the simulation reproduces the prescribed mean load level.

When modeling cutting forces, the appropriate independent variable is the distance traveled by the pick between the onset of two successive fracture events. The intensity of the random event flow, λ, is defined as the average number of chips formed per meter of cutting path. The literature reports λ values for coal and certain rock types [4]. As an approximation, λ can be estimated from cutting-force oscillograms by counting the average number of peaks per unit pick travel.

For modeling cutting forces in coal, the following relationship is proposed. It was obtained by generalizing our experimental results for cutting a coal-cement block with ZR4-80 picks mounted on the RKU-13 shearer and the data reported in [4]:

λ = 584000/(2920 + Z_mean), m⁻¹, (2)

where Z_mean is the mean cutting force, N.

To reduce the number of model parameters requiring identification, event flows were considered in which the time intervals, or equivalently the path intervals Δl, follow one-parameter distributions. For a Poisson event flow, the interval between events is a random variable with an exponential distribution described by the density function

 f(Δl) = λ exp (−λΔl).  (3)

Another one-parameter distribution considered here is the Rayleigh distribution, written as

where σ₀ is the distribution parameter determined from the expression σ₀ = 0.798/λ.

To assess whether these distributions are suitable for the problem under study, cutting-force oscillograms obtained in experimental tests of the RKU-13 shearer on a test rig with a coal-cement block were analyzed [20, 21]. Fig. 2 shows, as an example, the histogram of the interval distribution between successive chip formation events during cutting of the coal-cement block with a ZR4-80 pick mounted on the RKU-13 shearer.

Fig. 2. Histogram of intervals between successive chip formation events based on full-scale experimental data (n = 77)

Statistical analysis of the experimental data shows that the distribution of intervals between elementary events is fairly complex and is not well described by standard theoretical distributions. To adapt the exponential distribution to real cutting conditions, the possible range of the random variable Δl must be bounded from below and above in such a way that its mean value, equal to 1/λ, remains unchanged. The lower bound of the interval range, l_lower, is determined from l_lower = 1/(λK_lower), where K_lower is the ratio of the mean interval to the minimum distance between successive events. This parameter can be estimated from the oscillogram as the minimum spacing between cutting-force peaks, assuming the cutting speed is approximately constant. For cutting a coal-cement block, K_lower = 2.0.

The upper bound, l_upper, is determined by requiring the expected value of the truncated distribution to be equal to 1/λ. The probability density function of a random variable bounded both below and above is given by

Setting the expected value of this random variable equal to 1/λ yields the following transcendental equation for determining the upper bound, l_upper:

l_lower exp (−λl_lower) = l_upper exp (−λl_upper). (6)

The trivial solution l_upper = l_lower is disregarded, and the equation is solved numerically. For example, λ = 68 m⁻¹ and K_lower = 2.0 yield l_lower = 7.4 mm and l_upper = 25.9 mm. As shown in Fig. 2, these boundary values cover virtually the entire range of variation of the random variable Δl.

Fig. 3 presents the probability density functions of the random variable Δl for the three cases considered above.

Fig. 3. Probability density functions of the intervals between successive chip formation events:
1 – exponential; 2 – truncated exponential; 3 – Rayleigh

For the first case, the interval between successive chip formation events in the numerical simulation is generated as the random variable² [22]

Δl_i = −λ⁻¹ l_n (1 − ν), (7)

where n is a random variable uniformly distributed over the interval [0, 1].

In the second case, only values Dl_i obtained from Eq. (7) that satisfy the inequality

l_lower < Δl_i < l_upper, (8)

are retained.

In the third case, the interval between successive chip formation events follows a Rayleigh distribution. In the numerical simulation, this interval is determined from the following expression [21]

where ξ₁ and ξ₂ are realizations of independent normally distributed random variables with zero mean and unit variance.

For all cases, the pick position at the onset of the i-th elementary event is determined from

l_i = l_i – 1 + Δl_i. (10)

The elastoplastic deformation of the rock volume involved in a given elementary loading event is then

Δx_i(t) = l(t) – l_i,

where l(t) is the distance traveled by the pick by time t, that is, the pick coordinate.

The maximum deformation of the elementary rock volume at the moment of brittle fracture is ΔX = P_chip,i/C_ps,i. To simplify the model, C_ps,i and P_chip,i are assumed to be deterministic. For a prescribed mean cutting force, these quantities are taken to be the same for all elementary loading events and are denoted by C_ps and P_chip.

The mean cutting force, which can be determined by established methods [2–4], is used as the reference quantity for determining the required model parameters.

At any instant, the cutting force acting on the pick is defined as the sum of the forces generated in the individual elementary fracture events

where N₁ and N₂ are the sequence numbers of the first and last elementary fracture events that have started but have not yet ended. The summation procedure for the forces generated in the individual elementary fracture events is illustrated in Fig. 4.

Fig. 4. Summation of the forces generated in the individual elementary fracture events

For the proposed cutting-force model, the expected value of a realization of the random process over an arbitrary length L is determined as the ratio of the area enclosed by the polygonal line shown in Fig. 4 to the length L. This area is equal to the sum of the areas of the elementary triangles, whose number N is equal to the number of elementary loading events that occur over the trajectory segment under consideration. The area of each triangle is S_tri = 0.5 P_chip²/C_ps. Accordingly, the expected value of the cutting force is

Since λ = N/L, the final expression for determining the parameters of an elementary fracture event from the known mean value, that is, the expected value, of the cutting force can be written as

The generalized flowchart of the algorithm used to simulate rock cutting force is shown in Fig. 5.

Fig. 5. Generalized flowchart of the rock cutting force simulation algorithm

In Fig. 5, blocks 1 and 9 mark the beginning and end of the algorithm, respectively. Block 2 specifies the pseudo-stiffness coefficient С_ps and the coefficient K_lower. Block 3 calculates the mean cutting force using established methods. Block 4 determines the chip formation intensity λ and the maximum chipping force P_chip. In block 5, the random variables ν, ξ₁, and ξ₂ are generated, and Eqs. (7) and (8) or Eq. (9) are then used to obtain the random interval Δl_i, and, subsequently, the pick coordinate l_i. Block 6 calculates the elastoplastic deformation of the rock element being fractured, Δx_i(t). Block 7 uses Eq. (1) to determine the force associated with each chip formation event that has started but is not yet complete at the given instant. Finally, block 8 computes the cutting force at the time under consideration by summing these contributions according to Eq. (11).

¹ The prefix pseudo- is used to emphasize the elastoplastic nature of the resulting deformations.

² Bakalov V. P. Digital Simulation of Random Processes: A Textbook. Moscow: MAI Publishing House; 2001. 84 p.

Results

The proposed simulation model includes two main input parameters, λ and C_ps, which can be determined either through dedicated experiments or through parametric identification of the mathematical models developed for specific mining machines. In addition, the second variant of the model, which uses a truncated exponential distribution for the intervals between successive chip formation events, includes one additional parameter, K_lower.

Experimental studies of the cutting process were carried out on an RKU-13 shearer equipped with strain-gauge instrumented holders for measuring the three components of the interaction force between a radial pick and the rock [20]. Measurements were taken during cutting of a full-scale coal-cement block with a mean cutting resistance of 210 N/mm at different machine travel speeds.

To reproduce the experimental conditions as closely as possible in the simulation, the mean cutting force was set equal to its experimental value. During a single cut, this mean cutting force varied from zero to a maximum and then back to zero according to a sinusoidal law. With a total duration of 0.75 s for one cut, segments of the realizations from the time interval 0.175–0.575 s were analyzed. Oscillograms from five successive cuts, both in the full-scale and in the computational experiment, were joined together to form a composite realization 2 s in duration. For spectral analysis, the composite realizations were assembled from the original records after trend removal. The resulting composite realizations possess the property of stationarity and are sufficiently long in time to make the ergodicity assumption more justified. These realizations were subjected to statistical analysis, which yielded estimates of the expected value, standard deviation, distribution histograms, and spectral density of the random process.

Fig. 6. Fragments of oscillograms from the full-scale experiment (a) and the computational experiment (b)

Table

Cutting force, kN, in computational and full-scale experiments

Fig. 7. Histograms of the cutting-force in the full-scale experiment (a) and the computational experiment (b) at a mean load level of 5.67 kN

Fig. 8. Estimated normalized spectral densities of composite cutting-force realizations for the computational model and the full-scale experiment

The cutting force was simulated using the algorithm described above for three variants of the proposed model, corresponding to different distributions of the random variable representing the intervals between successive chip formation events: exponential in variant 1, truncated exponential in variant 2, and Rayleigh in variant 3. The closest agreement with the full-scale experimental data was obtained for variant 2. Fig. 6 shows fragments of oscillograms from the full-scale and computational experiments for this variant.

The table presents the statistical results for variants 2 and 3 obtained from the full-scale and computational experiments under two operating conditions with substantially different mean cutting forces.

Fig. 7 compares cutting-force histograms obtained in the full-scale and computational experiments for variant 2 of the model.

Fig. 8 shows plots of the estimated normalized spectral densities of composite cutting-force realizations obtained in the full-scale and computational experiments.

As shown in Fig. 8, the variance of the random cutting-force process is distributed similarly over the 0–100 Hz frequency range in the full-scale and computational experiments.

Discussion and conclusions

Analysis of the table shows that the smallest errors are obtained for variant 2 of the model, not exceeding 7% for the expected value and standard deviation and 15% for the maximum value. Variant 3 of the model also provides acceptable accuracy, with the corresponding errors reaching 18% and 29%, respectively.

Variant 2 of the developed model shows good agreement between the simulated and experimental cutting-force histograms, as illustrated in Fig. 7. As shown in Fig. 8, the spectral composition of the cutting force obtained in the computational experiment is also in close agreement with the experimental data.

The model validation shows that the developed stochastic model of cutting-force formation, specifically variant 2, agrees well with the experimental data in terms of the statistical characteristics, the distribution histogram, and the spectral composition at two representative mean load levels. Variant 3 yields somewhat larger errors, but it has an important advantage in that it requires only two parameters, λ and C_ps.

The modeling error of the developed mathematical model does not exceed the corresponding errors reported for different finite-element-based approaches and expensive commercial software packages, most of which have been developed abroad. For example, according to [16], the error of 3D DEM simulation for the standard deviation of sandstone cutting force exceeds 22%.

Given the simplicity of the proposed mathematical model and its acceptable level of agreement with real processes, the validated approach appears suitable for developing a stochastic mathematical model of cutting in hard soils by the working units of earthmoving machines, as well as for modeling the operating process of crushing machines [22].

Thus, three variants of a stochastic mathematical model of rock cutting have been developed. They involve only a small number of parameters, namely two or three, do not require expensive commercial software, and provide an adequate representation of the real process.

Of the three variants, the simulation model represented by variant 2 best satisfies the adequacy requirements. In this variant, the intervals between elementary loading events follow a truncated exponential distribution. The error in modeling the cutting force acting on an individual pick does not exceed 7% for the expected value and standard deviation and 15% for the maximum value. Good agreement is also observed between the distribution histograms and the spectral density plots obtained from the processing of full-scale and computational experimental data.

The proposed mathematical model of rock cutting can be used as a component of the integrated mathematical model, or digital twin, of a mining shearer for engineering analysis, design calculations, process simulation, and structural and parametric optimization at the design stage. The proposed representation of rock fracture as a flow of random events, namely elementary fracture events, also appears suitable for mathematical modeling of cutting in hard soils by the working units of earthmoving machines and for modeling the operating process of crushing machines.