Abstract
The soil-structure interface problem is an important part of soil-structure interaction research. These problems are mostly three-dimensional space problems, which is more complex to solve. In this paper, reduced stress and strain rate vectors are incorporated into the explicitly granular hypoplastic model by considering the plane strain state precisely. In addition, considering the important influence of roughness on the mechanical properties of contact surface, an improved hypoplastic model is established by incorporating the influence of roughness into the hypoplastic model, and the applicability of the new improved model is validated by comparing with the simulation results of the Mohr–Coulomb model, the explicitly granular hypoplastic models, and the experimental data. The results indicate that the improved model can be utilized to reflect the nonlinearity of the mechanical properties of the contact surface, which is in good agreement with the experimental data.
1. Introduction
The soil-structure contact surface is an important medium of transferring stress and strain between soil and structure. The unique mechanical properties of the contact surface have great influence on the strain and interaction between soil and structure. The study of mechanical properties of soil-structure interface has always been one of the important topics in civil engineering and geotechnics. In particular, interface modeling is an important part of simulating soil-structure contact behavior. Therefore, many researchers have proposed a large number of constitutive models.
In recent years, many interface models have been available, using elastic models, elastoplastic models [1–4], damage models [5, 6], and so on. Lashkari [1] proposed a plasticity model for the mechanical properties of sand-structure interfaces based on the critical state soil mechanics and the bounding surface plasticity theory. Boulon et al. [2] developed a strain-hardening elastoplastic model for the behavior of sand-structure contact surface under monotonic and cyclic loadings based on two plastic surfaces. Based on the disturbed state concept (DSC), Desai et al. [6] proposed a plastic model for liquefaction and stress-strain response of sand-concrete interfaces under dynamic conditions. Karalar and Cavusli [7] have used special material models (a WIPP-creep viscoplastic material model and a burger-creep viscoplastic material model) and studied the influences of the normal and shear stiffnesses on viscoplastic creep behavior of a concrete-faced rockfill (CFR) dam-foundation interface. Based on WIPP-creep viscoplastic material model and burger-creep viscoplastic material model, Karalar and Cavusli [8] have conducted the creep analyses of the concrete slab and rockfill materials-foundation, respectively. Karalar and Cavusli [9] have modelled the Oroville dam using Flack3D and studied the effects of far fault earthquakes on the mechanical properties of dam body-foundation interface. However, these models are described by complex formulation and many parameters. Meanwhile, traditional elastoplastic theories such as flow rule, destruction rule, and hardening rule have to be taken into account, and many assumptions must be selected.
Therefore, the hypoplastic models have attracted wide attention [10–14]. Fuentes et al. [10] developed a new constitutive model for sands by incorporating a loading surface with two hardening variables into a “Karlsruhe” hypoplastic model to reproduce memory effects and simulate plastic accumulation under cyclic loading. Based on a hypoplastic model with intergranular strains, Hleibieh et al. [11] presented a numerical model for a centrifuge experiment on tunnel embedded in sand under seismic loading. Li et al. [13] proposed a new macroelement developed within the framework of hypoplasticity to simulate the mechanical property of single vertical piles in sand. Based on the hypoplastic model by Gudehus [15] and Bauer [16], von Wolffersdorff [14] developed a new improved hypoplastic model with the asymptotic behavior. The models mentioned above are based on the nonlinear tensor function and written in few variables and a simple formulation. These models do not require elastoplastic theories such as yield surface, plastic potential, flow, and hardening laws and were developed without separating the deformation of the elastic and plastic parts. Therefore, compared to elastoplastic models, the hypoplastic models are more consistent with the actual situation and have obvious advantages in simulating the geotechnical mechanical properties.
In general, there are many factors affecting the mechanical properties of soil-structure interface, especially normal stress, relative density of soil, normal stiffness, structural surface roughness and hardness, and particle breakage, which have important effect on the mechanical property of the soil-structure contact surface [17]. Many researchers have conducted numerous experiments and simulations to prove the effect of roughness on the mechanical properties of contact surface [18–24]. Martinez and Frost [18] have performed laboratory experiments for friction sleeves sheared against sands to study the effect of surface roughness on the behavior of sand-structure contact surface. Haruyama [20] has conducted drained triaxial compression tests for the shear characteristics of uniform spheres of steel having different surface roughnesses to study the influence of surface roughness on the development of failure zone in the triaxial specimen. Wang et al. [23] developed the silicon piezoresistive sensor and used it to study the effect of roughness on the mechanical properties of the silty clay-concrete contact surface. Jensen et al. [24] have proposed an improved discrete element method (DEM) to simulate the mechanical property of granular media-structure interface and prove the effects of roughness and particle shape on the shear response of contact surface. In order to provide a consistent interface behavior model, the hypoplastic model takes account of the frictional contact between soil and structure with different roughnesses. However, roughness has been taken into account in different ways [25–28]. Arnold and Herle [25] have proposed a hypoplastic model by incorporating the parameter κr to describe the contact surface roughness and the parameter κr was depended on the soil-soil friction angle and the interface friction angle. Stutz et al. [27] have considered the parameter κr to depend on the contact surface roughness and established an enhanced hypoplastic interface model. In addition, hypoplastic models are mostly three-dimensional models that deal with space problems and require a lot of time and work. Therefore, some researchers have proposed two-dimensional hypoplastic interface models to reduce the CPU load [25, 27]. Arnold and Herle [25] have proposed the hypoplastic model for 2-D interface conditions and Stutz et al. [27] have developed an enhanced hypoplastic model by the incorporation of in-plane stresses into the model.
The research methods for the interface problems mentioned above are reviewed from the aspect of macromechanics. However, micromechanical analysis is also significant. Therefore, many researchers have conducted atomic scale simulation methods to study the interface problems. Based on the molecular dynamic modeling, Ghaffari et al. [29] carried out micronumerical simulations of mechanical properties of lubricant-sliding solids interface to study the effects of temperature and n-alkanes on the mechanical response of contact surface. Bai et al. [30] have used molecular dynamics model to simulate the friction and wear between diamond-like carbon (DLC) films and diamond tips with respect to micromechanical view and studied the effects of load, velocity, and surface roughness on the micromechanical response of contact surface.
Based on the existing granular hypoplastic model [14], in this paper, a new improved model has been developed by taking into account the relative surface roughness. In addition, the improved model was introduced with the reduced stress and strain tensors obtained by considering the plane strain condition precisely. Finally, the reliability of the improved interface model was validated by comparing with simulation results of the previous models and experimental data for different roughnesses.
2. Hypoplastic Interface Model
2.1. Stress and Strain Rate Tensor
Based on 2-D interface conditions, Arnold and Herle [25] proposed a 2-D hypoplastic interface model under the assumption of an isotropic stress state at the interface, while Stutz et al. [27] developed an enhanced model under the assumption of a transverse isotropic stress state at the interface.
In this paper, the new improved model was developed by considering the plane strain state precisely.
Figure 1 shows the soil-structure contact surface and the coordinate system. , , and are the displacement components in -, -, and -directions, respectively. As shown in Figure 1, the size in one direction is much larger than two sizes in the other two directions, while the geometrical shape and size in the lengthwise direction do not change. The external force (normal stress ) also acts parallel to the cross-section and does not change lengthwise. Therefore, the geometrical characteristics and external force of the soil-structure interface are the same as those of the plane strain problem. If we assume that any cross-section is the plane, any line vertical to it is the -axis, and the -direction is infinitely long, then , , , , , and do not change along the -direction and are only functions of and . As any cross-section can be considered as the plane of symmetry, then and the displacement vectors of all points are parallel to the plane.

