Modeling of Adhesive Particles Using a Combination of the Two-Body Interaction and Phase-Field Methods

Takuya Uehara

Department of Mechanical Systems Engineering, Yamagata University, Yonezawa, Japan.

DOI: 10.4236/ojmsi.2020.82003                

Abstract

Discrete materials such as powders and granular materials have been widely used due to their specific characteristics. The precise evaluation is accordingly becoming important, and various numerical schemes have been developed. However, the interactions among the constituent particles are still difficult to model precisely. Especially, contact conditions, which vary with material properties and circumstances, are difficult to formulate. In this study, a computational model for simulating adhesive particles on contact in a many-particle system is proposed. The interaction between the particles was represented by a two-body repulsive force that depends on the distance between particles and an additional adhesive force at the contact point. A phase-field variable was introduced to express the surface of each particle, and the adhesive force was formulated using the phase-field distribution. As a result, the adhesion of particles was properly expressed. For a mono-particle system, neighboring particles adhered and uniformly aggregated, while for a dual-particle system, several characteristic patterns were obtained depending on the initial arrangement of the particles. Repulsive contact was also considered as a specific case, and the corresponding results were obtained.

Keywords

Phase Field Model, Two-Body Interaction, Multibody System, Adhesive Force, Particle Aggregation, Computer Simulation

1. Introduction

Solid materials can be categorized into two types: continuum materials, which include metals, ceramics, and polymers, and discrete materials, which include sands, powders, and other granular substances. In engineering, continuum materials are used more often than discrete materials because discrete materials cannot sustain mechanical force. However, the application of discrete materials has recently expanded due to some of their specific characteristics both in macroscopic properties and microscopic structures [1] [2]. For example, the vacant space among their constituent elements or particles enables reduction in weight and noise- and vibration-absorbing properties, and the flexibility in shape-forming and in blending component is applied to materials processing [3] [4]. Accordingly, the precise evaluation and prediction of the behavior of discrete materials are becoming increasingly important.

The mechanical behavior of continuum materials can be described by conventional continuum mechanics, and the properties specific to the objective material, such as elasticity, plasticity, viscosity, and their combinations, can be characterized by their constitutive laws. On the other hand, the mechanics and dynamics of discrete materials are difficult to describe with continuum mechanics, even though their constituent elements are continuum materials. Their fluid-like behavior also makes it difficult to describe in the continuum mechanical framework [5]. One possible way to describe these materials is to treat each element as a particle and to solve the equations of motion for every particle. There are several variations in such computational methods [6], one of which is the discrete element method or distinct element method (DEM in both cases) [7] [8] [9]. Application of this method to practical engineering problems was once difficult because of the overwhelming number of particles but has become possible owing to rapid progress in computer technology. However, as this method is widely used, its accuracy has been raised as a critical issue. One of the most important points in this respect is how to treat the contact condition [10] [11]. The distance between two particles is usually applied to judge their contact; when the distance is below a threshold value, the two particles are considered to be in contact, and a repulsive force is imposed. In particular, a rigid- or soft-core potential is applied to calculate a repulsive force. Similarly, an adhesive force can also be imposed using the two-body term as a function of the distance between the particles, and we have demonstrated DEM simulation for investigating the stability of particle packing structure [12]. However, particular contact conditions, such as those in a case depending on temperature and humidity, are difficult to formulate.

In this study, we were thus motivated to express the contact state using a phase-field model, which is able to solve various interfacial problems [13] - [17]. One of the successful fields is the phase transformation, such as that between solid and liquid, and the interface may be mobile, depending on thermodynamic variables. We have also demonstrated phase-field simulation on various problems such as metal-foam structure [18] and domain tessellation [19]. Application of the phase-field model to the contact problem is now a prospective approach and may be possible by way of suitable numerical methods. The contact of two particles is similar to the solid-liquid interface or grain boundary in a polycrystalline material represented in the conventional phase-field models, for which a coupling with continuum stress/strain field was proposed before [20]. However, the contact face is different in the point that the target is not a continuum material but separated particles, and the mechanical force should also be calculated independently. As a first step to model this problem, we propose a simple expression of the contact state in a many-particle system using a phase-field variable. In particular, an adhesive force in the interfacial region is assumed as a function of the phase-field value. Simulations are then demonstrated, and the validity of the function is presented.

