2.1 Basic principle

In general, the homogenization theory works on a representative volume element (RVE) with a finite size, where the material is considered to be homogeneous at the macroscale, and its mechanical properties are unknown and will be directly derived from the heterogeneous microstructures, whose properties are known in advance. At this point, the RVE includes statistically sufficient information on the microstructure. Therefore, the macroscopic mechanical behavior of the heterogeneous material can be equivalent to those of the RVE on average. In principle, the larger RVE size can enhance the accuracy of the effective properties for a disordered microstructure since more microstructures will be counted, which in turn involves a higher computational cost. Therefore, a finite RVE size is usually selected to make a trade-off between accuracy and efficiency. Besides, RVE is supposed to satisfy the scale requirement. Specifically, RVE will be small enough compared with the scale of the macroscopic objective, while being adequately larger than the dimension of microstructures (Elmasry et al., 2022).

Microcracked rock can be considered as a specific example of matrix-inclusion heterogeneity. Randomly disordered and oriented microcracks are present in a given RVE domain with its boundary , as shown in Figure 2. The shape of microcracks is assumed to be ellipsoidal, with a small thickness (Eshelby, 1957).

               

2.2 A brief overview of micromechanics of rock

At present, considerable studies are being conducted to quantify the macroscopic property of rocks from micromechanics. Based on the Eshelby inclusion solution (Eshelby, 1961; Mura, 1987), various homogenization variants mentioned above are available for different geomaterials, for example, rock, soil, and concrete. These models are becoming increasingly robust as more influential factors are taken into consideration. Here, this paper aims to present a brief review of the homogenization scheme instead of the difference between these variants. According to the homogenization scheme, the multiscale model can be divided into single-level homogenization and multi-level homogenization.

In terms of microcracked rocks, homogenization is usually carried out with a single-level form as shown in Figure 2. The initial focus is on the estimation of effective properties in the absence of interaction between random uniformly distributed open microcracks (Eshelby, 1961; Kachanov, 1994; Mura, 1987; Walsh 1965). However, this deviates from reality, where compression is dominant in geoengineering; even in the damaged areas caused by human activity, closure–frictional microcracks account for a large proportion. Later, continuous improvements have been made to take more factors into account. Systematic studies on brittle and quasi-brittle geomaterials have been carried out by Shao and Zhu (Zhu, 2016, 2017; Zhu, He et al., 2016; Zhu, Hu, et al., 2008; Zhu, Kondo, et al., 2008; Zhu, Liu, Wang, et al., 2015; Zhu & Shao, 2015, 2017; Zhu et al., 2018). For example, a damage–friction-coupled model was developed based on the Eshelby inclusion solution with a same-orientated family of microcracks to describe the nonlinearity, unilateral effect, damage and friction, and dilatancy of anisotropic rocks (Zhu, Hu, et al., 2008; Zhu, Kondo, et al., 2008). However, the model fails to predict the material strength. To this end, a refined damage–friction model based on MT homogenization and a thermodynamics framework has been proposed (Zhu, Liu, Wang, et al., 2015; Zhu & Shao, 2015). It enables the model not only to predict the material strength and describe the hardening and softening behavior but also to account for the interaction of microcracks and their spatial distribution. Furthermore, the strength criterion of quasi-brittle rocks has been established by analyzing the damage–friction model based on MT homogenization, in which the three possible failure types (i.e., tension failure, shear–tension failure, and compression–shear failure) are identified and a smooth transition between them is theoretically verified. The analysis has provided the linkage of mesoscale parameters and macroscale counterparts and an analytical solution of the stress–strain relationship (Zhu, 2016, 2017; Zhu & Shao, 2015). On the basis of the plastic deformation attributed to the frictional sliding of microcracks, a plastic–damage-coupled model has been derived based on MT homogenization (Zhu, He et al., 2016). Later, this model has been refined to interpret the effect of a wide range of confining compressive stresses (Zhao et al., 2018). To enable the model to be used in engineering applications, a multiscale time-dependent damage creep model has been put forward to investigate the long-term strength of hard brittle rocks (Zhao, Lai, et al., 2022). In addition, attempts to develop a hydro–mechanical coupling model have already been made to predict the mechanical behavior under undrained triaxial compression tests of sandstone (Zhu et al., 2018). Apart from the prediction of rocks, this model has also been successfully tested in high-performance reactive power cement pastes (Zhang et al., 2019) and ice-saturated frozen silt (Zhao, Lai, et al., 2021). More detailed data on these multiscale models can be found in the previous work (Zhu & Shao, 2017).

