A Discrete Approach to Meshless Lagrangian Solid Modeling

The author demonstrates a stable Lagrangian solid modeling method, tracking the interactions of solid mass particles rather than using a meshed grid. This numerical method avoids the problem of tensile instability often seen with smooth particle applied mechanics by having the solid particles apply stresses expected with Hooke’s law, as opposed to using a smoothing function for neighboring solid particles. This method has been tested successfully with a bar in tension, compression, and shear, as well as a disk compressed into a flat plate, and the numerical model consistently matched the analytical Hooke’s law as well as Hertz contact theory for all examples. The solid modeling numerical method was then built into a 2-D model of a pressure vessel, which was tested with liquid water particles under pressure and simulated with smoothed particle hydrodynamics. This simulation was stable, and demonstrated the feasibility of Lagrangian specification modeling for fluid–solid interactions.


Introduction
Computational solid modeling is an incredibly valuable tool for todays engineers [1][2][3][4][5][6], thankfully due to the powerful computers available at low costs.Finite element analysis (FEA) is just one of many meshed solid modeling techniques used in countless industries to study mechanical stresses in numerous different components today.Computational fluid dynamics (CFD) is constantly evolving and is an important tool in engineering design [7][8][9].These numerical techniques enable a detailed and robust study of complicated geometries and nonlinear behavior with far less need for costly and time-consuming experimental studies.
The vast majority of numerical methods in engineering, including finite element, are meshed, analyzing stresses and mass flows in a fixed region of space.This approach to solving continuum mechanics is considered the Eulerian specification.This has the advantage of keeping the relative distances between discrete particles and elements fairly constant.One of the main disadvantages of this method is that the entire domain needs to be modeled, and if there are large voids and empty spaces, there is a waste in computational effort as empty domains are studied repeatedly at each time step.
Another approach to avoid this wasted computational effort is to use meshless numerical methods, studying the stresses and forces in the Lagrangian specification [10][11][12][13].The Reynolds transport theorem is most often used to convert a continuity equation from the Eularian to the Lagrangian domain, where v (m/s) is the velocity vector, V (m 3 ) is the volume, F is an arbitrary function, ∂/∂t is the partial derivative, and D/Dt is the total derivative, When characterizing continuum mechanics in the Lagrangian specification, rather than studying the mass that travels in and out of a specific discrete meshed domain, the mass itself is simulated as a set of discrete particles.Each particle has its own unique location, velocity, and stresses.At every time step, every particle is affected by all neighboring particles; because of this it is necessary to calculate the changes in relative distances between individual particles.A big challenge of meshless methods is the need to constantly determine the relative distances between all of the neighboring particles, a task that grows exponentially with increasing particle quantity.This can be mitigated with link-listing and other optimizing techniques [14,15].Meshless methods, however, have the advantage of not requiring a grid space for empty regions of the domain.This can be advantageous for simulations such as large deformations and explosion studies, the study of planetary motion, or extremely small nanoparticles floating through a nanofluid.
While meshless methods form a minority of numerical methods in practical applications, a few have been investigated.One popular method is smoothed particle hydrodynamics (SPH) [16][17][18][19][20]. SPH utilizes a host of different smoothing functions in order to determine the magnitude of the force impacts from each of the neighboring particles, to discretely solve the Lagrangian Navier-Stokes equations [10][11][12]21], where v i (m/s) is the velocity in the i-direction, B i (Newtons) is the body force (such as gravity), P (Pa) is the pressure, and µ (N•s/m 2 ) is the dynamic viscosity.At every time-step, it is necessary to calculate the separation distance and smoothing function for each neighboring pair of particles, and thus the computational resources grow exponentially with particle quantity; link-listing is often used to mitigate this.SPH also only models fluids, often but not exclusively using a modified Lennard-Jones (LJ) potential force to represent a fixed solid boundary [16,17,19,22].
Attempts have previously been made to branch off SPH to apply to solid modeling; this has been referred to as smoothed particle applied mechanics (SPAM) [23][24][25][26].SPAM has its origins as an effort to model planetary interaction in space, as well as statistical mechanics [27].The Lagrangian-based method is quite practical when dealing with relatively small particles (such as planets and stars) in a large domain filled with empty space; it is computationally wasteful to model vast sums of empty space with no activity in it.At smaller scales, however, SPAM and SPH suffer from tensile instability [28][29][30][31][32], as particles under tensile stress eventually become unstable regardless of the size of time-integration steps.Efforts to use artificial viscosity approaches have only had limited success [23,26].If these limitations could be overcome, SPH and SPAM could be used together for fluid-solid interaction (FSI) simulations [33][34][35][36][37][38][39], a capability that many commercial FEA and CFD software packages lack.
This effort is to investigate meshless modeling of solid mechanics, with the goals of eventually evolving into a practical SPAM-like tool for modeling solid mechanics.The eventual goal is to have a tool that can interact with SPH fluid particles, and perhaps become a useful FSI tool for fluid flows over a large domain.With an accurate model consisting of particles of liquids and solids interacting, an engineer could investigate eventual material failures and crack propagation in real time for highly dynamic fluid flow within a moving solid boundary.

