HOSS: an implementation of the combined finite-discrete element method

Nearly thirty years since its inception, the combined finite-discrete element method (FDEM) has made remarkable strides in becoming a mainstream analysis tool within the field of Computational Mechanics. FDEM was developed to effectively “bridge the gap” between two disparate Computational Mechanics approaches known as the finite and discrete element methods. At Los Alamos National Laboratory (LANL) researchers developed the Hybrid Optimization Software Suite (HOSS) as a hybrid multi-physics platform, based on FDEM, for the simulation of solid material behavior complemented with the latest technological enhancements for full fluid–solid interaction. In HOSS, several newly developed FDEM algorithms have been implemented that yield more accurate material deformation formulations, inter-particle interaction solvers, and fracture and fragmentation solutions. In addition, an explicit computational fluid dynamics solver and a novel fluid–solid interaction algorithms have been fully integrated (as opposed to coupled) into the HOSS’ solid mechanical solver, allowing for the study of an even wider range of problems. Advancements such as this are leading HOSS to become a tool of choice for multi-physics problems. HOSS has been successfully applied by a myriad of researchers for analysis in rock mechanics, oil and gas industries, engineering application (structural, mechanical and biomedical engineering), mining, blast loading, high velocity impact, as well as seismic and acoustic analysis. This paper intends to summarize the latest development and application efforts for HOSS.


Introduction
The combined finite-discrete element method (FDEM) was first conceived in 1989 by Munjiza while working at Tohoku University in Sendai, Japan. Initially FDEM was implemented only in 2D, with the first software platform being subsequently developed in 1990s at the University of Wales, Swansea and at the Massachusetts Institute of Technology (MIT). At the time, the FDEM implementation was written in C ++ and the code was called RG [1][2][3]. The code was initially intended to model fracture and fragmentation of cementitious materials and therefore was designed to handle finite displacements and finite rotations while incorporating large-strain-based deformability models.
Shortly thereafter RG was followed by the first commercial FDEM package called ELFEN, which was developed in 1994 for Rockfield Software Ltd [4]. From 1992 onwards proprietary FDEM developments were also done at MIT at their Intelligent Engineering Systems Lab by different researchers and Ph.D. students under the supervision of Professor Williams [5].
Starting in the late 1990s FDEM implementations were also introduced at Lawrence Livermore National Laboratory [6]. Munjiza, meanwhile, developed the open source Y-FDEM software package, which included both 2D and 3D FDEM [7][8][9][10][11][12][13][14]. In the early 2010s this code was extended to 2.5D, i.e., shells, by Lei, who developed the world's first 1 3 FDEM solution for fracture and fragmentation of layered shells and plates [15,16]. In the meantime, Lei also incorporated further enhancements to Y-FDEM by introducing a 3D crack model which was one of the first 3D fracture and fragmentation solutions in FDEM [15,17].
During this timeframe, FDEM developments also began at Imperial College London [18]. Initially the team effort turned Y-FDEM into what was a separate FDEM software platform called VGeST (Virtual Geoscience Simulation Tools) and more recently named Solidity. The Imperial College team has since coupled their FDEM with an inhouse fluid code called Fluidity [19].
Other efforts were also undertaken at Queen Mary University of London to simulate problems that were dominated by fluid [20,21]. The result was the CgLes-Y code (Conjugate Gradient Large Eddie Simulation) software platform, and it has now been made an open source code. The team members have continued with further developments at Tianjin University [22,23]. FDEM has also captured the attention of the Chinese Academy of Science, and a major effort has been put into further development of FDEM simulation tools [24,25]. Work at Wuhan University started with the open source Y-FDEM code by improving its capabilities in terms of contact, fracture, thermal effects and fluids [26,27].
In parallel, a group at the University of Toronto used the Y-FDEM code as a starting point and developed it for commercial applications via a code named Y-geo [28]. Their efforts continued as they adopted Y-geo to GPU (Graphics Processing Unit) achieving an initial speed-up of more than 30 times [29]. This was a pioneering work in porting a FDEM software tool to GPU. Recently, ambitious FDEM developments for GPU were done jointly in Japan and Australia where speed-ups of about 300 times have been reported [30].
In 2003 a FDEM group was established at Los Alamos National Laboratory. The initial aim of this group was to add sophisticated FDEM contact libraries into one of the Laboratory's proprietary multi-physics codes. It was quickly realized that a stand-alone FDEM platform was needed as communication requirements for the contact search and detection needed an independent and robust database. This new FDEM platform was developed using state-of-the-art software design approaches coupled with extensive verification and validation efforts. Initially the platform was called MUNROU (MUNjiza-ROUgier) and development efforts led to sophisticated material and fracture libraries as well as an efficient parallelization scheme. Over a decade later, a custom-developed 2D fluid solver was implemented for hydro-fracture analysis, where the fluid interacts within a contained solid domain and, instead of coupling, the solver was naturally integrated, i.e., the fluid and the solid domains were simulated using the same spatial discretization or mesh [31][32][33]. At this point the platform was renamed as the Hybrid Optimization Software Suite, a.k.a. HOSS [34,35]. Laboratory mission requirements demanded a more robust fluid approach and recent efforts have led to the incorporation of the Fluid-Solid Interaction Solver (FSIS) [36]. Given the team's experience with the earlier approach, an explicit computational fluid dynamics solver has now been fully integrated into the HOSS mechanical solver.
In terms of algorithms, robustness, capability, parallelization, etc., HOSS has proven itself to be one of the leading FDEM software platforms and, as such, it is used by researchers at LANL, commercial entities and designated universities. The results obtained by HOSS have been independently validated against numerous experiments, some of which are shown in this paper. Current LANL application spaces range from hydraulic fracturing, to hypervelocity impact (bolides), earthquake rupture, fully coupled THM (Thermal-Hydraulic-Mechanical) reservoir processes, just to name a few, while other more specialized applications such as underground test containment, nuclear weapons effects (cratering, pyroclastic flow, etc.), high explosive performance, and weapons penetration are also considered.
The rest of the paper is dedicated to highlighting the latest FDEM enhancements incorporated into HOSS, many of which are available to the community at large, together with succinct summaries of some of the most relevant application examples. It should be noted that the theoretical basis of HOSS has been partially disclosed in numerous papers, three patents, and in three books [14,37,38].

