Influence of Soil Heterogeneity on the Contact Problems in Geotechnical Engineering

Contact problems are widely encountered in geotechnical engineering, such as the contact between soils and concrete used in earth and rockfill dams, tunnels and coastal levees. Due to the unknown contact region and contact forces, the contact problems have strong boundary nonlinearity. In addition, soils have been recognized as heterogeneous materials in geotechnical engineering. The existence of the soil heterogeneity increases the nonlinearity of the contact problems. Currently, the contact problems are mostly analysed without considering the soil heterogeneity, which may not reflect the contact behavior well. In order to investigate the influence of soil heterogeneity on the contact problems, in this paper, a simple plane-strain contact problem is analysed as an example. In this example, Young’s modulus is taken to be a spatially variable. The local average subdivision (LAS) is used to model the heterogeneity of Young’s modulus. The penalty method is utilised to determine the contact behavior. By the first use of linking the penalty method with the LAS, the proposed approach can be used to analyse the contact problems considering soil heterogeneity. The results show that the influence of soil heterogeneity on the elastic contact problems is significant. The contact forces of the heterogeneous case present apparent variation compared to the results of the homogeneous case. The distribution of the contact force at a specific point is also normal when Young’s modulus is normally distributed, moreover, the coefficient of variation (COV) and the horizontal scale of fluctuation of Young’s modulus affect the extent of variation of the normal contact forces. The standard deviation of the normal contact force increases with the increase of the COV and decreases with the increase of the horizontal scale of fluctuation of Young’s modulus. From the analyses, to better predict the deformation/stress in the contact problems, heterogeneity needs to be considered.


Introduction
In geotechnical engineering, contact problems are commonly seen, such as the contact between soils and concrete used in earth and rockfill dams [1,2], tunnels [3] and coastal levees. In order to accurately predict the performance of these engineering projects, the contact between soils and concrete needs to be analysed. For example, the contact force applied on the concrete lining of tunnels is important for checking whether the strength of the concrete is qualified or not. If it is not qualified, the breakage of the concrete could result in water leakage and even subsequent failure of these projects. However, due to the unknown contact region and contact forces, the contact problems have strong boundary nonlinearity, which makes the problem difficult.
Currently, the frequently used method in solving contact problems contains the Goodman element method [4], the thin-layer element method [5], the Lagrange multiplier method [6][7][8][9] and the penalty method [10,11]. The Goodman element method is the most widely used method to simulate the interface contact between two materials with different modulus. In this method, the thickness of the Goodman element is zero [2,12]. In contrast, in the thin-layer element method, the interface is replaced by an equivalent solid element with a small thickness. In addition, the element is given a relatively low modulus and can experience large deformation [13]. The Lagrange multiplier method introduces additional variables (the Lagrange multiplier), e.g., the contact force, to enforce directly and exactly the contact constraints. The penalty method enforces the contact constraints by applying virtual springs with zero rest lengths to surfaces in contact. A point is penalized for penetrating another project. Qian et al. [14] compared the Goodman element method, the thin-layer method and the penalty method in a numerical simulation of a concrete-faced rockfill dam. The comparison results indicated that the penalty method can reflect the contact stress along the interface between the soil and the concrete slab well, especially when there is a gap which emerges between the soil and the concrete slab. In addition, compared to the Lagrange multiplier, the penalty method is easy to implement and there are no additional efforts required to solve the multiplier [15]. Therefore, it is worth investigating the contact behavior by using the penalty method.
In the previous research of contact problems, parameters are usually considered to be homogeneous. However, heterogeneity has proven to be influential in many fields, such as Zeng et al. [16], Amato et al. [17], etc. In geotechnical engineering, soils have already been recognized as heterogeneous materials. The effect of soil heterogeneity has been investigated in many geotechnical problems, such as bearing capacity [18,19], seepage [20][21][22] and slope stability/displacement [23,24]. In Fenton and Griffiths [18] and Cho and Park [19], the cohesion and friction angle are considered as random variables. The results showed that the bearing capacity of heterogeneous foundations presented significant variation compared to the homogeneous foundations. The variation was strongly related to the mean, standard deviation and the correlation length of the soil properties. Griffiths and Fenton [20] investigated the seepage beneath a water retaining structure. In the study, the hydraulic conductivity was taken to be a random variable. The heterogeneity of hydraulic conductivity influenced the exit gradient, the flow rate and the uplift force significantly. Liu et al. [21] demonstrated the influence of the heterogeneity in hydraulic conductivity on the pore pressure distribution and the subsequent slope stability of an embankment. Yetbarek et al. [22] numerically examined the effects of the heterogeneity in hydraulic parameters on subsurface water movement which accounts for root water uptake. It was observed that the maximum difference in soil moisture values could reach 30% near the surface soil layers, due to an increase in soil heterogeneity on irrigation days. Hicks and Li [23] studied the effect of the heterogeneity in undrained shear strength on a 3D slope stability. ElGhoraiby and Manzari [24] evaluated the combined effect of the spatial variability in the soil density and base motion on the lateral spread in the numerical simulation of centrifuge modeling. All this research has shown that the influence of soil heterogeneity should be considered.
Therefore, in order to investigate the effect of soil heterogeneity on the elastic contact problems, a simple contact problem of an elastic body on a rigid foundation is analysed in this paper and the penalty method is taken to solve the contact problem. In addition, Young's modulus of the elastic body is considered to be heterogenous. The local average subdivision (LAS) is used to simulate the heterogeneity of Young's modulus. Compared to other methods of generating random field, such as the turning bands method and the fast Fourier transform, the LAS was shown to be more efficient and reliable [25]. Moreover, it is ideally suited for use with the finite element method (FEM). By linking the penalty method with the LAS, the disadvantage of the conventional analysis of the contact problems mentioned above can be overcome. In this paper, firstly, the penalty method and the LAS are briefly introduced. Then, an illustrative example is presented. The results of the homogeneous and heterogeneous cases are compared with each other. Finally, the influence of the statistics of Young's modulus, i.e., the standard deviation and the scale of fluctuation, on the contact force is investigated.