Mechanics
In a Lagrangian solid modeling algorithm such as SPAM, each neighboring particle's relative distance would be defined by a smoothing function W (m −3 ), such as, where λ is the dimensionless ratio of the absolute distance between two specific particles (m) and the smoothing length h (m), The density is calculated as, where ρ a (kg/m 3 ) is the density of particle a, and m b (kg) is the mass of neighboring particle b.
The Lennard-Jones potential [16,22,23] acceleration δv LJ /δt (m/s 2 ) is calculated as, where D (m/s 2 ) is a constant of acceleration proportional to the particle velocity squared, r 0 (m) is the specified width of a solid particle, r ab (m) is the total distance between particle a and b, and x ab (m) is the directional distance between particle a and b.The values of M and N are arbitrary coefficients; M = 12 and N = 4 often have the best results in practical applications of SPH [16,17].
The acceleration δv a /δt (m/s 2 ) of a particle a by a neighboring particle b is thus, where P a (Pa) is the pressure of particle a, and can be calculated by the stress tensor σ ij (Pa), and ∆W ab (m −1 ) is the gradient of the normalized smoothing function between particles a and b.
One of the biggest challenges of Lagrangian solid modeling is overcoming the tensile instability [16,23,26].If the fluid or solid is under pressure (P > 0), these Lagrangian specification methods work well.If the solid is under tension, however, where (P < 0) and the particles of mass are attracted to each other, there is a tendency for all of the mass particles to clump together due to the fact that the derivative of most smoothing functions grows exponentially as the particles become closer together.This is a challenge to Lagrangian modeling of solid mechanics that must be overcome.

