Numerical Modeling of Soil Erosion with Three Wall Laws at the Soil-Water Interface

In the area of civil engineering and especially hydraulic structures, we find multiple anomalies that weakens mechanical characteristics of dikes, one of the most common anomalies is erosion phenomenon specifically pipe flow erosion which causes major damage to dam structures. This phenomenon is caused by a hole which is the result of the high pressure of water that facilitate the soil migration between the two sides of the dam. It becomes only a question of time until the diameter of the hole expands and causes destruction of the dam structure. This problem pushed physicist to perform many tests to quantify erosion kinetics, one of the most used tests to have logical and trusted results is the HET (hole erosion test). Meanwhile there is not much research regarding the models that govern these types of tests. Objectives: In this paper we modeled the HET using modeling software based on the Navier Stokes equations, this model tackles also the singularity of the interface structure/water using wall laws for a flow turbulence. Methods/Analysis: The studied soil in this paper is a clay soil, clay soil has the property of containing water more than most other soils. Three wall laws were applied on the soil / water interface to calculate the erosion rate in order to avoid the rupture of such a structure. The modlisitation was made on the ANSYS software. Findings: In this work, two-dimensional modeling was carried of the soil.in contrast of the early models which is one-dimensional model, the first one had shown that the wall-shear stress which is not uniform along the whole wall. Then using the linear erosion law to predict the non-uniform erosion along the whole length. The previous study found that the wall laws have a significant impact on the wall-shear stress, which affects the erosion interface in the fluid/soil, particularly at the hole's extremes. Our experiment revealed that the degraded profile is not uniform.


