Comparing the force due to the Lennard-Jones potential and the Coulomb force in the SPH Method

D.A. Bros F.P. Piccoli

a Octopus Dofleini: En vironmental Physics/Laboratory of Simulation of Flows with Free Surface - Federal University of Espírito Santo, Brazil

b Laboratory of Simulation of Flows with Free Surface - Federal University of Espírito Santo, Brazil

Abstract Seventy years after the Smoothed Particle Hydrodynamics (SPH) method was created, there was the need to use the solid limit to keep the domain particles of a certain region. This boundary's treatment is usually done with the use of force derived from potential Lennard-Jones.In our paper, a new boundary treatment technique is presented. Now, we use the Coulomb force to exert a repulsive force on the border to avoid outgoing particles leaving the simulation domain. To test this new technique, the problem of dam break and sloshing was simulated.

Keywords: Repulsive force; SPH; Coulomb Force; Lennard-Jones.

1.Introduction

1.1. SPH background

Lucy [1] and Gingold and Monaghan [2] propounded SPH independently in 1977, which presented a technique based on the concept of kernel reproductive. It was initially intended to do astrophysical phenomena simulation and has subsequently been extended to various problems within the scope of continuum mechanics, in which it is preliminarily applied to compressible fluids.

When SPH method was created, use of solid limits for keeping particles in the domain was not needed; the subsequent necessity of a boundary's treatment, as a result of the increase of the SPH's applications, was fulfilled with the use of a force derived from the L -J potential.

SPH method was originally developed in the study of stellar interactions with domains of infinite dimensions. Therefore, the boundary is not a crucial feature in astrophysical applications. Afterwards, the SPH evolved to the study of confined flows, leading to the need of creating boundary treatment techniques, since in these problems the boundaries actively influence the domain solution. Among these techniques,the force resultant from the L -J potential, in which the domain particles are prevented from leaving the domain, is widely used. However, in this force, the variables are hard to calibrate, since the solution of these variables involves the solution of 4th order polynomials or more. Therefore, in this research, we have implemented a new technique, the Coulomb Force Boundary Treatment, in which the variables are solved through linear analysis, decreasing the computation cost. To evaluate this technique, the sloshing and the dam break tests were used.

1.2. Lennard-Jones and Coulomb force

Most Smoothed Particle Hydrodynamics (SPH) models usually treat the domain boundary with the use of virtual particles (or ghost particles). These particles are "placed" on the domain limits. Thus, a repulsive force is attributed to them to avoid interpenetration. Classically, the repulsive force used is derived from the Lennard-Jones' (L-J) potential, widely employed in molecular dynamics.

It is important to note that although L -J potential naturally has a physical explanation, the nature of this force is purely numerical during the boundary treatment in SPH. It is due to L -J potential being a mathematical model that represents the interaction between neutral molecules, which are composed of the same substance [3] . However, in our case, we are using a liquid interacting with a solid, wherein real situations, the specific masses and viscosities and other properties between the contour and fluid are different. Therefore, the molecules of the solid and the liquid do not physically correlate with this type of energy.

Eq. (1a -c ) describes the force (Y) derived from the L -J potential. The coefficientDis the well potential or the depth of the well potential andr0is the (finite) distance at which the potential between two molecules is null. The potentialrab=corresponds to the stable equilibrium position.The negative energies correspond to bound atoms forming molecules (discrete energies) whereas positive energies form a continuum and correspond to dissociated molecules (separate atoms) [3] . The termp1represents the repulsion, while the negative termp2describes the attraction.

Thus, there is difficulty in implementing the L -J method in calibrating the parametersDandp. Although Liu and Liu[5] state that this technique is widely used in molecular dynamics, we have not yet found a direct study of these variables in the SPH method. Monagham [4] has presented a relationshipD=gH/h, whereHandhare, respectively, the height of the fluid column and the smoothing length. However, this expression cannot be used in configurations other than that presented in its work for the problem of a dam break. A simple test for the cavity problem can show that this relationship is without effect.

An alternative to the L -J force presented in this paper is the Coulomb force (FCoulomb), Eq. (2) , which, like the classical L -J technique, has the sole and exclusive purpose of keeping the particles inside the domain delimited by the virtual particles. Similar to the previous case, the polarization(or neutralization) has the same nature of the virtual particles, i.e., virtual and numerical, and therefore does not allow physical discussion.

wherekis the Coulmb's constant,QaandQbare the signed magnitudes of the charges.

The advantage over L -J lies in three points:

Fig. 1. Schematic of the configuration of the physical model of the dam breaking problem. Font: Cruchaga et al. [7] .

(1) Easy calibration and there is a physical relationship with the fluid for the case of the electric constant (despite the numerical nature). It facilitates the determination of this parameter, where the value variation is independent of the geometry of the problem. It does not occur directly in the case of L -J, that, in addition to the well potential needs calibration in some problems, the associative terms are high order polynomials

(2) is within the scope of classical mechanics, which makes the electric constant a deterministic and nonprobabilistic parameter, which is the case of the well

(3) has equivalence with the field forces of classical mechanics, and can be easily implemented.

2.Mathematical modelling

The SPH method and its basic formulation are widely disseminated in the literature [5] , so we will not worry about raising this question in this work and, to be more objective,we will deal directly with the physical-mathematical model.

The mathematical formulation applied to the sloshing problem, as most of the naturals phenomena, is based on fundamentals principles of physics and the conservation laws of classical mechanics are applied here. Fortunately, the laws governing hydrodynamic phenomena can be expressed concerning mathematical equations, which generally are partial differential equations. In this case, the equations governing the phenomenon are the conservation of mass, Eq. (3) , and conservation of momentum, Eq. (4) .

where τij= 2μ( ?ui/ ?xj+ ?uj/ ?xi) ,Pis pressure, μ is kinematic viscosity, andFiis the external force. The second term of Eq. (4) is the surface force due to pressure per unit mass,and the third term is the force due to surface shear stress (or due to viscosity) per unit mass. In the equations above,tis the time, andxiis the coordinate in the Cartesian plane or in three-dimensionali= 1 , 2, 3 .

Fig. 2. Comparison between the experimental profile and the SPH model. On the right is the simulation with the Coulomb force at the centre of the experimental flow Cruchaga et al. [6] and on the left is the simulation with the strength of L -J.

2.1. State equation of artificial compressibility

In the SPH method, an artificial compressibility technique is used to model the incompressible fluid as a quasicompressible fluid. According to Liu and Liu [5] , there are two ways to impose the incompressibility: by Eq. (5) used for cases which involve low Reynolds Numbers incompressible flows using SPH; and by Eq. (6) - Tait equation, which is applied to model free surface flows.

where β =c2ρ0/γ and γ are algebraic constant, which are obtained from the Riemman function [6] . ρ0is the reference density and βis a problem dependent parameter, which sets a limit for the maximum change of the density. In most circumstances, βcan be taken as the initial pressure [ 5 , p. 137].

Fig. 3. Geometric description of the experiment by Delorme et al. [8] , corresponding to the sloshing physical tank of the Naval Architecture Department(ETSIN), Technical University of Madrid (UPM). Font: Gotoh et al. [9] .

3.Results

The comparisons was be made with already consolidated experimental and numerical works. The first comparison is made with the experimental work of Cruchaga et al [7] .,Fig. 1 .

Fig. 2 shows the results of Coulomb force (right side) and the L -J Force (left side), both with 560 virtual particles on the boundaries. The speed of sound used was 10 m/s. The simulating fluid has the properties of water. The fluid domain consisted of 35 ×70 particles, with a water column length of 0.114 m and height of 0.228 m. It was observed that in both treatments the numerical flow profile approached the experimental flow profile.

In the second comparison, we use the work of Delorme et al. [8] ., Fig. 3 . In the previous test, we carried out a qualitative study on the flow kinematics. Now, in this second test, we evaluate the dynamic aspect ( Fig. 5 ). One blank line must be included above and below each equation. Delorme et al. [8] carried out a physical experiment to study a tank of 138,000 m 3of Liquefied Natural Gas, built at the Navantia Shipyard in Spain. A prototype on a scale of 1:50 whose dimensions are 90 × 58 × 5 cm was filled with water at 9.3 cm in height under an excitation subject for a period of 1.92 s and simulating time step ( ?t) of 5e-5s. A sinusoidal oscillatory motion is imposed on the tank ( Fig. 4 ). In this new analysis, the peaks follow the trend of the experiment; but they diverge from values for both SPH -Coulomb and SPH-LENNARD-JONES. However, Coulomb again appeared to be closer to the results of expropriation, see the first periodic cycles Fig. 6 .

4.Discussion

Both treatments presented good approximations when compared to the experimental data. However, SPH -Coulomb showed more coherence under the tendency to increase pressure peaks.

The first peak in each cycle results from the impact of the wave - acceleration. And the random behaviour can be interpreted as a highly nonlinear effect where the duration of the impact is very short, and the pressure ends up being defined by the waveform at the eminence of the shock and immediately after the collision, which results in the second peak (deceleration).

Fig. 4. Post-impact moments against the left wall. The "E" frames correspond to the experiment shown in the literature in time ?t /t = 0. 021 .

Fig. 5. (a) corresponds to the record of the pressure diagram found in Delorme et al. [8] ; (b) the boundary treatment with LENNARD-JONES force and; (c)the boundary treatment with Coulomb force.

Fig. 6. Pressure diagram with the L -J treatment (gray line) and the Coulomb force treatment (dark line).

5.Conclusion

We found that the force derived from L -J potential for boundary treatment is widely used. This method, although efficient, is difficult to calibrate and has great fluctuation when there is a change in the problem or change in the number of particles of the domain and/or contour, as we saw in the dam break. It is necessary to point out that the use of that force has only numerical character, with the unique and exclusive purpose of retaining the particles inside the domain. In this same line, we produce a new form of contour treatment using Coulomb force. Also, to present better results, as well as the morphology of the flow profile and the dynamics, the SPH -Coulomb is easy to calibrate. In this new treatment, we used the same values for their constants in all cases studied,and it is not necessary to change them when we modify the number of particles. Finally, even the flow profile approaching the physical flow, the SPH -Coulomb when tested for the flow dynamics, presented better approximations about SPH with LENNARD-JONES. However, we believe that with the appropriate values for well potential, we could achieve results with less noise.

Responsibility notice

The author(s) is (are) the only responsible for the printed material included in this paper.

Conflict of interest

This study was funded by Agência Nacional do Petróleo e Gás (ANP -Brazil) (grant number N/A). Funda??o de Amparo àPesquisa e Inova??o do Espírito Santo (FAPES) (grant number N/A). The authors declare that they have no conflict of interest.

Acknowledgements

Agência Nacional do Petróleo e Gás (ANP -Brazil), Programa Institucional da Universidade Federal do Espírito Santo em Petróleo e Gás (PRH-29-ANP) and Funda??o de Amparo àPesquisa e Inova??o do Espírito Santo (FAPES). FORA TEMER.