A semi-Lagrangian reproducing kernel particle method with particle-based shock algorithm for explosive welding simulation

The explosive welding process is an extreme-deformation problem that involves shock waves, large plastic deformation, and fragmentation around the collision point, which are extremely challenging features to model for the traditional mesh-based methods. In this work, a particle-based Godunov shock algorithm under a semi-Lagrangian reproducing kernel particle method (SL-RKPM) is introduced into the volumetric strain energy to accurately embed the key shock physics in the absence of a mesh or grid, which is shown to also ensure the conservation of linear momentum. For kernel stability, a deformation-dependent anisotropic kernel support update algorithm is proposed, which is shown to capture excessive plastic flow and material separation. A quasi-conforming nodal integration is adopted to avoid the need of updating conforming cells which is tedious in extreme deformations. It is shown that the proposed formulation effectively captures shocks, jet formation, and smooth-to-wavy interface morphology transition with good agreement with experimental results.


Introduction
High-rate materials processing technologies often entail combination of shocks, extreme plastic deformation, and material separation in the form of jetting. Explosive welding is one of these manufacturing processes which cannot be effectively modeled by the Lagrangian or Arbitrary Lagrangian-Eulerian (ALE) finite element methods due to the jetting formation resulting from the excessive shock-induced plastic flow behind the collision point. Explosive welding is widely used as a type of high velocity impact welding to join a wide variety of similar and dissimilar metals with relatively small heat affected zones (HAZ). As shown in Fig. 1, the detonation process produces an oblique impact angle of the flyer plate relative to the base plate, resulting in material stagnation at the collision point. A metal jet is ejected due to the high-pressure field induced by the impact and material stagnation. The jetting is considered necessary for quality welding since it removes the non-metallic films (such as oxide films) and other contaminants initially attached to the metal surface [1]. It is observed, with certain processing parameters, that a wavy morphology is formed at the interface of the two bonded materials (see Fig. 2). Accurate capturing of interface wavy patterns and the metal jetting phenomena in the numerical simulation is critical to the reliable assessment of weldability in the explosive welding processes.
Key physics in explosive welding to be captured in numerical modeling includes strong shock wave propagation, fluid (plastic flow)-structure interaction, contact kinetics and kinematics between the flyer plate and base plate, material fragmentation (jetting), and extreme plastic deformation under high strain rates. Such complex processes make the simulation of explosive welding challenging. In many cases, the solid structure is modeled by a traditional Lagrangian method, such as the Lagrangian finite element method (FEM). However, this formulation breaks down due to the severe mesh distortion caused by large deformations at the impacting interfaces. Grignon et al. [5] modeled the explosive welding of aluminum plates using RAVEN, a twodimensional explicit Eulerian finite difference code. The ability to capture jetting was shown, but the interfacial waves observed in experiments was not. In more recent studies, both jetting and interfacial waves in magnetic pulse welding were captured by Eulerian FEM [6,7]. However, in the Eulerian description, the interface is smeared, and an interfacial tracking technique is required, which is non-trivial. An arbitrary Lagrangian-Eulerian (ALE) finite element simulation [8] was also performed. This approach alleviated the mesh distortion issue, but the jetting was not captured, and the magnitude of interfacial waves obtained were lower than experimental observations. To model jetting in Lagrangian or ALE FEM, one can consider non-physical erosion which is one common technique to model fragmentation [9], but this introduces errors in mass conservation and wave reflection on the new artificially generated free surfaces. Overall, jetting is very difficult to model using techniques such as conventional Lagrangian or ALE FEM, finite difference methods, and isogeometric analysis (IGA).
In meshfree methods, the approximation functions are constructed based on a set of scattered nodes instead of the mesh. Therefore, they do not suffer issues in mesh alignment, mesh entanglement, and time-consuming mesh refinement associated with mesh-based methods. In particular, fragmentation can be modeled naturally, governed by physical laws.
Smoothed-particle hydrodynamics (SPH) has been widely used for modeling impact welding [10][11][12][13][14][15] and captured jetting and interfacial wave formation. Efficient computation of shape functions is one of SPH's advantages. Easy accessibility in commercial software such as LS-DYNA and ABAQUS is another strength of SPH. However, the lack of high order consistency, tensile instability, undesired numerical fracture, and inaccurate gradient estimates in stress calculations [16] are common difficulties in SPH methods that can lead to unstable and inaccurate solutions in highly deformed materials and require a highly refined domain discretization and additional treatments/modifications to achieve desirable accuracy and stability and avoid the aforementioned issues. In explosive/impact welding simulations, such issues can produce unstable and unphysical particle distributions when large interfacial waves are formed [11]. Also, when the flyer plate is under severe bending, instabilities in displacement and stress can be shown throughout the plate [14,15]. Wave at the metal interface, a copper-steel [2]; b titanium-stainless steel [3]; c aluminum-aluminum [4] The reproducing kernel particle method (RKPM) was developed as a correction of SPH in recovering the consistency conditions in the conventional SPH method [17][18][19], and Chen et al. [20] extended it to nonlinear and large deformation problems. Thereafter, RKPM has been applied to contact-impact in metal forming [21,22], impact-induced plugging failure, spall fracture, and penetration [23][24][25][26][27]. Chen and Wu [28] proposed semi-Lagrangian (SL) RKPM that does not involve visibility criteria to effectively model problems with extremely large deformation and fragmentation, and SL-RKPM has been successfully applied to fragmentation problems [25][26][27][29][30][31].
Strong shocks are produced throughout processes of explosive welding both in explosives and solids, and it is one of the key physics to capture to avoid unphysical stress oscillation and overly dissipated solutions. Artificial viscosity [32,33] is widely used for its simplicity and ease of implementation. However, this technique requires user-tunable problem-dependent parameters to capture shocks with desirable accuracy. To embed shock physics in the RKPM formulation, Roth et al. [34] proposed Riemann-SCNI, that introduces a Godunov flux term into the volumetric energy under the strain-smoothing framework of stabilized conforming nodal integration (SCNI) developed by Chen et al. [35]. The Riemann-SCNI were verified and validated on high strain rate impact problems [36]. This approach, however, requires conforming cell updates, which is not practical for impact-fragmentation problems such as explosive welding (which entails jetting). In this work, a new particle-based Godunov-type algorithm for shock modeling that eliminates the need of conforming cells is proposed, and the weak form is consistently integrated with a stabilized nonconforming nodal integration (SNNI) [25].
For investigation of the weldability of a pair of metals, jetting has to be accurately modeled, which requires sufficient resolution in the domain discretization of the impacting plates through the plate-thickness direction. Also, the region close to the material interface experiences extreme material stretch and the resulting deformation is highly anisotropic. To maintain desirable locality and kernel stability in the zone of highly anisotropic plastic deformation, a deformation-dependent anisotropic kernel support update scheme is introduced.
While the SNNI technique can circumvent the need for conforming cells, relaxation of the conforming condition results in loss of accuracy and convergence. In this work, a strain-smoothing cell update scheme is proposed to attain high conformity and is shown to maintain high accuracy and optimal convergence under severe distortion. A variationally consistent integration correction that Chen et al. [37] developed is selectively applied in the areas of low conformity.
The remainder of this paper is organized as follows. The governing equations and basic formulations of the semi-Lagrangian (SL) reproducing kernel (RK) approximation are briefly described in Sect. 2. In Sect. 3, a particlebased shock algorithm is introduced, followed by Sect. 4 on domain integration methods. In Sect. 5, a deformation-dependent update of kernel supports and strain smoothing cells is proposed. In Sect. 6, the proposed method is demonstrated through a set of benchmark problems, then applied to simulations on impact/explosive welding, where the convergence of the jetting and the interfacial waves is examined, and the numerical models are validated in terms of plate velocities, impact angles, and interfacial waves, by using sets of experimental results. Concluding remarks are given in Sect. 7.