Overview of HOSS capabilities
HOSS is an original parallel 2D/3D FDEM simulation platform designed with the main purpose of providing applied scientists or engineers with the ability of utilizing Computational Mechanics procedures for the analysis of complex continuum and/or discontinuum physical systems [39][40][41].
The key advantage of FDEM is the utilization of finite displacements, finite rotations, and finite strain-based deformability combined with suitable constitutive material laws which are then merged with discrete element based transient dynamics, contact detection, contact interaction solutions, and objective discrete crack initiation and crack propagation solutions [14]. At LANL, over the course of the past few years, extensive efforts have been undertaken to improve the accuracy and efficiency of some of the aforementioned portions of FDEM. Examples of this are the innovative solutions implemented in the FDEM solid solver, such as non-locking finite elements [42,43], smooth contact algorithms [44] and a unified cohesive zone model [45], and in the FDEM-based fluid-solid interaction solver [36].

Non-locking finite element formulation
When using lower order 2D or 3D finite element formulations, such as constant strain triangles (CSTRI) and tetrahedrons (CSTET), an issue that often appears is that these types of elements will numerically lock. This is due to their full integration scheme and is often recognized by a checkerboard pattern in the stress distributions. To alleviate this deficiency composite triangular (COMPTRI) and tetrahedral (COMPTET) finite elements were developed. Both COMPTRI and COMPTET formulations use selective integration of stresses in order to avoid this artificial stiff response or locking [42,43]. The developed finite elements utilize a unified hypo-/hyper-elastic approach that allows linkage to user-defined (isotropic or anisotropic) material models.
For the purpose of calculating the kinematics of deformation both of these formulations are based on Munjiza's multiplicative decomposition approach [38,42,43]. Furthermore, the FDEM model used in HOSS was designed in a way that the anisotropic properties of the solid can be specified in a cell-by-cell basis, therefore enabling the user to "seed" anisotropic properties that follow any type of spatial variation. It has been demonstrated that these composite finite element formulations do not exhibit problems related to volumetric locking as is often encountered when using CSTRI or CSTET finite elements.
A demonstration of HOSS' generalized anisotropic deformation kinematics is shown in Fig. 1. Figure 1a presents the setup of the problem under consideration [43], where a cross section of a mountain with four different material layers is shown. A high strain rate source is located inside of material layer #2. Each material layer has different sets of elastic properties. Two different cases were run: (1) isotropic elastic properties and (2) anisotropic elastic properties. The results shown in Fig. 1d and e show two snapshots of the particle speed for the isotropic case, while Fig. 1b and c shows the corresponding snapshots for the anisotropic case. Fracture patterns obtained are also indicated in the figure.
The results clearly demonstrate the effect of the material anisotropy on both the wave and fracture propagation. In the anisotropic case the shape of the wave propagation front is strongly affected by a combination of the material axes orientation and the anisotropic nature of the elastic properties of the materials.

Smooth contact algorithm
From its inception, FDEM has used a distributed potential contact force algorithm to resolve interaction between finite elements. The contact interaction algorithm relies on evaluation of the contact force potential field. The problem with existing algorithms is that the potential field introduces artificial numerical non-smoothness in the contact force, i.e., the contact forces calculated experience a jump (in amplitude and/or direction) when the contact points move from one finite element to another finite element, as demonstrated in Fig. 2. To overcome this issue, a new solution, which is named the Smooth Contact Algorithm (SCA) [44], was developed in HOSS.
In SCA, a smooth potential field is introduced according to the global geometry information of each discrete element. In particularly, SCA calculates contact potential at nodes of the finite element mesh by taking into account nodal  1 3 connectivity and existing discrete element boundaries. The contact force is then calculated as a function of the gradient of the potential field. Thus, a smooth contact evolution for a smooth surface is recovered. In addition, SCA features a generalized contact interaction law which significantly simplifies the contact force calculation and thus greatly improves the computational efficiency. Moreover, for fracture and fragmentation simulations, SCA can be naturally coupled with a Unified Cohesive Zone Model (UCZM) [45] to further improve the computational efficiency, wherein the solution concentrates the efforts of resolving contact only around the areas of interest, i.e., where fracture processes within the material are ongoing.

