Empirical formulae for electric double-layer repulsion between two arbitrarily inclined clay particles

Xingyu Shng,Juming Lu,Linfei Kung,Chen Yng,Guoqing Zhou

aState Key Laboratory for Geomechanics and Deep Underground Engineering,China University of Mining and Technology,Xuzhou,221116,China

bSchool of Mechanics and Civil Engineering,China University of Mining and Technology,Xuzhou,221116,China

cChina Railway Engineering Equipment Group Co.,Ltd.,Zhengzhou,450016,China

Keywords:Clay particles Electric double-layer repulsion Numerical analysis Empirical formula

A B S T R A C T To understand the mesoscopic mechanism of clayey soil in view of macroscopic behavior,it is essential to quantitatively calculate the electric double-layer repulsion between arbitrarily inclined clay particles.However,suitable calculation methods with high efficiency and accuracy are still rare at present in literature.Based on a great number of numerical calculations of the repulsion between two inclined platy clay particles,explicit empirical formulae for estimating electric double-layer repulsion between clay particles are put forward.Comparison between the empirical solutions and corresponding numerical results shows that the proposed formulae have a reasonable accuracy,and application of the presented formula is easy and efficient.

1.Introduction

Surfaces of clay particles are usually charged negatively due to ismorphous substitutions,and the electric double-layer forms surrounding the clay particles when immersed in electrolyte solution.If these clay particles approach each other,double electric layers overlap,and then physical interactions including doublelayer repulsion between them arise(Jiang et al.,1989).It is significant to accurately calculate the double-layer repulsion for understanding the macroscopic mechanical behaviors of clayey soil such as compressibility(Tripathy and Schanz,2007;Bayesteh and Mirghasemi,2013),expansibility(Laird,2006;Katti et al.,2009;Wang et al.,2017),and lateral earth pressure coefficient(Shang et al.,2015a).In addition,calculation of double-layer repulsion is important for getting insight into sedimentary process(Lu et al.,2008).In fact,how to rapidly and accurately calculate the doublelayer repulsion between clay particles has been one of concerns in research field related to soil mechanics for many decades(Verwey and Overbeek,1948;Hogg et al.,1966;Anandarajah and Lu,1991;Lu and Anandarajah,1992;Anandrajah and Chen,1994;Shang et al.,2015b).method(FEM) (Anandarajah and Lu,1991;Shang et al.,2015b)or finite difference method(FDM)(Anandrajah and Chen,1994),the numerical electric potential field in discretized space domain,which is occupied by charged clay particle and electrolyte solution,is first solved,and then the potential in relevant areas is numerically integrated to obtain the double-layer repulsion between clay particles.Generally, fine mesh is indispensable to the target precision in numerical methods(Anandarajah and Lu,1991;Anandrajah and Chen,1994;Shang et al.,2015b).Alternatively,calculation of repulsion between inclined clay particles can be approximated by the summation of repulsions between a dozen of ideal parallel particles with smaller size(Anandrajah and Chen,1994).Discrete element method(DEM)has become a powerful tool that relates macroscopic to mesoscopic mechanical behaviors of clay,in which macroscopic clay medium is regarded as an assembly of discrete interacting charged plates(Anandarajah,1994;Lu et al.,2008;Katti et al.,2009;Bayesteh and Mirghasemi,2013).

The existing calculation methods of double-layer repulsion between charged clay particles include analytical,numerical and empirical methods.Among them,the analytical method is the most accurate and simple one,but it has limitation upon application due to the simplified hypotheses.For example,analytical solutions are available for only some special cases such as doublelayer repulsion between two infinite parallel charged plates(Verwey and Overbeek,1948),and two charged spheres(Hogg et al.,1966).Using numerical methods including finite element However,since a large number of contact forces have to be updated in every time step,it is almost impractical to use the above time-consuming numerical method in claydiscrete element code.

