The Material Point Method (MPM) is widely used to analyse coupled (solid-water) problems under large deformations/displacements. However, if not addressed carefully, MPM u-p formulations for poromechanics can be affected by two major sources of instability. Firstly, inf-sup condition violation can arise when the spaces for the displacement and pressure fields are not chosen correctly, resulting in an unstable pressure field when the equations are monolithically solved. Secondly, the intrinsic nature of particle-based discretisation makes the MPM an unfitted mesh-based method, which can affect the system's condition number and solvability, particularly when background mesh elements are poorly populated. This work proposes a solution to both problems. The inf-sup condition is avoided using two overlapping meshes, a coarser one for the pressure and a finer one for the displacement. This approach does not require stabilisation of the primary equations since it is stable by design and is particularly valuable for low-order shape functions. As for the system's poor condition number, a face ghost penalisation method is added to both the primary equations, which constitutes a novelty in the context of MPM mixed formulations. This study frequently makes use of the theories of functional analysis or the unfitted Finite Element Method (FEM). Although these theories may not directly apply to the MPM, they provide a robust and logical basis for the research. These rationales are further supported by four numerical examples, which encompass both elastic and elasto-plastic cases and drained and undrained conditions.

The Material Point Method (MPM) has been shown to be an effective approach for analysing large deformation processes across a range of physical problems. However, the method suffers from a number of spurious artefacts, such as a widely documented cell crossing instability, which can be mitigated by adopting basis functions with higher order continuity. The larger stencil of these basis functions exacerbate a less widely discussed issue - small cuts. The small cut issue is linked to the arbitrary interaction between the physical body and the background mesh that is used to assemble and solve the governing equations in the MPM. There is the potential for degrees of freedom near the boundary of the body to have very small contributions from material points, which causes two problems: (i) artificially large accelerations/displacements at the boundary and (ii) ill conditioning of the global linear system. This paper provides a new mesh Aggregated MPM, or AgMPM, that mitigates the small cut issue by forming aggregated elements, tying the ill-behaved degrees of freedom to well posed interior elements. Implicit quasi-static and explicit dynamic formulations are provided and demonstrated through a series of numerical examples. The approach does not introduce any new numerical parameters and can be applied to implementations that adopt a lumped mass matrix. Aggregation is shown to significantly improve the stability of implicit implementations of the MPM, often at a lower computational cost compared to standard, non-aggregated, implementations. The technique improves the energy conservation and the stress field of explicit dynamic MPMs.

Accurate and robust modelling of large deformation three dimensional contact interaction is an important area of engineering, but it is also challenging from a computational mechanics perspective. This is particularly the case when there is significant interpenetration and evolution of the contact surfaces, such as the case of a relatively rigid body interacting with a highly deformable body. The numerical challenges come from several non-linear sources: large deformation mechanics, history dependent material behaviour and slip/stick frictional contact. In this paper the Material Point Method (MPM) is adopted to represent the deformable material, combined with a discretised rigid body which provides an accurate representation of the contact surface. The three dimensional interaction between the bodies is detected though the use of domains associated with each material point. This provides a general and consistent representation of the extent of the deformable body without introducing boundary representation in the material point method. The dynamic governing equations allows the trajectory of the rigid body to evolve based on the interaction with the deformable body and the governing equations are solved within an efficient implicit framework. The performance of the method is demonstrated on a number of benchmark problems with analytical solutions. The method is also applied to the specific case of soil-structure interaction, using geotechnical centrifuge experimental data that confirms the veracity of the proposed approach.

The ability to model rigid body interaction with highly deformable solids is a very useful tool in geoengineering, including the modelling of drag anchors on seabeds and seabed ploughing [7, 1]. However, these simulations entail several numerical challenges, such as modelling frictional contact, and incorporating inertia forces for analyses whose simulated time is considerable. Here the implicit Generalised Interpolation Material Point Method (GIMPM) is adopted to model the highly deformable solid, whilst a rigid body is used to model the significantly stiffer engineering object, such as an anchor. The whole system is integrated in time with the Newmark method with interaction between the two bodies occurring through a normal penalty contact and a penalty enforced Coulomb stick-slip friction law.

Modelling the interaction between rigid and deformable bodies holds significant relevance in geotechnical engineering, particularly in scenarios involving stiff engineering objects interacting with highly deformable material such as soil. These processes are challenging due to the combined nonlinear mechanisms including large deformation, elasto-plasticity, and contact with friction. For highly deformable material, the Material Point Method is a natural choice over the Finite Element Method due to its ability to handle large deformations without remeshing by carrying material information at points. This paper uses the Implicit General Interpolation Material Point Method (GIMPM) to demonstrate a new approach for modelling this type of interaction, and exploits the GIMPM's inherent definition of the boundary of a deformable domain to formulate a consistent contact formulation, negating the need for boundary reconstruction. The formulation is demonstrated through validations and comparisons to alternative methods for simulating contact. The combination of the contact formulation with an implicit framework is shown to be an efficient method for modelling geotechnical problems. The proposed method exhibits optimal convergence for contact problems, accurately captures stick-slip Coulomb friction, and ensures consistent stress fields at the contact surface of a rigid body.

