Generalized coordinate smoothed particle hydrodynamics with an overset method in total Lagrangian formulation

This study proposes a generalized coordinate smoothed particle hydrodynamics (GCSPH) method coupled with an overset method using the total Lagrangian formulation for solving large deformation and crack propagation problems. In GCSPH, the physical space is decomposed into multiple domains, each of which is mapped to a generalized space to avoid coordinate singularities as well as to flexibly change the spatial resolution. The SPH particles are then non‐uniformly distributed in the physical space (e.g., typically in a boundary‐conforming manner) and defined uniformly in each generalized space, similar to the standard SPH. The SPH particles in the generalized and physical spaces are numerically related by coordinate transformation matrices. The use of non‐uniform particle distributions decreases the total number of particles, thus significantly reducing the simulation cost. Three numerical cases: three‐dimensional brittle crack propagation, Taylor impact, and plugging failure, are presented to validate the proposed method. The numerical results are sufficiently accurate compared with the standard total Lagrangian SPH. Furthermore, these results also show the benefits of GCSPH, which not only effectively reduces the computational cost but also eliminates the effects of particle arrangement. These advantages allow GCSPH to perform high‐resolution three‐dimensional problems, which are otherwise costly to perform with the standard SPH.

methods, such as the finite-element method (FEM), are effective and have been widely employed in these problems. However, mesh-based methods often suffer from a significant drawback in mesh generation for complex and largely deformed boundaries owing to mesh distortion and entanglement. Hence, several meshfree methods have been proposed to overcome this drawback. 1 Smoothed particle hydrodynamics (SPH) is a well-established meshfree method originally developed by Gingold and Monagan 2 and Lucy 3 for astrophysics applications. SPH was subsequently extended to free-surface-flow problems, 4 large strain solid mechanics, including crack propagation and impact problems, [5][6][7][8] and fluid-structure interaction problems. [9][10][11][12][13] Especially in the case of fluid/structure interaction problems, the Lagrangian nature of SPH has the benefit of easily defining the solid/fluid interface even when the structure is highly deformable. For these reasons, various new algorithms and techniques for high-accuracy SPH simulations have been developed, 14,15 and further applications of SPH in various fields are expected in the future. Other meshfree methods, such as the element-free Galerkin method, 16,17 meshless local Petrov-Galerkin method, 18 moving particle semi-implicit method, 19,20 and material point method, 21 have also been developed in recent decades; in particular, the reproducing kernel particle method adds a correction function in the kernel representation to improve its accuracy and efficiency, especially in the case of impact and large deformation problems. 22 SPH can is the basis of many other meshless methods and is being continuously developed in many practical engineering problems. Specifically, crack initiation, propagation, and failure are basic and important issues in solid mechanics, and numerous analyses have been performed using SPH. [23][24][25][26][27][28] Nevertheless, the original SPH has several shortcomings in its formulation, for example, inconsistency, 29,30 tensile instability, 31 and rank deficiency. 32 In particular, tensile instability is the most critical shortcoming in solid mechanics problems. The source of this problem has been identified as the SPH formulation in the deformed current configuration, that is, the updated Lagrangian formulation. 31 Based on various studies on tensile instability, 29,[33][34][35] total Lagrangian SPH (TLSPH) 30,36 was proposed so that the governing equation, including the kernel function, could be formulated in the reference or initial configuration. Rausch et al. 37 modeled soft tissue damage and failure using TLSPH. Islam et al. 26,38 demonstrated deflection of a beam under impact, the crack propagation in a notched beam, Taylor impact test, and the Kalthoff-Winkler experiment using TLSPH. Young et al. 39 proposed the coupling scheme between total Lagrangian and Eulerian formulated SPH for high velocity impact on an aluminum plate. Following these studies, this study employs a total Lagrangian formulated SPH to solve impact and crack propagation problems.
In standard applications of SPH to solid mechanics, the particles are initially aligned uniformly in an analysis domain so that the interaction between any two particles can be effectively controlled by the kernel function that has an isotropic shape in the physical space. Owing to the nature of this isotropic kernel function, particles must be initially placed uniformly for accurate computation of the kernel approximation, and the locally controlled resolution for reducing the computational cost often used in FEM is not allowed in conventional SPH ( Figure 1). However, it is often beneficial to locally control the resolution, especially when considering complex geometries and crack propagation problems. Such non-uniform resolution methods have been extensively developed for fluid dynamics problems, wherein an adaptive particle refinement technique, that is, particle splitting and merging based on refinement criteria, are extensively used. For instance, Feldman and Bonet 40 proposed a particle-splitting method for adaptive refinement in dynamic fluid dynamics problems and investigated the numerical errors. Furthermore, Vacondio et al. 41 introduced a variable-resolution technique by splitting a parent particle into child particles and merging them for refinement and coarsening, respectively. Subsequently, several adaptive particle methods have been proposed for controlling the resolution in flow problems. [42][43][44][45][46] However, these techniques often encounter a problem in satisfying the conservation of mass, 42 and generally involve cumbersome implementation, especially for parallel computation. In solid mechanics, the particle distribution does not change significantly with time, unlike fluid dynamics. Therefore, static refinement 47 can be primarily considered, wherein the refinement zone with particle splitting is static and defined beforehand depending on the physical properties of the computational targets. [40][41][42]44 Such a static refinement with particle splitting is effective if the physical properties that require a high resolution are predetermined. Nevertheless, the particle splitting technique requires some empirical criteria for changing the resolution and often complicated implementation.
Yashiro and Okabe 48 proposed SPH in generalized coordinate systems, that is, generalized coordinate SPH (GCSPH). In the original GCSPH, particles are nonuniformly distributed in the physical space (e.g., typically in a boundary-conforming manner), whereas the governing equation is formulated in the generalized space with a uniform particle distribution. Notably, in GCSPH, the governing equations in the generalized space are derived by conducting generalized tensor analysis, which employs the covariant derivatives and makes the coordinate transformation F I G U R E 1 Schematic of locally controlled resolution in FEM and SPH applicable to high-order derivatives frequently required in solid mechanics. The capability of the abovementioned method 48 has been well demonstrated in the case of a quasi-static three-point bending problem of a thin plate and a high-velocity impact problem, where the intervals between the particles are varied in both the thickness and the in-plane directions, leading to a more efficient simulation compared to the standard SPH.
Although the governing equations have been successfully formulated in the generalized space, their demonstrations have been limited to non-curved geometries; thus, the application of GCSPH to curved geometries has not been attempted for solid mechanics problems thus far. Note that for fluid dynamics problems, Tavakkol et al. 49 proposed a similar method as curvilinear SPH with coordinate transformation using the chain-rule relation between the Cartesian and curvilinear coordinates. Meanwhile, the use of such coordinate transformations often involves a coordinate singularity: for instance, the center axis of the cylindrical coordinate cannot be uniquely defined in the generalized space. The cylindrical coordinate is often useful in solid mechanics problems, such as the Taylor impact problem, as will be demonstrated later in this article.
Thus, this study aims to extend the original GCSPH 48 to more complex geometries, including curved boundaries, for which the governing equations with covariant derivatives are reformulated on the basis of total Lagrangian SPH using the components of a Cartesian coordinate. The present GCSPH is inspired by the isoparametric elements developed in the FEM. Hence, it is possible to use the well-established knowledge on FEM meshing strategies. Furthermore, to overcome the aforementioned coordinate singularity and provide a more flexible refinement technique, we propose an overset method for GCSPH.
The remainder of this article is organized as follows. Section 2 describes the methodology of GCSPH in the total Lagrangian formulation with an overset method. Section 3 presents representative numerical examples (crack propagation and impact problems) to validate GCSPH and demonstrate its capability in reducing the computational cost. One of the examples simulates brittle crack propagation of a notched beam, while the other two target the Taylor impact test and plugging failure of a Weldox 460 E steel plate. Finally, Section 4 concludes the article.

