## 3 Numerical Model
### Load Parameters  
As described in the following Equations (1) and (2), during vibratory pile driving, the force on the pile \(F_{\text{P}}\) can be divided into three parts: (a) the vibrating force \(F_{\text{V}}\) developed by a set of counter-rotating eccentric mass in a vibrator, (b) the static force \(F_{\text{S}}\), which is the resultant force of supporting force and pile-vibrator deadweight, and (c) the resistance of pile penetration \(F_{\text{r}}\), which is the reaction force from the soil.  
\[F_{\text{P}}=F_{\text{V}}+F_{\text{S}}+F_{\text{r}} \tag{1}\]  
\[F_{\text{v}}=F_{\text{vm}}\sin(2\pi ft) \tag{2}\]  
where \(F_{\text{vm}}\) is the vibrating force amplitude, \(f\) is the vibrator frequency, and \(t\) is the driving time.  
The vibrator frequency \(f\) was 20 Hz, which is in line with the experiment. \(F_{\text{vm}}\) and \(F_{\text{s}}\) were assigned as 4.0 \(G\) and 0.2 \(G\), respectively, where the pile deadweight \(G\) was 34.5 N in this paper. The penetration resistance \(F_{\text{r}}\) was updated by summing up the pile-soil contact force during pile penetration.  
The dynamic analysis in FLAC3D is based on an explicit finite difference scheme, in which the model is solved numerically using an explicit finite difference formulation in time. The maximum stable time increment \(\Delta t_{m}\) is calculated by Equations (3) and (4) [21]:  
\[\Delta t_{m}=0.5 min\left(\frac{V}{C_{p}A_{max}^{f}}\right) \tag{3}\]  
\[C_{p}=\sqrt{\frac{E(1-v)}{\rho(1-2v)(1+v)}} \tag{4}\]  
where \(C_{p}\) is the compression wave speed, \(V\) is the tetrahedral subzone volume, \(A_{max}^{f}\) is the maximum face area associated with the tetrahedral subzone, the \(min\) ( ) function is applied to all zones.  
A smaller \(\Delta t_{m}\) brings more computation time. In this paper, the original \(\Delta t_{m}\) was around 1.6 \(\times\) 10\({}^{-8}\) s due to the excessive elastic modulus (206 Gpa) of the pile and the tube. In order to decrease the time cost, the density scaling method [6; 11] was used. In this method, referring to Equations (3) and (4), densities of the pile and the tube were  
\begin{table}
\begin{tabular}{c c c c} \hline \hline \multirow{2}{*}{**Type**} & **Normal Stiffness** & **Shear Stiffness** & **Friction Angle** \\  & \(k_{s}\) [GPa-m] & \(k_{s}\) [GPa-m] & \(\delta\) [\({}^{\circ}\)] \\ \hline Pile–soil & 0.2 & 0.2 & 10 \\ Tube–soil & 8 & 8 & 0 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Contact parameters in the numerical model.  
scaled down to increase \(\Delta t_{m}\), and the \(F_{\rm P}\) was also enlarged by the same scale to avoid the abnormal acceleration, see Equation (5) [6]. The alternative to decrease computation time is a mass scaling method, which can be used in finite element software like ABAQUS [25]. The difference between the two methods is that the pile acceleration will remain normal by enlarging the \(F_{\rm P}\) in the density scaling method.  
\[a_{\rm P}=\frac{ngF_{\rm P}}{nG} \tag{5}\]  
where \(n\) is the scaling factor, \(g\) is the gravity acceleration of 10 m/s\({}^{2}\).  
To determine the scaling factor (SF), a series of SF were tested. Figure 4 presents the relation between the SF and the \(\Delta t_{m}\). It is clear that the \(\Delta t_{m}\) stops increasing when the SF is greater than \(10^{5}\). In the present study, SF = \(10^{6}\) was applied, and correspondingly the \(\Delta t_{m}\) was \(3.1\times 10^{-6}\) s.