Apart from the single-level homogenization of microcracked rocks, a series of multiscale models with multiphase and multilevel homogenization have also been developed to further describe the heterogeneous rocks, for example, shale, sandstone, and concrete. As shown in Figure 3, shale can be regarded as homogeneous at the macroscale level (Level III) and heterogeneous at the millimeter scale (Level II), composed of quartz, calcite, dolomite, and the porous matrix at the lower scale (Level I). Specifically, the interfacial transition zone between inclusions and the matrix is also considered. At Level I, the porous matrix can be further divided into pores and elementary clay particles at Level 0. Homogenization starts from the lower scale to upscale, and upscale homogenization will be inherited from the lower-scale equivalent matrix. At the same level, homogenization will be performed for each inclusion one by one. (Chen et al., 2016). Based on this homogenization scheme, estimations of macroscopic effective properties have also been carried out for shale (Chen et al., 2016; Zhang et al., 2022), sandstone (Li et al., 2021), mudstone (Graham et al., 2022), unsaturated concrete (Yan et al., 2013), saturated concrete (Chen et al., 2019), and fiber-reinforced concrete (Jiang et al., 2019; Zhang, Ju, et al., 2020). However, inclusions are restricted to regular shapes, for example, circle and ellipse in 2D, and sphere and ellipsoid in 3D.

               

3.1 Multiscale continuum method

To begin with, the multiscale continuum refers to a collection of multiscale methods with different continuum-based methods, for example, the multiscale FEM (Ms-FEM) (Zhang et al., 2009; Zhang, Lu, et al., 2015), the multiscale finite volume method (Ms-FVM) (Durlofsky et al., 2007; Jenny et al., 2005, 2006), and the multiscale extended FEM (Ms-XFEM) (Xu, Hajibeygi, et al., 2021). Since these multiscale methods are respectively formulated from their origins (i.e., FEM, FVM, and XFEM), and almost the same upscaling scheme is utilized, Ms-FEM can be taken as an example for illustration.

Ms-FEM was developed by Hou and Wu (1997) on the basis of previous work (Babuška & Osborn, 1983; Babuška et al., 1994); the idea is that the numerical base function is extracted from the solution of the equilibrium of the local unit cell rather than predefinition by the macroscopic constitutive law. As shown in Figure 4, the computational domain is discretized into coarse- and fine-scale mesh, and the fine-scale unit cell with its coarse-scale mesh as the boundary condition will be locally solved to provide the numerical base function for the coarse mesh. Here, the coarse mesh is regarded as homogeneous, while the fine mesh is utilized to describe the mesoscale heterogeneity. In this way, mesoscale heterogeneity can be reflected in the macroscopic base function. In general, the base function can be solved by an analytical or numerical solution. In most cases, the analytical solution is inaccessible or difficult to obtain. In this sense, the numerical solution has prevailed. This multiscale method has advantages in capturing the influence of microscale parameter changes at the macroscale, which can not only save the computational cost but also ensure solution accuracy.

               

Relevant numerical examples mainly focus on the investigation of concrete (Yu et al., 2023) and porous media (Lu et al., 2016, 2017; Zhang et al., 2009; Zhang, Lu, et al., 2015), for example, flow analysis (Dostert et al., 2008), elastic consolidation (Zhang et al., 2009), elasto–plastic consolidation (Zhang, Lu, et al., 2015), multiphase coupling problem (Zhang, Lu, et al., 2015), crack propagation (Lu et al., 2016), and strain localization (Lu et al., 2017). Meanwhile, this multiscale method has also been extended to engineering applications, for example, pavement engineering (Allen et al., 2017; Kim et al., 2013; Sun et al., 2019), dam engineering (Lu et al., 2017, 2022; Zhang, Lu, et al., 2015), and foundation engineering (Lu et al., 2017, 2022). Apart from FEM as the multiscale matrix, other continuum methods can also be utilized as alternatives due to their respective advantages, such as XFEM and FVM. For example, Ms-XFEM is utilized to study the fracturing of heterogeneous media (Xu, Hajibeygi, et al., 2021). Representative examples are listed in Figure 5.

               

3.2 Hierarchical continuum/discrete coupling method

In this section, the continuum/discrete coupling method in a hierarchical form will be reviewed. With this kind of method, the continuum is treated as homogeneous, whose constitutive model can be derived from the discrete model in real time during the whole simulation instead of the prescribed phenomenal macroscale model in advance.

FEM/DEM can be taken as an example. As is shown in Figure 6a, FEM aims to resolve macroscale boundary value problems, while RVE represented by the heterogeneous mesoscale DEM model is embedded into each local Gauss point to interpret mesoscale information. Among them, numerical homogenization is utilized to bridge them together. At the outset, the macroscale strain at the Gauss point calculated by FEM serves as the boundary condition for its corresponding mesoscale DEM. The boundary applied to the RVE will be periodic, including displacement and rotation. This can overcome overestimation of the homogeneous deformation and the zero rotation boundary and underestimation of the uniform force and the free rotation boundary (Meng et al., 2019). In turn, the stress tensor and tangent operator calculated by numerical homogenization will be returned back to the Gauss point for the next iteration.

               

