# An optimal design of blasting

Q. Zeng, Pierre Navidi, Joseph Zarka

► To cite this version:

Q. Zeng, Pierre Navidi, Joseph Zarka. An optimal design of blasting. Rock Fragmentation by Blasting, 1996, Paris, France. 10.1201/9781003078104-8. hal-00116435

HAL Id: hal-00116435

https://hal.science/hal-00116435v1

Submitted on 28 Nov 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# An optimal design of blasting

Q. Zeng, P. Navidi & J. Zarka

Laboratoire de Mécanique des Solides, Ecole Polytechnique, Palaiseau, France

**Abstract:** Rock blasting is often considered to be still an art performed by only a few good experts. Numerical simulations are one of the means considered to help them. However, many aspects of such numerical simulations need to be improved such as the explosive model, the rock's response and of course, their coupling. This paper provides a review of blasting simulations. Some development are presented, especially in the rock damage model. It is also shown that an optimal design of blasting can be reached thanks to automatic learning techniques coupled to optimization tools which are discussed.

## 1 Introduction

Blasting is one of the basic operations in mining and quarrying industry. An efficient blasting for given geology conditions, rock volume to break and security conditions, results from a correct choice of explosive types, explosive quantity and the blasting planning. One can think of diameter, and length of boreholes, drilling pattern, initiation modes, firing sequence, delay time... Efficiency is then determined by the quality of rock fragmentation and a minimum cost. In this sense it is an optimization problem. At present, almost all of the blast design is based on experiences, at first because the behavior of the explosive and the surrounding fractured rocks are still not well known, secondly because of the complex coupling effects. However, the computer aided design techniques, numerical blasting models and expert systems can considerably reduce both the cost and duration of blast design while increasing the quality of prediction. This is the motivation of our research.

The key to a good blast design is a precise prediction. The first part of our work focus on the prediction by numerical simulations, in which the modelling is the discussion center.

The second part is on to the optimization of blast design. After introducing the general methodology of our approach we show an application example for blast design.

## 2 Simulation of blasting

Many computational models of rock blasting are based on data obtained from experiments or/and empirical relations derived from experimental data and analytical results for very simple cavities problems [27],[7],[13].... They are either too much simplified or very specific. Recently, more fundamental models based on the material constitutive laws, usually implemented in a finite element or finite difference code, have been attempted [16],[10],[12]....

From a mechanical point of view, the blasting problem couples two difficult fields : the detonation of explosives and the rock mechanics. At present, the two fields are still in course of development. This section begins by a brief introduction of an existing 3D computer code we used for our simulations. An explicit scheme algorithm enable us to handle the non-linear aspects and the fluid-structure coupling. Then a review of the models for the explosive and the rock is given. After, we introduce a new damaged rock constitutive relation. Finally some numerical results based on this new model are presented to show its tractability.

### 2.1 Computer code description

The program used can simulate most of the nonlinear dynamic transient phenomena [23]. This code can handle solid and fluid mechanics problem and the coupling between them. The main features used for our specific problems are :
* explicit central-difference time-integration algorithm

* Arbitrary Lagrangian-Eulerian (ALE) formulation [3] for fluid/structure interaction

* artificial viscosity in order to smear shock fronts[29]

* Incremental form of constitutive relations

## 2.2 Detonation model

Computation of the detonation propagation in explosives is a complex problem. It results from strong chemico-physics interaction during the detonation process[14],[6].... When a coupling with a structure is considered, the numerical modeling of the detonation is often simplified. This is because the response of the structure is the main interest. Such models are often of continuum type and admit the following hypotheses :

1. reduction of the reactive fluid in a fluid with two components, an initial substance and a final substance (explosive reaction products);

2. each of the two substances is separately in thermo-chemical equilibrium;

3. quasi Chapman-Jouguet detonation (a quasi-sonic detonation velocity).

These hypotheses reduce the complex computation to a simple fluid flow one with discontinuities which are then often smeared as follows :

* an artificial viscosity is used for the thermophysical variables of discontinuity (shock wave);

* a decomposition law is used for the thermochemical variables of discontinuity.

Therefore, the behavior of the explosive is represented by an equation of state and a decomposition law.

