## Abstract

Appropriate simulation set parameters are the precondition to obtain accurate results; while the simulation results are affected by multiple factors, it is thus crucial to investigate the sensibility of different factors. This paper first analyses the application situation of numerical simulation software in the field of geotechnical engineering and finds that Fast Lagrangian analysis of continua in three dimensions (FLAC3D) has been widely used on roadways or tunnels. Then, taking the roadway excavation process as the engineering background, FLAC3D was used to create 171 schemes of different simulation parameters and analyze the influence of different factors on the simulation results. The findings show that there is a considerable difference in the degree of effect of different parameters on the simulation results. Most of the factors have a remarkable effect on the numerical simulation results (displacement and stress), and only some factors (parameter uniformity and density) have almost no effect on the results. Meanwhile, the trend of displacement and stress is opposite in most cases. In addition, some neglected factors can also have a considerable effect on the simulation results, such as the zone amount; therefore, it is necessary to avoid the variation of nonstudy factors as possible when carrying out the numerical simulation. This study may significantly assist concerned engineers and technicians in developing a more organized and thorough grasp of the impacts of various parameters on simulation outcomes.

## 1. Introduction

The challenges of mining underground mineral resources have grown increasingly difficult and dangerous due to the increasing depth of mining. Numerous researchers have conducted studies to address these challenges that limit the safe and effective production of mines, from the appearance [1, 2] to the essence [3, 4], from the Macro [5, 6] to the Microscopic [7, 8] to the Micro [9, 10] structures, and from the single physical field [11, 12] to multiple physical coupling field [13, 14], and there have been many outcomes. Nevertheless, the complex and variable environment of underground roadways makes it difficult for traditional theoretical analyses [15, 16] to resolve a specific complex engineering problem, and it is laborious and time-consuming to conduct scaled physical simulation tests [17, 18], and it is difficult to reproduce overly complex scenarios, and the accuracy of the obtained results cannot be guaranteed. As the understanding of the properties of geotechnical materials grows and computers continue to develop, computational mechanics [19] is in a flourishing stage. Considering the mutual coupling relationship between various fields, computational mechanics, an emerging interdisciplinary discipline, has significant advantages in processing practice engineering problems. Numerical simulation code based on computational mechanics can simulate approximate object comprehensions for almost any complex operating conditions. At present, the common numerical simulation methods are Finite Element Method [20] (FEM), Finite Difference Method [21] (FDM), Discrete Element Method [22], Boundary Element Method [23], and so on. Among them, due to its precision and speed, the Fast Lagrangian analysis of continua in three dimensions (FLAC3D) [24], a representative of the FDM, has emerged as one of the most popular numerical simulation software.

Researchers and technicians have conducted a large number of simulation studies on roadway/tunnel excavation using FLAC3D. Unlu et al. [25] used FLAC3D to determine the variation rule of the radial boundary displacement along the longitudinal direction of a circular tunnel located in the initial stress field, and based on the deformation behavior of linear elastic materials, the expression for the radial displacement of the tunnel excavation surface was obtained by nonlinear curve fitting. Xiao et al. [26] based on the internal stress field distribution law of coal rock in the process of mine roadway excavation obtained by FLAC3D, combined with the force-electricity coupling relationship equation between the electromagnetic radiation intensity generated in the process of compression and deformation and rupture of the coal rock and the internal stress of the coal rock, and researched the change rule of the electromagnetic emission (EME) signals generated in the process of roadway excavation. Meng et al. [27, 28] used FLAC3D to simulate the rheological characteristics of a deep high-stress soft rock tunnel, analyzed the rheological deformation law of the top plate, bottom plate, and two gangs of the surrounding rock of the soft rock tunnel under the action of high stress, and obtained the depth of burial of the tunnel-time creep curve, lateral pressure coefficient-time creep curve, elasticity modulus-time creep curve, and hysteresis coefficient-time creep curve. Suo et al. [29] applied FLAC3D to analyze the plastic damage of the tunnel when the inner, overlapping, and outer staggered arrangement of the working face of the coal seam tunnel under the very close coal seam group, and the plastic damage of the tunnel, the vertical stress of the top plate, and the sinking displacement of the top plate. Niu et al. [30] established a strength attenuation model for the ruptured perimeter rock of a deep tunnel and implanted the model into FLAC3D to verify the reasonableness of the established model. Pongpanya et al. [31] used FLAC3D to study the damage behavior of the roadway at different depths of overburden and found that factors such as the depth of excavation, the bearing capacity of the support system, and the depth of overburden are closely related to the plastic zone of the roadway. Zuo et al. [32] used FLAC3D to obtain the distribution law of roof stress and displacement under anchor support and equal-strength beam support and found that the use of the concept of equal-strength beam support can significantly optimize the distribution of stress in the roadway so that the local stress in the roof plate shows a uniform state, and the deformation of the surrounding rock is effectively controlled, which verifies the feasibility of the equal-strength beam support technology. Xue et al. [33] proposed a comprehensive mining roadway over-support program using automatic support and anchor combination unit based on FLAC, described the structure and working principle of the support robot, and proposed a method for determining the working resistance of the over-support bracket based on the over-support bracket peripheral rock mechanics coupling model. Dibavar et al. [34] focused on the effect of longitudinal and transverse spacing on the stability of tunnels with “umbrella arch” support and proposed to keep the longitudinal and transverse spacing more than 2.5 times the diameter of the tunnel. Wang et al. [35] compiled a command flow for energy dissipation to realize the secondary development of FLAC3D strain softening constitutive model, which extends the energy calculation module of FLAC3D. Mahmoudi et al. [36] used FLAC3D to model tunnels adjacent to caverns and investigated the effect of location, size, distance, shape, and arrangement of the caverns on the displacements, bending moments, and axial forces in the tunnel lining. Wu et al. [37] used FLAC3D software to conduct stability analysis of five different section-shape tunnel models under the conditions of no support and different base-angle support angles and carried out simulation verification under the actual working conditions of the Sanshandao gold mine.