Considerable studies have been carried out on FEM/DEM coupling for sandstone, limestone, soil and rock mixture, and granular materials. Numerical examples mainly cover localization of porous sandstone (Wu et al., 2018, 2019; Wu, Liu, et al., 2020; Wu, Zhao, Liang, 2020a, 2020b), localization of porous limestone (Wu, Papazoglou, et al., 2020), localization of granular materials (Guo & Zhao, 2014; Zhao & Guo, 2014), and engineering applications, such as the retaining wall problem (Guo & Zhao, 2016a), the footing problem (Guo & Zhao, 2016a), and the intrusion problem (Liang et al., 2022). Several numerical examples are listed in Figure 6.

Under the assumption of small strain, FEM has the disadvantage of mesh distortion in the large deformation state. Although there are solutions to tackle this issue, for example, the remeshing technique, reliability, and accuracy cannot be guaranteed. To eliminate the mesh dependence of this coupling method, a potential solution is to replace FEM with the mesh-free method. One available method is the coupling material point method (MPM)/DEM method (Liang & Zhao, 2019; Liu et al., 2017), as shown in Figure 7a. MPM represents a continuum by Lagrangian particles embedded into a fixed background Eulerian mesh. Unlike other mesh-free methods, the field variables of MPM are solved on a fixed mesh with its encapsulated particles. In this coupling scheme, the fixed background mesh makes the coupling method mesh independent and suitable for large deformation. Besides, the Gauss point in FEM is replaced by the material points as carriers of RVE. For engineering applications, it has already been applied for modeling of penetration of a rigid footing into the soil foundation (Liang & Zhao, 2019; Liang, Zhao, Wu, Zhao, 2021), soil–pile interaction (Liang & Zhao, 2019), soil column collapse (Liang & Zhao, 2019; Liu et al., 2017; Zhao, Zhao, et al., 2022), granular silo discharge (Zhao, Chen, et al., 2022), anchor pullout in sand (Liang, Zhao, Wu, Soga, 2021), and slope failure (Wang et al., 2022; Zhao, Chen, et al., 2022).

               

As for the coupling method for large deformation, another alternative to FEM is smoothed particle FEM (SPFEM), which is an extension of particle FEM and smoothed FEM (Guo et al., 2021). In fact, smoothed particle FEM (SPFEM) inherits the advantages of both particle FEM and smoothed FEM. Instead of the Gauss point in FEM, the Delaunay triangular node associated with RVE can be utilized to store relevant field variables in SPFEM. In this way, interpolation can be avoided while remeshing after large deformation. Although MPM enables large deformation, the accuracy of MPM at the boundary is unreliable due to its local nature. In addition, an explicit integration algorithm in MPM requires a small time step to reach equilibrium. Third, it suffers from mesh/discretization sensitivity. Fortunately, these issues have been addressed in SPFEM. Furthermore, the SPFEM/DEM coupling method has been numerically demonstrated by soil footing failure and slope failure (Guo et al., 2021). Specific numerical methods and their corresponding representative examples are shown in Figure 7.

3.3 Concurrent approach

In recent years, a family of concurrent multiscale modeling methods has been developed to make full use of discontinuum- and continuum-based methods in microscale/mesoscale and macroscale modeling, respectively. Compared with the hierarchical approach, the concurrent approach is a strong-coupling form. The distinct feature of this kind of method is that the solution domain is decomposed into local and nonlocal subdomains that run simultaneously. Specifically, the local subdomain is mainly concerned with the local field of interest and can be solved by micromechanics, while the remaining nonlocal subdomain is determined by macro-mechanics. Besides, a two-way information transfer will occur between the micro–macro coupling domains. Due to the feature of this domain decomposition, this kind of method is also known as the partitioned-domain method.

Generally, three domain partition schemes are always utilized as shown in Figure 8, where the continuum and discontinuum domains are denoted by and , respectively. Type I refers to an idealized scenario where there is an interfacial boundary with zero thickness (i.e., ) as the compatible and equilibrium boundary. Type II requires a finite range of the interfacial zone (i.e., ) to reduce/eliminate the ghost force induced at the interface and enhance the accuracy. In these two schemes, subdomains of continuum and discontinuum are separately independent. By contrast, the continuum domain exists throughout the entire domain in Type III, with only a part of the discontinuum domain of high interest overlapping it. For clarity, these domain partition schemes will be subsequently discussed in detail along with the specific multiscale methods.

               

In terms of the type and scale of interest of the discontinuum-based method, the concurrent multiscale modeling methods to be summarized are classified into the following three categories: the atomistic/continuum coupling method, the DEM/continuum coupling method, and the LSM/continuum coupling method. The atomistic method aims at microscale investigation, e.g., QM and MD, and the other two methods have been developed for mesoscale modeling. Here, the continuum is generalized, and the current multiscale method allows the continuum to be FEM (Tadmor et al., 1996), the FDM (Cai et al., 2007), SPH (Liu & Liu, 2003), peridynamics (PD) (Tong & Li, 2016), MPM (Guo & Yang, 2006), NMM (Zhao et al., 2011), the radial point interpolation-based mesh-free method (RPIM) (Gu & Zhang, 2006), and so on. It should be noted that the focus of this work is multiscale modeling. Therefore, coupling methods using different numerical methods at the same scale are not included here.

