1. Introduction
The Eulerian view is historically the most commonly used in modelling of physical problems involving fluids and phenomena transport. This modelling requires the use of numerical methods employing grids or meshes, such as finite differences, finite volumes and finite elements, for the solution of partial differential equations (PDEs) or integral formulations [^{1}]. The greatest difficulties in using meshes are in ensuring consistency. Employing meshes can lead to a number of difficulties when dealing with fluid flow problems with free surface, mobile interfaces, complex geometries and topological changes. If there are these features, generating a mesh, a prerequisite for the numerical simulation of quality, becomes a difficult process, both timeconsuming and expensive.
Lagrangian modelling is based on the hypothesis that the problem domain can be divided into a finite number of particles that do not interact with one another (oneparticle model). Each particle center of mass receives the coordinates that define its position in space and the physical properties of each Lagrangian element are found at each instance in time. The Lagrangian particle model is meshless and has been an alternative applied in research, pointing to the use of more effective computational methods for solving complex problems, providing stable and accurate numerical solutions for PDEs or integral equations.
Some of the advantages of the use of Lagrangian modelling regarding the Eulerian are the simplicity and lower computational cost in complex geometries, not needing the use of meshes, nonproduction of numerical oscillations, conservation of mass at the Lagrangian element, and finally the graphical visualisation of the results, that allows a better understanding of the spatiotemporal evolution of the fluid and its properties.
Until a few years ago, Lagrangian particle modelling was not used in the simulation of problems in fluid mechanics and transport phenomena. In 2011, Pritchard ([^{2}]) reported that describing the movements of individual fluid particles would be impractical, advising the use of the control volume analysis in modelling for fluids. For a long time, Lagrangian modelling was put aside in favour of Eulerian methods due mainly to the processing time: in general, higher than that required by mesh methods.
Currently, however, the progress of the processing techniques, the rapid spread of multicore processors, and parallelisation, made the use of particle methods possible. Graphics processing units (GPUs), using Compute Unified Device Architecture (CUDA), developed by NVIDIA, allows the use of more than one million particles in the discretisation of the fluid domain [^{3}].
Most graduate and undergraduate students still study physical laws according to the Eulerian view, basically due to the academic domain of numerical methods that use meshes to solve the differential equations used in problem modelling. However, nowadays, the meshless particle methods have become an alternative for the solution of the physical problems, with a promising future, due mainly to the advance of the processing techniques. In this context, the learning of the Lagrangian modelling and numerical particle methods is a current need in the courses of Physics. In addition, there is an increase in researches that aim at the development of Lagrangian mathematical models and new computational tools to solve physical problems.
This paper aims to present the employment of a Lagrangian particle method (Smoothed Particle Hydrodynamics) in the implementation of a computer code for solving certain physical problems. Thus, a valuable new tool for the Physics teaching will be presented to the students and researchers.
The remainder of this article is organised as follows. Section 2 brings a briefly presentation of the SPH Lagrangian method fundamentals and numerical aspects. The computer code developed is presented in Section 3. The validation of the computer code (with numerical simulations performed and results analysis) is shown in Section 4. In Section 5 conclusions are presented.
2. SPH Lagrangian Particle Method Fundamentals
SPH is based on the valid mathematical identity for a scalar function
where
By replacing the Dirac Delta function with the smoothing function,
where
The essence of the SPH method is in the discretisation of the continuum domain (
Searching for the neighbouring particles can be performed directly or using grids, which leads to a lower number of mathematical operations and reduces the computational cost ([^{4}, ^{5}]). Figure 2 shows the direct search (in which all pairs of particles ”
Different kernels can be used, and so that they are considered suitable for interpolation, it is necessary that each one follows certain properties: smoothness, positivity, symmetry, convergence, decay, compact support, and normalisation within the domain of influence [^{4}]. Figure 3 shows the kernel used by Lucy in his pioneer SPH work [^{6}], and the cubic spline kernel.
SPH approaches physical properties (such as density), gradients (pressure gradients), divergence (velocity divergence) and Laplacians (temperature and velocity Laplacians) related to fluid flow, with a secondorder error.
Table I presents the kernels employed in this work.
KERNEL 










2.1. SPH Approximations for Differential Equations in Fluid Mechanics and Transport Phenomena
Modelling of fluid flows and energy transport is performed by the conservation (mass and energy) and momentum balance equations. The evolution of the field densities, velocities and energy is defined by equations (3)(5). Table II shows the conservation and momentum balance equations, according to the Lagrangian frame of reference, and their approximations by SPH to a viscous and incompressible fluid.
Differential Equation (Continuum domain)  SPH Approximation (Domain discretised by particles)  

Mass conservation: 

(3) 
Momentum balance: 

(4) 
Energy conservation: 

(5) 
where
A complete presentation of the SPH method can be found in [^{4}]. Fluid mechanics physical laws (mathematical Lagrangian modelling), computational implementations and results achieved until this moment (summarized in the this paper) can be found in [^{8}].
2.2. Errors in the SPH approximations
The approximations to a function and its derivatives taken by the SPH method have errors. This subject will be discussed below.
Applying the expansion of the Taylor series for the function
where
Due to the kernel symmetry, equation (6) results in:
From the analysis of equation (7), it is seen that the approximation to the function value has a secondorder error. Similarly, the approximations for derivatives of the first and second orders also have errors.
2.3. Artificial Viscosity
The transformation of kinetic energy into heat takes place in problems involving mainly shock waves and needs to be measured. That energy transformation can be represented as a form of viscous dissipation, called artificial viscosity. Its application in the simulations aims to avoid numerical instabilities and the interpenetration between particles. The formulation used in the modelling of artificial viscosity is shown in equation (8), as presented in [^{4}].
where
The linear and quadratic coefficients are parameters controlling the strength of the artificial viscosity. The viscosity associated with
The term related to artificial viscosity
2.4. Particle Inconsistency
Particles near the boundaries do not have complete domains of influence, which causes truncation of the kernel and consequent lower accuracy in the interpolations [^{1}]. This phenomenon is known as particle inconsistency. Figure 4 shows the occurrence of the truncation of the kernel in the coordinate polar system.
In the border regions, as shown in Fig. 5 (b), the kernel consistency is not guaranteed. There is not a full domain of influence and the kernel is truncated. A correction method based on the Taylor series called the Corrected Smoothed Particle Method (CSPM) is used in the correction of the inconsistency [^{9}].
3. Computer Code
Figure 5 shows the flowchart of the implemented computer code implemented using FORTRAN programming language. The algorithm has been written for serial computation.
Below follows a brief presentation of the code and its routines.
Initial conditions: the initial positions, velocities, densities, temperatures, support radius and other fluid particles’ physical properties are set to the beginning of the simulation.
Boundary conditions definition: the technique to be employed for the boundary treatment is defined.
Searching for the neighbouring particles: particles within the domain of influence of a reference particle may vary over time, and the search must be performed at each numerical iteration.
Pressure calculus: in this routine, updating of the pressure field acting on the particles is performed.
Kernel definition: smoothing function used in interpolations is chosen.
Mass conservation equation solution: the particles’ density calculus for each particle is conducted from solution of equation (3).
Obtaining of surface forces: approximations to surface forces acting on the particles are obtained (the first two terms of equation (4), pressure and viscous forces, respectively).
External forces application: the last term on the right side of equation (4) is the sum of all external forces. In the class of external forces (third term on the right side of the equation (4)) is the gravity. Other forces can also be considered such as the repulsive forces exerted by the particles on the contours of fluid particles (LennardJones forces) and free surface forces, if they exist in the studied problem.
Energy conservation equation solution: SPH approximations are employed for solution of equation (5).
Momentum balance equation solution: particles’ accelerations field is obtained from solution of equation (4).
Temporal integration: prediction of the particles’ properties to the next iteration. The CourantFriedrichsLewy (CFL) stability criterion was applied in the timestep setting to ensure convergence of results [^{10}].
Output files: these files are obtained at the end of each numerical iteration and, from their data, graphical representations of the fluid's physical properties are generated.
Accuracy: this routine verifies whether the desirable accuracy has been achieved or if a new iteration will have to be executed. In the latter case, we must return to step 2 (boundary conditions definition) and follow the remaining steps until a new verification of the accuracy of numerical results be performed.
In steps 6 and 7, the correction of the particles’ physical properties (density and pressure gradient) located near the boundaries can be made using the Corrected Smoothed Particle Method (CSPM), [^{9}].
In the present work, the simulations have been performed using a computer with 8GB of RAM memory, Intel Core i52450M 2
4. Code Validation
The numerical code has been used for solution of some physical problems. Before the beginning of simulations, an examination of the physical phenomenon was carried out to select the routines to be used (Fig. 5 shows all available routines). From the numerical results analysis, made through comparisons with analytical and results reported in the literature, it has been possible to validate the computational code for application in studied cases. In this section, the set of problems, their results, and a discussion of their results are presented and discussed.
4.1. Heat Diffusion in a Homogeneous Flat Plate
The equation that governs the energy transfer in the absence of pressure fields and energy dissipation is presented below (equation (9)):
where
The heat flux is obtained by Fourier's law. For an isotropic and homogeneous material, thermal conductivity is constant, and this law is given by equation (10):
where
In the absence of heat source, in a thermal steady state, the energy conservation equation is written as follows (equation (11)):
The calculation of the SPH Laplacian of temperature for each iteration was done explicitly, in a polar coordinate system, by solving equation (12), [^{11}]:
where
The technique of the separation of variables was employed to obtain the solution by series for the temperature distribution in an isotropic, homogeneous and flat plane, at the stationary state, as shows equation (13), [^{12}].
where
In this work the first ninety terms of the series were used.
The dimensions of the plate were 1.0 m x 1.0 m. The whole domain had temperature
The particles have been placed in the domain separated by a lateral distance
Simulations were carried out starting from the transient until the steady state was reached. Different combinations of numbers of particles (50, 60, 70, 80, and 90 per side of the flat plate) / kernels have been used.
Figure 6 shows the initial distribution of 50 ×50 particles in the domain (and the prescribed temperatures at the boundaries, attributed to the virtual particles) and the temperature distribution when the steady state was reached (t
For internal regions of the domain, the approximation of the Laplacian of the temperature by the SPH method gave results that were in agreement with the analytical solution.
The largest relative errors found in the steady state (
Figure 7 shows these larger relative errors for different combinations of particle numbers and kernels.
Smaller differences in temperatures in regions close to the boundaries were observed when a kernel of higher degree was utilized (quintic spline), except for the distribution with 3600 particles. The complete analysis of simulations’ results is in [^{11}, ^{13}].
From the results analysis, we can conclude that to achieve better results it is necessary to perform a correction of the temperature values in regions near the contours (due to the kernel truncation problem).
Figure 8 shows the CPU time for different combinations number of particles/kernel employed in simulations.
Figure 9 presents the algorithm implemented in a flowchart scheme (wherein the enabled routines in the computer program is shown).
4.2. Incompressible, uniform and isothermal fluid at rest within an immobile reservoir
This is one of the first problems that must be understood and solved in Fluid Mechanics. This hydrostatic problem consists of a tank that is at rest, open to the atmosphere, and filled with a liquid. The origin of the frame of reference, with a positive vertical orientation of pointing up, is fixed at the bottom of the reservoir, as shown in Fig. 10.
The solution for a Newtonian, incompressible and isothermal fluid at rest within an immobile reservoir is well known. The application of the macroscopic view and the continuum theory (where the molecular motion of the fluid are not considered) allows the use of partial differential equations in the problem mathematical modelling. The integration of the momentum balance equation in the vertical direction leads to equation (14):
The physical laws of mass conservation and momentum balance (equations (3) and (4)), and the concept of the modified pressure, presented in [^{14}], have been employed in the problem solving. The modified pressure concept is presented mathematically by equation (15):
where
Taking into account that
where
The use of the modified pressure changes the form of the momentum balance equation (equation (4)), as shown in equation (17).
An SPH approximation for the modified momentum balance equation is:
where
The initial distribution of the particles in the domain has been performed on a regular disposition with 50 particles per side of the tank with dimensions of 1.00 ×1.00 m, or 2500 water particles in total, at 20° C and at sea level (
After the solution of equation (18) and updating from the particle positions and velocities by the numerical temporal integration, for each numerical iteration, those remained unchanged, and the hydrostatic equilibrium has been maintained. The modified pressures field were restored with 0 at the start of each new iteration, since the particles were at rest (and free surface has not changed). After starting the simulation, at the time instance 5.00 s (50,000 numerical iterations), the coincidence of the positions of the centers of mass at that moment and their initial positions was verified (regardless of the kernel used) in agreement with the analytical solution (Fig.12).
In previous literature studies, standard SPH method has been employed and results have been obtained using the momentum balance equation in its traditional form, presented in equation (4) ([^{15}^{18}]). When the standard SPH was applied to the obtaining of approximations for the pressure and viscosity forces, terms which were present in the momentum balance equation and which were added to the gravitational force, made small mistakes appear. Despite of the small orders of magnitude from these errors, they have spread and increased during the numerical simulation. The particles did not remain in their initial positions, and their movements did not cease with time progress, which is not physically correct, considering the continuum hypothesis applied for fluids (which does not take into account the molecular motion). The particle and the continuum approaches, employed in SPH studies, can agree since the microscopic fluctuations, which are absent in continuum mechanics, can be ignored [^{19}].
In this research, by using the modified pressure concept, the errors from the approximations provided by the SPH method for the physical laws (written mathematically in the form of partial differential equations) were eliminated. In consequence the continuum theory was obeyed and the particles oscillations haven't appeared.
The processing times were 6 h 37 min 6 s and 6 h 57 min 23 s, when the cubic and quintic spline kernels were used, respectively.
The flowchart, in Fig. 13, shows the enabled routines in the computer program used for solution of the hydrostatic problem.
4.3. Dam Breaking
This is a classic problem of fluid dynamics. Dam breaking studies are of great importance to the prevention of accidents, which can cause serious environmental consequences and also damage to inhabitants located in the vicinity.
The fluid is considered incompressible, uniform and isothermal; the NavierStokes equations, mass conservation and momentum balance (equations (3) and (4)), must be solved.
Numerical simulations have been performed for two different heights of dammed water (
Figure 15 presents the 2D computational domains simulated and the initial distribution of the particles.
The absolute pressure (
where
The term
where
For the Tait equation to be applied in predicting the dynamic pressure, it must be ensured that the Mach number is not greater than 0.10. In simulations, the term
The free surface particles were marked and the constant pressure condition (zero) was applied first and renewed in every numerical iteration. The time step employed was 1.00 ×10
In the boundary treatment technique, the geometric boundary conditions have used. The collisions of the fluid particles against the solid walls (considered as welldefined planes) have been studied.
SPH is a model of one particle; that is, collisions between particles are not considered. An algorithm for the study of collisions and trajectories of the particles that collide against the solid walls of the domain has been implemented. After carrying out the temporal integration and obtaining the position of the centre of mass of each fluid particle by the SPH method (in the situation in which there are no walls delimiting the domain), the collision algorithm, based on mathematical and geometry fundamentals, brings the particles back into the domain. This is based on the fact that the fluid particles could not escape the domain due to the presence of the predefined geometrical planes. The degree of elasticity of the collisions was measured by a coefficient of restitution of kinetic energy (
After the solution of equation (4), temporal integration for the position and velocity has been solved using the Euler method.
Figures 16 to 18 show the experimental and numerical results for both simulated domains.
Table III shows the percentage differences between the abscissas of the wave fronts obtained by SPH and experimentally, calculated using the equation (21):
First Domain Simulated  Second Domain Simulated  











0.20  0.10  19.00  18.07  4.89  0.20  0.10  22.00  19.00  13.64 
0.20  29.00  30.28  4.41  0.20  36.00  37.00  2.78  
0.30  0.10  19.00  17.58  7.47  0.30  0.10  22.00  18.07  17.86 
0.20  29.00  28.70  1.03  0.20  36.00  35.16  2.33 
where
From the Table III analysis, it can be seen that the greatest differences between the wave fronts occurred in the initial stages of dam breaking. With the time progress, the differences decreased, as shown at time 0.20 s. The smaller differences between wave fronts (1.03% in the first and 2.33% in the second geometry) have occurred when it was used the coefficient
In the first geometry simulations, CPU processing times were 2 h 21 min 24 s (
In Figure 19, a flowchart shows the numerical code implemented for solution of dam break problem.
5. Conclusion
Due to the advances in computer processing techniques, there has been an increase in the applications of the particle methods for solution of problems in fluid mechanics and transport phenomena areas. SPH is a particle method that has advantages compared to traditional methods using meshes, e.g. better visualisation of the spatiotemporal evolution of the flow, and a lower computational cost in the study of complex geometries with topological changes or free surfaces beyond the control of undesirable numerical diffusion. From the graphical analysis of the results, by viewing the properties of particles (which may be fields of positions, velocities, accelerations, pressures, among others) a better understanding of the physical phenomenon that occurs is reached.
The learning of the Lagrangian modelling and numerical particle methods is a current need in the Physics courses. The recent progress of processing techniques, the rapid spread of multicore processors, and parallelisation, make possible the use of particle methods for physicals problems solving. In addition, there is an increase in researches aimed at the development of Lagrangian models applicable to physical problems and the need to development and use new computational tools.
The computational code presented, and validated for the study of certain physical problems, is a valuable tool and points to the future, as an instrument for the application of particle modelling in Physics teaching.