In the plane strain state, the reduced stress tensor can be defined as
The degenerated vectorial is written as
As shown in equation (2), the first input of the stress vector is the stress normal to the interface , the second and third inputs of the stress vector are the stresses normal to the interface in - and -directions and , respectively. The shear stress is . As mentioned in equation (1), the shear stresses are and ; therefore, these are not incorporated into the vectorial notation.
The reduced strain rate tensor is written as
The vectorial form is defined:
In the plane strain state, the strain in the -direction is , so the third input of the strain rate vector is . In order to match the stress vector to the strain vector, this is used in the strain rate vector. These reduced stress and strain rate tensors will be used to simulate the mechanical property of soil-structure interface.
2.2. Constitutive Formulation
The constitutive formulation used was defined by von Wolffersdorff [14]:withwhere and are the objective stress rate and stretching tensor, respectively; is a deviator stress and the coefficient is defined by critical friction angle ; is simply estimated from the angle of repose of loose material. According to the Matsuoka–Nakai [31] failure condition, the stress coefficient is written bywith the lode angle ,where is the deviator of ; is the second-order unity tensor and ; is the barotropy factor which controls the effect of the mean stress and is given as
The pycnotropy factors and control the influence of the relative density and are defined as
Pressure-dependent void ratios are defined bywhere , , and are limiting void ratios; , , and are maximum void ratio, critical void ratio, and minimum void ratio at zero pressure, respectively. Under increasing mean pressure , the limiting void ratios , , and decrease until the limiting values , , and are reached (Figure 2). In Figure 2, . The granular stiffness is the only parameter with the dimension of stress and is considered as the reference pressure. The exponent reflects the pressure-sensitivity of a grain skeleton. In addition, controls the relative density to peak friction dependency and controls the relative density to the soil stiffness dependency.

