Analytical solution of coupled soil erosion and consolidation equations by asymptotic expansion approach

Applied Mathematical Modelling xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

Analytical solution of coupled soil erosion and consolidation equations by asymptotic expansion approach M.A. Parron Vera a, F. Yakhlef a, M.D. Rubio Cintas a, O. Castillo Lopez a, P. Dubujet b, A. Khamlichi c,⇑, M. Bezzazi c a

Civil and Industrial Engineering Department, High Polytechnic School of Algeciras, University of Cadiz, Ramon Puyol Avenue, Algeciras 11202, Spain Tribology and Dynamics of Systems Laboratory, National School of Engineers at Saint-Etienne, 58 Jean Parot Avenue, Saint-Etienne 42100, France c Modelling and Analysis of Systems Laboratory, Faculty of Sciences at Tetouan, BP. 2121 M’hannech, Tetouan 93002, Morocco b

a r t i c l e

i n f o

Article history: Received 21 July 2011 Received in revised form 11 January 2014 Accepted 3 February 2014 Available online xxxx Keywords: Porous media Internal erosion Consolidation Damage Asymptotic expansion Finite differences

a b s t r a c t In this work, one-dimensional approximation of internal erosion taking place in a soil made from sand and clay mixture was considered. The clay phase that is susceptible to experience erosion under water flow discharge was assumed to be small. A new erosion law fixing the initiation threshold of erosion and integrating the effect of soil consolidation on internal erosion was proposed. Conversely, the effect of erosion on elastic soil deformation was also integrated through damage mechanics concepts. Asymptotic expansion of the coupled equations in terms of a perturbation parameter linked to the total amount of internal erosion that is likely to occur has been performed. This has enabled to view the internal erosion phenomenon occurring inside the soil as a perturbation affecting the classical soil consolidation equation, and further to evaluate the critical discharge gradient for which internal erosion starts. Equations at order zero that are provided by the asymptotic expansion were exactly integrated while an adequate finite difference scheme was introduced to solve the equations at order one. A parametric study was conducted after that in order to assess effects of the main factors on internal erosion and soil deformation. Ó 2014 Elsevier Inc. All rights reserved.

1. Introduction Stability of soils under buildings or hydraulic infrastructures such as levees or dams is of great importance in practice. Even if initially design of foundations meets well safety criteria, problems of durability can arise. One of the most important sources of degradations yielding to instabilities of foundations is internal erosion. This phenomenon is related to the migration of soil particles through the soil matrix by suffusion or piping. This phenomenon is generated by water seepage as the hydraulic gradient reaches a critical value. It starts as suffusion of soil fines that takes place inside the foundations of hydraulic retaining structures such embankment dams and levees, then progresses until reaching the advanced stage of piping. This last step consists of regressive erosion which develops from downstream of the critical seepage line towards its upstream with the formation of a continuous pipe, [1]. The risk threatening hydraulic infrastructures due to internal erosion is aggravated by the sudden character of this phenomenon. No external sign enabling its early detection is visible before piping occurrence. A dam may breach in fact only a ⇑ Corresponding author. Tel.: +212 600769960; fax: +212 539994500. E-mail addresses: [email protected] (M.A. Parron Vera), [email protected] (F. Yakhlef), [email protected] (M.D. Rubio Cintas), [email protected] (O. Castillo Lopez), [email protected] (P. Dubujet), [email protected] (A. Khamlichi), [email protected] (M. Bezzazi). http://dx.doi.org/10.1016/j.apm.2014.02.006 0307-904X/Ó 2014 Elsevier Inc. All rights reserved.

Please cite this article in press as: M.A. Parron Vera et al., Analytical solution of coupled soil erosion and consolidation equations by asymptotic expansion approach, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.02.006

2

M.A. Parron Vera et al. / Applied Mathematical Modelling xxx (2014) xxx–xxx

