Numerical manifold method with joint elements for jointed tunnel stability analysis

Abstract

When conducting simulations of jointed rock mass tunnels, the conventional Numerical Manifold Method (NMM) relies on contact algorithms entailing intricate contact detection and open-close iteration processes, thereby leading to restricted computational efficiency. To overcome this constraint, this study puts forward a Joint Element-based Numerical Manifold Method (JE-NMM), which aims to improve the existing situation. The method replaces traditional contact algorithms with joint elements characterized by cohesive constitutive models and incorporates a Newton–Raphson iterative convergence strategy, thus establishing a static analysis framework that is suitable for initially closed joint systems. Results from numerical benchmark tests demonstrate that the new method attains approximately 20-fold higher computational efficiency compared to the classical NMM while preserving a high degree of consistency with the results of physical model test results in stability analysis of jointed tunnel. Findings from engineering application studies show that a locally optimized support scheme, by precisely identifying instability zones, can achieve reinforcement effects comparable to full-face bolting while significantly decreasing the number of bolts used. This approach effectively overcomes the efficiency-related constraints imposed by traditional contact algorithms, providing a novel numerical tool for stability analysis and support optimization in jointed rock tunnels.

Graphical abstract

Keywords

NMM / cohesive zone model / jointed tunnel / stability / localized bolt support

1 Introduction

As a pivotal geological feature controlling the mechanical behavior of tunnel surrounding rock, the spatial distribution characteristics of rock mass joint networks have a decisive influence on the stability of the surrounding rock [1–3]. Current mainstream numerical analysis methods are largely based on continuum mechanics assumptions, incorporating joint effects into complex constitutive relationships through equivalent continuum models [4,5]. While this approach ensures computational efficiency, it poses difficulties in parameter determination, which heavily depends on engineering experience. Discontinuous methods such as the discrete element method (DEM) and discontinuous deformation analysis (DDA) can accurately characterize the discontinuous deformation features of joint surfaces [6–8]. However, their computational efficiency is restricted by the complexity of contact detection algorithms [9], and the mapping relationship between micro-parameters and macro-mechanical responses still lacks a universal calibration theory. Consequently, continuum–discontinuum numerical methods have emerged as a more promising alternative [10,11].

The Numerical Manifold Method (NMM), based on topological and differential manifolds, integrates continuous and discontinuous deformation calculations using dual cover systems [12]. Drawing on the contact theory of DDA for handling discontinuous problems [13], NMM decomposes contacts into vertex–vertex and vertex–edge types and incorporates constraint inequalities into the global equilibrium equations via the penalty function method, thereby enabling the resolution of various complex contact problems [14]. However, the penalty function method requires the introduction of spring stiffness parameters lacking physical significance, and their values significantly impact solution accuracy. As a result, numerous improved contact methods have been developed. Yang and Zheng [15] avoided the penalty function method by directly solving the tangential and normal contact forces for each contact pair, but the algorithm demands substantial memory and is unsuitable for handling multiple contacts. Fan et al. [16] introduced a cone complementary formulation into NMM to compute friction and contact, achieving favorable results. Yan et al. [17] employed virtual thin-layer elements in NMM that transmit only pressure, not tension, to address discontinuous contact problems in tunnel linings. While these improvements have achieved progress, none have fundamentally resolved the trade-off between algorithmic complexity and computational efficiency.

After years of development, NMM has been widely applied to tunnel stability analysis [18,19], slope stability analysis [20,21], rock crack initiation [22–24], hydraulic fracturing [25,26], and seepage-consolidation problems in rock and soil [27–29]. However, owing to the computational inefficiency arising from the numerous complex contacts in jointed rock masses, the application of NMM in jointed tunnels has been restricted.

In tunnel engineering, the contacts needing to be dealt with primarily initially closed rock mass joints or fractures [30,31], and most cases only require calculating small displacements and deformations of the tunnel, making static analysis feasible [32,33]. Therefore, abandoning the classical NMM contact theory can avoid tedious open-close iterations and convergence issues in complex contacts, while introducing Newton–Raphson (N–R) iteration ensures convergence, addressing the limitation of classical NMM static calculations that cannot guarantee convergence [34]. Based on this, this paper incorporates cohesive joint elements into NMM to replace the contact theory and algorithms, thereby forming the Joint Element-based Numerical Manifold Method (JE-NMM), which efficiently solves static tunnel stability problems through iterative convergence. This improvement significantly enhances NMM’s computational speed and is particularly suitable for handling initially closed joint systems frequently encountered in tunnel engineering, providing a new numerical tool for tunnel stability analysis under complex geological conditions.

The main structure of this paper is as follows. Section 2 briefly introduces the fundamentals of NMM, Section 3 presents the application of joint elements in NMM. Section 4 validates the effectiveness of JE-NMM through multiple numerical examples. Section 5 provides conclusions and future prospects.