GENERALIZED COORDINATE SMOOTHED PARTICLE HYDRODYNAMICS
This section describes the total Lagrangian formulation of GCSPH coupled with an overset method. First, the coordinate transformation between the physical and generalized spaces is introduced based on tensor analysis. 50,51 Second, the SPH approximation in the generalized space is described with its kernel representation and discretization using a particle approximation. Subsequently, the governing equations in the physical space can be reconstructed using the generalized coordinate-based SPH approximation and the coordinate transformation matrix. In addition, an overset method is introduced to enhance the present GCSPH. Finally, the numerical algorithm for the proposed method is presented with the governing equations in the total Lagrangian formulation.
Hereafter, the Einstein summation convention is used for these indices. The covariant and contravariant bases g i and g i are defined as where { , , , … } and {i, j, k, … } are varied as {1, 2, 3} and will be used for the physical and generalized coordinate systems, respectively. Note that the physical space is represented by the Cartesian coordinate system; thus, the covariant and contravariant components are undifferentiated, so that only the subindices { , , , … } are used. The Cartesian basis can be written as We are concerned with the coordinate transformation for the derivatives of an arbitrary physical quantity (scalar, vector, and tensor) between the Cartesian and generalized coordinate systems while keeping the components in the Cartesian coordinate system. According to the notations presented above, an arbitrary vector is written in each coordinate system as follows: Let us consider the coordinate transformation for the covariant gradient of an arbitrary vector as where ∇ = g i ( ∕ i ) and Equation (1) has been used. Similarly, the coordinate transformation of a divergence operation is given by Furthermore, the gradient vector of a scalar is transformed as Similarly, an arbitrary second-order tensor is written in each coordinate system as Accordingly, the coordinate transformations of ∇ ⋅ and ⋅ ∇ are as follows: All the derivatives of the vectors and tensors in Equations (6)- (12) are represented by the derivatives of the Cartesian components with respect to the generalized coordinate system. This enables us to avoid using the Christoffel symbol, which often burdens the formulation and computation owing to its complexity and high computational cost.
The present formulation ensures that the number of particles inside each support is sufficient if the generalized coordinate system is defined such that its shape region is appropriately deformed in the physical space. We describe the procedures and guidelines to define the generalized coordinate system later. Notably, the present transformation eventually corresponds to a chain-rule transformation of the Cartesian components, which has been generally adopted for the body-fitted coordinate system in mesh-based schemes for computational fluid dynamics, 52,53 and also for the isoparametric formulation in FE analyses. Equations (6)- (12) are derived only from the definitions of covariant and contravariant bases instead of relying on the chain-rule transformation. Therefore, the present transformation is directly applicable to higher-order tensors that are sometimes required in nonlinear constitutive laws, such as for viscoelasticity. Such a tensor-analysis-based coordinate transformation was adopted in the original GCSPH by Yashiro and Okabe 48 without curved coordinates. The present study directly extends this transformation to the curved coordinates for solid dynamics. Furthermore, the proposed GCSPH is enhanced by an overset method to flexibly vary the particle distribution and avoid coordinate singularities.

