Commercial mining / quarrying / civil-engineering operations

Abstract

A method for predicting solid material dislodgement due to explosive blasting, the method including: representing pre-blast solid material with a set of discrete particles, each with a 3D location and one or more mechanical properties of the solid material at that 3D location before the explosive blasting; representing at least one explosive blast with a blast time and a set of explosive particles, each with a 3D blast location in the pre-blast solid material and an initial internal energy; repeatedly, for a plurality of time steps, defining a 3D grid surrounding the discrete particles (e.g., plus two empty grid cells in each direction), and numerically determining / calculating updated velocities and updated 3D locations of the discrete particles based on the mechanical properties and the 3D locations, using the 3D grid; for at least one of the time steps corresponding to the or each blast time, including the explosive particles with the discrete particles to numerically determine the velocities and updated 3D locations of the discrete particles, and numerically determining / calculating the updated velocities and the updated 3D locations at the time steps after the or each blast time using the initial internal energies; and determining post-blast locations of the solid material using the repeatedly updated 3D positions of the discrete particles and the pre- blast solid properties of those discrete particles.

Description

- I -COMMERCIAL MINING / QUARRYING / CIVIL-ENGINEERING OPERATIONSRELATED APPLICATION

[0001] The present application is related to U.S. Provisional Patent Application No. 63 / 612810, filed on December 20, 2023, the originally filed specification of which is hereby incorporated by reference herein in its entirety.TECHNICAL FIELD

[0002] The present disclosure relates to explosive blasting planning and execution based on numerical simulations for commercial mining / quarrying / civil-engineering operations, including open-cut mining and underground mining.BACKGROUND

[0003] In commercial mining / quarrying operations, valuable rock is identified by geological surveying prior to blasting. When the rock is blasted, the valuable rock moves, and can be difficult to find during excavation. In civil-engineering operations, it may be difficult to predict where fragments of solid materials (e g., concrete) from a blast (or shot) will travel, which can make safe planning difficult.

[0004] Prior attempts to track material movement (thus dislodgement) during explosive blasting are either geometrical (referred to as "geometric approaches") or based on physics simulations (referred to as "physics-based approaches").

[0005] Geometric approaches, like "Smart Vectors" from ORICA, rely on heuristics to permutate the positions of post-blast materials from their initial positions to their final positions based on the location of empty cells between the floor (which is substantially notblasted) and the post-blast topography (which can be measured using non-contact surveys, e g., using LIDAR); however, these heuristics cannot account for complex mechanical behaviors during blasting, e.g., deformations, friction, rolling and sliding within the solid material. Experiments have shown incorrect distribution of the post-blast material density, requiring other geometric post-processing — like numerical relaxation methods — to attempt to distribute the post-blast blocks evenly; however, such relaxation methods tend to converge extremely slowly or not at all in complex cases with dense cluster of rocks or thin areas between the two topographies.

[0006] Physics-based approaches can account for more mechanical specificities of the solid material (e.g., rock) during a blast, but prior approaches are either restricted to two dimensions (2D) or long-running three-dimensional (3D) processes (taking several hours or even days of computation times) due to the high velocity and stiffness of the approach. None of these prior approaches result in satisfying solutions, e.g., due to insufficient swell in the result, and / or insufficient power through on the back of the blast, and / or insufficient matching to the floor and post-blast topographies. For example, "I-Blast Ultimate" is a simulation engine for blasting from Thierry Bernard Technology and DNA Blast Group (Nice, France) that simulates rock movement due to blasting in order to estimate where pre-blast rocks are in post-blast material; however, this simulation engine is undesirably slow and / or the simulated volume undesirably small / imprecise for at least some commercial mining / quarrying / civil-engineering operations.

[0007] It is desired to address / ameliorate at least one difficulty / limitation of the prior art, or to at least provide a useful alternative.SUMMARY

[0008] Disclosed herein is a method, performed with one or more physical processors, for predicting solid material dislodgement due to explosive blasting, the method including: a. representing, with the one or more physical processors, pre-blast solid material with a set of discrete particles, each with a 3D location and one ormore mechanical properties of the solid material at that 3D location before the explosive blasting; b. representing, with the one or more physical processors, at least one explosive blast with a blast time and a set of explosive particles, each with a 3D blast location in the pre-blast solid material and an initial internal energy; c. repeatedly, , with the one or more physical processors and for a plurality of time steps, defining a 3D grid surrounding the discrete particles (e.g., plus two empty grid cells in each direction), and numerically determining / calculating updated velocities and updated 3D locations of the discrete particles based on the mechanical properties and the 3D locations, using the 3D grid; d. for at least one of the time steps corresponding to the or each blast time, including the explosive particles with the discrete particles to numerically determine the velocities and updated 3D locations of the discrete particles, and numerically determining / calculating the updated velocities and the updated 3D locations at the time steps after the or each blast time using the initial internal energies; and e. determining, with the one or more physical processors, post-blast locations of the solid material using the repeatedly updated 3D positions of the discrete particles and the pre-blast solid properties of those discrete particles.

[0009] The repeated numerically determining of the velocities and the updated 3D locations may include performing numerical calculations of the velocities and the updated 3D locations using a plurality of graphics processing units (GPU) in parallel.

[0010] The repeatedly numerically determining of the velocities and the updated 3D locations may include using a quadratic B-spline function to combine the mechanical properties of ones of the discrete particles adjacent to each discrete particle.

[0011] The repeatedly numerically determining may include a Galerkin-style moving least squares (MLS) discretization.

[0012] The mechanical properties may include a Poisson coefficient value and a Young modulus value, and the numerically determining / calculating may include determining / calculating the velocities and updated 3D locations using the Poisson coefficient values and the Young modulus values.

[0013] The mechanical properties may include a tensile strength value, and the method may include repeatedly, for some or all of the plurality of time steps: a. determining a peak effective energy release rate for each discrete particle to determine a fragmentation status for that discrete particle, wherein the fragmentation status is determined to be fragmented if the peak effective energy release rate for that discrete particle is above a predetermined threshold value; and b. numerically determining / calculating the velocities and the updated 3D locations of the discrete particles based on their respective fragmentation statuses, including reducing / zeroing the tensile strength value in the mechanical properties to model fluid-like behavior when the fragmentation status represents fragmented.

[0014] The method may include determining, with the one or more physical processors, the fragmentation status of one or more of the discrete particles in the pre-blast solid material based on a fragmentation input (e.g., a user-defined volume representing prefragmented material that are affected by the eigenerosion model in the first numerical determination / calculation to account for previous blasts that would have already broken part of the solid material).

[0015] The mechanical properties may include a friction angle and friction hardening coefficients, and the numerically determining / calculating may include determining / calculating the velocities and updated 3D locations using the friction angle and the friction hardening when the fragmentation status represents fragmented.

[0016] The method may include increasing, with the one or more physical processors, a particle volume represented by each discrete particle to represent swelling of the solid material due to the explosive blasts, including none or more of the following: a. increasing the particle volume based on its updated velocity; b. increasing the particle volume based on its updated spin tensor; and / or c. increasing the particle volume based on its determinant of rate of deformation tensor.

[0017] The method may include removing, with the one or more physical processors, the explosive particles from the discrete particles after a pre-defined duration to emulate gaseous release into the air through cracks in the solid material.

[0018] The method may include representing the explosive particles with explosive material properties before and during ignition.

[0019] The method may include representing the explosive particles in their gaseous states by a Jones-Wilkins-Lee (JWL) model.

[0020] The 3D grid can include a grid spacing between 0.1 meters and 10 meters in each direction, e.g., substantially 1 meter.

[0021] The time steps may each have a duration of substantially 10 microsecond to 1000 microseconds.

[0022] The numerically determining / calculating the updated velocities may include explicit integration of forces

[0023] The method may include defining values of the 3D locations and the mechanical properties by discretizing, with the one or more physical processors, a block-model in the form of a set of non-overlapping cuboids representing the pre-blast solid material with material attributes attached to each cuboid.

[0024] The method may include automatically extending, with the one or more physical processors, the block-model using null blocks based on a user-defined floor and / or a measured pre-blast topography of the blast area to form the block-model to discretize.

[0025] The method may include defining the 3D blast locations by discretizing one or more explosive volumes in respective blast holes in the pre-blast solid material, including defining the 3D blast locations by a random / quasi-random uniform distribution in the or each explosive volume.

[0026] The method may include representing an unbreakable floor under the pre-blast solid material as a heightfi eld, and using the heightfield as a boundary condition in the repeatedly numerically determining.

[0027] The method may include repeatedly, , with the one or more physical processors and at one or more of the time steps, one or more of: a. reducing / eliminating the velocity values of the discrete particles below the unbreakable floor; b. reducing / eliminating downward velocity values of discrete particles adjacent to the unbreakable floor; and c. reducing / eliminating tangential velocity values of discrete particles adjacent to the unbreakable floor depending on the normal velocity and the distance of the discrete particle from the unbreakable floor.

[0028] The method may include fitting, with the one or more physical processors, the discrete particles to a measured post-blast topography, including by: a. applying an attraction force to the discrete particles to move the discrete particles into a post-blast volume of the measured post-blast topography; and / or b. applying a suction force to the discrete particles to move the discrete particles to fill the post-blast volume of the measured post-blast topography.

[0029] Disclosed herein is computer-readable storage having stored thereon computer- readable instructions configured to cause a system that includes the one or more physical processors in the form of a CPU and one or more GPUs to perform the method above.

[0030] Disclosed herein is a system including the one or more physical processors and the computer-readable storage above.|0031| Disclosed herein is a mining method, including the method above; and physically measuring the pre-blast solid material to determine the one or more mechanical properties of the solid material before the explosive blasting, and / or physically excavating and / or routing the solid material after the explosive blasting based on the determined post-blast locations of the solid material.BRIEF DESCRIPTION OF THE DRAWINGS

[0032] Some embodiments of the present invention are hereinafter described, by way of example only, with reference to the accompanying drawings, in which: a. FIG. 1 A is a flow chart of a method disclosed herein; b. FIG. IB is a schematic diagram of a system configured to perform the method of FIG. 1A showing memory transfers during an initialization subprocess of the method; c. FIG. 1C is a schematic diagram of the system showing memory transfers during a main simulation loop of the method; d. FIG. 2 is a schematic side view of a block-model used in the method; e. FIG. 3 is a screen shot of a topography mesh of the bench used in the method; f. FIG. 4 is a schematic side view of the block-model of FIG. 2 with additional null blocks added by the method,g. FIG. 5 is a screen shot of a defined volume defined by the topography mesh of FIG. 3 and discretized with a block-model, wherein dark and light blocks have been provided by a user, and medium-grey blocks are the additional null blocks added by the method; h. FIG. 6 is a schematic side view of the block-model of FIG. 4 following discretization by the method; i. FIG. 7 is a screen shot of the block-model discretized into particles by the method, wherein each particle has a diameter smaller or equal to a cellwidth; j . FIG. 8A is a vertical cross-sectional view of a discretized block-model with discretized explosive particles representing bulk explosive material in a blast hole, k. FIG. 8B is an enlarged view of a portion of FIG. 8A; l. FIG. 9 is a screen shot of the defined volume of FIG. 3 including a portion selected and marked as pre-fragmented, e g., due to a previous blast; and m. FIG. 10 is a vertical cross-sectional schematic view of discrete particles being fitted into a post-blast volume by artificial attraction and suction forces.DETAILED DESCRIPTIONOverview

[0033] Described herein is method 100 for predicting solid material movement (e g., dislodgement of natural solid materials such as rock and / or manufactured solid materials such as concrete / brick) due to explosive blasting. Also described herein is a system 150 configured to perform (or "execute") the method 100. The explosive blasting may be incommercial mining / quarrying / civil-engineering operations, including open-cut mining and / or underground mining.

[0034] The method 100 provides a physics-based approach for predictive and reactive movement simulations applied to explosive blasting. The method 100 can track / simulate solid material dislodgement more rapidly / accurately than prior methods / systems, e.g., tracking / simulating substantially 100k particles in substantially 10 seconds (s) using at least one Graphic Processing Unit (GPU) with an accuracy of about 1 meters (m) at all points. The method 100 can combine both speed (e.g., sub-minute execution times) and accuracy (e g., an average distance error sub 2 m between the actual / measured post-blast topography and the simulated post-blast topography). In experimental examples, the method 100 has shown exceptional computation times (sub-minutes execution times) while delivering sufficient simulation results (sub 2 m of average height difference error), including when the simulation parameters were selected, e.g., manually, e.g., using a set of blast models from the same operational site (e g., a mine) with a known post-blast topography (measured after the blast though photogrammetry), including with a human operator adjusting the selected values of a subset of the simulation parameters (including a coefficient multiplying an explosive strength in a Jones-Wilkins-Lee (JWL) model, a rock friction angle in a rock plasticity model, a rock crack threshold for an eigenerosion model, and / or a swell deformation rate multiplier for artificial swelling) and adjusting the cellwidth to keep the simulation times around 1 minute on the biggest model from the dataset of the operational site. In optimized experimental examples, the human operator ran the method 100 on the selected models and observed the simulated result; then, based on the errors observed (e.g. too much swell, not enough rock breaking, not enough movement, etc.), he or she adjusted the values of the simulation parameters (i.e., their selected values) and re-ran the simulations of the method 100 repeatedly until only a reasonably small difference (or "error") persisted between the simulation result and the real post-blast topography, e.g., less than substantially 1.5 meters, e.g., automatically quantified by measuring an average absolute height difference between the simulated post-blast topography and the real-world (measured) post-blast topography for that operational site.

[0035] The method 100 may differ from one or more pre-existing approaches by: using improved physics simulation subprocesses, applying a material point method (MPM) fornumerical calculations, running on one or multiple GPUs in parallel, applying artificial swelling, and / or enabling post-processing for reactive movement, among other features described hereinafter.

[0036] In mining / quarrying operations, the method 100 can predict rock movement due to typical explosive blasting operations, e.g., representing bulk explosive material filled into blastholes by a plurality of explosive material portions (specifically: gas particles, sa described hereinafter) distributed in specific 3D locations (in explosive volumes in the blastholes) in / across a pre-blast volume and activated (specifically: inserted) with selected timing, e.g., sub-ms timing, e.g., corresponding to blasting equipment provided by ORICA's Blasting Services products.

[0037] The method 100 can reduce energy usage during excavation, maximize ore extraction during excavation because ore rock locations are shown in the muck pile, and / or minimize ore dilution during post-blast processing (which includes excavation, trucking and ore extraction) The method 100 may allow for optimization of a planned blast (or "shot") before blasting (e.g., by selection of the bulk explosive type, the blasthole dimensions, the blasthole locations, and / or the selected blast timings) to obtain a desired movement direction / distance and / or fragmentation level. The method 100 may also allow for optimization of the post-blast material collection based on the knowledge of the locations of the pre-blast rocks within the post-blast muck pile, e.g., 3D maps of ore locations in the post-blast muckpile. The method 100 may be integrated with a commercially available blast movement modelling system, e g., OREPRO 3D from ORICA, that models blast movement to enable improved situational awareness, and grade control, safety, blast control, and economic outcomes.Method 100