few hours after first evidence of piping has been stated. Understanding the early stages of internal erosion during the suffusion phase is of huge interest from the practical point of view. Predicting internal erosion can give in this way the possibility to perform in better way health monitoring of hydraulic infrastructures through assessing active protection of these vital facilities. As pointed out in [2], suffusion occurs when a soil satisfies the following criteria: it has a structure that possesses both a coarser fraction containing voids and a finer fraction; the loose particles constituting the finer fraction are smaller than the constrictions formed by the coarser fraction, and finally water flow velocity is strong enough to transport the fines through the voids. Several physical models that describe internal erosion have been derived in the last decades [3–6]. In most of the proposed models, a continuum approach is adopted. The soil is viewed as a three-phase medium with the representative volume element consisting of solid particles, discharging fluid and fluidized grains which are extracted by erosion from the soil skeleton [3]. Internal erosion is activated locally by a seepage flowing with a supercritical gradient and was recognized to be strongly coupled with the hydro-poro-mechanical behaviour of the soil. The effect of internal erosion on the soil was integrated by introducing a tangent elastic plastic porosity matrix that affects soil deformation and by taking into account variations of permeability as function of porosity [7]. Based on experimental evidence that states that mass production of eroded particles tends to decrease over time, an erosion law had been suggested by assuming that the production rate of fluidized grains is proportional to the hydraulic gradient discharge [6–8,3]. Under this form, the equation of erosion was coupled to Darcy equation and thus to consolidation phenomenon taking place in the soil. Most of the basic ideas existing in the previous literature dealing with modelling internal erosion are carried out in the actual study. But, since focus will be done on the early stage of erosion, the soil skeleton is assumed to deform elastically in the vicinity of an equilibrium point that has been reached during the pre erosion stage. In addition to that, variations of porosity resulting from erosion are supposed to yield damage of the soil. Even if the determination of the thermodynamical potential associated to this damage mechanism is still an open problem, the soil Young’s modulus is allowed to vary according to the classical law of damage mechanics [10,11]. Using an extension of the erosion law that was first formulated in [7], the authors of the present work have given one-dimensional approximation of internal erosion equations [11]. This was performed in the special case of soil samples made from controlled mixtures of sand and clay but containing a small proportion of clay that is susceptible to experience erosion. The equations governing internal erosion problem were found to have the general form of coupled equations that include erosion, poro-elastic consolidation and transport of eroded particles. The big problem with the obtained equations is that they are instable and their direct numerical integration through the standard finite element method has given rise to poor results. This has been well assessed earlier, since numerical models based on coupled formulations even in the simplest soil consolidation problem without any erosion were recognized to become ill conditioned when the soil permeability is small or soil deformability is high. To overcome these numerical difficulties special treatment was proposed in the literature [12] by defining a new element through an averaging technique. This is why even for the decoupled problem analytical approaches are still extensively used to deal with soil consolidation. Xie et al. [13] have developed a general approximate analytical solution for one-dimensional consolidation with consideration of the threshold gradient under a time-dependent loading, and the approach was found to be quite reasonable. Considering large strain consolidation, Xie et al. [14] have given an explicit analytical solution for the one-dimensional problem. Comparisons that were made with the classical small strain theory have shown that the average degree of consolidation is modified. Using Laplace transform, Qin et al. [15] have presented an analytical solution to one-dimensional consolidation in unsaturated soils with the applied vertical loading varying exponentially with time. Due to the complexity of the non-linear consolidation of soft clay, other robust approaches were proposed like the differential quadrature method which was used in [16] to solve one-dimensional nonlinear consolidation problem for a multi-layered soil. Zheng et al. [17] have used also this method to implement solution of one dimensional nonlinear consolidation equation with various boundary conditions, and the obtained numerical results compared satisfactorily with analytical solutions. In the framework of fully coupled equations that comprise both internal erosion and consolidation of soils, literature dedicated to this subject is actually rare. Among the few works in this field, Wang and Wan [18] have observed numerical difficulties while using the standard finite element method to integrate the problem. These instabilities had taken the form of wiggles that had arisen in the case of high gradient conditions that follow large variation of the local field variables. In order to capture the associated local sharp changes, the field variables that enter the differential equations modelling the problem were averaged around material points under consideration. In a previous work [11], the authors have introduced an asymptotic expansion approach in order to solve in a stable way the one dimensional soil coupled internal erosion-consolidation equations. Within the context of this approach, the internal erosion consolidation problem is viewed to be as the classical consolidation problem perturbed by porosity variations. This suggested performing asymptotic expansion of the obtained equations in terms of a scale parameter chosen to be the magnitude of expected porosity variations. The classical consolidation problem appeared like this to coincide with the degenerated erosion free problem. Using the asymptotic expansion, the problem was shown to have a closed form analytical solution. To reproduce experimental evidence stating that erosion appears only when the seepage velocity is supercritical, a further extension is considered here for the erosion law which was first introduced in [11]. A threshold criterion for minimum gradient required to initiate erosion is formulated. The asymptotic expansion technique is used with some modifications that Please cite this article in press as: M.A. Parron Vera et al., Analytical solution of coupled soil erosion and consolidation equations by asymptotic expansion approach, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.02.006

M.A. Parron Vera et al. / Applied Mathematical Modelling xxx (2014) xxx–xxx

3

account for the new erosion law. Simulations are conducted after that in order to study the influence of key factors on the obtained erosion pattern, and in particular the influence owing to the critical gradient initiating internal erosion. 2. One dimensional modelling of coupled internal erosion-consolidation problem in soils Considering the special case of a soil sample that consists of a mixture made from sand and a small proportion of clay that is susceptible to be removed by erosion, a fully coupled internal erosion model is proposed. The model is considered in the case of one-dimensional approximation of the problem. Part of the equations that are used in the subsequent for modelling internal erosion is based on relevant previous works that were performed in this field such as those presented in [5] and [7]. A modification of these equations was introduced by the authors in [11]. It consisted of an extension of the erosion law and the consideration of damage resulting from erosion in the elastic soil behaviour. In the following further extension is proposed for the erosion law by introducing a critical threshold under which erosion cannot occur. 2.1. Quantitative description of the three-phase soil medium The internal erosion model considered in the following assumes the three-phase continuum approach [3]. For a given volume of soil, the three phases occupying it are: - the skeleton consisting of sand and clay grains that are not separated from it by erosion, its volume is denoted by Vs; - the fluid which is here water having volume Vw; - the fluidized solid, that are the grains of clay detached from the skeleton by the effect of discharging flow forces and which has volume Ver. - The relative volumetric composition of the three-phases is defined by means of soil porosity u and concentration of the fluidized solid c. Soil porosity is defined as

u¼

Vv V

ð1Þ

with Vv = Vw + Ver representing the volume of voids.Concentration of fluidized grains is given by,

c¼

V er Vv

ð2Þ

The apparent densities of respectively solid, water and fluidized grains phases are,

qs ¼

ms ; V

qw ¼

mw ; V

qer ¼

mer V

ð3Þ

where ms, mw and mer are, respectively the masses of solid, water and fluidized grains phases that are contained in the considered soil volume V. Denoting qs, qw and qer the real densities of, respectively solid, water and fluidized grains phases, one can derive easily the following relationships,

qs ¼ qs ð1  uÞ; qw ¼ qw uð1  cÞ; qer ¼ qer cu

ð4Þ

and shows that the soil density is given by,

q ¼ qs þ qw þ qer ¼ qs ð1  uÞ þ qw u þ ðqer  qw Þcu

ð5Þ

2.2. Kinematics of the three-phases In modelling kinematics of eroded particles it is assumed that the detached grains move immediately at the velocity of adjacent water particles. This is only an approximation since in reality kinematics of fluidized grains is complex and their velocity is in general small than that of the fluid [19]. Thus

v er ¼ v w

ð6Þ

er

w

where v is the velocity of an eroded grain and v is the velocity of adjacent fluid particles. The concept of Darcy velocity which represents a fraction of the true fluid particle speed and which is associated to the apparent percolation rate is used in the following [20]. The partial volume discharge flows of respectively water and fluidized grains (Darcy velocities) are given by,

qw ¼ v w uð1  cÞ;