In conclusion, although FLAC3D has been widely applied, most of the literature focuses on the use of the software to solve a specific problem but seldom focuses on the proper use of the tool itself. For beginners, how to set up various parameters in the simulation (not only geotechnical parameters) to obtain accurate and reliable simulation results has always been a problem for them. Therefore, in this paper, we designed 171 simulation scenarios based on the engineering background of roadway excavation (FLAC3D-based) and comprehensively and systematically explored the influence of twelve factors on the numerical simulation results, which can provide a careful reference for many FLAC3D users to design and adjust the parameters.

### 1.1. Basic Information on FLAC3D

#### 1.1.1. Software Overview

FLAC3D was created by Cundall et al. and is widely applied in geotechnical and mining engineering analysis and design at present [24]. It is a 3D numerical analysis code that was developed based on continuous medium theory and explicit FDM. FLAC3D is especially well suited for dealing with FEM difficult-to-solve complex geotechnical subjects, typically such as complex multiconditions, large deformation, nonlinear material behavior, occurrence, and development of destabilization damage [38]. Therefore, FLAC3D is a well-suited numerical simulation software for underground subterranean works.

The numerical simulation section of the China Geotechnical Forum (aka YanTuBBS) was surveyed, and the more active sections (those with more than 1000 posts) were statistically analyzed (Figure 1). FLAC is the most active section of them all (over 40% of the total). The literature search engines Google Scholar (GS), Web of Science (WOS), China National Knowledge Infrastructure (CNKI), and WanFang Academic Search System (WFASS) were used for keyword retrieval. Where WOS, an English literature search engine, only used the keywords “roadway” and “tunnel,” whereas CNKI and WFASS, Chinese literature search engines, only used the keywords “巷道” and “隧道” (i.e., “roadway” and “tunnel” in Chinese). GS is a literature search engine with mixed, and the outcomes are derived by applying logical operators like “AND/OR” to the search terms. Due to the different operating principles of each search engine, GS obtained the highest number of documents, while WOS obtained the least (Figure 2). The Large general numerical simulation software such as ANSYS, ABAQUS, MIDAS, and so on retrieves the papers in GS totaling more than 200,000. FLAC3D (about 34,500) is larger than FLAC2D (about 5950), so it can be seen that the 3D version is currently being used more frequently than its 2D counterpart in FLAC. Figure 3 demonstrates how FLAC is widely utilized in the field of “roadway” or “tunnel.”

In summary, FLAC3D has been widely used in the field of geotechnical engineering, especially in the field of “roadway” and “tunnel” Therefore, it is of great engineering significance and theoretical value to study the influence of various factors (in FLAC3D) on the stability of the tunneling process and to provide guidance and help for the beginners of many simulation software.