SPH approximation
Based on the set of the physical and generalized coordinates defined in the previous subsection, an arbitrary location in the physical space can be mapped to the generalized space as = ( 1 (x), 2 (x), 3 (x)). GCSPH performs a particle approximation in the generalized space; thus, a physical quantity f at an arbitrary point in the computational domain Ω can be represented as where a smoothing parameter h represents the radius of a radially symmetric compact support. Here, W is called the kernel function, and the following cubic B-spline function is adopted: where Suppose that N particles exist in the compact support of the kernel function around particle a at = a . Then, Equation (13) for particle a is approximated as where W ab = W( a − b , h), and m b and b represent the mass and mass density of the neighboring particle b, respectively. In the remainder of this article, the subindices a and b represent the values of particles a and b, for which the Einstein summation convention is not taken. The derivative of f in the i direction at particle a can be approximated in the following two forms: Furthermore, the present study adopts the corrective smoothed particle method 54 to improve the consistency of the particle approximation near the boundary, which aids us to avoid inconsistent representation of the spatial derivatives owing to truncation of the kernel function at the boundary.

Discretization of the governing equations in total Lagrangian GCSPH
In the total Lagrangian framework, the current position of a material point is written as x = x e , and the reference position of the same point is expressed as X = X e , where = {1, 2, 3} denotes the index of the physical coordinate.
Here, (x 1 , x 2 , x 3 ) and (X 1 , X 2 , X 3 ) represent the physical coordinates of the material point in the current and reference configurations, respectively. The deformation gradient F and its rateḞ are defined as The governing equations are established in the reference configuration as where , v, and e are the density, velocity, and specific energy, respectively, and the subscript 0 denotes the values in the reference configuration. Here, J = det |F| represents the determinant of the deformation gradient F. Furthermore, P is the first Piola-Kirchhoff stress tensor, defined as where is the Cauchy stress tensor, which will be represented later in Equation (35). In GCSPH, the governing equations are Equations (20)- (22), where the stress tensor is expressed in terms of the physical coordinate X in the reference configuration, while the spatial derivatives are taken with respect to the generalized coordinate i . Based on the SPH approximation, semi-discrete forms of the governing equations, that is, Equations (20)-(22), at particle a are respectively given by Here, Π ab; represents the components of the artificial viscosity in the physical space and is given by with the definitions of where c a and h a are the speed of sound and the smoothing parameter at particle a, respectively. The SPH approximations of the deformation gradient F and its rateḞ at particle a are expressed as Equations (25), (32), and (33) contain i ∕ X for the coordinate transformation between the physical and generalized coordinates in the reference frame. The coordinate transformation matrix at particle a is computed as Hence, i ∕ X at particle a can be numerically computed as the inverse of this matrix without using the analytical coordinate transformation. Let us consider a two-dimensional quarter cylinder as an example. If the cylindrical coordinate is adopted as a generalized coordinate in the entire domain, the particles in the vicinity of the origin (X 1 , X 2 ) = (0.0, 0.0) in the Cartesian coordinate show a highly dense distribution, while the particles in the far-field show a dilute distribution as shown in Figure 2. In the generalized coordinate, the particles are placed uniformly. For example, in the computation of the momentum conservation Equation (25), the construction of the neighboring list and the spatial derivatives of the spherical kernel function are computed in the generalized coordinate. This kernel gradient in the generalized coordinate can be F I G U R E 2 Schematic of solving the momentum conservation using the Cartesian and generalized coordinate systems transformed into the Cartesian coordinate using the chain rule, This means that the spherical kernel in the generalized coordinate is deformed in the Cartesian coordinate. This procedure allows us to solve Equation (25) in the Cartesian coordinate systems with arbitrary particle arrangements, while always using spherical kernels in the generalized coordinate systems with uniformly arranged particles. Meanwhile, the center axis of the cylindrical coordinate cannot be uniquely defined in the generalized space, that is, coordinate singularity. To overcome this problem, we introduce the overset method.