Implementation of the Penalty Method in the FEM
In this paper, the commonly used method, i.e., the penalty method, is utilised to solve the contact problems. In the penalty method, a small amount of penetration is allowed, and the contact force is proportional to the amount of penetration. In order to combine the penalty method with the FEM, the slave-master concept is introduced in the FEM.
In the slave-master concept, one of the two contact bodies is designated as a slave body and the other one is a master body, as shown in Figure 1. In FEM, one node from the slave body and two nodes from the master body are defined as one contact pair. In one contact pair, the first objective is to check whether the slave node X c is in contact with the master segment X 1 X 2 , i.e., the gap g n is positive (separated) or negative (in contact). If it is in contact, the contact forces at the slave node and the master nodes will be calculated by the penalty method. In addition, the contact stiffness at the three nodes will also be calculated.
influence of the statistics of Young's modulus, i.e., the standard deviation and the scale of fluctuation, on the contact force is investigated.

Implementation of the Penalty Method in the FEM
In this paper, the commonly used method, i.e., the penalty method, is utilised to solve the contact problems. In the penalty method, a small amount of penetration is allowed, and the contact force is proportional to the amount of penetration. In order to combine the penalty method with the FEM, the slave-master concept is introduced in the FEM.
In the slave-master concept, one of the two contact bodies is designated as a slave body and the other one is a master body, as shown in Figure 1. In FEM, one node from the slave body and two nodes from the master body are defined as one contact pair. In one contact pair, the first objective is to check whether the slave node Xc is in contact with the master segment X1X2, i.e., the gap gn is positive (separated) or negative (in contact). If it is in contact, the contact forces at the slave node and the master nodes will be calculated by the penalty method. In addition, the contact stiffness at the three nodes will also be calculated. After the contact search is finished, the stiffness and forces calculated by the penalty method will be added into the matrix equation as [26]: where, Kc is the contact stiffness matrix in the global coordinate, Fc is the contact force in the global coordinate, while K is the global stiffness matrix in the FEM and Fres is the residual force in the FEM generated from boundary conditions. In each contact pair (one slave node with one master segment), the contact stiffness matrix kc can be calculated as: in which, kN and kT are the stiffness matrices due to normal contact and consideration of friction, respectively.
For the calculation of kT, there are two different situations: the stick situation and the slip situation. For the stick condition, After the contact search is finished, the stiffness and forces calculated by the penalty method will be added into the matrix equation as [26]: where, K c is the contact stiffness matrix in the global coordinate, F c is the contact force in the global coordinate, while K is the global stiffness matrix in the FEM and F res is the residual force in the FEM generated from boundary conditions. In each contact pair (one slave node with one master segment), the contact stiffness matrix k c can be calculated as: in which, k N and k T are the stiffness matrices due to normal contact and consideration of friction, respectively.
For the calculation of k T , there are two different situations: the stick situation and the slip situation. For the stick condition, in which, ξ is the local coordinate of the projection point of the slave node on the master surface (see Figure 1). The superscript ' 0 ' denotes the previous time step. l 0 is the length of the master segment. For the slip condition, In addition, the contact force of each contact pair for the stick condition can be calculated as: f c = −ω n g n C n − ω t g t C t (7) The contact force of each contact pair for the slip condition can be calculated as: f c = −ω n g n C n − µω n g n sgn(g t )C t (8) In the above equations, the vectors C n , C t , P, Q are defined as follows: where, ω n and ω t are the penalty parameters, e n is the unit normal vector and e t is the unit tangential vector. The detailed derivation of K c and F c is referred to in Chapter 5 in Kim [26].

