###### Abstract
Additionally, the accuracy of such contact algorithms is affected by the determination of particle size, and they usually require relatively high particle resolution to obtain satisfactory accuracy and do not allow large inconsistency of particle spacing for different contact objects [10]. [PERSON] et al. [11] proposed a contact algorithm which uses a contact potential to impose contact conditions. The algorithm was implemented in the SPH method and applied to the simulation of various problems, including fluid sloshing and hypervelocity impact. It can preserve the meshless feature of the SPH well. However, it did not consider the friction effect. Meanwhile, it may produce contact force according to the contact potential when objects are sufficiently close, even if they are not actually in contact. [PERSON] et al. [12] developed a frictional contact algorithm to simulate the interaction between soils and rigid or deformable structures. [PERSON] et al. [13] presented a hybrid contact method to simulate frictionless and frictional contacts between soil and a rigid body and contacts between rigid bodies. The above two algorithms can be considered particle-to-surface contact algorithm. They both used the penalty method to enforce contact conditions. They defined surfaces for partial objects (generally rigid or small deformation objects) with boundary particles in the initial discretization, but without considering a possible change of topology of surfaces during the deformation process. Thus, they are mainly suitable for problems where at least one of the two contact objects is undeformable or does not experience very large deformation and material failure. In addition to this, they need to assume the offset from boundary particles (essentially particle size) to locate the actual surface positions. Therefore, their accuracy has an effect of the determination of the offset. [PERSON] et al. [14] presented a mortar segment-to-segment contact algorithm for frictionless contact problems and implemented it for two-dimensional (2D) problems. In this algorithm, surfaces are defined and located in a similar way in [12] for all objects, and the contact conditions are imposed by using the perturbed Lagrange method. This algorithm can provide accurate results for small and moderately large deformation problems. However, it is implemented in the total-Lagrangian SPH method and does not address significant surface changes, such as surface fragmentation, which may invalidate the surfaces initially defined. Thus, it is not applicable to problems with extremely large deformation and material failure. [PERSON] et al. [15] presented an algorithm with the artificial density concept. The algorithm was designed for the treatment of rigid-rigid contacts. Furthermore, a series of algorithms [16; 17; 18] was proposed for the SPH method to treat fluid-rigid/flexible body interactions. These algorithms are not suitable for treating contacts between solids, especially when solid surfaces have intense deformations. In our previous work [10], a versatile penalty-based surface-to-surface contact algorithm was presented. It can treat both frictionless and frictional problems and deal with problems with arbitrary deformation, from small deformations to extremely large deformations. It can reach higher accuracy and better stability than the particle-to-particle contact algorithm. The algorithm reconstructs surfaces automatically instead of predefining surfaces with connectivity between boundary particles in the initial phase of simulation and, therefore, preserves the advantages of the SPH in modeling large deformations. However, the algorithm presented in the previous work is limited to 2D problems. It is not satisfactory for many practical applications, since few practical problems can be simplified as 2D problems.