### 1.2. Fundamental Principles

FLAC3D is a finite difference numerical program, which mathematically uses the fast Lagrangian method. Among them, the FDM is a method for finding numerical solutions to definite problems of partial differential equations and systems of equations [39, 40]. The basic idea of FDM is to discretize the problem’s defining domain into a mesh (zone) and then, at the gridpoints, replace the differential quotient in the definite problem with the difference quotient according to the appropriate numerical discretization formulas, thus discretizing the problem into a different format, which in turn leads to a numerical solution. This method is widely used because it is easy to implement on computers [41]. The fast Lagrangian method [42] is a stepwise solution based on explicit differencing to obtain all the equations of motion and constitutive equations of the model, whose constitutive equations are derived from the basic stress-strain definitions and Hooke’s law, while the equations of equilibrium of motion are directly applied to the Cauchy equations of motion (which are derived from Newton’s law of motion). Its computational model is generally composed of several different shapes of three-dimensional units, that is, the dissected spatial unit mesh area, and each unit is further divided into tetrahedra consisting of four nodes in the computation, and the stress-strain of the tetrahedra is only transferred to the other tetrahedra through the four nodes, which is then transferred to other units. When a load is applied to a node, the load acting at that point only affects several surrounding nodes (neighboring nodes) for a tiny period. Using the equations of motion, the relative displacements between the units can be calculated based on the change in velocity and time of the unit nodes, which in turn leads to the unit strains, and then using the constitutive model of the unit, the unit stresses can be calculated. In the process of calculating the strain, the Gaussian integral theory is utilized to simplify the three-dimensional problem by transforming it into a two-dimensional problem. In the equation of motion, the viscosity of the geotechnical body is also fully considered, which is regarded as damping attached to the equation.

### 1.3. Simulation Step Recommends

The “roadway” and “tunnel” (the follow-up is all called: roadway) excavation simulation can be broken down into the following phases, as indicated in Figure 4, depending on the “FLAC3D user help manual” and the author’s practical experience with the software. The simulation process of statics is roughly divided in the “FLAC3D user help manual” into fourteen phases. The “Project Planning and Setup” section contains another seven steps, and the “Tips and Advice” section also gives thirteen tips and advice. Without going into the details of each section, readers can refer to the “FLAC3D help” section of FLAC3D for more specific information.

### 1.4. Building of Model and Simulation Schemes

#### 1.4.1. Numerical Simulation Model

The models are generated by the built-in command flow of FLAC3D, and some special section-shaped roadway models are generated by Rhino. To study the influence of different factors on the stability of the roadway excavation process, a certain initial model and its initial parameters are used as the baseline group (control group). To avoid the interference of irrelevant factors, the following assumptions are made for the model of the benchmark group: (1) the phenomenon of rock mass layering is not considered; (2) rock mass defects, such as fracture zones, joints, and fissures, are ignored, and the model is a homogeneous geological body; and (3) the geostress is applied perpendicularly on the surface of the model. The bottom, left, right, front, and rear faces of the model were fixed, and a compressive stress of 30 MPa was uniformly applied to all the faces except the bottom, and the displacement and plastic zones were zeroed after the model converged. Subsequently, according to the content of the study, to carry out the corresponding excavation simulation, the maximum unbalanced force in the model is 1 × 10 − 5 (the default convergence conditions of the software) to stop the calculation (part of the simulation program under the model cannot be converged, the calculation of 5000 steps).

### 1.5. Simulation Schemes

A total of twelve primary types and 171 groups of simulation test scenarios were created based on the baseline group by changing various factors. Due to the large number of scenarios, this section does not provide a detailed description of each set of scenario parameters. The parameters of the baseline group model and the specific simulation scenarios are shown in Table 1, and the scenario parameters that are not specifically described are the baseline group model parameters for variable control in this study. In addition, the displacement and stress results of the analysis section are the maximum value in the model (or in different directions). The analysis of the results is carried out in terms of displacement and stress.

Type | Factors | Specific scheme (the emphasized part is the baseline group) | |||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

Nongeotechnical parameters | Model size (cube size, m) | 10 | 20 | 30 | 40 | 50 | 60 | 70 | |||||

Zone amount(Mesh gradient coefficient = 1.0) | 52,000 | 320,000 | 512,500 | 972,000 | |||||||||

Mesh gradient coefficient | 0.8 | 0.9 | 1.0 | 1.1 | 1.2 | 1.3 | 1.4 | ||||||