3.4 Atomistic/continuum coupling method

The purpose of the atomistic/continuum coupling method is to bridge the gap between the atomistic and macroscopic scales and reduce the computational cost. The earliest attempt to explore micro–macro multiscale modeling dates back to the 1970s (Gehlen et al., 1972), where the crack propagation of a bcc iron lattice was simulated by coupling the nonlinearity of the crack tip at the atomistic level with the linear continuum at the far field. Over decades, several multiscale atomistic/continuum coupling techniques have been developed. Although these diverse variants with different names vary in the governing formulation, the coupling boundary at the interfacial region and the treatment of the continuum, the main idea of coupling atomistics (e.g., MD) with the continuum (e.g., FEM), remains unchanged. Therefore, they can be unified and compared in a single framework. A previous study of Miller and Tadmor (2009) provides a thorough examination of similarities and differences for 14 variants.

In terms of the governing formulation, the energy-based method, for example, quasicontinuum (QC), shares the same energy potential in different models to seek equilibrium (Tadmor et al., 1996). The commonly congeneric methods include the coupling of length scales (CLS) (Abraham et al., 1998, 2000; Rudd & Broughton, 2000), the bridging scale method (BSM) (Wagner & Liu, 2003), and the bridging domain method (BDM) (Xiao & Belytschko, 2004). Despite preserving physical energy conservation, a nonphysical ghost force will occur and further lead to spurious wave reflection at the coupling boundary for the dynamic condition, which has already been demonstrated by researchers (Adelman & Doll, 1976; Aubertin et al., 2009; Celep & Bažant, 1983; Doll & Dion, 1976). To eliminate or alleviate this superfluous phenomenon, a direct remedy, called the deadload ghost force correction, is used to compensate the relevant atoms with a negative deadload (Eidel & Stukowski, 2009; Li & Ming, 2014; Ortner & Zhang, 2016; Shenoy et al., 1999). In this way, ghost-force corrected QC (QC-GFC) with a constant deadload proposed by Shenoy et al. (1999) achieves mitigation to a certain degree, but at the cost of loss of energy conservation. Later, the accuracy was improved by applying this idea to a QC variant, known as cluster-based QC (CQC(m)-E) (Eidel & Stukowski, 2009). The force-based method is an alternative to the energy-based method, in which the force equilibrium at the interfacial boundary is a requirement for compatibility. In this way, the ghost force can be efficiently mitigated. The drawbacks lie in the energy loss due to the absence of a well-defined energy potential and the difficulty in achieving equilibrium (Miller & Tadmor, 2009; Tu et al., 2014). The representative studies include the finite element/atomistic method (FEAt) (Kohlhoff et al., 1991), coupled atomistic and discrete dislocations (CADD) (Shilkrot et al., 2002, 2004), and cluster-force QC (CQC[m]-F) (Knap & Ortiz, 2001).

Another key distinction between different methods lies in the treatment of the interfacial region. As shown in Figure 9a, the interfacial region can be subdivided into the handshake region and the padding region . Specifically, the handshake region is a mixture of atomistic and continuum domains, and its field variable will be a weighted function of both atomistic and continuum domains. In this way, it allows a relatively smooth shift of the field variable, rather than an abrupt vibration from one domain to the other. Furthermore, the padding domain is, in theory, a continuum, which will offer boundary conditions to those atoms in and . Due to the nonlocal nature of the atomistics, that is, MD, the thickness of is highly linked to the interaction range of those atoms in , which will provide a full complement to the boundary atoms. The second strategy is to remove the handshake region as illustrated in Figure 9b, for example, QC (Tadmor et al., 1996), in which coincidence of the atom and node near the interface is necessary to provide a two-way boundary path. This requires a downscaling mesh scheme on the continuum side and in turn leads to difficulty in mesh generation. Generally, a fine-to-coarse mesh, rather than a fully refined mesh, is used as it is far from the boundary.

               

To couple atomistics and continuum properly and find a two-way compatible variable field of to each other, two compatible approaches are commonly utilized, that is, strong compatibility and weak compatibility. In the case of strong compatibility, the movement of padding atoms remains consistent with that of finite nodes in a two-way form. Meanwhile, it necessitates the coincidence of at least one-layer finite nodes and padding atoms near the boundary. In this way, it enables higher accuracy, but the complication lies in the refined downscaling mesh generation near the boundary. The representative method covers QC (Tadmor et al., 1996). To circumvent this downside, the weak-form compatible condition was developed, where the displacement boundary condition between different models is imposed by an averaging or penalty function (Gooneie et al., 2017; Miller & Tadmor, 2009). By contrast, the accuracy is relatively low, which can be exemplified by BDM (Xiao & Belytschko, 2004).