2.3. The Contact Formulation considering Plane Strain Condition and the Surface Roughness
Arnold and Herle [25] assumed that the three normal stresses acting on the interface are equal, i.e., soil has isotropic properties, which do not match actual condition of the soil. Natural soils generally exhibit transverse isotropy with an axis of rotational symmetry perpendicular to the deposition plane. As a result, equations (1) and (2) were redefined as follows:
The following formulation is defined by Dove and Jarrett [33] and Jin et al. [34]:where and are respectively the maximum friction coefficient of the sand-structure contact surface and sand under the same normal stress and is the surface roughness.
Based on Rowe’s stress dilatancy theory, Bolton [35] proposed a relationship between the maximum friction angle , maximum dilatancy angle , and critical friction angle :
Combining equations (13) and (14), the critical friction angle of contact surface is written as
By substituting equation (15) in equation (6), it can be seen that the surface roughness has been included in the formulation of granular hypoplastic model.
3. Simulation Method
3.1. Simulation Model of Direct Shear
The direct shear simulation model was used to verify the new model by comparing with previous models. The size of the soil sample is 40 10 10 cm (Figure 3). The dimension of the structural part is 54 18 5 cm with the following characteristics: Young’s modulus of E = 1 GPa and Poisson’s ratio of . The structures of different roughnesses were used. The roughness R is defined as the relative height between the highest peak and the lowest along a profile of structure surface. Figure 4 shows the profile of the applied roughness. As shown in Figure 4, the structural surface consisted of regular triangular asperities in the numerical simulation. Each regular triangle is an isosceles triangle and is exactly the same to each other; it contains two geometric parameters: the asperity height (R) and the asperity width . In this test, the width was set at 2 mm and height (roughness) was set according to the condition of experiment and simulation. The assembly was meshed by a linear interpolation of eight noded elements (C3D8). The entire loading process is divided into two stages; the first stage applies normal stress on the top surface of the soil and, in the second stage, the soil is horizontally displaced by u = 7 cm in the positive direction. The friction characteristics of the contact surface are specified by the frictional subroutine (FRIC). The simulation was compared with the standard Mohr–Coulomb model and the explicitly granular hypoplastic models.

(a)

(b)

(c)