qer ¼ v w cu:

ð7Þ

Thus, the total discharge is

q ¼ qw þ qer ¼ v w u

ð8Þ

Please cite this article in press as: M.A. Parron Vera et al., Analytical solution of coupled soil erosion and consolidation equations by asymptotic expansion approach, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.02.006

4

M.A. Parron Vera et al. / Applied Mathematical Modelling xxx (2014) xxx–xxx

2.3. Continuity equations Denoting j the rate associated to loss of solid skeleton mass per unit time due to erosion and thus entering in the water domain, conservation of the solid phase mass writes

@ qs @ðqs v s Þ þ ¼ j @x @t

ð9Þ

where vs is the solid velocity, t time and x the axial coordinate. It is assumed that the whole clay grains that are extracted by the discharging water flow are added to the fluidized solid phase. Considering mass conservation for this phase, one obtains the following equation,

@ qer @ðqer v w Þ þ ¼j @x @t

ð10Þ

Finally, mass conservation for water phase writes

@ qw @ðqw v w Þ þ ¼0 @x @t

ð11Þ

Eqs. (4), (9), and (10) yield

@ðuðg  cÞÞ @ðcqÞ @u @ u @2u  þg  gð1  uÞ ¼0 @t @x @t @x @x@t where u is the solid phase displacement which is related to solid phase velocity Using Eqs. (4), (11), and (12) one arrives at

ð12Þ

vs = o u/ o t and g = qs/qer.

@q @2u @u @ u @u þ gð1  uÞ g þ ð1  gÞ ¼0 @x @x@t @t @x @t

ð13Þ

Assuming that qs ’ qer which is not far from reality, one obtains the following approximation g ’ 1 and Eqs. (12) and (13) reduce to,

@q @ 2 u @u @ u þ ð1  uÞ  ¼ 0; @x @x@t @t @x ! @c @c @ u @u @ u @2u ; u þ q ¼ ð1  cÞ þ  ð1  uÞ @t @x @t @x @x@t @t

ð14Þ

ð15Þ

The two Eqs. (14) and (15) have the four following unknowns: u, /, q and c. To perform closure of the problem, Darcy equation and soil equilibrium condition are considered. 2.4. Darcy equation The experimental law of Darcy gives water volume discharge that flows across the porous soil sample per unit time. Since the homogenized density of water and fluidized solid phases is given by qf = qsc + qw(1  c), Darcy velocity writes [20]

q¼

k @p g qf @x

ð16Þ

where p is the pore pressure, k the hydraulic permeability and g the constant of gravity. In practice c  1 and qf ’ qw constitutes a good approximation of fluid density. Variation of soil permeability as function of porosity is assumed to be governed by the Kozeny–Carman law [7]

k ¼ k0

u3 ð1  uÞ2

ð17Þ

where k0 is the initial soil permeability.Using Eqs. (16) and (17) with the approximation qf ’ qw one obtains

q¼

k0 u3

@p g qw ð1  uÞ2 @x

ð18Þ

2.5. Soil equilibrium Soil displacement is assumed to remain quasi-static. The stress state inside the soil sample is determined by considering Terzaghi principle which introduces the concept of effective stress [20]. According to this principle, the total stress filed is given by r = r’  p where r’ is the effective stress that governs soil behavior and p the pore pressure. Please cite this article in press as: M.A. Parron Vera et al., Analytical solution of coupled soil erosion and consolidation equations by asymptotic expansion approach, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.02.006

M.A. Parron Vera et al. / Applied Mathematical Modelling xxx (2014) xxx–xxx

5

Erosion phenomenon can be viewed to be associated to a degradation mechanism that happens inside the soil medium and which affects essentially its porosity. This last increases following grain detachment from the solid skeleton. As the purpose of this work is to predict the early stage of internal erosion, one may assume that the mechanical behaviour of the soil skeleton is almost elastic. Local linearization of soil constitutive equations around the actual equilibrium point can be performed for that.Applying the generalized Hook’s law to a linear isotropic homogeneous material yields

ex ¼

 1 rx  tðry þ rz Þ ; E

ey ¼

 1 ry  tðrz þ rx Þ ; E

ez ¼

 1 rz  tðrx þ ry Þ E

ð19Þ

where E is the Young’s modulus, t Poisson’s coefficient, e. the stains and r. the stresses. Considering a cylindrical soil sample having a configuration that obeys axial symmetry about x axis with confined lateral deformation for which ey = ez = 0 and ry = rz one obtains by using Eq. (19)

ry ¼ rz ¼

t 1t

rx ; rx ¼

Eð1  tÞ ex ð1 þ tÞð1  2tÞ

ð20Þ

Erosion taking place inside the soil evolves causing the deterioration of stiffness. This deterioration is normally referred to as damage [9] and [10]. Isotropic damage is generally assumed to simplify numerical calculations. To account for damage mechanism resulting from internal erosion, the net soil Young’s modulus is assumed to decrease linearly as function of porosity. Applying Eq. (20) to the soil sample according to Terzaghi principle amid damage yields

r0 ¼ E1 ð1 þ u0  uÞ

@u @x

ð21Þ

with

E1 ¼

Eð1  tÞ ð1 þ tÞð1  2tÞ

ð22Þ

the modulus of soil compressibility where E designates the pre-erosion soil Young’s modulus as measured for the initial porosity u0. As the local equilibrium of soil sample writes o r/ o x = 0, it follows then from Terzaghi postulate and Eq. (21) that

  @ @u @p ¼ E1 ð1 þ u0  uÞ @x @x @x

ð23Þ

Integrating Eq. (23) with respect to x yields the following equation

@u 1 ¼ ðp þ r0 Þ @x E1 ð1 þ u0  uÞ

ð24Þ