Parameter uniformity(Gaussian distribution, %) | 0 | 1 | 3 | 6 | 10 | ||||||||

Excavation mode | Relax excavation | Assign null | |||||||||||

Roadway area (rectangle, m^{2}) | 4 | 6.25 | 9.0 | 12.25 | 16 | 20.25 | 25.0 | ||||||

Roadway shape | Rectangle | Horseshoe | Three-center arch | Trapezoid | Round | Straight wall semi-circular arch | |||||||

Cyclic footage | 2 | 5 | 10 | 20 | |||||||||

In situ stress (MPa) | 5 | 10 | 15 | 20 | 25 | 30 | 35 | 40 | 45 | 50 | 55 | 60 | |

Lateral stress coefficient | 0.4 | 0.6 | 0.8 | 1 | 1.2 | 1.4 | 1.6 | 1.8 | |||||

Constitutive model | Elastic | Mohr–Coulomb | Mohr–Coulomb strain softening | ||||||||||

Density (kg/m^{3}) | 2200 | 2300 | 2400 | 2500 | 2600 | 2700 | 2800 | 2900 | 3000 | ||||

Geotechnical parameters (M–C) | Elasticity modulus (GPa) | 0.001 | 0.005 | 0.01 | 0.05 | 0.1 | 0.5 | 1 | 1.25 | 1.5 | 1.75 | 2 | 2.25 |

2.5 | 2.75 | 3 | 5 | 10 | 50 | 100 | |||||||

Tensile strength (MPa) | 0.001 | 0.005 | 0.01 | 0.05 | 0.1 | 0.5 | 1 | 1.25 | 1.5 | 1.75 | 2 | 2.25 | |

2.5 | 2.75 | 3 | 5 | 10 | 50 | 100 | |||||||

Dilatancy angle (°) | 0 | 2 | 5 | 10 | 15 | 20 | 25 | 30 | |||||

Cohesion (MPa) | 0.001 | 0.005 | 0.01 | 0.05 | 0.1 | 0.5 | 1 | 1.25 | 1.5 | 1.75 | 2 | 2.25 | |

2.5 | 2.75 | 3 | 5 | 10 | 50 | 100 | |||||||

Internal friction angle (°) | 5 | 10 | 15 | 20 | 25 | 30 | 35 | 40 | 45 | 50 | 55 | 60 | |

Poisson’s ratio | 0.05 | 0.1 | 0.15 | 0.2 | 0.25 | 0.3 | 0.35 | 0.4 | 0.45 | 0.5 |

Type | Factors | Specific scheme (the emphasized part is the baseline group) | |||||||||||
---|---|---|---|---|---|---|---|---|---|---|---|---|---|

Nongeotechnical parameters | Model size (cube size, m) | 10 | 20 | 30 | 40 | 50 | 60 | 70 | |||||

Zone amount(Mesh gradient coefficient = 1.0) | 52,000 | 320,000 | 512,500 | 972,000 | |||||||||

Mesh gradient coefficient | 0.8 | 0.9 | 1.0 | 1.1 | 1.2 | 1.3 | 1.4 | ||||||

Parameter uniformity(Gaussian distribution, %) | 0 | 1 | 3 | 6 | 10 | ||||||||

Excavation mode | Relax excavation | Assign null | |||||||||||

Roadway area (rectangle, m^{2}) | 4 | 6.25 | 9.0 | 12.25 | 16 | 20.25 | 25.0 | ||||||

Roadway shape | Rectangle | Horseshoe | Three-center arch | Trapezoid | Round | Straight wall semi-circular arch | |||||||

Cyclic footage | 2 | 5 | 10 | 20 | |||||||||

In situ stress (MPa) | 5 | 10 | 15 | 20 | 25 | 30 | 35 | 40 | 45 | 50 | 55 | 60 | |

Lateral stress coefficient | 0.4 | 0.6 | 0.8 | 1 | 1.2 | 1.4 | 1.6 | 1.8 | |||||

Constitutive model | Elastic | Mohr–Coulomb | Mohr–Coulomb strain softening | ||||||||||

Density (kg/m^{3}) | 2200 | 2300 | 2400 | 2500 | 2600 | 2700 | 2800 | 2900 | 3000 | ||||