Conservation laws
Let = ( , t) be the mapping between the initial configuration of a body, Ω ∈ ℝ n d , and the configuration of the body at time t , Ω ∈ ℝ n d , with the space dimension n d , where and are the material and spatial coordinates of Ω and Ω , respectively.
The conservation equations for a continuous, non-heat conducting media described in Ω are where (0) and are the densities in the initial and current configurations, respectively, = − is the displacement vector, is the Cauchy stress tensor, E is the total energy per unit mass, and = ∕ is the deformation gradient. In Eqs. (2) and (3), d(⋅)∕dt denotes the material time derivative; = d ∕dt and = d 2 ∕dt 2 are the velocity and the acceleration, respectively, and ∇ is the gradient operators with respect to . The body force is neglected as the welding process finishes in a short time so that the body force does not play a significant role.
The total energy E is given as, where e is the internal energy per unit mass. (1)

Semi-Lagrangian reproducing kernel approximation
The reproducing kernel (RK) approximation h ( , t) of a function ( , t) defined in a domain discretized by a set of NP nodes, I | I ∈ Ω NP I=1 (see Fig. 3), is constructed as [17]: where Ψ I ( ) is the RK shape function of node I and I (t) is the nodal coefficient of node I in the approximation of h ( ) . Ψ I ( ) is obtained by correcting a kernel function Φ a − I with compact support measure a as follows: are basis functions. The coefficients b ( ) | |≤n are determined such that the following reproducing conditions are satisfied up to | | = n: By solving b ( ) | |≤n from Eq. (7), the shape functions Ψ I ( ) are obtained as follows: where −1 ( ) is the inverse of the moment matrix ( ) , and − I is the vector of basis functions: In the RK approximation, kernel function Φ a − I defines the order of smoothness and the locality, and the basis function vector − I controls the order of completeness of the approximation. These attributes allow the RK approximation more flexibility in constructing approximation functions with desirable completeness, smoothness, and compactness than the interpolation-type approximations used in FEMs.
In the Semi-Lagrangian (SL) RK formulation [28], suitable for classes of problems which entail fragmentation, such as jets in explosive welding, the RKPM nodes follow the material motion, while the distance measure, − I , in the kernel function Φ a − I is defined in the current configuration. Figure 4 shows the difference between Lagrangian and SL kernels. The Lagrangian kernel follows the material deformation (Fig. 4b), while the SL kernel does not (Fig. 4c). Therefore, an SL kernel support may cover a different set of material nodes throughout the simulation, while Lagrangian kernel support covers the same set of nodes. Note that, in Fig. 4c, a smaller number of nodes are covered by a SL support under a deformation if supports are not adjusted, and may result in the failure to satisfy the reproducing conditions, Eq. (7), under severe deformation. This point will be revised later in the text.