Often, in case of the detonation/structure simulation, the equation of state has a pseudo-potential form $p(v, E)$ (pressure-volume-energy) to avoid the thermo-chemical iterative calculation. The decomposition law is of «bulk burn» type due to the homogeneity consideration. A variety of explosive behavior ranging from the simple to complex models are available, we use the following models[30]:

* The Jones-Wilkins-Lee (JWL) equation of state for the explosive products pressure :

$$p = A(1 - \frac{\omega}{R_1 V})e^{-R_1 V} + B(1 - \frac{\omega}{R_{21} V})e^{-R_2 V} + \frac{\omega E}{V}$$

where $V = \frac{v}{v_0}$ : relative volume

$p$ : explosive products pressure

$E$ : detonation energy per unit volume

$A, B, \omega$ :

,R_1, R_2$ : constant parameters

* if $A = B = 0$, we refind the simple «gamma-law» with

$$\omega = \gamma - 1 (\gamma = \Gamma = \frac{\partial \ln \rho}{\partial \ln v})$$

* the «C.-J volume burn» chemical decomposition law :

$$m = \frac{1 - V}{1 - V_{CJ}}$$

where the burn fraction $m$ is the indicator of decomposition progress.

* the current pressure is calculated by :

$$P = m \cdot p$$

The detonation velocity, the pressure at Chapman-Jouguet point $P_{CJ}$ and the initial energy $E_0$ should also be known.

After having done several simulations (explosion in water, explosion in an elastic medium), we noticed some interesting results [31] :

* computation with the equation of state «gamma law» has a similar result compared to the computation with the equation of state JWL.

* the explosive energy partition (shock energy/bubble energy) can be characterized by the coefficient $\gamma$ of the «gamma law», the greater $\gamma$ is, the greater is the shock energy.

* the modeling of blast loading by using defined pressure-time history is only a crude approximation specially for the case of the point initiation.

## 2.3 Rock mechanical model

Another essential point for the blasting simulation is a proper representation of the rock behavior. The approach which involves the formulation of appropriate material models (constitutive relation) is preferred. The main phenomena that must be accounted for is crack initiation and propagation. Two alternatives are then possible. One is the discrete model which takes into account the physical presence of the cracks in a direct approach[11],[8].... For the computation, this one is often very costly and suitable for the post-period of the blasting processus. Another way is to take a continuum damage model. We have chosen
the last method which is more consistent with a finite element approach. A continuum model [31] is therefore developed and implemented into the explicit code. This is an elastic anisotropically damaging behavior based on the damage mechanics theory introduced by Kachanov [20].

The constitutive relation is :

$$ \sigma = \bar{C}(D) : \epsilon $$

where $\bar{C}$ is the elastic effective stiffness tensor and D is the damage tensor.

In an incremental form :

$$ d\sigma = \bar{C}(D) : d\epsilon + d\bar{C}(D) : \epsilon \quad (1) $$

The main hypotheses are the following:

1. the dominant mode of failure is mode I

as proven by numerous experiments [21],[22],[15]....

2. the material non linearity is totally governed by damage and thus plastic deformations are disregarded

3. The expanding plan directions of the penny-shaped cracks are perpendicular to the tensile principal stress directions.

this concept is deduced from the impact experiment results [9].

## 2.3.1 Damage description

The nature of the rock flaws and the possible complex loadings require an anisotropic damage description. Generally, an eight order tensor should be introduced to describe the anisotropic damage, but the difficulty to identify all the parameters is obvious. It is thus desirable to reduce the order of that tensor. Based on the third hypothesis above, a symmetric second-order damage tensor D is introduced, reducing the model to an orthotropic one.

We can consider each of the eigenvalues $D_i$ as a group of parallel cracks perpendicular to the associated eigenvector. As if :

$$ D_i = \frac{4}{3} \pi \cdot f \quad (2) $$

where $f = N_\alpha \cdot c_\alpha^3$ is the cracks ($N_\alpha$ cracks of radius $c_\alpha$) group density.

This definition is the same as the one for a homogenization method [25], [19] so that the results of homogenization can be used to determine the effective elastic matrix which is going to be shown lately.

$D_i$ can be also defined with the main parameters as introduced by Grady and Kipp[16] :