Algorithm
In previous studies [23,26], the solid particles were treated similarly to liquid particles with SPH algorithms, with various techniques to avoid the tensile instability.Rather than using the traditional smoothing algorithms seen in SPAM, this model works by applying different steps for liquid-liquid, solid-solid, and liquid-solid particle interactions.All of the particles move freely (unless specified as fixed) within the domain, but accelerations and stress calculations are specific to the different classes of particles.The liquid-liquid particle interactions were all studied using the established Lagrangian CFD method of SPH [16][17][18].
The solid-solid particle interactions, however, used a much different approach from traditional SPH.The solid particles are given a specific cubic shape (that undergoes elastic strain), and particles of mass connected in tension are linked together in a specific contact matrix.Only particles linked together in this matrix can experience tension stress with a displacement apart from each other.Non-linked particles exert no stress on each other when they are proximate to each other unless they are close enough to be within the particle's shaped boundary; in this circumstance they are under compressive stress.While this bears similarities to a meshed approach, it is not Eulerian as the contact matrix merely relates to each mass particle which can travel freely within the domain.
This approach avoids the concerns of the tensile instability entirely.Some of the solid particles are not fused together, and thus have no attractive forces; they only experience a repulsive force when in close proximity to each other (such as during an impact).For particles pre-defined as in contact and fused together, they only experience a tensile force when pulled apart; this force reverses itself entirely when the particles are close enough that they cross each other's solid boundaries.By using this approach, the issue of tensile instability, where the attractive forces of particles in tension grow exponentially and result in solid particle clumping, is avoided entirely.
This algorithm models solid-solid forces with a numerical approach to Hooke's law [13,40] where σ (Pa) is the tensile stress, E Y (Pa) is Young's modulus, and is the dimensionless tensile strain, where δL (m) is the material deflection of an object with a natural length of L (m).This strain can occur from particle pairs that are in both compressive and tensile contact.Each solid particle represents a cubic particle of mass, with a predetermined length, width, and height of δx (m).If the length between two particles δx ab (m) is ever less then the particle length, the particle is then under a compressive acceleration repulsing the two particles apart; if particles linked in contact are ever farther apart than the particle length, then an attractive acceleration is applied to each particle.While the particle cubic dimension δx is constant, the true particle shape is rarely a cube, as it is subjected to strain.When the stress is determined, the particle strain is calculated, and thus the particle dimension δx u in the u direction, where E uu is the component of the strain tensor, where σ uu (Pa) is from the stress tensor.
The equivalent strain between two particles is thus, which can be used in Hooke's law to find the equivalent acceleration δv ab (m/s 2 ) of particle a due to compressive or tensile forces from particle b, where m a (kg) is the mass of the particle a, and v and w represent dimensions that are not u; the dimension of contact between the two particles being u.
When there is liquid-solid particle interaction, the particles exert a repulsive force on each other not dissimilar to Lennard Jones with some unique differences.The traditional LJ approach (Equation ( 7)) was originally developed to describe and model the interactions of molecules subjected to Van der Waal forces; at times using this approach for macroscopic fluid solid interactions results in liquid particles passing through a boundary of solid particles, and at other times the exponential nature of the equation results in unrealistic repulsive accelerations.Rather than use LJ, this effort approaches this problem by applying a force on the particles in contact proportional to what is needed to stop the particle from crossing the boundary and stop it from further travel.A variation of this approach was previously taken with SPH and discrete multi-hybrid physics (DMP) models to study the effects of blood and tissues in the human body [33][34][35][36][37][38][39].By using the combination of these approaches, a stable Lagrangian simulation of solid mechanics and fluid-solid interactions can be achieved.
The first step in this algorithm is to establish the model; this is a list of unique mass particles with their own locations, velocities, stresses, and (in the case of solids) strains.The next approach is to refresh the link-list; particles move and therefore the link-list database needs to be updated regularly to avoid errors.The link-list, by definition, is a quasi-grid that segregates particles into neighboring regions, so computer resources are not wasted calculating the impacts of neighboring particles far too separated to have a meaningful impact.The next step is to determine how much acceleration each particle can expect as a result of the tension (from linked solid particles), impacts, repulsive forces from solid-solid interactions, liquid accelerations based on SPH forces defined in Equation ( 9), and liquid-solid interactions.When the accelerations are determined, the velocities are updated to reflect such an acceleration over the time step, followed by a change in location depending on the individual velocity and time step duration of each particle.The final step of a given time step is to recalculate the new tensions, stress, and strains (for solid particles) based on the new particle locations.

Tension and Compression
A series of studies was conducted to verify this Lagrangian algorithm as a valid means of studying solid mechanics.First, a simple 3-D bar was pulled in tension, to determine if the tensile stress would match Hooke's law (Equation ( 11)).In addition, the simulation aimed to verify that the necking that occurred would properly match the Poisson's ratio ν, where L (m) and d (m) are dimensions of the solid under stress.
For this case study, the material parameters of steel will be used.For steel, the Young's modulus is E Y = 207 GPa, and the Poisson's ratio is ν = 0.3.The 3-D steel block was represented by a series of cubic particles, 23.4 mm in length per particle, with the particle dimensions of 5 × 11 × 5 in the x, y, and z direction, respectively, for a total length of 25 cm in the y direction.The steel bar was pulled 1.7132 µm in tension along the axial y direction.All of the particles were in tensile contact with their nearest neighbor.The peak tensile stress of 1.4758 MPa was observed at the center of the bar, identically matching to the analytical stress predicted with Hooke's law.In addition, the width in the X and Z direction decreased at a distance of 0.1329 µm, yielding an equivalent Poisson's ratio of 0.34125, well within the numerical error that can be expected with a 5 × 5 cross-section resolution.
This study was then reversed for the model, but in compression rather than tension.The 5 × 11 × 5 block with a Young's modulus of 207 GPa and a Poisson's ratio of 0.3 was compressed to have 1.7132 µm of deflection, and as expected, the magnitude of the stress was identical to the tension case, with 1.4758 MPa of compression stress.In addition, the width in the X and Z direction increased the exact same distance of 0.1329 µm; the only difference was the direction.As is observable in Figure 1, the tensile pulling of the bar resulted necking in the bar in the X and Z direction, whereas compression results in a bulging in the X and Z direction.