Particle-based Godunov-type shock algorithm
Explosive welding processes often involve strong shocks in the impacting metals. Effective shock algorithms are essential in capturing shock physics, as their strong tie with the colliding velocity determines characteristics of jetting formation in the explosive welding processes. To start, the weak form of the conservation of linear momentum in Eq. (2) is to find ∈ H 1 , = on Γ g , such that ∀ ∈ H 1 , = 0 on Γ g : where Γ g and Γ h denote the essential and natural boundaries, respectively, is the traction vector on Γ h , and ∇ ( ) ∶ ≡ u i,j ij where (⋅) ,j denotes a derivative with respect to x j in the current configuration. Throughout the text, repeated lowercase indices " i " and " j " imply summation unless otherwise stated. Let h and h be the RK approximation of and . Then, the corresponding (10) where the volumetric-deviatoric decomposition of the Cauchy stress has been employed: A standard Galerkin method developed directly based on Eq. (12) without effective shock algorithms could lead to incorrect modeling of shock speed and jump conditions, failure to meet entropy conditions, and severe non-physical oscillations in the presence of strong shock waves due to the Gibbs phenomenon at the shock front. To remedy this issue, artificial viscosity [32,33] is widely used for its simplicity and ease of implementation. However, this technique requires user-tunable problem-dependent parameters. Roth et al. [34,36] developed a shock algorithm which embeds the Godunov flux into Galerkin meshfree methods under the framework of SCNI. In this work, an extension of the Godunov-type shock algorithm called Riemann-SNNI is proposed such that the algorithm is purely particle-based, and thus does not require conforming integration cells for modeling fragmentation problems. This algorithm embeds the shock physics via the solution of a local Riemann problem for pair-wise interactions, and as a result, does not require tunable parameters. With such advantages, Riemann-SNNI is adopted in this study for both fluid and solid modeling.
In Riemann-SNNI, shock physics is embedded into the volumetric part of the internal force. By substituting Eq. (5) into Eq. (12), the deviatoric internal force f dev Ii and the volumetric internal force f vol Ii of node I are derived as In order to embed shock physics, f vol Ii is further modified by using the properties of the partition of unity of the RK shape function and the partition of nullity of the RK shape function derivatives in Eq. (7), following Hietel et al. [38]: where IJ−vol ij is the volumetric stress contribution from the interaction of node pair I-J , which is assumed constant in Eq. (15), and the coefficient vector IJ j is defined as As shown in Fig. 5, the coefficient vector IJ can be interpreted as an equivalent geometric coefficient by drawing a comparison to finite-volume methods and Riemann-SCNI, with the surface normal IJ , multiplied by the surface area of the virtual interface between node I and J . To investigate this geometric property of IJ in the proposed formulation, the following quantity is defined: where W L is the nodal weight of RK node L . With the above in hand, Eq. (15) can be expressed as where (15)  where h vol i is the volumetric stress contribution to the traction. Note that the last term of Eq. (20) is present because the stress divergence smoothing does not perform the integration over the natural boundary Γ h . Under this interpretation, Riemann-SNNI can be considered to be an extension of Riemann-SCNI as a particle-based approach. In addition, it is easily shown that IJ = − JI using the aforementioned interpretations, meaning that the exchange of stress between node I and J is equal and opposite, which maintains conservation of linear momentum in the shock algorithm.
Note that, in the SL-RK formulation, Ψ I ( )∕ x i can directly be obtained without requiring any mapping. Also, IJ j is computed by domain integration instead of integrating via the contours of a physical cell as done in the finite volume method. Such cell-free properties make Riemann-SNNI suitable for problems with large material distortion and fragmentation.
To introduce shock physics into vol I , a Godunov-type scheme [39] is introduced. In order to capture the shock discontinuity, an one-dimensional Riemann problem is locally defined between nodes I and J along the direction of IJ (see Fig. 5) with the left state P IJ L , v IJ L , IJ L and the right state The state variables are defined as (20) where P IJ * is the pressure solution of the local Riemann problem obtained from a Riemann solver. As such, the shock-enhanced volumetric internal force f vol * Ii is computed as follows: In order to solve the local Riemann problem, in this work, the Dukowicz Riemann solver [40], an approximate solver, is adopted for its robustness and efficiency: it requires approximately 20 scalar algebraic operations to solve a Riemann problem. Now, Eq. (12) can be expressed in the following semidiscrete equation: In this work, a central difference time integration method is utilized with a lumped mass matrix.
The conservation of mass, Eq. (1), is nodally solved by employing the forward Euler method with a time step size Δt as follows: where the superscripts n and n + 1 denote time steps.

Riemann-solution enriched shock algorithm for the energy equation
The Galerkin formulation of the energy equation, Eq. The approximation w h ( ) is constructed as follows:

Then, Eq. (26) can be written as
With the same procedure performed in Eqs. (15)- (22), the first term in the right-hand side of Eq. (28) can be approximated as follows: If the RK shape function Ψ I ( ) that is used for displacement is utilized, IJ j in Eq. (29) becomes the same as IJ j defined in Eq. (16).
is the velocity solution of the local Riemann problem defined in Sect. 3.1, To further simplify the equation, let E h ( ) be defined as where Ψ E I ( ) and E I are the energy shape function and the nodal coefficient, respectively. By applying the row-sum technique, the left-hand side of Eq. (28) is rewritten as follows: Now we have where M I = ∫ Ω Ψ I dΩ . Note that M I is the lumped mass matrix components used for the balance of linear momentum in explicit dynamics. For inviscid flow, (27) In this work, the Kronecker delta property is assigned to Ψ E I ( ) , i.e., Ψ E I J = IJ , and Eq. (31) becomes which leads to The forward Euler method is used to solve Eq. (36). This allows the energy to be nodally computed and directly used in the material laws without interpolation under the nodal integration framework. If the energy at an arbitrary location is of interest, the transformation methods [20,41] can be considered to recover the Kronecker delta property from standard RK shape functions.