Constitutive equations
The Jaumann rate ∇ of the Cauchy stress is employed 55 in the present study, such thaṫ where C e is the elasticity matrix. Furthermore, D and W are the rate of deformation and spin tensors, defined as which can be computed using Equations (32) and (33).

Rankine criterion
Three-dimensional crack propagation problems are presented in Section 3.1, where the Rankine criterion is used to handle brittle crack propagation. 38 The interaction factor f ab based on the damage index D ab is defined to characterize the interaction state between particles a and b as The damage evolution 38 is determined according to the following criterion: where (r ab ) t and (r ab ) 0 are the distances between particles a and b in the current and reference configurations, respectively, and max is the failure parameter. At the beginning of the simulation, there is no damage in the model, that is, D ab = 0 and f ab = 1, for all the interaction pairs. When damage initiates, that is, D ab = 1 and f ab = 0, a crack surface is generated by removing the interaction between particles a and b so that a brittle failure is modeled in the simulation.

Thermo-visco-plastic behavior
In this study, the Johnson-Cook model is used to consider plastic hardening, rate dependency, and thermal softening 26,56 in the Taylor and ballistic penetration tests. The yield stress y of this model is expressed as where A is the initial yield stress of the material, and B, C, m, and n are the hardening parameters. Here, n pl is the accumulated plastic strain, anḋ * pl is defined aṡ * wherėp l anḋ0 are the effective plastic strain rate and reference strain rate, respectively. Furthermore, the nondimensional temperature is defined as where T, T r , and T m are the current, room, and melting temperatures of the material, respectively. The increase in temperature caused by plastic work is estimated as where Δw p , C p , and are the increment in plastic work, specific heat capacity, and empirical constant set to = 0.9, respectively. The von Mises yield criterion y f = √ 3J 2 − y is employed to determine whether the stress state exceeds the yield surface, where J 2 = S ∶ S∕2 is the second invariant of the deviatoric stress tensor S = − (1∕3)tr( )I. The Wilkins criterion S n = c f S is used for the return mapping when the trial elastic stress state exceeds the yield surface, where c f = min( y ∕ √ 3J 2 , 1), and S n is the corrected deviatoric stress tensor. Finally, the following equations are used to compute the increment in the plastic strain, increment in the effective plastic strain, and accumulated plastic work density, 26 respectively:

Damage model
An appropriate damage model should be employed to handle fractures in the steel-plate penetration problem, which will be presented in the numerical test discussed in Section 3.3. In this study, the Johnson-Cook model coupled with a damage model is adopted with the following modified yield stress: 57,58 where r is the damage-accumulated plastic strain, given bẏr = (1 − D)̇p l . The damage variable D is determined by the Johnson-Cook criterion as where Δ pl and f are the incremental effective plastic and fracture strain, respectively. In this study, f is calculated as where * = m ∕ eq is the stress tri-axiality with m being the mean stress. Furthermore, D i (i = 1, … , 5) are the material constants.