The treatment of the continuum including the constitutive model and the numerical method varies in different coupled methods. In some cases, a simple linear elastic FEM is adopted, for example, CLS (Rudd & Broughton, 2000). In QC (Tadmor et al., 1996), the Cauchy–Born-based FEM is the basis of the continuum. Furthermore, FEAt adopts nonlinear and nonlocal elasticity to separately address the continuum and the padding region (Kohlhoff et al., 1991). Besides, plasticity is also considered in the continuum in CADD to describe the lattice dislocation (Shilkrot et al., 2002, 2004). Detailed information can be found in the previous review work (Miller & Tadmor, 2009). Indeed, other numerical methods have also been used to couple with atomistics, which will be discussed later.

As one of the most well-known methods, QC bridges a gap between the atomistic and continuum scale by coupling MD and FEM (Tadmor et al., 1996). As shown in Figure 10, representative atoms (abbreviated as repatoms) are selected from all systematic atoms, where fully resolved atoms are kept in the concerned region, while selectively coarse repatoms are collected in the remaining region. To obtain the total potential energy of the domain, each unselected atom will be subjected to linear interpolation according to its affiliated mesh (e.g., point A in Figure 10b). This greatly reduces the computational cost on the one hand and the advantage of coincidence of the element node and the atom with the same degrees of freedom is achieved on the other hand. Coupling can be completed by minimizing the total potential energy. QC was initially designed for the emergence of dislocations and lattice defects (Tadmor et al., 1996) and has subsequently been extended to cope with crack and fracture problems (Mikeš & Jirásek, 2022; Miller, Ortiz, et al., 1998; Miller, Tadmor, et al., 1998), rate-dependent fracturing (Ghareeb & Elbanna, 2021), nanoindentation (Moslemzadeh et al., 2019; Tadmor et al., 1999; Zhu, Zhao, Shao, 2016), and nano-cutting (Yang et al., 2021). Several representative numerical examples are listed in Figure 11. The latest development of QC can be found at www.qcmethod.com (Tadmor & Miller, 2022).

               

               

Owing to both the governing equation and the disparate features between the different methods in the coupled methods, for example, mismatch of the degrees of freedom or the time step, the ghost wave reflection exists at the interface between different scales, especially for the dynamic problems. The solution to this incompatible problem mainly consists of the usage of the force-based governing equation and the treatment of the interfacial region. As stated before, the force-based method always results in difficulty in seeking equilibrium and causes energy loss. At this point, the focus on the interfacial region is an alternative. For example, in the handshake region, overlap can be introduced to weigh the field variables of both atomistics and the continuum. Pioneering work involved the development of the generalized Langevin method (Adelman & Doll, 1976; Doll & Dion, 1976), in which the compatible condition is the average of the Hamiltonian energy of fine and coarse scales in the handshake region. Later, this conception was introduced into the coupled continuum/molecular/quantum method (Abraham et al., 1998; Broughton et al., 1999), BSM (Wagner & Liu, 2003), BDM (Xiao & Belytschko, 2004), and the Arlequin method (Dhia & Rateau, 2005) to mitigate the spurious wave reflection. However, these methods can only filter out low-frequency waves, but cannot work at the high frequency of dynamics. Yet, this shortcoming has been overcome by the damping technique in an enhanced BDM (Tabarraei et al., 2014).

Apart from the spatial scale, the temporal scale is also a key point to consider. Generally, the atomistic scale is smaller than the continuum scale. If the same time step is adopted in both atomistics and the continuum, it will result in the use of extensive computational resources. Therefore, the multiple-time step algorithm is an alternative, which allows a larger time step in the continuum than those in the atomistics, for example, BSM (Wagner & Liu, 2003) and BDM (Xiao & Belytschko, 2004).

Due to the high dependence of FEM on the mesh, it is time-consuming and complex for mesh generation when dealing with the refined downscaling mesh in strong compatibility, large deformation, and nonlinear problems. To eliminate this dependence, the mesh-free method has been developed and further used to serve as the continuum to couple with atomistics, such as the MD/SPH coupling method (Liu & Liu, 2003), the MD/RPIM-based mesh-free coupling method (Gu & Zhang, 2006), the MD/MPM coupling method (Guo & Yang, 2006), the MD/PD coupling method (Tong & Li, 2016), and other meshless coupling methods (Liang et al., 2006). Representative numerical examples, including wave propagation and dynamic impact, are presented in Figure 12. Although the workpiece dimensions of the above numerical examples are in the tens of micrometers and still short of the macroscopic level, these multiscale methods mathematically achieve the connection between micromechanics and macromechanics, which is beneficial in terms of saving time in terms of extensive computations (Ghareeb & Elbanna, 2020).

               

3.5 DEM/continuum coupling method