2.6. Erosion law Classical laws for internal erosion were proposed by putting either the erosion rate to be proportional to fluid discharge [6], [8] and [3] or by putting the particle discharge to be proportional to porosity gradient [19]. The advantage of this last is that, for small eroded particle concentration, erosion resulted in a porosity diffusion law which was found to uncouple the porosity from the coupled set of equations describing elastic–plastic soil consolidation, yielding the possibility of simple integration of equations. The problem is that this erosion law was found to be insufficient to explain all the observed experimental results [19]. In the subsequent an extension of the erosion law that was formulated first in Papamichos et al. [7] and modified in [11] is proposed under the following form

  @ eer 1 u a ¼ qS 1 hjqj  qer i uc @t ser

ð25Þ

where a is a positive real exponent, ser is a soil parameter measuring soil resistance to erosion, qer the critical Darcy velocity which indicates the onset of erosion, uc the maximal porosity that could occur owing to the fact that the soil skeleton contains a non erodible part and h  i = max (  , 0). In the erosion law as stated by Eq. (25), the erosion process initiation is not free and is controlled by the qer parameter. This equation gives rise to further nonlinearities as compared with its counterpart, Eq. (6) in Ref. [11]. The total time variation of local soil porosity is equal to volume time variation induced by erosion plus strain volume time variation, that is to say

@ u @ eer @ ex ¼ þ @t @t @t

ð26Þ

where eer is the volume variation produced by internal erosion. Please cite this article in press as: M.A. Parron Vera et al., Analytical solution of coupled soil erosion and consolidation equations by asymptotic expansion approach, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.02.006

6

M.A. Parron Vera et al. / Applied Mathematical Modelling xxx (2014) xxx–xxx

Substituting Eqs. (18) and (25) into Eq. (26) one obtains

+   *  @u q k0 u3 u a @p g qw ð1  uÞ2 @2u þ ¼ S 1   q er   2 3 uc @x @x@t @t g qw ser ð1  uÞ k0 u

ð27Þ

3. Asymptotic expansion of coupled internal erosion and consolidation equations Deriving Eq. (24) with respect to time yields,

  @2u 1 @p @ r0 p þ r0 @u þ ¼ þ @x@t E1 ð1 þ u0  uÞ @t @t E1 ð1 þ u0  uÞ2 @t

ð28Þ

Substituting after that Eq. (28) into Eq. (14) gives

  @q ð1  uÞ @p @ r0 ð1  uÞðp þ r0 Þ @ u @u @ u  ¼ þ þ @x E1 ð1 þ u0  uÞ @t @t @x @t E1 ð1 þ u0  uÞ2 @t

ð29Þ

Deriving now Eq. (18) with respect to x and substituting Eq. (24) into the obtained result gives

@2p 1 @p 1 @ r0 1 @u 3  u @ u @p E1 ð1 þ u0  uÞ @u @ u ¼ þ þ p  þ @x2 C v ðuÞ @t C v ðuÞ @t C v ðuÞð1 þ u0  uÞ @t uð1  uÞ @x @x C v ðuÞð1  uÞ @t @x r0 @u þ C v ðuÞð1 þ u0  uÞ @t

ð30Þ

with

C v ðuÞ ¼

k0 E1 u3 ð1 þ u0  uÞ

qw gð1  uÞ3

ð31Þ

Eqs. (24), (27), and (30) form a system of three independent equations having the three unknowns u, u and p. When solution of these three equations is obtained, Eq. (15) can be used to compute the concentration of the fluidized grains and Eq. (18) to compute the Darcy velocity q. After extensive investigation trying to integrate directly system of Eqs. (24), (27), and (30) by means of the standard finite element method, numerical instabilities have always emerged yielding to systematic divergence of results. This is known to be inherent to this kind of coupled equations because they are ill-conditioned as was recognized in some previous works [12–14] and [21]. Hence, asymptotic expansion technique is considered in the following in order to transform the erosion consolidation equations into a more stable form that is suitable for numerical integration. Examining Eq. (30) one can notice that it consists of the classical consolidation equation that appears to be perturbed by erosion taking place inside the soil sample: variation of porosity u. Assuming that the porosity is constant u = u0 and the case where r(t) = r0 is constant, Eq. (30) yields the classical consolidation equation

@2p 1 @p ¼ @x2 C v 0 @t

ð32Þ

where Cv0 = Cv(u0) is the consolidation coefficient. For further details about consolidation theory based on Terzaghi theory [15] or the nonlinear theory, the reader can consult references [14], [21] and [22]. The above remark invites us to choose as the scale parameter that governs the perturbed consolidation problem associated to Eqs. (24), (27), and (30), the following expression which is linked to the maximal porosity that could be reached due to the erosion



e¼ 1

u0 uc

a ð33Þ

where a is a scalar constant and uc the maximum porosity than could result from erosion. Parameter e is equal to zero when there is no possibility of internal erosion inside the soil sample, that is to say uc = u0. In order to simplify more the model, the initial porosity is assumed to be uniform over the whole soil sample, supposed to be free from any erosion, thus the initial gradient of porosity vanishes and

@ u0 @ u0 ¼ ¼0 @t @x

ð34Þ

As an approximation, the asymptotic expansion in terms of e is considered to be truncated at the first order, thus

p ¼ p0 þ ep1 u ¼ u0 þ eu1

ð35Þ

u ¼ u0 þ eu1 Please cite this article in press as: M.A. Parron Vera et al., Analytical solution of coupled soil erosion and consolidation equations by asymptotic expansion approach, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.02.006

M.A. Parron Vera et al. / Applied Mathematical Modelling xxx (2014) xxx–xxx

7

By substituting the asymptotic expansion (35) into Eqs. (24), (27), and (30), it can be shown easily that only the case where a = 1 will yield a polynomial expansion. One obtains in this case the following equations at orders 0 and 1.

@ 2 p0 1 @p0 ¼ ; C v 0 @t @x2

ð36Þ

@u0 p0 þ r0 ¼ ; @x E1

ð37Þ

@ u1 q k0 u30 ¼ s @t g qw ser ð1  u0 Þ2

*  !+  u0  u1 @p0 g qw ð1  u0 Þ2 @ 2 u1 þ  þ qer ; u0 @x @t@x k0 u30

ð38Þ

