## 3 3D SPH Surface-to-Surface Contact Algorithm Using Surface Reconstruction  
### 3D Surface Reconstruction of SPH Objects  
In SPH simulation, objects are discretized by particles without explicit connectivity as that defined for nodes in the FEM, so surfaces of objects are not directly available. To retrieve the surfaces of 3D objects for contact treatment, surface reconstruction is adopted. The general procedure of surface reconstruction is illustrated in Figure 1. Firstly, Delaunay triangulation is performed for the particles representing an object, and a series of tetrahedrons, which connect the particles, is obtained. Then, the circumradius for each tetrahedron is calculated, and the tetrahedrons with circumradius larger than a critical distance are identified and dropped. Finally, the four triangular faces of each remaining tetrahedron are checked to determine whether they are on the surfaces of the object. If a triangular face is shared by only one tetrahedron, then it is a boundary face and becomes a part of the object surfaces. All these triangular faces shared by only one tetrahedron compose the reconstructed surfaces of the object. Here, the Delaunay triangulation is used in surface reconstruction because of its distinct features and easy availability of its construction algorithm and program module. It is known that the Delaunay triangulation is unique for a point set of general position (no 4 points on the same plane and no 5 points on the same sphere in 3D space). Meanwhile, the regularity of tetrahedrons obtained with the Delaunay triangulation is generally good. Thanks to these features, reconstructed surfaces based on the Delaunay triangulation will be not significantly affected by the triangulation process. Rather, they will consist of triangular faces of good shapes. The Delaunay triangulation algorithms have been widely investigated, and some program modules are available for use. In this work, the \"delaunayTriangulation\" function in MATLAB [24] is embedded in the main program coded in Fortran to achieve a robust 3D Delaunay triangulation. The Delaunay triangulation can only produce the convex hull of a point set. It may contain tetrahedrons connecting points with large distances which are not part of the volume of an object. In the surface reconstruction, such tetrahedrons, along with some skinny tetrahedrons, are eliminated by checking that the circumradius of a tetrahedron is larger than a critical distance specified according to local particle spacing. When considering whether a tetrahedron is dropped during the surface reconstruction process, it is required to determine the circumradius of the tetrahedron and the critical distance. The circumradius of a tetrahedron \(I\) is calculated as  
\[R_{I}=S_{I}/6V_{I}, \tag{8}\]  
where \(S_{I}=\sqrt{p(p-l_{12}l_{34})(p-l_{13}l_{24})(p-l_{14}l_{23})}\), \(p=(l_{12}l_{34}+l_{13}l_{24}+l_{14}l_{23})/2\), \(l_{12}\), \(l_{13}\), \(l_{14}\), \(l_{23}\), \(l_{24}\) and \(l_{34}\) are the lengths of the edges of the tetrahedron, 1, 2, 3 and 4 denote the local numbers of particles at the four vertices of the tetrahedron, and \(V_{I}\) is the volume of the tetrahedron which is calculated as  
\[V_{I}=\frac{1}{6}\mathrm{det} \begin{pmatrix}1&x_{1}&y_{1}&z_{1}\\ 1&x_{2}&y_{2}&z_{2}\\ 1&x_{3}&y_{3}&z_{3}\\ 1&x_{4}&y_{4}&z_{4}\end{pmatrix} . \tag{9}\]  
The critical distance is determined as \(\alpha\overline{\Delta}_{I}\), where \(\alpha\) is a constant and \(\overline{\Delta}_{I}\) is the local particle spacing. In this work, the constant \(\alpha\) is 1.0, and the local particle spacing \(\overline{\Delta}_{I}\) is estimated by the smoothing lengths of the four particles at the vertices of the tetrahedron as \[\overline{\Delta}_{I}=(h_{1}/\alpha_{\rm h}+h_{2}/\alpha_{\rm h}+h_{3}/\alpha_{\rm h }+h_{4}/\alpha_{\rm h})/4. \tag{10}\]  
It should be mentioned that care must be given when the circumradius is determined from Equation (8).