2. Method

2.1. Two-Body Interaction between Particles

The main force acting on each particle can be represented by a two-body interaction, with Newton’s equation of motion being numerically solved for every particle. Generally, a strong repulsive force at a small distance and a weak attractive force at a large distance are represented by Lennard-Jones potential. In this study, only a repulsive force was assumed in the two-body interaction to make evident the effect of the adhesive term using a phase-field variable. The power index was chosen as n = 16 so that the force would converge smoothly to zero in an adequate range as follows:

$f_{rp}=\frac{E}{d}{\left(\frac{D_{rp}}{d}\right)}^{16}$, (1)

where, d is the distance between the center of particles, D_rp is the standard distance, and E is a parameter representing the strength of the force.

2.2. Phase-Field Model at Contact Region

To express the adhesive force, a phase-field model was applied to the surface area of the particles. The phase-field variable f was set as 1 in the particle and 0 in the outside region, with the value varying smoothly but steeply from 1 to 0 over the surface area of a certain width w, as schematically illustrated in Figure 1. The variation was calculated using the following equation:

$\tau \stackrel{\dot{}}{\varphi }={\epsilon }^2{\nabla }^2\varphi +\varphi \left(1-\varphi \right)\left(\varphi -0.5+M\right)$, (2)

where $\tau$ and $\epsilon$ are parameters and M is a term that controls transition of the phase field value $\varphi$. For example, if M is set as positive, the boundary line, which is denoted by a contour line of $\varphi$ = 0.5, moves toward domain represented as $\varphi$ = 0 and vice versa for a negative M value. The equilibrium distribution is approximated as

$\varphi =0.5+0.5\tanh \left(x/\delta \right)$, (3)

where δ is a parameter representing the width of the transient region, termed interfacial thickness, and x denotes the middle point in the interfacial region. Figure 1(b) represents Equation (3) with δ = 2, demonstrating that the value is distributed within a width of approximately 5δ.

Figure 1. Schematic illustration of the phase field on a particle surface and the initial distribution of $\varphi$ approximated by a hyperbolic tangent function.

In this study, the center of the interface was set to r = R_pf, where r is the distance from the center of particle. The parameter M was set as +M_pf for r < R_pf − w/2 and converges smoothly to 0 within a range of r < R_pf + w/2 so that the phase-field distribution would maintain its position on the surface.

2.3. Determination of Contact

When the two particles come close to one another, the phase-field variable centered in each particle overlaps in the middle range, and the phase fields then tend to merge. This phenomenon represents the adhesion of the particles; from a mechanical viewpoint, an adhesive force is created. In the present model, this force was represented in the two-body interaction as an attractive force, according to the following formula:

$f_{at}=k_{at}{\varphi }_{mid}{\left(r_{mid}-w\right)}^2{\left(r_{mid}+w\right)}^2/w^4$, (4)

where f_mid is the phase-field value at the middle point between two particles, r_mid is the deviation in distance from the equilibrium two-particle distance D_pf(= 2R_pf, i.e., r_mid = d − D_pf, where d is the distance between two particles), w is the width of the range in which the adhesive force functions, and k_at is a parameter. The function is a double-well type, has a peak value at r_mid = 0, and is distributed in a range of ±w, as illustrated in Figure 2(a). When the interparticle distance falls within this range, an attractive force is created according to Equation (4), and the absolute value is controlled by the phase-field value at the middle point. If the phase field is merged and has a value close to 1, the attractive force works, while the force is neglected if the phase field remains at a low value. The tendency of phase-field coalescence would depend on the phase-field parameter representing the interfacial properties, which generally depend on various material and environmental parameters. Detail on this point will be discussed in a future work.