2 Brief review of numerical manifold method

The NMM employs two cover systems: the mathematical cover and the physical cover. The mathematical cover consists of arbitrarily shaped patches that must completely cover the entire solution domain. These patches are referred to as mathematical patches, and weight functions are defined over them. The physical cover is generated by intersecting the mathematical cover with the computational domain, subdividing mathematical patches into physical patches. When different physical patches overlap (i.e., intersect), they give rise to manifold elements, which are exclusively used solely for integration. This describes the generation process of NMM’s dual cover system.

In Fig. 1, regular triangles are used as the mathematical cover, where the weight functions can simply take the finite element shape functions of triangular elements. Each mathematical patch corresponds to the set of six triangles surrounding any given vertex. In Fig. 1, the mathematical patch is cut by a discontinuity within the computational domain, yielding two physical patches. Due to the presence of the discontinuity, eight manifold elements are formed for integrating the physical patches. Crucially, while physical patches are split by discontinuities, mathematical patches remain intact. This separation between the integration domain (physical patches) and the interpolation domain (mathematical patches) ensures that discontinuities do not compromise computational accuracy. As a result, the NMM can smoothly handle both continuous and discontinuous problems.

NMM will establish a local displacement function for each physical patch in two-dimensional (2D)

(1)u_i(x)={\mathit{p}}^{\text{T}}\mathit{d},

where p is the polynomial basis function matrix and d is the corresponding displacement degree of freedom vector. The approximation function of the 0th-order local displacement is used in this paper

(2){\mathit{p}}^{\text{T}}=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right],

(3)\mathit{d}=\{d_{i1}{d_{i2}\}}^{\text{T}}.

The global displacement approximation function is weighted by the weight function to the local displacement approximation function

(4)u\left(x\right)=\sum _{i=1}^n{w}_i\left(x\right)u_i\left(x\right)=\sum _{i=1}^n{w}_i\left(x\right){\mathit{p}}^{\text{T}}\mathit{d}=\sum _{i=1}^n{N}_i\left(x\right)\mathit{d}\mathbf{=}\mathit{N}\mathit{d},

where w_i(x) is the weight function, n is the total number of physical patches, and N_i(x) is the shape function.

Applying the principle of minimum potential energy yields the time-integrated governing equation of NMM as

(5)\left({\mathit{K}}_{\text{e}}+{\mathit{K}}_{\text{d}}\right)\mathit{D}\mathbf{=}\mathit{F}-\beta {\mathit{F}}_{\text{d}},

where K_e is the global stiffness matrix, K_d is the inertia matrix, D is the displacement degree-of-freedom vector, F is the system force vector, F_d is the inertial force and β is a control parameter. When β = 1, the formulation performs dynamic analysis; when β = 0, it reduces to static analysis. In the static case, the system effectively introduces damping with stiffness K_d, requiring prolonged computation to achieve stabilization.

3 Implementation of joint elements in numerical manifold method

As a partition-of-unity method [35], NMM does not necessitate explicit node and element generation like the finite element method (FEM). Instead, akin to meshless methods, it only requires nodes to be located within mathematical patches (analogous to influence domains in meshless methods). This approach makes joint element construction highly convenient and independent of mesh generation.

During the generation of NMM’s dual cover system, in the presence of discontinuities, two manifold elements will necessarily form on either side. These can arise from discontinuities cutting through physical patches or coinciding with patch boundaries. Thus, after generating the dual cover system, all discontinuities are traversed to identify corresponding manifold element pairs sharing the same edge.

For any manifold element pair and its joint element (Fig. 2), the joint element is derived by virtually offsetting the shared edge along local coordinates son by distances d_n and d_s (They are usually taken as 0.001 to 0.01 times the average manifold element length). These offsets are numerically introduced for convenience and have negligible impact on results when small [36]. The joint element has a thickness of 2d_n and length l, where d_n << l. The nodal displacements of the joint element can be expressed as

(6){\mathit{u}}_{\text{e}}=\mathit{T}{\mathit{d}}_{\text{e}},

(7){\mathit{u}}_{\text{e}}=[{\begin{array}{cccccccc}{u}_{i^{+}}& v_{i^{+}}& u_{i^{-}}& v_{i^{-}}& u_{j^{+}}& v_{j^{+}}& u_{j^{-}}& v_{j^{-}}]\end{array}}^{\text{T}},

(8)\mathit{T}\mathbf{=}\left[\begin{array}{cc}\begin{array}{c}{N}_{i^{+}}\\ N_{i^{-}}\\ 0\\ 0\end{array}& \begin{array}{c}0\\ 0\\ N_{j^{+}}\\ N_{j^{-}}\end{array}\end{array}\right],