$D_i = N\nu$; with $\nu = \frac{4}{3}\pi \cdot c^3$ ($c$ and $N$ are mean values)

## 2.3.2 Damage evolution law

Generally, the second thermodynamic principle is used to determine the damage evolution law as in an elastoplastic problem. However the dissipative energy is not easily established. Grady and Kipp have developed an evolution law based on the activation and growth of an initial Weibull distribution of fracture-producing flaws, for an isotropic behavior [16]. It can be extended to an orthotropic material :

$$ \dot{D}_i = \begin{cases} \left[ \frac{8\pi(m+3)^2}{(m+1)(m+2)} \right]^{1/3} c_g n\left(\frac{\sigma_p^p - \sigma_u}{E}\right)^{\frac{1}{3}} D_i^{\frac{2}{3}} & \sigma_i^p > \sigma_u \\ 0 & \text{else} \end{cases} \quad (3) $$

with $i = 1,3 \quad n(x) = kx^m$

where

$c_g$ : constant fracture growth velocity

$\sigma_i^p, \sigma_u$ : principal stress and a constant depending on the material cohesion

$E$ : Young's modulus

$k, m$ : constant Weibull parameters

## 2.3.3 Damage criterion

The damage criterion is usually derived from the Griffith theory [17], [4], with the damage evolution law (3) shown above, we just take the following simple form :

$$ \sigma_i^p > \sigma_u $$

## 2.3.4 Effective elastic matrix

Based on the results of the homogenization methods [25], [19], the symmetric effective elastic matrix takes the form :

$$ \begin{aligned} [\bar{C}] &= \begin{bmatrix} [\bar{C}^{(1)}] & 0 \\ 0 & [\bar{C}^{(2)}] \end{bmatrix} \\ \bar{C}_{11}^{(1)} &= \frac{E}{\alpha}(1 - \nu^2 d_2 d_3) d_1 \quad \bar{C}_{22}^{(1)} = \frac{E}{\alpha}(1 - \nu^2 d_3 d_1) d_2 \\ \bar{C}_{33}^{(1)} &= \frac{E}{\alpha}(1 - \nu^2 d_1 d_2) d_3 \quad \bar{C}_{12}^{(1)} = \frac{E}{\alpha}\nu(1 + \nu d_3) d_1 d_2 \\ \bar{C}_{13}^{(1)} &= \frac{E}{\alpha}\nu(1 + \nu d_2) d_1 d_3 \quad \bar{C}_{23}^{(1)} = \frac{E}{\alpha}\nu(1 + \nu d_1) d_3 d_2 \\ \bar{C}_{11}^{(2)} &= 2\mu \frac{(1 - \beta D_2)(1 - \beta D_3)}{2 - \beta D_2 - \beta D_3} \\ \bar{C}_{22}^{(2)} &= 2\mu \frac{(1 - \beta D_1)(1 - \beta D_3)}{2 - \beta D_1 - \beta D_3} \\ \bar{C}_{33}^{(2)} &= 2\mu \frac{(1 - \beta D_2)(1 - \beta D_1)}{2 - \beta D_2 - \beta D_1} \end{aligned} $$

$\bar{C}_{ij}^{(2)} = 0 \quad (i \neq j)$

with

$d_i = 1 - D_i \quad (i = 1,3)$

$\alpha = 1 - \nu^2(d_1 d_2 + d_2 d_3 + d_1 d_3 + 2\nu d_2 d_3 d_1)$

