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], 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 Mechanics (CFD) is constantly evolving and is an important tool in engineering design.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 being simple to implement, and avoids the need for complex link-listing, as the neighboring particles or elements are fixed in space throughout the study.The disadvantage 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 [5][6][7][8].The Reynolds Transport Theorem is most often used to convert a continuity equation from Eularian to the Lagrangian domain, where ∂/∂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 [9] and other optimizing techniques.Meshless methods, however, have the advantage of not requiring a grid space for empty regions of space.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 are a minority of numerical methods in practical applications, a few have been investigated.One popular method is Smoothed Particle Hydrodynamics (SPH) [10][11][12][13].SPH utilizes a host of different smoothing functions in order to determine the magnitude of the force impacts from each neighboring particles, to discretely solve the Lagrangian Navier Stokes Equations [5][6][7]14] ρ ), where i and j represent the Einstein notation for the dimensional direction, v i (m/s) is the velocity, B i (Newtons) is the body force (such as gravity), P (Pa) is the pressure, and µ 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 incompressible fluids, using a modified Lennard-Jones (LJ) potential force to represent a fixed solid boundary [10,11,13,15].
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) [16][17][18][19].SPAM has its origins as an effort to model planetary interaction in space, as well as statistical mechanics [20].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 suffers from tensile instability [21][22][23][24][25], 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 [16,19].If these limitations could be overcome, SPH and SPAM could be used together for fluid-solid interaction (FSI) simulations, a capability that most 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 would be characterized by a smoothing function, 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 [10,15,16] is calculated as where D (Newtons) is a constant of force 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 has the best results in practical applications of SPH.The acceleration of a particle is thus where P a (Pa) is the pressure of particle a, and can be calculated by the stress vector where ii represents Einstein notation One of the biggest challenge of Lagrangian solid modeling is overcoming the tensile instability [10,16,19].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 grow 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 [16,19], the solid particles were treated similarly to liquid particles with SPH algorithms, with various techniques to avoid the tensile instability.Rather than use 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 [10][11][12].
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 linked 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 particles 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.Solid particles that have not always been fused together but comes into contact (such as during an impact) do not experience any attractive forces; only a repulsive force that only exists when the particles are within the cubic boundary.When outside of the boundary, no repulsive force is possible.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.
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 biggest difference in this approach is that the force is proportional to what is needed to stop the particle from crossing the boundary and stop it from further travel.The traditional LJ approach (Eq.7) was originally developed to describe the interactions of molecules; 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.By using the combination of these approaches, a stable Lagrangian simulation of solid mechanics and fluid solid interactions can be achieved.

Tension and Compression
A series of studies was conducted to verify this Lagrangian algorithm as a valid means of studying solid mechanics.First, an extremely simply 3-D bar was pulled in tension, to determine if the tensile stress would match Hooke's Law [8,26] where σ (Pa) is the tensile stress, E Y (Pa) is Young's modulus, and is the dimensionless tensile strain In addition, the simulation aimed to verify that the necking that occured 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 23.4 mm 3 cubic particles, with the particle dimensions of 5•11•5, pulled 0.8566 µm in tension along the axial 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, matching near identically 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 for a Y increase of 0.8566 µm, yielding an equivalent Poisson's Ratio of 0.34125; within 13% of the analytical Poisson's Ratio of 0.3.
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 0.8566 µm of deflection, and as expected, the magnitude of the stress was identical to the compression 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 for a Y increase of 0.8566 µm; the only difference was the direction.As is observable in Fig. 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 23.4 mm 3 cubic particles.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 where τ (Pa) is the shear stress, and γ is the dimensionless shear strain where δL (m) is the tangential deflection, and H is the height 90 • of the object 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.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
Now that this Lagrangian solid modeling effort was demonstrated effective in tension, compression, and shear of a simple 3-D steel bar, it will be further validated by simulations of Hertz contact [27][28][29][30] between a large disk and a flat elastic plate.The 2-D model comprises of an 81 by 10 flat plate of steel, with a particle dimension of 25 cm.The 2-D disk is represented as a series of 25 particles of identical dimensions assembled to give a one-particle disk with a radius of 2 meters (Fig. 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.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.
What is the clear 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 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.As 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 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 Hertz theory would therefore predict that the maximum stress at the point of contact would be The model described was run for five simulations of forced deflections, ranging from 26 to 49 µm.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 (Eq.19).As observed in Fig. 3, the load and peak stress matched remarkably with the analytical Hertz predicted stresses and loads, and the displacement (Fig. 4) and velocity (Fig. 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-5 inches as well as 20 to 50 meter, and in all cases, the predicted load and maximum stress always matched the Hertz predicted values remarkably.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 (Fig. 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).
The simulation was ran for 1000 time steps, averaging 146.8 nanoseconds in duration.The model was stable, with little 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 all pushed outwards (Fig. 7-a).The liquid particles fluctuated, as they were trapped within the pressure vessel, evenly shifting up and down radially as the liquid settled in the pressure vessel (Fig. 7-b).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.

Conclusion
This effort demonstrated a stable method to numerically model solid mechanics with a Lagrangian specification.Rather than using Eulerian mechanics which utilizes 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 reduced computational expense for models with large regions of open space.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 using a Lagrangian solid mechanics method to study fluid solid interactions simultaneously with the modeling of the fluid and the solid, which often is 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

PreprintsFigure 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 δ.

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