Radom Field Generation
In this paper, the LAS method [27] is undertaken to model the heterogeneity of soil parameters. The implementation of the LAS requires the information of soil parameters, which is the distribution of parameters (i.e., the mean and standard deviation). The LAS can generate standard normal random fields first, then these fields can be transformed to the required fields, which follows the input distribution.
In addition to the distribution of the soil parameters, in order to model the correlation of parameter values of two points, the scales of fluctuation in the vertical and horizontal directions are also indispensable. The scales of fluctuation are used to build the correlation function used in the LAS. In this paper, an exponential Markov correlation function is taken as follows: where, ρ(τ 1 , τ 2 ) is the correlation coefficient between the values of two points, τ v and τ h are the lag distance between two points in the vertical and horizontal directions in a random field, respectively. θ v and θ h are the vertical and horizontal scales of fluctuation, respectively. The reader refers to Samy [28] and Spencer [29] for the detailed implementation of the LAS in FEM.

Illustrative Example
In this paper, a simple plane-strain contact problem is analysed as an example. Figure 2 shows the flowchart of using the proposed method to analyse the illustrative example. The LAS is first used to generate multiple realisations of the random field, and then each realisation is imported to the contact analysis using the penalty method. Finally, all the results from all the realisations are collected and analysed. In this example, an elastic block is pressed against a rigid foundation, which is extracted from the contact between the overlying soil and the concrete lining in tunnels. The elastic block is 4 m long and 2 m high. On the top of the elastic block, there is a uniformly distributed load. The uniformly distributed load is 200 kN/m. The finite element mesh for this problem is shown in Figure 3, and the size of the mesh is 0.2 m by 0.2 m. Each element has 4 nodes. In this example, the Poisson's ratio ν is 0.3. Young's modulus is taken to be a random variable, which follows a normal distribution. overlying soil and the concrete lining in tunnels. The elastic block is 4 m long and 2 m high. On the top of the elastic block, there is a uniformly distributed load. The uniformly distributed load is 200 kN/m. The finite element mesh for this problem is shown in Figure  3, and the size of the mesh is 0.2 m by 0.2 m. Each element has 4 nodes. In this example, the Poisson's ratio ν is 0.3. Young's modulus is taken to be a random variable, which follows a normal distribution.

Validation of the Numerical Results
In order to validate the program used in this paper, the results computed by the program are compared with the results obtained from the widely used commercial software Abaqus. In this comparison, two scenarios have been selected: in the first scenario, there is pressed against a rigid foundation, which is extracted from the contact between the overlying soil and the concrete lining in tunnels. The elastic block is 4 m long and 2 m high. On the top of the elastic block, there is a uniformly distributed load. The uniformly distributed load is 200 kN/m. The finite element mesh for this problem is shown in Figure  3, and the size of the mesh is 0.2 m by 0.2 m. Each element has 4 nodes. In this example, the Poisson's ratio ν is 0.3. Young's modulus is taken to be a random variable, which follows a normal distribution.

Validation of the Numerical Results
In order to validate the program used in this paper, the results computed by the program are compared with the results obtained from the widely used commercial software Abaqus. In this comparison, two scenarios have been selected: in the first scenario, there