Based on numerous numerical results of double-layer repulsion,an explicit empirical formula of the repulsion between two inclined clay particles that contact with each other at one end was proposed by Lu and Anandarajah(1992),and the inputs included the mandatory geometric information such as length,thickness and orientation of clay particle,and some physical parameters including surface charge density,ion concentration and so on.Later,Anandarajah(1994)generalized the above empirical formula using an assumed distribution of electric potential to make his new formula applicable to two separate clay particles,but its rationality was not validated.Therefore,verification of empirical formula of the double-layer repulsion for the case of two arbitrarily inclined clay particles with reasonable precision is still unavailable up to now.In addition,the above empirical formulae have the limitation of low computational efficiency because their application still involves finding the solution of the potential on the mid-plane between two infinite parallel charged plates,which needs to solve nonlinear integral equation(Katti et al.,2009).Consequently,it is necessary to put forward a new empirical formula of the repulsion between two arbitrarily inclined clay particles with high computational efficiency as well as reasonable accuracy.

The aim of this study is to obtain a reasonable and efficient formula applicable to two arbitrary charged clay particles.The paper is organized as follows.First,the existing calculation methods of electric double-layer repulsion are briefly reviewed.Then numerical calculations of the repulsion between two clay particles are carried out based on the calculating procedure presented by Shang et al.(2015b).Finally,based on the numerical results,an empirical formula that can reasonably and efficiently estimate the double-layer repulsion between two arbitrary charged clay particles is proposed.

2.Review of the existing calculation methods for doublelayer repulsion

In a two-dimensional(2D)discrete element analysis of clayey soil,clay particles are usually treated as charged plates,as shown in Fig.1(Anandarajah,1994;Katti et al.,2009;Bayesteh and Mirghasemi,2013).For two infinite parallel clay particles presented in Fig.1a,the dimensionless double-layer repulsive stress can be estimated by(Verwey and Overbeek,1948):

whereφdis the dimensionless potential at the mid-plane as shown in Fig.1a,andris the repulsive force per unit length of the particle.However,due to the lack of an analytical solution to the mid-plane potential which is applicable to wider ranges of inter-particle distance and surface potential in the past(Shang et al.,2015c),numerical calculation is still needed in application of Eq.(1).As an example,the above mid-plane potential can be obtained by numerically solving the nonlinear integral equation below(Katti et al.,2009):

Fig.1.Calculating models for repulsive force between two platy clay particles.

Fig.1b shows the most general 2D model of two arbitrary clay particles which are finitely long,non-parallel as well as asymmetrical.However,it is reasonable to assume that only the interaction between the parts of particle surfaces directly facing each other makes significant contribution to repulsive force due to the localized nature of the repulsive force(Anandarajah,1994;Bayesteh and Mirghasemi,2013).In other words,the model illustrated in Fig.1b can be reduced to that of two symmetrical but non-parallel charged particles with finite length,as presented in Fig.1c.For this symmetrical model,Anandarajah and Lu(1991)indicated that the corresponding dimensionless repulsive force can be estimated by the following equation:

whereΓmdenotes the mid-plane shown in Fig.1c,andRis the dimensionless double-layer repulsive force acting along the direction perpendicular to the mid-plane.However,the error of the repulsion estimated by Eq.(3)can be large.A modified equation has been presented by Shang et al.(2015b)as follows:

whereζdenotes the unit vector along the axis parallel to the above mid-plane.For the case of two non-parallel clay particles,the distribution of potential along the mid-plane(Eqs.(3)and(4))can only be obtained by solving Poisson-Boltzmann(PB)equation numerically(Anandarajah and Lu,1991;Anandrajah and Chen,1994;Shang et al.,2015b).As noted earlier,this is almost impracticable in actual discrete element analysis of clay.A feasible approximation method pioneered by Derjaguin(1934)is used to divide the inclined clay particles shown in Fig.1c into smaller particles,and the double-layer repulsion between smaller nonparallel particles can be approximated by that between parallel particles with the same size(Anandrajah and Chen,1994).Consequently,calculation of repulsive force between one pair of inclined particles is converted into that between dozens of pairs of parallel particles.It should be noted that estimation of repulsive force between parallel clay particles still involves numerically solving nonlinear equation such as Eq.(2).Although a new closed-form solution of mid-plane potential between two infinite charged plates(Shang et al.,2015c)eliminates the need to solve Eq.(2)numerically,the PB equation has to be solved in order to obtain the gradient of potential in Eq.(4).As a result,the computation cost of clay discrete element analysis using the above approximation method will be expensive.

In addition,for the above approximation method,it can be deduced that the reasonable number of subdivision of charged non-parallel particles with large inter-particle angle may be a few dozen times that of particles with small angle.Therefore,it is difficult to predict a reasonable subdivision number prior to the calculation of double-layer repulsion.

To simplify the calculation of double-layer repulsion between clay particles further,after solving PB equation by means of FEM,Lu and Anandarajah(1992)obtained a great number of numerical repulsive forces between two non-parallel symmetrical clay particles contacting each other at one end,as shown in Fig.1d,according to Eq.(3).Based on this,they put forward the following empirical formulae of repulsive forceRand its acting positionl,i.e.the distance from the contact point to the location at which the resultant repulsive force is applied to one particle:

whereLis the length of clay particles;θis the angle between two particles;f1andf2are the fitted explicit functions of the length,potential and angle;andR∞denotes the double-layer repulsive stress between two parallel charged particles,which is the product of the repulsive stress calculated by Eq.(1)and the length of the particle.

Because the arrangement of clay particles in actual clay discrete element analysis is primarily random,the model shown in Fig.1d is not representative.Therefore,application of Eqs.(5)and(6)is limited.For the general model presented in Fig.1c,Anandarajah(1994)assumed that the potential at the mid-plane between two infinite clay particlesφdis the function of particle surface potential φ0and inter-particle distanceD,i.e.φd= φ0e-D/2，and that the hyperbolic cosine function of potential at the mid-plane can be approximated by its first-order Taylor series expansion,i.e.coshφd≈1+φ2d/2.Based on these,the empirical formulae of repulsive forceRand the distancelfor the model shown in Fig.1c are given as follows:

whereDis the minimum distance between two symmetrical particles;andR0andl0are determined according to Eqs.(5)and(6),respectively.It shows that the assumptions used to deduce Eqs.(7)and(8)are likely to limit the scope of their application.In addition,the rationality of these new formulae was not verified(Anandarajah,1994).For example,using numerical integration after solving PB equation by means of FEM,the numerical dimensionless double-layer repulsionRand the above-mentioned distancelare obtained for the model shown in Fig.1c with given parameters,i.e.inter-particle angle of 30°,dimensionless particle length of 0.5,and dimensionless surface potentials of 6 and 8.Fig.2 presents the numerical results by varying the minimum inter-particle distance as well as the corresponding empirical results determined according to Eqs.(7)and(8).

It can be seen from Fig.2a and b that the empirical repulsive force estimated by Eq.(7)agrees well with the corresponding numerical results when the distance between particles is large,but in a global sense the former is significantly larger than the latter.This is due to the above-mentioned Taylor series approximation of midplane potential,which only has reasonable accuracy when the potential is small.In the case of large inter-particle distance,the mid-plane potential is small because the influence of particle surface potential on the mid-plane potential is weak.The empirically estimated position of repulsive force is independent of interparticle distance,as shown in Fig.2c and d,according to Eq.(8),whilst the corresponding numerical position is not the case.This is due to the oversimplified assumption used in Eq.(8)that the interparticle distance has no influence on the acting position.Therefore,it is necessary to explore new empirical formulae to calculate repulsive force and its position accurately for the general situation as shown in Fig.1c.

3.New formulae of double-layer repulsive force and verification

To facilitate the following calculation,all the relevant parameters are dimensionless following the previous studies(Anandarajah and Lu,1991;Lu and Anandarajah,1992;Anandrajah and Chen,1994;Shang et al.,2015b).Two key parameters,dimensionless lengthLand surface potentialφ0,are given below:

Fig.2.Comparison of numerical and empirical solutions of value and position of electric repulsive force in the case ofθ=30°and L=0.5.

whereL0is the actual length of the particle,nis the electrolyte concentration,vis the valence of cation in the liquid phase,eis the elementary electric charge,ε is the dielectric constant of the liquid phase,κis the Boltzmann’s constant,Tis the absolute temperature of the system,CECis the cation exchange capacity for clay,FCis the Faraday’s constant,andSis the specific area of clay.The ranges of the dimensionless variables adopted in the numerical analyses in Lu and Anandarajah(1992)include surface potential of 4-8,angle between particles from 0°to 120°,constant minimum distance of 0,and length of particle from 1 to 8.In order to obtain more targeted formulae,referring to the typical values of particle length,specific area,CECand electrolyte concentration of kaolinite,illite and montmorillonite clay(Mitchell and Soga,2005),and using Eqs.(9)and(10),the ranges of the parameters are enlarged,such as the dimensionless surface potential of 4-10,the minimum distance of 0.1-10,and the particle length of 0.5-10.The details of the parameters used in the following numerical analyses are shown in Table 1.Accordingly,1512 numerical simulations are carried out.

The numerical calculation of repulsion forces needs the potential distribution which is obtained by solving PB equation using thenonlinear FEM(more details can be seen in Shang et al.(2015b)).Fig.3 shows some typical potential contours calculated.Once the above potential field is obtained,the repulsive force can be estimated using Eq.(4).

Table 1 Parameters for numerical calculation of electric repulsive force between charged clay particles.

The numerical results obtained are fitted using MATLAB nonlinear regression.With reference to the function form presented in Lu and Anandarajah(1992),the following empirical formulae are put forward to estimate the repulsive force and its acting position between the symmetrical charged platy particles shown in Fig.4:

Fig.3.Typical numerical potential contours.

Fig.4.Schematic of the numerical model.

where

It is shown that the convergent squared correlation coefficients,which can be used to illustrate how well the above empirical equations represent the numerical data overall,for Eqs.(11)and(12)are 0.9657 and 0.9612,respectively.

Comparing with the previous empirical formulae(Eqs.(5)and(6)),the effect of inter-particle distance has been considered in the presented formulae(Eqs.(11)and(12)).For example,there is a term related to the inter-particle distanceDin Eq.(17).TheR≈in Eq.(11)denotes the repulsion between two parallel particles with finite length.It includes the effect of the distanceDas well as the effect of potential gradient,as shown in Eq.(4).

It should be noted that determination ofR∞in Eq.(5)involves evaluation of the potential at the mid-plane between two infinite charged plates,which needs to solve the nonlinear integral equation(Eq.(2))numerically(Lu and Anandarajah,1992;Katti et al.,2009).In order to improve the computation efficiency,an explicit formula forR≈is proposed,which is applicable to a wide range of dimensionless inter-particle distance from 0.1 to 10 and dimensionless surface potential from 4 to 10.This formula is given as follows:

Eq.(19)is obtained by means of nonlinear curve fitting method similar to that used by Shang et al.(2015c)and has compact form and high accuracy with the squared correlation coefficient larger than 0.992.The combination of Eqs.(11),(12)and(19)gives explicit empirical formulae applicable to evaluating the inter-particle repulsive force and its acting position for the case shown in Figs.4 and 1c,thereby the general case of Fig.1b according to the interpretation presented in Section 2.It can be seen that any numerical process is not needed by using the above empirical formulae.

Fig.5.Comparison of the numerical and presented empirical solutions of value and position of electric repulsive force.

Fig.5 shows the comparison of numerical and empirical solutions for three kinds of typical cases:(i)L=4,and θ=15°;(ii)L=6,and φ0=6;and(iii)φ0=8,and θ=90°.It should be noted once again that all the parameters involved in Fig.5 are dimensionless,and all the calculated results correspond to the symmetrical model shown in Fig.4.It can be seen that the calculated repulsive force decreases but its acting position increases with the increasing interparticle distance.Both of them tend to be constant as the interparticle distance increases.

In general,the agreement between numerical and empirical results is quite good in some cases but to some extent less in others.The left graph in Fig.5a shows that the repulsive force increases with increasing surface potential and decreasing interparticle distance,and the empirical solution is consistent with the corresponding numerical solution except that the former is slightly smaller than the latter when the inter-particle distance falls into the range of 0.5-1.The right graph in Fig.5a shows that the acting location increases with decreasing surface potential and increasing inter-particle distance,and the empirical solution is larger than the numerical one when the inter-particle distanceDis smaller than 1,while the reverse is the case forD＞2.It can be seen from the left graph in Fig.5b that the repulsive force decreases with increasing inter-particle distance and angle,and good agreement exists except whenDis in the range of 0.5-1 andθis nonzero,where the empirical solution is slightly smaller than the numerical one.The right graph in Fig.5b shows that the acting location increases with increasing inter-particle distance and decreasing inter-particle angle,and the empirical solution is slightly smaller than the numerical one whenDis in the range of 0.5-1 while the reverse is the case forD＞ 2 and 15°＜θ＜ 30°.The left graph in Fig.5c illustrates that the variation of repulsive force with increasing particle length is relatively complicated when θ =90°,and good agreement exists except when the inter-particle distance is in the range of 0.5-1,where the empirical solution is slightly smaller than the numerical one;while the reverse is the case forD=0.1 andL＞2.The right graph in Fig.5c shows that the acting location increases with increasing particle length,and the empirical solution is larger than the numerical one when the interparticle distance is smaller than 1,while the reverse is the case forL＞2 andD＞4.It should be noted that the above differences between empirical calculations and numerical solutions are the results of curve fitting method as the curve fitting method is designed to produce an approximating function in the sense of the minimum total deviation.

The discrepancy between the empirical and numerical solutions for all other cases,which are not presented herein for the sake of brevity,was found to be in the same order as shown in Fig.5.Based on the above analyses and considering that the interparticle repulsive force and the acting location for the general model presented in Fig.1b can be approximated by those for the symmetrical model shown in Fig.1c,the empirical formulae presented herein,i.e.Eqs.(11),(12)and(19),have the advantages of extensive application range,high computational efficiency and reasonable accuracy.It can be expected that the present empirical solution is particularly useful for the discrete element analysis of clay aiming at investigating the macro-micro mechanical links of clayey soil.

4.Conclusions

According to the new calculation method for the repulsive force between two inclined charged plates and the compact explicit solution of potential at the middle of the two parallel charged plates recently proposed by the authors,a great number of numerical simulations of repulsive force and the acting location are carried out for two charged clay particles with various particle surface potentials,inter-particle angles and inter-particle distances.Based on the numerical results,new empirical formulae of repulsive force and its acting location for two arbitrarily inclined platy particles are put forward.Using the presented formulae,the repulsion between charged clay particles can be explicitly estimated with a reasonable accuracy using the particle surface potential,pore fluid physical properties,relevant geometric information,and so on.The proposed formulae in this paper,as a basic calculation tool,will be useful for quantitatively exploiting and understanding the macroscopic mechanical behavior of clayey soil from the scale of clayey particle.

Conflicts of interest

The authors wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgements

The authors gratefully acknowledge the financial support from“The Fundamental Research Funds for the Central Universities”(Grant No.2017XKQY052).