The kinematic behaviour of drag embedment anchors has become a recent research focus due to the increase in offshore renewable energy devices. This is due to their potential use as an anchoring system for future floating wind applications, in addition to the need to understand their penetration behaviour as a part of the Cable Burial Risk Assessment. Studies on the behaviour of anchors typically consist of field scale or model centrifuge tests, where such facilities are not readily available to all and can result in significant cost. In addition to this, measuring the load–penetration behaviour of an anchor has proven to be a significant challenge, as any contact-based methods are likely to influence the penetration behaviour of the anchor. In this paper, a novel wireless method of recording the inclination of the anchor and calculating the penetration depth is presented. A comparison of the penetration behaviour of a Class F (AC-14) anchor has been made in sand using centrifuge and 1g model scale testing. The results indicate that the 1g testing can match the behaviour of the anchor testing in the centrifuge in terms of both the position of the anchor and its orientation during the dragging event.

Phase field approaches are an increasingly popular method for modelling complex fracture problems and have been applied to a number of real-world settings. In some applications, pressure forces, or more generally traction terms, must be considered on the crack faces. However, application of appropriate boundary conditions to represent these tractions is non-trivial, since phase field models do not include a direct description of the fracture surface due to their diffuse nature. This paper summarises one- and two-domain approaches to including immersed traction boundary conditions and states the authors’ intention to implement and evaluate these methods in an hp-adaptive discontinuous Galerkin finite element framework.

The phase field method is becoming the de facto choice for the numerical analysis of complex problems that involve multiple initiating, propagating, interacting, branching and merging fractures. However, within the context of finite element modelling, the method requires a fine mesh in regions where fractures will propagate, in order to capture sharp variations in the phase field representing the fractured/damaged regions. This means that the method can become computationally expensive when the fracture propagation paths are not known a priori. This paper presents a 2D hp-adaptive discontinuous Galerkin finite element method for phase field fracture that includes a posteriori error estimators for both the elasticity and phase field equations, which drive mesh adaptivity for static and propagating fractures. This combination means that it is possible to be reliably and efficiently solve phase field fracture problems with arbitrary initial meshes, irrespective of the initial geometry or loading conditions. This ability is demonstrated on several example problems, which are solved using a light-BFGS (Broyden–Fletcher–Goldfarb–Shanno) quasi-Newton algorithm. The examples highlight the importance of driving mesh adaptivity using both the elasticity and phase field errors for physically meaningful, yet computationally tractable, results. They also reveal the importance of including p-refinement, which is typically not included in existing phase field literature. The above features provide a powerful and general tool for modelling fracture propagation with controlled errors and degree-of-freedom optimised meshes.

The Carbon Trust (2015) 'Cable Burial Risk Assessment (CBRA) Methodology' document is widely used in the offshore subsea cable industry to define the cable burial Depth of Lowering (DoL). To-date published work on anchor penetration depths has focused on single homogeneous soil units, offering limited information on the response of different soil layering combinations, and associated contrasting geotechnical properties between soil units. By interrogating >11,000 shallow cores from the entire UK North Sea area, we demonstrate that 'layered' soil combinations (e.g., 'sand over clay') are statistically common across the North Sea study area. The results also highlight the importance of updating current CBRA approaches to include 'layered' soils, and associated changes in geotechnical properties (e.g., strength and density) between single and layered soil units. In addition, we collated geotechnical data for input into physical and numerical modelling undertaken by the University of Dundee and Durham University respectively (see Sharif et al., 2023; Bird et al., 2023 a, b), to assess the implications for the current CBRA Methodology. Ultimately the goal is to create a new CPT-based tool for better constraining the DoL, as part of the EPSRC research grant 'Offshore Cable Burial: How deep is deep enough?'.