Consequently, the total force acting between two particles can be summarized as follows:

Figure 2. A function for representing adhesive force in the contact region (a) and corresponding two-body force (b).

$f=\{\begin{array}{l}{f}_{rp}\text{for}d<D_{pf}-w\hfill \\ f_{rp}+f_{at}\text{for}{D}_{pf}-w<d<D_{pf}-w\hfill \\ 0\text{for}d>D_{pf}+w\hfill \end{array}$. (5)

Note that f_at does not always work, even though the distance falls in the contact range if the phase-field value is zero. The total force for f_md = 1 is illustrated in Figure 2(b).

2.4. Dual-Particle System

The phase-field model elucidated in Section 2.2 is used to distinguish whether the point is inside or outside a particle, the phase field being merged when two particles come into contact. In this model, however, a combination of particle types could not be dealt with. Therefore, the phase-field model was extended to consider a system containing two types of particles by assigning two different values to the phase field as $\varphi$ = +1 and −1. The fundamental behavior of the model was the same as that explained in Sections 2.2 and 2.3, while the center of interfacial area was represented by an contour line of $\varphi$ = 0. The attractive force was provided in Equation (4), by modifying the sign of k_at. This means that the attractive force does not work between the different types of particles, as the phase-field value at the middle point f_md = 0. As a result, it was expected that two particles of the same type would adhere, while particles of the different type would repel due to the repulsive force.

2.5. Computational Procedure

Simulation code was programmed according to the following procedure: 1) initial setting, 2) calculation of the two-body force based on the interparticle distance and the phase-field distribution, 3) movement of the particles by solving the equation of motion, 4) update of the phase-field distribution by numerically solving Equation (2) using the finite difference method (FDM), and 5) repeating steps (2) - (4) to obtain sufficient time steps. This flow is summarized in Figure 3.

Particles were initially arranged on the points of a regular square lattice, and the initial phase-field distribution was provided according to Equation (3). The

Figure 3. Flow chart of the computational procedure.

velocity-Verlet algorithm was used for numerical integration of the equation of motion, and equally spaced grids were used for FDM.

3. Model and Conditions

A two-dimensional model was applied in this study. A square region with dimension L × L and with periodic boundary condition was considered. Particles were arranged at equally spaced N_p × N_p lattice points, and their initial deviation was randomly determined. The initial velocity of each particle was also determined randomly, while the total kinetic energy was controlled at a certain level, K_E0, to prevent extreme motion. The parameters are listed in Table 1, where ΔL is the gird interval for FDM calculation and Δt is the time increment for integration. The length of the calculation region L was set as (D₀ + δ_L) N_p= N_LΔL, where N_L was the number of grids in the x and y directions. When δ_L = 0, the particles aligned along the x and y directions to contact each other at the equilibrium distance D₀. The parameters marked with asterisks, w and δ_L, were varied in every case, and the standard value is given in the table. Note that all variation in length is standardized by particle radius (i.e., R₀ = D₀/2 = 1.0) and that energy and time are also non-dimensional.

The first 1000 time steps, which were sufficient to obtain an equilibrium state in the present model, were devoted to relaxation of the phase-field distribution; the particles were fixed at their initial position and only the phase-field equation was solved. Then, the motion of the particles as well as the phase-field distribution was simulated.

4. Simulation

4.1. Mono-Particle System

A mono-particle system, in which all particles had the same radius and surface properties, was first simulated. Figures 4(a)-(c) show the results for δ_L = 0.0, 0.1, and 0.2, respectively. When the model volume was set to δ_L = 0.1, particles were moderately packed, and neighboring particles were in contact with one another after the initial relaxation without motion of particles, as shown in Figure 4(b) (ii). The particles were then made mobile at the 1000th time step by