$\beta = (1 - \nu)(2 - \nu) \quad \nu : \text{poisson's ratio}$

The result is in $0(f^2)$ order to the constant strain or constant stress method and the formulae are easy to calculate.
### 2.3.5 Incremental algorithm

During each time step the following algorithm is used for the determination of the stress and the damage tensor values : given $\sigma(t)$, $\epsilon(t)$, $D(t)$, $\dot{\epsilon}(t + dt)$:

$$ 
\begin{aligned}
D(t + dt) &= D(t) + \dot{D}(\sigma(t)) \\
d\tilde{C} &= \tilde{C}(D(t + dt)) - \tilde{C}(D(t)) \\
\sigma(t + dt) &= \sigma(t) + \tilde{C}(t) : \dot{\epsilon}(t + dt)dt + d\tilde{C} : \epsilon(t)
\end{aligned}
$$

### 2.3.6 Dominant fragment size

The calculation of the dominant fragment size (fragment size corresponding to the largest volume fraction of material) is also derived from the result given by Grady and Kipp which is strain-rate dependant. For sake of simplicity the formula is now fracture stress dependant :

$$ L_M = \eta(k, m) \left( \frac{\sigma_c}{E(m+3)} \right)^{-\frac{m}{3}} $$

where $\eta$ is a constant and $\sigma_c$ the tensile fracture stress.

## 2.4 Simulation results

The damage evolution law based on the Weibull distribution leads to predictions which are in good agreement with the strain-rate effects discussed by different authors[16],[2][10]...). In order to show the capability of the model proposed to take into account anisotropic effects, the following numerical experiment is done: we consider several cases of constant strain rate loading on a single element. The first loading is carried out up to a value of $D_1 = 0.4$ (largest eigenvalue of the damage tensor). After changing the loading direction, a second loading is applied to the same specimen, the curves in dot lines correspond the results for various directions of the second loading respectively 0°, 30°, 45°, 90° Fig(1).

### 2.4.1 Calculations in plane strain with an imposed pressure function

We use the model presented above to carry out the analysis on various geometries. One of them is a circular cylindrical cavity problem. Due to symmetry a quarter of the geometry is considered. The configuration and the loading, a sine function pressure imposed inside the cavity, are given in Fig 2 and 3:

The surrounding rock is an oil shale which has the properties :

<table><thead><tr><th>ρ(1/m³)</th><th>E (MPa)</th><th>ν</th><th>m</th><th>k(1/m³)</th><th>ε<sub>g</sub>(m/4)</th></tr></thead><tbody><tr><td>0.002</td><td>1.065 · 10<sup>4</sup></td><td>0.4</td><td>8</td><td>1.7 × 10<sup>18</sup></td><td>1300</td></tr></tbody></table>

In Fig 4 and 5, the distributions of the largest principal value of the damage tensor are respectively

<figure><img src="image_2.png" /><figcaption>Figure 1: Anisotropic effect</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Figure 2: Configuration</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3: Pressure function</figcaption></figure>

shown at time $t = 0.25ms$ et $t = 0.5ms$. This result underlines that :

* the damage begins around the cavity and progresses with decreasing value.

* a second damage area appears when the stress wave reaches the free surface and is reflected.

* the eigenvector directions of $D$ tells us that the
<figure><img src="image_5.png" /><figcaption>Figure 4: First principal value of D at $t = 0.25ms$</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5: First principal value of D at $t = 0.5ms$</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Configuration of two blasts</figcaption></figure>

cracks in the first area are radial and those in the second are parallel at the free surface.

To see the influence of blasting with delay time, two calculations are carried, one with simultaneous detonations and the second with a delay $\tau = 0.6ms$. The configuration is shown in Fig 6 with $T = 0.05(ms)$.