The complexity of the physics of rock blasting is a longstanding modelling challenge. This work presents in detail a three-dimensional, material non-linear finite element based model for wave propagation, combined with a postprocessing procedure to determine the fracture intensity caused by blasting. The rock is described with the Johnson–Holmquist-2 constitutive model, an elastoplastic-damage model designed for brittle materials undergoing high strain rates and high pressures and fracturing; it is also combined with an instantaneous tensile failure model. Additionally, material heterogeneity is introduced into the model through variation of the material properties at the element level, ensuring jumps in strain. A detailed algorithm for the combined Johnson–Holmquist-2 and tensile failure model is presented and is demonstrated to be energy-conserving, and is complemented with an open-source MATLABTM implementation of the model. A range of sub-scale numerical experiments are performed to validate the modelling and postprocessing procedures, and a range of materials, explosive waves and geometries are considered to demonstrate the model's predictive capability quantitatively and qualitatively for fracture intensity. Fracture intensities on 2D planes and 3D volumes are presented. The mesh dependence of the method is explored, demonstrating that mesh density changes maintain similar results and improve with increasing mesh quality. Damage patterns in simulations are self-organising, and form thin, planar, fracture-like structures that closely match the observed fractures in the experiments. The presented model is an advancement in realism for continuum modelling of blasts as it enables fully three-dimensional wave interaction, handles damage due to both compression and tension, and relies only on measurable material properties.

Open source codes are a key ingredient to greater research integrity and accountability in computational science and engineering. However, many of these codes have not been developed with modification of the base code as their primary consideration. Existing codes may provide an environment for researchers to quickly test out their ideas under different physical conditions in a high level way but they are not always ideal for those interested in the development of numerical methods. The majority of existing open source discontinuous Galerkin finite element codes are written in C++ and there is a significant learning curve for junior researchers to adopt, un-derstand and modify the underlying code/routines. This paper presents an open source hp-adaptive discontin-uous Galerkin finite element code written in MATLAB that has been explicitly designed to make it easy for users, especially MSc/PhD-level researchers, to understand the method and implement new ideas within the core code. Although the code is focused on solving problems in linear elasticity, it is straightforward to modify it to solve other physical equations.

ABSTRACT: Simultaneous detonation of charges in closely spaced boreholes is a commonly used blasting technique for frag-mentation and construction. Blasting is a challenging phenomenon to model, due to the complexity of the mechanical deformation, fracturing, and fragmentation, and the spatial scales involved. Blast wave models must consider the possibility of constructive interference between separate waves in three-dimensions. A three-dimensional finite element method is used to study the induced damage in a general hard rock tunnel blast setting. The Johnson Holmquist-2 elastoplastic-damage model is used to quantify shear and tensile failure. In simulations with two blastholes separated by up to one meter, damage patterns emanating from boreholes interact to form self-organizing ‘fracture-like’ structures. Mechanical interaction between the two blasts is a function of the input charge wave properties, blasthole separation, and distance along the charge, with high concentrations of damage at the free boundary representing the tunnel wall. Constructive interference between the two blast waves is not shown to directly induce additional damage zones, and instead, interaction results from the overlap and slight extension of each blasthole’s damage zone towards the other. 1 INTRODUCTION Drilling and blasting continues to be a key method for mining and excavation (Lu et al., 2012). Many blasting applications, such as bench blasting or tunnel excavation, involve drilling and blasting multiple holes in close proximity to one another in order to generate overlapping fractured zones (Cho and Kaneko, 2004). Efficient use of explosives reduces costs, improves sustainability, and ensures that the extent of the excavation disturbed zone is minimized, the latter being an important constraint for the construction of radioactive waste disposal facilities (e.g., Tsang et al., 2005; Kwon et al., 2009) and many civil engineering projects (e.g., Sharafat et al., 2019). Due to the complex physics of blasting and rock fragmentation, a fully mathematical description of the process is not possible, and blast operation designs often make use of empirical formulas (Liqing and Katsabanis, 1997). In particular, the dynamic weakening (damaging) of rock introduces challenges for models, and upscaling of damage is necessary to avoid resolving the complex geometry of fracturing and fragmentation. Numerical approaches are therefore used widely to understand the process of blasting (e.g., Liqing and Katsabanis, 1997; Yilmaz and Unlu, 2013; Yang et al., 2015; Hajibagherpour et al., 2020; Ji et al., 2021; Pu et al., 2021).

This article presents a hpr-adaptive crack propagation method for highly accurate 2D crack propagation paths which requires no a priori knowledge of the tip solution. The propagation method is designed to be simple to implement, only hr-adaptivity is required, with the propagation step size independent of the initial mesh allowing users to obtain high fidelity crack path predictions for domains containing multiple cracks propagating at different rates. The proposed method also includes a crack path derefinement scheme, where elements away from the crack tip are derefined whilst elements close to the crack tip are small so capture the fidelity of the crack path. The result is that the propagation of cracks over an increasingly larger distances has negligible increased computational effort and effect on the propagation path prediction. The linear elastic problem is solved using the hp discontinuous Galerkin symmetric interior penalty finite element method, which is post-processed to obtain the configurational force at each tip to a user defined accuracy. Several numerical examples are used to demonstrate the accuracy, efficiency, and capability of the method. Due to the method's high accuracy crack path solutions of benchmark problems that are prolifically used in the literature are challenged.

