Coupling of IGA and Peridynamics for Air-Blast Fluid-Structure Interaction Using an Immersed Approach

We present a novel formulation based on an immersed coupling of Isogeometric Analysis (IGA) and Peridynamics (PD) for the simulation of fluid-structure interaction (FSI) phenomena for air blast. We aim to develop a practical computational framework that is capable of capturing the mechanics of air blast coupled to solids and structures that undergo large, inelastic deformations with extreme damage and fragmentation. An immersed technique is used, which involves an a priori monolithic FSI formulation with the implicit detection of the fluid-structure interface and without limitations on the solid domain motion. The coupled weak forms of the fluid and structural mechanics equations are solved on the background mesh. Correspondence-based PD is used to model the meshfree solid in the foreground domain. We employ the Non-Uniform Rational B-Splines (NURBS) IGA functions in the background and the Reproducing Kernel Particle Method (RKPM) functions for the PD solid in the foreground. We feel that the combination of these numerical tools is particularly attractive for the problem class of interest due to the higher-order accuracy and smoothness of IGA and RKPM, the benefits of using immersed methodology in handling the fluid-structure coupling, and the capabilities of PD in simulating fracture and fragmentation scenarios. Numerical examples are provided to illustrate the performance of the proposed air-blast FSI framework.


Introduction
In [1,2], the authors developed an immersed fluid-structure interaction (FSI) formulation for air blast. The formulation made use of the weak forms of the fluid and structural mechanics equations and two meshes, background and foreground.
The paper is outlined as follows. In Section 2, the continuum-and discrete-level FSI formulations, which include the fluid and solid subproblems in the weak form and their coupling, are presented. In Section 3, numerical examples are presented to verify the accuracy of the coupled IGA-PD FSI formulation and to demonstrate the ability of the new formulation to simulate fracture and fragmentation in brittle and ductile solid structures subjected to blast loading. Section 4 makes concluding remarks and outlines future research directions to address the remaining numerical challenges.

Mathematical formulation
In this section, we develop an immersed FSI formulation that accommodates compressible flow and an inelastic solid or structural model coming from a standalone weak formulation and its discretization. We first provide a general continuum variational framework and then specialize to the coupling of IGA on the background and PD on the foreground grids at the discrete level. The resulting methodology resembles the classical Immersed Boundary Method [36] and a more recent Immersed Finite Element Method [37], but without the use of ad hoc smoothed delta functions and with the force transfer motivated by variational arguments. In what follows, all the vectors are column vectors unless otherwise noted. Bold symbols denote tensors of rank one (i.e., vectors) or higher. Underscored terms indicate the PD field variables at the bond level.