Stabilized non-conforming nodal integration (SNNI) with particle-based shock algorithm
Chen et al. [35,42] proposed Stabilized conforming nodal integration (SCNI) as an alternative to Gauss quadrature to meet the linear integration constraint (the requirements to pass a linear patch test) in the domain integration of Galerkin meshfree methods, and to remedy the rank instability in direct nodal integration. Instead of using the direct nodal gradient, SCNI employs divergence with a smoothed gradient ∇ (⋅) in each nodal representative domain Ω L by where W L is the integration weight associated with node L . The strain smoothing cell Ω L is conforming such that Ω = ⋃ NP L=1 Ω L and Ω I ∩ Ω J = ∅ with I ≠ J (see Fig. 6a). Chen et al. [35] showed that, by introducing the gradient smoothing, the linear integration constraint is satisfied and optimal convergence is achieved for RKPM with linear bases. Because SCNI requires conforming cells, it is suitable for Lagrangian formulations where the smoothed gradients are computed only once in the initial configuration. However, for the semi-Lagrangian approximation which is constructed on the current configuration, the regeneration of conforming cells at each time step is impractical. In order to address this issue, Guan et al. [25] proposed stabilized non-conforming nodal integration (SNNI). SNNI applies the same gradient smoothing concept as introduced in SCNI but the smoothing cells can be non-conforming as shown in Fig. 6b. The shape of the smoothing cell can be chosen as a cube or sphere for simplicity and ease of implementation. The disadvantage of relaxing the conforming construction is suboptimal convergence rates, and lower accuracy is obtained compared to SCNI. Nevertheless, SNNI is considered for the semi-Lagrangian formulation due to the balance of stability, efficiency, and accuracy, as the need for constructing conforming cells is avoided. In Sect. 5.2, an algorithm for the update of strain smoothing cells for better accuracy will be introduced, which allows SNNI to converge at the optimal rates for linear bases. The accuracy of SNNI can be further enhanced by a variationally consistent correction discussed in the next section.
With SNNI, the internal force terms are calculated as follows:

Variationally consistent correction of SNNI
During a simulation of explosive welding, highly non-conforming strain smoothing cells can be generated at the metal interface caused by the extreme plastic deformation, and this leads to low solution accuracy. Chen et al. [37] proposed a variationally consistent (VC) integration to increase solution accuracy and convergence by satisfying the integration constraints, and this method can be considered to address the low accuracy issue due to the highly non-conforming strain smoothing cells. While the VC correction is applicable to arbitrary numerical integration methods with an arbitrary order of the RK basis to achieve desirable convergence rates, in this study, the correction is made to SNNI with the linear RK basis, such that the assumed gradient of the test shape function, ∇ Ψ I ( ) satisfies the first order integration constraint: where c Ii is the nodal coefficient, the operator ∫ represents numerical integration, and R I ( ) is defined as where supp Ψ I ( ) denotes the kernel support of the RK shape function Ψ I ( ) . The unknown coefficient Then, the internal force terms in Eqs. (38) and (39) is corrected as follows: Due to the property in Eq. (41), For a kinematically admissible function space, we have The fact that Eqs. (47) and (52) are identical indicates that the RKPM Riemann-SNNI formulation yields linear exactness with the VC correction introduced in Eqs. (44) and (45).

Deformation-dependent updates
of kernel support and smoothing zone

Update of kernel supports
During the numerical simulation of the explosive welding processes, the distribution of RKPM nodes changes extensively due to large plastic deformation. Because the nodal points in SL-RKPM follow the material points' motion, kernel supports with a small normalized support size can result in numerical fracture, or unphysical distribution of nodes. On the other hand, the employment of very large support size (large locality) is inaccurate in modeling highly local phenomena such as strain localization, shearbands, and fracture. Liu et al. [12] applied an isotropic kernel support update formula to avoid unphysical distribution of particles in the SPH computation of explosive welding. However, for the explosive welding process, extremely large shear deformation produced along the material interface leads to excessive material stretch in the direction nearly parallel to the interface, where the deformation is highly anisotropic. Also, for investigation on the weldability of a pair of metals, the metal jet must be accurately modeled, which requires sufficient resolution in the domain discretization of the impacting plates through the plate-thickness direction. As such, anisotropic domain discretization should be introduced for efficient modeling of the jet, and correspondingly, an anisotropic kernel support update scheme is needed to maintain the locality without numerical fracture and non-physical distribution of nodes.
In this work where plane strain is assumed, the kernel supports are updated by employing the following equations: where as shown in Fig. 7, a x ′ and a y ′ are the dimensions of a kernel support in the x ′ -and y ′ -directions, respectively, and is the orientation angle of the support. v x ′ and v y ′ are the velocities in the x ′ -and y ′ -directions, respectively. ( ) represents the time derivative, and (⋅) ,x � and (⋅) ,y � denote the quantity is the spatial derivative of (⋅) with respect to x ′ and y ′ , respectively. In this work, the local x ′ -y ′ coordinate is defined such that x ′ -direction is initially parallel to the material interface. The support dimensions a x ′ and a y ′ of a node are no longer updated once the material is fully damaged at the node (see Sect. 6.1.1 for the damage model used in the work). As such, the SL-RK approximation allows the material to be naturally separated.