Engineers require accurate determination of the configurational force at the crack tip for fracture fatigue analysis and accurate crack propagation. However, obtain- ing highly accurate crack tip configuration force values is challenging with numer- ical methods requiring knowledge of the stress field around the crack tip a priori. In this thesis, the symmetric interior penalty discontinuous Galerkin finite element method is combined with a residual based a posteriori error estimator which drives a hp-adaptive mesh refinement scheme to determine accurate solutions of the stress field about about the crack. This facilitates the development of a novel method to calculate the crack tip configurational force that is accurate, requires no a priori knowledge of the stress field about the crack tip with, its error bound by an error estimator which is calculated a posteriori. Benchmark values of the crack tip con- figurational force are presented for problems containing multiple mixed mode cracks in both isotropic and anisotropic materials. Additionally, the hp-adaptivity is com- bined with a mathematical analysis of the stress field at the crack tip to critique the convergence and limitations of other methods in the literature to calculate the crack tip configurational force. Two methods for staggered quasi-static crack prop- agation are also presented. An rp-adaptive method which is simple to implement and computationally inexpensive, element edges aligned with the crack propagation path with the exploitation of the discontinuous Galerkin edge sti↵ness terms exist- ing along element interfaces to propagate a crack. The second method is denoted the hpr-adaptive method which combines the accurate computation of the crack tip configuration force with r-adaptivity to produce a computationally expensive but accurate method to propagate multiple cracks simultaneously. Further, for indeter- minate systems, an average boundary condition that restrains rigid body motion and rotation is introduced to make the system determinate.

Engineers require accurate determination of the configurational force at the crack tip, and corresponding stress intensity factors, for fracture fatigue analysis and accurate crack propagation. However, obtaining highly accurate crack tip configuration force values is challenging with methods requiring knowledge of the stress field around the crack tip a priori. This paper proposes a method which aims to remove the necessity of knowing the stress field a priori whilst producing very accurate values of the configurational force at a static crack tip. The proposed method is demonstrated to be path independent and is combined with a robust a posteriori residual error estimator which is indicative of the accuracy of the configurational force calculation. This makes it possible to generate accurate values for the configurational force acting both perpendicular and parallel to the crack edges. Accuracies are achieved which are at least 10(4) times more accurate than other numerical methods which make no assumption about the local tip stress field. Therefore accurate benchmarks are determined in this paper for inclined, split and tree crack problems. In addition the new method is shown to obtain very similar values for the configurational force compared to results obtained using other methods which require knowledge of the stress field at the crack tip. The techniques presented in this paper open the door to configurational force-based methods being used for fatigue analysis.

MATLAB is adept at the development of concise Finite Element (FE) routines, however it is commonly perceived to be too inefficient for high fidelity analysis. This paper aims to challenge this preconception by presenting two optimised FE codes for both continuous Galerkin (CG) and discontinuous Galerkin (DG) methods. Although this has previously been achieved for linear-elastic problems, no such optimised MATLAB script currently exists, which includes the effects of material non-linearity. To incorporate these elasto-plastic effects, the externally applied load is split into a discrete number of loadsteps. Equilibrium is determined at each loadstep between the externally applied load and the arising internal forces using the Newton–Raphson method. The optimisation of the scripts is primarily achieved using vectorised blocking algorithms, which minimise RAM-to-cache overheads and maximise cache reuse. The optimised codes yielded maximum speed gains of ×25.7 and ×10.1 when compared to the corresponding unoptimised scripts, for CG and DG respectively. It was identified that with increasing refinement of the mesh, the solver time begins to dominate the overall simulation time. This bottleneck has a greater disadvantage on the DG code, predominantly due the asymmetric nature of the global stiffness matrix. The implementation of an efficient solver would see further improvement to the overall run times, particularly for large problems.

This paper presents a novel combination of two numerical techniques to produce a method for solving fracture mechanics problems. A weak form meshless method, the cracking particles method, forms the basis of the mechanical model while crack propagation direction is calculated using configurational forces. The combined method is presented here for 2D quasi-brittle crack propagation. The configurational force approach has the advantage that it provides a prediction of the crack propagation direction which does not require decomposition of the stress and displacement fields for mixed-mode crack problems. The use of a meshless method removes the need for remeshing and it is therefore eminently suitable for multiple crack problems. The paper includes a discussion on the configurational force calculations via contour integration and domain integration and results are presented that show both approaches to be path independent when the integrations over the two crack surfaces cancel out, with domain integration generally providing better accuracy than contour integration. The contribution from the crack surfaces to the configurational force is discussed, and shown to have little influence on the final result while being easily affected by the oscillations around the crack tip. In addition, the relationship between the configurational force and the J-integral is explained. The proposed method is demonstrated on several examples, including multiple crack propagation, where good agreements with results from the literature are obtained.