Shear
The next step in the validation of this numerical model would be to test it in shear.This model would also be in steel, with a 11 × 3 × 3 block of cubic particles, 23.4 mm in length, for an overall length of 25 cm.With the Young's modulus E Y of 207 GPa and a Poisson's ratio of 0.3, the shear modulus G Y (Pa) can be found simply as, which ultimately results in a shear modulus of G Y = 79.6154GPa.The shear modulus can be used to find the shear stress as a linear function of shear strain, where τ (Pa) is the shear stress, and γ is the dimensionless shear strain, where δL (m) is the tangential deflection, and H is the height of the object 90 • tangent to the deflection.
The model used was subjected to a shear deflection of 18.3429 µm in the X direction.The dimensionless shear strain γ is found by taking the deflection over the height (9.3619 cm) of the bar (Equation ( 20)) and the shear strain can be used in Equation ( 19) to find the shear stress τ, τ = (79.6154×10 9)×(1.9593×10−4 ) = 15.5992×10 6 .
The numerical maximum shear stress was observed to be 15.8993MPa, an error of less than 2%.This close match further validates this Lagrangian numerical method as a reasonable approach to solving simple solid mechanics.

Hertz Contact Simulation
After this Lagrangian solid modeling approach was demonstrated effective in tension, compression, and shear of a simple 3-D steel bar, it was then further validated by simulations of Hertz contact [41][42][43][44] between a large disk and a flat elastic plate.The 2-D model comprises a 51 by 10 flat plate of steel, with a particle dimension of 159 µm, to form a 8-mm by 1.5-mm plate.The 2-D disk is represented as a series of 53 µm particles of identical relative dimensions assembled to give a disk, three particles thick, with a radius of 12.7 mm (Figure 2).This disk will experience no elastic deflection due to having an infinite Young's modulus; the disk would be forced down up to 49 µm into the elastic plate in five discrete increments.The elastic plate particles are free to move, except for the final bottom row (opposite side of the Hertz contact); these particles are fixed.The model will be validated by comparing the deflection and stresses in this elastic plate, and comparing the numerical results to analytical Hertz predictions.The control parameter in this simulation is the actual deflection down at the center of the point of contact δ (m); this deflection can be used within the Hertz contact model to determine the maximum stress.The first step is to determine the reduced Young's modulus E Y (Pa) and reduced radius R'.As the disk does not deflect, it is treated as an infinitely stiff material E Y = ∞, and thus, where E Y (Pa) is the Young's modulus, and ν is the dimensionless Poisson's ratio.The disk is infinitely stiff, E Y,disk = ∞, and thus the reduced Young's modulus (in MPa) is, The reduced radius is easily found as, 1 Because the disk is of infinite stiffness and does not deflect elastically, the width a (m) of the contact region can be found with simple trigonometry, where δ (m) is the deflection distance.Hertz theory would therefore predict that the maximum stress at the point of contact would be, One significant challenge of using the Lagrangian methods to model Hertz contact is to make sure the time step is short enough so that the model is stable.The maximum time step δt max (s) allowed in these simulations, where δx (m) is the minimum particle width, and C s (m/s) is the speed of sound; this approach has been used in successful applications of SPH [16,17].This is the longest permitted time step; to better ensure stability, a variable time step δt (s) was determined at each time step, where V max (m/s) is the maximum particle velocity at a given time step.The benefit of Equation ( 27) is that it ensures that at no point will a particle ever move further than 10% of a particle-width in a given time-step, ensuring the model remains stable.The value of δt is found with each time step, and if it is ever greater than the value of δt max found with Equation ( 26), then δt max is used as the time step.While this approach is extremely stable, it also results in extremely short time steps; the simulation under discussion uses time steps ranging from 15 to 56 nanoseconds.The model described was run for five simulations of forced deflections, ranging from 26 to 49 µm.In a given simulation, the disk would be forced down a given increment, and the plate particles would have time to settle; after they settled the disk would be forced down to the next displacement increment.According to Hertz contact mechanics equation, the analytical load to accomplish this would be 2.66 to 5.98 kN, with the peak stress ranging from 7.3 to 10 GPa (Equation ( 25)).As observed in Figure 3, the load and peak stress matched remarkably with the analytical Hertz predicted stresses and loads, and the displacement (Figure 4) and velocity (Figure 5) all settled into place after a period of time.The error percentages for the load varied from 2% to 10.4%; the error percentages for the peak stress were even lower, ranging from 0.2% to 5.4%.This close match was repeated for models of identical dimensions, with small disk radii ranging from 2 to 5 inches as well as 20 to 50 meters, and in all cases, the predicted load and maximum stress always matched the Hertz predicted values remarkably (Figure 3).This Hertz study strongly demonstrates the feasibility of this Lagrangian numerical model to predict the stresses and strains in a solid model.