@ 2 p1 1 @p1 u1 ð3  u0 þ u20 Þ @p0 qw gð1  u0 Þ3 ðp0 þ r0 Þ @ u1 3  u0 @p0 @ u1 qw gð1  u0 Þ3 ¼  þ  þ 3 C v 0 @t C v 0 u0 ð1  u0 Þ @t @x2 @t u k0 u30 k0 E1 u0 0 ð1  u0 Þ @x @x @u0 @ u1  ; @t @x @u1 1 ¼ ðp1 þ u1 ðp0 þ r0 ÞÞ: E1 @x

ð39Þ

ð40Þ

When r0 is specified as well as the initial and boundary conditions, Eqs. (36) and (37) can be integrated analytically to give explicitly the basic terms p0 and u0 of the asymptotic expansion appearing in Eq. (35). To get the first order terms of the asymptotic expansion, numerical integration of the partial differential equations (38) and (40) should be considered. This is performed in the following via a finite difference scheme. In order to reduce the number of unknown variables, one can notice that the variable u1 can be eliminated from Eqs. (38) and (40). Thus, substituting Eq. (40) into Eq. (38) and multiplying the result by E1 gives

*  !+    p0 þ r0 @ u1 qs k0 u30 u0  u1 @p0 g qw ð1  u0 Þ2 1 @p0 1 @p1 þ 1 ¼  þ qer u þ E1 @t 1 E1 @t E1 @t g qw ser ð1  u0 Þ2 u0 @x k0 u30

ð41Þ

4. Analytical solution of coupled erosion and consolidation equations Solution of Eqs. (36)–(40) is developed in the following in the situation where the left hand side of the soil sample x = 0, Fig. 1, is suddenly subjected to an increase of pressure from the initial value pei to the value pe while the right hand side x = L is assumed to be fixed against the axial displacement and subjected to the constant pressure ps. At the initial time t = 0, the pressure is assumed to be linearly distributed over the domain of the soil sample, while the total stress r0 is taken to be constant. The initial conditions associated to Eqs. (36)–(41) are,

x L

u0 ðx; t ¼ 0Þ ¼ ui ; p0 ðx; t ¼ 0Þ ¼ pei þ ðpsi  pei Þ; r0 ðx; t ¼ 0Þ ¼ r0 ;

ð42Þ

p1 ðx; t ¼ 0Þ ¼ u1 ðx; t ¼ 0Þ ¼ 0 The boundary conditions applying at the left hand extremity x = 0 which correspond to a given pressure, vanishing strain as well as vanishing pressure perturbation gradient are

p0 ðx ¼ 0; tÞ ¼ pe ;

p1 ðx ¼ 0; tÞ ¼ 0;

r0 ðx ¼ 0; tÞ ¼ r0

ð43Þ

The boundary conditions applying at the right extremity x = L that correspond to given pressure, and vanishing displacement and vanishing pressure perturbation gradient are,

p0 ðx ¼ L; tÞ ¼ ps ; u0 ðx ¼ L; tÞ ¼ 0;

p1 ðx ¼ L; tÞ ¼ 0;

u1 ðx ¼ L; tÞ ¼ 0:

ð44Þ

Eqs. (36) and (37) along with the initial and boundary conditions (42) to (44) have the following closed form solutions

 2 2  1 npx

x

x 2X pe  pei  psi ð1Þn n p Cv0t sin p0 ðx; tÞ ¼ pe 1  exp  þ psi  2 L L p n¼1 L n L

0

ð45Þ

x L

Fig. 1. Geometry of the soil sample.

Please cite this article in press as: M.A. Parron Vera et al., Analytical solution of coupled soil erosion and consolidation equations by asymptotic expansion approach, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.02.006

8

M.A. Parron Vera et al. / Applied Mathematical Modelling xxx (2014) xxx–xxx

2 2

n p Cv0 t 1 exp  h npx

i p ðL  xÞ2 p 2L X L2 u0 ðx; tÞ ¼ 2 ðpe  pei  psi ð1Þn Þ cos  si  ððpe  pei Þð1Þn  psi Þ  e 2 L p E1 n¼1 E1 2L E1 n 

L2  x2 r0 ðL  xÞ  E1 2L

ð46Þ

Eqs. (39) and (41) can now be integrated numerically by means of a finite difference scheme to calculate p1 and u1. Denoting Dt and Dx the time and spatial discretisation steps and giving a grid point having coordinates (mDx, nDt), finite approximations of the partial derivations appearing in Eqs. (39) and (41) are considered under the following form.

ðu Þnþ1  ðu1 Þnm @ u1 ðmDx; nDtÞ ¼ 1 m ; @t Dt

ð47Þ

@p1 ðp Þnþ1  ðp1 Þnm ðmDx; nDtÞ ¼ 1 m ; @t Dt

ð48Þ

n

n

ðu Þ  ðu1 Þm @ u1 ðmDx; nDtÞ ¼ 1 mþ1 ; @x Dx n

ð49Þ

n

n

ðp Þ  2ðp1 Þm þ ðp1 Þm1 @ 2 p1 ðmDx; nDtÞ ¼ 1 mþ1 ; @x2 Dx2

ð50Þ

where ðÞnm ¼ ðÞðmDx; nDtÞ 2

u0 Þ 0 As seen from Eq. (27), for the erosion to take place, the condition  @p > gqwkð1 qer should be satisfied. Deriving Eq. (45) @x u3 0

0

with respect to x and substituting the result in the previous condition, one arrives at the following necessary condition that specifies the lower bound of the applied pressure pe for which the erosion occurs

pe >

g qw ð1u0 Þ2 L qer k0 u30

2 2

P n n p Cv0 t þ psi þ 2 1 n¼1 ðpei þ psi ð1Þ Þ exp  L2

P1 n2 p2 C v 0 t 1 þ 2 n¼1 exp  L2

ð51Þ

Substituting Eqs. (47)–(50) into Eqs. (39) and (41) gives at time tn+1 = (n + 1) Dt the following solutions for the interior points n

n

ðu1 Þnþ1 ¼ m