Update of SNNI strain smoothing cells
The strain smoothing cells employed for SNNI can be proportionally adjusted by using the updated kernel supports without additional computational cost as: where D x ′ and D y ′ are the updated dimensions of a smoothing cell in x ′ -and y ′ -directions, respectively, and is the orientation angle of the smoothing cell. The superscript 0 on parameters denotes their values at the initial configuration. Figure 8 schematically illustrates the effect of using the updated smoothing cells under rigid body motion, bending, and shear deformation. Although the cells are not strictly conforming, they maintain a high level of conformity for the deformation modes produced during explosive welding simulations. Herein, the updated strain smoothing cell algorithm is termed adaptive rectangular cell. Similar to the kernel support update, the smoothing cell of a node is no longer updated once the material is fully damaged at the node. To further investigate the effect of the adaptive rectangular cell on the solution accuracy, convergence studies are carried out with the following two-dimensional Poisson equation: where Ω = (−1, 1) × (−1, 1) . The details of the problem considered in this work are listed in Table 1.
Linear bases with a normalized support size of 1.8 is employed for the RK approximation. The domain is discretized with 120, 435, 1653, and 6441 nodes. To construct discretization, uniformly distributed reference nodal points * I = x * I1 , x * I2 are considered. Then, the nodal points, I = (x I1 , x I2 ) , are defined in the following manner to resemble the wavy interfacial morphology of a specimen fabricated by explosive welding:   Figure 9 shows four different levels of refinement used for the convergence study. As shown in Fig. 10, five different types of strain smoothing cells in SNNI are considered for domain integration with different levels of conformity of the smoothing cells. A conformity parameter is defined as follows: where NP is the number of nodal points and n d is the space dimension, e.g. n d = 2 for a 2-dimensional problem, and e Ii is a normalized deviation from the first order integration constraint in Eq. (41): where a I1 and a I2 are the support sizes of node I in the x 1 -and x 2 -directions, respectively, and R I is defined in Eq. (42). In Eq. (58), repeated indices do not imply summation. The conformity of each type of cell considered in this study is specified in Fig. 10. As shown in Fig. 11, the convergence performance of the RKPM integrated with SNNI methods with adaptive rectangular smoothing cell of Fig. 10b and perturbed adaptive rectangular smoothing cell I ( = 3.13 ) of Fig. 10c is comparable to the solution of SNNI with the conforming Voronoi smoothing cells (equivalent to SCNI). In Fig. 11, the numbers in the legends denote the average convergence rates. The accuracy and convergence rates tend to decrease as the conformity decreases.
Although RKPM-SNNI with the perturbed adaptive rectangular smoothing cells (Fig. 10c, d) show superior convergence performance compared to that with the nonupdated non-conforming smoothing cells (Fig. 10e), the convergence rates can become suboptimal when smoothing cell with low conformity is employed. This implies that the numerical integration may need to be corrected to obtain desirable accuracy in the region where the SNNI Due to the additional computational cost involved in Eq. (59), in this study, effective plastic strain p is alternatively used as a criterion to determine the nodes for which VC-SNNI is applied. Figure 12 shows a distribution of rectangular smoothing cells adaptively updated during the simulation of a magnetic pulse welding of Titanium and Copper (see Sect. 6.5 for details.). The nodes with I ≤ 3.0 are coded blue in Fig. 12b. It is shown in Fig. 12c that the zone with p ≥ 0.7 includes most of the smoothing cells with low conformity. As a conservative approach, p ≥ 0.5 is chosen in the explosive welding simulation presented in Sect. 6. As a remark, I could also be calculated only at certain time steps, say, every 100 time steps, to determine the nodes selected for the adaptive VC correction.

Numerical examples
In this section, the accuracy of particle-based shock algorithm is first examined by solving high velocity impact problems and a one-dimensional detonation problem. Then, numerical examples of impact/explosive welding processes are presented to validate the proposed numerical methods. For all the examples of explosive/impact welding, plane strain is assumed as the out-of-plane dimension of the plates are sufficiently larger than the in-plane dimensions. For all the examples, RKPM equipped with the cubic B-spline kernel function, linear monomial bases, and SNNI is selected. Normalized support sizes of 1.25 and 2.0 are used for the metal plates and explosives at initial stage, respectively. The algorithms for the anisotropic update of the kernel supports and strain smoothing cells are applied to the metal plates. The upper limit of the sizes of the kernel supports and strain smoothing cells is set to three times their initial sizes. Because the hydrostatic expansion is the major deformation mode of the detonation gas, an isotropic kernel support update, a i = 0 a 0 i , is utilized for the explosives where a i and are the support size in the i-th direction and the mass density, respectively, with superscript 0 denoting the quantity is at t = 0 . Similarly, the support size that is three times the initial support size is set as the upper limit for the explosives. For effective modeling of material fragmentation, the Quasi-linear approximation [29] is employed (See "Appendix 1"). To impose essential boundary conditions, the boundary singular kernel method [41] is used.

Material models
The material models employed in this work is summarized. The specific model parameters are provided in the individual numerical examples.

Constitutive equation and failure criterion
The Johnson-Cook constitutive model [43] and failure criterion [44] are used to describe the strain-rate dependent behavior of the metals.
In the Johnson-Cook models, the accumulation of damage D is influenced by plastic strain, strain rate, stress state, and temperature, and expressed as where is the equivalent plastic strain and f is the failure strain, where T * = T − T r ∕ T m − T r is the normalized temperature at temperature T used to characterize thermal softening with the reference temperature T r and the melting temperature T m , and D 1 -D 5 are the material constants and the stress triaxiality ratio * is given by Following [45], the damage is not updated when * ≤ −1∕3 . The triaxiality ratio-based law allows to capture realistic damage evolution. The damage rapidly grows when the material ejects as a jet due to the positive triaxiality ratio. Once the material is fully damaged, the reevaluation of SL-RK shape functions in the current configuration allows the damaged material to be naturally separated from the remaining body. In contrast, when the material is subjected to high pressure, the triaxiality ratio is negative and the damage hardly develops, which The equivalent stress eq , is defined as where ̇ * =̇∕̇0 represents the plastic strain rate normalized by a reference strain rate. The adiabatic condition is assumed in this work to calculate the temperature T . The Johnson-Cook constitutive equation is solved using the radial return mapping algorithm in this work.