Overset method
Based on the coordinate transformation between the physical and generalized spaces, GCSPH can handle nonuniform particle distributions while keeping the standard SPH discretization in the generalized space for uniform particle distributions. However, such a coordinate transformation often suffers from coordinate singularities, for example, singularity on an axis of the cylindrical coordinate, which limits the applicability of GCSPH to more complex geometries. We propose the overset method for GCSPH to overcome this problem. A schematic of the overset method is shown in Figure 3. Let us consider a two-dimensional quarter cylinder again. If the cylindrical coordinate is adopted as a generalized coordinate in the entire domain, the particles in the vicinity of the origin (X 1 , X 2 ) = (0.0, 0.0) show a highly dense distribution, and the particle at the origin of the physical space cannot be uniquely mapped to the generalized space. Therefore, in Figure 3, the computational domain is decomposed into two subdomains to differentiate between the inner and outer parts of the cylinder, where the red and blue particles are mapped to Generalized Spaces 1 and 2, respectively * . The cylindrical coordinate is applied to the blue particles, which are defined as the rectangular domain in Generalized Space 2. Meanwhile, the distributions of the red particles are the same for the physical space and Generalized Space 1; thus, the coordinate singularity is removed at the origin (X 1 , X 2 ) = (0.0, 0.0). The calculations of Equation (25) at the red and blue particles are conducted in the physical space with the corresponding coordinate transformation matrix, whereas the kernel function and its derivative are calculated in each generalized space. For communication between the two generalized spaces, transition particles are introduced in the green dashed line, which are mapped to both Generalized Spaces 1 and 2. The interactions between red particles are computed using the coordinate in Generalized Space 1. Similarly, the interactions between blue particles are computed based on Generalized Space 2. However, the interactions between red and blue particles, that is, between red and transition particles, are computed based on Generalized Space 1. This maintains F I G U R E 3 Schematic of particle distributions in an overset method for GCSPH the continuity of the velocity and density fields. In this study, it is assumed that generalized Space 1 is more comprehensive in that the particles related to Generalized Space 2 (cylindrical coordinate) can also be mapped to Generalized Space 1, but not vice versa.
Finally, the black solid closed curves in Figure 3 indicate examples of the domains of influence, each of which is also called a kernel shape in this study, in the physical and generalized spaces. The domain of influence is made circular around a particle in the generalized space for this two-dimensional example. The domain of influence in Generalized Space 2 changes in shape according to the distribution of the blue particles in the physical space. Meanwhile, the shape of the domain of influence remains a circle in the physical space for the red particles, as Generalized Space 1 is the same as the physical space in this example. As such, the domain of influence can be varied such that the resolution of the particles changes adaptively and efficiently depending on the geometry and associated particle distributions.

Algorithm
Algorithm 1 summarizes the computational process of GCSPH in the total Lagrangian formulation with the overset method. The velocity Verlet 59,60 method is adopted for the time integration in this study. In Algorithm 1, the subscripts t, t + Δt, and t + Δt∕2 denote the times of each quantity.

REPRESENTATIVE NUMERICAL EXAMPLES
Three numerical examples are presented using the proposed GCSPH in the total Lagrangian formulation with the overset method. The first example is crack propagation in a pre-notched beam, where the brittle crack is modeled using the Rankine criterion. The second and third examples are the Taylor impact test and ballistic penetration of a Weldox 460 E steel plate, respectively, in which thermo-plastic deformations and plugging failures are examined. In the second and third examples, two generalized spaces are defined in the overset framework to efficiently reduce the number of particles and remove a coordinate-transformation singularity. Algorithm 1. GCSPH in the total Lagrangian formulation with the overset method 1: Define initial conditions, including particle distributions X and in the physical and generalized spaces, respectively; 2: if the number of generalized spaces ≥ 2 then 3: Define transition particles with two neighboring generalized coordinates; 4: end if 5: Find interaction particles in the domain of influence and calculate the derivative of the kernel function W∕ i in each generalized space; 6: Calculate the coordinate transformation matrix between the physical and generalized coordinates and store i ∕ X ; 7: for n = 1 to n end do 8: Calculate the velocity at t + Δt∕2 as v t+Δt∕2 = v t +̇v t Δt/2;

14:
Output the field variables at t + Δt. 15: end for F I G U R E 4 Pre-notched beam subjected to impact load.

Crack propagation in pre-notched beam
As shown in Figure 4 . The failure parameter for the Rankine criterion, max , is set to 0.03. By comparing GCSPH with TLSPH, the validity of the proposed total Lagrangian formulated GCSPH was verified. Basic numerical parameters of TLSPH are the same as those of GCSPH, as mentioned below. Three simulation cases, as shown in Table 1, are conducted. Case (i) is the standard TLSPH with a particle interval of 3.8 mm. Case (ii) is GCSPH with the same particle arrangement as case (i), and it is essentially the same as case (i), although it involves a coordinate transformation between Cartesian and generalized coordinates. Case (iii) is GCSPH with a smaller particle interval in the x and z direction and with the same particle interval in the y direction as cases (i) and (ii). This pre-notched beam example is essentially a two-dimensional problem, and hence, case (iii) is a typical example of GCSPH with high resolution in the xz plane for crack propagation and low resolution in the depth direction. The parameters of the artificial viscosity given by Equation (28) are set as 1 = 1 and 2 = 1, and the time increment is set at Δt = 2 × 10 −5 (ms).
The crack propagation patterns in the final stage are shown in Figure 5. In all the cases, the beam subjected to the impact load shows a bending deformation with crack paths deviating from the vertical plane toward the impact area, in which a mixed mode of crack propagation emerges. These results are consistent with the numerical results in the literature. 23,38 The propagation paths obtained in GCSPH are almost the same as that obtained in TLSPH. This implies that the generalized coordinate transformation with the total Lagrangian formulation is successful and that the proposed method can simulate such a crack propagation even with anisotropic particle distributions. Note that the present example is essentially a two-dimensional problem and involves a single generalized space without the overset method; thus, three-dimensional problems with the overset method will be presented in the following examples.