Currently, the atomistic/continuum coupling method is preferred for solid mechanics of metals, nanomaterials, or composites. Due to the complexity of rocks and the limitation of computational resources, it is difficult and unfeasible to study the mechanical behavior of rocks at the microscale. To make a trade-off between accuracy and efficiency, grain-scale modeling is a viable alternative. Similar to the idea of atomistic/continuum coupling, only the anticipated domain is modeled by the mesoscale model, while the remainder is sufficiently modeled by the continuum model based on the macroscale law. In this way, various multiscale DEM/continuum coupling techniques have been devised for rock mechanics. Due to the intrinsic difference between DEM and atomistic methods, for example, the scale, the computational scheme, and the particle/atom arrangement, the concept of the coupled atomistic/continuum method is unsuitable for direct extension to DEM/continuum coupling (Tu et al., 2014). Specifically, both the temporal and spatial scales of the DEM are larger than those of the atomistic method. For example, the DEM particle is of finite size, contrary to the material point of the atom. Moreover, the atomistic method is nonlocal, while the DEM is local. Besides, irregular particle packing is always utilized in DEM, as opposed to the regular arrangement in the atomistic method.

In this section, mesoscale modeling mainly focuses on particle-based DEM while varying in the continuum numerical method. The idea of this coupling method was developed in the last century (Lorig & Brady, 1984; Munjiza et al., 1995). After decades of development, this kind of method has already been used as an effective tool to explore rock mechanics problems. Like the atomistic/continuum coupling method, the two-way information transition between different models is also a key issue in the DEM/continuum coupling method. To address this compatibility problem, there are three common categories, which separately correspond to schemes in Figure 8.

The first one is the contact model. It is also referred to as the master–slave coupling method (Rabczuk et al., 2006), the edge-to-edge coupling method (Tu et al., 2020), the triangular wall coupling method (Bai et al., 2022), and the surface coupling model (Cheng et al., 2023). The concept behind coupling is that the various models directly interact with each other at the interface, where the equivalent variables (i.e., force and torque) are transmitted between the contact point of the DEM and the finite element, for example, triangular mesh. As an illustration, Figure 13a shows the contact model. The key lies in the determination of the contact point, the penetration depth, and the contact normal. Furthermore, if the finite element is regarded as the slave contact, these field variables from DEM can be equivalently projected to the element node by the shape function, and vice versa (Dang & Meguid, 2013; Oñate & Rojek, 2004). Meanwhile, the damping scheme at the interface is also utilized to stabilize the computational system. In this way, the finer element size at the interface is not required, which thus enhances computational efficiency. This coupling scheme is always utilized in cases of rock/soil–structure interaction, for example, rock cutting and impact of the projectile against the rock (Oñate & Rojek, 2004), pile–soil interaction (Elmekati & Shamy, 2010), and soil–tunnel lining interaction (Dang & Meguid, 2013) by using the DEM/FEM coupling method. Besides, the well-known computational framework developed by Itasca provides an integration of its their own methods PFC and FLAC, also known as DEM and FDM. Numerical examples cover underground excavation (Cai et al., 2007; Qu et al., 2019; Yu & Chen, 2021), soft soil reinforced by stone columns (Indraratna et al., 2015; Xu, Zhang, et al., 2021), and dynamic compaction in granular soils (Jia et al., 2018). Some of these numerical examples are illustrated in Figure 13. However, this scheme will contribute to the occurrence of spurious wave reflection at the interface, especially for the high-frequency wave in dynamics. This will be alleviated by the volume coupling method, which will be discussed next.

               

The second method is the overlapping domain method (Bai et al., 2022), also known as the volume coupling method (Cheng et al., 2023). This method is an extension of the Arlequin algorithm (Dhia & Rateau, 2005) and the BDM (Xiao & Belytschko, 2004). This method allows overlap between the DEM domain and the continuum domain, as shown in Figure 14a. In this overlapping domain, a smooth two-way information transition between different models will be achieved by the weighted work contributed by both DEM and the continuum, which is typically done in a linear way. Besides, compatible coupling necessitates additional kinematic restrictions imposed by either Lagrange multipliers or the penalty function method (Li et al., 2015; Rojek & Oñate, 2008; Tu et al., 2014). In this manner, wave propagation problems and the false wave reflection at the interfacial region can be suitably minimized. With various weighting functions, a generalized overlapping domain method is proposed, in which four common coupling methods are compared (Tu et al., 2014, 2020). Among them, it has been demonstrated that the separate domain coupling method is the best choice for allowing the high-frequency wave to travel. In the above methods, either the displacement compatibility or the force interaction condition is adopted at the interface, which will cause stress discontinuity. This has been overcome by a mixed displacement and stress compatibility condition (Tu et al., 2020). Another issue of coupling is the alignment between micro–macro material properties. A common solution is based on the homogenization of RVE (Cheng et al., 2023; Rojek & Oñate, 2008). This coupling scheme has already been successfully applied to investigations of rock cutting (Rojek & Oñate, 2008), the impact on concrete (Zhou et al., 2022) and wave propagation (Tu et al., 2014) by DEM/FEM coupling, and tunnel excavation (Bai et al., 2022) by DEM/FDM coupling; selective numerical examples are presented in Figure 14.

               