[0038] The method 100 is performed (or "executed") by the system 150, which includes (as shown in FIGs. IB and 1C) a combination of at least one central processing unit (CPU) 154 including corresponding random-access memory (RAM), forming a CPU side 10A of the system 150, and one or more graphic processing units (GPUs) 156, each including corresponding video random-access memory (VRAM), forming a GPU side 10B of thesystem 150, and computer-readable storage 152, e.g., disk memory and / or cloud memory, which is readable by at least the CPU 154, having stored thereon computer-readable instructions configured to cause the or each CPU 154 and the GPUs 156 to perform (or "execute") the method 100, including the subprocesses thereof described hereinafter. As shown in FIGs. IB and 1C, the storage 152 is configured to communicate with the CPU 154 (including to transfer data representing 200, 202, 204, 206 as described hereinafter), and the CPU 154 and GPUs 156 are configured and mutually connected to communicate with each other (including to transfer data representing 204B, 202B, 600, 700, 702, 704, 706 as described hereinafter).

[0039] The method 100 includes the following subprocesses, as shown in FIG. 1 A: a. an initialization subprocess 102 performed by both the CPU side 10A, and the GPU side 10B; b. an update explosives subprocess 104 performed by the CPU side 10A and partially by the GPU side 10B when uploading / removing particles, following and responsive to the initialization subprocess 102; c. an update sparse grid subprocess 106 performed by the GPU side 10B, following and responsive to the update explosives subprocess 104, d. an estimate substep length subprocess 108 performed by the GPU side 10B, following and responsive to the update sparse grid subprocess 106; e. a particle and grid update subprocess 110 performed by the GPU side 10B, following and responsive to the estimate substep length subprocess 108; f. a force integration and boundary handling subprocess 112 performed by the GPU side 10B, following and responsive to the particle and grid update subprocess 110; g. a more substeps determination subprocess 114 performed by the CPU side 10A, including reading a substep value from the GPU side 10B andcomparing that substep value with the remaining time, following and responsive to the force integration and boundary handling subprocess 1 12; h. a convergence / di vergence check subprocess 116 performed by the CPU side 10A, following and responsive to the more substeps determination subprocess 114 determining that no further substeps are required in a GPU simulation loop 126 (also known as the "simulation step", which loops if more substeps are determined to be required in the more substeps determination subprocess 114); i. a convergence / divergence decision subprocess 1 18 performed by the CPU side 10A, following and responsive to the convergence / divergence check subprocess 1 16; j. a reactive movement correction subprocess 122 performed by the GPU side 10B (mostly) and by the CPU side 10A (the CPU side 10A only deals with initialization of the post-blast solid and a final projection step), following and responsive to the convergence / divergence decision subprocess 118 determining that the main simulation loop 120 has converged / diverged; and k. an output subprocess 124 performed by the CPU side 10 A, following and responsive to the reactive movement correction subprocess 122.Initialization subprocess 102Input data - blocks, pre-blast topography, floor, explosives, fragmentation

[0040] The initialization subprocess 102 commences with the CPU side 10A receiving preblast input data (also referred to as "simulation parameters") including a block-model 200 (including global and / or local mechanical values and explosive values, as described hereinafter), a pre-blast topography 202 and a floor mesh 204, e.g., as shown in FIG. 2. As shown in FIG. IB, the block-model 200, the pre-blast topography 202, the floor mesh 204, and a plurality of non-overlapping blocks 206 (forming the block-model 200) are transmitted from the storage 152 to the CPU 154 in the initialization subprocess 102. As shown in FIG. IB, in the initialization subprocess 102, the CPU 154 generates a GPUversion of the pre-blast topography 202 for each GPU 156 (in the form of one or more heightfields 202B of respective portions of the pre-blast topography 202 for the respective GPUs 156) and transmits the GPU versions to the GPUs 156A,156B (selected based on a split or partition of the simulated 3D space between the GPUs 156, as described hereinafter) and their respective VRAMs; similarly, the CPU 154 generates a GPU version of the floor mesh 204 for each GPU 156, in the form of a heightfield 204B (or portions thereof corresponding to respective portions of the floor mesh 204) and transmits the GPU versions to the GPUs 156A,156B (selected based on the split or partition of the simulated 3D space between the GPUs 156, as described hereinafter) and their respective VRAMs.

[0041] The block-model 200 represents the pre-blast solid material (i.e., the solid material before it is blasted), including 3D locations of the solid materials and one or more mechanical properties of the solid materials at those 3D locations. The block-model 200 includes the plurality of non-overlapping blocks 206 (e.g., cuboids) forming the blockmodel 200, e.g., as shown in FIG. 2. The block-model 200 is typically provided by a user, with attributes (grade, density, etc.) attached to each block 206 to represent the mechanical properties of the pre-blast solid material in the block 206, and 3D locations attached to each block 206 to define the 3D locations of the pre-blast solid material. The attributes attached to each block 206 include mechanical property values that are used to simulate / model mechanical movement of the solid material in the method 100, and these mechanical property values include, or represent, for the set of 3D locations in the solid materials: (i) a corresponding set of Poisson coefficient values, and a corresponding set of Young modulus values for co-rotated linear elastic constitutive modelling (such that each 3D point in the block-model 200 has a defined Poisson coefficient value and Young modulus value); (ii) a corresponding set of tensile strength values for fragmentation modelling (such that each 3D point in the block-model 200 has a defined tensile strength value); and (iii) a corresponding set of the selected rock friction angles, thus friction angle coefficient values, and a corresponding set of friction hardening coefficient values for sand plasticity modelling (such that each 3D point in the block-model 200 has a defined friction angle coefficient value and friction hardening coefficient value). An example block data structure may include the attributes in the form of global parameters for all rocks involved in the block-model, e.g., the selected mechanical properties mentioned hereinbefore (thecoefficient multiplying the explosive strength in the JWL model, the rock friction angle in the rock plasticity model, the rock crack threshold for the eigenerosion model, and / or the swell deformation rate multiplier for the artificial swelling). Another example block data structure may include the attributes per-block, generated by a blast planning software package, e.g., with manual annotation of the block-model (and tracking through time), and / or generated automatically from geological surveys, e.g., resulting from the material extracted from drill holes.

[0042] The pre-blast topography 202 is a plurality of 3D topography locations representing a pre-blast upper surface of the solid material, e.g., in the form of a 3D topography mesh 302 as shown in FIG. 3. The pre-blast topography 202 is measured from an actual surface of the solid material, and may be obtained from a photogrammetry system, e.g., a prior-art photogrammetry system, e.g., from PROPELLER AERO. The pre-blast topography 202 is represented in a data structure, generated by the CPU side 10A, including a topographic heightfield 202B, which allows for efficient collision-detection on the GPU side 10B. The GPU 156 accesses and uses the heightfield data structure 202B instead of the mesh data structure 202 (from which the heightfield 202B is determined) during the GPU simulation loop 126. The topography heightfield 202B is initialized on the CPU side 10A, and then uploaded to the VRAM of the GPU side 10B so the GPU 156 can access it. Example preblast topography data structures includes meshes sampled every 2 meters, with the extents depending on the extend of the blast, e.g., between 200m x 200m, and 500m x 500m.

[0043] The floor mesh 204 (also referred to as the "floor") defines as an unbreakable static topography, and includes a plurality of 3D locations representing a pre-blast lower surface of the solid material. This floor mesh 204 may be computed from the pre-blast topography, where all triangles inside of a user-defined perimeter are set at a given depth value. The floor mesh 204 is represented in a data structure, generated by the CPU side 10A, that includes the heightfield 204B, which allows for efficient collision-detection on the GPU side 10B. The GPU 156 accesses and uses the heightfield data structure 204B instead of the mesh data structure 204 (from which the heightfield 204B determined) during the GPU simulation loop 126. The heightfield 204B is initialized on the CPU side 10A, and then uploaded to the VRAM of the GPU side 10B so the GPU 156 can access it. The GPU 156 accesses the heightfield data structure 204B instead of the mesh data structure 202 fromwhich the heightfield 204B is computed. The floor mesh 204 may be selected to correspond to a vertical limit of the explosive material, or may be selected to correspond to a vertical level of a subsequent bench (e.g., when pre-conditioning simulation is needed). The floor mesh data structure, and thus the floor heightfield 20A, can have the same extents / precision of the pre-blast topography data structure.

[0044] The pre-blast input data also represent a selected fragmentation threshold value (e g., selected as a global coefficient, or selected per-block) that is selected based on the desired ease of fragmentation, e.g., chosen manually for each mine or material type based on measurements of pre-existing blasts.

[0045] The pre-blast input data may include a user-defined volume representing prefragmented material, e.g., a selected pre-fragmented volume 902 (e.g., manually drawn by the user, or identified as fragmented by a previous simulation run of the method 100), as shown in FIG. 9. All discrete blocks (and their corresponding simulation particles) below the surface of the selected pre-fragmented volume 902 (down to the floor 204) are marked / recorded as pre-fragmented, accounting for a previous blast that occurred on the front of this one

[0046] The pre-blast input data also represent an explosive blast plan for the solid material, including with a set of blast times (e g., a delay time for each explosive volume 804 in the blast plan), corresponding 3D blast locations in the pre-blast solid material, and corresponding explosive attributes for each explosive material type. The blast plan typically includes a plurality of mutually separate blast locations, each representing a blasthole with explosive material therein, e.g., holes 802A - 802E in FIG. 8 A. The separate blast locations are in volumes ("explosive volumes", or "blasting volumes"), e.g., substantially cylindrical, in the defined volume between the pre-blast topography 202 and the floor 204. The explosive volumes 804 can be substantially cylindrical volumes in the blastholes (e.g., as defined in the blast plan). In other words, an explosive blast can include a plurality of explosive blastholes 802A, 802B, 802C, 802D and 802E, and each explosive volume 804 can represents the explosive material (e.g., bulk explosive material, such as ANFO or ANE) loaded into the blast holes 802A, 802B, 802C, 802D and 802E. The explosive attributes attached to each explosive volume 804 include an initial internalenergy value represented by internal energy per unit initial volume (e.g., substantially 2 GJ / m3) and (in some embodiments) explosive mechanical property values to simulate / model mechanical movement of the explosive material in the method 100. The explosive attributes are represented in the input data. The mechanical values of the explosive material may include the same mechanical attributes associated with the blockmodel 200 (including Poisson coefficient values, Young modulus values); however, unlike the solid material, the explosive material does not typically require tensile strength values or friction coefficients for fragmentation modelling because the explosive material is typically represented by fluid behavior in the blast hole. The explosive blast plan data structure can have the same extents / precision of the pre-blast topography data structure, although the blast plan extent will generally be smaller than the topography extent (to leave sufficient room for the particles to spill on if they get pushed on a non-discretized area). The precisions are generally in the same order of magnitude; however, the precision of the topography meshes can be completely independent from the blast precision as far as the main simulation loop 120 is concerned. x' tension by null blocks

[0047] Following the receiving of the pre-blast input data, the CPU side 10A automatically extends the block-model 200 to fill the defined volume defined by the pre-blast topography 202 and the floor mesh 204 to form an extended block-model 400. In this process, as shown in FIGs. 4 and 5, the block-model 200 is extended automatically by null blocks 402, which are blocks (e.g., cuboids) with assumed attributes defining mechanical properties, e.g., representing assumed rock types where no geological information has been provided in the pre-blast topography 202 (the null blocks 402 are not air). The mechanical properties of the null blocks 402 can be global values, e.g., set to have the same values as other ones of the blocks in the block-model 200. The null blocks 402 can have the same formats, sizes and shapes as the other blocks and as each other, or different sizes / shapes, because all blocks in the method 100 (including the null blocks 402 and the block-model 200 blocks) can have mutually different selected sizes, within the spatial limits of the floor 204 and the pre-blast topography 202, so long as no blocks are mutually overlapping: thus the null blocks 402 may be sized adaptively, i.e., the sizes of the null blocks 402 can be selected for convenience since they will be broken down to smaller particles in thesubsequent discretization substep; e.g., the null blocks 402 may be generated down to a user-defined threshold (e.g., Ixlxl meter), and sized adaptively, including being smaller along the edges, and bigger if there is enough room between features of the pre-blast topography 202. As shown in FIG. 4, the automatic extending fills gaps inside the defined volume that is defined by, and between, the floor mesh 204 and the pre-blast topography 202. In the example of FIG. 5, example null blocks 502 substantially or completely surround the example blocks 504,506 to form an example extended block-model 500.

[0048] As shown in FIG. 4, the null blocks 402 at least fill the defined volume and may extend beyond the defined volume because they are cuboids: using cuboids for the blockmodel 200 and the null blocks 402 allows the extended block-model 400 to be composed entirely or substantially of un-rotated, non-overlapping cuboids, and this allows for efficient discretization in the following phase of the initialization subprocess 102.

[0049] As shown in FIG. IB, the automatic generation of the null blocks 402 (and the example null blocks 502) and formation of the extended block-model 400 (and the example extended block-model 500 and the example blocks 504,506) with the null blocks 402 (and the example null blocks 502) is performed by and in the CPU 154 and its RAM in the initialization subprocess 102.Discretization of (extended) block model

[0050] Following the formation of the extended block-model 400 (if required), the CPU side 10A automatically forms a discretized block-model 600, as shown in FIG. 6, by discretizing the extended block-model 400 (or by discretizing the block-model 200 if no extension by the null blocks 502 was necessary, which could be rare).

[0051] The discretizing of the extended block-model 400 (or the block-model 200) includes subdividing the extended block-model 400 (or the block-model 200) into a set of discrete particles, which may also referred to as "discretization points" or "material points". The discretized block-model 600 thus fills (or substantially fills) the defined volume, between the topology 202 and the floor 204, with the discrete particles. As shown in FIG.6, the set of discrete particles may extend beyond the defined volume because the blocks are non-overlapping cuboids, as mentioned hereinbefore.