Taylor impact test
As discussed in the previous subsection, the coordinate transformation of the pre-notched beam from the physical to the generalized space could be readily conducted, although only one set was carried out. However, performing GCSPH with only one generalized space for complex structures is difficult, as demonstrated in Section 2.5. Indeed, the particles near the axis show a coordinate singularity even for a solid cylinder if a cylindrical coordinate is adopted, although those far away from the axis can be converted into a cuboid domain in the generalized space. The overset method proposed in this study overcomes this problem.
To demonstrate the efficiency of GCSPH with the overset method, the Taylor impact test, which has been widely studied through experiments 61 and simulations, 26,62 is discussed in this section. As shown in Figure 6, a cylinder with initial velocity v 0 = 181 (m/s) impacts a rigid wall, where the length and diameter of the cylinder are h = 37.97 (mm) and 2r = 7.595 (mm), respectively. To capture the salient metal behaviors, such as plastic hardening, rate dependency, and thermal softening, the Johnson-Cook model is used. The relevant parameters are listed in Table 2.
In this study, only a half model is considered, and the symmetrical condition is applied to the planes of symmetry of the model. In GCSPH with the overset method, the cylinder is decomposed into two domains to avoid the coordinate singularity on the axis when using the cylindrical coordinate transformation. In Figure 7, the entire model is divided into read and blue domains, which are mapped to Generalized Spaces 1 and 2, respectively, as mentioned in Figure 3. The blue particles in the physical space are distributed in a cylindrical coordinate manner to reduce the number of particles far away from the axis and in the cylinder-axis direction, which corresponds to a cuboid domain with uniform spacing in Generalized Space 2 ( Figure 7D). Meanwhile, the red particles in the physical space are distributed using the Cartesian coordinate in the cross-section to avoid the cylindrical coordinate singularity and have the same distribution in the cross-section in Generalized Space 1 ( Figure 7C). Note that particles in the interface layers between the blue and red particles are defined as transition particles with both Generalized Coordinates 1 and 2. These transition particles correspond to the blue particles surrounding the red particles in Generalized Space 1 ( Figure 7C) and the blue particles in the two bottom layers in Generalized Space 2 ( Figure 7D).
The following five simulation cases are conducted to verify the validity of the GCSPH with the overset method. Case (i) is standard TLSPH with a low resolution, that is, a particle interval Δp = 0.47 (mm). Case (ii) is standard TLSPH with a middle resolution, Δp = 0.38 (mm). Case (iii) is standard TLSPH with a high resolution, Δp = 0.29 (mm). Case (iv) is GCSPH with a low resolution, where the number of particle in the cross section layer n lay = 486, and the number of  particles in the axis direction n axi = 16. Case (v) is GCSPH with a high resolution, where n lay = 1066 and n axi = 21. The detailed condition can be seen in Table 3. Case (iv) has a significantly small number of axial particles, n lay . This reduces the total number of particles, n total , to less than 1/4 of that in the highest resolution TLSPH (case iii), even though the cross-sectional resolution in case (iv) is higher than that in case (iii). It should be noted that this is one of the benefits of GCSPH with the overset method. By allowing for non-uniform particle arrangement, high-resolution calculations can be performed while dramatically lowering computational costs. In the Appendix, the results of GCSPH for various particle arrangement conditions are shown. The accumulated plastic strain and temperature distribution at 0.05 (ms) are shown in Figure 8. It can be seen that in all the cases, the end of the cylinder is deformed like a mushroom shape by plastic deformation due to the impact on the plate, and the temperature is increased due to plastic work. These results are consistent with those in References 26,62. For a more detailed comparison, the cross-sectional shapes of the impacted part at the same time are compared, as shown in Figure 9. Note that the range of the color legend is different from that of Figure 8. In the standard TLSPH, particles are arranged in a rectangular manner, thus the outer edge shape is not exactly circular. Consequently, there is non-uniform  27 Nevertheless, in GCSPH, the particles are arranged using cylindrical coordinate manners, and thus uniform strain distribution along the outer edge was observed. This is also one of the benefits of GCSPH. GCSPH allows particle arrangement along with curved geometries, thereby eliminating the effects of the arrangement, such as the non-uniform strain distribution. Furthermore, GCSPH has numerous cross-sectional particles and a small number of axial particles, enabling high-resolution analysis with a realistic total number of particles. Figure 10 shows the changes in the length, mushroom radius, effective plastic strain, and maximum temperature over time obtained in the lowest and highest resolution TLSPH (cases i and iii) and GCSPH (cases iv and v). Although there is a slight difference in the final length, h, TLSPH and GCSPH provide almost identical results. The results here are validated with the reference data in Table 4. It can be seen from the comparisons that the maximum effective plastic strain and temperature increase in the present study are in good agreement with the previous works of FEM 62 and standard TLSPH; 26 in particular, the final length h obtained by GCSPH is much closer to the experimental result. 61 An important aspect in formulations of solid mechanics by SPH is energy conservative features. Figure 11 shows the changes in the energy over time obtained in the lowest and highest resolution TLSPH (cases i and iii) and GCSPH (cases iv and v). In all the cases, the total energy is fluctuated at the time of the collision. This has been reported in a previous work. 63 The energy conservation is not so good in the low-resolution cases ((i) TLSPH and (iv) GCSPH), and the increment of total energy in GCSPH is larger than that in TLSPH. This might be because performing coordinate transformations lead to loss of accuracy. However, in the high-resolution cases ((iii) TLSPH and (v) GCSPH), GCSPH showed a better conservation feature of total energy, closer to that of TLSPH, despite the presence of oscillations and a slight increase in total energy. To improve energy conservation, Islam et al. proposes a formulation of particle approximation that is highly energy conservative. 63 Furthermore, especially for GCSPH, the accuracy loss due to the chain rule would be improved by describing the derivative in the strong conservative form, which is used in computational fluid dynamics. 64 These methods will be implemented into GCSPH in the near future.