Pressure Vessel
The last step in this effort is to join the Lagrangian specification solid modeling efforts with smoothed particle hydrodynamics; this was demonstrated with a 2-D pressure vessel.A model of a 2-D circular pressure vessel was built; the solid pressure vessel was five particles thick, whereas the liquid was comprised of a 2-D circular region with a radius of twenty particles (Figure 6).Each particle, both fluid and solid, had a particle length of 23.4 mm.The solid particles were steel, with the same parameters as the earlier studies; the Young's modulus is E Y = 207 GPa, the Poisson's ratio is ν = 0.3, and the default density ρ = 7800 kg/m 3 .The liquid water has a default density of ρ = 1000 kg/m 3 , and a bulk modulus of K bulk = 2.15 GPa.The water is initially set in the pressure vessel at a pressure of 1 atmosphere (101,135 Pa).One fundamental distinction was the effective shape of the particles; the solid particles were cubic (as was the case in all the prior examples), whereas the liquid particles were effectively spherical.The reason for this is that the solid particles undergo shear in all three dimensions, and the cubic approach is necessary to accurately determine the tensions between fused solid particles.The liquid particles are shapeless (liquid inherently has no shape), but as SPH uses the smoothing length (Equation ( 4)), the impacts of neighboring particles are all proportional to the radial separation distance and the smoothing length.
The simulation was ran for 400 time steps, averaging 146.8 nanoseconds in duration.The model was stable, with few dramatic shifts in particle position; this is expected as in the pressure vessel, only elastic expansions are expected as the particles settle into place.It was observed that the steel particles are all pushed outwards (Figure 7a) as the liquid pressure stabilizes within the pressure vessel (Figure 7b).This model demonstrated the feasibility of merging the Lagrangian liquid numerical method of SPH with Lagrangian solid modeling, with applications in the study of fluid-solid interactions.

Conclusions
This effort demonstrated a stable method to numerically model solid mechanics with a Lagrangian specification.Rather than using Eulerian mechanics which utilize a meshed grid, this model used a Lagrangian approach that tracked individual particles of the mass, and how these particles interacted with each other.This approach has the advantage of not being confined to a fixed domain, as well as having reduced computational expense for models with large regions of open space, and has models with large deformations.A solid bar in compression, tension, and shear was simulated, and it matched the stresses and strains remarkably with the analytical results expected with Hooke's law.A much more detailed simulation of a rigid disk being applied to a flat plate was then studied, and the model matched the analytical loads, strains, and stresses expected with Hertz contact theory; the Hertz model matched the numerical model for a large host of loads, particle dimensions, and disk sizes.Finally, this Lagrangian solid model was implemented with the Lagrangian computational fluid dynamics method of smooth particle hydrodynamics to build a stable model of water under pressure in a 2-D pressure vessel.This demonstrated the ability of the Lagrangian solid mechanics method to study fluid-solid interactions simultaneously with the modeling of the fluid and the solid, which are often modeled separately.By using this model, a host of solid mechanics applications can be better modeled, resulting in overall better engineering design.

Figure 1 .
Figure 1.The steel bar (a) with no stress, (b) in tension, and in (c) compression, with 10 5 exaggerated deflection.

Figure 3 .
Figure 3.Comparison of Hertz and numerical solid modeling, for both (a) the numerically integrated total load, and (b) the maximum stress, calculated as a function for a given displacement δ.

Figure 4 .
Figure 4.The displacement, where each line represents the displacement of a specific mass particle in the (a) y and (b) x direction, for the top five center-most particles.The black dashed line in (a) represents the forced displacement of the infinitely-stiff disk.

Figure 5 .
Figure 5.The velocity change, where each line represents the velocity of a specific mass particle in the (a) y and (b) x direction, for the top five center-most particles.

Figure 6 .
Figure 6.The pressure vessel, at the start of the simulation.Green circles represent water particles, and blue dots represent steel particles.

Figure 7 .
Figure 7.The pressure vessel steel particle radial movement as a function of time, both for the (a) solid particles and the (b) liquid water particles.