Equation of state
For modeling the detonation product, the following Jones-Wilkins-Lee (JWL) equation of state [46] is employed: where P , , 0 , and e are pressure, density, initial density, and internal energy per unit mass, respectively; A , B , R 1 , R 2 , are the model constants.
For modeling the volumetric behavior of metals, the Mie-Gruneisen equation of state [47] is used to capture the sharp increase of pressure when the density and the internal energy increase. In the Mie-Gruneisen equation of state, the pressure P is described as follows: where , 0 , and e are the density, the initial density, and the internal energy per unit volume; = ∕ 0 − 1 ; and 0 , C , and S are model constants. For the shock algorithm introduced in Sect. 3.2, the equations of state are employed in solving the local Riemann problems.

Impact of elasto-plastic bars
This problem is analyzed to examine the effectiveness of the proposed node-based Riemann-SNNI shock algorithm in capturing elasto-plastic shock waves. Consider a one-dimen-  Figures 13 and 14 show the solutions at t = 0.2 obtained by Riemann-SNNI and uncorrected SL-RKPM, respectively. Compared to the analytical solution, Riemann-SNNI accurately captures the stress and velocity profiles including the elastic precursors followed by the plastic shock waves without oscillation in contrast to the uncorrected counterpart where oscillation is present.

Two-dimensional high-velocity plate impact with rarefactions
A symmetric impact of 8 mm × 2 mm aluminum plates is modeled as shown in Fig. 15a. The impact velocity of 1000 m/s is imposed as the initial condition of the upper  Table 2.
As shown in Fig. 15, rarefaction waves are produced at the impact plane and propagate toward the opposite boundaries. The maximum pressure obtained using Riemann-SNNI is 7.95 GPa which agrees well with the compressive pressure of 8.0 GPa experimentally observed by Marsh [48]. Compared to the numerical results of uncorrected RKPM shown in Fig. 16 where severe oscillation is produced with maximum pressure of 15.67 GPa, the proposed particle-based shock algorithm accurately captures the sharp pressure profile without oscillation.

Modeling of explosive detonation
In simulations of explosive welding, modeling the detonation process is required. In this work, high explosives (HE) are chosen as the focus. As shown in Fig. 17a, detonation of HE typically generates three distinct zones: unreacted high explosive, reaction zone, and detonation product. If the detonation process is assumed to be steady, the Chapman-Jouguet (C-J) hypothesis can be utilized. Under the C-J hypothesis, a plane shock front propagates with a constant detonation velocity D , and the fully reacted detonation product right behind the C-J plane is characterized with C-J pressure P CJ and C-J specific volume V CJ = 1∕ CJ where CJ is C-J density. P CJ and V CJ are determined such  that P CJ , V CJ becomes a point of tangency at which the Rayleigh line meets the Hugoniot curve of the detonation product, as shown in Fig. 17b. The mass density, CJ , and the particle velocity, v CJ , at the C-J point are expressed as CJ = 0 ( + 1)∕ with = 0 D 2 ∕P CJ − 1 and v CJ = 1 − 0 ∕ CJ D , respectively, with initial mass density of explosive, 0 .
For numerical modeling of detonation process, the total energy per unit mass at C-J point, E CJ = e CJ + 0.5v 2 CJ , is strongly imposed to a node when the node is placed within a thin initiation zone which propagates with a constant detonation velocity D . The internal energy per unit mass at C-J point, e CJ , is determined such that the JWL equation of state, Eq. (64), is satisfied given P CJ and CJ . In this study, the width of the initiation zone is set to the nodal spacing, h.
As a demonstrational example, a one-dimensional detonation of 0.2-m-long TNT initiated at the center of the explosive is simulated. Due to the symmetry of the problem, Ω = {x|0 ≤ x ≤ 0.1} with an essential boundary condition u(0) = 0 is considered. The explosive is modeled with 2001 equally spaced RK nodes. The JWL equation of state is used and the model constants are listed in Table 3. The pressure distribution, particle velocity distribution, and density ratiopressure profile obtained by the RKPM simulation are shown in Fig. 18a-c, respectively. The peak pressures and peak particle velocities match the corresponding C-J values. The density ratio-pressure profile agrees well with the reference solution under adiabatic expansion provided by Lee et al. [46].

Geometrical change of HE detonation-driven flyer plate
In this section, RKPM-predicted geometrical change of Aluminum plates driven by HE is compared to experimental data: plate velocity and impact angle, two key parameters that are used to estimate weldability in a practical  application. For the experiments presented in Sect. 6.4, the explosion was initiated by the line wave generator that creates a planar detonation front, hence the plane strain assumption is valid.  Figure 20 shows the vertical plate velocity versus vertical plate displacement profile. In the legend, "RKPM" and "Test" denote the numerical results and the experimental measurements, respectively. "C1", "C2", "C3", and "C4" represent values computed or measured at x C1 , x C2 , x C3 , and x C4 , respectively. The velocity-displacement curve denoted by "RKPM: Unconfined" is a case computed in the absence Fig. 18 One-dimensional detonation problem: a pressure distribution, b particle velocity distribution at 2 μs to 14 μs with an interval of 2 μs, and c density ratio-pressure profile   The next example is the comparison of numerically predicted and experimentally measured impact angles. Four aluminum to aluminum welding cases listed in Table 6 are considered. The 76.2 mm long flyer plate and base plate are modeled with 16,850 and 13,817 RK nodes, respectively. The explosive is modeled with 3592 nodes for Case I and II, and with 7184 nodes for Case III and IV. The numerical results are compared with the experimental data in Fig. 21.
The experimental values were measured with X-ray flash radiography techniques. The numerical values were measured after the impact angle reached a constant value. The impact angles measured from the numerical simulation show a good agreement with the experimental observation with an average error of 5%.