ða0 Þnm ðd0 Þm ðs1 Þnm  ðd0 Þm ðr1 Þnm  n n n ðd0 Þm þ ða0 Þnm ðe0 Þnm  ðd0 Þm ðb0 Þm Dt

ðp1 Þnþ1 ¼ m

ðe0 Þnm ðr1 Þnm Dt þ ðd0 Þm ½1  ðb0 Þm Dtðs1 Þnm  n n n ðd0 Þm þ ða0 Þnm ðe0 Þnm  ðd0 Þm ðb0 Þm Dt

n

n

ð52Þ

n

ð53Þ

n

where the quantities ða0 Þnm ; ðb0 Þm ; ðd0 Þm ; ðe0 Þnm ; ðr1 Þnm and ðs1 Þnm are given in the Appendix at the actual time tn = nDt as function of known quantities. For the boundary points x = 0 and x = L, the following expressions are obtained by applying the boundary conditions in terms of p1

ðu1 Þnþ1 ¼ 1

ðr 1 ÞnMþ1 ðr 1 Þn1 and ðu1 Þnþ1 n n Mþ1 ¼ 1  ðb0 Þ1 Dt 1  ðb0 ÞMþ1 Dt

ð54Þ

Eqs. (37) and (40) can now be used in order to calculate the axial strain as function of time and axial coordinate. The 0 Appendix gives the analytical expression of @u . @x 2

u1 Now, assuming on one hand that the soil sample is relatively rigid, so as the term @@t@x can be neglected, and on the other hand that no erosion occurs, Eq. (27) shows that if the initial porosity is uniform then it will remain constant in time. Eq. (30) degenerates in this case to that of a simple linear consolidation problem which is governed by Eq. (36). This particular case is viewed here as a limit problem that can be recovered directly as a particular case from the asymptotic expansion Eqs. (36)–(40).

5. Results and discussion A homogeneous soil sample occupying a cylindrical domain having the length L = 5 m, Fig. 1, is considered. The soil is assumed to be at a static initial equilibrium state with a uniform initial porosity ui = 0.35, when its left side extremity is subjected to pressure pei and its right side extremity to pressure ps. At the instant t = 0 a pressure increment of magnitude Dp = 3

3

pe  pei is applied at the left side extremity of the sample. The other soil data are: qw ¼ 1000 kg:m , qs ¼ 2650 kg:m , E0 ¼ 2  108 Pa, t = 0.25 and uc = 0.38. Please cite this article in press as: M.A. Parron Vera et al., Analytical solution of coupled soil erosion and consolidation equations by asymptotic expansion approach, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.02.006

M.A. Parron Vera et al. / Applied Mathematical Modelling xxx (2014) xxx–xxx

9

Porosity

0.365 0.36 0.355 3000

0.35

2000

0

1

1000

2

3 Coordinate x (m)

4

0

5

Time (s)

Fig. 2. Variations of porosity as function of axial position and time.

5

x 10

Pressure (Pa)

6 4 2 0 3000 -2 0

2000 2

4

Coordinate x (m)

1000 6

0

Time (s)

Fig. 3. Variations of pressure as function of axial position and time.

0.365

x= 0.5m

Porosity

0.36 x= 1m

0.355 x= 1.5 m

x= 2 m

0.35 0

500

1000

1500 Time (s)

2000

2500

3000

Fig. 4. Time history of porosity at different locations.

Please cite this article in press as: M.A. Parron Vera et al., Analytical solution of coupled soil erosion and consolidation equations by asymptotic expansion approach, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.02.006

10

M.A. Parron Vera et al. / Applied Mathematical Modelling xxx (2014) xxx–xxx 5

4

x 10

x=1 m

3.5

Pressure (Pa)

3 2.5

x= 2 m

2 x= 3 m

1.5 1

x= 4 m

0.5 0 0

500

1000

1500

2000

2500

3000

Time (s) Fig. 5. Time history of pressure at different locations.

0.365

0.3504

0.3501 0.35 0

Porosity

Porosity

0.3502

3000 1

2

3

5 0

0.355

0

1000 4

0.36

1

2

3

1000 4

5 0

Coordinate x

Coordinate x

(a) q er = 5 × 10−6 m.s −1

(b) q er = 9 ×10−7 m.s −1

0.365

0.36

0.36

0.355 3000

0.35

2000 1

2000

Time

0.365

0

3000

0.35

2000

Porosity

Porosity

0.3503

2

3

1000 4

5 0

Time

0.355 3000

0.35 0

2000 1

2

3

Time

1000 4

5 0

Coordinate x

Coordinate x

(c) q er = 7 ×10−7 m.s −1

(d) q er = 5 ×10−7 m.s −1

Time

Fig. 6. Influence of the critical gradient on the erosion pattern.

Please cite this article in press as: M.A. Parron Vera et al., Analytical solution of coupled soil erosion and consolidation equations by asymptotic expansion approach, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.02.006

M.A. Parron Vera et al. / Applied Mathematical Modelling xxx (2014) xxx–xxx

11

-3

x 10 -1.04 -1.06

Strain

-1.08 -1.1 -1.12

No erosion Partial erosion Full erosion

-1.14 -1.16 -1.18 -1.2 0

2000

4000 6000 Time (s)

8000

10000

Fig. 7. Influence of erosion on soil consolidation strain as obtained for x = 3.75 m; partial erosion corresponds to qer ¼ 7  107 m:s1 and full erosion is associated to qer = 0 ms1.