Unified cohesive zone model
In FDEM, Cohesive Zone Models (CZMs) usually are used to simulate the fracture and fragmentation of solids. CSMs can be classified into intrinsic and extrinsic. In intrinsic CZMs (ICZMs) the interfaces between the finite elements are active from the beginning of the simulation, while in extrinsic CZMs (ECZMs) the interfaces between finite elements are dynamically introduced when the stress state at the interface reaches some specified criteria. ICZMs are the most common type of CZM used within the FDEM community. However, ICZMs are more computationally expensive because all the interfaces are always active and therefore cohesive interaction has to be resolved in all the FDEM domains used in the simulation. Due to the same reason, ICZMs also introduce some artificial stiffness in the system which could influence the accuracy of the solutions. On the other hand, ECZMs are more efficient, since the interfaces are introduced only when fracture processes are about to start. One disadvantage of ECZMs is that, when the interfaces are introduced, this sudden change on the local topology of the model (i.e., mesh) can lead to numerical noise or chatter in the form of oscillations. This is commonly known as the "time discontinuity" problem described in [46]. In the UCZM model featured by HOSS [45], damage surfaces are dynamically inserted into the model according to the local stress state, much like the ECZM approach. However, in order to avoid sudden jumps in the simulations, the state variables are smoothly transitioned from continua to discontinua through an algorithm that properly balances the nodal forces during the process. Moreover, within the UCZM framework the point of transition from continua to discontinua (i.e., the point at which the interfaces between finite elements are dynamically inserted) is controllable through the introduction of a threshold parameter α. In Fig. 2 The contact potential field together with its contour defined by previous solutions (a) [37], and by the new solution (b). Both the jump in the amplitude and direction of the contact force calculated in the previous work can be observed when the contact point moves from one finite element to another finite element. The contact force calculated with the new solution is smooth both in terms of amplitude and direction as the contact point moves from one finite element to another. Figure modified from Lei et al. [44] Fig. 3 Comparison of 2nd Generation algorithms and the Next Generation as seen in HOSS. Figure modified from Lei et al. [31] essence, an extrinsic CZM can be retrieved from UCZM if α = 1.0; while an intrinsic CZM is obtained when α = 0.0. For an α between 0.0 and 1.0, a combined intrinsic-extrinsic UCZM is formulated.
The above enhancements not only improve the method's solid physics capability but can significantly improve speed and performance in any code that implements them. As can be seen in Fig. 3, when all three are combined, one can effectively eliminate the sharp jump seen in Y-code simulations for the artificial stiffness in the fracture model, the potential field definition in contact operations and the CST limitation for stress calculations are all avoided.