High velocity impact welding
The validation problem considered in this section is a high velocity impact problem using magnetic pulse without explosion which is an experiment conducted by Vivek et al. [51]. This is a well-controlled experiment that allows systematic validation of the predicted interfacial morphology as well as jet formation. The experiment was magnetic pulse welding of commercially pure Titanium grade 2 and Copper, where a magnetic pulse acts on the bottom surface of the flyer plate to accelerate the plate upward as shown in Fig. 22. The impact angle between the plates is 24°. As listed in Table 7, the thicknesses of the flyer plate and the base plate used in the simulation are 0.508 mm and 0.762 mm, respectively. The impact velocity of the titanium flyer plate driven by the magnetic pulse is set as an initial velocity v 0 p , of 770 m/s, which produces a weld velocity of approximately 1900 m/s. An essential boundary condition, u y = 0 , is imposed on the top surface Ω g of the upper base plate. The material models parameters are listed in Tables 8 and 9.
The plates are discretized with the initial nodal spacings h x ′ and h y ′ (Table 10) along the normal and tangential directions, respectively, with respect to the contact surfaces, i.e., the top surface of the flyer plate and the bottom surface of the base plate. The nodal spacing gradually increases as they move away from the contact surfaces with maximum spacing of h max in both directions given in Table 10. In this work, three different levels of refinements are utilized to investigate the solution convergence. The details of the discretizations are listed in Table 10. Smaller nodal spacing h y ′ is chosen compared to h x ′ to better capture the jet formation.
All models ran on a single compute node with two CPUs (16 cores per CPU, 2.3 GHz per core) and 256 GB of RAM. The run-times of Model 1, Model 2, and Model 3 were 1.18 × 10 4 s, 3.29 × 10 4 s, and 1.09 × 10 5 s, respectively. Also, the run-times per time step of each model were 0.406 s, 0.692 s, and 1.753 s, respectively. Note that the computational efficiency can be enhanced by dynamic load balancing and the employment of coupling of Lagrangian and semi-Lagrangian RKPM [53]. Figure 23 shows the progressive deformations during the welding. The collision point moves from left to right with the weld velocity of V w = 1900 m/s. A hump is generated  at the collision point by the stagnation of material, which initiates the jet formation. The jetting starts approximately at t = 0.35 μs as shown in Fig. 23a. The amount of jetting gradually increases until it reaches a steady state. The interfacial wave is also gradually built up and reaches a steady state at approximately half of the total plate length. Figure 24 shows the jet formation predicted by three different discretizations. The amounts of the jets from the different models converge to a similar range. As shown in Figs. 23 and 24, the jet is mainly ejected from the titanium flyer plate. It is a reasonable prediction as the density of copper is twice the density of titanium. Figure 25 shows the steady-state interfacial wave obtained numerically and experimentally. Model 2 and Model 3 results show convergent wave pattern, wavelength, amplitude of wave, and It is worth noting that, as shown in Figs. 23, 24 and 25, the very fine welding features such as the jet layers and the interface waviness and tails are well captured. These excessive anisotropic deformations are properly modeled with RK shape function with only 1.25 normalized support size at the initial stage. This feature is hard to achieve throughout the entire welding processes without a proper update of the kernel supports and smoothing cells in semi-Lagrangian meshfree methods presented in Sect. 5.

Explosive welding
In this section, an explosive welding experiment conducted by Bahrani and Crossland [2] is used as further validation of the welding simulation in this study. As shown in Fig. 26 and Table 11, a 3.2-mm-thick PETN-based sheet explosive is used for accelerating a 0.8-mm-thick stainless steel flyer plate which is placed above a 1.2-mm-thick mild steel with 6° of an oblique angle. The rubber buffer used in the original experiment to protect the flyer plate from explosive-induced damage is ignored in this work. An essential boundary condition = 0 is imposed on Ω g . The explosive, flyer plate, and base plate are discretized with 21,973, 67,529, and 75,573 RK nodes, respectively. The initial nodal spacing at the metal surfaces to be in the contact region is 4.1 μm in the plate-thickness direction and 8.2 μm in the longitudinal direction. The nodal spacing gradually becomes larger with increasing distance, as described in Sect. 6.5. The material model parameters are listed in Tables 4, 12, and 13. Figure 27 shows the predicted pressure field and equivalent plastic strain field produced during the welding process. As shown in Fig. 27a-d, the propagation of a sharp shock through the flyer plate generated by the detonation and a radially propagating shock through both plates by impact is well captured by the shock algorithm introduced in this work. It is shown in Fig. 27b, c that the impact induced bending in the flyer plate yields an impact angle larger than the initial plate angle , which favors jet formation. Compared to the sound speed of steel, approximately 6000 m/s, the numerical weld velocity is approximately 4000 m/s at steady state, which causes the radial shock propagation. Note that, Cowan and Holtzman [55] established, when the elastic strength of material is exceeded, a subsonic weld velocity generates a jet which is necessary for welding. As can be seen in Fig. 27e, the maximum equivalent plastic strain of around 1.5 shows the excessive plastic deformation along the weld interface where the elastic strengths of the Fig. 22 Simulation setup of the plate dimension, the initial angle, and the initial plate velocity. In the experiment conducted by Vivek et al. [51], a magnetic pulse propels the flyer plate  1 3 materials are well-exceeded. The resulting jet is captured by the RKPM simulation as shown in Fig. 28a. Figure 28a also shows the gradual build-up of the interfacial wave, including the smooth-wavy transition, evolution of interfacial wave, and a steady state region, which all qualitatively agree with the experimental observations. The interfacial wave reaches the steady state slower than the magnetic pulse welding presented in Sect. 6.5. This is attributed to the gradually increased impact angle and vertical plate velocity v p plotted in Fig. 28b. The edge effect of the detonation and the small initial standoff distance between the flyer and the base plates lead to smaller and v p at the beginning of welding as shown in Fig. 27a-c. With the larger initial standoff distance between the flyer and base plates with larger horizontal distance from the left end of the plates, as shown in Fig. 26, and with the decreasing edge effect in the plate axial direction, and v p increase gradually and approach a steady state where ≈ 13 • and v p ≈ 1000 m/s. The interfacial wave pattern also reaches a steady state. Figure 28c shows a snapshot of the fully developed wave in the numerical simulation. Compared to the experimental result (Fig. 28d), the overall magnitude and the shape of the wave are well-captured. Figure 28e presents the equivalent plastic strain field. Within the narrow region along the material interface, highly localized plastic strain is accumulated, demonstrating the likelihood of grain refinement induced by adiabatic strain localization. The temperature distribution obtained by the proposed method is shown in Fig. 28f. The upper limit of the color range is the melting temperature Jet formation at approximately t = 0.9 μs with different discretizations: the initial nodal distance at the interface, h x , h y , is a (12, 6), b (8,4), and c (6,3) in micrometer of the mild steel, 1811 K. The material points with white color indicate the potential molten area. The tail and head regions contain a large amount of molten material points, which is consistent with the experimental observation ( Fig. 28d) with molten zones formed near the tails and heads. Also, the numerical results demonstrate potential melted area along the material interface. This implies that the processing parameters used in the example could result in a considerable amount of heat affected zones that can negatively affect the advantage of solid-state welding techniques.