[0052] The discretizing of the extended block-model 400 (or the block-model 200) includes the CPU side 10A using a selected particle size that defines a size for the discrete particles. The selected particle size may be automatically selected by the CPU side 10A based on input data representing user-defined parameters for the discretization; these user- defined parameters includes selected values for a simulation cell-width, and a desired number of particles per cell-width, thus the size for the discrete particles can be automatically selected by dividing the simulation cell-width by the desired number of particles per cell-width. The cell-width is one of the main parameters affecting accuracy and computation times. Smaller cell-width means more accuracy, but higher computation times. The cell-width may be selected to be as small as possible based on a pre-selected time budget, which depends on the blasting operation. Example cell-width values include substantially 1 m and substantially 2 m, and values between substantially 1 m and substantially 2 m, or between substantially 0. I m and substantially 10 m in some applications. Fractional sizes of the cell-width are addressed by creating smaller particles, e.g., if 2.5 cells fit into a block, instead 3 particles are fit into the block with the discretized block size divided by 3 as diameter. Accordingly, the discrete particles in the discretized block-model 600 can include two or more mutually different diameters.|0053| As shown in FIG. 7, an example discretized extended block-model 700 includes a set of example discretized null blocks 702 surrounding a set of example discretized blocks 704, 706. The example discretized blocks 704, 706 have the same area / volume as the respective example blocks 504,506 with less blocky boundaries. Compared to FIG. 5, FIG. 7 shows a greater area of the discretized blocks 704, 706 because the particle rendering allows for seeing through small gaps between adjacent points (since they are not rendered with their real simulated radius). So, in FIG. 7, some of the example discretized null blocks 702 are not visible because they are just a thin layer on top of the discretized blocks 704, 706

[0054] Each of the discrete particles includes mechanical property values that correspond to or are equal to the mechanical property values of the block-model 200 at that particle's 3D location. In other words, the set of discretized particles includes or represents: (i) a corresponding set of Poisson coefficient values and a corresponding set of Young modulus values (such that each 3D point in the discretized block-model 600 has a defined Poissoncoefficient value and Young modulus value); (ii) a corresponding set of tensile strength values for fragmentation modelling (such that each 3D point in the discretized block-model 600 has a defined tensile strength value); and (iii) a corresponding set of friction angle coefficient values and a corresponding set of friction hardening coefficient values (such that each 3D point in the discretized block-model 600 has a defined friction angle coefficient value and friction hardening coefficient value). In some implementations, some of these values are assigned to the mechanical model bound to the particle instead of the particle itself, in other examples, these values are assigned to the particle itself.

[0055] As shown in FIG. IB, in the initialization subprocess 102, the discretization of the extended block-model 400 to form the discretized block-model 600 occurs in the CPU 154 and its RAM in the initialization subprocess 102, and the CPU 154 transmits the discretized block-model 600 (or the example discretized extended block-model 700, including its discretized null blocks 702 and its discretized blocks 704, 706), or split portions thereof, to the one or more respective the GPUs 156A,156B (the split portion of the discretized block-model 600 being sent to the GPU assigned to that partition of the simulated 3D space, as described hereinafter) and their respective VRAMs.Discretization of explosive volumes 804

[0056] As shown in FIG. 8A, the initialization subprocess 102 automatically discretizes each of the explosive volumes 804 by defining a set of explosive particles (also referred to herein as "blast points") for each 3D blast location. The generated set of 3D locations of the explosive particles in each blast location can be generated to form a uniform distribution, or a random or quasi-random uniform distribution, in the or each explosive volume 804, thus substantially filling the explosive volume 804 with the explosive particles. The filling method may include the initialization subprocess 102 dividing the or each blast volume 804 into a plurality of pieces, e.g., NZ = H / (2R) cylindrical pieces, where cylinder’s height is H, and N explosive particles with radiuses R are required, then discarding any particle outside the shape of the blast volume 804 until N / NZ particles are preserved on each piece.

[0057] The explosive particles each include their 3D location and an (initial) internal energy value (from the initial internal energy of the corresponding explosive material type), e.g., the selected coefficient multiplying the explosive strength. The explosive particles also each include a plurality of specific blasting parameters (or "JWL" parameter, described hereinafter), which include: A, B, and omega (which are three dimensionless linear coefficients) and Rl, R2 (which are two dimensionless nonlinear coefficients).

[0058] The explosive particles each include (which may include having in its data record, or being uniquely associated with, e.g., by a data pointer): (i) an ignition time that is based on, and / or equals, the delay time (or blast time) for the corresponding explosive volume 804 in the blast plan, at which time the corresponding explosive particles are inserted into the simulation loop as gas particles, and (ii) a lifetime that represents how long after the ignition time the gas particles of that explosive volume 804 remain in the simulation loop.

[0059] As mentioned hereinbefore, in some embodiments, the explosive particles may each include the explosive mechanical property values (mentioned hereinbefore) that are used for elastic constitutive modelling and / or sand plasticity modelling. The mechanical values for the elastic constitutive model properties of the explosive materials define the explosive material properties in their solid state. In other embodiments, the solid elastic explosive material before and during ignition are not simulated. Whether to model the solid elastic explosive material before and during ignition is selected dependent on the spatial accuracy of the method 100: if the cell-width is an order of magnitude smaller than the drill-hole’s diameter, deformation of the solid elastic explosive material may be simulated. If, on the other hand, when the selected grid-cell size (e.g., 1 to 2 m for subminute simulations) is an order of magnitude larger than the explosive’s diameter (e g., the diameter of the explosive volume 804), the absence of rock inside the explosive volume 804 need not be simulated, so the explosive particles are inserted in between the other discretized particles, e.g., as shown in FIG. 8B. Unlike the other discretized particles (representing the block-model 200), the explosive particles need not require tensile strength values for fragmentation modelling, e.g., representing fluid explosive material types in contract to relatively brittle rock types.

[0060] The explosive particles each include (which may include having in its data record, or being uniquely associated with, e.g., by a data pointer): (i) an ignition time that is based on, and / or equals, the delay time (or blast time) for the corresponding explosive volume 804 in the blast plan, as which time the corresponding explosive particles are inserted into the simulation loop as gas particles, and (ii) a lifetime that represents how long after the ignition time the gas particles of that explosive volume 804 remain in the simulation loop.

[0061] In some examples, the CPU 154 instantiates the explosive particles based on the JWL model for sending to the GPU 156, including using a plurality of JWL models if there is a plurality of explosive types For example, as shown in Code Appendix N (abbreviated hereinafter to "Apx N"), e.g., the CPU 154 can call the 'JWLExplosive::constitutive_moder function, and for each explosive type, a’ : constituti ve model can be instantiated, and added to a list that is then sent to the GPU 156.Upload to GPU side 10B

[0062] As shown in FIG. 1A, in the initialization subprocess 102, after the CPU side 10A forms the discretized block-model 600 (e.g., 700, comprising "rock particles"), the CPU side 10A sends the discretized block-model 600 to the GPU side 10B, which in turn receives the discretized block-model 600 for use in the GPU simulation loop 126.

[0063] As shown in FIG. IB, in the main simulation loop 120, after the CPU side 10A forms the or each discretized explosive volume 804 (e.g., comprising "explosive particles"), the CPU side 10A sends the or each discretized explosive volume 804 to the GPU side 10B, which in turn receives the or each discretized explosive volume 804 for use in the GPU simulation loop 126.

[0064] For the multi-GPU case, the simulated 3D space is split into as many partitions as there are GPUs 156. The split is operated in such a way that each of these partitions contains roughly the same number of rock particles, thus each partition is associated to one of the GPUs 156. The set of particles is split into two or more subsets, with one subset per partition (thus two particles located in the same partition are put in the same particle subset). Each particle subset is then uploaded to the corresponding GPU 156. Whenexplosive particles are generated (which includes adding them as the gas particles), their positions are checked against the partitions to determine on which one of the GPUs 156 the explosive particles are to be uploaded. Thus one particle is only ever present on a single GPU at a time, and is never moved from one GPU to the other, even if the particle's simulated trajectory is such that it moves outside of the partition associated to this GPU (since the partitioning system is only used for the first time a particle is initialized / uploaded to a GPU).Main simulation loop 120

[0065] As shown in FIG. 1 , the main simulation loop 120 includes: the update explosives subprocess 104 (which starts the main simulation loop 120), the GPU simulation loop 126, the convergence / divergence check subprocess 1 16, ad the convergence / divergence decision subprocess 118.

[0066] The main simulation loop 120 repeatedly performs the update explosives subprocess 104 and the GPU simulation loop 126 until the method 100 determines a main loop end condition, specifically by the convergence / divergence check subprocess 116 detennining that either the simulation has diverged (which means the simulation has either entered an unrecoverable state, e.g., NaN in the simulation, or has entered obviously undesirable state, e.g., if some particles are sent flying very far in the air, or outside of the reasonable perimeter) or that the simulation has converged (which means the simulation has reached an equilibrium that will no longer make any notable progress, e.g., if all particles have nearly stopped moving).

[0067] As shown in FIG. 1C, simulation outputs from the output subprocess 124, including replay files and statistics, are transmitted from the CPU 154 to the storage 152 after the main simulation loop 120.Update explosives subprocess 104

[0068] As shown in Code Appendix A (abbreviated hereinafter to "Apx A"), which includes pseudocode implementing the update explosives subprocess 104, the update explosives subprocess 104 includes the CPU side 10 A:a. starting a simulated time at commencement of the main simulation loop120, e.g., starting at 0 seconds — all hole timings may be shifted so that the first hole time is almost zero to avoid spend a long time simulating a still model before the first explosion happens; b. for each of the explosive particles that are still live, i.e., explosive particles that have been previously added to the main simulation loop 120 as respective gas particles, and have not yet expired and been removed, wherein the live explosive / gas particles represent the remaining explosive forces at the current time of the main simulation loop 120 (e g., see Apx A, line 1 — abbreviated hereinafter to "Apx A:l"): i . determ i n e i f th e si mul ated ti m e i s n ow greater th an th e expi rati on time for the gas particles (wherein each expiration time is the sum of the explosive particle's ignition time and the explosive particle's lifetime), i.e., determine if the explosive particle has been live (thus in its gas form) in the main simulation loop 120 for a duration greater than its lifetime (e g., see Apx A:2), and ii. if the simulated time is greater than the expiration time of that explosive particle, then remove that explosive particle from remainder of the main simulation loop 120, e.g., by setting the status to not alive (e.g., see Apx A:3), e.g., by deleting each of these expired particles from the GPU memory — in other words, the explosive particles detected as non-live are removed from the simulation memory, they no longer exist in VRAM, and are no longer part of the lists of active or waiting explosives; c. for each of the explosive particles that are still to be initiated since commencement of the main simulation loop 120, that is, for the explosive particles that have not yet been inserted as gas particles (e.g., see Apx A:4): i. determine if the simulated time matches an ignition time of the explosive particle (thus the hole timing controls when the explosiveparticles are added to the simulation) (e.g., see Apx A:5) — this matching is performed in-between two simulation steps, and a simulation step is at most equal to a user-defined value (e.g., 14ms in example implementations), or equal to the time before the next explosive liveness event (either an explosive must be removed (based on the lifetime expiring), or added (based on the ignition time occurring)), whichever is smaller, ii. if the simulated time matches the ignition time for that explosive particle, then generate a gas particle (also referred as a "JWL particle") substantially at the 3D location defined by the explosive particle (in other words, insert a particle — the gas particle — at that location, e.g., by initializing the material properties of the added gas particles using the or each JWL model for the or each explosive type, e.g., see the JWLExplosive::particle function in Apx N), thus generating many gas particles for the drill hole, as in FIG. 8B (e.g., see Apx A: 6), and iii. send the generated gas particles (in the explosive volume 804) to the GPU side 10B (e.g., see Apx A:7) in a memory transfer from the CPU 154 to the GPU(s) 156, as shown in FIG. 1C; and d. determine and set the next simulation step duration to be less than a lowest of the ignition times / timings of to-be-exploded explosive particles, representing the explosives that have not yet been initiated, and the next of the expiration times of the remaining particles, such that the duration of the next simulation step is less that the next time to add or remove the gas particles representing the explosive particles (e.g., see Apx A:8).

[0069] In FIG. 8, the blastholes are depicted in a side cut view of a blast, and the rightmost hole 802E has reached its ignition time and has thus been discretized into the gas particles, whereas the holes on the left 802A - 802D haven’t reached their ignition times so no gas particle have been added for them yet. The discretized particles (e g., rock)overlapping with the initiated blasthole 802E are retained in place for the GPU simulation loop 126, thus overlapping rock particles are not removed

[0070] The explosive volumes 804 generally have a plurality of mutually different ignition times / delay times defined in the blast plan. Everything in the simulation loop is affected by the same simulation step / substep lengths (or "timestep"), which may varies depending on the explosive lifetimes or the next ignition time.

[0071] In the update explosives subprocess 104, for each explosive blast time (or "delay time"), the CPU side 10A includes (or inserts) the explosive particles as gas particles with (or into) the set of discrete particles for the GPU simulation loop 126 to operate (which includes, as described in more detail hereinafter, numerically determining the velocities and updated 3D locations of the discrete particles, and numerically determining / calculating the updated velocities and the updated 3D locations at time steps (each with a time step duration) at and after the or each explosive blast time using the initial internal energies of the included / inserted explosive particles). Thus, the numerical determination of the explosive particles' velocities or equivalently forces is done at least at each explosive blast time among the plurality of time steps. The explosive ignition times and lifetimes control when the explosive particles are inserted or removed from the simulation loop. As long as an explosive particle is involved in the simulation, its velocity, position, etc. are naturally updated at the same time, in the MPM pipeline (described hereinafter), as the other (rock) particles involved in the simulation.

[0072] With reference to Apx A: 5-6, the gas particles are added to the simulation only at the time indicated by the corresponding blasthole’s timing because the update explosives subprocess 104 generates and inserts the gas particle substantially at the 3D locations of the explosive / gas particles in the explosive volume 804 when the simulated time is substantially equal to the ignition time of that explosive volume 804.

[0073] With reference to Apx A: 3, the explosive / gas particles are removed when the simulated time is equal to the lifetime of their explosive volume 804; in other words, the update explosives subprocess 104 models the transient nature of each explosive blast by removing the explosive particles from the discrete particles after the lifetime. The lifetimeis how long the gas particle exist in the simulation, after which they are removed, thus the gas particles are always active in the simulation loop or deleted from the relevant GPU's memory (the VRAM).

[0074] The lifetime may be a pre-defined duration, set for all the explosives involved in the blast. The deletion of the gas particles after the lifetime effectively models the explosive products of the initiation as gaseous products that no longer contribute to the simulated volume, e.g., the explosive products may be regarded as being released as gases into the air through cracks in the solid material. In embodiments with a large cell-width (e g., when the selected grid-cell size is an order of magnitude larger than the explosive’s diameter, mentioned hereinbefore), the method 100 effectively removes the explosive particles from the simulation after the pre-defined lifetime duration to simulate their release into the air through the rock cracks; for example, a cell-width of, say, 1 meter would require a crack of at least 3 meters for the gas to escape freely, which may represent an unrealistic amount of movement and fragmentation in typical operations.

[0075] The gas particles may be simulated in the particle and grid update subprocess 110, FIG 1 according to a Jones-Wilkins-Lee (JWL) model that defines mathematical relationships between gas particle pressures, the internal energies and particle stresses, allowing determination of new particle velocities after integration of forces, e.g., as described in Xiong Zhang et al, "The Material Point Method: A Continuum-Based Particle Method for Extreme Loading Cases" Elsevier Inc. (978-0-12-407716-4), Chapter 6.