3.2. Simulation Scheme
During this research, we use the frictional subroutine (FRIC keyword in ABAQUS) to compile the new model into the FRIC subroutine; the aim is to simulate the mechanical property of the contact surface. Figure 5 shows the operation flow of the FRIC subroutine. The parameters of the improved hypoplastic model , hs, n, ed0, and so on, are inputted into the main program. After that, the contact state of soil-structure interface is estimated. If the interface is in contact, next stage (ABAQUS calculation kernel) is done, and if the interface is not in contact, it is returned to first stage. The stresses and displacements are calculated from the ABAQUS calculation kernel stage. In input FRIC stage, the reduced tensors (Section 2) are inputted to FRIC subroutine. In addition, the in-plane stress and void ratio e are inputted to this stage as additional state variables. In the next stage (Transformed to UMAT input), in order to access UMAT, input format of FRIC is formatted and transformed to UMAT subroutine. The stress is calculated and updated by Jacobian matrix. After the UMAT call, if the relative error “err” is greater than the tolerance “TOL,” it is gone to next stage (Newton–Raphson iteration), and if the relative error “err” is smaller than the tolerance “TOL,” it is gone to Update FRIC stage. A new value of must be estimated by Newton–Raphson iteration until err TOL. In Update FRIC stage, the stress tensors TUMAT are transformed back to TFRIC for each increment.

4. Results
4.1. Comparison with Mohr–Coulomb Model
The Mohr–Coulomb model and the improved model are compared using the condition of direct shear simulation model (Section 3.1). The parameters for the Hostun sand used in the simulation are defined in Table 1. The contact of the Mohr–Coulomb model is characterized only by the friction coefficient; it is decided to consider the roughness of the contact surface [36]. The two models were verified using normal stresses of 100 kPa and 200 kPa. The surface roughness was 46 . The relation between the shear stresses and shear displacements simulated by two models is shown in Figure 6. The default friction model is the Coulomb model of linear elasticity-ideal plasticity, so it cannot represent the hardening/softening properties of the contact surface, whereas the improved model is a typical nonlinear model; the results of the subroutine simulation are in good agreement with the actual condition of the contact surface and represent the effect of roughness and normal stress on the shear stress to the shear displacement. Consequently, it is proved that the new improved model is correct and feasible.

(a)

(b)
4.2. Comparison with the Explicitly Granular Hypoplastic Models
The model from Arnold and Herle [25], the model from Stutz et al. [27], and the improved model are based on the hypoplastic model from Wolffersdorff [14]. As described above (Section 2), the difference of three models lies in the interface conditions under the assumption of reduced stress and strain tensors. Therefore, the new improved model is compared to the model from Arnold and Herle [25] and the model from Stutz et al. [27].
In the following, the three models are simplified as(i)Model proposed by Arnold and Herle [25]: AH(ii)Model proposed by Stutz et al. [27]: STZ(iii)New improved model proposed in this paper: NEW
The three models use the same hypoplastic parameters, as shown in Table 1. The roughness of the contact surface for the three models was considered in the same way as that defined in Section 2.3. The applied normal stresses are 100 kPa and 200 kPa; the applied roughnesses of the contact surface are 2.4, 20.5, and 46.0 . Figure 7 shows the relation between the shear stresses and shear displacements of the three models according to the surface roughness with 100 kPa and 200 kPa applied normal stresses. The AH model gives the lowest shear stress at different roughnesses of the contact surface and different normal stress levels. The three models have the similar trend but the shear stresses obtained vary according to the roughness and the normal stress. In addition, in the cases of three models, the shear stresses increase as the surface roughnesses increase. This is the reason why the restraint effect of structural surface on soil in contact zone increases to raise the resistance of the sliding with the increase of surface roughness. This is also the reason that by increasing the surface roughness, the interface coefficient of friction increases so that the restriction on the soils by structures increases. As shown in Figures 7(a) – 7(f), in all cases, the difference between the improved model (NEW) and the STZ model is smaller compared to that between the new improved model and the AH model. This reason is the same as that described in Section 2: the improved model (NEW) and the STZ model have used the assumption of the transverse isotropy at the contact surface; however, the AH model has been proposed under the assumption of the isotropy at interface. As a result, the new improved model (NEW) and the STZ model show a small difference in their shear stress response.