Validation of the Numerical Results
In order to validate the program used in this paper, the results computed by the program are compared with the results obtained from the widely used commercial software Abaqus. In this comparison, two scenarios have been selected: in the first scenario, there is no friction along the interface, and in the second one the coefficient of friction of the interface is 0.3. In both cases, Young's modulus is assumed to be uniform across the whole research domain and equal to the mean of Young's modulus. interface is 0.3. In both cases, Young's modulus is assumed to be uniform across the whole research domain and equal to the mean of Young's modulus. Figure 4 shows the comparison between the results of the program and the results of Abaqus. The light blue line with the blank circle and the grey line with the blank triangle show the results of the normal contact force for the first scenario. The light red line with blank square and the orange line with cross show the normal contact force for the second scenario. The blue line with the solid circle and the green line with solid diamond show the frictional contact force for the second scenario. From Figure 4, the results computed by the program are close to the results of Abaqus, which proves the validity of the program.

Comparison between the Homogeneous Case and Heterogeneous Case
In order to investigate the influence of heterogeneity on the contact problems, a homogeneous case is analysed first. As in the previous subsection, Figure 4 presents the results of the homogeneous case, which are the computed contact force applied on the bottom of the elastic block along the interface when the coefficient of friction is 0.0 and 0.3. From Figure 4, it can be seen the normal contact force is symmetrical when the Young's modulus is uniform. In addition, Figure 5 shows the deformation of the elastic block under loading for both scenarios.