[0076] The particle and grid update subprocess 110 may include exemplary functions to solve pressure and energy and to update internal energy and pressure, e.g., as shown in Apx H. The resulting stress tensor may be determined using an exemplary function to determine the Kirchhoff stress, e.g., as shown in Apx H.

[0077] In the update explosives subprocess 104, the CPU side 10A automatically determines what GPU time step duration (or "substep length") to use in the next for the GPU simulation loop 126, e.g., to be between substantially 10 microsecond (ms) and 1000 ms, e g., substantially 100 ms. As typical blasting operations involve extremely high velocities (due to the explosive gas particles), as well as very high material stiffnesses (dueto rocks) and geometric discontinuities (due to fractures), generally extremely small time steps (around l e4seconds) are used approach satisfying Courant-Friedrichs-Lewy (CFL) stability conditions. The use of the extremely small time steps (around le4seconds) allows use of explicit integration in the force integration and boundary handling subprocess 112 that is significantly more computationally efficient per time step than an implicit scheme that would involve the expensive resolution of a complex sparse linear system The determined GPU time step duration is determined based on: the particle velocities, rock stiffness, and explosive parameters (when explosive particles present), and a relationship representing the CFL stability conditions.

[0078] In exemplary examples, the update explosives subprocess 104 may be implemented using functions from G1THUB, e.g.: a. the timestep limit for each particle may be computed using the following routine: https: / / github.eom / dimforge / sparkl / blob / v0.2.l / src_kemels / cuda / timestep.rs #L81-L112 b. the call to "estimate_particle_timestep_length" may implement modelspecific timestep limits, including for rock particles using the following functions: http s : / / github . com / dimforge / sparkl / bl ob / vO .2.1 / src_core / dy nami cs / model s / e lasticity corotated linear. rs#L 105-L 113 http s : / / github . com / dimforge / sparkl / bl ob / vO .2.1 / src_core / dy nami cs / timestep / elasticity sound speed timestep bound. rs#L28-L33

[0079] For explosive particles, the update explosives subprocess 104 may be implemented using functions timestep bound and adiabatic sound speed in Apx H.GPU simulation loop 126

[0080] The GPU simulation loop 126 includes repeatedly, for a plurality of substeps (each with the determined GPU time step duration) performing (or "executing") the following subprocesses, as shown in FIG. 1 : the update sparse grid subprocess 106, the estimatesubstep length subprocess 108, the particle and grid update subprocess 110, the force integration and boundary handling subprocess 1 12, and the more substeps determination subprocess 114.

[0081] The GPU simulation loop 126 repeatedly performs its subprocesses until the method 100 determines a sub loop end condition, specifically by the more substeps determination subprocess 114 determining that no more substeps are required in the current main simulation loop 120.

[0082] The GPU simulation loop 126 includes the GPU side 10B making CUDA API calls to the CPU side 10A for orchestrating execution of GPU kernels in the GPU side 10B, e.g., for all of the subprocesses in the GPU simulation loop 126.Update sparse grid subprocess 106

[0083] The update sparse grid subprocess 106 generates the "grid" onto which the particle values are mapped. The update sparse grid subprocess 106 may be understood with reference to Xinlei Wang et al., "A Massively Parallel and Scalable Multi-GPU Material Point Method" in ACM Trans. Graph. (SIGGRAPH), Vol. 39, No. 4, Article 30 (July 2020), and Ming Gao et al., "GPU Optimization of Material Point Methods" in ACM Trans. Graph. (SIGGRAPH), Vol. 37, No. 6, Article 254 (November 2018).

[0084] As shown in Apx B, which discloses pseudocode implementing the update sparse grid subprocess 106, the update sparse grid subprocess 106 includes the GPU side 10B: a. identifying 3D chunks by iterating on each particle, dividing that particle’s position by the cell-width, rounding, and subtracting 2, which gives the cell / particle association — each cell coordinate is divided by 4 to identify the chunk., e.g., cubes, of cells, e.g., 4x4x4 cubes of cells, wherein each chunk contains at least one particle — multiple particles can be mapped to the same cell (and thus the same chunk as well) — the grid-cell size is the size of each cell — and the grid-cell size is a global parameter of the whole simulation, selected by the user (e.g., see Apx B:l);b. resizing a vector of the 3D chunks based on the required number of chunks in order to update the size of the GPU buffer that will contain all the 3D chunks required — in this context, “vector” means GPU buffer (i.e., not a linear algebra “vector”) — this is done to provide the right amount of memory available for the chunks involved in this GPU step (e.g., see Apx B:2) in a memory allocation, as shown in FIG. 1C; c. sorting particle indices based on their containing chunk indices (e g., see Apx B:3); d. allocating the chunks and populating a hash-map as an inventory of the existing chunks — the hash-map is a data structure that maps keys (here, chunk coordinates) to some data (here, a chunk pointer to the chunks buffer), and the inventory is populated so that (1) during the particle-to- chunk mapping, the same chunk is not created twice, and (2) the location of a chunk in memory can be easily located based on its coordinates, e g., needed by the G2P2G as described hereinafter (e.g., see Apx B:4) in a memory allocation from the CPU 154 to the GPU(s) 156, as shown in FIG. 1C; and e. if there is a plurality of GPU's (forming a multi -GPU) in the GPU side 10B, identifying halo-blocks as blocks in the radii of effect particles originating from the two GPUs — after each GPU updates its sparse grid based on the particles it is responsible for simulating, the GPU will have a known set of blocks, and each GPU will then gather the block coordinates of all the blocks existing on the other GPUs, and any block coordinate a given GPU has in common with the block coordinates gathered from the other GPUs is marked as halo-block (e.g., see Apx B:6).

[0085] In exemplary examples, the update sparse grid subprocess 106 may be implemented using functions available from G1THUB, e.g.: a. https: / / github.eom / dimforge / sparkl / blob / v0.2.l / src_kernels / cuda / grid_updat e.rsb. https: / / github.eom / dimforge / sparkl / blob / vO.2.l / src_kemels / cuda / prefix_su m.rs c . http s : / / gi thub . com / dimforge / sparkl / bl ob / vO .2.1 / src_kemel s / cuda / sort.rs

[0086] The update sparse grid subprocess 106 includes defining a 3D grid (also referred to as a "topologically invariant sparse grid" or a "grid") substantially / wholly surrounding the discrete particles. The grid surrounds the particles, plus 2 empty grid cells in every direction. The grid is composed of regular cells, each with a width referred to as the "cellwidth" or "grid-cell size". The grid-cell size defines the grid's grid spacing. As mentioned hereinbefore, grid-cell size may selected to be within an order or magnitude of the particle size, e.g., between 0.1 meters and 10 meters in each direction, e.g., substantially 1 meter. This grid-cell size controls a spatial resolution of the main simulation loop 120 by controlling a radius of force propagation from one particle to the other particles within the simulated continuous material because of the compact support described hereinafter: for example, if the kernel function has compact support (e g., is a quadratic b-spline) over a selected number of grid cells (e.g., 3x3x3), then any ones of the particles that are closer than the selected number of grid cells (e g., 3) have their velocities and locations updated due to mechanical properties and velocities of the surrounding particles within that selected number of grid cells, thus are treated by the particle and grid update subprocess 110 as being part of the same continuous material, with cohesive forces growing stronger as the particles get closer to one another. The method 100 avoids direct particle-particle interaction because the particle mechanical quantities get transferred to the grid, and then the grid mechanical properties (updated velocities) are transferred to the particles, thus the grid is always the intermediary between two particles interacting.Estimate substep length subprocess 108

[0087] As shown in Apx C, which discloses pseudocode implementing the estimate substep length subprocess 108, the estimate substep length subprocess 108 includes the GPU side 10B: a. setting or re-setting the substep to be equal to the simulation step length, which initializes the loop — the simulation step length is the user-decidedmax acceptable simulation step, and the rest of the substep length is iteratively updated from this base value (by always taking the minimum with other calculated substep length values in the substep length calculation loop) (e.g., see Apx C:l); b. for each particle in the model, thus all particles, including any explosive particles (i.e., the inserted gas particles) from the update explosions subprocess 104 (e g., see Apx C:2): i. setting or re-setting the substep to be equal to the lessor of: the substep, and the cell-width divided by the particle-velocity because, within a substep, a particle should not travel by more than one cellwidth to preserve stability (e.g., see Apx C:3); and ii. setting or re-setting the substep to be equal to the lessor of: the substep, and a value based on the particle’s constitutive model to satisfy the CFL stability conditions, because CFL stability conditions keep the numerical stability of the simulation loop in check for explicit integrators (without the limitation on substep length, the simulation may blow up due to the numerical integrations diverging) (e g., see Apx C:4); and c. reading the minimum substep length from the GPU side 10B (e.g., see Apx C:5), to the CPU side 10A, in a memory transfer from the GPU(s) 156 to the CPU 154, as shown in FIG. 1C.

[0088] The estimate substep length subprocess 108 is used because the numerical integration involves the approximations (linearizations) of the various non-linear equations in the model, and larger substeps result in larger approximations, resulting in potential over-estimation of forces, resulting in divergence.Particle and grid update subprocess 110

[0089] The particle and grid update subprocess 1 10 is repeated for each iteration of the GPU simulation loop 126 in the GPU side 10B.

[0090] The particle and grid update subprocess 110 uses the discretization particles that hold / access the mechanical property values of the simulated solid material, thus representing the solid material as the mechanical property values in the 3D locations of the discretization particles.