Weak form of the coupled FSI problem
We take, as the starting point, the coupled FSI formulation proposed in [1], where the fluid mechanics part of the problem is governed by the Navier-Stokes equations for compressible flow and the structural mechanics part is modeled as an isothermal large-deformation inelastic solid. Let Ω denote the combined fluid and solid domain, and let Ω f and Ω s denote the individual fluid and solid subdomains in the spatial configuration, such that Ω f Ω s = Ω and Ω f Ω s = ∅. Both the fluid and solid problems are formulated in a variational setting, where we define the following semilinear forms and linear functionals that comprise the weak forms of the fluid and solid subproblems Here, Y denotes a set of pressure-primitive variables [38,39], where p is the pressure, v is the material-particle velocity vector, and T is the temperature. Y and W, the vector-valued trial and test functions, respectively, are the members of S and V, the corresponding trial and test function spaces, respectively, defined on all of Ω. The superscripts f and s correspond to the fluid and solid quantities, respectively. Γ f H and Γ s H are the subsets of the fluid-and solid-domain boundaries where natural boundary conditions are imposed, and H f and H s contain the prescribed values of the natural boundary conditions. A 0 and A i are the so-called Euler jacobian matrices, F p i , F d i , and F σ i are the pressure, viscous/thermal, and stress fluxes, respectively, and S is the volume source. (See [1] and references therein for the details.) The subscript ω on the semilinear forms and linear functionals denotes the domain of integration, comma denotes partial differentiation, and i = 1, . . . , d, where d = 2, 3 is the space dimension.
With the above definitions, following [1] the coupled FSI formulation may be written as: Find Y ∈ S, such that ∀W ∈ V, where the integration over the fluid mechanics domain is replaced by the integration over the combined domain minus that over the solid domain. This form of the coupled problem statement is convenient for the application of an immersed approach to the discretization of the coupled FSI equations (see, e.g., [40]). It was shown in the original reference [1] that the above coupled formulation satisfies a.) The fluid and solid governing equations in their respective subdomains; and b.) The compatibility of the kinematics and tractions at the fluid-solid interface. The formulation given by Equation (8) may be employed directly in the case the weak form of the solid problem can be written in terms of the test and trial functions coming from the spaces V and S defined over Ω. However, in order to have more flexibility for the solid problem formulation and discretization, the above formulation needs to be modified, which we do as follows. First, we assume that the solid formulation is given in the weak form that is defined over a different set of test and trial functions 2 : FindỸ ∈S, such that ∀W ∈Ṽ, We also assume that we can develop a relationship between the functions in S andS, and the functions in V andṼ, by defining the mappingsỸ = ΠY In the discrete setting these mappings may correspond to interpolation or projection, and will be detailed in what follows. These assumptions allow us to re-write the solid problem as: Find Y ∈ S, such that ∀W ∈ V, As a consequence, the coupled FSI formulation given by Equation (8) may be modified as: Find Y ∈ S, such that ∀W ∈ V, which, as will be shown in what follows, enables us to flexibly incorporate a much larger class of solid and structural formulations in our immersed FSI framework.

Correspondence-based PD framework for solid modeling
In this section we focus on the PD formulation of the solid. The rate form of the energy balance law on Ω s may be expressed as [10,41]: Here, the first and second terms on the left hand side denote the rate of the total kinetic and internal energy, respectively, while the right hand side gives the total supplied power. Here, ρ s is the solid mass density in the current configuration, v is the velocity field, and s is the volume source term. The overdot symbol indicates the time derivative. In the PD formulation the strain-power densityU at a material point x is regularized as [42,43]U Here, H(x) is the family set (also known as the horizon) of x defined as where δ is the horizon size, x − x denotes a PD bond between x and x in the horizon, and the strain-power densityU x − x is given as a bond-level quantity. The normalized weighting function α x − x is constrained to satisfy the equation For simplicity of notation, to distinguish the material-point-and bond-level quantities, the bondassociated fields are marked using underscores in the remainder of this paper. For example, Equation (15) may be re-written asU whereU =U(x), α = α x − x ,U =U x − x , and H = H(x). The bond-associated strain-power density is given byU where σ and L are the bond-associated power-conjugate Cauchy stress and velocity gradient tensors, respectively. The total stress power, using Equations (14), (18) and (19), may be re-expressed as where T = T x − x is the so-called PD force state (with units of force per unit volume squared). The PD velocity state v is defined as In the correspondence PD framework, the spatial gradients are computed using an integral form. Using the computed velocity gradient, a classical constitutive law is employed to evaluate the Cauchy stress tensor. As a result, the way L is calculated determines the relation between the force state and Cauchy stress at the bond level, which we detail in the Appendix. Re-writing Equation (14) using Equations (20) and (21) yields where T = T x − x . Noting that Equation (22) must hold for all choices of v enables us to recast the PD linear-momentum balance equation in the weak form given by Equation (12) with the following definitions of the semi-linear and linear forms: and It is important to note that no higher regularity than L 2 is required forṽ andw in order for the weak form of the PD linear-momentum balance equation to be well defined. Also note that the solid is assumed to be isothermal and the mass balance is satisfied in the Lagrangian description at each material point, resulting in the PD formulation that only makes use of the momentum-balance PDE written in terms of the velocity unknowns. The extension to a thermally coupled solid does not present a conceptual difficulty; however, it is not pursued here.