Table 1. Simulation parameters.

Figure 4. Simulation results for mono-particle system. The initial 1000 steps were devoted to relaxation without particle motion, and the particles were made mobile at the 1000th time step. The color indicates the phase field value in (i) - (iii), and the contour line of $\varphi$ = 0.5 is drawn in (iv).

solving the equation of motion; however, as shown in Figure 4(b) (iii), the particles did not move at all. This was because the contact distance at the initial configuration was equal to the equilibrium distance of the balance between repulsive and adhesive forces.

When the volume of the system was made smaller, the particles were packed more densely. This meant that the repulsive force was present at the initial configuration, and particles changed their position by moving alternately, as shown in Figure 4(a) (iii). It should be noted that when we performed similar simulations using a conventional Lennard-Jones potential, a hexagonal arrangement was obtained. In this case, local adhesion was formed among the four initially neighboring particles, and the contact was maintained. As a result, an arrangement different from the hexagonal one was obtained. In contrast, when the density was low, some of the contacts between particles were lost and a clear cleavage formed, as shown in Figure 4(c) (iii). The particles then moved freely before finally settling into equilibrium in a hexagonal arrangement.

4.2. Dual-Particle System

A dual-particle system was next considered. The phase-field parameters were set differently depending on the particle types so that adhesive force would work only between the same types of particles and no additional force would work between two different types. The radius was identical for both types. Figures 5(a)-(c) show the results for δ_L = 0.2, 0.1, and 0.1, respectively. The color indicates the phase-field value; red and blue represent $\varphi$ = 1.0 and −1.0, respectively, corresponding to the two types of particles, while green indicates $\varphi$ = 0 outside the particles. When the system volume was large (δ_L = 0.1, 0.2), particles of the same type came close and adhesive contacts were generated. Consequently, particle pairs were formed, as shown in Figure 5(a) and Figure 5(b). Initially, each particle had four primary neighbors of different particle type and four secondary neighbors of the same type. The primary neighbors were moved away by repulsive forces and adhesion occurred with one of the secondary neighbors. In contrast, when the volume was small and particles were closely packed, particles of the same type tended to form a linear connection, and a stripe pattern was formed, as shown in Figure 5(c) (iv). From the initial state, pairs of particles first formed at the 2000th time step, as shown in Figure 5(c) (iii). Then, multiple pairs attracted one another and arranged in a line, except for the partial unconnected pairs, which remained in their original places. Note that the vertical line formation was by chance. A horizontal stripe pattern was obtained in another trial, as will be shown later.

4.3. Trials with Different Initial Conditions

The effects of the initial conditions and computational parameters on the results of the dual-particle system were next investigated. As explained earlier, random numbers were used to determine the initial position and velocity of particles. The random series were changed, and simulations were carried out several times for each set of initial conditions. Figure 6 shows some of the typical results for the conditions presented in Figure 5, where the variable Rndrefers the maximum deviation from the complete lattice points in the initial configuration, and “Trial 2” denotes a trial with the same conditions but with different random series. When the initial configuration deviated from the lattice points, the symmetry was broken. When the volume was large and there was significant space for the particles to move around, the influence of the initial condition was negligible, and a similar pattern was observed, as shown in Figure 6(a). In contrast, when the volume was small and particles were closely packed, it became difficult for the particles to form a symmetrical pattern. Instead, contact with more than

Figure 5. Simulation results for a dual-particle system. The color indicates the difference in types of particles. Red represents $\varphi$ = 1 and blue represents $\varphi$ = −1.

Figure 6. Simulation results for a dual-particle system with varied initial conditions, representing the effect of random series used in the initial condition.

two particles occurred more frequently, and a linear connection tended to form. The linear shape was intermittent when δ_L = 0.1, as shown in Figure 6(b), whereas a complete stripe pattern through the periodic boundary condition was obtained when δ_L = 0.0, as shown in Figure 6(c).

Figure 7 presents the results for different initial configurations. When the same types of particles were initially aligned linearly, this arrangement was maintained throughout the simulation, while the connection was broken when δ_L = 0.2. In all cases, the particle position shifted slightly by keeping the linear connection so that the particle position became alternate. Similarly, when particles of the same type were arranged in a block, as shown in Figure 7(b) (i), the block pattern was maintained, while the local positions of the particles were rearranged depending on the system volume.

4.4. Effect of the Phase-Field Parameter

In the current model, the adhesive properties are represented by phase-field parameters. The real situation should be investigated carefully, but a proposal of the modeling is the objective of this paper; therefore, only a sample is presented here.

The interface thickness w was varied, and the results when w = 0.1, one half of the value used in the case presented in Figure 5, are shown in Figure 8. The width of the surface area, represented by a gradation in color, obviously became thinner. In addition, the adhesive tendency was affected; no adhesion pair was observed for δ_L = 0.2, and no linear connection or stripe pattern was observed for δ_L = 0.0, as shown in Figure 8(a) and Figure 8(c), respectively. No specific effect was observed for moderate density (δ_L = 0.1), as shown in Figure 8(b). This feature reveals the effectiveness of the present model not only at the rigid contact but also in a wide range of adhesive zones.

4.5. Repulsive Particles

Finally, a repulsive force was assumed in the interfacial area as a special case. This condition was easily realized by changing the parameter K_pf to a negative

value. Figure 9 shows the results for K_pf = −0.1 with w = 0.1 and δ_L = 0.2. For the regularly alternate initial configuration, an adhesive connection did not occur, and the alternate order was maintained, as shown in Figure 9(a). In addition, the local position was rearranged into a hexagonal configuration with all particles being dispersed equally. In the case that a stripe pattern was initially assigned, the repulsive interface detached the connected particles, and consequently, the same alternate pattern was formed, as shown in Figure 9(b).

These results, presented in Sections 4.4 and 4.5, reveal that the surface properties and connective characteristics of particles can be varied by changing the phase-field parameters; thus, it can be concluded that the effectiveness of the present model has been adequately verified.

5. Conclusions

In this study, a computational model for simulating particle adhesion in a many-particle system was proposed. Interaction was represented by a combination of two-body interaction and phase-field models. A prominent repulsive

Figure 7. Simulation results for a dual-particle system with varied initial conditions, representing the effect of the initial configuration.

Figure 8. Effect of the phase-field parameter with an interface thickness of w = 0.1 (in contrast with w = 0.2 in Figure 5).

Figure 9. Simulation results in the case of repulsive contact.

force was represented by a simple two-body force that depended on the interparticle distance. A phase-field variable was introduced to express the surface of each particle, and the adhesive force was formulated using the phase-field value at the contact region. As a result, the adhesion of particles was adequately expressed, and several aggregate patterns were obtained. The effect of the interfacial thickness was investigated, and it was confirmed that the adhesion condition can be controlled by varying the phase-field parameters. Repulsive contact was also assigned as a specific case, and the corresponding results were obtained. These results demonstrate the effectiveness of the present model for solving further practical problems.