[0091] As shown in Apx D, which discloses pseudocode implementing the particle and grid update subprocess 110, the particle and grid update subprocess 110 includes the GPU side 10B: a. first, transferring the grid values to the particle values (referred to as "G2P") using an MLS approximation, and calculating new / update values (e.g., see https: / / github.eom / dimforge / sparkl / blob / v0.2.l / src_kemels / cuda / g2p2g.rs# L353-L474, and https: / / github.eom / dimforge / sparkl / blob / v0.2.l / src_kernels / cuda / g2p2g.rs# L193-L245), including: i. transferring the velocities (velocity values or velocity gradient values) from the grid cells to the particles (e.g., see Apx D:2, and http s : / / github . com / dimforge / sparkl / bl ob / vO .2.1 / src_kernel s / cuda / g2p 2g.rs#L2 1, and Apx K, line below “Update particle velocity”) — wherein each particle velocity is forcefully capped to a selected value (e.g., cell_width / substep time', see Apx K, block below “Cap particle velocity”) to be 100% guaranteed that the particle does not move by more than one grid cell, ii. integrating the particle velocity to obtain its new position (e.g., see Apx K, line below “Integrate position”), iii. calculating the particle velocity gradients from the grid velocities (e.g., see Apx D:3, and http s : / / github . com / dimforge / sparkl / bl ob / vO .2. 1 / src_kernel s / cuda / g2p 2g.rs#L232, and http s : / / github . com / dimforge / sparkl / bl ob / vO .2.1 / src_kernel s / cuda / g2p 2g.rs#L234), andiv. calculating averaged positive elastic strain energy from the cells to the particles (e.g., see Apx D:4, and http s : / / gi thub . com / dimforge / sparkl / bl ob / vO .2.1 / src_kernel s / cuda / g2p 2g.rs#L233). b. second — although only during the reactive movement correction subprocess 122 — calculating an artificial suction force based on a mass difference between the cell containing the particle and adjacent cells (e g., see Apx D:5, and https: / / github.eom / dimforge / sparkl / blob / v0.2.1 / src_kemels / cuda / g2p2g.rs# L238-L244), c. third, updating each particle in the simulation (thus all particles — non-live particles are removed as soon as they become non-live) by, for each particle (e.g., see Apx D:7, and https: / / github.eom / dimforge / sparkl / blob / v0.2.l / src_kemels / cuda / g2p2g.rs# L238-L244): i. updating its fragmentation status based on the particle elastic strain energy (based on the selected rock crack threshold and eigenerosion (e.g., see Apx D:8, and Apx K: block below the “Modified Eigenerosion” comment), ii. integrating its velocity gradient to update its deformation gradient (e.g., see Apx D:9, and Apx K: block below “Deformation gradient update” ), iii. calculating its particle stress using the constitutive model and deformation (e.g., see Apx D: 10, and Particle stress calculation in http s : / / github . com / dimforge / sparkl / bl ob / vO .2.1 / src_core / dy narni cs / models / elasticity_corotated_linear.rs#L28-L75, and https: / / github.eom / dimforge / sparkl / blob / v0.2. l / src_core / dynamics / models / elasticity corotated linear.rs#L28-L75),iv. updating its particle positive elastic strain energy (e.g., see Apx D: 1 1 , and Strain energy calculation in http s : / / gi thub . com / dimfbrge / sparkl / bl ob / vO .2.1 / src_core / dy nami cs / models / elasticity_corotated_linear.rs#L77-L90, andStrain energy update in Apx K, block below “Update Pos energy” ), v. only during the reactive movement correction subprocess 122, calculating the spring-like attraction force to the particle if it is outside the post-blast solid (e.g., see Apx D: 12, and Apx K, the block below “Artificial spring-like force”, in which there is a reversed concept of interior for the post-blast solid because the "GpuGridProj ectionStatus: Inside" state corresponds to a particle being outside of the post-blast solid); d. fourth, for each fragmented particle in the simulation loop (e.g., see Apx D: 13): i. projecting the particle stress based on the Drucker Prager plasticity model (e g., see Apx D: 14, and https: / / github.eom / dimforge / sparkl / blob / v0.2. l / src_core / dynamics / models / plasticity_drucker_prager.rs#L34-L95), and ii. applying artificial swell to the particle (e.g., see Apx D: 15, and Apx K, code below “Apply artificial swelling” ); and e. fifth, transferring the particle values to the grid values (referred to as "P2G") using the MLS approximation (e.g., see https: / / github.eom / dimforge / sparkl / blob / v0.2.l / src_kemels / cuda / g2p2g.rs# L476-L553, and https: / / github.eom / dimforge / sparkl / blob / v0.2.l / src_kemels / cuda / g2p2g.rs# L280-L350) including: i. transferring the momentum, the mass, and the forces from the particles to the grid cells (e.g., see Apx D: 16, andhtp s : / / github . com / dimforge / sparkl / bl ob / vO .2.1 / src_kernel s / cuda / g2p 2g.rs#L332-L333), and ii. transferring the positive elastic strain energy from the particles to the grid cells (e g., see Apx D: 17, and http s : / / github . com / dimforge / sparkl / bl ob / vO .2.1 / src_kernel s / cuda / g2p 2g.rs#L334-L335).

[0092] The example source code and pseudo code show the example variables are being used in each example calculation, e.g., see Apx H, I, J, L, M, N, and structure definitions on the open-source "sparkl" repository on GITHUB.

[0093] The particle and grid update subprocess 110 may be referred to as "G2P2G" because it includes the transferring the grid values to the particle values (referred to as "G2P") followed by the transferring the particle values to the grid values (referred to as "P2G"). By way of the G2P2G, the particle and grid update subprocess 110 numerically determines / calculates the updated particle velocities (which collectively form a "velocity field"), and updates the 3D locations of the particles based on their mechanical properties and their 3D locations in the 3D grid G2P2G is a more efficient approach than P2G2P, especially for the multi-GPU case.

[0094] The GPU side 10B maintains the 3D grid ("topologically invariant sparse grid") for the purpose of calculating a new velocity field from the particles’ mechanical properties and boundary conditions in each iteration of the GPU simulation loop 126, thus at each time step. The ephemeral nature of this grid (it is reset at each iteration of the GPU simulation loop 126), eliminates typical mesh inversion problems encountered in meshbased methods (like finite element analysis), avoiding expensive remeshing and numerical instabilities due to large deformations and dynamic topology changes (fractures). The GPU simulation loop 126 uses a double-buffering strategy, thus keeping two grids in memory: at substep N, grid A contains mechanical properties transferred from particles during the previous substep, and grid B is reset, which is then populated during G2P2G (which effectively reads from grid A and writes to grid B), then A and B are swapped. The 3D grid is composed of regular cells with a width equal to the cell-width. This cell-widthplays a key role in the spatial resolution of the simulation as it controls the radius of force propagation within the continuous material, as explained hereinbefore with reference to the update sparse grid subprocess 106. For example, because the compact support of quadratic b-splines covers 3x3x3 grid cells, particles that are closer than 3 cell-widths from one another will start being considered as being part of the same continuous material (with cohesive forces growing stronger as the particles get closer to one another). As a result, a compromise must be considered. A smaller cell-width results in a higher spatial resolution of the simulation: simulated cracks can be thinner, and movement is more detailed (less artificial stiffness), closer to ground truth. However, a small cell-width value results in more grid cells, and more particles to cover the same volume of the model, resulting in a significant increase in computational times (e.g., dividing the cell-width by 2 increase computation times by a factor of at least 8). In experimental embodiments, cell-widths of around 1 meter were used.

[0095] The particle and grid update subprocess 110 uses a kernel function, e g., see Kernel function implementation in https: / / github.com / dimforge / sparkl / blob / vO.2.1 / src_core / dynamics / solver / kernel.rs#L8- L136

[0096] In the particle and grid update subprocess 110, the determining of the velocities and the updated 3D locations includes using a quadratic B-spline function as the kernel function to combine the mechanical properties of ones of the discrete particles adjacent to each discrete particle. The particle-to-neighbor-particle velocity transfer is indirect, through the grid, and is essentially described by particle and grid update subprocess 110 (the G2P2G). The quadratic B-spline function is used to compute coefficients in G2P (e.g., see https: / / github.eom / dimforge / sparkl / blob / v0.2. l / src_kemels / cuda / g2p2g.rs#L224- L228) and P2G (e g., see https: / / github.eom / dimforge / sparkl / blob / v0.2. l / src_kernels / cuda / g2p2g.rs#L312-L316).

[0097] In the particle and grid update subprocess 110, the discretization operations use a moving least squares (MLS) approximation, e.g., a Galerkin-style moving least squares (MLS) discretization, e.g., as described in Yuanming Hu et al., "A Moving Least Squares Material Point Method with Displacement Discontinuity and Two-Way Rigid BodyCoupling" in ACM Trans. Graph. (SIGGRAPH), Vol. 37, No. 4, Article 150 (August 2018), or in Yuanming Hu's GTTHUB repository: https: / / github.com / yuanming- hu / taichi_mpm.Fragmentation

[0098] To model fragmentation of the solid material, such as internal failure / ripping of the rock portions due to the explosive blasting, the particle and grid update subprocess 110 includes the GPU side 10B: a. determining a peak effective energy release rate for each discrete particle to determine a fragmentation status for that discrete particle, wherein the fragmentation status is determined to be fragmented if the peak effective energy release rate for that discrete particle is above the selected fragmentation threshold value, and b. numerically determining / calculating the velocities and the updated 3D locations of the discrete particles based on their respective fragmentation statuses, including reducing / zeroing the tensile strength value in the mechanical properties to model fluid-like behavior when the fragmentation status represents fragmented.

[0099] The peak effective energy release rate for each discrete particle is calculated from a normalized weighted average of the positive (extension) elastic strain energy density of all the particles within a neighborhood of the discrete particle being processed: these elastic strain energies (e.g., see https: / / github.eom / dimforge / sparkl / blob / v0.2.l / src_core / dynamics / models / elasticity_corot ated_linear.rs#L79-L83) obtained from the particle’s deformation gradient (e g., see Apx K - lines below “Deformation gradient update” in the ! is_fluid case), e.g., based on the eigenerosion method described in Kun Zhang et al., "Dynamic brittle fracture with eigenerosion enhanced material point method" in International Journal for Numerical Methods in Engineering, Volume 121, Issue 17 p. 3768-3794. The selected fragmentation threshold values are an input to the simulation, as mentioned hereinbefore.

[0100] The particle and grid update subprocess 110 also determines the fragmentation status of one or more of the discrete particles in the pre-blast solid material based on a fragmentation input, e.g., a user-defined volume representing pre-fragmented material that are affected by the eigenerosion model (and the sand plasticity model mentioned hereinafter) in the first numerical determination / calculation (substep) to account for previous blasts that would have already broken part of the solid material.

[0101] Thus the discrete particles may be in the fragmented state (represented by a data flag) because they are in the user-defined volume (e.g., 902) representing pre-fragmented material, or because they are changed to fragmented by the particle and grid update subprocess 110 (e.g., Apx D:8).

[0102] To model the fragmented discrete particles, the particle and grid update subprocess 110 includes the GPU side 10B determining / calculating the velocities and updated 3D locations using the friction angle and the friction hardening when the fragmentation status represents fragmented (e g., see lines 7 and 8 of Apx D), thus representing the fragmented particles using a sand plasticity model, e.g., a Drucker-Prager sand plasticity method, e.g., as described in Gergely Klar et al., "Drucker-prager elastoplasticity for sand animation", ACM Trans. Graph. (SIGGRAPH), Vol. 35, No. 4, Article No. 103 (pp 1-12). The sand plasticity model eliminates the rocks' tensile strength but preserves friction, which improves on the "fluid-like" behavior represented merely by the tensile strength going to zero.Swelling

[0103] To model the swelling, the particle and grid update subprocess 110 includes increasing a particle volume represented by each discrete particle to represent swelling of the solid material due to the explosive blasts. This increasing of the particle volume can be referred to as "artificial swelling". Real-life blasts results in irregularly shaped rocks creating air gaps as they roll out of their tightly packed initial positions. The method 100 uses relatively large cell-widths, so the rock fragments are typically 1 -particle wide with an implicit ellipsoid shape due to anisotropic deformations, and this does not provide enough geometric precision to define air gaps between rocks; however, these air gaps generallyresult in a significant increase in the post-blast muck pile volume, so, in the particle and grid update subprocess 1 10, the GPU side 10B introduces artificial swelling to represent this volume increase.

[0104] For a desired change of volume dV, the artificial swelling is achieved by altering the deformation gradient F determinant of the particle to swell withinducing an isotropic dilation due to the model’s elasticity generating compressive forces to recover a unit determinant (e.g., see Apx K - code below “Apply the volume increase”).

[0105] Based on the velocity gradient (L) of a particle, the amount of artificial swelling added to a particle is configured so that: a. particles with larger linear velocities are affected by more swell, thus the swell can be calculated using a monotonic or linear relationship between swell and linear velocity (e.g., see Apx K, code below “Swell factor based on movement”); b. particles with larger spin tensor fl =(L - Lr), which is equivalent to the angular velocity vector’s magnitude, are affected by more swell, thus the swell can be calculated using a monotonic or linear relationship between swell and the magnitude of the angular velocity for that particle (e.g., see Apx K, code below “Swell factor based on rotations”, although in some implementations the parameter particle_swell.volume_increase_spin_multiplier may always be set to 0.0); and c. particles with larger determinant of their rate of deformation tensor D =(L + Lr) (faster deformations) are affected by more swell, thus the swell can be calculated using a monotonic or linear relationship between swell and the determinant of the rate of deformation for that particle (e g., see Apx K, code below “Swell factor based on deformations”).

[0106] The spin tensor is a skew-symmetric matrix, thus an equivalent 3D angular rotation vector can be extracted by taking three of its off-diagonal components (namely.(3, 2), (1, 3), (2, 1)]), and calculating the magnitude of this 3D vector to provide an estimate of the angular velocity.

[0107] Note that, while adding the artificial swelling, the amount of swell actually applied is clamped / limited (e.g., see the calls to “clamp” in Apx K) to some maximum to avoid overshooting or the swelling causing artifacts.|0108| Thus the artificial swelling includes: a. increasing the particle volume based on its updated velocity such that particles with larger linear velocities are affected by more swell; b. increasing the particle volume based on its updated spin tensor such that particles with larger spin tensor [Q] or equivalent angular velocity vector’ s magnitude are affected by more swell; and / or c. increasing the particle volume based on its determinant of rate of deformation tensor, e g., represented by the selected swell deformation rate multiplier, such that particles with larger determinant of rate of deformation tensor [D] (representing faster deformations) are affected by more swell.

[0109] In some embodiments, (a) and (b) may be unnecessary, so swell can be determined using only (c), the deformation.

[0110] The particle and grid update subprocess 110 may include limiting the volume increase (due to swell) to be above zero, and limiting the maximum particle radius (after swell) to be below a selected threshold The maximum particle radius after swell is the radius of effect of the kernel support function because any value larger than that would result in undesired numerical fractures. For the quadratic B-spline, this implies a maximum particle radius of 1.5 times the cell-width. In some embodiments, the limit is selected to be 1 times the cell-width, e.g., see Apx K, code below “Total volume increase amount”. The input data can define the swell cap, e.g., based on a user selection orprevious measurements. The minimum volume increase is 0.0 to not artificially shrink particles.MPM-MLS10111] The particle and grid update subprocess 110 includes numerically determining the grid using a numerical technique used to simulate continuum materials, e g., using a material point method (MPM), e g., as described in De Vaucorbeil et al. "Material point method after 25 years: Theory, implementation, and applications", in Advances in applied mechanics 53 (2020): 185-398. The numerical technique used in the particle and grid update subprocess 1 10 (e.g., MPM) is mesh-free, which is suitable for simulating large deformations and fractures. This numerical technique may be classified in the “mesh-free” class of numerical methods because the “mesh” (the grid) is only a temporary medium for carrying out the force calculation at each step so it doesn’t suffer from any of the “meshbased” issues, like element inversion, that are caused by the mesh being persistant and deforming.

[0112] In some embodiments, the numerical technique used in the particle and grid update subprocess 110 includes a moving least square (MLS) method of the MPM (referred to as the "MPM-MLS" method), e g., as described in Yuanming Hu et al., "A Moving Least Squares Material Point Method with Displacement Discontinuity and Two-Way Rigid Body Coupling" in ACM Trans. Graph. (SIGGRAPH), Vol. 37, No. 4, Article 150 (August 2018).

[0113] In some embodiments, the numerical technique used in the particle and grid update subprocess 110 includes the MPM-MLS method performed in parallel across the plurality of GPUs on the GPU side 10B (referred to as a "multi -GPU MPM-MLS method"), e g , as described in Xinlei Wang et al., "A Massively Parallel and Scalable Multi-GPU Material Point Method" in ACM Trans. Graph. (SIGGRAPH), Vol 39, No. 4, Article 30 (July 2020). In experimental examples, the multi-GPU MPM-MLS method allowed sub-minute calculation times for blast movement prediction.

[0114] In some embodiments, the MPM-MLS method and the multi-GPU MPM-MLS method include using a quadratic B-spline function as the kernel function.Multi-GPU

[0115] If there is a plurality of GPUs in the GPU side 10B, the particle and grid update subprocess 110 performs the numerical calculations of the particle velocities and the particle 3D locations using the plurality of GPUs in parallel with each other, i.e., such that the plurality of the GPUs are performing the numerical calculations in the particle and grid update subprocess 110 mutually simultaneously. Using the GPUs in parallel can allow for sub-minute calculation times

[0116] This parallel use of the GPUs may be regarded as parallelizing the complete physics simulation pipeline, including the numerical calculation of the velocities / positions, the particle deformations (deformation gradient), the particle stress tensors, and the pressure / intemal energy (for gas particles). The velocity update based on the boundary conditions (the floor) as well as the simulation substep determination are also parallelized. The grid-related bookkeeping, which includes updating its data structure at each simulation step because the grid is sparse, is also parallelized, e.g., see orchestration calls to the CUDA driver to execute the GPU kernels on particles or grid cells in parallel in https: / / github.com / dimforge / sparkl / blob / v0 2. l / src / cuda / cuda_mpm_pipeline.rs#L323- L638 (this is an example GPU loop where each iteration runs one sub-step).Force integration and boundary handling subprocess 112

[0117] As shown in Apx E, which discloses pseudocode implementing the force integration and boundary handling subprocess 112, the force integration and boundary handling subprocess 112 includes the GPU side 10B: a. for each grid cell (e.g., see Apx E: 1): i. performing the explicit integration by setting the cell velocity to be equal to: the momentum and the forces times the time step (dt), divided by the mass (e g., see Apx E:3 and http s : / / github . com / dimforge / sparkl / bl ob / vO .2. 1 / src_kernel s / cuda / gri d_update.rs#L62-L63),ii. handling the one or more solid boundaries represented by the floor 204 or the topography mesh 202 (converted to the heightfields 202B, 204B for use on the GPUs 156) by, on each boundary shape (e.g., see Apx E:5):1. calculating the closest point, i.e., the point on the solid boundary, that is closest to the cell center and normal from the cell center (e g., see Apx E:6, and http s : / / github . com / dimforge / sparkl / blob / vO .2.1 / src_kernel s / c uda / grid_update.rs#L65-L96); and2. projecting the cell velocity based on this point & the normal (e.g., see Apx E:7, and https: / / github.eom / dimforge / sparkl / blob / v0.2.l / src_kernels / c uda / grid update. rs#Ll 10-L149).

[0118] The numerically determining / calculating of the updated velocities (in the velocity field) in the force integration and boundary handling subprocess 112 includes the explicit integration of forces (e g , see Apx E:3). If rocks are given a high stiffness, explicit integration generally implies that the time step must be smaller (than with implicit integration schemes) in order to avoid numerical issues (divergence of force calculation resulting in velocity oscillation or explosions). The explicit integration scheme may be significantly more computationally efficient per time step than an implicit scheme that would involve the expensive resolution of a complex sparse linear system.

[0119] The force integration and boundary handling subprocess 112, and thus the method 100, represents the unbreakable floor 204 under the pre-blast solid material as the heightfield 204B, and using the heightfield 204B as a boundary condition in the repeatedly numerically determining (substep). The boundary condition is used to define the floor 204 of the simulation (otherwise particles would just fall down indefinitely due to gravity). The floor 204 of the simulation can be interpreted as modeling the rocks below the bench that are sufficiently far down in the ground to not be moved by the blast. The force integration and boundary handling subprocess 112 includes repeatedly, at one or more of the timesteps: reducing / eliminating the velocity values of the discrete particles below the unbreakable floor 204; and applying a corrective impulse to particles below the unbreakable floor 204 in order to bring them back above the floor 204. The method includes repeatedly, at one or more of the time steps: reducing / eliminating any downward velocity values of discrete particles adjacent to the unbreakable floor 204: wherein "downward" means "the parts of the velocity that goes in the direction opposite from the floor’s normal at the floor’s point closest to the rock particle." The force integration and boundary handling subprocess 112 includes repeatedly, at one or more of the time steps: reducing / eliminating tangential velocity values of discrete points adjacent to the unbreakable floor 204 depending on the normal velocity and the distance of the discrete point from the unbreakable floor 204, thus modelling the solid ground as a non- deformable, static, and rigid body.

[0120] Thus the solid ground, modeled as a non-deformable, static, rigid body, is represented in the force integration and boundary handling subprocess 112 (and thus effectively coupled with the MPM simulation) (e.g., see friction application in https: / / github.eom / dimforge / sparkl / blob / vO.2.1 / src_kernels / cuda / grid_update.rs#Ll 1 1 - L150) by: a. (optionally) eliminating any grid-cell velocity below the floor topography; b. eliminating any entering velocity from grid cells above the floor topography (e.g., see https: / / github.eom / dimforge / sparkl / blob / v0.2.l / src_kernels / cuda / grid_updat e.rs#L144-L145); and c. eliminating part of the tangential velocity of grid cells above the floor topography, depending on the normal velocity and the distance of the grid cell from the ground (e g., see https: / / github.eom / dimforge / sparkl / blob / v0.2.l / src_kemels / cuda / grid_updat e.rs#L127-L142).More substeps determination subprocess 114

[0121] In the more substeps determination subprocess 114, the CPU side 10A determines whether more substeps are required to cover the current simulation step — the CPU 154 checks if the total time covered by the already executed substeps is equal to the desired simulation step time. If the CPU 154 determines that more substeps are required, the GPU side 10B proceeds to repeat the GPU simulation loop 126 by commencing the update sparse grid subprocess 106 again. If the GPU side 10B determines that no more substeps are required, the GPU side 10B ends the GPU simulation loop 126 and instructs / signals the CPU side 10A to commence the convergence or divergence subprocess 116.Convergence / divergence check subprocess 116 and decision subprocess 118

[0122] As shown in Apx F, which discloses pseudocode implementing the convergence or divergence subprocess 1 16 and the convergence / divergence decision subprocess 1 18, the convergence or divergence subprocess 116 includes the CPU side 10A: a. the CPU side 10A reading the new particle positions and velocities from the GPU side 10B, as received from the GPU simulation loop 126 (e.g., see Apx F : 1 ) in a memory transfer from the GPU(s) 156 to the CPU 154, as shown in FIG. 1C; and b. for each particle, the CPU side 10A (e.g., see Apx F:2): i. determining if any of its coordinates or velocity is not a number (NaN) or is too large — e.g., determined by the particles leaving a user-defined perimeter (e.g., set to the perimeter of the floor 204), or reaching a point higher than a user-defined altitude (e.g., 100 meters above the highest particles at their initial position) from the original blast (e g., see Apx F:3); and, if so, exiting the main simulation loop 120, and considering the simulation to be "diverged" (e.g., see Apx F:4), otherwise ii. determining if the magnitude of its velocity is larger than a small user-defined value ("epsilon"), e g., 0.5m.s1(see Apx F:5); and, if so, marking the simulation as not converged (e.g., see Apx F:6).

[0123] As shown in Apx F, the convergence / di vergence decision subprocess 118 includes the CPU side 10A performing the following: a. if the simulation is marked (by the convergence or divergence subprocess 116) as not converged, or if some explosives have not been processed yet (determined if their ignition-time+lifetimes are greater than the simulated time) (e.g., see Apx F:7), continuing to the next step in the main simulation loop 120, i.e., repeating the update explosives subprocess (e.g., see Apx F:8), or b. if the simulation is not marked (by the convergence or divergence subprocess 116) as not converged, and if no explosives remain to be processed (determined if their ignition-time+lifetimes are all less than the simulated time) (e.g., see Apx F:9), exiting the main simulation loop 120 and considering the simulation to be converged (e.g., see Apx F: 10).Reactive movement correction subprocess 122

[0124] In the event where the user provides a post-blast topography, the predictive simulation can be followed by a post-processing stage that fits the fragmented particles into the post-blast volume enclosed by the floor 204 and post-blast topography. This involves two artificial forces: (1) a spring-like attraction force moving particles towards the post-blast solid interior if they are outside of it, and, (2) a suction force applied to particles near empty space, to even the particle distribution within the post-blast volume.

[0125] As shown in FIG. 10, the reactive movement correction subprocess 122, and thus the method 100, can include fitting the discrete particles to a measured post-blast topography 1002 (typically extracted from drone scans), including by: applying an artificial attraction force 1004 (using a spring stiffness coefficient) to ones of the discrete particles in an upper region / layer 1006 of the processed discrete particles to move the discrete particles to be within the measured post-blast topography 1002. The reactive movement correction subprocess 122, and thus the method 100, can include fitting ones of the discrete particles in the upper region / layer 1006 to fill a volumetric gap 1008 under the measured post-blast topography 1002 and above the upper region / layer 1006 by: applyingan artificial suction force 1010 (using a suction stiffness coefficient) to the discrete particles to move the discrete particles in the upper region / layer 1006 under the measured post-blast topography 1002 to fill the post-blast volumetric gap 1008. The (spring-like) attraction force is calculated as mass-spring-damper system.

[0126] Thus, given the floor topography 204 and a post-blast topography (typically extracted from drone scans), the method 100 simulates blasting of the rock and fits them into the post-blast volume. The method 100 may be regarded as preserving the relationship between pre-blast material and post-blast material to identify the nature / composition of the post-blast material in its locations based on the known nature / composition and locations of the pre-blast material.

[0127] The spring-like attraction is calculated as mass-spring-damper system with critical damping. With inputs including a spring stiffness parameter k, particle mass m, distance d, and normal n, — of which only stiffness parameter k, is an input value, where the other parameters (mass, distance, normal) are all dependent on the current particle position (and the mass is specific to each particle) — the attraction force is given by f = - kdn — 2\lkmv.

[0128] Note that k = 0 for particles located inside of the post-blast volume.|0129] Given a suction stiffness parameter k’, the suction force s applied to a particle with position p is calculated as a weighted sum of the mass-difference between the cell, with mass mt, containing the particle, and all the grid cells with positions cyincluded in the non-zero domain of a compact weighting function IR3-> IR:

[0130] In some implementations, is selected as the same kernel function as the rest of the simulation loop (quadradic B-spline) to ensure compatibility with the existing optimized grid-to-particle-to-grid transfer kernel in the particle and grid update subprocess 110, resulting in no performance penalty. There needs be only 1 kernel since the method 100 picks the same function, e g., quadratic B-spline.

[0131] As shown in Apx G, which discloses pseudocode implementing the reactive movement correction subprocess 122, the reactive movement correction subprocess 122 includes the CPU side 10A: a. computing the post-blast solid where the particles have to be fit (e.g., see Apx G: l); b. performing (also referred to as "running") the GPU simulation loop 126 but (e.g., see Apx G:2) i. without any explosive particles in the 3D grid (e.g., see Apx G:3), and ii. forcing any particle that didn’t move to remain static — these are the particles that moved (between their initial position before simulation, and final position after simulation) by a distance smaller than a selected threshold, e.g., 0.05m, or particles that moved by more than 0.05m, but less than another selected threshold (here 3.0m) and were not fragmented by the simulation, (e.g., see Apx G:4), and iii. applying the artificial spring and suction forces in the particle and grid update subprocess 1 10 (e g., see Apx G:5); and c. instead of performing the convergence / divergence check subprocess 116 and decision subprocess 118, performing the main simulation loop 120 for a selected number of times (or "iterations"), optionally with a plurality of mutually different stiffness coefficients for the artificial forces, e.g., running it 50 times with low stiffness coefficients for the artificial forces, and150 times with higher stiffness coefficients (e.g., see Apx G:6); and d. projecting any remaining non-static outliers into the boundary of the postblast solid — this effectively computes the point on the post-blast solid surface that is closest to an outlier particle, and teleports that outlier particle to that closest point directly (e.g., see Apx G:7).

[0132] The reactive movement correction subprocess 122 (referred to as a "postprocessing step") brings the simulation closer to the ground-truth by correcting movement and swell errors, allowing more accurate localization of post-blast materials within the muck pile.Output subprocess 124|0133| The output subprocess 124 includes: determining post-blast locations of the solid material using the repeatedly updated 3D positions of the discrete particles and the preblast solid properties of those discrete particles, which gives movement and fragmentation information.System 150

[0134] As shown in FIG. IB and FIG. 1 C, the system 150 includes: a. a CPU 154, including RAM, in the CPU side 10A; b. a machine readable storage, including a storage disk and / or cloud storage, accessible by the CPU 154; and c. one or more GPUs 156A,156B, . . . in communication with the CPU 154.

[0135] FIG. IB and FIG. 1C show memory transfers between the CPU 154 and the GPUs 156A,156B: e.g., "204B" and "202B" represent the GPU (heightfield) version of the meshes 204 and 202, and "Apx B:2" means Apx B, line 2.

[0136] In an example, the system 150 includes: a. as the CPU 154, an AMD Ryzen Threadripper 3970X; b. as each GPU 156, an NVIDIA RTX A6000; and c. as a motherboard, an MSI Creator TRX40.

[0137] The CPU 154 and the CPUs 156 are in the form of physical processors, and the storage 152 includes tangible, non -transitory, machine-readable memory that stores thecomputer-readable instructions (configured to cause the or each CPU 154 and the GPUs 156 to perform the method 100) described hereinbefore. These instructions, when executed by the physical processors, effectuate the method 100 and its subprocesses described herein.

[0138] Mining methodsThe hereinbefore-described method 100 is used in mining methods to improve mining outcomes. The mining methods can include physically measuring the pre-blast solid material to determine the one or more mechanical properties of the solid material before the explosive blasting. The mining methods can also include physically excavating and / or routing the blasted solid material (i.e., after the explosive blasting) based on the determined post-blast locations of the solid material. The mining methods can include optimising dig direction (including angle) during the excavating by controlling an excavator (or "digger") using the determined post-blast locations: this can improve efficiency because changing dig directions can be cumbersome. In structurally- controlled operations, changing dig direction can prevent serious dilution and increase grades. The routing of the excavated blasted solid material can be selected / controlled by the determined post-blast locations in order to select or minimize ore dilution, and select an optimal destination for each excavated load based on the determined ore content and the operation of the mine, e.g., routing loads to waste, stockpiles or the mill based on a determined ore grade from the determined post-blast locations (from the method 100). The determined 3D post-blast locations of the solid material can be used in mining plant equipment, such as the diggers, which can receive the post-blast locations to update their three-dimensional geological maps of the muck pile with the new ore boundaries during excavation.Interpretation

[0139] As used herein, the term “set” corresponds to or is defined as a non-empty finite organization of elements that mathematically exhibits a cardinality of at least 1 (i.e., a set as defined herein can correspond to a unit, singlet, or single element set, or a multiple element set), in accordance with known mathematical definitions (for instance, in a manner corresponding to that described in An Introduction to Mathematical Reasoning: Numbers, Sets, and Functions, "Chapter 11 : Properties of Finite Sets" (e g., as indicated on p. 140),by Peter J. Eccles, Cambridge University Press (1998)). Thus, a set includes at least one element. Tn general, an element of a set can include or be one or more portions of a system, an apparatus, a device, a structure, an object, a process, a procedure, physical parameter, or a value depending upon the type of set under consideration.

[0140] The FIGs. included herewith show aspects of non-limiting representative embodiments in accordance with the present disclosure, and particular structural elements shown in the FIGs. may not be shown to scale or precisely to scale relative to each other. The depiction of a given element or consideration or use of a particular element number in a particular FIG or a reference thereto in corresponding descriptive material can encompass the same, an equivalent, an analogous, categorically analogous, or similar element or element number identified in another FIG. or descriptive material associated therewith.

[0141] The presence ofin a FIG. or text herein is understood to mean "and / or" unless otherwise indicated, i.e., “A / B” is understood to mean “A” or “B” or “A and B”.

[0142] The recitation of a particular numerical value or value range herein is understood to include or be a recitation of an approximate numerical value or value range, for instance, within + / - 20%, + / - 15%, + / - 10%, + / - 5%, + / - 2.5%, + / - 2%, + / - 1%, + / - 0.5%, or + / - 0%. The term "substantially" can indicate a percentage greater than or equal to 50%, 60%, 70%, 80%, or 90%, for instance, 92.5%, 95%, 97.5%, 99%, or 100%.

[0143] Many modifications will be apparent to those skilled in the art without departing from the scope of the present invention. Reference to one or more embodiments herein, e.g., as various embodiments, many embodiments, several embodiments, multiple embodiments, some embodiments, certain embodiments, particular embodiments, specific embodiments, or a number of embodiments, need not or does not mean or imply all embodiments.

[0144] Throughout this specification and the claims or statements which follow, unless the context requires otherwise, the word "comprise", and variations such as "comprises" and "comprising", will be understood to imply the inclusion of a stated integer or step or groupof integers or steps but not the exclusion of any other integer or step or group of integers or steps.

[0145] Any reference in this specification to any prior publication (or information derived from it) or to any known matter (e g., the experiments and examples mentioned hereinbefore), is not, and should not be taken as, an acknowledgment or admission or any form of suggestion that the prior publication (or information derived from it) or known matter forms part of the common general knowledge in the field of endeavour to which this specification relates.CODE APPENDTCTESCODE APPENDIX AUpdate explosives

[0001] For each explosive particle live in the simulation:

[0002] If it existed for a duration greater than its lifetime:

[0003] Remove the particle from the simulation.

[0004] For each remaining explosive:

[0005] If the simulated time is equal to the hole timing:

[0006] Discretize the explosive into JWL particles.

[0007] Send these particles to the GPU.

[0008] Set the next simulation step length so it doesn’t exceed the lowest timing of the explosives that were not inserted yet.CODE APPENDIX BUpdate sparse grid

[0001] Identify chunks of 4x4x4 cells containing at least one particle.

[0002] Resize the vector of chunks based on the required number of chunks.

[0003] Sort particle indices based on their containing chunk indices.

[0004] Allocate and populate a hash-map as an inventory of existing chunks

[0005] / * Multi-GPU only: * /

[0006] Identify halo-blocks: blocks in the radiuses of effect particles originating from two GPUs.CODE APPENDIX CEstimate substep length10001] Set substep := simulation step length

[0002] For each particle:

[0003] Set substep := min(substep, cell-width / particle-velocity)

[0004] Set substep := min(substep, value based on the particle’s constitutive model to satisfy the CFL stability conditions)|0005| Read minimum substep length from GPUCODE APPENDIX DParticle and grid update (G2P2G)

[0001] / * [G2P]Using the MLS approximation: * /

[0002] Transfer velocities from grid cells to particles.

[0003] Calculate particle velocity gradients from grid velocities.

[0004] Calculate averaged positive elastic strain energy from cells to particles

[0005] Calculate the artificial suction force based on mass difference between the cell containing the particle and adjacent cells(only when running reactive movement correction).|0006| / * [Particle update] * /

[0007] For each particle:

[0008] Update fragmentation status based on the particle elastic strain energy (eigenerosion).100091 Integrate the velocity gradient to update the deformation gradient.

[0010] Calculate particle stress using the constitutive model and deformation.|00111 Update the particle positive elastic strain energy.

[0012] Calculate the spring-like attraction force to the particle if it is outside the post-blast solid (only when running reactive movement correction)|0013| For each fragmented particle:

[0014] Project the particle stress based on plasticity model.

[0015] Apply artificial swell to the particle.

[0016] / * [P2G] Using the MLS approximation: * /

[0017] Transfer momentum, mass, and forces from particles to grid cells.

[0018] Transfer positive elastic strain energy from particles to grid cells.CODE APPENDIX EForce integration & boundary handling

[0001] For each grid cell:

[0002] / * Explicit integration. * /

[0003] Set cell-velocity := (momentum + forces x dt) / mass

[0004] / * Solid boundary handling. * /

[0005] For each boundary shape:

[0006] Calculate the closest point and normal from the cell center.

[0007] Project the cell velocity based on this point & normal.CODE APPENDIX F(Periodically) Check for convergence of divergence

[0001] Read the new particle positions and velocities from the GPU

[0002] For each particle:

[0003] If any of its coordinates or velocity is NaN or is too far from the original blast:

[0004] Exit the simulation and consider it diverged.

[0005] If the magnitude of its velocity is larger than an epsilon:

[0006] Mark the simulation as not-converged.

[0007] If the simulation is marked as not-converged or some explosives have not been processed yet:

[0008] Continue to the next step.

[0009] Else:

[0010] Exit the simulation and consider it converged.CODE APPENDIX G(Optional) Reactive movement correction

[0001] Compute the post-blast solid where the particles have to be fit.|0002| Run a simulation loop similar to the above; however:

[0003] - Without any explosive involved.

[0004] - Forcing any particle that didn’t move to remain static.

[0005] - With the artificial spring and suction forces included in the Particle and grid update (G2P2G) step.

[0006] - Instead of a convergence check, the simulation loop can be run only a fixed number of times, e.g., run it 50 times with low stiffness coefficients for the artificial forces, and 150 times with higher stiffness coefficients.

[0007] Project any remaining non-static outliers into the boundary of the post-blast solid. NOTE: this step 0007 happens on the CPU only.CODE APPENDIX FT - JWL CONSTITUTIVE MODEL use crate: math:: {Matrix, Real, Vector}; / * 3D matrix, f32, 3D vector * / use sparkl3d_core: dynamics: models: :ActiveTimestepBounds;#[cfg(not(feature = "std"))] use na::ComplexField; pub struct JwlEos { pub detonation speed: Real, pub a: Real, pub b: Real, pub rl : Real, pub r2: Real, pub omega: Real, pub directional coeff: Vector<Real>,} / / can be obtained by solving a linear equation: p = coeff_a + coeff_b * E' let coeff a = ignition factor * (terml + term2), let coeff b = ignition factor * self.omega / rel v; let e = particle_internal_energy, / / NOTE: this formula is based on the corresponding book, except we multiply the numerator / / and denominator by 'parti cl e_volumeO' to avoid some divisions by very small values. let pressure =(coeff_a * particle_volumeO + coeff_b * e) / (particle_volumeO + coeff_b * dv); let energy = e - dv * pressure;(pressure, energy)} pub fn is fluid(&self) -> bool { true pub fn kirchhoff_stress(&self, particle_pressure: Real, particle_fluid_deformation_gradient_det: Real,) -> Matrix<Real> {Matrix: :from_diagonal(&(self. directi onal coeff* (-particle_pressure * particle fluid deformation gradient det)), )} pub fn update_intemal_energy_and_pressure(&sel , particle_velocity_gradient: &Matrix<Real>, particle volumeO: Real, particle volume fluid: Real, particle ignition time: &mut Real, particle internal energy: &mut Real, particle_pressure: &mut Real, dt: Real, cell width: Real,) { / / Note: we use ' -parti cle.ignition time', which is equal to(init_particle_ignition_time - elapsed_time)' / / which equals ' elapsed time - init_particle_ignition_time' . let ignition_factor = if *particle_ignition_time > 0.0 {0.0let (pressure, energy) = self.solve_pressure_and_energy( particle velocity gradient, particle volumeO, particle volume fluid,*particle_intemal_energy, dt, cell_width, ignition factor,);*particle_internal_energy = energy;*particle_pressure = pressure;*particle_ignition_time -= dt;} pub fn active timestep bounds(&self) -> ActiveTimestepBounds {ActiveTimestepBounds:: CONSTITUTIVE MODEL BOUND| ActiveTimestepBounds: :PARTICLE_VELOCITY_BOUND| ActiveTimestepBounds: :PARTICLE_DISPLACEMENT_BOUND| ActiveTimestepBounds: : SINGLE PARTICLE STABILITY BOUND pub fn timestep_bound(&self, parti cl e densityO: Real, particle_density_fluid: Real, particle_fluid_deformation_gradient_det: Real, particle internal energy: Real, particle_pressure: Real, cell width: Real,) -> Real { cell width / self.adiabatic sound speed( particle densityO, particle density fluid, particle fluid defor ation gradient det. particle internal energy, particle pressure,)}CODE APPENDIX I - CONSTITUTIVE MODELS ENUMERATION#[cfg(not(feature = "std"))] use na::ComplexField; use sparkl_blast_core::{JwlEos, Parti cl eExplosive}; use sparkl_core::dynamics::{ParticleData, Parti cl eVolume}; use sparkl_core::math:: {Matrix, Real}; use sparkl core: :prelude: : { ActiveTimestepBounds, CorotatedLinearElasticity,ParticleVelocity}; use sparkl_kemels: :DevicePointer;#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]#[derive(cust_core::DeviceCopy, Copy, Clone, PartialEq)]#[repr(C)] pub enum GpuBlastConstitutiveModel {CorotatedLinearElasticity(CorotatedLinearElasticity, DevicePointer<ParticleData>), JwlEos(JwlEos, DevicePointer<ParticleExplosive>), / * Appendix H * / } impl GpuBlastConstitutiveModel { pub fn is fluid(&self) -> bool { match self {Self: :CorotatedLinearElasti city (m, _) => m.is_fluid(),Self::JwlEos(m, _) => m.is_fluid(),pub unsafe fn pos_energy(&self, particle_id: u32, particle volume: &ParticleVolume,_particle_phase: Real,) -> Real { match self {Self::CorotatedLinearElasticity(m, d) => { let hardening = (*d.as_ptr().add(particle_id as usize)).elastic_hardening; m.pos energy (particle volume, deformation gradient, hardening)}Self::JwlEos(..) => 0.0,}} pub unsafe fn kirchhoff_stress(&self, particle id: u32, particle volume: &ParticleVolume, particle_phase: Real, velocity gradient: &Matrix<Real>,) -> Matrix<Real> { match self {Self: :CorotatedLinearElasti city (m, d) => { let hardening = (*d.as_ptr().add(particle_id as usize)).elastic_hardening; m.kirchhoff_stress( particle_phase, hardening,&particle_volume.deformation_gradient,)}Self::JwlEos(m, d) => { let pressure = (*d.as_ptr() add(particle_id as usize)).pressure; m.kirchhoff_stress(pressure, particle_volume.fluid_deformation_gradient_det()) }}} pub unsafe fn update internal energy and_pressure(&self, particle id: u32, particle volume: &ParticleVolume, dt: Real, cell width: Real, velocity gradient: &Matrix<Real>,) { match self {Self::CorotatedLinearElasticity( .) => {}Self::JwlEos(m, d) => { let ignition time = &mut (*d.as_mut_ptr().add(particle_id as usize)).ignition_time; let internal_energy =&mut (*d.as_mut_ptr() add(particle_id as usize)).internal_energy; let pressure = &mut (*d.as_mut_ptr() add(particle_id as usize)). pressure; m.update_internal_energy_and_pressure( velocity gradient, particle volume.volumeO, particle volume.volume fluid(), ignition time, internal energy, pressure, dt, cell_width,);}pub fn active_timestep_bounds(&self) -> ActiveTimestepBounds { match self {Self::CorotatedLinearElasticity(m, _) => m.active_timestep_bounds(),Self::JwlEos(m, _) => m.active timestep boundsQ,pub unsafe fn timestep_bound(&self, particle_id: u32, parti cl e volum e : &Parti cl eV olum e, particle_vel: &Parti cl eVelocity, cell width: Real,) -> Real { match self {Self::CorotatedLinearElasticity(m, d) => { let hardening = (*d.as_ptr().add(particle id as usize)). elastic hardening; m . timestep b ound( particle volume. densityOQ,&particle_vel. vector, hardening, cell_width,)}Self::JwlEos(m, d) => { let internal_energy = (*d.as_ptr(').add(particle_id as usize)).internal_energy; let pressure = (*d.as_ptr().add(particle_id as usize)).pressure; m.timestep_bound( particle_volume.densityO(), particle_volume. density _fluid(), parti cl e vol um e . fl ui d deform ati on gradi ent_det(), intemal_energy, pressure, cell_width,)}}}}CODE APPENDIX J - PER-PARTICLE EXPLOSIVE AND SWELL ATTRIBUTES. use crate:: math:: Real;#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] #[derive(Copy, Clone, Debug, PartialEq, bytemuck: od, bytemuck: :Zeroable)] #[repr(C)] pub struct ParticleSwell { / / For particles with artificial volume modification. pub volume increase: Real, pub volume increase velocity: Real, pub volume increase minimal velocity: Real, pub volume increase bounds: [Real; 2], pub volume increase deformation multiplier: Real, pub volume increase movement multiplier: Real, pub volume increase spin multiplier: Real,} impl Default for ParticleSwell { fn default() -> Self {Self { volume increase: 0.0, volume increase velocity: 0.0, volume increase minimal velocity: 0.0, volume increase bounds: [0.0, 0.0], volume increase deformation multiplier: 1.0, volume_increase_movement_multiplier: 0.0, volume_increase_spin_multiplier: 0.0,}}}#[cfg_attr(feature = "cuda", derive(cust_core::DeviceCopy))]#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))] #[derive(Copy, Clone, Debug, PartialEq, bytemuck: od, bytemuck: :Zeroable)] #[repr(C)] pub struct ParticleExplosive { pub ignition time: Real, pub internal energy: Real, pub pressure: Real,} impl Default for ParticleExplosive { fn default() -> Self {Self { ignition_time: 0.0, intemal_energy: 0.0, pressure: 0.0,CODE APPENDIX K - PARTICLE UPDATE use crate: :GpuBlastConstitutiveModel; / * Appendix I * / use sparkl_blast_core:: {Parti cl eExplosive, ParticleSwell}; / * Appendix J * / use sparkl_core: : dynamics: :Parti cl eFracture; use sparkl_core::math::{Matrix, Real, Vector}; use sparkl_core::prelude::{ActiveTimestepBounds, ParticlePhase, ParticlePosition, ParticleStatus, ParticleVelocity, ParticleVolume,}; use sparkl_kemels::cuda: :{InterpolatedParticleData, Parti cleUpdater}; use sparkl kernels:: {DevicePointer, GpuColliderSet, GpuGridProj ectionStatus,GpuParti cl eModel } ;#[cfg(not(feature = "std"))] use na::ComplexField;#[derive(Copy, Clone, cust_core::DeviceCopy)]#[repr(C)] pub struct BlastParticleUpdater { pub models: DevicePointer<GpuParticleModel>, pub blast models: DevicePointer<GpuBlastConstitutiveModel>, pub parti cl e swell: DevicePointer<ParticleSwell>, / / For artificial swell. pub particle explosive: DevicePointer<ParticleExplosive>, / / For internal energy and pressure. pub particle fracture: DevicePointer<ParticleFracture>, / / For eigenerosion. pub artificial_pressure_stiffness: Real,} impl ParticleUpdater for BlastParticleUpdater { fn artificial_pressure stiffness(&self) -> Real { self, artifi ci l_pres sure stiffnes s}#[inline(always)] unsafe fn estimate_particle_timestep_length(&self, cell width: Real, particle id: u32, particle status: &ParticleStatus, parti cle volume : &Particl eV olume, particle_vel: &Parti cl eVelocity,) -> (ActiveTimestepBounds, Real) { let blast model = &*self.blast_models.as_ptr().add(particle_status.model_index); let active timestep bounds = blast_model.active_timestep_bounds(); if active_timestep_bounds.contains(ActiveTimestepBounds::CONSTITUTIVE_MODEL_B OUND) {( active timestep bounds, blast model timestep_bound(particle_id, particle_volume, particle_vel, cell width),)} else {(active timestep bounds, Real : :MAX)}}#[inline(always)] unsafe fn update particle and compute kirchhoff stress(&self, dt: Real, cell width: Real, colliders: &GpuColliderSet, particle id: u32, particle_status: &mut ParticleStatus, particle_pos: &mut ParticlePosition, particle_vel: &mut Parti cl eVelocity, particle volume: &mut ParticleVolume, particle_phase: &mut ParticlePhase, interpolated data: &mut InterpolatedParticleData,) -> Option<(Matrix<Real>, Vector<Real>)> { let parti cl e swell = &mut *self.particle_swell as_mut_ptr().add(particle_id as usize); let parti cl e explosive = &mut *self.particle explosive,as_mut_ptr(),add(particle id as usize); let particle fracture = &mut *self.particle fracture.as mut_ptr().add(particle id as usize); let model = &*self.models.as_ptr().add(particle_status.model_index); / / NOTE: we always use our own constitutive models (Appendix I). let constitutive model =&*self.blast_models.as_ptr().add(particle_status.model_index); / *particle_volume.deformation_gradient[(0, 0)] +=(interpolated data. velocity _gradient_det * dt)particle volume, dt, cell_width,&interpolated_data velocity gradient,);} / ** Apply artificial swelling.* / if particle_phase.phase == 0.0 { / / Swell factor based on deformations. let swell def velocity factor = (-(0.5* (&interpolated data. velocity gradient+ interpolated data.velocity gradient.transpose())),determinant()* particle swell.volume increase deformation multiplier).clamp(0.0, 10.0).abs(), / / Swell factor based on rotations. let spin_tensor = 0.5* (&interpolated_data. velocity gradient- interpolated_data.velocity_gradient.transpose()); let swell_spin_velocity_factor = Vector: :new(spin_tensor[(0, 1)], spin_tensor[(0, 2)], spin_tentsor[(l , 2)]).normO* particle swell.volume increase spin multiplier) .clamp(0.0, 1.0); / / Swell factor based on movement. let swell mov velocity factor = (particle vel.vector.xy().norm()* particle swell.volume increase movement multiplier).clamp(0.0, 1.0); let swell velocity factor = swell_def_velocity_factor + swell_spin_velocity_factor + swell mov velocity factor; / / Total volume increase amount. let new_volume_increase = (parti cle_swell.volume_increase+ (particle_swell.volume_increase_velocity * swell_velocity_factor).max(particle_swell.volume_increase_minimal_velocity)* dt).clamp( particle_swell.volume_increase_bounds[0], particle_swell.volume_increase_bounds[l],); / / Apply the volume increase. let delta volume = new volume increase - parti cle_swell.volume_increase; let j_factor = parti cle_volume.volumeO / (particle volume volumeO + delta_volume); particle_volume.volumeO += delta_volume; particle_volume.deformation_gradient *= j_factor.powf(1.0 / 3.0); particle_swell.volume_increase = new_volume_increase;} / ** Apply plasticity.* / if let Some(plasticity) = &model. plastic model { plasticity. update_particle(particle_id, particle_volume, particle_phase.phase); } let mut penalty force = Vector: :zeros(), / ** Artificial spring-like force.* / if particle_status.is_static { parti cle_vel . vector, fill (0.0); interpolated_data,velocity_gradient.fill(0,0);} else { if let GpuGridProjectionStatus::Tnside(collider_id) = interpolated_data.projection_status{ if let Some(collider) = colliders.get(collider id) { if collider.penalty stiffness > 0.0 { let depth = interpolated data.proj ection scaled dir.normQ; if depth > 0.0 { let effective depth = depth.min(lO.O); / / Spring with critical damping. penalty force += (-interpolated data.projection scaled dir)* (effective_depth / depth)* collider.penalty _stiffness- particle_vel.vector* 2.0 / ** Don’t crash the whole simulation if a particle diverged. an() 0.0 / / Isolated particles tend to accumulate a huge amount of numerical / / error, leading to completely broken deformation gradients. / / Don’t let them destroy the whole simulation by setting them as failed.|| (!is fluid && particle volume. deformation gradicnt[(O, 0)].abs() > 1.0e4){ particle status.failed = true; particle explosive.intemal energy = 0.0; particle_explosive.pressure = 0.0; particle_volume.deformation_gradient = Matrix: :identity(); return None; / ** Update Pos energy.* / { let energy = constitutive_model.pos_energy(particle_id, particle_volume, particle_phase. phase); particle_phase.psi_pos = particle_phase.psi_pos.max(energy);} / ** MPM-MLS: the AP1C affine matrix and the velocity gradient are the same. * / if I particle status.failed { let stress = constitutive model. kirchhoff stress( particle id, particle volume, particle_phase.phase,&interpolated_data. velocity gradient,);Some(( stress, penalty force))} else {Some((Matrix: :zeros(), Vector: :zeros()))CODE APPENDIX L - PARTICLE DEFINITION use crate: : {Parti cl eExplosive, Parti cl eSwell}; / * Appendix J * / use spark!3d_core:: dynamics::}ParticleData, ParticlePhase, ParticlePosition, ParticleStatus, ParticleVelocity,ParticleVolume, }; use sparkl3d_core::math:: {Point, Real}; use sparkl3d core: :prelude: :ParticleFracture;#[cfg(not(feature = "std"))] use na::ComplexField;#[cfg_attr(feature = "serde-serialize", derive(Serialize, Deserialize))]#[derive(Copy, Clone, Default, Debug)] pub struct BlastParticle { pub position: ParticlePosition, pub velocity: ParticleVelocity, pub volume: ParticleVolume, pub status: ParticleStatus, pub phase: ParticlePhase, pub data: ParticleData, pub swell: ParticleSwell, pub explosive: Parti cl eExplosive, pub fracture: ParticleFracture,} impl BlastParticle { pub fn new(model: usize, pos: Point<Real>, radius: Real, densityO: Real) -> Self { Self: :with internal energy(model, pos, radius, densityO, 0.0)} pub fn with internal energy} model: usize, pos: Point<Real>, radius: Real, densityO: Real, intemal_energy_per_volume_unit: Real,) -> Self {let volumeO = (radius * 2.0).powi(3);Self { status: ParticleStatus { model index: model,.. Default : : default()}, volume: ParticleVolume { volumeO, radiusO: radius, mass: volumeO * densityO,Default:: default})}, position: Parti cl ePosition { point: pos }, explosive: ParticleExplosive { intemal energy: intemal_energy_per_volume_unit * volumeO,..Default::default()},.. Default: : default()}CODE APPENDIX M - EXAMPLE OF JWL EXPLOSIVE VALUES use crate: :math:: {Point, Real, Vector}; use crate: :{BlastParti cl e, JwlEos}; / * Appendix L, Appendix H * / / / See https: / / www.researchgate.net / publication / 281440526_Determination_of_the_JWL_Consta nts_for_ANFO_and_Emulsion_Explosives_from_Cylinder_Test_Data / / "Determination of the JWL Constants for ANFO and Emulsion Explosives from Cylinder Test data / / See table 3 pub struct AU50; impl AU50 { pub const FACTOR: Real = 1.0; pub const D: Real = 3233.0; pub const DENSITYO: Real = 830.0 * Self::FACTOR; pub const JWL EO: Real = 1.562e9 * Self::FACTOR; pub const JWL A: Real = 216.044e9 * Self::FACTOR; pub const JWL B: Real = 1.838e9 * Self::FACTOR, pub const JWL Rl : Real = 7.162;pub const JWL R2: Real = 0.865; pub const JWL OMEGA: Real = 0.34; pub const BULK COEFF 1 : Real = 1.5; pub const BULK COEFF2: Real = 0.06; pub fn constitutive_model() -> JwlEos {JwlEos { detonation_speed: Self::D, bulk coeff 1 : Self: :BULK_COEFF 1 , bulk_coeff2: Self::BULK_COEFF2, a: Self::JWL_A, b: Self::JWL_B, rl: Self::JWL_Rl, r2: Self::JWL_R2, omega: Self: : J WL OMEGA, directi onal coeff: Vector: :repeat(1.0),}} / / NOTE: The model integer is the index of the instance of ' Self: constitutive model})' / / added to the ' CudaParti cl eModel Sef . pub fn particle_with_eO( model: usize, position: Point<Real>, radius: Real, eO: Real,) -> BlastParticle {BlastParticle: :with_internal_energy(model, position, radius, Self::DENSITY0, eO)} pub fn parti cl e(model: usize, position: Point<Real>, radius: Real) -> BlastParticle { BlastParticle: :with_intemal_energy(model, position, radius, Self::DENSTTY0,Self::JWL_E0)}}CODE APPENDIX N - GENERIC JWL EXPLOSIVE DEFINITION. use crate: :math:: {Point, Real, Vector}; use crate:: {BlastParticle, JwlEos}; / / See https: / / www.researchgate.net / publication / 281440526_Determination_of_the_JWL_Consta nts_for_ANFO_and_Emulsion_Explosives_from_Cylinder_Test_Data / / "Determination of the JWL Constants for ANFO and Emulsion Explosives fromCylinder Test data / / See table 3 pub struct JWLExplosive { pub detonation speed: Real, pub density: Real, pub jwl_eO: Real, pub jwl a: Real, pub jwl b: Real, pub jwl rl : Real, pub jwl_r2: Real, pub jwl omega: Real, pub bulk coeffl : Real, pub bulk coeff2: Real,} impl JWLExplosive { pub fn constitutive_model(&self) -> JwlEos {JwlEos { detonation_speed: self.detonation_speed, bulk coeffl : self.bulk coeffl, bulk_coeff2: self.bulk_coeff2, a: self.jwl a, b: self.jwl_b, rl: self.jwl_rl, r2: self.jwl_r2, omega: self.jwl_omega, directi onal coeff: Vector: :repeat(l .0),}} pub fn parti cl e(& self, model: usize, position: Point<Real>, radius: Real) -> BlastParticle{BlastParticle: :with internal energy(model, position, radius, self. density, self.jwl eO)}}

Claims

CLAIMS1. A method, performed with one or more physical processors, for predicting solid material dislodgement due to explosive blasting, the method including: representing, with the one or more physical processors, pre-blast solid material with a set of discrete particles, each with a 3D location and one or more mechanical properties of the solid material at that 3D location before the explosive blasting; representing, with the one or more physical processors, at least one explosive blast with a blast time and a set of explosive particles, each with a 3D blast location in the pre-blast solid material and an initial internal energy, repeatedly, with the one or more physical processors and for a plurality of time steps, defining a 3D grid surrounding the discrete particles (e g., plus two empty grid cells in each direction), and numerically determining / calculating updated velocities and updated 3D locations of the discrete particles based on the mechanical properties and the 3D locations, using the 3D grid; for at least one of the time steps corresponding to the or each blast time, including the explosive particles with the discrete particles to numerically determine the velocities and updated 3D locations of the discrete particles, and numerically determining / calculating the updated velocities and the updated 3D locations at the time steps after the or each blast time using the initial internal energies; and determining, with the one or more physical processors, post-blast locations of the solid material using the repeatedly updated 3D positions of the discrete particles and the pre-blast solid properties of those discrete particles.

2. The method of claim 1, wherein the repeated numerically determining of the velocities and the updated 3D locations includes performing numerical calculations of the velocities and the updated 3D locations using a plurality of graphics processing units (GPU) in parallel.

3. The method of claim 2, wherein the repeatedly numerically determining of the velocities and the updated 3D locations includes using a quadratic B-spline function to combine the mechanical properties of ones of the discrete particles adjacent to each discrete particle.

4. The method of claim 3, wherein the repeatedly numerically determining includes a Galerkin-style moving least squares (MLS) discretization.

5. The method of any one or more of the preceding claims, wherein the mechanical properties include a Poisson coefficient value and a Young modulus value, and the numerically determining / calculating includes determining / calculating the velocities and updated 3D locations using the Poisson coefficient values and the Young modulus values.

6. The method of any one or more of the preceding claims, wherein the mechanical properties include a tensile strength value, and the method includes repeatedly, for some or all of the plurality of time steps: determining a peak effective energy release rate for each discrete particle to determine a fragmentation status for that discrete particle, wherein the fragmentation status is determined to be fragmented if the peak effective energy release rate for that discrete particle is above a predetermined threshold value; and numerically determining / calculating the velocities and the updated 3D locations of the discrete particles based on their respective fragmentation statuses, including reducing / zeroing the tensile strength value in the mechanical properties to model fluid-like behavior when the fragmentation status represents fragmented.

7. The method of the preceding claim, including determining, with the one or more physical processors, the fragmentation status of one or more of the discrete particles in the pre-blast solid material based on a fragmentation input, e g., a user-defined volume representing pre-fragmented material that are affected by the eigenerosion model in the first numerical determination / calculation to account for previous blasts that would have already broken part of the solid material.

8. The method of any one or more of the preceding claims, wherein the mechanical properties include a friction angle and friction hardening coefficients, and the numerically determining / calculating includes determining / calculating the velocities and updated 3D locations using the friction angle and the friction hardening when the fragmentation status represents fragmented.

9. The method of any one or more of the preceding claims, including: increasing, with the one or more physical processors, a particle volume represented by each discrete particle to represent swelling of the solid material due to the explosive blasts, including none or more of the following: increasing the particle volume based on its updated velocity, increasing the particle volume based on its updated spin tensor; and / or increasing the particle volume based on its determinant of rate of deformation tensor.

10. The method any one or more of the preceding claims, including removing, with the one or more physical processors, the explosive particles from the discrete particles after a pre-defined duration to emulate gaseous release into the air through cracks in the solid material.

11. The method any one or more of the preceding claims, including representing the explosive particles with explosive material properties before and during ignition.

12. The method any one or more of the preceding claims, including representing the explosive particles in their gaseous states by a Jones-Wilkins-Lee (JWL) model.

13. The method of any one or more of the preceding claims, wherein the 3D grid includes a grid spacing between 0.1 meters and 10 meters in each direction, e g., substantially 1 meter.

14. The method of any one or more of the preceding claims, wherein the time steps each have a duration of substantially 10 microsecond to 1000 microseconds.

15. The method of any one or more of the preceding claims, wherein the numerically determining / calculating the updated velocities includes explicit integration of forces.

16. The method of any one or more of the preceding claims, including defining values of the 3D locations and the mechanical properties by discretizing, with the one or more physical processors, a block-model in the form of a set of non-overlapping cuboids representing the pre-blast solid material with material attributes attached to each cuboid.

17. The method of any one or more of the preceding claims, including automatically extending, with the one or more physical processors, the block-model using null blocks based on a user-defined floor and / or a measured pre-blast topography of the blast area to form the block-model to discretize18. The method of any one or more of the preceding claims, including defining the 3D blast locations by discretizing, with the one or more physical processors, one or more explosive volumes in respective blast holes in the pre-blast solid material, including defining the 3D blast locations by a random / quasi-random uniform distribution in the or each explosive volume.

19. The method of any one or more of the preceding claims, including representing an unbreakable floor under the pre-blast solid material as a heightfield, and using the heightfield as a boundary condition in the repeatedly numerically determining.

20. The method of the preceding claim, including repeatedly, with the one or more physical processors and at one or more of the time steps, one or more of: reducing / eliminating the velocity values of the discrete particles below the unbreakable floor; reducing / eliminating downward velocity values of discrete particles adjacent to the unbreakable floor; andreducing / eliminating tangential velocity values of discrete particles adjacent to the unbreakable floor depending on the normal velocity and the distance of the discrete particle from the unbreakable floor.

21. The method of any one or more of the preceding claims, including fitting, with the one or more physical processors, the discrete particles to a measured post-blast topography, including by: a. applying an attraction force to the discrete particles to move the discrete particles into a post-blast volume of the measured post-blast topography; and / or b. applying a suction force to the discrete particles to move the discrete particles to fill the post-blast volume of the measured post-blast topography.

22. Computer-readable storage having stored thereon computer-readable instructions configured to cause a system that includes the one or more physical processors in the form of a CPU and one or more GPUs to perform the method of any one of claims 1 to 21.

23. A system including the one or more physical processors and the computer-readable storage of claim 22.

24. A mining method, including: a the method of any one of claims 1 to 21 for predicting solid material dislodgement due to explosive blasting; and b. physically measuring the pre-blast solid material to determine the one or more mechanical properties of the solid material before the explosive blasting, and / or physically excavating and / or routing the solid material after the explosive blasting based on the determined post-blast locations of the solid material.

Images