Comparison between the Homogeneous Case and Heterogeneous Case
In order to investigate the influence of heterogeneity on the contact problems, a homogeneous case is analysed first. As in the previous subsection, Figure 4 presents the results of the homogeneous case, which are the computed contact force applied on the bottom of the elastic block along the interface when the coefficient of friction is 0.0 and 0.3. From Figure 4, it can be seen the normal contact force is symmetrical when the Young's modulus is uniform. In addition, Figure 5 shows the deformation of the elastic block under loading for both scenarios.    Figure 5a,b, it can be found that because of the friction, the shape of the zone with large horizontal displacement becomes significantly smaller. The horizontal displacement reduces much near the interface. Moreover, Figure 5c,d indicate that the shape of the zone with large vertical displacement changes as well. In addition, the maximum horizontal and vertical displacements of the case with no friction are larger than those of the case when the coefficient of friction is 0.3.   Figure 5 that the horizontal/vertical displacement is symmetrical. However, due to the emergence of friction, the distribution of the displacement across the whole domain is quite different. From Figure 5a,b, it can be found that because of the friction, the shape of the zone with large horizontal displacement becomes significantly smaller. The horizontal displacement reduces much near the interface. Moreover, Figure 5c,d indicate that the shape of the zone with large vertical displacement changes as well. In addition, the maximum horizontal and vertical displacements of the case with no friction are larger than those of the case when the coefficient of friction is 0.3.
Then, Young's modulus of the elastic block is assumed to be spatially variable. The mean and coefficient of variation (COV) (defined as the standard deviation divided by the mean) of Young's modulus are 100 MPa and 5%, respectively. Additionally, the scale of fluctuation of Young's modulus is chosen to be 1.0 m. Figure 6 shows one realisation of the random field. From Figure 6, each finite element mesh has been assigned with a different value. In this paper, the total number of realisations used in the computation is 500. In each realisation, the contact forces can be calculated.  that the horizontal/vertical displacement is symmetrical. However, due to the emergence of friction, the distribution of the displacement across the whole domain is quite different. From Figure 5a,b, it can be found that because of the friction, the shape of the zone with large horizontal displacement becomes significantly smaller. The horizontal displacement reduces much near the interface. Moreover, Figure 5c,d indicate that the shape of the zone with large vertical displacement changes as well. In addition, the maximum horizontal and vertical displacements of the case with no friction are larger than those of the case when the coefficient of friction is 0.3.
Then, Young's modulus of the elastic block is assumed to be spatially variable. The mean and coefficient of variation (COV) (defined as the standard deviation divided by the mean) of Young's modulus are 100 MPa and 5%, respectively. Additionally, the scale of fluctuation of Young's modulus is chosen to be 1.0 m. Figure 6 shows one realisation of the random field. From Figure 6, each finite element mesh has been assigned with a different value. In this paper, the total number of realisations used in the computation is 500. In each realisation, the contact forces can be calculated.   Figures 7 and 8, it can be found that the heterogeneity of Young's modulus has a significant influence on the results of contact force. First, due to the spatial variability of Young's modulus, the distribution of normal contact force along the contact surface is not symmetrical. Moreover, Figures 6 and 7 show that the COV of the contact force at a specific point changes with the x coordinate. Furthermore, it can be seen from Figure 8 that the variation of the frictional contact force is less than that of the normal contact force, which suggests that the influence of heterogeneity has more of an impact on the normal contact force. be found that the heterogeneity of Young's modulus has a significant influence on the results of contact force. First, due to the spatial variability of Young's modulus, the distribution of normal contact force along the contact surface is not symmetrical. Moreover, Figs. 6 and 7 show that the COV of the contact force at a specific point changes with the x coordinate. Furthermore, it can be seen from Figure 8 that the variation of the frictional contact force is less than that of the normal contact force, which suggests that the influence of heterogeneity has more of an impact on the normal contact force.    Later, to investigate the detailed variation of the contact forces at each point, the normal contact force at the middle bottom point A (shown in Figure 3) of the elastic block has been selected as an example. The following Figures 9-11 show the normal and frictional contact force at the point A of the elastic block from 500 realisations when the coefficient of friction is 0 and 0.3, respectively. In the homogeneous case, the normal contact forces at the point A for two scenarios (one has no friction, and the coefficient of friction is 0.3 in the other one) are 40,099.7 N and 35,283.8 N. From Figures 9a, 10a and 11a, it can be found that the distribution of the contact forces follows a normal distribution. When there is no friction, the mean and the standard deviation of the normal contact force are 40,068 N and 835.6 N. When the coefficient of friction is 0.3, the mean and the standard deviation of the normal contact force are 35,279 N and 677.2 N. Meanwhile, the mean and the standard deviation of the frictional contact force are −2.0746 N and 182.1 N. The contact force of the homogeneous case is quite close to the mean of the heterogeneous case. By comparing Figure 9a and 10a, it can be seen that the contact with no friction generally has larger normal contact forces, which is consistent with the homogeneous case. In addition, from Figures 9b and 10b, the probability of the contact force exceeding a certain value can be calculated. For example, in the scenario without friction, if the normal contact force of the homogeneous case (40,099.7 N) is taken as a limit, the probability of exceeding the limit is 1 − P (normal contact force < 40,099.7 N), i.e., 48.2% in this paper. Similarly, when the coefficient of friction is 0.3, if the normal contact force of the homogeneous case (35,283.8 N) is taken as a limit, the probability of exceeding the limit is 49.8%. When the coefficient of friction is 0.3, if the frictional contact force of the homogeneous case (0 N) is taken as a limit, the probability of exceeding the limit is 51%.
force of the homogeneous case (40,099.7 N) is taken as a limit, the probability of exceeding the limit is 1 − P (normal contact force < 40,099.7 N), i.e., 48.2% in this paper. Similarly, when the coefficient of friction is 0.3, if the normal contact force of the homogeneous case (35,283.8 N) is taken as a limit, the probability of exceeding the limit is 49.8%. When the coefficient of friction is 0.3, if the frictional contact force of the homogeneous case (0 N) is taken as a limit, the probability of exceeding the limit is 51%.   Since the influence of heterogeneity on the contact force is significant, the heterogeneity certainly has an influence on the corresponding displacement of the elastic block. Figure 12 shows the displacement of the elastic block from two realisations when there is no friction. Comparing Figure 12 with Figure 5a,c, it can be seen that the displacement of the heterogeneous case is significantly different from that of the homogeneous case. For the horizontal and vertical displacements, the variation of the horizontal displacement is much more obvious. Figure 13 shows the displacement of the elastic block from two realisations when the coefficient of friction is 0.3. From the figure, the variation of the horizontal and vertical displacements is not significant compared to the scenario with no friction. Since the influence of heterogeneity on the contact force is significant, the heterogeneity certainly has an influence on the corresponding displacement of the elastic block. Figure 12 shows the displacement of the elastic block from two realisations when there is no friction. Comparing Figure 12 with Figure 5a,c, it can be seen that the displacement of the heterogeneous case is significantly different from that of the homogeneous case. For the horizontal and vertical displacements, the variation of the horizontal displacement is much more obvious. Figure 13 shows the displacement of the elastic block from two realisations when the coefficient of friction is 0.3. From the figure, the variation of the horizontal and vertical displacements is not significant compared to the scenario with no friction.  Since the influence of heterogeneity on the contact force is significant, the heterogeneity certainly has an influence on the corresponding displacement of the elastic block. Figure 12 shows the displacement of the elastic block from two realisations when there is no friction. Comparing Figure 12 with Figure 5a,c, it can be seen that the displacement of the heterogeneous case is significantly different from that of the homogeneous case. For the horizontal and vertical displacements, the variation of the horizontal displacement is much more obvious. Figure 13 shows the displacement of the elastic block from two realisations when the coefficient of friction is 0.3. From the figure, the variation of the horizontal and vertical displacements is not significant compared to the scenario with no friction.