Discrete form of the coupled FSI problem
To develop a discrete version of the FSI formulation, we first define S h ⊂ S and V h ⊂ V, the background-domain finite-dimensional trial-and test-function spaces, respectively. The discrete trial and test functions may be expressed as and where Y B and W A are the vector-valued control-point degrees of freedom (DOFs) and weights, respectively, N (x)'s in this work are assumed to be the B-Spline basis functions defined everywhere in Ω, and N cp is the dimension of the B-Spline space. Note that all the components of the trial and test-function vectors are approximated using the same basis functions. We also define the finite-dimensional trial and test function spacesS andṼ, respectively, for the PD solid. The discrete trial and test functions may be expressed as where the solid domain is represented using a finite number N mp of material points or PD nodes,Ỹ Q andW P are the nodal DOFs and weights, respectively, and χ P (x) is a characteristic function of a PD node P that attains a value of unity at x P ∈ Ω s , the spatial location of the PD node, is zero at all other PD nodes, and satisfies Ω s where V P is the local volume of the PD node. Note that the exact shape of the characteristic functions is not necessary for our purposes because, ultimately, we will make use of nodal quadrature to evaluate the integrals in the PD formulation. We now define the mapping Π between the discrete background and foreground functions using nodal interpolation. Namely, we constrain the PD nodal DOFs and weights to the background test and trial functions, respectively, asỸ andW which yieldsỸ andW and defines the mapping Π.
With the above definitions of the test and trial functions, the spatially-discretized immersed FSI formulation may be stated as: Note that, in the space-discrete case, we introduce the SUPG stabilization (B st [48][49][50][51] operators to the compressible flow formulation to address the convective instability of the Galerkin technique in the regime of convection dominance and to provide additional dissipation in the shock regions. Both the stabilization and shock-capturing operators are proportional to the strong-form PDE residual of the compressible-flow Navier-Stokes equations, which retains the method consistency. More details on the definition of these operators that are employed in the present formulation of the compressible flow may be found in [38,52]. Introducing the definitions of the test and trial functions from Equations (26)-(27), (28)- (29), and (31)-(34) into the space-discrete weak form given by Equation (35) gives where I = 1, 2, . . . , d + 2 is the nodal or control-point DOF index and e I is a Cartesian basis vector of dimension (d + 2) with one in the I th entry and zero in the remaining entries. Because Equation (36) holds for all background-domain discrete weights W AI , we arrive at the following vector form of the semi-discrete FSI problem: Find the background-domain control variables where and Equation (37) states that at each control point the generalized force vectors coming from the fluid and solid parts of the problem balance; Equation (40) states that the PD solid nodal forces may be computed directly from the PD formulation on the foreground mesh provided the kinematic data can be interpolated on the PD nodes from the background solution field; and Equation (39) defines a consistent transformation of the force vector from the foreground to the background grid.
Introducing the definitions of the semi-linear and linear forms of the PD formulation given by Equations (23)-(25) into Equation (40) yields an explicit expression for the PD force vector on the foreground grid:R where H P is the horizon of the PD node located at x P and the integral is evaluated using nodal quadrature. Note that in the present model the PD solid will only contribute to the linear-momentum equation balance of the discrete FSI formulation given by Equation (37).
Remark 2.1. The formulation developed in this section presents a minimally intrusive approach to immersed FSI, and opens the door to using not just PD, but any other suitable discrete formulation of solid and structural mechanics as part of the FSI coupling framework. All that is required is the ability to develop a mapping of the background solution to the DOFs employed by the solid and structural mechanics formulation.

Time integration and the lumped mass matrix
The semi-discrete FSI equations are integrated in time using the lumped-mass explicit Generalized-α predictor-multi-corrector procedure [53] adopted for the immersed FSI formulation and detailed in [1]. This approach requires the computation of a lumped mass matrix for the coupled system. A consistent mass matrix for the coupled FSI problem is obtained by adding the contributions from the fluid and solid subdomains on the background mesh. Denoting byM s P Q the coefficients of the consistent mass matrix from the PD foreground mesh and using the relations developed in this section, one can show that its counterpart on the background mesh is obtained through the transformation The diagonal entries of the lumped mass L s A may be obtained using a classical row-sum technique, namely whereL s P are the diagonal entries of the lumped mass on the PD foreground mesh given bỹ Equation (43) states the diagonal entries of the lumped mass matrix transform in the same way from the foreground to the background mesh as the nodal forces (see Equation (39)). Because the present PD formulation only solves the momentum balance equation, L s A 's are added only to the diagonal entries 2, . . . , d + 1 of the (d + 2) × (d + 2) control-point block of the background lumped mass matrix.

Numerical examples
Three 2D numerical simulations are provided here to showcase the capabilities of the proposed IGA-PD framework in blast loading and fragmentation events. The first example serves as a verification case for the developed formulation in solving problems involving large deformation and plasticity. The next two cases are demonstrative examples that include fracture propagation and fragmentation in brittle and ductile solids subjected to blast loading. In all the presented computations, C 1 -continuous quadratic NURBS functions are used for the background solution. Unless otherwise noted, RK functions with quadratic consistency, rectangular support, and the BA stabilization technique [34,54] are employed in the PD formulation. The PD support size δ is chosen with respect to the discretization size ∆x and the order of the kernel function n is chosen such that δ/∆x = n + 1 [55]. The air properties are used for the fluid with constant viscocity µ = 1.81 × 10 −5 kg/(m s), Prandtl number 0.72, and adiabatic index γ = 1.4.

Chamber detonation
In this FSI example, a steel bar undergoes a detonation blast loading condition. The problem description is shown in Figure 1 The detonation condition is enforced by setting higher-than-ambient values of the pressure and temperature, i.e., we set p = 6.75 MPa and T = 1465 K in a semi-circular area with a radius of 6.1 mm, centered on the left wall. No-penetration and free-slip boundary conditions are prescribed at the chamber walls.
The IGA-PD simulations are compared with the reference immersed simulations in [1,2], which we refer to as IGA-IGA. Three different discretization levels are considered: coarse -Fluid: 40 × 40 elements; Solid: 60 × 30 elements (PD nodes); medium -Fluid: 80 × 80 elements; Solid: 120 × 60 elements (PD nodes); fine -Fluid: 120 × 120 elements; Solid: 180 × 90 elements (PD nodes). In the PD case, each foreground element is replaced by a meshfree node at its centroid with an equivalent volume. The time step used for the coarse, medium, and fine cases are 0.5 ms, 0.25 ms, and 0.166 ms, respectively.
The bond-associative, quadratic (n = 2) PD solid and the quadratic IGA solid are used in the computations. Air pressure at different time instants during the simulation is compared between the IGA-PD and reference simulations in Figure 2 using the finest mesh. The plastic strain contours in the final configuration of the specimen are shown in Figure 3. A good agreement between the results is obtained. Note the large deformation at the bar corners and the mushrooming at the left edge of the bar. Also note that the shock waves bounce off the right wall, impact the specimen, and leave a permanent indentation on the right edge of the bar.
The discretization effects are studied in this problem. The time history of the center-of-mass displacement of the bar, pressure at the detonation center (Point 1 in Figure 1), and pressure at the right-wall center (Point 2 in Figure 1) are shown in Figure 4 for different discretization levels of IGA-PD and IGA-IGA. The results indicate good convergence of the IGA-PD solution with mesh refinement, and a good agreement between the IGA-PD and IGA-IGA results.
We study the effects of the BA stabilization and order of accuracy on the PD solid behavior. The linear (n = 1) and quadratic (n = 2) PD variants with and without BA corrections are considered. The PD simulation results are compared with their quadratic IGA counterparts in Figure 5. The non-BA PD results are of poor quality due to the presence of unstable modes that deteriorate the solution and eventually lead to divergence of the simulation (around t = 0.25 ms here). The results of quadratic PD model equipped with BA stabilization are in good agreement with the reference solution (better than the linear + BA version), which confirms the conclusions of [55] that the combination of BA stabilization and higher-order corrections provides the best quality results in correspondence-based PD. In the remainder of this paper, we utilize the quadratic + BA model for the PD calculations.

Brittle material subjected to internal explosion
In this FSI example, we consider a hollow cylindrical block of elastic, brittle material subjected to internal detonation. As shown in Figure 6, the detonation is initiated at the center of the hollow cylinder with the inner radius of 7 cm and outer radius of 10 cm. The background domain is a 30 cm × 30 cm square enclosing the cylinder domain. The solid elastic material has Young's modulus E = 200 GPa, Poisson's ratio ν = 0.29, and initial density ρ s = 7870 kg/m 3 . To simulate brittle fracture, the von Mises stress failure criterion is used with the critical stress σ cr = 3 GPa. As noted in Appendix C, a PD bond is broken once its associated von Mises stress exceeds σ cr . Initially, the air is at rest with p = 0.1 MPa and T = 290 K. The detonation condition is enforced by setting the initial pressure p = 6.75 MPa and temperature T = 1465 K in a circular area with a radius of 3.5 cm, centered inside the hollow cylinder.
In this problem, the background domain and foreground solid are discretized into a uniform mesh and semi-uniform nodal setting (uniform along the θ-direction), respectively. We consider four discretization levels in this problem, D1-D4, with solid node spacings of 2.5 mm, 2 mm, 1.5 mm, and 1 mm, respectively. The fluid mesh size is set to three times of the solid node spacing in each case. The time step used for D1-D4 are 2.5 µs, 2 µs, 1.5 µs, and 1 µs, respectively. The results for D4 is shown in Figure 7, where the air speed is depicted on the background grid while the damage field is shown on the PD nodes. As expected, the impact of shock waves with the inner side of cylinder results in nucleation of several cracks at different locations along that side. Following the propagation of cracks through the solid structure, it is completely fractured and split into multiple fragments. The pressurized air passes through the empty space between fragments in the final snapshot in the figure. Because of the nature of the immersed methodology, fully-damaged solid nodes could still move within the computational domain without triggering mesh distortion issues that typically occur in a finite element foreground discretization. Although we do not present a comparison to experimental or other computational results, the predicted qualitative behavior appears to be physically reasonable. As shown in Figure 8, the fractures on the final configuration of the solid appear to converge toward a given state.

Ductile material subjected to internal explosion
The problem consists of a hollow square block of ductile material subjected to internal explosion. As shown in Figure 9, the detonation is initiated at the center of the hollow square with inner dimension of 10 cm and outer dimension of 16 cm. The background domain has the dimensions of 40 cm × 40 cm. Isotropic linear hardening rule is used for the solid ductile material, which has Young's modulus E = 200 GPa, Poisson's ratio ν = 0.29, yield stress σ Y = 0.4 GPa, hardening modulus H = 0.1 GPa, and initial density ρ = 7870 kg/m 3 . To simulate ductile fracture, a plasticity-driven failure approach given in Equation (C.2) is used with¯ P th = 0.18 and¯ P cr = 0.2. The air is at rest initially with p = 0.1 MPa and T = 290 K. The detonation condition is initiated by setting the pressure p = 6.75 MPa and temperature T = 1465 K in a circular area with a radius of 5 cm, centered inside the hollow square.
In this example, the background domain and foreground solid are discretized uniformly. Four discretization levels are considered here, D1-D4, with the solid node spacings of 2.5 mm, 2 mm, 1.5 mm, and 1 mm, respectively. In each case, the fluid mesh size is set to four times that of the solid node spacing. The time step used for the D1-D4 meshes is 1.8 µs, 1.2 µs, 0.9 µs, and 0.6 µs, respectively. Figure 10 shows the simulation results of D4 discretization, where the air speed is shown on the background mesh while the PD damage field is plotted on the foreground PD nodes. The stress singularity created by the re-entrant corner of the hollow square leads to crack nucleation in that location, as expected. The cracks initially propagate along the diagonal, then broaden and branch. Eventually, the solid structure is completely fractured and split into multiple fragments with the air occupying the newly open empty space. Note the permanent deformation of the fragments occurring due to the plasticity effects. As shown in Figure 11, the fragmented solid structure deformed configuration converges toward a unique shape, which is hard to achieve for ill-posed local fracture models due to mesh-dependency issues.

Discussion and conclusions
In this paper, a computational FSI framework for simulating air blast events based on an immersed IGA-PD approach is presented. B-Splines are used to construct the trial and test functions in the background domain for the coupled FSI problem while the foreground solid is discretized using correspondence-based PD, which is a meshfree methodology. The PD nodes are constrained to follow the background-grid kinematics, while the solid internal forces are computed directly using the PD discretization and then transferred to the background grid to complete the coupled discrete FSI formulation. This novel approach to immersed FSI has several advantages over the existing methods. In particular, using PD eliminates the need to explicitly track the crack interfaces and makes modeling fragmentation a relatively easy task. Likewise, we do not require to track the fluid-structure interface due to the use of an immersed method. Also, the use of higher-order and smooth B-Spline functions in the background and RKPM functions for the PD solid in the foreground results is higher-order accurate solutions for both the fluid and solid fields.
The results of the proposed IGA-PD formulation match well with the solutions using other conforming and immersed methods, which provides good verification of the present methodology. Although no experimental comparison is presented for the air-blast-structure interaction problems involving dynamic fracture, the IGA-PD simulations provide physically reasonable results that also show convergence under mesh refinement.
We also note some challenges using the current approach. While we obtain the desired strong coupling by constraining the solid to the background kinematics, simulating fragmentation becomes complicated since the smooth background discretion is not designed to appropriately capture the material discontinuities in the foreground solid solution. In practice, the damage zones, whose size scales with the background element size (see Figures 8 and 11), appear to be artificially thick leading to excessive structural damage. Refining the background mesh helps reduce the size of the damage zones, however, it considerably increases the computational costs. Relaxing the dependence of the foreground kinematics on the background solution, by leveraging such techniques as developed in [35], may be a way forward to address this challenge.
While we do not address this issue rigorously in the present work and leave it for the future developments, we explore a more practical way to simulate blast-induced fracture and fragmentation using the current formulation. In the blast events, the key factor that governs structural failure is the total impulse that is transmitted to the solid as a result of the explosion. Once enough loading is applied to the structure, it can nucleate cracks and potentially lead to complete disintegration. As a possible approach, we consider decoupling the foreground solution from the background discretization after damage starts growing in the solid material and letting the PD solid evolve on its own afterwards. We explore this idea using the last numerical example (ductile structure under blast loading). As shown in Figure 12, removing the constraint between the foreground and background fields once enough impulse is delivered to the structure and letting the standalone PD formulation govern the solid problem from that point on results in sharper cracks and thinner damage zones (bottom row in the figure) as opposed to the thicker damage regions in the fully coupled case (top row in the figure). While this approach involves some challenges in identifying when to decouple the solid, we observe that the unconstrained foreground solution is able to predict sharper fractures and better quality fragmentation. Future work will address this issue in a more rigorous way.

Appendix A. Semi-Lagrangian correspondence PD formulation of an inelastic solid
In the semi-Lagrangian PD model [34], the velocity gradient tensor is computed using only the current field variables and without a need for mapping to the reference configuration. As a prerequisite step for calculation of L, in this model, the point-level velocity gradient L is evaluated first using Here, Φ is a set of gradient kernel functions determined using a higher-order meshfree method that is asymptotically compatible with the local gradient operator [55]. We employ the reproducing kernel (RK) shape functions [56][57][58], as follows, to construct a meshfree differential operator equipped with an n-th order accuracy: where the so-called influence state ω is a weighting function dependent on the relative distance between material points with respect to the support size. We use a cubic B-spline kernel to define the radial influence function, i.e., whereη is defined as the current length of bond normalized with the support size, i.e., where x = x − x is called the PD position state. To satisfy the identity condition in Equation (17), the normalized weighting function α is defined as Example Appendix A.1. For a 2-dimensional system and a cubic kernel function (n = 3), the involved terms are: andΞ [3] = 1, 0, 0, 0, 0, 0, 0, 0, 0 0, 1, 0, 0, 0, 0, 0, 0, 0 . (A.10) To obtain a naturally-stabilized PD correspondence model, the bond-associative (BA) stabilization method of [34,54] is utilized to compute the velocity gradient at the bond level: where L = L(x). L is shown to be a second-order operator if L is equipped with a second or higher order of accuracy [55]. At the bond level, the velocity gradient is used to update the Cauchy stress using the classical theory, i.e.,ψ (∇L) = σ : L (A. 12) where σ is the power conjugate of L, and a classical rate-based constitutive law governs the strain-power density functionψ. The PD force state T x − x in this model reads [34] T = ω σ x ω |x| 2 + z Φ, (A. 13) where z(x) is given by z = H 0.5 ω + 0.5 ω ω σ I − x x |x| 2 dH (A.14) with the identity tensor I.

Appendix B. Solid constitutive equations
In this work, the rate-form constitutive model of [1] is adopted, which is based on the standard J2 flow theory with isotropic hardening [59] and thus suitable for modeling metal plasticity. Note that the constitutive relation is applied at the bond level. In this constitutive approach, the time history of the bond-associated velocity gradient derives the evolution of the bond-associated Cauchy stress. The rate-of-deformation tensor is additively decomposed into elastic and plastic parts, i.e., To maintain objectivity, the Jaumann rate of the Cauchy stress [60] is used, i.e., σ Y is the yield stress that depends on the bond-associated equivalent plastic strain¯ p , which evolves according to the associative plasticity rule [59], i.e., where H is the hardening modulus.
Appendix C. Bond-associative damage correspondence modeling The bond-associative failure model [34,61] is used here. In this correspondence-based approach, a continuum damage theory is utilized to compute the state of material damage in the solid material. The influence state, using this technique, is modified as where ω r is the conventional (radial) undamaged influence function as defined in Equation (A.3). ω D is the damage-dependent component of the influence state. A classical continuum damage model is incorporated to evolve the bond-associated damage parameter D ∈ [0, 1] based on the bond-associated internal variables such as strain, stress, stretch, and temperature. It is required that ω D (0) = 1 for an undamaged bond (D = 0) and ω D (1) = 0 for a fully damaged bond (D = 1). Damage irreversibility is enforced by constraining ω D as a non-increasing function of D.
In this work, to model the brittle fracture phenomenon, we break a bond once its associated equivalent (von Mises) stress exceeds a critical value. When the failure criterion is met, the bond breakage is enforced sharply to respect the immediate crack propagation nature of brittle materials.
For simulating the ductile fracture process, the gradual aspect of material degradation in malleable materials is considered. A plasticity-driven failure approach is adopted to decay the load-carrying capacity of a bond in its transition from the undamaged to fully-damaged state, i.e., where¯ P th and¯ P cr are the model parameters corresponding to the undamaged and completely damaged states, respectively.