The third method is a specific case of the overlapping domain, which is an extension of the BSM to the DEM/continuum coupling method. Different from the two above-mentioned coupling methods, in which an apparent separation between DEM and the continuum exists, the DEM-based BSM decomposes the whole computational domain into a coarse FEM mesh region covering the whole domain and a fine particle region. In this way, the whole DEM region is overlapped with the coarse FEM region, as shown in Figure 15a. The local displacement of particles consists of coarse- and fine-scale components. With the help of the projection operator, it allows the coarse- and fine-scale equations to work independently, which facilitates a multiple-time step scheme. For the dynamic problem, the nonreflecting boundary condition and introduction of virtual particles along the interface are necessary to allow high-frequency wave travel between the FEM and DEM domains. Numerical analyses for slope stability (Li & Wan, 2011), foundation pit (Li & Wan, 2011), hydro-mechanical behavior of granular materials (Wan & Li, 2015), and granular material fracturing (Wang & Li, 2015) have been carried out, as shown in Figure 15b. The main contributions to this coupling method come from the research group of Li (Li & Wan, 2011; Wan & Li, 2015; Wang & Li, 2015).

               

Like the atomistic/continuum coupling method, the DEM-based BSM will result in a larger time step in the continuum than that in the DEM because both the DEM and the continuum adopt the explicit time integration and the element size of the continuum is generally larger than the particle size in DEM. To circumvent the expensive computational cost in the case of adopting the same time step in both DEM and the continuum, the multiple-time step algorithm is a common solution, which allows different time steps in different models (Elmekati & Shamy, 2010; Tu et al., 2014).

3.6 LSM/continuum coupling method

Although LSM, similar to a coarse MD, has been refined over decades, it is still limited to small-scale simulations due to the high computational cost. One of the most well-known variants is DLSM, which realizes the micro–macro linkage in the parameter calibration (Zhao, 2010; Zhao et al., 2011). To further make the DLSM fully multiscale, it is coupled with NMM, which is a macroscopic model that transforms FEM and DDA into a unified framework. As shown in Figure 16a, the Particle Manifold Method (PMM) as a transition zone is developed to achieve the interaction between DLSM and NMM. The PMM element is a combination of the DLSM particle and the NMM manifold element, where PMM interacts with DLSM through particle contact and NMM via the manifold element. In this manner, PMM acts as a cushion to absorb the ghost wave reflection between the micro–macro interface. Another distinct feature is that the macroscopic PMM element can be automatically degraded into microscopic DLSM particles to address rock failure problems once it satisfies the user-defined criterion. Besides, it is noted that DLSM and NMM work simultaneously at the same small time step calculated by DLSM to guarantee numerical stability in the multiscale DLSM. Based on this multiscale framework, the dynamic collapse of a tunnel under rock blasting was numerically investigated, as listed in Figure 16b.

               

As another variant of the LSM, the four-dimensional LSM (4D-LSM) is also extended to a multiscale version by coupling NMM (Jiang et al., 2021), which also utilizes the same computational strategy as the multiscale DLSM. The enriched multiscale 4D-LSM has been numerically tested by large-scale rock-cutting tests.

3.7 Other concurrent methods

Given the natural divergence in the understanding of multiscale modeling among different researchers, the definition of multiscale methods in the academic community is extensive. In this sense, a few numerical methods are still ambiguous, aside from the aforementioned concurrent methods with clear and uncontroversial classifications. One is the mesh-free/continuum coupling method, such as SPH/FEM (Chuzel-Marmot et al., 2011), the element-free Galerkin (EFG)/continuum coupling method (Rabczuk et al., 2006), PD/FEM (Wang et al., 2019), and the PD/boundary element method (BEM) (Chen & Yu, 2022). Although these above-mentioned mesh-free methods have local characteristics, their governing equations or constitutive models are still based on the macroscopic scale. The purpose of these coupling methods is to leverage the mesh-free method to compensate for the shortcomings of the continuum-based method in coping with discontinuous problems. The other method is called the mesh refinement technique, where fine and coarse meshes coexist in the computational domain (Figure 17). The fine mesh is utilized in the position of interest, while the remainder is roughly described by the coarse mesh. Numerical examples cover FEM (Yu et al., 2013), extended FEM (XFEM) (Zhou & Yang, 2012), LSM (Jiang et al., 2021), and RFPA (Li et al., 2022). However, both of the aforementioned methods adopt the single-scale constitutive model with the same numerical framework, which makes it necessary to further discuss the rationality of this technique for multiscale modeling.

               

In the present work, the multiscale model necessitates a lower-scale local property to reveal the upscale mechanical responses or the influence of upscale loading on the lower-scale responses in turn. As a result, different subdomains should contain constitutive models of different scales. Therefore, these coupling numerical methods are not covered in detail in this work.

To further help readers quickly understand the focus of this study, the multiscale numerical methods used here are summarized in Table 2.

4.1 Large deformation of RVE