Conclusion
In this study, a semi-Lagrangian RKPM framework with a particle-based shock algorithm and deformation-dependent kernel support and gradient smoothing cell update algorithm for simulation of impact/explosive welding was proposed. A particle-based Godunov-type algorithm was introduced such that the shock physics in both fluids and solids were embedded into the SL-RKPM formulation. The formulation is shown to ensure conservation of linear momentum. Also, the proposed approach does not require conforming integration cell construction and use of user-tunable parameters, making it easy to implement for impact-fragmentation problems. A high velocity plate impact problem and a detonation problem demonstrated that the proposed shock algorithm captures sharp shock propagation without spurious oscillation. An algorithm for the update of kernel supports was also proposed to maintain locality, preclude numerical fracture, and unphysical distributions of nodes under the extreme deformation induced by impact/explosive welding processes. Based on the updated supports, the gradient smoothing cells were accordingly updated to obtain high conformity and achieve optimal solution convergence in RKPM with SNNI for domain integration. For accurate domain integration near the metal interface where high conformity is difficult to achieve, variationally consistent SNNI was adaptively employed. The effect of gradient smoothing zone conformity on the solution accuracy and convergence properties was studied in terms of conformity measure, and the regions where VC-SNNI was applied was selected accordingly. Further, the proposed SL-RK approximation was employed with the quasi-linear technique to efficiently model the jet formation, while keeping desirable solution accuracy.
The particle-based Godunov-type shock algorithm captured both the HE explosion-driven sharp shock propagation and the plate impact-driven radial shock propagation in metal pieces without oscillation in simulation of explosive welding. Stable nodal distributions at the region close to the metal interfaces, where the deformation was large and  highly anisotropic, were attained by the proposed kernel support and strain smoothing cell update algorithm. The explosive welding kinematics including the plate velocities and impact angles of flyer plates driven by a high explosive detonation were accurately predicted by the proposed methods when compared to experimental results. In both the magnetic pulse welding and the explosive welding simulations, metal jets were well-captured under a subsonic weld velocity. The interfacial waves in their steady state quantitatively and qualitatively agreed with the experimental results reported in literature in terms of wavelength, wave amplitude, and direction of tails, along with the overall buildup of the waves.
Compared to the standard SPH approaches, the proposed methods showed promising results where more stable and realistic nodal distributions along the metal interface with large interfacial waves was obtained, and the sharp shock propagation was accurately captured as shown in Figs. 13, 15, and 27, without any pre-defined tunable parameters. Further, the gradient smoothing cell update algorithm combined with the adaptive VC correction ensures better solution convergence with refinement. Although the calculation of the semi-Lagrangian RK shape functions at every time consumes additional CPU, this disadvantage can be improved by coupled Lagrangian and semi-Lagrangian RKPM techniques recently proposed in [53] and will be included in the future work.

Appendix 1. Quasi-Linear RK Approximation
When a problem entails material fragmentation, such as in an explosive welding simulation, the moment matrix ( ) in Eq. (9) can easily become singular due to the lack of enough support coverage in the SL approximation. Yreux and Chen [29] proposed the quasi-linear (QL) RK approximation as a remedy of the possible singular moment matrix, and this technique is employed in this study. The QL-RK approximation ensures a nonsingular moment matrix by introducing an additional nonsingular matrix * ( ) into the standard moment matrix as follows: Here, * ( ) is defined as, with a set of N S sampling points * k ( )  where S( ) = I |Φ a − I ≠ 0 . A set of six ( N S = 6) sampling points is employed here, defined as where 1 = kh x , 0, 0 , 2 = 0, kh y , 0 , 3 = 0, 0, kh z where h x , h y , and h z are nodal spacing in x -, y -, and z-direction, respectively. In this work, k is set to 0.1 and is set to 0.001 as suggested in Yreux and Chen [29]. The QL-RK shape function Ψ I ( ) is constructed as follows: