On the Calculation of Van der Waals Force between Clay Particles

: Complex physical–chemical interactions between clay particles including the van der Waals force control the macroscopic behaviors of clay. The calculation of van der Waals force is essential in the discrete element method (DEM) that has been widely used to enhance the understanding of the behavior of clay. Due to its high computational e ﬃ ciency, a plate-wall model developed in the literature was adopted to obtain the van der Waals force between the neighboring clay particles approximately in the published DEM simulations. However, di ﬀ erent choices of the ideal wall result in di ﬀ erent magnitudes as well as directions of total van der Waals forces. To investigate the e ﬀ ect, a new rigorous plate-plate model was put forward and solved using an optimized Cotes integration method. Based on the comparison between the calculated results based on plate-wall and plate-plate model, the above e ﬀ ect was analyzed for the cases of face-face, edge-edge and edge-face contact types. Then, necessary advice for the reasonable use of the ideal-wall model was given accordingly. rectangular Cartesian coordinate systems O 1 TSW and O 2 XYZ moving with plate1 and plate2, respectively, are deﬁned, was simpliﬁed into that between a ﬁnite plate and an inﬁnite plate which is referenced to as a wall as shown in Figure 2. According to the ideal plate-wall model shown in Figure 2 where Z-axes and W-axes are assumed to be parallel to each other, a close-form formula for van der Waals force was obtained [17], which has been used widely in the later DEM simulations of clays [13–15,18,19] due to the high e ﬃ ciency of close form solution. special case of the present plate-plate model. X.-Y.S.; Validation, X.-Y.S. and K.Z.; Formal Analysis, X.-Y.S. and K.Z.; Investigation, Q.-Y.Z.; Resources, G.-Q.Z.; Data Curation, K.Z.; Writing—Original Draft Preparation, X.-Y.S., K.Z. and W.-X.Q.; Writing—Review & Editing, X.-Y.S. authors


Introduction
van der Waals forces act between all atoms and molecules [1] and play a role in many important phenomena such as the aggregate stability of soil [2], adsorption behavior of clay [3,4], and macroscopic deformation of clay [5]. The well-known DLVO theory (Deryaguin-Landau-Verwey-Overbeek) considering the contribution of van der Waals force has been used extensively to predict the stability of colloidal suspensions including clay-water systems [6]. In fact, there exists a tight link between the macro mechanical behaviors of clays and the interactive forces between the clay particles [5,6]. Exploiting these links is necessary for the deep understanding of the mechanism of macro behaviors of clay. At the nanoscale, molecular dynamics simulations have been adopted successfully to exploit the physical properties of clay-water system [7][8][9]. In these simulations, van der Waals interactions between atoms are usually taken into account by means of the model of force field, e.g., CLAYFF model (CLAY Force Field) [8]. Based on molecular dynamics, it has been shown that the parameter related to van der Waals energy has profound influence on the macroscopic behaviors such as capillary forces on clay particle [7]. However, due to the limited computational resources, most of the molecular simulations of clay-water focused on idealized infinite clay lamellae [9]. To simulate multi-aggregate systems of clay, a proper meso-scale interaction potential of the clay aggregation, which is the essential part of coarse-graining upscaling method, has been developed based on full atomistic simulations [10]. The other coarse-graining strategy by which each clay particle is composed of hundreds of spherical sub-particles physically representing a great deal of atoms has also been proposed to simulate the atoms has also been proposed to simulate the macroscopic clay behavior [11]. In the field of soil mechanics, discrete element method (DEM) which simulates the macroscopic behavior of soil by considering the inter-particle forces quantitatively has been used to reveal the abovementioned links for nearly four decades [12]. DEM simulations of clays have also been developed to achieve a quantitative multiscale interpretation of macro mechanical behavior of clays [13,14]. In such simulations, it is essential to calculate DLVO forces including van der Waals as well as electrostatic forces accurately and efficiently. Although DLVO theory fails to predict short range repulsion which results from hydration force [1,6], it was satisfactorily used in DEM simulations of clays to reproduce the typical mechanical behavior of clayey soil, in which both electrostatic and van der Waals forces were treated as constants once the distance between neighboring clay particles become small enough [13,15]. It is well known that both van der Waals and electrostatic forces decrease with the increasing inter-particle distance [1,6], yet the relative magnitude of these two forces depends on many factors such as the relative orientation of the interacting particles, ion valence, concentration of the pore electrolyte and so on [13,14,16]. The calculation of van der Waals forces between two finite clay particles is the concern of this study.
Based on the well-known London theory for attractive energy between two-unit molecules, the theories for attractive energy between two macro-bodies have been developed significantly by Hamaker, de Boer, Casimir and Polder, and Lifshitz [1,6,17]. Following these classical theories, Anandarajah and Chen put forward a combined theory that was referred to as Hamaker-de Boer-Lifshitz theory in order to account for the retardation effect in a condensed system as well as the effect of geometry of macro-bodies [17]. To calculate the van der Waals forces between neighboring clay particles reasonably and efficiently, the geometry of the clay particle has been assumed to be a cubic plate. Furthermore, the calculation of van der Waals forces between two arbitrary plates as shown in Figure 1, in which two rectangular Cartesian coordinate systems O1TSW and O2XYZ moving with plate1 and plate2, respectively, are defined, was simplified into that between a finite plate and an infinite plate which is referenced to as a wall as shown in Figure 2. According to the ideal plate-wall model shown in Figure 2 where Z-axes and W-axes are assumed to be parallel to each other, a closeform formula for van der Waals force was obtained [17], which has been used widely in the later DEM simulations of clays [13][14][15]18,19] due to the high efficiency of close form solution.  The analytical formula of the van der Waals forces for the above ideal plate-wall model is [17] = 6 sin2 (−1) 4 − 1 − 12( + ) + = , = + , = + + sin , = + sin , = + + cos , = + sin , = + sin + cos , = + + sin + cos = + sin + cos , = + + sin + cos (1) where A is the Hamaker constant 1.59 × 10 −20 J, c is a constant dependent on the characteristic length of interaction (49.363 nm is taken throughout the present study following [17]), w1, d1 and L1 are the width, thickness, and length of plate1, respectively, d2 is the thickness of the wall (i.e., plate2 in Figure  2), and d and are the minimum distance and angle between plate1 and the wall, respectively. It was also shown that the force Fpw in Equation (1) is in the direction of Y-axes, i.e., perpendicular to the selected wall shown in Figure 2 [17]. Anandarajah and Chen further gave the similar formulas of the coordinates of action of Fpw on plate1, which are not presented here for space reasons.
Chen adopted a rigorous method of adding van der Waals forces between all the pairs of molecules over these two plates [16]. By this method, an expression of a multi integral instead of analytical formula like Equation (1) was obtained. To simplify the calculation, Chen assumed that Zaxes and W-axes shown in Figure 3 were parallel to each other. As a result, the relative inclination between these two plates can be determined by one angle α. The expression of van der Waals force acting on plate2 in Figure 3 where V1, V2 denote the volume of plate1 and plate2, r is the distance between an arbitrary molecule A(x,y,z) in plate2 and B(t,s,w) in plate1, denotes the direction of the vector AB with respect to The analytical formula of the van der Waals forces for the above ideal plate-wall model is [17] F pw = Aw 1 6π sin 2α where A is the Hamaker constant 1.59 × 10 −20 J, c is a constant dependent on the characteristic length of interaction (49.363 nm is taken throughout the present study following [17]), w 1 , d 1 and L 1 are the width, thickness, and length of plate1, respectively, d 2 is the thickness of the wall (i.e., plate2 in Figure 2), and d and α are the minimum distance and angle between plate1 and the wall, respectively. It was also shown that the force F pw in Equation (1) is in the direction of Y-axes, i.e., perpendicular to the selected wall shown in Figure 2 [17]. Anandarajah and Chen further gave the similar formulas of the coordinates of action of F pw on plate1, which are not presented here for space reasons. Chen adopted a rigorous method of adding van der Waals forces between all the pairs of molecules over these two plates [16]. By this method, an expression of a multi integral instead of analytical formula like Equation (1) was obtained. To simplify the calculation, Chen assumed that Z-axes and W-axes shown in Figure 3 were parallel to each other. As a result, the relative inclination between these two plates can be determined by one angle α. The expression of van der Waals force acting on plate2 in Figure 3 are where V 1 , V 2 denote the volume of plate1 and plate2, r is the distance between an arbitrary molecule A(x,y,z) in plate2 and B(t,s,w) in plate1, β i denotes the direction of the vector AB with respect to O 2 XYZ coordinate system, D is the minimum distance between the two plates as shown in Figure 3, d 2 is the Minerals 2020, 10, 993 4 of 15 thickness of plate2, and l x and l z are the coordinates of O 1 in the O 2 XYZ system. The expressions of the locations of total force on the two plates C 1 and C 2 as shown in Figure 3 were also obtained [16], which are not present here for space reasons.
Minerals 2020, 10, x FOR PEER REVIEW 4 of 15 O2XYZ coordinate system, D is the minimum distance between the two plates as shown in Figure 3, d2 is the thickness of plate2, and lx and lz are the coordinates of O1 in the O2XYZ system. The expressions of the locations of total force on the two plates C1 and C2 as shown in Figure 3 were also obtained [16], which are not present here for space reasons. The calculation related to Equation (2) is too time consuming [16] to be used in practical DEM simulation. To reduce computation effort, Chen put forward an approximation method of superposition [16]. However, the method of superposition involves the solutions of nine equations, among which there is a nonlinear one. The reason why the abovementioned analytical solution based on the ideal plate-wall model instead of the approximation method is widely used in DEM simulations is most likely its simplicity.
However, there is a dilemma of how to choose a proper ideal wall in the two plates shown in Figure 1. It seems that either of the two plates can be chosen as the ideal wall. Nevertheless, it can be deduced that different choices result in different magnitudes as well as directions of total van der Waals forces between neighboring clay particles. Therefore, the present study aims to give reasonable advice for the proper use of the ideal plate-wall model by comparing its solutions and those of a more rigorous plate-plate model.
The remainder of this paper is organized as follows. Firstly, the integral expression of a rigorous plate-plate model in which two plates are arbitrarily located in three-dimensional space is given and its numerical solution is verified. Then, the effects of choice of the ideal wall in two plates on the calculated van der Waals force are investigated systematically for various contact types including face-face, face-edge contact, and so on. Finally, according to the corresponding analysis, necessary advice for the reasonable use of the ideal-wall model is given.

van der Waals Force for the New Rigorous Plate-Plate Model
For the model of two arbitrarily located plates as shown in Figure 3, the coordinates of any molecule B(xB, yB, zB) in plate1 with respect to the O2XYZ system are related to the counterpart B(t,s,w) in the O1TSW system by a rotational matrix instead of only one angle α as described in Chen's plateplate model. The calculation related to Equation (2) is too time consuming [16] to be used in practical DEM simulation. To reduce computation effort, Chen put forward an approximation method of superposition [16]. However, the method of superposition involves the solutions of nine equations, among which there is a nonlinear one. The reason why the abovementioned analytical solution based on the ideal plate-wall model instead of the approximation method is widely used in DEM simulations is most likely its simplicity.
However, there is a dilemma of how to choose a proper ideal wall in the two plates shown in Figure 1. It seems that either of the two plates can be chosen as the ideal wall. Nevertheless, it can be deduced that different choices result in different magnitudes as well as directions of total van der Waals forces between neighboring clay particles. Therefore, the present study aims to give reasonable advice for the proper use of the ideal plate-wall model by comparing its solutions and those of a more rigorous plate-plate model.
The remainder of this paper is organized as follows. Firstly, the integral expression of a rigorous plate-plate model in which two plates are arbitrarily located in three-dimensional space is given and its numerical solution is verified. Then, the effects of choice of the ideal wall in two plates on the calculated van der Waals force are investigated systematically for various contact types including face-face, face-edge contact, and so on. Finally, according to the corresponding analysis, necessary advice for the reasonable use of the ideal-wall model is given.

Van der Waals Force for the New Rigorous Plate-Plate Model
For the model of two arbitrarily located plates as shown in Figure 3, the coordinates of any molecule B(x B , y B , z B ) in plate1 with respect to the O 2 XYZ system are related to the counterpart Minerals 2020, 10, 993 5 of 15 B(t,s,w) in the O 1 TSW system by a rotational matrix instead of only one angle α as described in Chen's plate-plate model.
r xt r xs r xw r yt r ys r yw r zt r zs r zw is the coordinates of O 1 with respect to O 2 XYZ system, R 12 is the rotational matrix, and the element of R 12 , r ij denotes the directional cosine between i axes and j axes. Accordingly, the distance between an arbitrary point A(x,y,z) in plate2 and B(t,s,w) in plate1 is and the direction of the vector AB with respect to O 2 XYZ coordinate system is The expression of van der Waals force acting on plate2 in the present model can be obtained by substituting Equations (4) and (5) into Equation (2). As can be seen, Chen's plate-plate model can be regarded as a special case of the present plate-plate model.

Numerical Solution of van der Waals Force of the New Model
It can be noticed that multiple integral is involved in the numerical solution of Equation (2). Due to the highly nonlinear integrand which is a function of Equations (4) and (5), a reasonable numerical procedure has to be elected to obtain the accurate solutions. It should be noted that Monte Carlo [20] and quasi-Monte Carlo [21] methods have been adopted by the author to address this problem. However, plenty of calculations indicated that the integral values are sensitive to the selection of random sample points, and the solutions are often unable to converge with the increasing number of sample points. An optimized Cotes integral method [22] is elected on a trial and error basis.
According to the used integral method [22], the numerical precision of multiple integrals over a six-dimensional cubic domain can be guaranteed by continuous division of the integral interval in the following way: the interval [a i ,b i ] is divided into 4m i equal subintervals, then each subinterval can be divided into five equal parts repeatedly until the desired precision has been reached. The corresponding numerical integral formula is Minerals 2020, 10, 993

Verification of the Obtained Numerical Solution
An example of the present model is taken to verify the numerical solution. In this example, the two plates have the same length and width of 0.1 µm and the same thickness of 0.01 µm, they are parallel to each other, and the distance between them is 0.001 µm. The example can equally be described in terms of necessary quantities introduced previously, This example is solved numerically in both frameworks of Chen's model and the present model using the abovementioned Cotes integral method. The integral intervals along the length and width of the both plates are divided into 4m subintervals, and those along the thickness of the both plates are divided into 0.4m subintervals. Nine values of m ranging from 5 to 35 are taken. At the same time, the example is also solved using the ideal plate-wall formula Equation (1). We used F ppy and F pwy to denote the calculated force based on the plate-plate and plate-wall model, respectively. All the forces are along the Y-axes since the two plates are parallel to each other. Figure 4 presents the numerical results. It can be seen that the solutions based on Chen's plate-plate model and those on the present model are exactly the same since the former is a special case of the latter, and the integral converges when the number of division m exceeds 27.5. The converged ratio of F ppy to F pwy is about 0.3, which means that the ideal plate-wall model overestimates the total van der Waals force between neighboring clay particles.

Verification of the Obtained Numerical Solution
An example of the present model is taken to verify the numerical solution. In this example, the two plates have the same length and width of 0.1 μm and the same thickness of 0.01 μm, they are parallel to each other, and the distance between them is 0.001 μm. The example can equally be described in terms of necessary quantities introduced previously, d1 = d2 = 0.01, D = 0.001, l1 = l2= w1 = w2 = 0.1, lx = lz = 0, α = 0 (the unit μm is omitted for brevity).
This example is solved numerically in both frameworks of Chen's model and the present model using the abovementioned Cotes integral method. The integral intervals along the length and width of the both plates are divided into 4m subintervals, and those along the thickness of the both plates are divided into 0.4m subintervals. Nine values of m ranging from 5 to 35 are taken. At the same time, the example is also solved using the ideal plate-wall formula Equation (1). We used Fppy and Fpwy to denote the calculated force based on the plate-plate and plate-wall model, respectively. All the forces are along the Y-axes since the two plates are parallel to each other. Figure 4 presents the numerical results. It can be seen that the solutions based on Chen's plateplate model and those on the present model are exactly the same since the former is a special case of the latter, and the integral converges when the number of division m exceeds 27.5. The converged ratio of Fppy to Fpwy is about 0.3, which means that the ideal plate-wall model overestimates the total van der Waals force between neighboring clay particles. The possible reason for the overestimation is that the ideal plate-wall model treats one of the two finite plates as an infinite plate. To verify the speculation, more calculations in which the size of plate1 is fixed and that of plated2 increases gradually while the two plates are parallel to each other with the line linking the centers of the two plates along the normal of the plate were carried out. Table  1 presents the calculated results. It can be seen that the calculated force based on the present model increased with the increasing size of plate2. The ratio of Fppy to Fpwy is about 0.93 when the area of plate2 is 10 6 times that of plate1.  The possible reason for the overestimation is that the ideal plate-wall model treats one of the two finite plates as an infinite plate. To verify the speculation, more calculations in which the size of plate1 is fixed and that of plated2 increases gradually while the two plates are parallel to each other with the line linking the centers of the two plates along the normal of the plate were carried out. Table 1 presents the calculated results. It can be seen that the calculated force based on the present model increased with the increasing size of plate2. The ratio of F ppy to F pwy is about 0.93 when the area of plate2 is 10 6 times that of plate1.

Effect of the Choice of Ideal Wall
Three basic contact types are shown schematically in Figure 5, among which two plates are parallel to each other in the face-face type. The effect of the choice of the ideal wall will be investigated separately in this section.

Face-Face Type
For this type, we only focus on the effect of the choice of ideal wall on the magnitude of the calculated van der Waals force since the force is always along the direction perpendicular to the face, i.e., Y-axis shown in Figure 6 in which two specific cases of face-face type are illustrated. All the endpoints of particle1 are projected on the particle2 in case (a), while only part of them are projected on the partilce2 in case (b). The geometric parameters related the two plates are: l1 = w1 = 0.1, l2 = w2 = 0.2, d1 = d2= 0.01. lx and lz are shown in Figure 6. The distance between the two plates, D, increases from 1nm to 10 nm.

Face-Face Type
For this type, we only focus on the effect of the choice of ideal wall on the magnitude of the calculated van der Waals force since the force is always along the direction perpendicular to the face, i.e., Y-axis shown in Figure 6 in which two specific cases of face-face type are illustrated. All the endpoints of particle1 are projected on the particle2 in case (a), while only part of them are projected on the partilce2 in case (b). The geometric parameters related the two plates are: l 1 = w 1 = 0.1, l 2 = w 2 = 0.2, d 1 = d 2 = 0.01. l x and l z are shown in Figure 6. The distance between the two plates, D, increases from 1 nm to 10 nm.
i.e., Y-axis shown in Figure 6 in which two specific cases of face-face type are illustrated. All the endpoints of particle1 are projected on the particle2 in case (a), while only part of them are projected on the partilce2 in case (b). The geometric parameters related the two plates are: l1 = w1 = 0.1, l2 = w2 = 0.2, d1 = d2= 0.01. lx and lz are shown in Figure 6. The distance between the two plates, D, increases from 1nm to 10 nm.  Figure 7 presents the calculated ratios of Fpwy to Fppy of the two specific cases shown in Figure 6. For each case, two different choices of the ideal wall were made. It can be seen that the calculated Fpwy corresponding to the choice of treating plate1 as the ideal wall is larger than that of plate2 as the ideal wall, and the ratio of Fpwy to Fppy corresponding to case (a) is larger than that to case (b). The reason why the calculated Fpwy corresponding to plate1 being the wall is much larger is that plate2 has a larger area. Therefore, it seems like the choice of the ideal wall has a remarkable effect on the calculated force depending on the geometric configuration of the two plates. To reduce the overestimation of the ideal plate-wall model as much as possible, it is reasonable to choose the plate with a smaller area as the ideal wall. For the specific cases shown in Figure 6, or two-dimensional cases, the following procedure is put forward. After choosing one plate i as the ideal wall, the other plate j of the two plates is projected on the chosen plate. Then the projected length and width are regarded as the new counterparts of the plate j. The new sizes are denoted by the variables with an apostrophe shown in Figure 8 and Figure 9 that present the projection for the two specific cases, as plate1 and plate2 are chosen as the ideal wall, respectively.  Figure 7 presents the calculated ratios of F pwy to F ppy of the two specific cases shown in Figure 6. For each case, two different choices of the ideal wall were made. It can be seen that the calculated F pwy corresponding to the choice of treating plate1 as the ideal wall is larger than that of plate2 as the ideal wall, and the ratio of F pwy to F ppy corresponding to case (a) is larger than that to case (b). The reason why the calculated F pwy corresponding to plate1 being the wall is much larger is that plate2 has a larger area. Therefore, it seems like the choice of the ideal wall has a remarkable effect on the calculated force depending on the geometric configuration of the two plates.  Figure 7 presents the calculated ratios of Fpwy to Fppy of the two specific cases shown in Figure 6. For each case, two different choices of the ideal wall were made. It can be seen that the calculated Fpwy corresponding to the choice of treating plate1 as the ideal wall is larger than that of plate2 as the ideal wall, and the ratio of Fpwy to Fppy corresponding to case (a) is larger than that to case (b). The reason why the calculated Fpwy corresponding to plate1 being the wall is much larger is that plate2 has a larger area. Therefore, it seems like the choice of the ideal wall has a remarkable effect on the calculated force depending on the geometric configuration of the two plates. To reduce the overestimation of the ideal plate-wall model as much as possible, it is reasonable to choose the plate with a smaller area as the ideal wall. For the specific cases shown in Figure 6, or two-dimensional cases, the following procedure is put forward. After choosing one plate i as the ideal wall, the other plate j of the two plates is projected on the chosen plate. Then the projected length and width are regarded as the new counterparts of the plate j. The new sizes are denoted by the variables with an apostrophe shown in Figure 8 and Figure 9 that present the projection for the two specific cases, as plate1 and plate2 are chosen as the ideal wall, respectively.  To reduce the overestimation of the ideal plate-wall model as much as possible, it is reasonable to choose the plate with a smaller area as the ideal wall. For the specific cases shown in Figure 6, or two-dimensional cases, the following procedure is put forward. After choosing one plate i as the ideal wall, the other plate j of the two plates is projected on the chosen plate. Then the projected length and width are regarded as the new counterparts of the plate j. The new sizes are denoted by the variables with an apostrophe shown in Figures 8 and 9 that present the projection for the two specific cases, as plate1 and plate2 are chosen as the ideal wall, respectively. It can be seen from Figure 10 that the calculated forces are relatively small and the same no matter which plate is chosen as the ideal wall. Therefore, for the specific three-dimensional cases in which the projection of one plate on the other plate is a rectangle or two-dimensional case, the proposed projection procedure can eliminate the adverse effect of the choice of the ideal wall. To reduce the overestimation of the ideal plate-wall model as much as possible, it is reasonable to choose the plate with a smaller area as the ideal wall. For the specific cases shown in Figure 6, or two-dimensional cases, the following procedure is put forward. After choosing one plate i as the ideal wall, the other plate j of the two plates is projected on the chosen plate. Then the projected length and width are regarded as the new counterparts of the plate j. The new sizes are denoted by the variables with an apostrophe shown in Figure 8 and Figure 9 that present the projection for the two specific cases, as plate1 and plate2 are chosen as the ideal wall, respectively.
(a) Schematic diagram corresponding to case (a) in Figure 6 (b) Schematic diagram corresponding to case (b) in Figure 6      It can be seen from Figure 10 that the calculated forces are relatively small and the same no matter which plate is chosen as the ideal wall. Therefore, for the specific three-dimensional cases in which the projection of one plate on the other plate is a rectangle or two-dimensional case, the proposed projection procedure can eliminate the adverse effect of the choice of the ideal wall.  Figure 11 illustrates two specific cases of edge-edge type. The geometric parameters related to the two plates are the same as those shown in face-face type: l1 = w1 = 0.1, l2 = w2 = 0.2, d1 = d2 = 0.01. Two specific cases of relative positions of the two plates are shown in Figure 11 where various distance between plates is along Y direction in case (a) and X direction in case (b). According to the previous observation, the calculated force for these cases will be small since the corresponding   It can be seen from Figure 10 that the calculated forces are relatively small and the same no matter which plate is chosen as the ideal wall. Therefore, for the specific three-dimensional cases in which the projection of one plate on the other plate is a rectangle or two-dimensional case, the proposed projection procedure can eliminate the adverse effect of the choice of the ideal wall.  Figure 11 illustrates two specific cases of edge-edge type. The geometric parameters related to the two plates are the same as those shown in face-face type: l1 = w1 = 0.1, l2 = w2 = 0.2, d1 = d2 = 0.01. Two specific cases of relative positions of the two plates are shown in Figure 11 where various distance between plates is along Y direction in case (a) and X direction in case (b). According to the previous observation, the calculated force for these cases will be small since the corresponding projection area is small.  Figure 11 illustrates two specific cases of edge-edge type. The geometric parameters related to the two plates are the same as those shown in face-face type: l 1 = w 1 = 0.1, l 2 = w 2 = 0.2, d 1 = d 2 = 0.01. Two specific cases of relative positions of the two plates are shown in Figure 11 where various distance between plates is along Y direction in case (a) and X direction in case (b). According to the previous observation, the calculated force for these cases will be small since the corresponding projection area is small. Figure 10. Comparison of the calculated force of different choices based on projection procedure. Figure 11 illustrates two specific cases of edge-edge type. The geometric parameters related to the two plates are the same as those shown in face-face type: l1 = w1 = 0.1, l2 = w2 = 0.2, d1 = d2 = 0.01. Two specific cases of relative positions of the two plates are shown in Figure 11 where various distance between plates is along Y direction in case (a) and X direction in case (b). According to the previous observation, the calculated force for these cases will be small since the corresponding case (a): lx = 0.2, lz = 0.05 case (b): ly = −0.01, lz = 0.1 Figure 11. Two specific cases of edge-edge contact. D/nm Figure 11. Two specific cases of edge-edge contact. Figure 12 presents the ratios of the plate-plate solutions of the two specific cases of edge-edge type, F ee , to that of case (a) of face-face type shown in Figure 6, F ff , while keeping the same distance between two plates. It can be seen that the van der Waals forces in the case of the edge-edge type are much smaller than those of face-face type under the same conditions of plate size and distance between plates. So, it seems that van der Waals force for the case of edge-edge type can be ignored in DEM simulations. Therefore, the effect of the choice of the ideal wall will not be analyzed further.  Observations on the aspect of the effect of ideal wall choice are similar as those in the case of face-face type. Therefore, for the two specific cases shown in Figure 13, the abovementioned projection procedure has been performed before calculating the relevant plate-wall solutions as illustrated in Figures 14 and 15. Observations on the aspect of the effect of ideal wall choice are similar as those in the case of face-face type. Therefore, for the two specific cases shown in Figure 13, the abovementioned projection procedure has been performed before calculating the relevant plate-wall solutions as illustrated in Figures 14 and 15.      Figure 14. Schematic diagram of the projection step treating plate1 as the wall.

Edge-Edge Type
It can be seen form Figure 16 that unlike the results of the face-face type, the calculated forces are different for the different choice of the ideal wall even after the abovementioned projection procedure has been performed. For the case (a) shown in Figure 13, the calculated force corresponding to plate2 being the wall is larger than that to plate1 being the wall, and the difference between the calculated forces corresponding to the different choices for α = 30 • is smaller than that for α = 45 • . However, for the case (b) shown in Figure 13, the calculated force corresponding to plate2 being the wall is smaller than that to plate1 being the wall, while the difference between the calculated forces corresponding to the different choices for α = 30 • is smaller than that for α = 45 • . The above differences reflect the effect of the new length and distance after projection as shown in Figures 14 and 15. Comparing the above observations to the geometric parameters shown in Figures 14 and 15, it can be concluded that the change in the distance instead of length due to the different choice of ideal wall has the critical influence on the calculated forces.
(a) Schematic diagram corresponding to case (a) in Figure 13 (b) Schematic diagram corresponding to case (b) in Figure 13    It can be seen form Figure 16 that unlike the results of the face-face type, the calculated forces are different for the different choice of the ideal wall even after the abovementioned projection procedure has been performed. For the case (a) shown in Figure 13, the calculated force corresponding to plate2 being the wall is larger than that to plate1 being the wall, and the difference between the calculated forces corresponding to the different choices for α = 30° is smaller than that for α = 45°. However, for the case (b) shown in Figure 13, the calculated force corresponding to plate2 being the wall is smaller than that to plate1 being the wall, while the difference between the calculated forces corresponding to the different choices for α = 30° is smaller than that for α = 45°. The above differences reflect the effect of the new length and distance after projection as shown in Figures 14  and 15. Comparing the above observations to the geometric parameters shown in Figures 14 and 15, it can be concluded that the change in the distance instead of length due to the different choice of ideal wall has the critical influence on the calculated forces.  Figure 15. Schematic diagram of the projection step treating plate2 as the wall. Therefore, to reduce the overestimation of the ideal plate-wall model as much as possible, a reasonable choice of the ideal wall means a smaller distance between the two plates for the case of edge-face contact type.

Direction and Location of the Calculated Forces
The case (a) shown in Figure 13 has been chosen to investigate the difference in the direction and location of the calculated forces based on the ideal plate-wall and plate-plate wall model. As Therefore, to reduce the overestimation of the ideal plate-wall model as much as possible, a reasonable choice of the ideal wall means a smaller distance between the two plates for the case of edge-face contact type.

Direction and Location of the Calculated Forces
The case (a) shown in Figure 13 has been chosen to investigate the difference in the direction and location of the calculated forces based on the ideal plate-wall and plate-plate wall model. As illustrated in the last subsection, for the studied case, choosing plate1 as the wall led to a smaller calculated force. Therefore, the calculations based on the ideal plate-wall model treat plate1 as the ideal wall while the distance between the two plates is kept 0.001 µm, and the angle α increases from 15 • to 75 • .
According to the ideal plate-wall model [17], the direction of the total van der Waals force is perpendicular to the ideal wall. For the studied case, the direction is just along the S-axes shown in Figure 13a. To show the direction of the calculated force based on the plate-plate model, the coordinate components of the force based on plate-plate model are presented in Table 2. The last column is the direction of the calculated force based on the plate-wall model, with plate1 being the wall. When plate2 is treated as the wall, the direction of the calculated force is just along the Y-axes, shown in Figure 13a. It can be seen from Table 2 that the direction corresponding to plate-plate model is somewhere in between the direction corresponding to plate-wall model, with plate1 being a wall and plate2 being a wall. In addition, it can be seen that the difference in the directions between plate-plate and plate-wall model changes nonlinearly with the increasing angle between two plates. However, comparing with the traditional plate-wall model widely used in DEM [13][14][15][17][18][19], it would be reasonable if the direction of the total van der Waals force is taken as the average of the directions corresponding to plate-wall model with plate1 being a wall and plate2 being a wall. Table 3 shows the difference between the location of the calculated force based on the plate-plate model and that of the plate-wall model. Plate2 in Figure 13a is treated as the wall and the location is the acting point of the calculated force on the plate2. We use X pp and X pw to denote the X-coordinate of the acting point based on the plate-plate and plate-wall model, respectively. The other variables in Table 3 denote the relevant quantity in a similar way. It can be seen that the difference in the X-coordinate between the plate-plate and plate-wall model is most significant, the Z-coordinate is smallest, and the Y-coordinate is in between. This difference in the X-coordinate increases nonlinearly with the increasing angle between two particles. However, the projection length and the calculated van der Waals force becomes small with the increasing angle between two particles. Therefore, the influence of the above difference in the location of calculated force on the calculated moment will not be as large as the difference itself shown in Table 3.

Conclusions
Review on the literature related to the calculation of van der Waals force between two clay particles shows that the ideal plate-wall model has been widely used in clay DEM simulations. This study developed a plate-plate model without any assumption about the relative location of the two plates, and an optimized Cotes integration method has been used to calculate the van der Waals force numerically. Through the comparison of the calculated results based on the plate-wall model and the plate-plate model, the effect of the choice of the ideal wall has been investigated.
The calculated results show that for the face-face contact type, it is reasonable to choose the plate with a smaller area as the ideal wall in general. For the three-dimensional cases which can be simplified into two-dimensional cases, the projection procedure is advised before using the plate-wall model. This procedure chooses one plate i as the ideal wall, the other plate j of the two plates is projected on the chosen plate. Then the projected length and width are treated as the new counterparts of the plate j. The aim of the above is to reduce the overestimation of the ideal plate-wall model as much as possible.
The van der Waals forces in the case of edge-edge type are much smaller than those of face-face type. It seems reasonable that van der Waals force for the case of edge-edge type can be ignored in DEM simulations.
In additional, for the edge-face contact type, a reasonable choice of the ideal wall is to choose one plate as the ideal wall so that there is a smaller distance between the two plates. It is reasonable that the direction of the total van der Waals force is taken as the average of the directions corresponding to the plate-wall model, with plate1 being a wall and plate2 being a wall. The difference in the location of the calculated force based on the plate-wall model and plate-plate model increases nonlinearly with the increasing angle between clay particles. It has limited effect on the calculated moment.

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