Case (i) TLSPH (ii) TLSPH (iii) TLSPH (iv) GCSPH (v) GCSPH
These results confirm the validity of GCSPH with the overset method and show the benefits of GCSPH, such as high-resolution analysis and particle arrangement along with the curve shape. In the next section, we describe the three-dimensional ballistic penetration simulation, which is a significant effective example of the use of GCSPH.

Ballistic penetration of steel plates
The final example is a penetration test of the Weldox 460 E steel plate, 58,65 which is schematically shown in Figure 12. Such penetration tests have been addressed using two-dimensional SPH simulations with axisymmetric idealization in a previous work. 28 However, impact tests on anisotropic materials, such as carbon-fiber-reinforced-plastic laminates, require three-dimensional analysis. 66,67 Therefore, in this study, a penetration test on a steel plate, which has been conducted in previous studies, is analyzed using three-dimensional GCSPH simulation to validate GCSPH and to extend it to  three-dimensional analysis. Here, the plate has a diameter of 500 (mm) and two thicknesses of 8 and 10 (mm). It is modeled using the damaged thermo-visco-plastic material model. The material parameters are listed in Table 5. A projectile with a diameter of 20 (mm) and a length of 80 (mm) is modeled as a rigid body with the same density. As the material undergoes considerable plastic deformation, high strain, and thermal softening, the Johnson-Cook model with damage, as discussed in Section 2.4.3, is used to describe the dynamic behavior of the steel plate, where the fully damaged particles (D = 1) are removed during the simulation. The material constants are set as follows: 26  For the impact problem, the plate's contact zone and the projectile are of greater importance than the region away from the contact zone. In the standard SPH, the space of neighboring particles is kept uniform to obtain accurate results, which often results in a significant increase in the computational cost. Thus, GCSPH with the overset method is expected to be effective for this plate penetration problem with a smaller number of particles. In the physical space, numerous particles are distributed in the vicinity of the impact zone on the plate, while a few particles are distributed far away from the impact zone. Analogous to the Taylor impact test, only a half model is considered, and symmetrical conditions are applied to the plane of symmetry. Two generalized coordinates are used for the plate, while only one generalized coordinate is used for the projectile. As shown in Figure 13A,B, the plate with a thickness of 8 (mm) is discretized with 40,821 particles, in which 41 layers of particles (red set, 27,511 particles) in the radius direction are placed within a radius of 30 (mm), and 10 layers of particles (blue set, 13,310 particles) are placed far away from the contact zone. Furthermore, 11 layers of particles are placed in the thickness direction. The particle distributions in the generalized spaces are shown in Figure 13C,D. By taking advantage of the overset method, the plate with a thickness of 10 (mm) can also be discretized with the same scheme. If TLSPH is performed with the same particle spacing as the contact zone of GCSPH, 1.9 million particles are required, which is practically impossible in standard computing resources.
The steel plate penetration problem has been widely studied by many researchers, and predictions of the residual velocity of a projectile and the velocity of a plug are common topics of interest. The responses of a plate subjected to an impact load are local bending deformation, thermal softening, accumulated damage, and adiabatic shear band, which eventually lead to plugging failure. Figure 14 shows the changes in the residual and plugging velocities over time in the two impact velocity cases, and the velocity distributions obtained with various impact velocities v 0 are shown in Figures 15 and 16. Figures 15 and 16 show the snapshots at the time at which the projectile fully penetrates, that is, at which the bottom of the projectile becomes lower than that of the remained plate. It can be seen from these figures that the bending deformation and plugging failures of the plate are caused by the projectile; a larger impact velocity v 0 leads to a shorter time to full penetration. Consequently, larger impact velocity v 0 leads to minor bending deformation. The simulated projectile residual velocity and plugging velocity are compared with the experimental results 65 in Tables 6 and 7. Although there are some over/underestimations, GCSPH with the overset method can reproduce three-dimensional impact fracture efficiently.