Influence of the Statistics of Young's Modulus on the Contact Force
The previous subsection has demonstrated the influence of the heterogeneity on the results of contact force and displacement. The random field generation is controlled by the statistics of Young's modulus, i.e., the mean, the COV and the scale of fluctuation. Therefore, in this subsection, the purpose is to investigate the influence of different standard deviation and different scale of fluctuation.
First, the influence of different standard deviation is investigated. The COV of Young's modulus is chosen to be 5% (which has been used in the previous subsection), 10% and 20%, respectively. In the investigation, the coefficient of friction remains 0.3 and the scale of fluctuation is assumed to be 1.0 m. Figure 15 shows the PDFs and CDFs of the normal contact force at the point A when the COV is 5%, 10% and 20%. From Figure 15, it proves that with the increase of the COV of Young's modulus, the mean of the normal contact force does not have significant variation, whereas the standard deviation of the normal contact force increases. It is because that with the increase of the COV, the range of the value of Young's modulus becomes wider. When the COV is 5%, the mean and the standard deviation of the normal contact force are 35,279 N and 677.2 N. When the COV is 10%, the mean and the standard deviation of the normal contact force are 35,266 N and 1371.3 N. When the COV is 20%, the mean and the standard deviation of the normal contact force are 35,255 N and 2878.4 N.

Influence of the Statistics of Young's Modulus on the Contact Force
The previous subsection has demonstrated the influence of the heterogeneity on the results of contact force and displacement. The random field generation is controlled by the statistics of Young's modulus, i.e., the mean, the COV and the scale of fluctuation. Therefore, in this subsection, the purpose is to investigate the influence of different standard deviation and different scale of fluctuation.
First, the influence of different standard deviation is investigated. The COV of Young's modulus is chosen to be 5% (which has been used in the previous subsection), 10% and 20%, respectively. In the investigation, the coefficient of friction remains 0.3 and the scale of fluctuation is assumed to be 1.0 m. Figure 15 shows the PDFs and CDFs of the normal contact force at the point A when the COV is 5%, 10% and 20%. From Figure 15, it proves that with the increase of the COV of Young's modulus, the mean of the normal contact force does not have significant variation, whereas the standard deviation of the normal contact force increases. It is because that with the increase of the COV, the range of the value of Young's modulus becomes wider. When the COV is 5%, the mean and the standard deviation of the normal contact force are 35,279 N and 677.2 N. When the COV is 10%, the mean and the standard deviation of the normal contact force are 35,266 N and 1371.3 N. When the COV is 20%, the mean and the standard deviation of the normal contact force are 35,255 N and 2878.4 N.