The results in Fig 7) and (8 are given for the largest eigenvalues of the damage tensor at $t = 1.5ms$ where the damage results are convergent.

We observe that the detonation with delay modifies the damage distribution and increases the damage volume (10% for this example).

<figure><img src="image_8.png" /><figcaption>Figure 7: simultaneous detonation</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 8: with a delay time $\tau = 0.6ms$</figcaption></figure>

### 2.4.2 3-D calculation coupled with detonation

The configuration is shown in Fig (9) :

<figure><img src="image_10.png" /><figcaption>Figure 9: Configuration in 3-D</figcaption></figure>

The explosive is idealized. Its brisance is near the ANFO. The equation of state used is the «gamma law» where the parameters are as follows :

<table><thead><tr><th>γ</th><th>ρ₀(1/cm³)</th><th>E₀(MPa·cm³/cm²)</th><th>p<sub>cj</sub>(MPa)</th><th>D<sub>cj</sub>(cm²/µs)</th></tr></thead><tbody><tr><td>2.8</td><td>0.9</td><td>770</td><td>2.765×10³</td><td>0.4</td></tr></tbody></table>

The space between the explosive and the surrounding rock is occupied by air. the decoupling ratio (borehole diameter to explosive diameter) is 1.2, the confinement is perfect (the stemming is also occupied by
rock). A point detonator is put at the bottom of the borehole. The rock is the same as in the first example (4).

In Fig (10), the damage result ( largest eigenvalue of the damage tensor) is shown at $t = 0.375$ ms,

<figure><img src="image_11.png" /><figcaption>Figure 10: damage result in 3-D</figcaption></figure>

# 3 Blast optimization

The blasting simulation based on the advanced modeling allows us to predict more generally and more precisely the blasting results. It can be used also to better understand fragmentation mechanisms, the roles played by the different design parameters [31]. However, It can not be yet a perfect tool for the blast design mainly because:

* it is difficult to simulate the entire blasting process with a single model often much simplified.

* the strong non linearity and the local/global problems in blasting make the computation too expensive. Even if it is technically possible to simulate the entire blasting process by more sophisticated modeling (e.g.: a combined finite element/discrete element model [8]), it is practically impossible to simulate the real sequence blasting problem.

In this context, the use of the automatic learning technics for searching rules (or models) seems to be a rational solution, and opens a realistic way to the blasting optimization.

In this section we will first give a brief introduction to the automatic learning technics. The methodology of blasting optimization is then discussed. Finally an example of automatic learning and optimization are given.

## 3.1 Automatic learning

In many technical domains, experts are aware that, very often, there is no a complete solution for a real world problem. But there is a possibility to build

a data base of examples. Each sample of this base is obtained experimentally or numerically with sometimes some fuzzy or missing informations. The automatic Learning Expert Systems generator (LEES), developed at LMS during 1986-1990, can build a set of rules based on the description of such a examples base. Several applications have been shown [26], [24], [28], [18], [5].... The main difficulty is to provide a good description of the examples base. The initial descriptors (let us call descriptor each single field of the data base) appear then as variables for the final model.

Suppose that the expert defines a set of these initial descriptors describing the observed and/or measurable phenomena, ($x_i$) generally in a limited number (e.g. the radius of the borehole, explosive energy etc...). Indeed only a subset of this initial set is sufficient to make the description of one example but it is not known a priori. These descriptors are then integrated within our actual limited knowledge to determine a set of more complex ones which are the intelligent descriptors, ($X_i$) (e.g. results of a numerical simulation or the ftime-history of the radius of the bubble in the water blast test, etc...). These descriptors are function of the initial descriptors and some of them may be considered among the conclusions ($C_k$) and then also be learned (e.g. the damaged rock volume).

In order to build the full knowledge of the problem or the model, the expert must prepare a file containing for each example :

The input descriptors:

$x_1$ | ... | $x_6$ | $X_1$ | ... | $X_{40}$ | ... |

The output descriptors:

$C_1$ | $C_2$ | ...

The symbolical option of LES will give a set of rules while its numerical option will generates a polynomial based model which is often easier to handle.

### 3.1.1 Generating rules

For a simple problem, the data base, with the initial descriptors ($x_i$) and the conclusions $C_i$, can be presented to the system. The system generates then a mathematical expression based on a polynomial expression :

$$C_I = \text{POL}(x_1, x_2, \dots, x_n)$$

For more complex problems, it is better to use the data base with the intelligent descriptors ($X_i$) in order to help the system to find a more reliable model for the conclusion $C_I$ :

$$C_I = \text{POL}(X_1, X_2, \dots, X_k).$$

Once the model is created, new unknown examples may be introduced for evaluating the conclusions.

The building of the set of ($X_i$) requires making them intrinsic in order to be able to reproduce them
for various types of problems independently of the nature or type of problem considered.

## 3.2 Optimization

Many optimization techniques are available. But we have the following points in mind:

* Design variables can be discrete and also take their value in a finite set.

* From a technological point of view many solutions can lead to the same type of results. This means that if $C_I$ is the cost function then in the design variable space POL can have flat domains (with no gradient). Many optima can exist simultaneously that is POL does not define a convex.

* Constraints generally varies depending on the domain in which the design set is looked for.

Blasting optimization by the traditional methods based on the mathematical programming theory such as quasi-Newton, conjugate gradient... seems to be impossible. On the other hand, the evolution methods such as the genetic algorithms are more suitable for a general optimization problem.

Contrary to the traditional methods based on the continuum concept, the genetic algorithms make evolve a population of solutions with a manner suggested by the genetic. Their advantages compared with the traditional methods are as follows:

* the probabilistic characteristics can handle local optima.

* there is no other restriction for the functions (cost functions, constraint function), excepted to be calculable.

* the algorithm is extremely general.

* design variables can be of any type.

The main disadvantage is costly computation time. Optimization by coupling the genetic algorithms with automatic learning is then a good strategy because the generation of the specific rules (for the cost function, the constraint functions) by automatic learning allows to compensate the handicap of the genetic algorithms method.

It is interesting to observe that the fact of expressing all the relations in a simple polynomial way enable us to carry out genetic optimization very effectively.

The other interesting points of this approach can be summarized as follows:

* It gives us a way of dealing with inverse problems, that cannot be handled by classical schemes.

* It lets us understand the essential descriptors which drives the desired conclusion. This can be a starting point for experts to build new PHYSICAL or MATHEMATICAL MODELS.

* By a clever choice of the intelligent descriptors, it is possible to extrapolate models resulting from a reduced set of simple samples to more complex situations.

## 3.3 Generation a damaged volume rule by LES

In this section we use the above methodology in a simple blasting design problem. Suppose that the explosive type and the rock site are given. We search a damaged volume $v_d$ (broken rock) expression in function of the drilling pattern parameters in order to carry out next the optimization calculation.

The database is built by numerically simulating blasting in plane strain. The chosen descriptors are:

borehole diameter $d$; burden $B$; decoupling ratio $R_d$ and two ratios $d/R_d$, $d/B$.

The result in Tab(1) is given by eleven simulations

<table><caption>Table 1: damaged volume database</caption><thead><tr><th>d(m)</th><th>B(m)</th><th>R<sub>d</sub></th><th>v<sub>d</sub></th><th>d/R<sub>d</sub></th><th>d/B</th></tr></thead><tbody><tr><td>0.050</td><td>0.600</td><td>1.0</td><td>1.952774</td><td>0.05</td><td>0.0833</td></tr><tr><td>0.050</td><td>0.5066</td><td>1.677</td><td>0.79336</td><td>0.0298</td><td>0.0987</td></tr><tr><td>0.050</td><td>0.600</td><td>1.172</td><td>1.467875</td><td>0.0427</td><td>0.0833</td></tr><tr><td>0.050</td><td>0.600</td><td>1.667</td><td>0.675639</td><td>0.0298</td><td>0.0833</td></tr><tr><td>0.040</td><td>0.500</td><td>1.3333</td><td>0.94768</td><td>0.03</td><td>0.08</td></tr><tr><td>0.040</td><td>0.600</td><td>1.333</td><td>0.871909</td><td>0.03</td><td>0.0667</td></tr><tr><td>0.100</td><td>0.750</td><td>1.257</td><td>1.520870</td><td>0.0796</td><td>0.133</td></tr><tr><td>0.100</td><td>0.750</td><td>1.474</td><td>1.265204</td><td>0.0678</td><td>0.133</td></tr><tr><td>0.200</td><td>1.000</td><td>1.2067</td><td>1.570654</td><td>0.1657</td><td>0.2</td></tr><tr><td>0.200</td><td>0.750</td><td>1.2067</td><td>2.957848</td><td>0.1657</td><td>0.267</td></tr><tr><td>0.200</td><td>1.000</td><td>1.843</td><td>1.182392</td><td>0.1085</td><td>0.2</td></tr></tbody></table>

First we ask the LES to build a polynomial model. Limiting the system to terms of order 2 with cross terms. the LES output is:

$$v_d =$$

$$0.85B \cdot d - 2.44B \left(\frac{d}{R_d}\right)^2 + (24.134 - 45.6R_d + 23.17R_d^2) \frac{d}{B} \quad (5)$$
(and $e_m = 0.00477$.)
3.4 Optimization problem

The problem can be set as following :

For a given explosive type and rock in-site, we
search the design parameters values (d,R<sub>d</sub>,B) which
maximize the specific charge

$$
E_s = \frac{v_d}{v_e (\text{explosive volume})}.
$$

Taking back the result (equation (5) ) divided by
($\frac{d}{R_d}^2$) which is in proportion to the explosive volume,
the problem becomes :

$$
\begin{align*}
\max f = & \quad 4.256 B \frac{R_d}{d} - 88.95 B + (90.39 - 92.644 R_d + \\
& \quad 25.54 R_d^2) \frac{R_d^2}{Bd}
\end{align*}
$$

with the constraints :

$$
1.0 \le R_d \le 2.0
$$

$$
0.4 \le B \le 1.0
$$

$$
0.5 \le d \le 0.2
$$

$d/R_d \ge 0.03$ (the explosive cartridge diameter is
supposed superior to $0.03m$)

This "simple" problem is firstly resolved by one
of the traditional methods, then by the genetic al-
gorithm.

• Conjugate Gradient

With a starting point : {d = 0.04, B = 0.8, Rd = 1.0}, a solution has
be found such that :

$$
\{d = 0.04, B = 0.5, Rd = 1.09358\}
$$

$$
\Rightarrow f = 1187.1658
$$

With another starting point
: {d = 0.1, B = 0.8, Rd = 2.0}, a different solution
has be found such that :

$$
\{d = 0.06, B = 0.5, Rd = 1.09358\} \Rightarrow f = 995.79
$$

This means that their are at least two local so-
lution. In fact the solution found depends highly
on the starting point.

• Genetic Algorithm

With several different starting populations, the
following solution is found :

$$
\{d = 0.04, B = 0.5, Rd = 1.094118\}
$$

$$
\Rightarrow f = 1187.1654
$$

Now, suppose that borehole diameter *d* and the ex-
plosive cartridge diameter can only take some fixed
discrete values :

$$
\frac{d}{R_d} \in \{0.048, 0.05, 0.066, 0.084, 0.1, 0.15, 0.2\}
$$

$$
\frac{d}{R_d} \in \{0.03, 0.036, 0.04, 0.05, 0.64\}
$$

The solution found by genetic algorithm and chang-
ing R<sub>d</sub> by d/R<sub>d</sub> is:

$$
\{d = 0.048, B = 0.5, d/R_d = 0.045 \rightarrow R_d = 1.0667\}
$$

$\Rightarrow f = 980.96$

One can notice that these points where not in the ini-
tial database. These values can then be used for a
new numerical simulation (or experiment). The dam-
aged volume obtained is then compared with the one
predicted by the model. If they agree, the model gen-
erated by the LLS is consistent. Otherwise we can
add a new example to our data base and build a new
model. This work can be done iteratively until the
model is consistent.

4 Conclusion

In this paper, it has been shown that the modelling of
micro-craks by a continuum damage model may im-
prove the quality of the prediction. For the continuum
rock model considered, the mixed phenomenological
- microphysical approach seems to be tractable; the
microphysical approach allows to build a consistent
model thus limiting the number of parameters, it also
allows quantitative results for the fragmentation.The
phenomenological approach simplify the model and
decrease the computing time.

However, even if the blast prediction by simula-
tion appears to be a powerfull prediction tool, it can
only be an element of the industrial blasting design.
Indeed, many other relevant parameters and techni-
cal aspects cannot directly enter the mathematical
model. The use of the automatic learning lets us
overcome this step efficiently, as we showed it on a
simple example. The coupling between genetic algo-
rithms and automatic learning seems to be a powerful
potential approach to blasting optimization.

References

[1] *First International Symposium on Rock Fragmentation By Blasting*, Luleå, Sweden, 1983.

[2] L.G Margolin & T.F Adams. Numerical simulation of fracture. In *First International Symposium on Rock Fragmentation By Blasting* [1].

[3] C.W. Hirt A.A. Amsden and J.L. Cook. An arbitrary lagrangian-eulerian computing method for all flow speeds. *J. of Computational Physics*, 14:227-253, 1974.

[4] H.D. Bui. *Mécanique de la Rupture Fragile*, volume 1. MASSON, 1978.

[5] Michale Bulik. *Sur l'Optimisation de la Protection Parasismique*. PhD thesis, Université Paris VI, 1994.

[6] R. Chéret. *La Détonation Des Explosifs Condensés*, volume 1 of *Scientifique*. MASSON, 1988.
[7] W.I. Duvall. Vibration associated with a spherical cavity in an elastic medium. *U.S. Bureau of Mines*, RI(4692), 1950.

[8] A. Munjiza et al. On a rational approach to rock blasting. *Computer Methods and Advances in Geotechnics*, pages 857-862, 1994.

[9] D.A. Shockley et al. Fragmentation of rock under dynamic loads. *Int. J. Rock Mech. Min. Sci. & Geomech. Abstr.*, 11:303-317, 1974.

[10] P.J. Digby et al. Computer simulation of blasting-induced vibration, fracture and fragmentation processes in brittle rock. In *First International Symposium on Rock Fragmentation By Blasting* [1].

[11] R. Hamajima et al. Analysis for discontinuous medium considering elemental deformation. *Computer Methods and Advances in Geotechnics*, pages 877-880, 1994.

[12] T.F. Adams et al. Simulation of rock blasting with the shale code. In *First International Symposium on Rock Fragmentation By Blasting* [1].

[13] R.F. Favreau. Generation of strain waves in rock by an explosion in a spherical cavity. *Journal of Geophysical Research*, 74(17):4267-4280, 1969.

[14] W. Fickett and W.C. Davis. *Detonation*. University of California Press, 1979.

[15] D. E. Grady and M. E. Kipp. *Dynamique Rock Fragmentation*, volume 1 of *Academic Press Geology*, section 10. Academic Press Inc, second edition, 1987.

[16] D.E. Grady and M.E. Kipp. Continuum modelling of explosive fracture in oil shale. *Int. J. Rock Mech. Min. Sci. & Geomech. Abstr.*, 17:147-159, 1980.

[17] A. A. Griffith. The phenomenon of rupture dans flow in solids. *Philosophical Transactions of the R. Soc. A*, 221:163-170, 1920.

[18] J.M. Hablot. *Construction de solutions exactes en élastro-plasticité. Application à l'estimation d'erreur par apprentissage*. PhD thesis, Ecole Nationale des Ponts et Chaussées, 1990.

[19] A. Hoenig. Elastic moduli of non-randomly cracked body. *Int. J. Solids Structures*, 15:137-154, 1979.

[20] Mark L. Kachanov. Time of the rupture process under creep condition. *Izvestiya Akademii Nauk;SSSR, Otd Tekhn Nauk*, 8:26-31, 1958.

[21] U. Langefors and B. Kilhlström. *The Modern Technique of Rock Blasting*, volume 1. JOHN WILEY and SON, 1963.

[22] S. McHugh. Computational simulation of dynamically induced fracture and fragmentation. In *First International Symposium on Rock Fragmentation By Blasting* [1].

[23] MECALOG. *RADIOSS*. FRANCE, 1986.

[24] P. Navidi and J. Zarka. Clever optimal design of materials and structures. In *Second French-Korean Conference on Machine Learning*, 1993.

[25] S. Nemat-Nasser and M. Hori. *Micromechanics overall properties of heterogeneous materials*, volume 37 of *Applied Mathematics And Mechanics*. North-Holland, 1993.

[26] M. Schoenauer and M. Sebag. Incremental learning of rules and meta-rules. In *7th Int. Conference on Machine Learning*, 1990.

[27] J.A. Sharpe. The production of elastic wave by explosive pressure, part i - theory and empirical field observations. *Geophysics*, 7(2):144-154, 1942.

[28] Michel Terrien. Systèmes experts par apprentissage en contrôle non-destructif. In *Les systèmes experts & leurs applications*, volume 3, pages 439-447. EC2, 1991. Onzièmes Journées Internationales "Les systèmes experts & leurs applications", Avignon, 27-31 Mai 1991.

[29] J. Richtmayer Von Neumann. A method for the numerical calculation of hydrodynamical shock. *J. of applied physics*, 1950.

[30] Mark L. Wilkins. Calcul de détonations mono et bidimensionnelles. In *Symposium High Dynamic Pressures*, Paris, 1967. I.U.T.A.M., Dunod.

[31] Q. Zeng. *Optimisation et l'utilisation des ex-plosifs en Génie civil*. PhD thesis, Ecole Nationale des Ponts et Chaussées, 1995.