CONCLUSIONS
This study proposed GCSPH with an overset method using the total Lagrangian formulation, in which the coordinate transformation between the physical and generalized (local coordinate) spaces was properly employed to control the spatial resolution and reduce the number of SPH particles. The main results can be summarized as follows.
1. The proposed GCSPH allows non-uniform particle distributions in the physical space to vary the spatial resolution, whereas the SPH particles are uniformly distributed in the generalized space, as in the standard SPH. The governing equations are solved in the physical space with the corresponding coordinate transformation matrix, whereas the kernel function and its derivative are calculated in the corresponding generalized space. The formulation adopted in this method is an extension of the original GCSPH 48 to curved coordinates by means of tensor analysis. Furthermore, this study is the first to generally represent the vector and tensor formulations with their gradients in the SPH discretization for curved coordinate systems, which is useful for implementing higher-order tensors for more general continuum mechanics. 2. The proposed GCSPH is enhanced by the overset method, where the computational domain is decomposed into multiple subdomains with local coordinates. The overset method enables us to flexibly control the spatial resolution so that, for example, the non-contact zone in ballistic penetration of a steel plate contains a smaller number of particles compared to the contact zone. Furthermore, coordinate singularities, which are inevitable in the cylindrical or spherical coordinate system, can be avoided by using the singularity-free coordinate system in the local regions. 3. To alleviate tensile instability, total Lagrangian formulations were implemented for finite deformation problems for the first time in the SPH formulation with coordinate transformation techniques. 4. Numerical examples of three-dimensional brittle crack propagation, the Taylor impact test, and plugging failure of the Weldox 460 E steel plate were studied using GCSPH. By comparing it to standard TLSPH, it was confirmed that GCSPH can perform stable computations for crack propagation and impact problems and has the same energy conservation feature as TLSPH. The numerical results were in good agreement with literature. 5. The strong two benefits of GCSPH were confirmed: the first is eliminating the effects of the arrangement such as the non-uniform strain distribution on the outer edge caused by particle arrangement and curved geometries. The second is the effective reduction of computational cost by changing the particle distribution in the physical space. These benefits enable GCSPH to analyze three-dimensional problems with high resolution, which is very costly to perform with the standard SPH.
In summary, this is the first study to not only formulate the SPH in the generalized coordinate system with the total Lagrangian framework but also enhance it with the overset method for solid mechanics problems, which allows more flexible static resolution control and stable numerical simulation. However, the inconsistency in zero-order approximation with an arbitrary distribution of particles has been resolved by the so-called corrected SPH (CSPH), 29,30 which achieves the first-order consistency in the particle approximation and ensures the conservation of the linear and angular momenta in the governing equations. We believe that CSPH will be applied to generalized (curvilinear) coordinates in the future. As a further extension, we are planning to combine GCSPH coupled with the simplified marker and cell (SMAC) algorithm to analyze the large deformation of hyperelastic materials 68 using a system that has a local high resolution and low computational cost with the overset method.

APPENDIX
To demonstrate the accuracy of GCSPH, the number of particles in the axial direction and each layer (cross section) is varied, where the smoothing factor and artificial viscosity parameters are set as = 1.1, 1 = 1, and 2 = 1. The accumulated plastic strain and temperature distribution in the final stage are shown in Figures A1 and A2, with the number of particles increasing in the axial direction and each layer, respectively. It can be observed that for all cases, the end of the cylinder is deformed like a mushroom shape by plastic deformation due to the impact on the plate, and the temperature is increased due to plastic work. In Figure A1, the length and mushroom radius are varied between 34.63 (mm) ≤ h ≤ 34.73 (mm) and 5.26 (mm) ≤ r ≤ 5.34 (mm), respectively, by varying n axi between 16 and 31. The number of particles in the axial direction does not significantly affect the results, and the differences in h and r are lower than 2% of their largest values. A similar discussion arises for Figure A2 in assessing the effect of n lay , that is, the cross-sectional resolution, on the final state, where n lay was varied from 317 to 1066.