Marine clay is viewed as an assembly of clay platelets in this study. Each clay platelet is assumed to be extruded based on a planar shaped polygon with multiple edges. The extrusion distance is the thickness of the clay platelet. This assumption implies that the numerical clay in this study is applicable to clay minerals, such as kaolinite and illite, that do not present a significant spatially curving shape. The interactions between two clay platelets include double layer repulsive interaction, van der Waals attractive interaction, and mechanical repulsive interaction. The former two interactions are closely related to the shapes and relative alignment of two neighboring clay platelets of interest. Although the pore water is not explicitly simulated in the DEM model, the pore water effects are implicitly considered in the interaction law between clay platelets as outlined below.
2.1. Calculation Scheme of Double Layer Repulsive Interaction
Figure 1 presents an arbitrary alignment of two clay platelets. Without loss of generality, the planar shapes of clay platelet 1 and clay platelet 2 are assumed to have eight and six edges, respectively. As shown in
Figure 1a, the normal directions of the two platelets are denoted by
n1 and
n2. A mid-plane is defined as the plane that bisects the acute angle made by the two platelets. A local coordinate system is defined with the normal direction of the mid-plane as
z axis, and the direction of
n2 ×
n1 as
y axis, while the origin can be put anywhere in the mid-plane. Since the clay platelet surfaces are negatively charged, then in
Figure 1a, the lower facet of clay platelet 1 and the upper facet of clay platelet 2 constitute a pair of inclinedly aligned, negatively charged facets. The electric interaction of the two clay platelets can be attributed to the interaction of the two facets.
In
Figure 1b, projecting the two facets on the mid-plane gives two intersecting polygons, as outlined by the thin dashed line segments. The intersection area is outlined by solid line segments. A side-view along the positive-
y direction is presented in
Figure 1c for easy understanding. The intersection area is the range of effective double layer repulsive interaction, while other parts of the two facets’ projections contribute little to the interaction and are omitted.
Following the method in [
11], the repulsive force can be calculated in a simplified way as follows. First, as presented in
Figure 1b,c, the intersection area is divided into bands along the
x-axis direction with a band width of
w. The two facets (lower facet of clay platelet 1 and upper facet of clay platelet 2) are also divided in this way, as depicted in
Figure 1c. The interaction between the two inclinedly aligned facets is now simplified as the summation of the interactions between a series of oppositely positioned band pairs with respect to the mid-plane, one from clay platelet 1 and the other from clay platelet 2 in a pair. Take band
i in
Figure 1b as example. The band has a length of
li in the
y-axis direction. The center distance of the two corresponding bands from the upper and lower platelet is
hi in the
z-axis direction, as shown in
Figure 1c. According to the work in [
11], the double layer repulsive force between the two inclinedly aligned bands,
Fe,i, can be approximated by that between two parallelly aligned bands with the same distance
hi, that is
where
n is ion concentration in the pore water,
k = 1.38 × 10
−23 J/K is the Boltzmann constant,
T is Kelvin temperature, and
is the normalized mid-plane electric potential of two parallel plates at a distance of
hi. In Equation (1),
Fe,i’s direction is normal to the mid-plane.
The normalized mid-plane electric potential can be obtained by numerically solving the Poisson–Boltzmann equation and the solutions are presented in
Figure 2. The normalized mid-plane electric potential is defined as
where
v is the ion valence,
e* = 1.6 × 10
−19 C is the charge of an electron, and
is the mid-plane electric potential.
In
Figure 2, the normalized distance is defined as
where
is the relative dielectric constant of pore water and
= 8.8542 × 10
−12 C/(V·m) is the dielectric constant in vacuum condition.
Figure 2 shows that the normalized mid-plane electric potential decreases non-linearly with the normalized distance. The curves also depend on the normalized surface electric potential of the clay platelet defined as
, with
being the surface electric potential.
By the summation operation, the double layer repulsive force,
Fe,1, and moment,
Me,1, applied on clay platelet 1 by clay platelet 2 can be expressed in the local coordinate system as
where
Nx is the number of divided bands in
Figure 1b, and
xi and
yi are the coordinates of band centers in the local coordinate system.
The above derivations apply to any convex planar shape of clay platelet. It is noted that the core idea of the above simplified calculation scheme of a double layer repulsive interaction is to replace the inclinedly aligned facets with a series of staircase-like, parallelly aligned bands. The finer division of the bands leads to the higher accuracy. To address this issue, three clay platelets are used as illustrative examples. As shown in
Figure 3, without loss of generality, the planar shapes of clay platelets A, B, and C are assumed to have four, seven, and five edges, respectively. The planar extent of the platelets is between 0.5 and 0.6 μm, and the thickness is 0.06 μm (typical size range of kaolin platelets). The alignments of the three clay platelets are presented as insets in the following related figures. If not stated otherwise, the parameters are
T = 293 K,
= 50 mV,
n = 0.0015 mol/L,
= 80.1, and
v = 1.0.
Figure 4 presents the variations of double layer repulsive force with the inclination angle
θ between clay platelets B and C under different band division number
Nx. When the inclination angle
θ is changed, the net distance of the two clay platelets,
dn, is kept constant at 1 nm.
Figure 4 shows that the double layer repulsive force decays very quickly as the inclination angle increases. When
θ is greater than 7°, the repulsive force is lower than 1‰ of that when
θ = 0° and it can be omitted. This conclusion is intrinsically consistent with that in [
16] for clay platelets with rectangular planar shape.
Figure 4 also indicates that when
Nx exceeds 8, the repulsive force will only change quite slightly as
Nx further increases. It can be concluded that
Nx = 8 is fine enough to divide the bind for the accurate calculation of double layer repulsive force.
Figure 5 presents the variations of double layer repulsive force with the relative inclination angle and the net distance between clay platelets A and B under various surface electric potentials. The division number
Nx takes 8 here.
Figure 5a shows that the repulsive force can be omitted when the inclination angle exceeds 10°, similar to conclusions from
Figure 4.
Figure 5b shows that the repulsive force decays very quickly as the net distance
dn increases. When
dn is greater than 10 nm, the repulsive force is lower than 1‰ of that when
dn = 1 nm and can be omitted.
2.2. Calculation Scheme of Van Der Waals Attractive Interaction
The van der Waals interaction between clay platelet 1 and clay platelet 2 is the summation of all interactions of molecule pairs respectively belonging to the two platelets. The van der Waals attractive energy,
u, between two molecules are first given by in [
17] as
where
B is the London constant,
is a vector connecting the centroids of two molecules, and thus
is the centroid distance of the two molecules.
The van der Waals attractive force between two molecules can be obtained by differentiation of Equation (6) with respect to
as
In
Figure 6a, the outlined solid volume of clay platelet 1 corresponds to the intersection area in the mid-plane (see
Figure 1) and is denoted as clay platelet 1′. The same applies to clay platelet 2 and clay platelet 2′. Since the van der Waals attractive force decays with the distance very fast (to a power of −7), the intersecting clay platelets 1′ and 2′ would make the most contribution to the interaction while the other non-intersecting parts’ contribution can be omitted. Similar to the band division in
Figure 1b, the clay platelet 1′ is divided into bands and the van der Waals interaction between clay platelet 2′ and every divided band is calculated. In
Figure 6b, taking band
I of clay platelet 1′ as example, its van der Waals interaction with clay platelet 2′ can be approximated by the interaction with an in-plane infinitely extended platelet with the same thickness as platelet 2′, which has an analytical solution. The extended volume has little influence on the interaction since it is located far away from band
i.
Integrating Equation (7) over the volume of an in-plane infinitely large platelet with a thickness of
t gives its van de Waals attractive force with a point as
where
h is the distance of the point to the upper surface of the platelet, and
ρ1 is the molecule density in the platelet.
Integrating Equation (8) over the volume of band
I gives the van der Waals attractive force applied on band
I as
where
A = is the Hamaker constant,
is the acute angle made by the two clay platelets, and
d1–
d4 are the distances of the four vertices of band
I to the upper surface of clay platelet 2′, as illustrated in
Figure 6b.
Then, the van der Waals attractive force,
Fv,1, and moment,
Mv,1, applied on clay platelet 1 by platelet 2 can be determined by
Similar to the results in
Figure 4, although not given here for conciseness, it is found that
Nx = 8 is fine enough to divide the clay platelet in
Figure 6 for the accurate calculation of van der Waals attractive force.
When two clay platelets are aligned with the acute angle made by them being large, say greater than 70°, the calculation scheme presented above can be replaced by a simpler scheme as discussed below. As presented in
Figure 7, the
z-axis of the local coordinate system is aligned with the normal direction of clay platelet 2. Denote the closest point on clay platelet 1 to clay platelet 2 by
E and
EA is the extrusion edge of clay platelet 1. Among the two lateral facets sharing edge EA, i.e., facet
EABF and facet
EADH, facet
EABF makes a smaller angle with plane
oxy than that facet
EADH makes. A thin cuboid,
EABF-IJKL, is constructed by taking the facet
EABF as one of its lateral facets and extending in the direction away from clay platelet 2. The extending distance takes the average edge length of the planar polygon of clay platelet 1. In fact, this average edge length is an ideal choice since the thus constructed cuboid would cover a large portion of the volume of clay platelet 1. A slight deviation from this extending distance has a very limited influence since the van der Waals interaction is mostly contributed by clay platelet 1′s molecules located near the clay platelet 2. An infinitely large clay platelet 2′ is constructed the same way as in
Figure 6b.
With the above considerations, the van der Waals interaction between clay platelet 1 and platelet 2 can be approximated by the interaction between the extended thin cuboid and the infinitely large clay platelet 2′, which has an analytical solution.
Integrating Equation (8) over the volume of the thin cuboid gives the van der Waals attractive force applied on clay platelet 1 as
where
hA–hL are the distances from the thin cuboid’s vertices to the upper facet of clay platelet 2, and
,
, and
are the
z components of the unit vectors along directions
,
, and
. Similarly, the moment can be integrated.
Figure 8 presents the variations of van der Waals attractive force with the relative inclination angle and the net distance between clay platelets A and B under various Hamaker constants. The division number
Nx takes 8 here. In
Figure 8a, when the inclination angle is less than 73°, Equation (9) is used, and otherwise Equation (12) is used. Here, 73° is just for illustration and other delimiting value around 70° can be used as well without changing the conclusion below. In
Figure 8a, the smooth transition across 73° indicates that Equations (9) and (12) are both properly developed.
Figure 8a shows that the van der Waals attractive force first quickly decreases as the inclination angle increases and then presents a slow decaying rate until the transitional angle about 45° is reached. After that, an increase of the attractive force is observed with the inclination angle and the increasing rate soars after the inclination angle exceeds 80°. Again, this conclusion is intrinsically consistent with that in [
16] for clay platelets with rectangular planar shape.
Figure 8a also shows that the van der Waals attractive force can be omitted when the inclination angle is between 10° and 80°, with the attractive force less than 1% of that in parallelly aligned platelets case.
Figure 8b shows that the attractive force decays very quickly as the net distance
dn increases. When
dn is greater than 10 nm, the attractive force is lower than 1‰ of that when
dn = 1 nm and can be omitted.
2.3. Total Interaction between Two Clay Platelets
In DEM simulations, a contact model is applied to describe the interaction between two discrete objects at the contact. The objects in three-dimensional DEM simulations can be balls, rigid blocks, ball clumps, walls, etc. In traditional DEM, two objects can interact only when they are in physical contact, i.e., dn = 0 nm. However, for the DEM simulation of clay, the total force/moment received by a clay platelet from surrounding platelets is a summation of traditional mechanical contact interaction (present when two objects being in physical touch) and the double layer repulsion and van der Waals attractive interactions (considered in certain net distance and relative inclination angle ranges as discussed above).
The double layer repulsion and van der Waals attractive interaction calculation scheme were implemented into DEM software PFC3D [
18]. The built-in linear rolling resistance contact model in PFC3D [
18] is used as the physical touch interaction, where the normal, tangential, and rolling interactions can be expressed as
where
Fn and
Fs are normal and tangential forces, respectively,
Mr is the rolling resistance moment,
Kn,
Ks,
Kr are the corresponding contact stiffnesses,
un is the overlap,
nc is the unit contact normal vector,
is the incremental tangential displacement,
is the incremental rolling angle, μ is the frictional coefficient, μ
r is the rolling resistance coefficient, and
R is an equivalent contact radius.
In Equations (13)–(15), the contact stiffnesses are
where
E is the contact modulus and
ξk is the normal-tangential stiffness ratio.
Note that in Equations (9) and (12), when the net distance between two clay platelets is small, say less than 1 nm, the van der Waals attractive force would be very large. In fact, in net distance range smaller than 1 nm, the repulsive van der Waals interaction (not discussed above due to irrelevance) will start to dominate and the absorbed water molecule layers on clay platelet surfaces will start to overlap to generate repulsive force. The two issues are abstracted by setting a mechanical cut-off gap
dm and when the net distance
dn is less than
dm, the following two aspects are considered: (1) van der Waals attractive and double layer repulsive interactions are kept constant, i.e., no change with decreasing
dn; (2) the starting point of traditional mechanical contact interaction is shifted from
dn = 0 to
dn =
dm [
11].
Figure 9 visualizes the interactions between clay platelets A and C (see
Figure 3), where the mechanical cut-off distance
dm = 1.0 nm and the normal contact stiffness
Kn = 60 N/m. In
Figure 9, the net force is a summation of the double layer repulsive force, van der Waals attractive force and mechanical contact repulsive force, with positive value representing repulsion. It is interesting to note that the net force is repulsive when
dn is less than 1.25 nm and a local maximum repulsion is observed at
dn = 2.0 nm. The net force becomes attractive when
dn is less than 1.25 nm and reaches a local maximum attraction at
dn =
dm (=1.0 nm here). When
dn is less than
dm, the mechanical contact repulsion starts to dominate.
2.5. DEM Simulation of Sedimentation and Consolidation Processses
The commercial software PFC3D [
18] was used for the numerical simulation. In PFC3D, clay platelets are represented by rigid blocks generated from the above constructed nine templates. Hence, 500 clay platelets are generated with random positions and directions in a slender box composed of six boundary walls, as shown in
Figure 11a. The corresponding void ratio
e is 21.0, avoiding any overlap between platelets. These clay platelets are randomly drawn from the nine templates. The solid density is set to 1650 kg/m
3 (considering the buoyancy). Note that 500 platelets are enough to capture the behavior focused on in this study, and the related computation cost is acceptable, which will be discussed later.
The model parameters are listed in
Table 1. A Hamaker constant of 1.0 × 10
−19 J is typical for a kaolinite-NaCl-water system. The surface electric potential may change in a wide range for clay minerals and a typical value of 200 mV is applied here. Numerous experimental data indicate that pore water chemistry can significantly change the physical and mechanical behavior of clay. In our simulations, the ion concentration is changed to investigate such effects. Salt concentration in marine clay pores poses large variation, depending on the specific location [
19]. For example, near river mouth where a river meets the ocean, the sea water can be diluted to a large extent. In
Table 1, the ion concentration of 0.64 mol/L is a typical up-bound value for most ocean area while the reduced concentrations can be representative of diluted cases. The physical contact model parameters are selected based on common DEM simulation practice.
The clay platelets in the box are allowed to settle down under gravity to simulate the sedimentation process under static water condition. The settled clay platelet assembly is then subjected to one-dimensional compression to simulate the consolidation process due to the weight of subsequent sediments lying above. One-dimensional compression is numerically achieved by moving the top boundary wall downward. The moving rate is very low to allow full arrangement adjustment of clay platelets.