Optimization of Fairhurst-Cook Model for 2-D Wing Cracks Using Ant Colony Optimization (ACO), Particle Swarm Intelligence (PSO), and Genetic Algorithm (GA) ()
1. Introduction
The study of rock failure under compression is a matter of importance in rock mechanics. Two main mechanisms, including axial splitting and shear failure, have been identified for this type of rock failure [1]. Failure by splitting mechanism can be related to the spalling phenomenon [2]. This phenomenon always occurs at high stress contact points, for example in a ball bearing. Spalling occurs where the maximum shear stress occurs below the rock surface.
Splitting parallel to the direction of maximum compression is a common type of macroscopic fracture of a brittle rock in the vicinity of its surface. Modelling of brittle failure is one of the greatest challenges of material failure analysis which results from the irreversible and very rapid propagation and connection of cracks, in a process called fracturing.
Materials like natural ice and rock are heterogeneous and crystalline with different behaviors under variation of applying forces. Generally, their behavior can be categorized in two main types, including ductile behavior under high confining pressures, and brittle behavior under low confining pressures [3]. Because most rocks are brittle at low temperature and low confining pressure, virtually most of the rocks at or near the Earth’s surface exhibits brittle failure mechanism. In other words, failure in these rocks occurs as deformation-induced loss of cohesion [4] [5] [6] [7].
One of the main mechanisms proposed for modelling of axial splitting is the sliding crack or “wing crack” model which was originally proposed for a 2D crack in a plate [3] [8] [9] [10]. Splitting failure begins with a primary crack (pre-crack) and then sliding mechanism creates secondary cracks or so called, wings, at the edges of those primary cracks. The macroscopic failure occur when a series of cracks extend and link together to split the material. The wedged crack can be modelled as a straight representative crack, where the wedging forces on the pre-crack area are opening the crack. This is often called the Fairhurst-Cook model [3]. The whole process of creating a 2D wing crack with Fairhurst-Cook model and an optimization process on this model is done by Dyskin et al. (1994). The basis of this optimization was the crack semi-length (l) and a critical point for this parameter found in this study which results in an unstable mode for the crack [11]. Another Optimization on the Fairhurst-Cook model is the purpose of this paper which is the optimization of wedging force based on the crack angle to find its critical points.
2. Ultimate Rock Strength
The ultimate rock strength and the orientation of the macroscopic failure plane depend on the confining pressure applying on the rock. Some models for micro-cracking under compression (e.g. [8] and [12]) have been proposed based on the fact that frictional sliding on pre-existing cracks results in the formation of tension cracks at their tips. It has been shown that under triaxial compression, the micro-crack distribution is almost uniform until the applied stress reaches the ultimate strength, and at the ultimate strength a region of high crack-density (tension cracks) emerges along a plane which eventually becomes a macroscopic shear failure plane [1] [13]. Variations of the “ultimate strength” and the orientation of shear failure with confining pressure using this model are already known and presented by Nemat-Nasser and Horii (1984) [1].
As discussed before, under uniaxial compression, the tension cracks grow at a sharp angle relative to the orientation of overall failure plane with confining pressure. Variation of this “crack angle” changes the amount of wedging force with specific but different trends, resulting different amounts of opening in wing cracks and “ultimate strengths” for rocks to fail under compression. Excluding these trends is another purpose of this paper.
3. Fairhurst-Cook Model
As it can be seen in Figure 1, under uniaxial compression, crack opens in the shape of a wing with a specific crack angle (α) and also a specific crack amplitude (a) (This parameter is also called “half-length” or “semi-length” in other studies; but here it is called “crack amplitude” which is more proper and differentiates it from the other “semi-length” or “l”).
The amount of this wedging force (f) which creates the wing crack can be determined from the following equation:
f = 2.a.σ.β(α) (1)
where “σ” represents the maximum principal stress, “a” represents the crack amplitude, “f” represents the wedging force and “β(α)” is a function of crack angle which can be determined from the following equation:
β(α) = sin2(α).cos(α).(1-tan(α).tan(μ)) (2)
In the above equation, “μ” represents the internal frictional angle of the rock.
4. Sensitivity Analysis
In this study, a complete sensitivity analysis performed to find the optimized amounts of parameters in Fairhurst-Cook Model. Three main parameters
![]()
Figure 1. Schematic of 2-D wing crack growth according to Fairhurst-Cook Model [14].
including crack amplitude (a), crack angle (α) and maximum principal stress (σ) and a fixed amount of 89 degrees for internal friction angle (μ) have been considered in this sensitivity analysis. According to this model, these parameters have the main roles in determining the amount of wedging force which creates a two-dimensional wing crack, while the internal friction angle (μ) should be determined according to the rock sample which is being cracked.
Three different stages have been considered in the sensitivity analysis process. In the first stage of this process, crack amplitude analyzed and this process repeated for the second and third stages by analyzing crack angle (α) and maximum principal stress (σ). Performing this sensitivity analysis in Matlab environment, results in Figure 2 as the output.
As it can be seen in Figure 2, it includes two different parts which are differentiated in a specific point. By focusing on this figure, the trend of changing the wedging force creating wing crack by variation of crack angle in different crack amplitudes can be excluded. This process is shown in Figure 3.
As it is shown in Figure 3, trends change in a certain point (α = 60.7˚). Before this point, in any crack angle, the amount of wedging force increases by increasing the crack amplitude, but this trend changes after this point (α = 60.7˚) and by increasing the crack amplitude, wedging force will decrease. This amount of crack angle (α = 60.7˚) virtually acts like a border line and specifies the maximum amount of crack angle which causes the increase of wedging force. In other words, this specific crack angle (α = 60.7˚) has considered to be the “peak point” for crack angle in analyzing the change of wedging force by variation of crack amplitude.
![]()
Figure 2. Change of wedging force (f) by variation of crack angle (α) in different crack amplitudes (a) by σ = 60 (Mpa).
Another fact which can be excluded from Figure 1 is the reversion of wedging force sign from positive to negative. The reason of this conversion is thought to be the conversion of applying wedging force from tensile to compression as a result of increasing the crack angle and decreasing the angle between force direction and its applying area, resulting in closing of the crack lips, instead of opening. In is thought that the reason for the variation of trends is also the same, indicating that the bigger wedging force still stays bigger after peak point angle, but with a different type of force and in an apposite way, applying compression force instead of tensile and closing the crack instead of opening it. This fact highlights the importance of finding the critical angle or so called, the “peak point”.
In the next stage of analyzing, the change of wedging force creating wing crack by variation of crack angle has been examined. For this purpose, one of the curves in Figure 2 has been chosen (a = 11 mm) and by focusing on it, the trend of changing the wedging force by variation of crack angle has been excluded. As it can be seen in Figure 4, this trend changes in a specific point (α = 37.7˚). Before this point, by fixing the amount of crack amplitude (a) and maximum principal stress (σ), the amount of wedging force increases as the amount of crack angle increases. In other words, there is a straight proportion between these parameters. After that specific point (α = 37.7˚) which is considered to be the “peak point” in variation of wedging force by changing the crack angle, that trend changes and in a reverse relation, the amount of wedging force decreases as crack angle increases. Just like the first stage, this specific amount of crack
![]()
Figure 3. Different trends of changing the wedging force (f) by variation of crack amplitude (a).
![]()
Figure 4. Different trends of changing the wedging force by variation of crack angle.
angle (α = 37.7˚) acts like a border line which limits the increase of wedging force by increasing the amount of the relevant parameter.
For the third stage of this sensitivity analysis which includes analyzing the wedging force variation by changing the maximum principal stress (σ), three different amounts of maximum principal stress considered in a fixed amount of crack amplitude (a). In this stage, the variation of wedging force by changing the crack angle in different amounts of maximum principal stress has been investigated and by its excluded trends, plotted in Figure 5 and Figure 6.
As it can been seen in Figure 6, trend of changing the wedging force by changing the crack angle in different amounts of maximum principal stress changes in a specific amount of crack angle (α = 60.7˚), like the previous stages. This point which considered to be the “peak point”, differentiates the trends and returns it from straight relation to the reverse. Just like previous stages. This peak point limits the increase of wedging force by increasing the relevant parameter.
According to this sensitivity analysis, optimized amounts of each parameter by fixing other relevant parameters have been excluded, but finding the optimum amounts of parameters without fixing others is much more important. For this reason, using some of the efficient algorithms of artificial intelligence for optimization is suggested. In this research three algorithms including Ant Colony Optimization (ACO), Particle Swarm Intelligence (PSO) and Genetic Algorithm (GA) were used and their results have been compared in order to find the most efficient method in general or separate parts (before and after the peak
![]()
Figure 5. Variation of the wedging force (f) by changing maximum principal stress (σ).
![]()
Figure 6. Different trends of changing the wedging force (f) by variation of maximum principal stress (σ).
point angle) and also to achieve the more reliable results.
5. Ant Colony Optimization (ACO)
Ant colony optimization (ACO) is a population-based metaheuristic algorithm which can be used to find approximate solutions to hard and discrete optimization problems. ACO is an algorithm for finding optimal paths that is based on foraging behavior of some ant species searching for food [15]. In this algorithm, a set of software agents called “artificial ants” search for good solutions to a given optimization problem [16]. At first, the ants distribute randomly to start searching for food. When an ant finds a food source, it walks back to the colony depositing pheromone on the ground in order to mark its favorable path which should be followed by other ants. In this situation, other ants are more likely to follow that path. As more ants find a path, it gets stronger and this process repeats so many times until there are a couple streams of ants traveling to various food sources near the colony. However, the pheromone trail starts to evaporate over time, thus its attractive strength reduces. The more time it takes for an ant to travel down the path and back again, the more time the pheromones have to evaporate [17]. Because the ants drop pheromones every time they bring food, shorter paths are more likely to be stronger, hence optimizing the “solution”. “Positive Feedback” eventually leads to all the ants following a single path. Ant colony optimization exploits a similar mechanism for solving optimization problems. To apply ACO, the optimization problem is transformed into the problem of finding the best path on a weighted graph. The artificial ants incrementally build solutions by moving on the graph. The solution construction process is stochastic and is biased by a pheromone model, that is, a set of parameters associated with graph components (i.e. nodes) whose values are modified at runtime by the ants. This process is repeated many times until the stopping condition is being satisfied and a satisfactory result is achieved. The most famous application of this algorithm is finding near-optimal solutions to the traveling salesman problem (TSP) [15] [18] [19].
6. Particle Swarm Intelligence (PSO)
Particle Swarm Intelligence originally formed based on the movement of organisms in animal swarms such as bird flocks or fish schools to simulate their social behavior [20]. A basic version of PSO algorithm starts by having a population (called a swarm) of candidate solutions (called particles). These particles are moved around the search-space to find the optimized solution (called global best solution) [21]. In this process, each particle determines its movement through the search-space by combining some aspects of the history of its own current and best locations with those of other members of the swarm, with some random perturbations. When improved positions discovered, these positions will guide the movements of the swarm. The process is repeated many times trying to find a satisfactory solution and is expected to find a global or near-global optimized answer at the end of all algorithm cycles, although it is not guaranteed.
This algorithm can be summarized in four main steps, which are repeated until the stopping condition is satisfied:
・ Assigning initially random positions and velocities for all of the particles in the search-space
・ Evaluation of the fitness of each individual particle
・ Updating the individual and global best positions
・ Updating the velocity and position of each particle [22]
In past several years, PSO Algorithm has been successfully applied in many researches and different application areas. It is demonstrated that PSO algorithm gives better results in a faster, cheaper way compared with other methods [23].
7. Genetic Algorithm (GA)
Genetic Algorithm (GA) is a heuristic optimizer algorithm based on the evolutionary ideas of natural selection and principles of “survival of the fittest” from “Charles Darwin”. Genetic Algorithms are commonly used to generate proper solutions to optimization problems by relying on bio-inspired operators such as “crossover” and “mutation”.
GA simulates the natural selection process among individuals (also called “phenotypes”), each of which representing a possible solution. Each possible solution has a set of properties and is known as a “chromosome” (also called “genotype”). Representing the solutions as binary strings of 0 s and 1 s are more common, but other types such as “real strings” are also applicable. The initially created individuals then lead through the process of evolution [24].
Starting with a randomly generated population of chromosomes, the evolution occurs as a process of fitness based selection and recombination to produce a better population in each iteration called a “generation” [25]. This process will be done by defining a proper “fitness function” and improving the initial population through repetitive application of the mutation, crossover and selection operators.
In each generation, the fitness of every individual in the population which is usually the value of the objective function in the problem, is being evaluated and the GA creates a new population by a new group of chromosomes with resulted fitness values. In this process, firstly, parents are selected to mate, based on their fitness, producing “offspring”, so better solutions with more fitness are given better chance to reproduce by crossover operation. The offspring inherit characteristics from both parents, but not equally. As parents mate and produce offspring, some new rooms must be freed for the newly generated chromosomes.
Since the population contains more information than each individual fitness, GA combines the good information hidden in a solution with good information from another one in the mating pool, in order to produce new solutions with good information inherited from both parents [24]. After the GA mates the new individuals and mutates some of them, the population undergoes a complete generation change. Some of the individuals in the population are replaced by new ones, since the size of the population has to remain fixed and shouldn’t increase. The population will then consist of offspring plus a few of the older individuals, which the genetic algorithm allows to survive to the next generation, because they are the best members of the population, called “elite individuals”. Following this procedure, it is hoped, but not guaranteed, that over many generations, better solutions will remain while the least fit solutions die out.
Each generation will contain, on average, more good genes than the previous one. Once the population is not producing much better solutions than previous generations, the algorithm is said to have converged to a specific set of solutions for the problem. Eventually, the algorithm terminates when either it converged to a proper solution with satisfactory fitness level or a specific number of generations has been produced [26] (Figure 7).
8. Optimization Process
As it has been discussed earlier, three different optimization algorithms including Ant Colony Optimization (ACO), Particle Swarm Intelligence (PSO) and Genetic Algorithm (GA) have been used in order to find the optimum or near optimum parameters influencing the wedging force which creates a wing crack. This optimization has been performed in several cases including general maximizing and minimizing of the wedging force and also separate maximizing and minimizing of this parameter, isolated in before and after the peak point angle (37.7˚).
The First algorithm was ACO and it has been performed by a colony of 46 ants trying to find the optimum amounts of crack angle (α), crack amplitude (a), maximum principal stress (σ) and frictional angle (μ) in 100 iterations. Two specific upper and lower limits have been considered for each parameter in order to specify their valid range of variation. This range was from 0 to 120 degrees for
![]()
Figure 7. Schematic of Genetic Algorithm.
crack angle, from 0 to 89.999 degrees for frictional angle, from 0.001 to 0.046 meters for crack amplitude and finally 10 to 100 mega-Pascal for maximum principal stress. The results of this optimization are illustrated in Table 1.
The second optimization algorithm was Particle Swarm Intelligence (PSO) and it has been performed in 50,000 iterations with 10,000 swarms and 4 particles for each swarm, representing crack angle (α), crack amplitude (a), maximum principal stress (σ) and frictional angle (μ). Just like the previous algorithm, two specific upper and lower limits with the same amounts have been considered for each parameter in order to specify their valid range of variation. The results of this optimization are also illustrated in Table 1.
The third and final optimization algorithm was Genetic Algorithm (GA). In this study, real genetic algorithm used and it has been performed with 200 chromosomes and 4 genes for each chromosome. These genes represented the related parameters including crack angle (α), crack amplitude (a), maximum principal stress (σ1) and frictional angle (μ) and each chromosome represented a resulting amount of wedging force (f). The optimization process repeated for 4000 iterations with the same upper and lower limits for different parameters and like the previous stages, its results is illustrated in Table 1.
These steps repeated several times, in order to achieve better results of wedging force without exiting the valid range for each parameter and as it can be seen in this table, Particle Swarm Intelligence (PSO) algorithm achieved better results with more amount of wedging force at the end. So these amounts of parameters will be more proper to use as the optimum amounts and this method is more efficient in general maximizing of wedging force in wing crack model comparing to the ACO and GA algorithms. These steps have been repeated for the general minimization stage with the same ranges for the parameters. The results are shown in Table 2 which indicate that GA is the optimum method in this stage of the optimization.
![]()
Table 1. Results of the optimization algorithms in general maximization of the wedging force.
![]()
Table 2. Results of the optimization algorithms in general minimization of the wedging force.
As it can be understood from Table 1 and Table 2, the GA and PSO algorithms both had proper results with near optimum amounts of wedging forces, but ACO didn’t show any acceptable results; so this algorithm has been neglected in the next stages of the optimization. The results of these stages of optimization are shown in Table 3 and Table 4.
For better understanding and easier access to the best algorithms and optimum amounts of parameters, a summary of these methods and results are shown in Table 5. Using these optimum amounts helps to accelerate or decelerate (in different cases of maximizing and minimizing) the opening of a wing crack by more applying wedging force in the tensile or comparison modes, avoiding the waste of time and resources.
9. Conclusions
In the first part of this study, a complete sensitivity analysis performed on the main parameters of Fairhurst-Cook Model to determine their relations and
![]()
Table 3. Results of the optimization algorithms in separate maximization of the wedging force.
![]()
Table 4. Results of the optimization algorithms in separate minimization of the wedging force.
![]()
Table 5. Summary of best algorithms and optimum results in different cases.
trends of change. Concluded results from this sensitivity analysis are summarized as below:
・ The trend of changes for each of these parameters reversed in a specific amount of crack angle (α) as a “peak point”.
・ This peak point is about 60.7 degrees in the analysis of crack amplitude and also maximum principal stress and it is about 37.7 degrees in crack angle analysis.
・ In all of above cases, straight relations between that specific parameter and the amount of wedging force observed before the peak point, but it has changed to reverse relations after that point.
・ The peak point virtually acts like a border line and limits the increase of wedging force by increasing the crack angle.
・ The reason for the reversion of wedging force sign is thought to be the reversion of force from compression to tensile with the same magnitude order; which means the biggest force still remains the biggest one, but in a reverse direction and it happens because of the changes in crack angle, making an angle of more than 90 ̊ between the direction of wedging force and its applying area.
Knowing the amount of peak points can be helpful for understanding the ultimate rock strength under uniaxial compression and also to apply optimum wedging force on a rock sample and create a wing crack. In order to achieve this optimum wedging force, a set of optimum amounts of related parameters should be used. For this purpose, three different optimizer algorithms (e.g. ACO, PSO and GA) were used and compared to identify the best algorithms and optimum results. Although it is not guaranteed that these resulted amounts of parameters are complete optimum, but they are expected to have proper applications as the local and near-global optimums with satisfactory resulting wedging force to create a wing crack.
Conflicts of Interest
The authors declare no conflicts of interest regarding the publication of this paper.
Nomenclature
a = crack amplitude
f = wedging force
l = crack semi-length
α = crack angle
β(α) = model function depending on crack angle
σ = maximum principal stress
μ = internal friction angle