where d_e represents the degrees of freedom of physical cover nodes from adjacent manifold elements, and N denotes shape functions.

Given negligible thickness, linear displacement distribution is assumed across the joint element. The horizontal displacement difference between its upper and lower surfaces is

(9)\mathrm{\Delta }\mathit{u}=\left\{\begin{array}{c}\mathrm{\Delta }\mathit{u}\\ \mathrm{\Delta }\mathit{v}\end{array}\right\}={\mathit{B}}^{\prime }{\mathit{u}}_{\text{e}},

(10){\mathit{B}}^{\prime }=\left[\begin{array}{cccccccc}-{\mathit{\text{z}}}_1& 0& -{\mathit{\text{z}}}_2& 0& {\mathit{\text{z}}}_2& 0& {\mathit{\text{z}}}_1& 0\\ 0& -{\mathit{\text{z}}}_1& 0& -{\mathit{\text{z}}}_2& 0& {\mathit{\text{z}}}_2& 0& {\mathit{\text{z}}}_1\end{array}\right],

(11)\{\begin{array}{l}{\mathit{\text{z}}}_1=1-\xi ,\\ {\mathit{\text{z}}}_2=1+\xi ,\end{array}-1\le \xi \le 1,

where ζ typically represents 2D Gaussian integration points.

As joint relative displacements are defined in local coordinates son while Eq. (9) operates in global coordinates xoy, transformation is required via

(12)\mathit{R}\mathbf{=}\left[\begin{array}{cc}\cos \alpha & \sin \alpha \\ -\sin \alpha & \cos \alpha \end{array}\right],

where α is the angle between the global x-axis and local s-axis.

The joint element stiffness matrix is then formulated as

(13){\mathit{K}}_j={\int }_{-\frac{1}{2}}^{\frac{1}{2}}{\left(\mathit{R}{\mathit{B}}^{\prime }\mathit{T}\right)}^{\text{T}}{\mathit{D}}_j\left(\mathit{R}{\mathit{B}}^{\prime }\mathit{T}\right)\text{d}\xi ,

where D_j represents the tangential constitutive relation of the joint element. The classical bilinear traction-separation law [37] (Fig. 3) fails to consider friction under compression and only simplifies the joint behavior as a simple spring constraint. Therefore, to more accurately describe the mechanical behavior between joint surfaces during compressive loading, this study incorporates the Willner friction theory [38] when friction effects manifest on the joint surfaces, addressing the limitation of the original law by explicitly considering frictional interactions under compression.

During the initial elastic phase, the tangent stiffness matrix is

(14){\mathit{D}}_j=\left[\begin{array}{cc}{\displaystyle \frac{G_0}{2d_{\mathrm{n}}}}& 0\\ 0& {\displaystyle \frac{G_0}{2d_{\mathrm{n}}}}\end{array}\right],

where G₀ = E/(2(1 + v)), E₀ = E/(1 − v²), with E and v being the Young’s modulus and Poisson’s ratio of adjacent manifold elements.

Under the bilinear traction-separation law, the damage variable Damage for cohesive elements is computed as follows, noting that frictional behavior under compression does not consider element damage

(15)Damage=\frac{{\delta }_{\text{m}}^{\text{f}}({\delta }_{\text{m}}^{\text{max}}-{\delta }_{\text{m}}^{\text{0}})}{{\delta }_{\text{m}}^{\text{max}}({\delta }_{\text{m}}^{\text{f}}-{\delta }_{\text{m}}^{\text{0}})},

where {\delta }_{\text{m}}^{\text{0}}, {\delta }_{\text{m}}^{\text{f}}, and {\delta }_{\text{m}}^{\text{max}} denote initial, failure, and maximum historical strains, respectively. The stiffness matrix K_j is updated iteratively based on damage evolution.

When the joint element is in a compressed state, the tangential displacement mismatch is w_s calculated and decomposed into recoverable elastic deformation and unrecoverable plastic deformation

(16){\mathit{w}}_{\text{s}}={\mathit{w}}_{\text{s}}^{\text{el}}+{\mathit{w}}_{\text{s}}^{\text{pl}},

In the static friction state, the frictional force is calculated from the micro-sliding displacement, and the tangential frictional force t_s is:

(17){\mathit{t}}_{\text{s}}=k_{\text{s}}({\mathit{w}}_{\text{s}}-{\mathit{w}}_{\text{sa}}^{\text{pl}}),

where k_s is the tangential spring stiffness, and {\mathit{w}}_{\text{sa}}^{\text{pl}} represents the unrecoverable tangential displacement mismatch from the previous time step.

Next, considering the Mohr–Coulomb sliding friction law, we first assume the object is in a static friction state to obtain the trial frictional force {\mathit{t}}_{\text{s}}^{\text{tst}}

(18){\mathit{t}}_{\text{s}}^{\text{tst}}=k_{\text{s}}({\mathit{w}}_{\text{s}}-{\mathit{w}}_{\text{sa}}^{\text{pl}}).

The sliding condition is derived as

(19)F_{\text{s}}=-c-\mu \left|{\mathit{k}}_{\text{n}}{\mathit{w}}_{\text{n}}\right|+\left|{\mathit{t}}_{\text{s}}^{\text{tst}}\right|,

where c is the tangential cohesion, μ is the friction coefficient, w_n is the compressive embedding amount, and k_n is the normal spring stiffness.

When F_s ≤ 0, the joint is in the static friction state, and {\mathit{w}}_{\text{s}}^{\text{pl}}={\mathit{w}}_{\text{sa}}^{\text{pl}} at this time.

When F_s > 0, the joint enters the sliding friction state. Using a non-associated slip criterion, the tangent stiffness matrix D_j^{\mathrm{\prime }} for the joint element material is derived as

(20){\mathit{D}}_j^{\mathrm{\prime }}=\left[\begin{array}{cc}{k}_{\text{n}}& 0\\ 0& k_{\text{s}}\left(1-{\displaystyle \frac{\mathrm{\Delta }\gamma k_{\text{s}}}{\left|t_{\text{s}}^{\text{tst}}\right|}}\right)\end{array}\right],

(21)\mathrm{\Delta }\gamma =\frac{-c-\mu \left|t_{\text{n}}\right|+\left|t_{\text{s}}^{\text{tst}}\right|}{k_{\text{s}}}.

Detailed derivations can be found in Ref. [38].

With joint elements implemented, classical NMM time integration for static analysis becomes unnecessary. This work employs N–R iteration to handle contact nonlinearity, applying incremental loading where the governing equation for JE-NMM at the kth iteration of the nth increment step is

(22)({\mathit{K}}_{\text{e}}+{\mathit{K}}_{{}_j}^{(k)})\mathrm{\Delta }{\mathit{D}}^{(k)}={\mathit{F}}^{\mathrm{n}}-{\mathit{F}}_{\mathrm{i}\mathrm{n}\mathrm{t}}({\mathit{D}}^{(k)}),

(23){\mathit{D}}^{(k+1)}={\mathit{D}}^{(k)}+\mathrm{\Delta }{\mathit{D}}^{(k)},

where F_int represents the internal force vector from joint elements.

4 Numerical examples

This section initially assesses the computational efficiency of JE-NMM through a benchmark example. Subsequently, two numerical examples are conducted to assess the performance and accuracy of cohesive joint elements in the NMM. A layered tunnel case is then simulated and compared with experimental results to demonstrate the engineering applicability of JE-NMM. Finally, a jointed tunnel with randomly generated joint surfaces is analyzed for stability, and local bolt are applied for support to exhibit the practical value of JE-NMM.

4.1 Two contact blocks

A simple two-block contact problem is selected to compare the computational efficiency between classical NMM (C-NMM) and JE-NMM. Both blocks share identical material parameters: elastic modulus E = 30 GPa, Poisson’s ratio ν = 0.25, and density 100 kg/m³. Only gravitational loading is considered, and the contact surface is frictionless. In SNMM, a static solution scheme is employed with 400 analysis steps and a time increment of 0.1. For consistency, JE-NMM uses the same 400-step increment.

It should be noted that the two methods employ fundamentally different solution strategies, and thus adopt distinct convergence criteria reflective of their respective frameworks. C-NMM simulation is considered converged when the maximum displacement increment between consecutive steps remains below 10⁻⁶ m for 5 consecutive steps. In contrast, JE-NMM adopts a true static formulation based on the N–R iterative scheme, with convergence defined by the relative residual norm below 10⁻⁶ of the global equilibrium equation. Although the mathematical forms of the stopping criteria differ, both ensure a comparable level of solution accuracy.

Both methods were run on a laptop with a 13th Gen Intel(R) Core(TM) i5-13500H CPU. The computational efficiency under different mesh sizes is shown in Fig. 4, and the results are presented in Fig. 5. Both methods produce comparable outcomes, while JE-NMM achieves an average 20-times speedup over C-NMM. The computational advantage of JE-NMM becomes more prominent as the element count increases. It is worth noting that the 400-step increment in JE-NMM is solely for incremental load application to ensure convergence. In practical scenarios, 100–200 steps usually suffice, which would double the computational efficiency compared to the results in Fig. 4(b).

4.2 Stacked beams

A cantilever beam composed of two identical stacked beams is modeled. Both beams share material parameters: elastic modulus E = 1500 MPa, Poisson’s ratio ν = 0.25, and are subjected to a concentrated force of 1500 N. Following the setup in Ref. [39], both upward and downward vertical forces are considered. The model dimensions and NMM mesh configuration are illustrated in Fig. 6. The cohesive elements’ maximum traction force ft and fracture energy Gc are set to zero.

The y-direction displacement contour plots from JE-NMM are compared in Fig. 7 (Magnify 1000 ×). The normal gap distance between the two beams under vertical concentrated forces is analyzed in Fig. 8. The results demonstrate high computational accuracy of the proposed method, aligning closely with Ref. [39], thereby validating its effectiveness.

4.3 Square plate with a crack

A square plate with elastic modulus E = 10 GPa, Poisson’s ratio ν = 0.3, and cohesive joint parameters ft = 20 MPa, Gc = 20 J/m² is modeled, as described in Ref. [40]. The model dimensions and JE-NMM mesh are shown in Fig. 9.

This case evaluates the compressive and frictional behavior of joint elements, which are nearly consistent with the results of Ref. [40] as shown in Fig. 10. This indicates that the constitutive law of the compressive-frictional joint elements used in this study can effectively capture the frictional behavior of joint surfaces.

4.4 Layered joint tunnel

Following the experimental setup in Ref. [41], a prototype model of a layered tunnel is established. The bedding dip angle is set to 10°, with rock material parameters: elastic modulus E = 10 GPa, Poisson’s ratio ν = 0.25, and density 2400 kg/m³. A confining pressure of 0.3 MPa is applied around the tunnel, and gravitational loading is considered. The model dimensions and NMM mesh are depicted in Fig. 11. Cohesive elements are assigned ft = 30 MPa and Gc = 10 J/m².

The x-direction and y-direction displacement contour plots of the layered tunnel model (Fig. 12) reveal distinct bedding behavior, with discontinuous displacements concentrated around the tunnel portal. The computed stress contour plot (Fig. 13) shows clear stress concentration zones, particularly at the upper-left corner of the tunnel portal. The failure on the right side of the portal stems from the collapse initiated at the upper-left corner.

A quantitative comparison was conducted by scaling the final displacement data from the physical model test according to the designated similarity ratio and comparing them with the engineering-scale displacement data from our simulation, as shown in Table 1. The results show a good agreement for the left tunnel haunch displacement in the y-direction, with a relative error of only 6.0%. However, discrepancies were observed for the x-direction displacements at the arch haunch and the y-direction displacement on the right side. These differences are primarily attributed to the fact that the experimental data represent the final displacements after failure, whereas our simulation, constrained by the small-deformation assumption, cannot accurately capture the post-failure state. Consequently, the computed displacements, particularly the vertical ones, are smaller than the experimental measurements.

Although JE-NMM focuses on static small-displacement analysis and cannot fully capture progressive failure, the results remain valid for engineering applications. Proactively controlling the initial failure location identified by JE-NMM can effectively mitigate progressive collapse during tunnel support design.

4.5 Cross-jointed tunnel analysis

A randomly generated jointed tunnel model is analyzed using JE-NMM to evaluate stability and optimize support strategies. The tunnel model spans a 100 m × 100 m domain, with a portal radius of 10 m. The model boundaries are fixed on both sides and at the bottom, with a free surface at the top. The rock mass is assumed to be homogeneous and isotropic moderately weathered sedimentary rocks such as sandstone or siltstone. Rock mass parameters include elastic modulus E = 5 GPa, Poisson’s ratio ν = 0.25, and density 2400 kg/m³, with gravity as the sole loading. Cohesive elements are defined with ft = 30 MPa and Gc = 10 J/m². The joint network was generated using the Monte Carlo method, comprising two sets whose dip angles follow normal distributions: the first set has a mean of 70° and a standard deviation of 10°, and the second set has a mean of 170° and a standard deviation of 40°. A tunnel with six intersecting joints is first generated, as shown in Fig. 14. Three monitoring points are established at the tunnel portal, as indicated by 1# and 2# in Fig. 14. The NMM demonstrates low sensitivity to mesh size, enabling accurate results even with coarse discretization.

The stability of the unsupported cross-jointed tunnel is depicted in Fig. 15. Significant vertical displacement occurs above the portal due to two long vertical joints, while displacement elsewhere remains minimal. To tackle this issue, local bolts are installed in the high-displacement zone above the portal. Initially, bolts were arranged within a 90° range, followed by incremental expansions of 90° per stage until full-pattern support (360° ) was achieved. Post-support displacement contour plots (Fig. 16) demonstrate comparable performance among the four schemes, all significantly reducing portal displacement. Notably, the full-pattern bolts exhibited the greatest reduction, followed by the 90° configuration, while other schemes showed slight differences.

Displacement data from three monitoring points around the tunnel portal are presented in Fig. 17. Figure 17 compares vertical displacements at Monitoring Points 1 and 2, revealing a nonlinear relationship between bolt coverage and displacement reduction. The 90° bolt arrangement demonstrated displacement control efficacy comparable to the 360° configuration at the tunnel crown (Monitoring Point 1#), significantly outperforming the 180° and 270° configurations. However, at the left tunnel section (Monitoring Point 2#), the 90° arrangement showed the poorest performance due to the absence of bolt reinforcement in this area. Despite its limitations, the 90° configuration requires only 7 bolts compared to 25 bolts for full-pattern support, resulting in reduced costs and installation time. Consequently, localized bolt support offers superior cost-effectiveness and practicality for this engineering scenario.

5 Comparison with existing approaches

The proposed JE-NMM builds upon the partition-of-unity framework of NMM and introduces cohesive joint elements to model initially closed discontinuities. Its key advantages become evident when compared with existing numerical strategies for jointed rock analysis.

Conventional FEM incorporating cohesive zone models (CZMs) require the mesh to conform precisely to the geometry of pre-existing joints or cracks, i.e., the mesh must be “adapted” to discontinuities. This constraint severely limits their applicability to complex, randomly distributed joint networks, as repeated remeshing becomes inevitable during propagation or reorientation. Although the eXtended Finite Element Method (XFEM) alleviates this limitation by enriching the approximation space to represent discontinuities independently of the mesh, it still faces significant challenges in handling multiple intersecting or branching joints. The enrichment functions in XFEM become increasingly complex and numerically unstable when discontinuities cross or coalesce, and integration over cut elements requires specialized quadrature schemes that hinder scalability. In contrast, JE-NMM inherits the intrinsic mesh-independent nature of the NMM: discontinuities are embedded directly within the dual-cover system without any requirement for mesh alignment or enrichment. Consequently, JE-NMM can naturally and robustly handle arbitrary joint configurations, intersecting, branching, or densely packed, using a fixed background mathematical cover, significantly enhancing modeling flexibility and pre-processing efficiency.

When compared with purely discontinuous methods such as the DDA, the DEM, or the Combined Finite–Discrete Element Method (FDEM), JE-NMM avoids the computationally intensive open–close iteration and explicit contact detection that dominate the cost in those frameworks. FDEM, while capable of capturing fracturing and block interaction, relies on frequent neighbor searches and contact resolution between discrete bodies, leading to high computational overhead in dense joint systems, similar to DEM and DDA. By replacing explicit contact constraints with cohesive joint elements governed by a traction–separation law and solved via a N–R iterative scheme, JE-NMM achieves a static equilibrium solution with guaranteed convergence and dramatically reduced computational cost, approximately 20 times faster than classical NMM (which shares similar contact-handling mechanisms with DDA). This makes JE-NMM particularly suitable for engineering-scale static stability assessments of jointed tunnels, where efficiency, robustness, and ease of setup are prioritized over full dynamic block kinematics.

It is important to note, however, that the current efficiency advantage of JE-NMM is contingent upon the assumption of small deformations and initially closed joints that remain in contact (i.e., no separation or large sliding). If the method were extended to large-deformation scenarios, where joints may open widely, undergo significant sliding, or cause block detachment, explicit contact detection and unilateral constraints would inevitably need to be reintroduced. In such cases, the computational cost would increase substantially, and the efficiency gap between JE-NMM and traditional contact-based methods (e.g., DDA or C-NMM) would narrow. Therefore, while JE-NMM offers a highly efficient solution for quasi-static, small-deformation tunnel stability problems, its applicability to large-displacement or progressive failure analyses remains a subject for future development.

Thus, JE-NMM strikes a favorable balance: it retains the geometric flexibility of meshfree and discontinuous methods, overcomes the multi-joint limitations of XFEM, avoids the contact bottlenecks of DEM/DDA/FDEM, and delivers computational efficiency and convergence stability akin to continuum-based approaches, yet without the mesh-conforming constraint of traditional cohesive FEM.

6 Conclusions

This study incorporates joint elements based on a CZM into the NMM, effectively circumventing tedious contact detection and open-closed iteration convergence processes. The N–R iterative scheme guarantees convenient convergence for static analyses, making the method suitable for stability assessments of tunnels with initially closed joint systems.

The developed JE-NMM achieves a 20-times speedup over classical NMM in static scenarios. Through two benchmark cases, the adopted cohesive model exhibits high accuracy. When applied to jointed tunnel stability analysis, the results of JE-NMM closely match the experimental data. Comparative studies on full-pattern and localized bolt support schemes reveal that targeted reinforcement of critical failure zones with fewer bolts achieves comparable stabilization to full-pattern support, highlighting cost-effectiveness.

It should be emphasized that the current JE-NMM formulation is developed under the assumptions of small deformations, initially closed joints, and purely elastic rock behavior. The method is not applicable to scenarios involving large joint opening, block separation, or post-peak plastic yielding of the rock mass. Future work will extend the framework to elastoplastic constitutive models and large-deformation kinematics.

References

[1] Zhang K S , Wu W , Zhang M , Liu Y , Huang Y , Chen B . New method to identify optimal discontinuity set number of rock tunnel excavation face orientation based on Fisher mixed evaluation. Underground Space, 2024, 17: 300–319

[2] Sun B J , Mei Y C , Li W T , Zhang C , Shao X , Li T , Li W , Zhao W , Wang L . Experimental study on the deformation and failure mode of surrounding rocks of a tunnel in a jointed rock mass. Engineering Failure Analysis, 2024, 163: 108470

[3] Lu S , Sun Z Y , Zhang D L , Liu C , Wang J , Huangfu N . Numerical modelling and field observations on the failure mechanisms of deep tunnels in layered surrounding rock. Engineering Failure Analysis, 2023, 153: 107598

[4] Chang L , Konietzky H . Application of the Mohr–Coulomb yield criterion for rocks with multiple joint sets using Fast Lagrangian Analysis of Continua 2D (FLAC2D) software. Energies, 2018, 11(3): 614

[5] Chen X , Bažant Z P . Microplane damage model for jointed rock masses. International Journal for Numerical and Analytical Methods in Geomechanics, 2014, 38(14): 1431–1452

[6] Zhu H H , Wu W , Chen J Q , Ma G , Liu X , Zhuang X . Integration of three dimensional discontinuous deformation analysis (DDA) with binocular photogrammetry for stability analysis of tunnels in blocky rockmass. Tunnelling and Underground Space Technology, 2016, 51: 30–40

[7] Gao J Y , Peng S Y , Chen G Q , Fan H . An extended discontinuous deformation analysis for simulation of grouting reinforcement in a water-rich fractured rock tunnel. Journal of Rock Mechanics and Geotechnical Engineering, 2025, 17(1): 168–186

[8] Cheng H , Zhao H B , Xie X K . Deformation characteristics and layout optimization of roadway in complex jointed rock mass: A case study based on discrete element method. Computational Particle Mechanics, 2024, 11(4): 1735–1754

[9] Wang X , Wu W , Zhu H H , Lin J S , Zhang H . Acceleration of contact detection between arbitrarily shaped polyhedra based on multi-cover methods in three dimensional discontinuous deformation analysis. International Journal of Rock Mechanics and Mining Sciences, 2020, 132: 104387

[10] Hu Y K , Yan C Z , Jiao Y Y , Wang L , Jia Y , Wang Y . Analysis and countermeasures of asymmetric failure in layered surrounding rock tunnels based on FDEM: A case study. Engineering Failure Analysis, 2025, 167: 109049

[11] Zhang J G , Li Y M , Li Y W , Wang M , Liu Q S , Du C L , Yu H G , Zheng X Q . 3D GPU-accelerated FDEM for fracturing and stability analysis of jointed rock masses due to tunnel excavation. Engineering Failure Analysis, 2025, 174: 109498

[12] Shi G H. Manifold method of material analysis. In: Transactions of the 9th Army Conf. on Applied Mathematics and Computing. Minneapolis, MN: Army Research Office, 1991, 57–76

[13] Shi G H. Discontinuous deformation analysis—A new model for the statics and dynamics of block systems. Dissertation for the Doctoral Degree. Berkeley, CA: University of California, 1988

[14] Shi G H . Contact theory. Science China Technological Sciences, 2015, 58(9): 1450–1496

[15] Yang Y T , Zheng H . Direct approach to treatment of contact in numerical manifold method. International Journal of Geomechanics, 2016, 17(5): E4016012

[16] Fan H , Huang D R , Wang G . Cone complimentary-based numerical manifold method modeling frictional and cohesive contact problems. Applied Mathematical Modelling, 2021, 89(2): 1341–1356

[17] Yan P F , Ren K B , Cai Y C . Numerical modeling techniques for shield tunnel lining structure using the numerical manifold method (NMM). Engineering Analysis with Boundary Elements, 2024, 162: 363–374

[18] Xu X Y , Wu Z J , Liu Q S . An improved numerical manifold method for investigating the mechanism of tunnel supports to prevent large squeezing deformation hazards in deep tunnels. Computers and Geotechnics, 2022, 151: 104941

[19] Xu X Y , Wu Z J , Weng L , Chu Z F , Liu Q S , Wang Z Y . Investigating the impacts of reinforcement range and grouting timing on grouting reinforcement effectiveness for tunnels in fault rupture zones using a numerical manifold method. Engineering Geology, 2024, 330: 107423

[20] Ghaffarpour Jahromi S , Bodaghi F . Slope stability analysis based on the numerical manifold method and the graph theory-case study evaluation. Geotechnical and Geological Engineering, 2020, 38(5): 5523–5534

[21] Zhu S , Yu Z J Z , Tan F , Lv J . Searching slope critical slip surface based on the NMM and equivalent plastic strain. Engineering Analysis with Boundary Elements, 2024, 162: 45–57

[22] Zhang Z S , Wang S H , Wang C G , Wang P . A study on rock mass crack propagation and coalescence simulation based on improved numerical manifold method (NMM). Geomechanics and Geophysics for Geo-energy and Geo-resources, 2021, 7(1): 5

[23] Jiang Y , Wang Y , Cai Z , Zhang Y , Liu Z , Zhang F . Modeling quasi-static crack propagation using preconditioned numerical manifold method. Engineering Analysis with Boundary Elements, 2024, 159: 138–149

[24] Ning Y J , Lu Q , Liu X L . Fracturing failure simulations of rock discs with pre-existing cracks by numerical manifold method. Engineering Analysis with Boundary Elements, 2023, 148: 389–400

[25] Li G , Wang K , Tang C A , Qian X . An NMM-based fluid-solid coupling model for simulating rock hydraulic fracturing process. Engineering Fracture Mechanics, 2020, 235: 107193

[26] Zhang G X , Li X , Li H F . Simulation of hydraulic fracture utilizing numerical manifold method. Science China. Technological Sciences, 2015, 58(9): 1542–1557

[27] Li X L , Zhang H . Analyzing unconfined seepage flow with corner singularity using an enhanced second-order numerical manifold method. Computers and Geotechnics, 2024, 167: 106101

[28] Lin S , Cao X T L , Zheng H , Li Y , Li W . An improved meshless numerical manifold method for simulating complex boundary seepage problems. Computers and Geotechnics, 2023, 155: 105211

[29] Wu W A , Zheng H , Yang Y T . Numerical manifold method for dynamic consolidation of saturated porous media with three-field formulation. International Journal for Numerical Methods in Engineering, 2019, 120(6): 768–802

[30] Wu W , Zhu H H , Lin J S , Zhuang X , Ma G . Tunnel stability assessment by 3D DDA-key block analysis. Tunnelling and Underground Space Technology, 2018, 71: 210–214

[31] Li L , Jiang Q , Liu Q , Fan Q , Zheng Q S , Zhao H R , Wei C . Study on the fracture characteristics of surrounding rock of underground cavern based on CZM-FDEM technique. Engineering Failure Analysis, 2026, 184: 110303

[32] Kawamoto T , Aydan Ö . A review of numerical analysis of tunnels in discontinuous rock masses. International Journal for Numerical and Analytical Methods in Geomechanics, 1999, 23(13): 1377–1391

[33] Li X B , Gong F Q , Tao M , Dong L , Du K , Ma C , Zhou Z , Yin T . Failure mechanism and coupled static-dynamic loading theory in deep hard rock mining: A review. Journal of Rock Mechanics and Geotechnical Engineering, 2017, 9(4): 767–782

[34] Luo S M , Zhang X W , Cai Y C . The variational principle and application of numerical manifold method. Applied Mathematics and Mechanics, 2001, 22(6): 658–663

[35] Lin J S . A mesh-based partition of unity method for discontinuity modeling. Computer Methods in Applied Mechanics and Engineering, 2003, 192(11–12): 1515–1532

[36] Cai Y C , Zhu H H . A locking-free nine-DOF triangular plate element with incompatible approximation. International Journal for Numerical Methods in Engineering, 2017, 109(7): 915–935

[37] Cornec A , Scheider I , Schwalbe K H . On the practical application of the cohesive model. Engineering Fracture Mechanics, 2003, 70(14): 1963–1987

[38] Eberhard PHu B. Advanced Contact Dynamics. Nanjing: Southeast University Press, 2003, 112–115

[39] Sun G , Yi Q , Sun Y , Wang J . The Goodman contact element in geotechnical engineering based on the virtual element method. Archive of Applied Mechanics, 2023, 93(4): 1671–1697

[40] Sun Y , Yi Q , Wang J , Sun G . The virtual element method for rock mass with frictional cracks. Engineering Analysis with Boundary Elements, 2021, 133: 255–268

[41] Li Y H , Zhao W Y , Zhu M G , Liu D Z , Yu H . The deformation and unstability mechanism of surrounding rock of tunnel in the soft and gently dipping layered strata. Journal of Mining and Safety Engineering, 2024, 41(6): 1148–1157

RIGHTS & PERMISSIONS

Higher Education Press