## 2 Materials and Methods
### Finite Element Modeling  
Agricultural soils experience plastic deformations after penetration induced by an engaging tillage tool. The behavior of the soil before yielding is elastic, but the elastic range is usually small and does not represent the real stress range [27] (i.e., the stresses applied by tillage tools are beyond the elastic range). Therefore, defining the soil material in an FEM model consists of importing both the elastic and plastic properties and selecting the yield and plastic strain behavior using constitutive models. The Drucker-Prager (D-P) model [28] was used in this study, it is a modified version of von [PERSON] model, and considers the influence of hydrostatic pressure in failure. The extended D-P models are used to model frictional materials, such as soils and rock, where material yield is associated with hardening (i.e., the material strength increases with increasing stress level). The D-P model is defined as follows:  
\[\mathrm{F}=\mathrm{t}-\mathrm{p}\times\mathrm{tan}\ \beta-\mathrm{d}=0 \tag{6}\]  
where F is the yield function, t is the deviatoric stress, \(\mathrm{p}\) is the normal stress, \(\beta\) is the internal friction angle, and d is the cohesion of the material. Moreover, \(\mathrm{p}\) and t are given as follows [29]:  
\[\mathrm{p}=\frac{1}{3}(\sigma_{1}+\sigma_{2}+\sigma_{3}) \tag{7}\] \begin{table}
\begin{tabular}{c c c c} \hline  & **Soil 1** & **Soil 2** & **Soil 3** \\ \hline Bulk density (g/cm\({}^{3}\)) & 1.98 & 1.53 & 1.22 \\ Moisture (\(\%\)) & 0 & 0 & 25 \\ Cohesion (Kpa) & 9.8 & 13.5 & 10.2 \\ Internal friction angle (\({}^{\circ}\)) & 30.8 & 28.2 & 25.5 \\ \hline \end{tabular}
\end{table} 
Table 2: Physical properties of experimental soils. Figure 3: (**a**) The soils used in experiments. (**b**) The devices in experiments. (**c**) The simulation model for analyzing penetration process which was built in ABAQUS. \[t=\frac{1}{2}q\left[1\right.+\frac{1}{k}-(1-\frac{1}{k})\left.(\frac{r}{q})^{3}\right] \tag{8}\]  
\[q=(\sigma_{1}-\sigma_{3}) \tag{9}\]  
where k (flow stress ratio) is the ratio of tension yield stress to compression yield stress in triaxial tests (0.778 \(<\) k \(<\) 1). k and \(\beta\) were calculated as follows:  
\[K=\frac{3-\sin\varphi}{3+\sin\varphi} \tag{10}\]  
\[\tan\beta=\frac{6\sin\psi}{3-\sin\psi} \tag{11}\]  
where \(\varphi\) is the internal friction angle and \(\psi\) is the Dilatancy angle. The elastic modulus for each soil was determined by unconfined uniaxial compression test [30]. The Poisson's ration \(\
u\) was calculated as follows:  
\[\
u=\frac{1-\sin\psi}{1+(1-\sin\psi)} \tag{12}\]  
The soil cohesion and the angle of internal friction were measured through direct shear test. The soil property parameters used in FEM model were shown in Table 3.  
The 3D models were developed in ABAQUS/Explicit as shown in Figure 3. The soil model part was 200 mm long, 200 mm wide, and 100 mm deep. The axis of each specimen was vertical in the center of the soil. Each specimen was considered as a rigid body and meshed with R3D4 element type. And the soil was meshed with C3D8R element type, which is commonly used for 3D stress-strain analysis of continuum materials [31]. As reported, the mesh size importantly affects the tillage force and calculation time [1]. The soil model had a fine mesh with 307,200 elements which can ensure the accuracy of the results. Meanwhile, Eulerian Lagrangian re-meshing coupled with mesh control was utilized to preserve a high-quality mesh of the soil model during the specimen penetration and prevent excessive mesh distortion [32]. This is an important feature in problems with high degree of material deformation. The interaction between specimen and soil was simulated with a general contact law and tangential behavior. A penalty contract method was used which searched for node-to-surface penetrations in the current configuration. The rigid body was selected as the master and the deformable part as the slave automatically by the software.