Geotechnical parameters (M–C) | Elasticity modulus (GPa) | 0.001 | 0.005 | 0.01 | 0.05 | 0.1 | 0.5 | 1 | 1.25 | 1.5 | 1.75 | 2 | 2.25 |

2.5 | 2.75 | 3 | 5 | 10 | 50 | 100 | |||||||

Tensile strength (MPa) | 0.001 | 0.005 | 0.01 | 0.05 | 0.1 | 0.5 | 1 | 1.25 | 1.5 | 1.75 | 2 | 2.25 | |

2.5 | 2.75 | 3 | 5 | 10 | 50 | 100 | |||||||

Dilatancy angle (°) | 0 | 2 | 5 | 10 | 15 | 20 | 25 | 30 | |||||

Cohesion (MPa) | 0.001 | 0.005 | 0.01 | 0.05 | 0.1 | 0.5 | 1 | 1.25 | 1.5 | 1.75 | 2 | 2.25 | |

2.5 | 2.75 | 3 | 5 | 10 | 50 | 100 | |||||||

Internal friction angle (°) | 5 | 10 | 15 | 20 | 25 | 30 | 35 | 40 | 45 | 50 | 55 | 60 | |

Poisson’s ratio | 0.05 | 0.1 | 0.15 | 0.2 | 0.25 | 0.3 | 0.35 | 0.4 | 0.45 | 0.5 |

## 2. Analysis of Simulation Results

### 2.1. Model Size

The model size is the overall size (length, width, and height) of the numerical model created, and its value mainly relates to the influence of boundary effects. Figure 5 shows the simulation results in fourteen schemes of model size (with seven model sizes in two conditions of considering gravity and did not consider). The roadway could not converge when the model size is 10 × 10 × 10 m due to the roadway outline being too near to the model boundary. The displacement of the roadway tends to gradually decrease with the increase in model size. On the contrary, the stress in the model tends to gradually increase as the model size increases. Additionally, when gravity is considered, the displacement and stress in the model are slightly higher than gravity, which is not considered.

### 2.2. Zone Amount

In FLAC3D, the model is mainly composed of zones and gridpoints, and the more zones there are, the more time is needed for computation. Figure 6 exhibits the simulation results in four schemes of zone amount. The displacement and stress in the model tend to gradually increase with the increase in zone amount. It should be noted that the vertical stress is higher than the horizontal stress in terms of the magnitude of the increase.

### 2.3. Mesh Gradient Coefficient

A common strategy used in simulations is to encrypt the mesh (zone) attached to the important region to improve computational accuracy. Here, the mesh gradient coefficient is the degree to which the mesh is gradually encrypted from the model boundary to the roadway outline. Figure 7 expresses the simulation results in seven schemes of mesh gradient coefficient. The displacement of the roadway tends to increase as the mesh gradient coefficient increases. The stress in the model is to decrease and then level off as the mesh gradient coefficient increases; the vertical stress is a little higher than the horizontal stress.

### 2.4. Parameter Uniformity

Metallic materials have a uniform distribution of physical and mechanical parameters, whereas rock mass is a nonhomogeneous material. However, the usual simulation treats the rock mass as a homogeneous material and assigns it to uniform parameters. Here, the Gaussian distribution command was utilized to assign nonuniform values to the rock mass material. Figure 8 demonstrates the simulation results in five schemes of Gaussian distribution rate. As the Gaussian distribution rate increases (within 10%), the displacement and stress in the model remain essentially unchanged.

### 2.5. Excavation Mode

There are two main ways to excavate a roadway in FLAC3D: one way is to use the “null” command to “excavate” the roadway, which is a kind of sudden change; the other way is to use the “relax excavate” command, which gradually reduces the physical–mechanical parameters in the region, which is a kind of slow change, which can avoid the interference of inertia effect on the result due to sudden change to a certain extent. In practical engineering, tunnel excavation is a “slow” process, but the “null” command is commonly used for simulation. Figure 9 displays the simulation results in two schemes of excavation mode. When “relax excavate” is used, the displacement in the model can be reduced by about 4.5%, and the stress can be increased by about 7.5% compared with “null.” The excavation mode has a slight effect on the simulation results.

### 2.6. Roadway Area

The cross-sectional area of the roadway is one of the direct factors concerning the stability of excavation. Therefore, many studies on the stability of the roadway directly focus on the roadway span and height (area) as the key considerations. Figure 10 manifests the simulation results in seven schemes of roadway area (rectangle). As the area of the roadway increases, the displacement in the roadway continues to grow in an approximately linear fashion. The stress is decreasing, the rate of stress decreases is also decreasing. In addition, the values of vertical and horizontal stresses converge as the area of the roadway increases.