In this study, the phase field was used to distinguish the contact state and was applied as a scaling factor to a two-body adhesive force. The function form, however, is not fixed yet, and the construction of a better formula should be explored. The effect on the tangential direction, or a frictional force, should also be considered to make the model more realistic. Improving these points will lead to a variety of further applications of this model, such as to mixtures of particles, aggregation in fluid, and the stability of particle packing structure.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Hinrichsen, H. and Wolf, D.E. (2004) The Physics of Granular Media. Wiley-VCH, Hoboken. https://doi.org/10.1002/352760362X
[2]	Cambou, B., Mgoariec, H. and Nguyen, N.-S. (2016) Granular Materials at Meso-Scale: Towards a Change of Scale Approach. ISTE Press and Elsevier, Amsterdam.
[3]	Fennis, S.A.A.M. and Walraven, J.C. (2012) Using Particle Packing Technology for Sustainable Concrete Mixture Design. Heron, 57, 73-102.
[4]	Tsukerman, S.A. (2013) Powder Metallurgy. Pergamon.
[5]	Andreotti, B., Forterre, Y. and Pouliquen, O. (2013) Granular Media: Between Fluid and Solid. Cambridge University Press, Cambridge.
[6]	Poeschel, T. and Schwager, T. (2010) Computational Granular Dynamics: Models and Algorithms. Springer-Verlag, Berlin.
[7]	Cundall, P.A. and Strack, O.D.L. (1979) A Discrete Numerical Model for Granular Assemblies. Geotechnique, 29, 47-65. https://doi.org/10.1680/geot.1979.29.1.47
[8]	O’Sullivan, C. (2011) Particulate Discrete Element Modelling: A Geomechanics Perspective. CRC Press, Boca Raton.
[9]	Sakai, M. (2012) Numerical Simulation of Granular Flows. Maruzen, Tokyo. (In Japanese)
[10]	O’Sullivan, C. and Cui, L. (2009) Micromechanics of Granular Material Response during Load Reversals: Combined DEM and Experimental Study. Powder Technology, 193, 289-302. https://doi.org/10.1016/j.powtec.2009.03.003
[11]	Mio, H., Shimosaka, A., Shirakawa, Y. and Hidaka, J. (2007) Cell Optimization for Fast Contact Detection in the Discrete Element Method Algorithm. Advanced Powder Technology, 18, 441-453. https://doi.org/10.1163/156855207781389519
[12]	Uehara, T. (2018) Modeling and Simulation of Particle-Packing Structures and Their Stability Using the Distinct Element Method. Open Journal of Modelling and Simulation, 6, 59-70. https://doi.org/10.4236/ojmsi.2018.64005
[13]	Emmerich, H. (2003) The Diffuse Interface Approach in Materials Science—Thermodynamic Concepts and Applications of Phase-Field Models. Springer-Verlag, Berlin, Heidelberg.
[14]	Provatas, N. and Elder, K. (2010) Phase-Field Methods in Materials Science and Engineering. Wiley-VCH, Hoboken. https://doi.org/10.1002/9783527631520
[15]	Biner, S.B. (2017) Programming Phase-Field Modeling. Springer International Publishing, Berlin. https://doi.org/10.1007/978-3-319-41196-5
[16]	Steinbach, I. and Shchyglo, O. (2011) Phase-Field Modelling of Microstructure Evolution in Solids: Perspectives and Challenges. Current Opinion in Solid State and Materials Science, 15, 87-92. https://doi.org/10.1016/j.cossms.2011.01.001
[17]	Nestler, B. and Choudhury, A. (2011) Phase-Field Modeling of Multi-Component Systems. Current Opinion in Solid State and Materials Science, 15, 93-105. https://doi.org/10.1016/j.cossms.2011.01.003
[18]	Uehara, T. (2015) Phase-Field Modeling for the Three-Dimensional Space-Filling Structure of Metal Foam Materials. Open Journal of Modelling and Simulation, 3, 120-125. https://doi.org/10.4236/ojmsi.2015.33013
[19]	Uehara, T. (2016) Numerical Simulation of a Domain-Tessellation Pattern on a Spherical Surface Using a Phase Field Model. Open Journal of Modelling and Simulation, 4, 24-33. https://doi.org/10.4236/ojmsi.2016.42003
[20]	Uehara, T., Fukui, M. and Ohno, N. (2008) Phase Field Simulations of Stress Distributions in Solidification Structures. Journal of Crystal Growth, 310, 1331-136. https://doi.org/10.1016/j.jcrysgro.2007.12.035