The accuracy of RVE is not relevant to the hierarchical coupling method. However, the shape of RVE will be severely distorted while dealing with large deformation problems, such as excessive shear, compression, or tension, as shown in Figure 18. At this point, RVE will lose fidelity and be unrepresentative since fewer number of particles remain in the narrow direction. To overcome this shortcoming, reinitialization is a drastic solution. However, the historical deformation will be reset, leading to an unacceptable error in modeling of history-dependent materials (Wang & Zhang, 2021). Another existing solution is use of the adaptive RVE scheme, where the deformed part will be mapped to its image counterpart by taking full advantage of periodic conditions. By doing so, the cuboid RVE shape remains unchanged and the historical deformation is also preserved (Qu et al., 2021; Wang & Zhang, 2021). However, only the case of large compression deformation is validated. In addition, a macroscopic constitutive model is also utilized to replace the distorted RVE by Wang et al. (2022), but additional parameter calibration is required. In this sense, ways to improve RVE at large deformation still need to be determined.

               

4.2 Mesh/element size dependence

The mesh-based method inherently has the disadvantage of mesh size dependence, and this is also true of the coupling methods. In the hierarchical methods, different mesh sizes have diverse strain localization, as shown in Figure 19, especially the width of the shear band, and a finer mesh will lead to a narrow shear band (Guo & Zhao, 2014; Lu et al., 2022). Since RVE is affiliated with the Gauss point or the material point, a significant influence of the element type on the numerical results can also be observed. Coarse mesh with a high-order element type can yield an obviously narrow concentration. A fine mesh or element is beneficial for improving accuracy but at a high computational cost. In this way, particular attention needs to be paid to balancing accuracy and efficiency, depending on the specific problem. This issue also exists in the MPM/DEM coupling method due to the presence of the background mesh in MPM (Liang & Zhao, 2019). Currently, several numerical techniques, including adaptive Ms-FEM (Lu et al., 2022), hierarchical coupled SPFEM/DEM method (Guo et al., 2021), and the second-order regularized FEM/DEM method (Desrues et al., 2019), are proposed to alleviate this drawback. Still, more validation and an advanced technique to overcome this issue are needed. Similar to the case of the hierarchical coupling method, the influence of mesh size at the transitional interface between different models has also been reported in the concurrent coupling method (Li et al., 2015). However, it is not as severe as that in the hierarchical coupling method.

               

4.3 Multiphase and multifield coupling

Rock mechanics problems generally involve multiphase (e.g., solid, fluid, gas, and critical state) and multifield (e.g., mechanics, hydraulics, thermal, chemical) issues. Various phenomenological macroscale constitutive models have been developed, but they are incapable of revealing the microscale/mesoscale mechanism. In recent years, several thermo-mechanical (Zhao, Zhao, et al., 2022) and hydro-mechanical (Guo & Zhao, 2016b) models have been embedded into multiscale methods. However, more influential factors (i.e., multiphase and multifield) should be taken into account in the future.

4.4 Full-scale spatiotemporal modeling

Currently, multiscale modeling in engineering applications mainly focuses on the coupling between the mesoscale and the macroscale due to the limitation of computer power, which lacks the description of microstructures. There are few studies on full-scale modeling from micro to macro scales in engineering applications. To overcome this issue, an advanced homogenization theory in the hierarchical multiscale coupling method or more mature concurrent multiscale techniques linking atomistics, mesoscale discrete modeling, and continuum at the nanometer to millimeter scale and even meter/kilometer scale should be developed. Apart from the multiscale feature in space, multiple temporal scales should also be considered. Generally, the creep constitutive model has always been adopted to carry out long-term analysis (e.g., several years or even longer) in both continuum- and discontinuum-based numerical methods (Zhang & Zhao, 2023; Jin et al., 2024). Due to the presence of damping, it is hard to realize the connection between discrete methods using time integration (e.g., Newton's Second Law) and physical time. In addition, few studies have included full-scale spatiotemporal modeling, which may become one of the key areas of research in the future.

4.5 High-performance computing

Generally, the bottleneck of numerical modeling in engineering applications, on the one hand, is computer power, which is beyond the user's control. To the authors' knowledge, the maximal computational model of the underground blasting problem with more than one billion 4D-LSM particles has been successfully tested at the TianHe-3 supercomputer as shown in Figure 20 (Fu & Zhao, 2022). In the model, the particle size is 0.135 m, but it is still too large to reveal the mesoscale mechanism. Although multiscale modeling reduces the computational cost to a certain degree, on the other hand, high-performance computing is still required to provide technical support for multiscale modeling in engineering applications. Currently, various numerical methods have already enabled CPU- (i.e., OpenMP and MPI) and GPU-based parallel computing. Different from the traditional numerical method, the multiscale numerical method involves different numerical techniques, making parallelization difficult. At present, attempts are being made to realize parallel computing for multiscale numerical methods (Desrues et al., 2019; Guo & Zhao, 2016b; Zhao, Zhao, et al., 2021), where parallelization works only at the RVE level. In addition, it is only tested on a personal computer, and thus, full-scale engineering applications are not yet possible. At this point, researchers should pay attention to taking advantage of multiscale numerical modeling and high-performance computing (e.g., supercomputers) to achieve full-scale engineering applications in the case of limited hardware conditions.