### 2.7. Roadway Shape

The shape of the roadway is a key factor affecting the state of geostress redistribution after excavation. Here, based on engineering experience, the study was carried out with typical six shapes of roadway sections. In addition, the model built in this section was generated by Rhino software (which also utilized the griddle plug-in to optimize the mesh), and the zone amount and so on were inevitably disturbed, so it was difficult to keep consistent with the baseline group. For this reason, the model size (changed to 50 × 50 × 5.5 m) and the cyclic footage (changed to 5.5 m) were also changed to speed up the simulation process. Figure 11 indicates the simulation results in six schemes of roadway shape. For the displacements, overall, the displacements of the sidewall of the roadway are generally larger than the other positions, while the horizontal displacements are generally slightly larger than the vertical displacements, in addition to the overall displacements: rectangle > trapezoid > straight wall semicircular arch > horseshoe > three-center arch > round. When the stress law is opposite to it, the stress of the roadway shape schemes with large displacement is relatively small.

### 2.8. Cyclic Footage

The cyclic footage is the length of the roadway or tunnel that is excavated at one time during the boring process. Figure 12 shows the simulation results for eight cyclic feed scenarios (four of the cyclic footage lengths under two stopping calculation conditions, the maximum imbalance force reaches 1 × 10^{−5}, or a total of 2000 steps are calculated). The displacements in the model are generally greater when the maximum unbalanced force is used as the convergence condition than when the computational step is used as the convergence condition. This is because the model needs to be calculated after each excavation until equilibrium, so the cumulative steps for the convergence condition of maximum unbalanced force are significantly higher than 2000. Hence, the displacement values are larger for the scheme with maximum unbalanced force as the convergence condition. It should be noted that the stress values in the model are very little affected by this factor. In addition, the effect of cyclic footage on the simulation results is “fluctuating, not simply increasing or decreasing, suggesting that there may be an optimal value for cyclic footage within the simulation range. At the same time, the effect of cyclic footage on the results is relatively slight.

### 2.9. In Situ Stress

The in situ stresses are natural stresses present in the earth’s crust that have not been disturbed by engineering, whose values are influenced by many factors (one of the main ones being the depth of burial) and are at the heart of underground engineering hazards, such as rockburst, collapse, and so on, which are also known as stress-induced hazards. Figure 13 illustrates the simulation results for twelve in situ stress values. The displacement of the roadway increases as the in situ stress increases (as the in situ stress increases, the increase in the rate of displacement is also growing), and the stress in the model also increases. A gentler phase exists at about 30 MPa and before a significant and sustained increase begins at about 40 MPa. In addition, as the in situ stresses exceed 40 MPa, the vertical stresses also begin to be significantly greater than the horizontal stresses.

### 2.10. Lateral Stress Coefficient

The formation of the in situ stress field is influenced by a variety of factors, such as self-gravity stress, tectonic stress, and residual stress, and presents an extremely complex state. Therefore, the values of horizontal and vertical stresses are often different. The lateral stress coefficient is defined as the ratio of horizontal stress to vertical stress. Figure 14 shows the simulation results for eight lateral stress coefficients. The overall displacement of the roadway roughly increases with the growth of the lateral stress coefficient. In particular, when the lateral stress coefficient is small, the displacement of the sidewall of the roadway is larger than that at the floor, and with the increase of the lateral stress coefficient ( >1.0), the displacement at the floor begins to be larger than that at the sidewall, the horizontal stress begins to be larger than that at the vertical stress, and the vertical stress and the maximum shear stress tend to be gradually stabilized.

### 2.11. Constitutive Model

A constitutive model, also known as the mechanical constitutive equation of a material, or the stress-strain model of a material, is a mathematical expression that describes the mechanical properties of a material (stress-strain-strength-time relationship). Currently, the M-C and H-B constitutive models are the most commonly used regarding the simulation of engineering scales (geotechnical). To be able to compare the influence of these constitutive models on the simulation results to a certain extent, we have chosen to analyze only the constitutive models with some of the same parameters (Elastic/E, Mohr–Coulomb/M-C, and Mohr–Coulomb strain softening/MSC). Figure 15 shows the simulation results for three constitutive models. It can be seen that the effect of the constitutive model on the simulation results is prominent. Especially for the MSC constitutive model (using the softening coefficients commonly used in the FLAC3D case), the values of the channel displacements are significantly higher than those of other constitutive models, even though the parameters before softening are consistent with those of the M-C constitutive model.