(a)

(b)

(c)

(d)

(e)

(f)
4.3. Comparison with Experimental Results
The first verification is done by our experimental data. A direct shear experiment on the glass sand-concrete structure interface is conducted with an improved direct shear apparatus to analyze the effects of the roughness of the contact surface and the normal stress on the shearing strength of the interface. The properties of the glass sand are shown in Table 2. The used concrete plates have had relative surface roughnesses R of 0.1, 1, 1.5, and 2 mm. The shear rate is 0.8 mm/min, and the cross-sectional area of the sample is 60 × 60 mm. The applied normal stresses are 50, 100, and 200 kPa.
The parameters of the glass sand were decided from previous study [32]. The parameters are given in Table 3. The new improved model (NEW) and the AH and the STZ models were compared with our experimental data. As described in Section 2.3, the roughness of the contact surface was considered in the same way. Figure 8 shows the relation between the shear stresses and shear displacements for the three models and our experimental data according to the surface roughness with 50 kPa applied normal stress.

(a)

(b)

(c)

(d)
As shown in Figure 8, the simulation result of the new improved model (NEW) is the nearest to the experimental data in shear stress response than other two models. This reason is that the new improved model (NEW) represents the mechanical property of the contact surface by considering the plane strain state precisely. As a result, the simulation results agree well with the experimental data.
Figure 9 shows the shear stresses-shear displacements of the three models and our experimental data according to the normal stress with 1 mm applied surface roughness. The three models show the similar response in shear behavior at different applied normal stress levels. However, the predictions obtained from the improved model (NEW) were better than those obtained from the AH and STZ models. This reason is that the improved model (NEW) has been proposed by considering the plane strain condition precisely and more exactly simulates the mechanical property of the contact surface than other two models.

(a)

(b)
The last verification is done by using experimental data by Shahrour and Rezaie [37]. They used modified direct shear box to study the effect of surface roughness ( = 1.0/Rmax = 46.0 ) on the relation between shear stress and shear strain under the constant normal stress ( = 100 kPa). The sample was Hostun sand (e0 = 0.95) and the parameters are given in Table 1. The new improved model (NEW) and the AH and the STZ models were used to reproduce the experimental data; the results are illustrated in Figure 10. As shown in Figure 10, three models simulate the hardening properties of the contact surface for loose sand. In particular, the simulation result of the new improved model (NEW) agrees better with the experimental data than those of other two models. This reason is that the improved interface model (NEW) considers the plane strain condition at the contact surface to represent the real state more accurately than the previous models. Therefore, it is found that the new model has high reliability.

5. Summaries and Conclusions
The new improved model is incorporated by new reduced stress and strain rate tensors. The stress and strain rate tensors have been reduced by considering the plane strain condition accurately. Considering the relation between the relative roughness R and critical friction angle , the roughness of the contact surface has been also introduced into the new improved model.
The new model was compared to the Mohr–Coulomb model. The result shows a better prediction with the new model. The new improved model is a typical nonlinear model, so it represents the hardening/softening properties of the contact surface. The new model was compared to the model developed by Arnold and Herle [25] (AH) and the model proposed by Stutz et al. [27] (STZ). The difference between the new model and the STZ model is smaller compared to that between the new model and the AH model.
The simulation results of the new improved model are close to the experimental data. These were demonstrated by the comparison of our test and the experimental data from Shahrour and Rezaie [37].
The new improved model represents the effects of roughness and normal stress on the shear stress to the shear displacement. The new model has advantage to reduce calculation time. The limitation of the improved model is that interface problem under monotonic loading can only be modelled. In the future, the proposed model may be widely used to solve boundary value problems such as interface problem under cyclic loading.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors appreciate the financial support of science and technology project of PowerChina Road Bridge Group Co., Ltd., and the National Natural Science Foundation of China (grant no. 51178044).