FDEM-based fluid-solid interaction solvers
In Computational Mechanics the quest to accurately resolve the interaction between fluid and solid domains has resulted in an extensive portfolio of research, methods and literature. The challenge of ensuring accurate physics is exacerbated when researchers must account for variations in time and length scales. For FDEM, the problem acquires yet another level of complexity as the solid material can break and fracture due to the loads imparted by the fluid.
Over a decade ago literature began to appear wherein FDEM community members were utilizing an enhanced fluid pressure boundary condition (BC) technique to account for fluid-solid interaction [47]. Essentially, the effect of fluid pressure on newly created fractures is accounted for by combining a source pressure-time history with an assumed pressure decay with distance from the source. The approach artificially represents, to a degree, the effects of the presence of pressurized fluid in cracks. That stated, users should apply this boundary condition carefully as the assumptions for it need to be thoughtfully established.
The Integrated Solid-Fluid Interaction Solver The BC mentioned above can be challenged by complex, real-world operations that require accuracy and reliability. To meet this challenges, in 2013 LANL unveiled its Integrated Solid-Fluid (ISF) Interaction Solver (US Patent #US10275551B2) [32,33]. The HOSS-ISF was created for simulating problems with a solid-dominated fluid interaction, such as hydraulic-fracture problems, enhanced geothermal systems, or carbon capture related. The HOSS-ISF accounts for fluid flow through fracturing porous solids, fluid flow through crack manifolds, pressure wave propagation through fluid and fluid-solid interaction.
One of the main advantages of the HOSS-ISF is that the fluid phase is described using the same grid as the solid phase via a modified Eulerian formulation. This eliminates the need of continuously mapping variables between the fluid and solid domains. Most importantly, the HOSS-ISF features an explicit time integration solver with an apertureindependent critical time step size [31]. Previous solutions have applied a similar approach but with the great shortcoming of having time step sizes that are aperture-dependent, thereby highly restricting the simulation's efficiency as aperture dimensions increase. This issue was resolved in HOSS by utilizing a pressure node at the end of each individual fracture and introducing a "fluid velocity node" at the middle of each individual fracture (Fig. 4). The pressure nodes are used to calculate the fluid pressure according to  an appropriate material law while the velocity node resolves the dynamic equilibrium of the flow system [31].
The Fluid-Solid Interaction Solver A different type of fluid-solid interaction solver is needed to simulate the effects of an explosively generated shock wave propagating through air and impinging on building structures. For this purpose the FSIS (Fluid-Solid Interaction Solver) [36] was developed and implemented into HOSS. The solver was designed to account for a fluid's compressibility and transient dynamics while also being highly robust and easy to parallelize (Fig. 5). Some of the main design parameters for the FSIS were: (a) to be explicit in terms of time integration and maintain a time step similar to the solid phase time step; (b) to allow for several fluid grids that could move with the solid; and (c) to allow for the resolution of multi-phase flow problems.
The FSIS combines a traditional FDEM Lagrangian framework (i.e., the solid solver) with fluid domains represented on Eulerian domains [36]. The solver accounts for fluid flow through fracturing porous solid, fluid flow through crack manifolds, and pressure wave propagation through fluid. The interaction between the solid and fluid domains is resolved via a modified version of the immerse boundary method [48].
An example of this newly developed FSIS capability is shown in Fig. 6. This 2D example consists on a high explosive source (TNT) emplaced inside a disk made of rock material. Between the source and the disk and outside the disk there is air at room temperature (denoted by the blue color in Fig. 6a). The explosive source is initiated at its center. It is worth noting that the FSIS solver also features a relatively simple internal energy-based reactive burn algorithm that simulates the ignition of the cells containing explosive after a certain threshold is reached. As can be seen, the FSIS can capture the breakage of the solid due to the loads imparted by the fluid. Computationally, the addition of the fluid solver to the existing HOSS framework adds less than 10% CPU overhead to the simulations.
Coupled THM Reservoir Modeling Coupled Thermo-Hydro-Mechanical (THM) simulators are necessary to explore the physical processes involved in complex subsurface operations such as unconventional fossil energy production and underground nuclear test detection. For example, hydraulic fracturing operations are designed to increase fluid pressures in specific areas of interest; as a result, this causes the rock to fracture under pressure, temperature and stress that are from ambient conditions. Shown in Fig. 7 is an example where HOSS was coupled to an Embedded Discrete Fracture Model based reservoir scale fluid solver [49]. The domain is a 100 m by 100 m square with Dirichlet boundary conditions on the left and right sides. Twenty years of injection are simulated, represented by a constant 10 MPa boundary condition on the left-hand side of the model. During the injection process the pressure at the right-hand side of the model was kept at 0 MPa. No flow boundary conditions were applied to the top and bottom of the model. The initial temperature of the rock material was 180 degrees Celsius, while the temperature of the fluid being injected was set to  50 degrees Celsius. The rock is assumed to contain two preexisting fractures in the shape of a + sign centered with the model. The effective stress in rock matrix due to the change in pressure and temperature after 20 years is shown in Fig. 7. It is clear from the figure that the presence of the fractures is perturbing the thermal stress distribution. This is a great advantage provided by HOSS, where preexisting discrete fractures in the rock media can be explicitly modeled and their influences are sensed by both the solid and fluid domains. In terms of the linking of the two codes, HOSS provides information such as porosity and permeability to the fluid solver, while the fluid solver passes fluid state variables such as pressure and temperature back to HOSS.

HOSS applications
In this section a number of numerical examples obtained with the HOSS software are presented.
Generalized Shape Dry Particle Interaction Rougier and Munjiza [50] presented a demonstration of FDEM's capability to model irregular shapes. Unlike with the Discrete Element Method (DEM), the representation of particles of general shape is very simple to achieve using FDEM.
A performance evaluation of contact interaction and contact detection algorithms was conducted utilizing the model shown in Fig. 8, which consists of a cubical, hollow raster of general shaped particles placed inside a rigid spherical container. For this case, four different types of particles were created and the raster was centered inside a rigid spherical container with no initial overlap between the particles. Each particle is given an initial velocity of 100 m/s, pointing toward the center of the container. The system is run starting from this initial state and a concentration of particles in the central area of the container is first observed, Fig. 8a. After that, a stress wave travels through the particles pushing them to move outwards until they encounter the rigid spherical container, Fig. 8c. Eventually the system reaches a steadystate point, Fig. 8d, where the particles keep moving and interacting with each other inside the container. This result is very similar to what is obtained in a molecular dynamics simulation of a gas inside a container. The main difference is that in the case shown in Fig. 8 the system is composed from solid pebbles of different shape instead of atoms.
ISF Solver Hydro-fracture As described previously, LANL researchers incorporated an integrated solid-fluid solver (ISF) within the framework of HOSS, [31][32][33]. A pressure-time boundary condition for the HOSS-ISF solver is specified on the inside face of the borehole, as shown in Fig. 9a, with the peak pressure being 20 MPa. The meshes for both cases are shown in Fig. 9b and e. The gray and red stripes oriented at 45º shown in Fig. 9e indicate the different layers of the rock. The material properties (i.e., elastic and strength properties) of the rock around the borehole and inside of the stripes for cases #1 and #2, respectively, are the same. The main difference is that the inter-layer strength properties used in case #2 are weaker that those corresponding to the rock. Snapshots for the results obtained for the two cases, taken at the same time,  Figure 9d and g display the maps of the maximum principal stress, i.e., the principal stress with the highest value (assumes tension is positive), for each case. The fracture patterns obtained in close proximity to the borehole are very similar for both cases. However, as the fractures propagate further, a significant difference is observed for case #2 when compared with case #1. This is due to the weak inter-layer properties used in the model which presents an "easier path" for the fractures to propagate (see Fig. 9g). One important aspect to highlight is that in the HOSS-ISF solver the time step size is not dependent on the fracture aperture, as is the case for other FDEM-based solvers. For the cases shown in Fig. 9, considering a maximum aperture of a max ≈ 0.012 m , the critical time step for traditional solvers is Δt fluid crit = 6.0 × 10 −12 s , while the critical time step for the HOSS-ISF solver is Δt = 5.0 × 10 −8 s . Therefore, this significant improvement introduced by the HOSS-ISF solver allows four orders of magnitude increase in the computational efficiency [31]. This example demonstrates that HOSS-ISF can simulate the effects of pre-existing interfaces on the outcome of hydraulic fracture operations.
Stick Slip Behavior in Sheared Granular Earthquake Fault Gouge Among the many numerical methods used to describe the evolution of granular systems, the DEM has been the most popular. However, DEM is incapable of capturing detailed deformation related to stress distributions and elastic wave propagations in the system under study.
The FDEM, on the other hand, provides a natural solution to modeling this system since it not only resolves the interaction between the different parts of the granular system but it also provides a solution for stress wave propagation inside the solid media. Gao et al. [51,52] used HOSS to explicitly simulate a 2D sheared granular fault system including both gouge and plate, as shown in Fig. 10. The first step in the modeling of this system is to apply a pre-defined normal stress on the gouge material. This was achieved by applying a confining stress to the bottom stiff bar while the top one was kept fixed. The system is allowed to reach equilibrium under this confining stress by using a velocity relaxation technique. The next step began from the start by allowing the system to slowly shear by moving the top stiff bar to the right with constant speed v. The shearing action is transmitted from the plates to the gouge via a set of teeth placed on the faces of the plates that are facing the gouge (see inset in Fig. 10). As the system was sheared, the gouge particles resisted the load via a number of force (or stress) chains, as shown in Fig. 11a where the map of Cauchy shear stress (τ xy ) is presented. However, at certain points in time, the gouge particles were not able to resist the load in the current configuration, giving way to the occurrence of slip events. After each slip event the force chain structure inside the gouge changed considerably as is demonstrated by comparing both sides of Fig. 11 [51,52].
The HOSS simulation results show good agreement with the physical experiment in terms of its capability to capture the stress chain development (Fig. 11) and the distribution of the size of seismic events [51,52]. In this case, the force or stress chains are emergent features of a complicated system composed of many different interacting parts. This example demonstrates that HOSS is a very powerful tool for the analysis of the gouge material failure and its consequences on earthquake rupture processes. The emergent property results (stress chains) achieved here clearly define how FDEM can be used as "virtual experimentation" tool. Instead of having to re-tool the experimental apparatus researchers can now simply adjust gouge dimensions or initial conditions to further analyze slip events in earthquake fault.
High Velocity Impact on Unconsolidated Beds In the context of the NASA InSight mission, which placed the first seismometer on the surface of Mars, interest was focused on a better understanding of seismic sources generated by impacts. HOSS was used to investigate the properties of impacts in regolith and to develop and validate a model of this unconsolidated target material under extreme loading.
This validation was conducted with data from an experimental campaign conducted in June 2012 at the NASA/ AVGR facility [53]. The setup consisted of a closed cylindrical tank, 1 m in radius and in height, filled with a controlled atmosphere and target bed. These experiments used a pumice sand target bed with grain size of 0.1-0.2 mm and porosity of 62%. The shot reproduced here used a 6.3 mm radius Pyrex bead of mass 0.29 g as impactor. The case simulated corresponded to a vertical impact with an initial velocity of 0.98 km/s. Since the model is symmetric, and with the aim of reducing the requirements of computational resources, a 30º wedge of the tank was simulated. In the FDEM framework the sand granular structure and grain displacements are accounted for by allowing the Lagrangian finite elements to move freely with respect to each other while they interact with each other via frictional contacts. A snapshot of the results obtained using HOSS is shown in Fig. 12 where the colors indicate the magnitude of the particle velocity. The top-left corner of the image shown in Fig. 12 shows the cratering process being created by the action of the projectile. Very good agreement was obtained between the calculated velocity-time histories for sensors located along a vertical line below the impact point and the observed records [54]. This example demonstrates that HOSS can be used to reproduce high velocity impacts on unconsolidated materials.
Modeling of Split Hopkinson Pressure Bar (SHPB). Rougier et al. [55] presented a classic example of FDEM laboratory replication for a 3D Split Hopkinson Pressure Bar (SHPB) simulation, as shown in Fig. 13.
For this particular type of experiment, FDEM can be used to model the entire experimental setup starting from the pressure wave traveling though the incident bar to its ultimate resolution in the transmission bar. The simulation results seen in Fig. 14 show great agreement with the physical experiment in terms of actual replication of breakage and fragmentation. In addition, the simulation also captures relevant experimental measurements as seen by the evolution of the calculated tensile stress as a function of time. Most computational mechanics simulations readily capture the "hardening" portion of the experiment, i.e., up to and including the peak, while little attention is payed to capturing the actual "softening" or post-peak failure behavior. However, as  1 3 a fracture and fragmentation tool, it is imperative for FDEMbased models to replicate the post-peak with the consequent fracture and fragmentation processes, as is seen here for the 3D solid mechanics simulation. The exceptional agreement between HOSS' results and the experimental observations demonstrate that models comprising of full experimental  [53] setups concentrating on rock (or other material type) failure can be properly represented in HOSS.
Uniaxial Compression Experiment Fracture Coalescence FDEM is an ideal tool for the numerical simulation of fracture processes; this provides an opportunity to better understand the underlying mechanisms of crack initiation, propagation, and coalescence in geologic materials. As stated by Euser et al. [53], fracture coalescence is a critical phenomenon that creates large, inter-connected fractures from smaller cracks, affecting fracture network flow and seismic energy release potential. An example of this is shown in Fig. 15 where the results obtained for a numerical study of a uniaxial compression test on a rectangular granite sample with a pre-existing flaw in its center are presented. Figure 15a presents a time evolution of the fractures developing inside of the sample. In this case, the amount of damage indicates the degree to which the bonds at the interface of the finite elements have degraded, i.e., damage equal zero means intact bonds between the finite elements, while damage equal one means bonds completely broken between the finite elements. According to the type of loading experienced by each interface the bonds can fail in three main ways: pure tension, pure shear, and mixed tension/shear, as is shown in Fig. 15b. Figure 15c shows the evolution of the maximum principal stress (i.e., the greatest value of the principal stress-tension being consider positive) with time as the fractures develop inside the sample. One important feature to highlight is the complex stress state developed inside the sample as the numerical experiment progresses. This is particularly evidenced by the stress concentrations around the tips of the fractures and at both sides of the pre-existing flaw, as shown in Fig. 15c.
The type of analysis presented by Euser et al. [53] can also be extrapolated to larger scales. For example, when developing mining operations there is a need to understand the behavior of mine pillars under complex loading conditions [56]. The mine pillars may contain pre-existing flaws or fractures that will degrade their loading response, making them more susceptible to experience failure in the form of strain burst or pillar burst [57]. A demonstration of HOSS for a problem similar to that of the mine pillar one is shown in Fig. 16, where a 25 m × 25 m granite block with two embedded in situ fracture sets is subjected to uniaxial loading imposed by keeping the bottom plate fixed while moving the top plate downwards with constant speed. Figure 16b displays the initial and final stage of the fracture map. As can be seen there is a high degree of complexity and the method appears to perform quite adequately in representing fracture coalescence for much larger scale problems. This demonstrates that HOSS can be used to analyze complex real-world size fractured rock systems.
Triaxial Direct-Shear Experimental Fracture Results obtained from research efforts aimed at improving the understanding of fracture initiation and propagation in shale rocks using HOSS is shown here. A key factor in any simulation of hydraulic fracturing processes is to correctly represent the effects of fluid pressure on the faces of the fractures. In order to address this an enhanced fluid pressure BC was implemented, where fluid pressure on the fractured surfaces was obtained as a function of both the time and the distance from the fluid source. The fluid pressure calculated from the solver is applied to the original free faces as well as the faces created by fracturing. This solver was used in the work presented by Carey et al. [58] and Lei et al. [47] where an experimental shale sample (Fig. 17) was loaded in direct shear conditions by the action of two opposing semi-circular, rigid anvils. The focus was to characterize the behavior of a generic shale sample and compare this with experimental observations.
Simulation results compared to triaxial core flooding experiments are presented in Fig. 17. The fractures that cross and propagate along horizontal bedding planes were successfully simulated using the experimental stress conditions and injected fluid pressure within the fractures. The simulation captures the primary fracture network (red lines) as well as bedding plane fractures and densely fractured regions.
The enhanced fluid pressure BC used in HOSS is therefore able to accurately replicate the main features of a fracturing process while showing the importance of fluid penetration into fractures as well as layering in determining final fracture patterns.   [59] As the complexity of the experiment increases, the method continues to perform more than adequately. The experimental and simulation objectives were to predict fracture-permeability behavior, i.e., fracture formation under stress, effect of fracture formation on fluid conductivity and the effect of reactive fluid-fluid and fluid-rock processes. Shown in Fig. 18 are a 2D view of stress propagation in a 3D simulation sample (left) and a 3D view of the damage propagation inside of the sample (right).
Flyer Plate Dynamics In the work presented by Chau et al. [59] HOSS is used to simulate the response of a flyerplate impact experiment in a Westerly granite sample that contains a randomized set of cracks. The model (Fig. 19) consists of 3 plates: aluminum flyer plate, granite target plate, and PMMA momentum trap. Within the sample, there are a number of pre-existing cracks each with initial lengths determined by a power law distribution, and randomly distributed orientations. A Velocity Interferometer System for Any Reflector (VISAR) is used as the interferometer system to measure the free rear surface particle velocity history and in the HOSS model a sensor to record the particle velocity is placed on the back of granite sample. As shown in Fig. 19, HOSS can accurately reproduce the VISAR plot and capture the spall region and stress peaks observed in flyer plate experiments in granite. The number and the distributions of pre-existing fractures was also studied to get better understanding of the effect of structural cracks on the mechanical behavior and the failure path of Westerly granite under high strain rate impact. Moreover, the simulation results for the spall strength are in quite good agreement with the actual experimental values (Table 1).
Extending the above technique to other materials like metals is a matter of incorporating appropriate material models. In this case, utilizing a relative simple metal plasticity model (Fig. 20) HOSS does an exceptional job of capturing the nonlinear dynamics of the problem. Additionally, as seen in Fig. 20 the ultra-high strain rate dynamics create a highly volatile region of fracture in the impacted metal target plate.
Mixed-mode Failure in Carrara Marble Boyce et al. [60] investigated the formation of hybrid fractures under mixed tensile and compressive stresses using the HOSS FDEM. There are two common types of fracture in brittle geo-materials, tensile and shear.
Modeling a series of indirect tension tests (see Fig. 21a) with confining pressures that increase from 30 to 150 MPa allows for the development of stress states adequate to generate a spectrum of mixed-mode failure conditions. In Fig. 21b the temporal propagation of a single mixed-mode fracture process is shown. The cooler colors show shear fracture domination for initial failure at higher confining pressures near the edge of the specimen. Once the break begins, the axial force acts on the reduced cross-sectional area thus increasing the tensile stress and causing tensile damage in the center of the specimen. Figure 21c shows a progressive change in fracture orientation as the confining pressure increases.

HOSS-FSIS High Explosives Modeling
A key feature of the fluid solver within the HOSS-FSIS, Munjiza et al. [36], is that the governing equations are resolved using different control volumes, resulting in discretization errors of either third, second or first order. Also, the fluid solver is fully explicit with respect to the time integration and therefore conditionally stable, with a time step that is synchronized with the FDEM solid time step. In this section three demonstration examples for high explosives (HE) modeling are   Fig. 22 shows the setups for four cases used to study the effect of porosity within a generic HE material that is placed inside a rigid container. Case I considers the solid energetic material with no porosity, Case II takes into account a 40% porosity within the energetic material, Case III: presents an HE charge that is composed of 2/3 solid and the top 1/3 with 40% porosity, and Case IV is the mirrored situation of Case III, i.e., energetic material with 1/3 with 40% porosity and 2/3 solid on the top. In all cases the amount of energetic material is the same. The initiation is achieved by the mechanical action of a slapper or flyer plate impacting the HE at high speed. Results are shown in the bottom portion of Fig. 22 where temperature maps for each case are shown for the same point in time. As expected, Case I exhibits a very uniform temperature map, while Case II presents a quite convoluted picture with hot and cold spots inside of the HE caused by the porosity. A similar effect can be observed in Cases III and IV's porous regions. The wave fronts propagating toward the opening of the container for Cases I and IV are quite uniform and flat. This is quite different to what is seen in Cases II and III where the HE porosity considerably distorts the wave front. These cases demonstrates the applicability of the FSIS for the study of pressure wave propagation in HE materials when they are interacting with solid parts.
HOSS-FSIS Demonstration #2. A similar type of problem, but with a different setup is presented in Fig. 23 where a 10% porosity HE-like material is emplaced inside a 5-mmthick copper disk, see Fig. 23a. The pores inside of the HE material, as well as the surroundings of the disk are filled with air. The HE material is detonated at the center of the disk. Similar to what was observed in the previous example, the presence of the pores generates a complex behavior inside of the HE charge, as it is demonstrated by the results shown in Fig. 23b. Even at later times during the simulation the pressure maps inside of the fracturing disk are quite varied, as is demonstrated by the picture shown in Fig. 23c. As the copper shell breaks, pathways are created for the high pressure, high temperature gases to escape into the surrounding atmosphere. This is evidenced by Fig. 23d, where a map of the velocity of the gasses is presented and the jetting process can be clearly seen. This example demonstrates the potential of HOSS-FSIS to study the evolution of HE charges (in many different shapes) and how the HOSS-FSIS Demonstration #3. Following the previous examples, this demonstration shows that HOSS-FSIS can also be used to simulate the propagation of acoustic signals due to explosive charges. In Fig. 24a a typical case where a source is emplaced inside a room that has rigid walls and rigid cylindrical obstacles (pillars) is presented. The results obtained are shown in Fig. 24b where a complex pressure map is observed inside the room, which will also determine HOSS-FSIS High Energy Blast Event Modeling Two idealized high energy blast event cases are shown here to illustrate that the FSIS solver can be used to capture structural response and debris flow. The model setup is similar for the two implementations in that there is a representative tunnel portal with a blast door emplaced some distance into the tunnel.
In Fig. 25a an explosive source is emplaced on a tunnel portal steel blast door that is 0.5 m thick. The results obtained are shown in Fig. 25b, where it evidenced that the use of HOSS-FSIS for this type of modeling can provide a plethora of information that would be unavailable from mechanical only analysis results or one-way coupling alternatives (i.e., no true fluid-solid interaction). First, with this event we see total obliteration of the blast door as well as a resultant 1 GPa pressure wave traveling down the tunnel itself. Just as important the particle velocity of the damaged material from the blast door and the rock is traveling at 1600 m/s when the calculation is stopped. The advantages of using a Lagrangian solver for damage assessment here is self-evident as it opens the door for having a capability that can better capture portal obliteration and pyroclastic flow effects.
The setup for a demonstration of the effects of a larger explosive source emplaced at a 10 m range and 15 m HoB from the portal entrance is shown in Fig. 26a. In the right panel of Fig. 26b, the rock material is being lofted because of material vaporization (Tillotson EOS [61] is being utilized in this calculation). In the future, the vaporization processes will be handled by specially designed donation algorithms which will allow the transfer of the vaporized, melted or rubbled finite elements from the Lagrangian to the Eulerian domain. This will enable researchers to better handle extreme loading conditions (leading to vaporization  and melting) which lead to highly distorted elements that usually present a challenge for regular Lagrangian-based codes.
Tunnel Blast Door Penetration In this section, HOSS' capability to capture high strain rate penetration debris flow is demonstrated. A cross section of the 3D model used for this demonstration is shown in Fig. 27.
The results obtained for the 3D multiple door impact event are shown here in Fig. 28. The calculation itself is stable and the HOSS not only captures the high velocity ejecta but also clearly shows the resultant rock damage caused by the high stresses the entrenched steel door experiences. The code's capability to capture debris laden flow is evident and it shows the effect the first door's ejecta flow has on the subsequent doors.
Hypervelocity Impact Traditionally, empirical laws have been used to determine the seismic efficiency of bolide or asteroid impacts. However, empirical laws are not always reliable since they provide seismic efficiency estimates that vary from 10 −2 to 10 −5 . Recent developments in FDEM made numerical modeling a promising approach to investigate the physics of impacts to provide answers to the outstanding question of seismic efficiency. The simulations shown here consider a 1 km diameter body impacting a 15 km diameter half space, which is based on a hypervelocity verification/validation paper from Pierrazzo et al., 2015 [62], see Fig. 29. Data are collected at a variety of locations and are then compared to data from other hydrocodes simulations of the same problem. The results reported in here concentrate on an impact velocity of 20 km/s for an oblique impact.
A snapshot showing the wave propagation and the material being ejected for the 45 degree angle impact with an initial projectile velocity of 20 km/s is shown in Fig. 30 along with comparisons of the calculated peak pressure with distance. There is good agreement between the results  Earthquake Dynamic Rupture Modeling Utilizing HOSS FDEM, a new physics-based dynamic earthquake rupture modeling framework to investigate the fundamental mechanisms of coseismic off-fault damage and its effect on rupture dynamics, radiation, and the overall energy budget has been developed by Okubo et al. [63] and Okubo [64]. FDEM was chosen to model the dynamic earthquake rupture with the coseismic off-fault damage simply because it allowed for the activation of both off-fault tensile and shear fractures based on prescribed cohesion and friction laws that could quantify the effect of the coseismic off-fault damage. Other state-ofthe-art numerical techniques generally used for earthquake rupture modeling were unable to describe the detailed offfault fracturing processes (Mode I tensile and Mode II shear fractures) due to their computational and model formulation limitation. Briefly shown here is one aspect of the cross validation of FDEM against community accepted earthquake models and then an application of the above framework to the dynamic earthquake rupture modeling of 2016's Kaikoura earthquake in New Zealand.
A cross validation of HOSS FDEM was performed [63] to assess its accuracy against other numerical schemes: the finite difference method (FDM), the spectral element method (SEM), and the boundary integral equation method (BIEM) as these approaches have been extensively verified in previous studies. Figure 31 shows the comparison of slip velocity history at 9 km from the center of the main fault where the results of HOSS are compared to FDM, SEM, and BIEM. The slip velocity history of HOSS is consistent with the other numerical schemes except for the peak slip velocity and the rupture arrival time is slightly faster than the others. These small differences are explained by the artificial viscous damping which coincidentally also significantly removes the high-frequency numerical noise.
Fracture damage patterns around faults, induced by dynamic earthquake rupture, provide an invaluable record to clarify rupture processes on complex fault networks. The 2016 Mw 7.8 Kaikoura earthquake in New Zealand was   29 General model setup for the hypervelocity impact simulations, a normal (90° with respect to horizontal) impact and oblique (45° with respect to horizontal) are considered one of the most complex earthquakes ever documented that ruptured at-least 15 faults. High-resolution optical satellite imagery provided distinctive profiles of fault displacement fields that aided in visualizing the off-fault damage pattern. In the work presented by Klinger et al. [65] this observational approach, coupled with the HOSS FDEM numerical approach captured the rupture and activation of off-fault damage, thereby allowing researchers to determine Kaikoura's most likely rupture scenario. Of the several nucleation options studied, Fig. 32 highlights the Papatea fault nucleation option wherein confirmatory seismic instrumentation compared to rupture simulation results show the earthquake started in the southern portion of that fault structure. Furthermore, as demonstrated by Klinger et al. [65] compelling simulation fracture damage results when compared to observational mapping confirm that the dynamic earthquake rupture modeling framework can be a mainstay in the earthquake community for analysis of these highly hazardous natural phenomena events.

Conclusions
It is incumbent upon the FDEM community to develop benchmark problems for discontinua fracture simulations that showcase the hybrid method's capability to capture relevant physics and demonstrate experimental replication; as stated by Nobel Laureate Richard P. Feynman: "It doesn't matter how beautiful your theory is, it doesn't matter how smart you are. If it doesn't agree with experiment,  it's wrong." FDEM is making great strides in establishing itself as a mainstay in the field of Computational Mechanics of Discontinua.
In this paper, a collection of examples combining previous research together with new results obtained using HOSS has been presented. The results demonstrate that a wide range of length (from sub-millimeter to several kilometers) and time (from fractions of a nanosecond to months) scales can be addressed with FDEM-based simulation technology. The power and versatility of FDEM to describe very disparate material behavior, starting from pure continuum, to fracturing solids, and all the way to full discontinuum systems while considering the effects of the presence of fluids has also been evidenced.
FDEM has an ever expanding user-base and an ever expanding field of application. Research groups who use the method are positioned in the six major continents and university usage is widely increasing. In this paper the latest HOSS developments have been summarized demonstrating that this numerical tool can facilitate multidisciplinary and interdisciplinary work.
Acknowledgements The authors would like to thank the Los Alamos National Laboratory's Feynman Center for Innovation for their support during the development and distribution of HOSS. HOSS is available for commercial licensing, more details can be obtained in the following link: https ://www.lanl.gov/proje cts//feynm an-cente r/ index .shtml . Also, for universities only, HOSSedu [34] is available for non-commercial licensing. HOSSedu is a parallel 2D/3D FDEM dry fracture platform with 2D fluid implementations. This research used resources provided by the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. Department of Energy National Nuclear Security Administration under Contract No. 89233218CNA000001. By acceptance of this article, the publisher recognizes that the U.S. Government retains a nonexclusive, royaltyfree license to publish or reproduce the published form of this contribution, or to allow others to do so, for U.S. Government purposes. Los Alamos National Laboratory requests that the publisher identify this article as work performed under the auspices of the U.S. Department of Energy. Los Alamos National Laboratory strongly supports academic freedom and a researcher's right to publish; as an institution, however, the Laboratory does not endorse the viewpoint of a publication or guarantee its technical correctness. This manuscript has been approved for unlimited release under Los Alamos National Laboratory's LA-UR-20-20646.

Compliance with ethical standards
Conflict of interest On behalf of all authors, the corresponding author states that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.