### 2.12. Geotechnical Parameters

The physical–mechanical parameters of a geotechnical body (hereafter referred to as geotechnical parameters) are one of the most crucial factors to be considered in engineering design. In numerical simulation, different constitutive models require different geotechnical parameters, and different constitutive models can be selected according to different engineering characteristics and engineering requirements. The M-C constitutive model is one of the most commonly applied constitutive models for geotechnical bodies. The required parameters are elasticity modulus, tensile strength, dilatancy angle, cohesion, internal friction angle, and Poisson’s ratio in addition to the density of the surrounding rock is also necessary. Figure 16 displays the simulation results for a total of seven influencing factors and ninety-four sets of simulation scenarios.

For the stress, some factors have essentially no effect on it, such as density, elasticity modulus, tensile strength, and dilatancy angle. The stress has increased gently and slowly with the Poisson’s ratio growth; additionally, the stress fluctuates but has a total rising trend with the cohesion and internal friction angle enlarging. For the displacement, the density also has essentially no effect on it. The growth of most factors (elasticity modulus, tensile strength, cohesion, and internal friction angle) causes the displacement to show a pattern of decreasing and then leveling off. In addition, an increase in dilatancy angle leads to an increase in displacement, while a growth in Poisson’s ratio causes the displacement to fluctuate within a small range.

## 3. Discussion

Numerical simulation has become an indispensable means to solve complex engineering problems, and reasonable setting parameters are the key to guaranteeing the reliability of numerical simulation results. In this paper, we take the tunnel excavation as the engineering background and FLAC3D as the research object to explore the effect of various factors on the simulation results. Although some studies have analyzed the influence of different factors on roadway excavation, for example, Peng et al. [43] studied the influence of the ground stress and the roadway area on the damage zone of the roadway excavation. Zhong-Cheng et al. [44] investigated the influence of different geotechnical parameters on roadway deformation and damage by carrying out orthogonal numerical simulation experiments and concluded that the influence of various factors on roadway deformation is in the following order of significance: cohesion > tensile strength > elasticity modulus > internal friction angle > Poisson’s ratio. In addition, there have been studies on the grid size effect [45-47], but they have not been extended to the engineering scale of the roadway excavation. Most of these studies focus on the effect of one or more factors, which may not provide readers with more comprehensive and intuitive help. To the best of our knowledge, nothing like this paper has been reported. Therefore, the results of this paper may provide researchers or engineers with more comprehensive guidance and reference for numerical simulation.

The factors discussed in this paper are not only limited to geotechnical and construction parameters but also involve some modeling parameters (Figure 17). The results of this study show that the impacts of roadway excavation on rock mass are not only affected by “direct” parameters such as geotechnical and construction parameters but also by modeling parameters such as model size and zone amount. Even the effects of these modeling factors are as strong as those of geotechnical and construction parameters, but they are often “neglected.” For example, when the model size was increased from 20 × 20 × 20 m to 70 × 70 × 70 m, the displacement in the model decreased by 21%, and the stress increased by 38%, while all other parameters were kept constant (Figure 5). When the number of meshes was increased by about three times (from 320,000 to 972,000 considering the self-gravity), the displacement in the model increased by 31%, and the horizontal stress increased by 55% (Figure 6). When the mesh gradient coefficient was increased from 0.8 to 1.4, the displacement in the model doubled by about six times, and the horizontal and vertical stresses were reduced by 50% and 58%, respectively (Figure 7). The effect of such a change on the results is significant. Therefore, when using numerical simulation to study a particular problem, parameter calibration should be performed first, and once all parameters are calibrated, no further changes should be made to parameters other than those of the study objective, or this may lead to significant deviations in the results.

The purpose of conducting numerical simulations is mostly to get regular conclusions (of course, there are some studies to get specific values, at which time it is necessary to try to guarantee that each simulation parameter is consistent with the actual parameter), so some parameters can be simplified appropriately. For example, in reality, the geotechnical body is a nonhomogeneous material, and its parameters cannot be as homogeneous as in common numerical models [48], but assuming its parameters as homogeneous, which has almost no effect on the simulation results (Figure 8). Similarly, the excavation of the roadway can be performed by using commands such as “assign null” and “relax excavation.” Theoretically, “relax excavation” is more in line with the actual excavation process, but “assign null” is more commonly used, and the results show that the effect of these two methods on the results is extremely slight, so the excavation pattern may not interfere with the results more significantly (Figure 9).