Influence of the Statistics of Young's Modulus on the Contact Force
The previous subsection has demonstrated the influence of the heterogeneity on the results of contact force and displacement. The random field generation is controlled by the statistics of Young's modulus, i.e., the mean, the COV and the scale of fluctuation. Therefore, in this subsection, the purpose is to investigate the influence of different standard deviation and different scale of fluctuation.
First, the influence of different standard deviation is investigated. The COV of Young's modulus is chosen to be 5% (which has been used in the previous subsection), 10% and 20%, respectively. In the investigation, the coefficient of friction remains 0.3 and the scale of fluctuation is assumed to be 1.0 m. Figure 15 shows the PDFs and CDFs of the normal contact force at the point A when the COV is 5%, 10% and 20%. From Figure 15, it proves that with the increase of the COV of Young's modulus, the mean of the normal contact force does not have significant variation, whereas the standard deviation of the normal contact force increases. It is because that with the increase of the COV, the range of the value of Young's modulus becomes wider. When the COV is 5%, the mean and the standard deviation of the normal contact force are 35,279 N and 677.2 N. When the COV is 10%, the mean and the standard deviation of the normal contact force are 35,266 N and 1371.3 N. When the COV is 20%, the mean and the standard deviation of the normal contact force are 35,255 N and 2878.4 N.    Table 1 has shown the mean and standard deviation of the normal contact force for different horizontal scale of fluctuation. The reason could be that with the increase of the horizontal scale of fluctuation, the value of Young's modulus along the horizontal direction varies less, which can be seen in Figure 16c, so less variation causes decreases in the standard deviation.  Table 1 has shown the mean and standard deviation of the normal contact force for different horizontal scale of fluctuation. The reason could be that with the increase of the horizontal scale of fluctuation, the value of Young's modulus along the horizontal direction varies less, which can be seen in Figure 16c, so less variation causes decreases in the standard deviation.

Conclusions
In this paper, a simple elastic problem has been studied to investigate the influence of heterogeneity on the contact forces and deformation. It is the first time that the linking the penalty method has been used with the LAS. In the example, the parameter Young's modulus has been considered as a spatial variable. The penalty method and the LAS have been linked together to solve the contact problem and model the heterogeneity of Young's modulus.
First, in the homogeneous case in which Young's modulus is taken to be uniform across the whole research domain, it can be found that the normal/frictional contact force and the deformation is symmetrical. Compared to the homogeneous case, the results of the case where Young's modulus is considered to be spatially variable (following a normal distribution) shows that the influence of soil heterogeneity on the elastic contact problems is significant. Due to heterogeneity, the contact force and deformation shows evident variation. The variation of the normal contact force at any specific point follows a normal distribution, and the COV changes with the x coordinate. When the soils are considered to be heterogenous, it can be found that there is a certain probability when the contact force exceeds a certain value (result calculated in the homogeneous case). This implies that the conventional analysis of the contact problems in which the soils are taken as homogeneous could induce risk. In addition, the variation of the frictional contact force is less than that of the normal contact force.
Later, the influences of different COVs and horizontal scales of fluctuation have been investigated. The results show that with the increase of the COV of Young's modulus, the standard deviation of the normal contact force increases. However, with the increase of the horizontal scale of fluctuation of Young's modulus, the standard deviation of the normal contact force decreases.

Conclusions
In this paper, a simple elastic problem has been studied to investigate the influence of heterogeneity on the contact forces and deformation. It is the first time that the linking the penalty method has been used with the LAS. In the example, the parameter Young's modulus has been considered as a spatial variable. The penalty method and the LAS have been linked together to solve the contact problem and model the heterogeneity of Young's modulus.
First, in the homogeneous case in which Young's modulus is taken to be uniform across the whole research domain, it can be found that the normal/frictional contact force and the deformation is symmetrical. Compared to the homogeneous case, the results of the case where Young's modulus is considered to be spatially variable (following a normal distribution) shows that the influence of soil heterogeneity on the elastic contact problems is significant. Due to heterogeneity, the contact force and deformation shows evident variation. The variation of the normal contact force at any specific point follows a normal distribution, and the COV changes with the x coordinate. When the soils are considered to be heterogenous, it can be found that there is a certain probability when the contact force exceeds a certain value (result calculated in the homogeneous case). This implies that the conventional analysis of the contact problems in which the soils are taken as homogeneous could induce risk. In addition, the variation of the frictional contact force is less than that of the normal contact force.
Later, the influences of different COVs and horizontal scales of fluctuation have been investigated. The results show that with the increase of the COV of Young's modulus, the standard deviation of the normal contact force increases. However, with the increase of the horizontal scale of fluctuation of Young's modulus, the standard deviation of the normal contact force decreases.
From the analysed results, it suggests that, in order to better predict the deformation/stress in the contact bodies (such as soil-concrete contact), the spatial variability needs to be considered. To achieve this, the proposed approach of linking the penalty method with the LAS can be applied.

Conflicts of Interest:
The authors declare no conflict of interest.

CDF
cumulative distribution function COV coefficient of variation FEM finite element method LAS local average subdivision PDF probability density function