An application is considered in the following for the particular case where k0 = 106, pe=5  105 Pa, pei = 4  105 Pa, ps = 0 Pa, r0 = 4  105 Pa, qer = 8  107 ms1 and ser = 1 kgm2. In these conditions C v 0 ¼ 2:883  103 m2 :s1 and E1 ¼ 2:4  108 Pa. Taking the number of spatial nodes N = 21, 90 terms in the series of Eqs. (45) and (46), the time interval duration T = 3  103 s and the time step Dt ¼ 5 s, the pressure and the porosity surfaces as function of time t and coordinate x are shown in Figs. 2 and 3. As shown in Fig. 2, the pressure gradient converges towards a constant value as pressure tends to a linear profile when time is sufficiently high. Fig. 3 shows that the alteration of porosity is high in the vicinity of the inlet extremity of the soil sample and it diminishes rapidly within a small distance from it. This indicates that the erosion process is more active at the beginning for points that are close to the extremities where high pressure is applied. But with elapsing time, elastic deformation increases because of damage mechanism reducing soil Young’s modulus. In the same time an increase of permeability follows any increase of erosion according to the Kozeny-Karman equation (17), which yields then to a decrease of erosion. Due to these concomitant mechanisms that reduce porosity, a stationary state is reached after some elapsed time and a uniform gradient of porosity extends on the whole soil sample length. The next plots, Figs. 4 and 5 depict the time history of porosity and pressure at different locations of the soil sample. Fig. 4 shows that porosity increases rapidly at the beginning for points located near from the left soil sample extremity and attains after that an asymptotic maximum. For the other points that are far from this extremity, porosity is almost unaffected. Fig. 5 shows that for all locations, the pressure increases with high rate at the first instants, then the rate decreases and the pressure tends towards a stationary value. The final pressure profile will be that for which a constant pressure gradient is reached. If we keep now all the parameters unchanged and modify only the critical Darcy velocity qer for which erosion starts to develop in the soil sample, then the various erosion patterns that are represented in Fig. 6 are obtained. No erosion takes place for the case (a) as the critical Darcy velocity is not reached in the soil sample. Porosity variations which remain in this case small are elastic and result from only consolidation. Decreasing the threshold for which erosion starts yields a net augmentation of erosion. The erosion pattern will localize at early time near from the inlet extremity for high values of qer, but for small values of this parameter, erosion will spread rapidly larger into the sample and finally attains a uniform gradient regime for qer < 5  107 m:s1 . Effect of erosion on consolidation phenomenon of the soil sample can further be emphasized through studying the stationary achieved consolidation strain as function of erosion amount that is likely to take place in the medium. With the same data than above, Fig. 7 shows the obtained stationary consolidation strain for three different erosion cases. This figure shows that erosion damages the soil sample; this reduces its effective rigidity which provokes a decrease of the reached strain settlement increment. Fig. 7 shows that settlement is modified when the soil undergoes erosion. This offers the possibility to use this information in practice as an indicator that enables to monitor hydraulic infrastructures. Using classical analytical formulas [20] that give the consolidation strain according to the linear theory one obtains exactly the same curve than that associated to the erosion free case as given in Fig. 7. 6. Conclusion An improved analytical modeling of internal erosion phenomenon occurring in a soil sample that consists of sand and clay mixture was performed under the assumption of one-dimensional approximation. An enhanced erosion law integrating eroPlease cite this article in press as: M.A. Parron Vera et al., Analytical solution of coupled soil erosion and consolidation equations by asymptotic expansion approach, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.02.006

12

M.A. Parron Vera et al. / Applied Mathematical Modelling xxx (2014) xxx–xxx

sion initiation threshold was proposed. The modeling uses the continuum approach and is based on the interpretation that the soil coupled equations involving both internal erosion and consolidation are those of a perturbed classical consolidation problem. The perturbation was directly linked to porosity variation and the amplitude of this variation set scale factor. An asymptotic expansion technique was used in order to derive the first order approximation of the coupled equations, for which the solution has been obtained in a closed form. This has enabled to conduct a parametric study regarding the influence of the critical gradient on the erosion pattern. Localized erosion regimes were obtained for moderate erodable soils. For low erosion resistant soils a uniform gradient of erosion over the soil sample was obtained. The proposed modelling enabled quantifying the influence of erosion initiation threshold on erosion outline and showed how this last is sensitive to the key intervening factors. However, to be valid this modelling requires that erosion should remain small in order that the first order asymptotic expansion could be acceptable. This can further be addressed through considering higher order asymptotic expansions. Acknowledgments Part of this study was financed by research grants accorded by the Spanish Agency AECID, the French CNRST and the Moroccan CNRST. The authors would like to thank them warmly. Appendix A The quantities p0 and u0 are given explicitly as function of t and x by Eqs. (45) and (46). The partial derivatives of these quantities that intervene in the problem are obtained explicitly under the following series form

 2 2  1 npx

 @p0 2p C v 0 X n p Cv0t n sin ðx; tÞ ¼ n p  p  p ð1Þ exp  e ei si L @t L2 L2 n¼1

ðA:1Þ

 2 2  1 npx

 @p0 p  pe 2 X n p Cv0t cos ðx; tÞ ¼ si  pe  pei  psi ð1Þn exp  2 L n¼1 L @x L L

ðA:2Þ

 2 2  1 h npx

i @u0 2C v 0 X n p Cv0t ðpe  pei  psi ð1Þn Þ cos  ð1Þn ðx; tÞ ¼  exp  2 L @t E1 L n¼1 L

ðA:3Þ

 2 2  1 npx p ðL  xÞ p x r @u0 2 X 1 n p Cv0t 0 ðpe  pei  psi ð1Þn Þ sin ðx; tÞ ¼  exp  þ si þ þ e 2 L L @x pE1 n¼1 n E1 E1 L E1 L

ðA:4Þ

Giving the solution at the time instant tn = nDt, the quantities ðp1 Þnm and ðu1 Þnm are known. Eqs. (45), (A.1), (A.2), (A.3), and (A.4) enable then to calculate the solution at time tn + 1 = (n + 1)Dt as specified in Eqs. (52)–(54) through using the following quantities

*

qS k0 u30 @p0 g qw ð1  u0 Þ2 wer ¼   qer @x k0 u30 g qw ð1  u0 Þ2 ser

+ ðA:5Þ

a0 ¼

1 E1  p0  r0

ðA:6Þ

b0 ¼

1 @p0 E1 wer  E1  p0  r0 @t E1  p0  r0 u0

ðA:7Þ

c0 ¼

E1 wer E1  p0  r0

ðA:8Þ