In addition, the roadway area (Figure 10), roadway shape (Figure 11), in situ stress (Figure 13), and lateral stress coefficient (Figure 14) all have a more significant effect on the displacements and stresses of the model. The type of constitutive model, on the other hand, determines the various parameters required for the simulation and has an extremely significant effect on the results (Figure 15). At the same time, different geotechnical parameters based on the M-C constitutive model have different effects on the tunneling process and may even lead to a certain arching of the stress (horizontal) distribution characteristics (Figure 18), but this phenomenon has not been reported to the best of our knowledge, perhaps because the cohesion of the rock mass in practical engineering is difficult to reach 10 MPa.

A variety of factors can have a significant impact on the simulation results, and the stress and displacement with the change of the influencing factors show different patterns, but in most cases, the change rule of displacement and stress is opposite (Figures 5,7,9-11,15,16(e), (f) and (g), and in a small number of cases, the change rule of displacement and stress is consistent (Figures 6, 13, 14 and 16(d)), and there are a few cases in which stress is not affected, only displacement changes (Figure 16(b) and (c)). This phenomenon still needs to be further analyzed.

In this paper, a large number of simulation experiments are carried out to visualize the influence of different factors on the simulation results, which may be able to provide a reference for the parameter design and adjustment in numerical simulation research. Generally speaking, the simulation results are affected by many factors; even the zone amount and model size, which are often neglected, will have a great influence on the simulation results. Therefore, when the research is carried out through numerical simulation, it is very necessary to carry out a large number of pretests to calibrate the model parameters to obtain realistic simulation results, and at the same time, it should be pointed out that the calibration is not only about the physical–mechanical parameters of the geotechnical body.

In addition, the design of the schemes in the study is based on the author’s experience, and to further reflect the influence of the parameters on the results, the designed parameter values may exceed the limits that can be observed in practical engineering. Moreover, due to the rich variety of factors studied in the article, an in-depth and detailed analysis of each set of simulation results is difficult to be carried out; thus, the paper only briefly discusses the pattern of the results.

## 4. Conclusions

Aiming at the challenges in the numerical simulation process, a brief survey and analysis of simulation software in the field of geotechnical engineering were carried out first, followed by designing more than 100 sets of research scenarios based on FLAC3D to analyze the influence of various factors on the stability of the roadway excavation with the roadway excavation as the engineering background. The main conclusions are as follows:

Through searching and analyzing the literature, we found that FLAC3D has been widely used in the field of geotechnical engineering, especially in the field of “roadway” or “tunnel.” A brief introduction to the software is given, and the suggested steps to carry out the hydrostatic simulation are given.

The simulation results show that there is a significant difference in the influence of different parameters on the simulation results. The ones with greater influence are model size, zone amount, mesh gradient coefficient, roadway area, roadway shape, lateral stress coefficient, constitutive model, elasticity modulus, dilatancy angle, cohesion, and internal friction angle; a little influence: excavation mode, cyclic footage, tensile strength, and Poisson’s ratio; trifling influence: parameter uniformity and density.

A variety of factors can have a significant effect on the simulation results, and the stress and displacement show different patterns with the changes of the influencing factors, but in most cases, the trend of displacement and stress is the opposite. In addition, even neglected factors (e.g., mesh density, model size, etc.) can have a significant impact on the simulation results. Therefore, when researching through numerical simulation, it is necessary to carry out a large number of pretests to calibrate the model parameters to obtain realistic simulation results, and it should be pointed out that the calibration is not only about the physical–mechanical parameters of the geotechnical body.

## Data Availability

The data that support the findings of this study are available on request from the corresponding author. The data are not publicly available due to privacy restrictions.

## Conflicts of Interest

The author(s) declare(s) that there is no conflict of interest regarding the publication of this paper.

## Acknowledgments

This work was financially supported by the Key Project of Education Department of Hunan Province (22A0293), the General Project of Education Department of Hunan Province (22C0235), and the Postgraduate Scientific Research Innovation Project of Hunan Province (QL20220213).