Introduction
Earthen hydraulic structures, which include dams, dikes or levees, are used worldwide to limit standing or flowing water in a given area. But the presence of water itself in such structures can cause significant damage that can lead to failure of the structure. This damage can be caused by three main mechanisms that are slope instability, overtopping and internal erosion. According to statistics from Foster et al. (2000), approximately 30-50% of dam failures are caused by internal soil erosion, which make internal erosion one of the major causes of dam's instability With global warming, which causes climatic disorder, the rise of lever waters, it is evident that natural disasters bound to extreme rainy events and afterward floods are more frequent and devastating. With climate disorder and rising lever waters caused by global warming, it is obvious that natural disasters related to extreme rainfall occurrences and subsequent floods are more frequent and devastating. Due to this, many dam failures have occurred around the world, some of them are reported by Foster et al. (2000) [1].
Soil erosion is a natural process caused by different natural events like wind, rain and moving water, excessive erosion is often caused by human activities, such as deforestation, road construction and intensive farming. It affects farming, construction projects, homeowners living near rivers, oceans and terrestrial slopes. The rate of soil erosion for a specific land area is determined from the loss of soil mass over time; you can calculate it by measuring the loss of soil mass over a specific time period. The vulnerability of geotechnical structures interface to fluids makes its maintenance an important area of research [2,3]. Understanding the underlying mechanisms and quantifying the effects of pertinent variables are both important. For sand erosion, its unconsolidated and weakly bonded particle assembly makes it highly vulnerable to flowing liquid. This can be modeled by different approaches, which are roughly divided into two types: continuum-based models and discrete models.
The common feature for all the continuum-based/empirical models is in the establishment of the conditions where sand failure will occur, and also in the use of some parameters which are calibrated with laboratory and/or field observations to predict the sand erosion rate. Various models for the prediction of sand erosion rate have been developed [2,4], the one developed by Bonelli et al. (2006Bonelli et al. ( , 2008 was for interpreting the whole erosion test with a constant pressure drop [4,5]. This type of model yielded a characteristic erosion time depending on the initial hydraulic gradient and the coefficient of erosion [6][7][8][9][10][11][12]; Another area of research is piping, which involves the formation and development of a continuous pipe within an earthen structure by the erosion of surrounding soil material driven by a hydraulic gradient. Failure of this type of structure can occur in one of three modes: piping through the fill, piping from the fill to the foundation and piping through the foundation.

Modelling Approach of the Interface Erosion
Turbulence models based on the hypothesis of high turbulence are not valid in the region close to the walls. There are then several solutions to take into account the presence of viscous effects: (1) Do not solve this region and use a wall law instead, (2) Introduce damping functions forcing the behavior of the model, henceforth called low-Reynolds model, or (3) Solve different equations for the main flow and the near-wall flow. The standard k-ε or RNG models are based on approach (1) [13][14][15][16].
Thus, in the first approach, the inner region is not resolved and a wall law is then used to connect the region affected by viscous effects to the fully turbulent part of the boundary layer (Figure 2.a). This more or less-evolved law makes it possible to estimate the average speeds as well as the turbulent quantities in the unresolved part but is nevertheless based on very constraining hypotheses. There are different levels of more or less evolved wall laws, the simplest one is the standard wall law (relation 3), valid in the region characterized by + = 30. The disadvantage of this law comes from the fact that it is based, among other things, on a very large Reynolds number of the flow, the absence of an adverse pressure gradient. To satisfy these points, other laws more elaborated are then available and built most of the time on a two-layer decomposition of the boundary layer (two-layer law). Depending on the calculation point placed in the boundary layer, the wall law will assign appropriate values. Thus, the use of this type of law allows finding information in the viscous under layer. The value of the meshing criterion must therefore satisfy + ≈ 1 [17]. It can therefore be noted that the wall law approach can be less consuming in terms of calculation time and storage space depending on the degree of complexity of the law, the size of the first element being more important [17][18][19][20][21].

Figure 2. Digital processing of the area near a wall. (a) Approach of wall laws; (b) Approach of low-Reynolds models
The value of the meshing criterion must therefore satisfy + ≈ 1. It can therefore be noted that the wall law approach can be less consuming in terms of calculation time and storage space depending on the degree of complexity of the law, the size of the first element being more important. Approaches (2) and (3) require in all cases much smaller mesh sizes, small enough to search for information in the viscous underlay ( Figure 2.b). The number of calculation elements then becomes higher than that of a calculation based on high-Reynolds models using standard wall laws [22][23][24][25][26][27].
With √ / homogeneous at a speed that is called speed of friction, noted ( ), and having no physical meaning (one prefers nevertheless to manipulate a speed instead of a surface tension). For average speeds, it is in fact better to focus on the / ratio, and not on U alone. Similarly, it is preferable to look at the ratio ′ ′ ̅̅̅̅̅̅̅ / 2 instead of ′ ′ ̅̅̅̅̅̅̅ alone. The friction speed is around 5% of the flow velocity. The transverse evolution of the longitudinal velocity in the inner region of the boundary layer shows several zones composing this inner region of the boundary layer: (1) the viscous sub-layer ( +< 5) in which the Reynolds tensor is negligible in front of the viscous stresses and in which the velocity evolves linearly with the distance to the wall: And Equation 2 the logarithmic zone (30 < +< 1000) dominated by turbulent stresses and in which a logarithmic evolution of speed is observed: Where; k and B being experimentally obtained model constants. The region between these two previous zones (5 < + < 30) is called the buffer zone. The inner region of the boundary layer therefore depends only on the parameters + and , and thus remain isolated from the main flow. Outside the logarithmic zone, the Equation 3 is no longer valid because the effects of the main flow are less and less negligible. The velocity relation which depended only on + must tend asymptotically towards a relation which this time is a function of / . We then enter a region characterized by a deficit velocity law and characterized by a relation of the type:  The CODE Fluent offers different types of wall treatment.

Standard Wall Function
It is an integration model connecting conditions at and near the wall based on the universal profiles of the turbulent boundary layer (parietal law and logarithmic law), these laws are based on the approach proposed by Launder and Spalding [8]: Where; * = The turbulent kinetic energy is calculated in the whole field. The wall condition imposed for k is then: / = 0. The calculation of k and ε, in a fluid cell adjacent to the wall, is performed with the assumption of local equilibrium between. The calculation of k and , in a fluid cell adjacent to the wall, is performed with the assumption of local equilibrium using scalable wall functions that ensure that the distance to the wall employed in the wall functions is such that y + ≥11.126 regardless of the level of grid refinement near the wall. This value of y + marks the intersection of the linear and logarithmic velocity profiles. It is clear that the evolving wall function will produce identical results to the standard wall function for values of y + ≥11.126.
The use of the standard equilibrium and non-equilibrium wall functions implies that the near-wall grid node lies entirely in the logarithmic region, but this can be difficult to achieve in overall applications given the varying geometric and velocity scales associated with the flow. Scalable wall functions can prevent erroneous modeling of the laminar and buffer regions of the boundary layer by effectively moving the near-wall mesh point to y + =11.126, regardless of the proximity to the wall. It is important to know that these functions can be activated for all high Reynolds number turbulence models that use equilibrium wall functions without surface roughness. They are not currently implemented for use with non-equilibrium wall functions or with the Reynolds stress transport model. Where; : is the kinetic energy calculated in the considered cell.

Enhanced Wall Treatment
In the case of improved walls, a single speed profile is used. This profile covers the entire boundary layer (inner and outer parts), in which case Fluent uses the following relationship Where is the transition function given by: With = 0.01 and = 5, and composed of a part; Laminar, + with: Or the α coefficient takes into account the parietal pressure gradient, As well as a turbulent part, + with: In this part also, the pressure gradient is represented by, through the coefficient ̀: Or we have: In the case where the above coefficients, α, þ and y are all equal to zero, the relation (3) admits as solution a classical logarithmic wall law [17][18][19][20][21].

Scalable Wall Function
Irrespective of the near-wall grid refinement's level, the wall distance applicated to the wall function which is + ≥ 11.126 was ensured by the scalability of the wall functions. Which marks the intersection of the velocity profile's logarithmic and linear. It's clear that a scalable wall function generates identical results to the standard wall function at the value of + ≥ 11.126.
Standard equilibrium and non-equilibrium wall function imply that the log-law region must entirely include the nearwall grid. It is difficult to achieve with the general application due to the velocity scales and varying geometrical with the flow. The erroneous modeling of the laminar and buffer regions of the boundary layer was be avoided by scalable wall function near to mesh point + = 11.126 no matter his nearest to the wall [9].

Analysis and Interpretation of Results
The fluid domain is assumed axisymmetric, which has 117 mm length in the axial direction (z) and 3 mm large in the radial direction (r). The inlet section and the outlet section are respectively in the left and the right side as the Figure 4 shown:     Figure 6 presented the value of y + deferent pressure used = 6000 , gives the three wall laws too as axial coordinate function. Figure 7 gives the function of the axial coordinate of the wall-shear stress for the three wall law's hydraulic gradients. Figure 7 shows the uniforms of the wall-shear along the hole length. The wall-shear stress can exceed 7 times and it has a permanent value inside the hole's plateau zone. Which is decline the habitual hypothesis that is used to derive the onedimensional approaches of piping flow. The wall-shear and erosion are maximal in the upstream extremity. The twodimension's results obtained of flow state clearly showing huge variations compared with the one-dimension's results. The initial stage can be predicted just in case of rigidity of the wall, and linearity of the erosions law was valid. So, the strong effect on the erosion rate was observed by the pressure gradient's application.     11 shows the erosion rate (of 10 -6 kg/s), corresponding to the amount of mass production per unit time because of the erosion. This amount is a result of the Integration of the erosion law along the whole length of the hole and the multiplication of the results and the initial circumference of the hole. The linear classical erosion law: ̇= ( − ) was used, with the erosion constants = 5.5 10 −4 / and τ = 7 as identified for a soil sample [10,11]. The results also indicate that the SWF (Scalable Wall Function) law is not compatible to improve the near wall calculation due to poor estimates.

Figure 11. Erosion Rate for 6000Pa by wall laws
By understanding and analyzing Figure 11, we notice that it gives us high shear stress at the entrance. This makes us conclude that the erosion of entrance in experiments is due to the excessive efforts of fluid. A process of the successive armoring may be reflected by increasing the initial period and decreasing the erosion rats as well. Even if the looser and smaller particles are directly deleted from the adjacent area to the axial hole, the course's soil matrix reset intact. Which making the surface rougher and more erosion-resistant over time. While the stress shears up to the critical value, which is able to remove and destroy the developed armoring, because of the excessive force's application to the soil underneath the specimen stat to speared along the whole hole.

Conclusion
Experimental studies were used to study the internal erosion wall of the soil and its effect on the hydraulic structure, Validation of those experimental results has no enough data at the numerical side which is far from our process. The main use of the Hole Erosion Test (HET) is to determine characteristics of the soil erosion, (erosion rate coefficient. In this study, we developed a new numerical model for modeling internal erosion of soil, in our model we determine the shear stress for three different pressures in order to compare the effects of three wall laws at the level of the solid/fluid interface. Using an interpretation of the hole erosion test with three constant pressure drops we obtained an analytical scaling law in the interface. The main advantage of this simulation is the ability to determine the shear stress in the wall, even if we have à turbulent flow inside the hole, which is not the case for the experimental study where we neglect the turbulent flow and consider only laminar flow. In this work, two-dimensional modeling was carried of the soil.in contrast of the early models which is one-dimensional model, the first one had shown that the wall-shear stress which is not uniform along the whole wall. Then using the linear erosion law to predict the non-uniform erosion along the whole length. The previous studies show that the consequence of the wall laws has not a negligible effect on the wall-shear stress, which affects the erosion's interface in the fluid/soil, thus, particularly at the extremity of whole the hole where is maximal. The eroded profile is not uniform, this is the result of our experiment.  Erosion Rate SWF-6000Pa Erosion Rate ScWF-6000Pa Erosion Rate EWT-6000Pa

Data Availability Statement
The data presented in this study are available in article.