d0 ¼

  1 E1 C v 0 E1  p0  r0

ðA:9Þ

 e0 ¼

f0 ¼

 p0 þ r0 ð3  u0 þ u20 Þ @p0 E1 ðp0 þ r0 Þ wer   C v 0 ðE1  p0  r0 Þ u0 C v 0 ðE1  p0  r0 Þ C v 0 u0 ð1  u0 Þ @t

gw @u0 K 0 @t



3  u0 @p0 u0 ð1  u0 Þ @x

ðA:10Þ

ðA:11Þ

Please cite this article in press as: M.A. Parron Vera et al., Analytical solution of coupled soil erosion and consolidation equations by asymptotic expansion approach, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.02.006

M.A. Parron Vera et al. / Applied Mathematical Modelling xxx (2014) xxx–xxx

g0 ¼

E1 ðp0 þ r0 Þwer C v 0 ðE1  p0  r0 Þ

ðr 1 Þnm ¼ a0 ðp1 Þnm  ðu1 Þnm  c0 Dt ðs1 Þnm ¼ ðp1 Þnm 

g Dt f0 Dt  Dt  ðu1 Þnmþ1  ðu1 Þnm  0 þ ðp1 Þnmþ1  2ðp1 Þnm þ ðp1 Þnm1 d0 Dx d0 d0 Dx2

13

ðA:12Þ ðA:13Þ ðA:14Þ

Appendix A. Supplementary material Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/ j.apm.2014.02.006. References [1] B.A. Kissi, M.A. Parron Vera, M.D. Rubio Cintas, P. Dubujet, A. Khamlichi, M. Bezzazi, L. El. Bakkali, Predicting initial erosion during the hole erosion test by using turbulent flow CFD simulation, Appl. Math. Model. 36 (2012) 3359–3370. [2] C.F. Wan, R. Fell, Assessing the potential of internal instability and suffusion in embankment dams and their foundations, J. Geotech. Geoenviron. Eng. 134 (2008) 401–407. [3] I. Vardoulakis, M. Stavropoulou, P. Papanastasiou, Hydromechanical aspects of sand production problem, Transp. Porous Media 22 (1996) 225–244. [4] I. Vardoulakis, Fluidisation in artesian flow conditions: hydromechanically stable granular media, Géotechnique 54 (2004) 117–130. [5] R.G. Wan, J. Wang, Analysis of sand production in unconsolidated oil sand using a coupled erosional-stress-deformation model, J. Can. Petrol. Technol. 43 (2004) 47–52. [6] M. Stavropoulou, P. Papanastasiou, I. Vardoulakis, Coupled wellbore erosion and stability analysis, Int. J. Numer. Anal. Meth. Geomech. 22 (1998) 749– 769. [7] E. Papamichos, Failure in rocks: hydro-mechanical coupling for erosion, Revue Française de Génie-Civil 8 (2004) 709–734. [8] R. Sakthivadivel, S. Irmay, A review of filtration theories, Technical Report of Hydraulic Engineering Laboratory, 15–4, University of California at Berkeley, 1966. [9] J. Lemaitre, How to use damage mechanics, Nucl. Eng. Des. 80 (1984) 233–245. [10] L. Kachanov, Introduction to Continuum Damage Mechanics, Martinus Nijhoff Publishers, Netherlands, 1986. [11] F. Yakhlef, A. Khamlichi, M. Bezzazi, M.A. Parron Vera, Ph. Dubujet. Solution of internal erosion equations by asymptotic expansion. MATEC Web of Conferences Volume 1, 2012, CSNDD 2012 – International Conference on Structural Nonlinear Dynamics and Diagnosis, Article number: 10007, 2012. http://dx.doi.org/10.1051/matecconf/20120110007. [12] P. Mira, M. Pastor, T. Li, X. Liu, A new stabilized enhanced strain element with equal order of interpolation for soil consolidation problems, Comput. Methods Appl. Mech. Eng. 192 (2003) 4257–4277. [13] K.H. Xie, K. Wang, Y.L. Wang, C.X. Li, Analytical solution for one-dimensional consolidation of clayey soils with a threshold gradient, Comput. Geotech. 37 (2010) 487–493. [14] K.H. Xie, C.J. Leo, Analytical solutions of one-dimensional large strain consolidation of saturated and homogeneous clays, Comput. Geotech. 31 (2004) 301–314. [15] A. Qin, D. Sun, Y. Tan, Analytical solution to one-dimensional consolidation in unsaturated soils under loading varying exponentially with time, Comput. Geotech. 37 (2010) 233–238. [16] R.P. Chen, W.H. Zhou, H.Z. Wang, Y.M. Chen, One-dimensional nonlinear consolidation of multi-layered soil by differential quadrature method, Comput. Geotech. 32 (2005) 358–369. [17] G.Y. Zheng, P.L. Li, C.Y. Zhao, Analysis of non-linear consolidation of soft clay by differential quadrature method, Appl. Clay Sci. 79 (2013) 2–7. [18] J. Wang, R.G. Wan, Computation of sand fluidization phenomena using stabilized finite elements, Finite Elem. Anal. Des. 40 (2004) 1681–1699. [19] E. Papamichos, I. Vardoulakis, Sand erosion with a porosity diffusion law, Comput. Geotech. 32 (2005) 47–58. [20] K. Terzaghi, P.B. Peck, G. Mesri, Soil Mechanics in Engineering Practice, Wiley-Interscience, New York, 1996. [21] J.G. Wang, G.R. Liu, P. Lin, Numerical analysis of Biot’s consolidation process by radial point interpolation method, Int. J. Solids Struct. 39 (2002) 1557– 1573. [22] P. Cornetti, M. Battaglio, Nonlinear consolidation of soil modelling and solution techniques, Math. Comput. Model. 20 (1994) 1–12.

Please cite this article in press as: M.A. Parron Vera et al., Analytical solution of coupled soil erosion and consolidation equations by asymptotic expansion approach, Appl. Math. Modell. (2014), http://dx.doi.org/10.1016/j.apm.2014.02.006