CC BY
21
УДК 519.7
Метод сглаженных частиц на основе потенциалов для оценки динамического разрушения компактированных и гранулированных материалов
R. Brighenti, N. Corbari
Пармский университет, Парма, 43100, Италия
Компьютерное моделирование твердых тел получило интенсивное развитие с появлением метода конечных элементов. Однако данный метод не позволяет адекватно решить ряд задач из-за сильного искажения сетки в расчетах по лагранжевой схеме в условиях больших смещений, высокоскоростного нагружения, фрагментации, измельченной твердой фазы, при взаимодействии жидкости с конструкцией и т.д. Дискретизация рассматриваемой области с помощью несвязанных узлов сетки может служить решением описанной проблемы. Кроме того, дискретная природа сплошной среды (наблюдаемая на микроуровне) позволяет применить дискретизацию, типичную для гранулированных материалов, что в свою очередь позволяет моделировать большие деформации и разрушение (растрескивание, измельчение, фрагментацию, образование кластеров) за счет изменений во взаимодействии между частицами.
Метод сглаженных частиц представляет собой бессеточный метод частиц в лагранжевой формулировке, широко применяемый во многих инженерных областях. В статье в рамках методов сглаженных частиц предложен унифицированный метод частиц на основе потенциалов для механического моделирования сплошных сред и гранулированных материалов в условиях динамического нагружения. Взаимодействие частица-частица и частица-граница моделируется с помощью силовых операторов, обусловленных природой изучаемого материала (сплошной, зернистый и т.д.). Показана возможность моделирования существенных изменений в геометрии механической системы, например, растрескивание, образование кластеров или гранулированный поток. На примерах демонстрируются возможности предлагаемого подхода.
Ключевые слова: метод частиц, метод сглаженных частиц, контакт, сплошная среда, гранулированная среда
A potential-based smoothed particle hydrodynamics approach for the dynamic failure assessment of compact and granular materials
R. Brighenti and N. Corbari
University of Parma, Parma, 43100, Italy
Computational simulation of solids has experienced a rapid development since the formulation of the finite element method. However a number of problems cannot be properly solved by using the finite element method because a severe mesh distortion in computations of Lagrangian scheme may arise for very large displacements, high speed impact, fragmentation, particulate solids, fluid-structure interaction, leading to lack of consistency between the numerical and the physical problem. The discretization of the problem domain with nodal points without any mesh connectivity would be useful to overcome this difficulty. Moreover the discrete nature of continuum matter— usually observed at the microscale—allows to adopt such a kind of discretization that is natural for granular materials and enables us to model very large deformations, handle damage—such as fracture, crushing, fragmentation, clustering—thanks to the variable interaction between particles.
In the context of meshless methods, smoothed particle hydrodynamics (SPH) is a meshfree particle method based on Lagrangian formulation that has been widely applied to different engineering fields. In the present paper a unified computational potential-based particle method for the mechanical simulation of continuum and granular materials under dynamic condition, is proposed and framed in the SPH-like approaches. The particle-particle and particle-boundary interaction is modelled through force functionals related to the nature of the material being analyzed (solid, granular, ...); large geometrical changes of the mechanical system, such as fracture, clustering, granular flow can be easily modelled. Some examples are finally proposed and discussed to underline the potentiality of the approach.
Keywords: particle method, smoothed particle hydrodynamics (SPH), contact, continuum solids, granular media
Nomenclature
a (rinfl) —function of the influence radius of a particle; Ay, AQy —cross section area of the truss assumed between
particles i and j and its reference value, respectively; bi —body force component in the i-th direction; C'—tangent elastic tensor of the material;
© Brighenti R., Corbari N., 2015
c(s)—function accounting for the no-penetration condition; d* —equivalent diameter of two particles in contact; dt = 2ri —diameter of the generic particle i; E'(s), E0 —elastic modulus of a linear element that represents the material connecting the two particles i, j and corresponding reference value for s >> 0; Ei —elastic modulus of the generic particle i; E —equivalent Young modulus of the elastic contact between two particles;
Etot(x) = n(x) —total energy of the particle system; Fij (s) —generic force acting between a couple of particles at the effective distance s;
F;, Fd, Fb, Fe —internal force, damping force, boundary and
external force vectors, respectively;
Ftot i —vector of the total force acting on the particle i;
h—smoothing length or support dimension;
K(s) —stiffness of the particles bonding depending on their
effective distance s;
Kn —normal stiffness of the particle-boundary contact; mi, M —mass of the particle i and mass matrix, respectively;
n, t—unit vectors normal and parallel to the tangential plane in the contact point between a generic particle and the boundary surface, respectively; Pi —force vector acting on particle i; q—unit vector identifying the direction connecting the two particles centers;
r—distance between the centers of a couple of particles; r —average radius of a couple of particles; rinfl —radius of influence (or maximum interacting distance) of a particle r^;
r0 —distance between the centers of the particles at which
F (r = ro)
s, se —effective distance between particle surfaces and at the equilibrium state, respectively; s' = r - (dj 2 + dj J 2) —distance between two particle surfaces;
t, At—time variable and time integration step amplitude, respectively;
Tn, Tt —particle-boundary contact surface forces normal and parallel to the tangential plane at the contact point, respectively;
Tt ij —tangential force between the colliding particles i, j; Vp —volume of the particle p;
w—displaced distance of one particle into another or into the contact surface;
W (| x - xi |) —Kernel or smoothing function for the smoothed particle hydrodynamics method; x0 —position vector identifying the equilibrium state; xi, xi, xi —position, velocity and acceleration vector of the particle i, respectively;
x, x, x —position, velocity and acceleration vector, respectively;
a—coefficient defining the maximum copenetration depth e j (x;-) —strain tensor components at the point identified by xi;
Y—thickness of a soft layer added to the boundary surface
in order to smooth the contact forces;
S(| x - xi |) —Dirac delta function placed at xi;
8 = ar —maximum copenetration amount between two
particles;
Oy —stress tensor components;
®(x), ®tot(x) —generic strain energy potential and potential of the particle system, respectively; Ad —damping coefficient; n—viscosity coefficient;
X(w)—smoothing function for the force particle-boundary contact evaluation;
—coefficient of dynamic friction between particles and boundaries;
M-md —coefficient of dynamic friction between particles; vi —Poissons ratio of the generic particle i; p—mass density.
1. Introduction
Numerical modelling of engineering problems, specially the mechanical ones, has attracted the interest of several researchers aimed at assessing the response of materials both during service and at their ultimate conditions. To this end plenty of computational tools have been developed, starting from the most diffused finite element method appeared in the middle of the XX century [1, 2], the finite difference method [3], the more recent boundary element [4], meshless [5] and natural element methods [6], just to cite the mains.
Problems dealing with mechanical systems under particular operating conditions, such as large displacements, high velocity impact, fracture, failure, fragmentation, clustering, fluid-solid interaction, and so on, typically do not allow a simple and reliable solution through classical numerical techniques in the context of the so-called Lagrangian approach, because the lack of consistency between the numerical and the physical problem due to the severe distortion of the discretized domain. This problem has been partially solved by introducing special remeshing techniques within the finite element method, even if the remeshing itself is not simple and moreover requires a proper mapping of the state variables from the old mesh to the new one, with the consequence of introduction of numerical errors.
The description of the domain of the problem with only nodal points but without the necessity of mesh connectivity, is a possible approach to overcome this drawback. In this context the so-called smoothed particle hydrodynamics has been one of the first particle meshless methods in computational mechanics [7, 8]. The original idea of the pioneering developers lay in the fact that physically the collective motion of particles in particle-wise systems (such
as the formation and evolution of stars or galaxies, the movement of a liquid, gas flow, particulate assemblies, etc.) may be all modelled by the governing equations of classical Newtonian hydrodynamics.
Such mesh-free Lagrangian method is characterised by variable nodal connectivity, where the material is discretized with particles interacting each other with interpolation functions. It operates by transforming a discrete field into a continuum one by using a localized kernel function that operates as a smoothing interpolation. The method can be naturally applied to discrete systems but also to continuum ones, where the nodal points are not directly associated to the single particle composing the matter [9-14]. The main fundament of the method is to get the system of discrete algebraic governing equations of a continuum problem, by localizing the partial differential equations through the use of a convoluted integration.
The noticeable advantages of such a method became suddenly evident after its first introduction and was subsequently extended and adopted for solving a wide range of physical and engineering problems, well beyond the initial classical hydrodynamics field. More recent researches have considered and extended smoothed particle hydrodynamics approach, such as the meshfree method, for solving various engineering problems [15-20].
If the discrete nature of solid matters can be recognised at the microscale or at the molecular level [21-23], materials at the macroscale are typically modelled as continuous media, suggesting the use of continuous computational techniques such as the finite element method. However discrete models can be conveniently applied even to macroscopi-cally continuum systems, where each representative particle can be associated to a cluster of continuum elements; discrete element methods or particle methods, usually identified as molecular dynamics methods at the nanoscale level, are numerical procedures for the solution of different engineering problems. The fundamental assumption for this class of methods consists in the description of the material through an assemblage of separate discrete elements, such as atoms, molecules, grains, particle solid elements, etc. [2427].
According to the discrete approach, complex nonlinear interactions between bodies and within bodies, are numerically simulated in order to get the motion of particles governed by nonlinear differential equations. Such discrete nature of solids is well evident for some classes of materials such as the granular ones that are constituted by several low deformable particles, usually interacting each other through elastic contact, cohesive and friction forces. Among the different problems requiring the mechanical simulation of materials, geomechanical and powders ones can be naturally studied by exploiting their well evident granular nature [24, 28-31].
Discrete methods have relevant applications in macroscopic scale simulations such as mineral processing, rock
blasting, crushing, sand mechanics, mechanical behavior and failure of compact or granular bodies [32-34]. Moreover materials like gases and liquids can also be simulated with the simple scheme of interacting particles [35, 36].
On the other hand, a generic solid body can be assumed to have a particle structure with a proper nature of the particle interaction forces. From this observation the mechanical behaviour of a material can be made to range from the very incoherent cases (characteristic of fluids, granular, powders, etc.) up to the compact materials ones (typical of polycrystalline materials), respectively.
The mechanical description of a body can be made through the cinematic and static properties of its particles, such as the displacements, velocity, accelerations, internal forces, stresses and so on: in a fluid matter the mechanical response is obtained by following the motion of interacting particles enclosing the physical properties of the flow, while in a continuum these properties are usually macroscopically averaged, e.g. the fields of density, momentum, etc. over a small reference volume are of main interest.
The above observations suggest the potentiality to use a discrete model for the simulation of different class of solids by properly choosing the nature of the interaction forces between their discrete constituent elements. Such an approach allows us to get the overall response of the material at the macroscale, which is the main interest of materials science and mechanics of materials. Moreover the particle discretization for a solid naturally allows to tackle the problem from a dynamic point of view, enabling the solution of high strain rate, impact, large displacements and large strains problems involving history-dependent behaviour, plasticity, etc. that are more naturally expressed in a Lagrangian computational framework.
As mentioned above, in problems involving large deformations, a purely Lagrangian approach applied in the framework of continuum mesh-based computational techniques cannot be appropriate due the requirement of complex remeshing and smoothing, too small time steps, and so on; all the above drawbacks have promoted the development in the last decade of particle methods in solid mechanics, achieving now a wide popularity. Moreover the possibility to deal with short and long range forces (such as the electrostatic ones), meets the physical nature of the materials at small or nanoscales, where interatomic forces obeying such no local interaction exist.
In the present paper, a discrete element approach for the mechanical simulation of both continuum and granular-like materials—formulated on the basis of the force potential interaction concept—is presented in the framework of the smoothed particle hydrodynamics theory.
After an introductory part related to the forces existing in materials characterised by a macroscopic particle nature, a discrete approach for continuum solids and discrete incoherent aggregates is presented by considering the dynamic nature and large strain characteristics of the problem. Fi-
nally, some examples demonstrating the capability to deal with the dynamic failure of solids and granular flow problems are presented to show the versatility and generality of the proposed approach.
2. Particle model for solid materials
According to the introductory section, the mechanics of continuous materials can be described through a particle model, once the discrete element interacting forces have been properly described.
As recalled above, such an approach is straightforward when the material is effectively composed by several discrete elements as at the nano/microscale for solids and for granular-like matters, where the particles can be immediately identified as the atoms/molecules or grains present in the solid. Continuous materials can also be idealized as a cluster of elements reciprocally joined by forces that modify when the particles relative position changes; if the particles tend to approach each other a repulsive force arises (such as in compact solids and granular materials), while an attractive force appears when they tend to move far away (such as in compact solids—typically characterised by elastic-like forces—or in cohesion granular matters).
The above recalled forces can be thought to have a local nature when they are associated to the contact between elements (for instance when granular materials are considered), while they assume a longer range characteristic in the case of continuous materials; this observation suggests to adopt a unified potential-based formulation for the assessment of the interparticle forces, by properly choosing the distance of influence for the transmission of forces.
2.1. Potential-based formulation for particle-particle interaction
Analogously to the atomic description of solids—where the material can be characterized by a potential energy functional n(x)—at the meso- or macroscale it can be assumed the existence of a potential function. The equilibrium equations referred to a single discrete element, corresponding to the configuration of minimum energy of the system, can be written:
d^tot dO tot( x)
P = 0 with
dx; dxi Etot(x) = n(x) = Otot(x)-Pi -x;.,
(1)
where O tot( x) and P; are the strain energy and the force applied to the particle i, respectively.
The stationary condition of the above functional identifies the equilibrium configuration:
d 2 Eo
dx 2
(x - x0):
d 2 O h
(x - xo):
= P -ÔEtot
9x
= K • ( x - xo) = P
(2)
that has been obtained by using the power series expansion of the potential energy function, starting from its equilibrium state identified by the position vector x0:
dEti
Etot( x) = ; Etot(xo) + -
tot
dx
(x -xo) +
+1( x-
2
d 2 Et
tot
ax2
(x - x 0) + ...>- Pi • x ;
(3)
and considering that 9Etot / 9x | Xo = 0 at the equilibrium. In the above relation the stiffness matrix of the system can be recognised; its components are
K.
d 2 Eu
d 2 O.
J 9xl 9x. 9xl 9x.
1 j 1 j
It is worth mention as at the atomic scale, various potential functions have been introduced to represent the mechanical behaviour of materials such as the Morse potential [37] and the Lennard-Jones (LJ) one [38].
2.2. Potential function for continuous and granular materials
The simplest mechanical stress-strain law describing the behaviour of a continuum linear elastic solid, can be represented by the generalized Hooke's relation that assumes the force between two closed material points to be proportional, through proper coefficients, to the strain value occurring in such a small region (Fig. 1).
By considering two particles i, j as representative of the corresponding material points in a continuum solid (Fig. 2, a), their reciprocal stretch—evaluated with respect to the equilibrium distance re —can be written as
£ = -
1
- + -2
1
- + — 2
(4)
where r =| x; - xj |, re = | x0; - x0j | are the distances between the particle centers in a generic case and at the equilibrium state, respectively, while s, se are the corresponding effective distances counterparts (Fig. 2) defined as
Fig. 1. Scheme of pair-wise interparticle forces
Fig. 2. Couple of particles at their equilibrium distance (a). Maximum copenetration corresponding to F (r = r0) ^ (b) and generic relative position of two particles with geometrical parameters (c)
s = r + 28-d + dj )/2 = r + 28-2r. In the above relation the maximum copenetration depth 81 = 82 = 8 = ar has been expressed as a fraction a of the average radius r of the particles. The above relation (4) has been written by considering large displacements effect, while the maximum copenetration distance enables the no-penetration condition between particles to be fulfilled.
By adopting a linear stress-strain relationship to calculate the force acting along the line joining the particle centers, the following expression holds:
Fj (s) =
AjE '
1
if s < rinfl +
0,if s > rinfl +
d{ + d;
- + 28,
(5)
+ d;
- + 28,
where Aj, E' = c(s)E0 are the cross-section and the elastic modulus of a linear element that represents the material connecting the two particles i, j. In the above relation the Young modulus has been assumed to vary with s since the no-penetration condition is represented by an increasing value of E when s ^ 0, i.e. c(s) = 1, if s >> 0, c(s) ^ if s ^ 0. A possible choice of the function c(s) having the above properties can be as follows: c (s) = [1 + (s +1)- m ] x x[1 - (s +1)-m ]-1 in which the exponent m controls the gradient of the elastic modulus versus s function when s ^ 0. Of course the no-penetration condition assumption (s > 0 or 0 < w < 28) entails an unlimited compressive strength of the particles when they approach each other more and more (up to the limit s = 0 or w ^ 28 or r = r0); the particles overlapping distance has been indicated with w, while the minimum particles distance is r0.
The above force-relative displacement relation (5) can be obtained from the functional:
0( s) = AijE 'x
2
s
-sse
A+ 1( s - se)3
(6)
while the corresponding stiffness is given by - ^ j A;E{
K ( s ) =
ds2
1 s-
+
(7)
In Eq. (5) the influence radius rinfl indicates the maximum interacting distance between a couple of particles; it represents the maximum distance beyond which two particles are not interacting, i.e. no bonding force exists. This assumption can be introduced for continuum matters since the internal forces between the material points exist even if they are not in reciprocal contact; this implies the existence of a nonlocal interaction enabling to represent long distance forces suitable to better describe the continuum nature of the material.
In Fig. 3 the distance dependence of the potential, of the force and of the stiffness is displayed for the linear elastic particles interaction by assuming the representative equilibrium distance se/ di = 1, while the same functions are reported in Fig. 3, b for the case of no-penetration constraint.
The consistent representation of the mechanical behaviour of the material connecting the particles, requires a proper and correct definition of the area Aj indicated in Eq. (5); since the influence distance rinfl determines the number of neighbour particles interacting each other, such an area must depend on this maximum distance, i.e. Aj = = a (rinfl) Ao ij, where the reference area can be evaluated as Ao, j = (di + dj )2 nl 4, that is to say such a reference area can be assumed as the cross-section of the cylindrical region of material placed between two particles. The above function a (rinfl) can be evaluated by a best fitting procedure once the particles arrangement in the space and their influence distance rinfl have been assumed.
In order to use the same potential formulation also for granular materials—where the internal forces appear only when the particles are in reciprocal contact and their motion make them to approach each other—it can be simply adopted an influence distance such that rinfl/( dj 2) <= 2, i.e. the interaction of particles occurs only when they are in contact under compression condition, while no forces exist when s' > 0; a tangential force can eventually assumed to be present between the colliding particles i, j when s' < 0. It
Dimensionless effective distance Dimensionless effective distance
between particle surface, s/dj between particle surface, s/dt
Fig. 3. Interparticle potential and corresponding forces and stiffness according to the linear elastic behaviour for an equilibrium distance se/ d = 1 (a); corresponding curves by considering the no-penetration condition (b)
can be evaluated as:
Tt,j(w) = -sign(ut?rei)mm(n 11 u^ ||, Mmd 1 Fj(w)|), (8) where n is the viscosity coefficient and ^md is the dynamic friction coefficient of the materials, while ut rel is the relative tangential velocity between the two particle surfaces in contact. In the above relation the normal force has been indicated with Fj (w) (instead of Fj (s)) to indicate that the friction occurs when the particles are in contact, i.e. an overlapping length arises. Equation (8) provides the tangential force as the minimum value between the dynamic friction and the viscosity action.
It has been determined as the function a(rinfl) is almost independent by the particle size once a specific particles arrangement is adopted, while the Poisson coefficient tends to the value 1/3 as the influence radius increases, irrespectively of the particles layout used to represent the domain of the problem of interest [39].
When the simulation of compact materials is considered, a reasonably choice of the influence radius value could be equal to about rinfl/(d l /2) = 4 (in this way the forces have a limited nonlocal character and exist for neighbour particles even if they are not in direct contact), corresponding to the correction function value equal to about a (rinfl) = = 0.4 for both cubic or tetrahedral arrangements. On the other hand, granular materials can be simulated by adopting rinfl/( d/2) = 2, i.e. the interaction of particles exists only when they are in reciprocal compressive contact, while no forces are present when s' > 0.
3. Smoothed particle hydrodynamics interpretation of the force potential approach
Smoothed particle hydrodynamics approximation of a function defined in the region Q at a point identified by its position vector x l, can be written as [9-14]:
f(x) = J f(x)8(|x - xi|)dQ, (9)
Q
i.e. the sought value of the function at x l is obtained through an interpolation of the known values f (x l) over the domain Q by mean of a convolution integral.
For numerical purpose the Dirac delta function S(| x - x 11) in Eq. (9) is usually approximated by adopting a proper decreasing function W(|x - x 11) of the distance |x - x (often called kernel or smoothing function, [13, 14]). The above approximation becomes:
f (xi) =| f (x )W (| x - x |)dQ,
Q
with W(| x - xl |) such that: 7 W (|x - xi |,h)dQ = 1,
9Q
-W(|x-xl |,h) = 0 if |x-xl |> h, (10)
lim W (| x - xl |, h) = S(| x - xl |)
h^0
in which the smoothing length (or support dimension h) defines a spherical region centred at x l, | x - x 11 < h, where the volume integral has to be evaluated.
Usually a Gaussian-like, W(| x - x 11, h) = e"|x-x'2/h x x(rch2)"n2, 1 < n < 3 or B-spline functions,
W (| x - x 11, h) = (1 - 3/2 x2 + 3/4 x3)/(fth3), 0 < x = | x - x l H h < 1,
W(| x - xi |, h) = [1/4 (2 - x)3 ]/(nh3), 1 < x < 2, are used for the kernel expression, where the parameter h is the smoothing length.
If the nth derivative of the above function is sought, it can be assessed through the smoothed particle hydrodynamics approximation by using the following expression:
f(n)(xi) = | f (x )W(" )(| x - x 11 )dQ. (11)
Q
The above relations (10), (11)—when evaluated numerically—can be rewritten by replacing the integrals with summations:
f (Xi) = £m,^^W(|x-Xi|, h),
k (12)
f(n)(Xi) = E W(n)(| x, - Xi|, h)
k Pk
that have to be made over the index k spanning all the particles; xk, mk, pk are the position, mass and density of the generic particle k, respectively, and f (x k) is the function value at x,. In practice the summation is made only over the neighbours particles of the particle placed at xi because of the compact support of the kernel function W(|x, -x-|, h).
As is well-known the dynamic equilibrium equations of an infinitesimal material volume can be written as:
do,-,- dv■
—j-p-i- + bi = 0 or o, , -px + bi = 0, dx; dt ' 11 -
(13)
where tfj, p, vt, bi are the stress tensor components, the mass density, the velocity and the body force in the i-th direction, respectively.
Adopting the above numerical smoothed particle hydrodynamics representation for the stress tensor,
(x p) = L mk
oij (xk) pk
W(|Xk -xI,h),
oij, j(x p) = L mk
(xk)
(14)
W - (| Xk - x p h)
- k Pk
the equilibrium equations (12) in the three dimensional space, written at the position of the particle p, can be expressed as:
a, (xk)
L mk—-W j(| xk - xp K h) - Ppxi(p) + bi = 0
k Pk i = 1, 2, 3, or
Vp L Vk Ojj (xk ) W, j (|xk - x p|,h )
k
-VpPp x(p) + Vpbi = 0.
(15)
The second equation of (15) has been obtained by multiplying each term of the first equation of (15) by Vp (that represents the volume of the particle p), and by posing
Pi = mk/Vk •
It can be recognised as the first term in the first equation of (15) represents the internal force vector, while the second term corresponds to the inertial term of Eq. (13) and the last one provides the body force vector components.
On the other hand, by considering the dynamic equilibrium equation written for the particle p:
Fi(p) - mp xp + Ftot(p) = 0 (16)
and by substituting the expression (5) for the internal force exerted on the particle p by the generic particle k, we get:
Vp L
a(h) AmkE'
k[ p]
Vp
-s,
e(kp)
se(kp )
1
+
2
(
xk xp \ se(kp)
se(kp)
mp xi(p) + ^tot, i(p)
J = 0,
qi(kp) '
(17)
where k[p] indicates the generic particle k falling in the neighbour of particlep, i.e. all the particles included in the domain of influence of particle p, according to the smoothing function adopted; se = | xk - x | e is the equilibrium distance between the particles k andp and q is the unit vector identifying the direction p-k. It can be recognised as: Vp E VkOy(Xk)Wy(|Xk -xp|, h) =
k[ p]
= Vp L Vk
k[ p]
(l. -
a(h) AmE '
VkVp
\xk - x ,
se(kp)
se(kp )
1
+
2
-s,
e(kp)
se(kp)
qi(kp ),
(18)
i.e. the particle equation written by expressing the forces through the potential function, provides an expression analogous to the smoothed particle hydrodynamics approximation, where qi(jp) is the i-th component of the versor connecting the particles k and p. In other words, it is instructive to note that the term a j (xk) Wj (| xk - x p |, h) in the smoothed particle hydrodynamics approximation corresponds to the force existing between the particles p-k divided by the product of their volumes, Vk Vp.
3.1. Strain and stress evaluation
According to the model for the particle reciprocal forces presented above, an uniaxial stress state acting in the direction of the line joining two interacting elements can be assumed to arise in the material. Since there are several particles within the influence distance of a specific particle i and given that the stress state at xt must be those of the medium representing the domain of the problem, it is straightforward to think the strain measured in the generic direction i-j to be expressed through the classical transformation law:
Ey (xi- j) = £(xi) nj,
eï (xi- j) =
\ xi x j \ se(ij)
se(ij)
(19)
1
+
2
(
xi x j \ se(ij)
se(ij
where ny is the unit vector identifying the i-j direction and Ey (xi-y ) is the strain along such direction itself. Having the strains Ey (xi), j = 1, 2, ..., p inp different directions
(where p indicates the number of particles that are within the region of influence of the particle of interest i), the above relations (19) can be used to determine the strain tensor components at xi, for example through a least square-like algorithm or others equivalent approaches.
The above procedure, in order to provide acceptable strain values, requires that p > 6 and that the particles surrounding the particle i are well distributed in the space.
Once the strain field in the matter is known, the corresponding stress field can be evaluated by using the proper stress-strain constitutive relation of the material; in incremental form it can be written:
do(x,.) = C'd £(x,), (20)
where C' is the tangent elastic tensor of the material and d£(Xj) is the strain rate tensor at xi.
The evaluation of the stress field within all the domain of the problem must be carried out by determining the stress tensor at all the particles point, ^(Xj). The stress tensor— known at discrete points—can be finally interpolated through a smoothing interpolation procedure such that expressed by Eq. (10) where the smoothing length can be assumed identical to the influence distance adopted for the particles interaction.
The force arising between two particles of the medium can be easily exploited to determine the failure condition between them; for brittle continuum materials it can be simply verified by checking if the maximum principal stress exceeds the tensile strength, while under compression the crushing failure can be represented through an abrupt reduction of the influence distance, i.e. by reducing the material to an assembly of incoherent grains that behave only when compressed each other, with (eventually) the limitation of the no-penetration condition.
3.2. Particle-boundary contact
The particle-boundary interaction can be modelled similarly to the particle-particle case presented above; in particular, the contact mechanics concepts can be adopted to particularise the force-particle overlapping law. In the case of an impact between an elastic spherical particle and an elastic plane the contact stiffness, the equivalent Young modulus, and the equivalent contact radius can be defined as [34]:
K (w) = ^TT = 2^/fwV2' (21)
2 1 ,,2 A-1 ,»
1 -v2 1 -v- * d di
E = -- +-- , r =— = —,
Ei E. 2 2
where Ei, vi, dt are the elastic modulus, the Poisson ratio and the diameter of the particle, while E-, v- are the elastic modulus and the Poisson ratio of the elastic boundary surface (assumed to be flat, i.e. with a local radius of curvature of the contact area tending to infinity). The normal force,
acting perpendicularly to the tangent plane at the particle-surface contact point, and the tangential one can be evaluated as:
T (w) = Kw3/2n,
X , (22) T =-MT„| = -(^K„w3/2)t, 1 xi,t1
where n is the unit normal to the boundary surface in the contact point and -t = - Xitj | xi 11 represents the unit vector opposite to the direction of the tangential velocity with respect to the boundary of the particle i.
In order to avoid possible numerical instabilities due to the high value of the contact force given by Eqs. (26), (27) when a very stiff plane is considered, the contact force co-penetration relation, Tn (w), can be attenuated by introducing a smoothing function x(w) for the contact stiffness:
Kn (w) = X( w) Kn , (23)
0 < x(w) = (w/y)1 m <1 , m >1.
The above relationship corresponds to assume that the interaction with the boundary takes place when the particle touches a layer having a small thickness y superposed to the surface. The contact force increases progressively up to the final copenetration depth w.
It can be observed that the potential function concept developed for the particles interaction is suitable also to describe the particle-boundary forces. In fact it can be expressed through a suitable potential with an influence distance (measured between the particle surface and the surface boundary) for the impacting particle i equal to rinfl = = ri +y (see Eq. (5)).
Moreover the case of continuum solid or granular materials colliding with others elastic solids or boundaries can be treated exactly in the same way by simply using the proper force potential.
In the examples presented below a soft layer with relative thickness equal to y/(di/ 2) = 0.1 and m = 4 (see Eq. (23)) is adopted for the numerical simulations.
4. Numerical implementation
The dynamic governing equations of the discretized problem are given by:
-MX + Fi + Fd + Fe + Fb = 0 or MX = Ftot, (24)
where M is the mass matrix, X is the vector of the particle center accelerations (note that the rotation degrees of freedom are neglected since a small rotational inertia of the particles is assumed), while Fi, Fd, Fe, F6 are the internal force vector (obtained by assembling the forces F.- (i)q, see Sect. 2.2), the damping force vector, the external force vector and the vector of actions due to the collision of the particle with the elastic boundaries, respectively.
For numerical stabilization purpose, the above dynamic equations can written by adding some damping forces. Note
that in the above described model no velocity-related damping has been adopted for the particles and so it is convenient to introduce an artificial numerical damping. It can be thought to be opposite to the particle current velocity, i.e., according to Cundall and Strack [24], the following equation for the numerical damping force Fd can be established:
-XdFsign (Fx),
(25)
where Xd is the damping coefficient, sign(-) is the sign function and F, X are the actual force on the particle of interest and its velocity, respectively.
An explicit numerical integration scheme can be adopted for solution of the equations of motion, allowing a separate evaluation for each particle once the force vector acting on it is known. The solution of a large equation system is therefore avoided. Explicit time integration methods, also if conditionally stable, are often adopted in particles simulations due to the reduction of memory requirement thanks to the absence of matrices such as the stiffness one.
In the present particle method, the explicit leapfrog integration scheme—or Verlet integration approach—is adopted [40]; it provides a satisfactory numerical stability, time-reversibility and preservation of the symplectic form on phase space. Its name recall that even time derivatives of position are known at on-step points, whereas odd derivatives are available at mid-step points. The method operates by determining the particle acceleration at the current time step n and consequently the velocity at the time step n +1/2 as:
xi,n = FTinlmi, xi,n+12 = xi,n-1/2 + At Xi,n > (26)
where Xi n-1/2 = (Xi n - Xi,n-1)/At is the mean velocity vector across the time step n + n +1. The position of particle i at the time step n + 1 is given by the vector:
Xi, n+1 _
Xi, n + At ( xi, n-1/2 +A¿ xi, n ) or
vi ,n+1
_ Xi,n + i, n+1/2 •
(27)
It can be noted that the positions of particles are obtained at time instants tn = nAt, while the velocities are known at tn+1j2 = (n ± 1/2)At.
The upper limit of the time integration step At is a crucial parameter to ensure stable results. This value is generally obtained from the smallest natural oscillating period of a couple of particles, At = a^J m/K, where the minimum ratio mlK among all the particles constituting the dis-cretized solid must be considered and a < 1 is a proper constant to be introduced [41].
5. Numerical applications
5.1. Elastic cantilever beam under harmonic load
The capability of the present discrete approach to model continuum elastic bodies is verified in the present section by studying the dynamic response of a cantilever beam under a harmonic shear force applied at its free extremity (Fig. 4).
The beam is supposed to have an elastic modulus equal to E = 5 • 109 Pa and mass density equal to 2400 kg • m3, while its dimensions are as follows: length l = 0.6 m, height b = 0.06 m and width a = 0.03 m; the harmonic concentrated load F is assumed to vary according to the relation F (t) = F0sin(2nft), where F0 = 100 N and the frequency is f = 7.77 Hz; a negligible damping has been considered.
The force potential described in Sect. 2.2 has been used to describe the interaction forces inside the beam by adopting an influence distance of the particles equal to double of the particle diameter, rinfl = 2 d, while the beam is discre-tized through about 1700 particles arranged in cubic fashion (Fig. 4, a).
2D plane stress finite element analysis of the beam (discretized by 216 quadrilateral quadratic elements) has been performed for comparison of the results with those provided by the particle method. The time history of the vertical displacement and velocity of point A, located at the free edge of the beam, is represented for both the finite element and the discrete model (Fig. 4, b).
It can be observed that the discrete approach results are in satisfactory agreement with the finite element ones, both in term of magnitude of displacement and velocity, as well as of period of oscillation, thus confirming the good capability of the developed particle methods approach to correctly represent the dynamic behaviour of an elastic body.
5.2. Impact of a cubic block on a cantilever beam
This example considers the impact of a cubic elastic body on a cantilever beam. The falling cube has an initial velocity equal to v0 and is located above the beam at a distance H, its horizontal position with respect to the beam
Clamped section
2 ^ l
AV(t)
0.004
0.3 b
-0.004
-0.3
0.00
0.04 0.08 Time t. s
0.12
Fig. 4. Cantilever beam under harmonic load (a); time histories of the free edge vertical displacement and velocity (b). FE—finite element method, PM—particle method
is given by the distance d, while it is rotated by an angle a with respect to the horizontal plane. The material of the beam is supposed to be brittle with an elastic modulus equal to E = 3-106 Pa and a tensile strength equal to ft = 2-105Pa, while the force potential described in Sect. 2.2 has been used to describe the interaction forces inside the solids and between the different parts of the beam when failure takes place, by adopting an influence distance of the particles equal to rinfl = 2d and rinfl = d, respectively.
Below the beam an elastic flat horizontal boundary surface n lying on the x-y plane (with elastic modulus E = = 3-1010 Pa) is assumed to be placed at a distance h from the bottom of the beam (Fig. 5). In the numerical particle approach the beam, assumed to have an elastic modulus E = 3 -106 Pa and a tensile strength of 2 -105 Pa, is modelled through 1200 spheres (arranged in a cubic lattice with diameters assumed to be normally distributed with a mean value equal to 5 mm and variance 0.1 mm), while the falling body is supposed to have an elastic modulus equal to E = = 3 -107 Pa and a very high tensile strength to prevent its failure. The geometry of the system is characterised by a = = 0.02 m, L = 10a, H = a and h = 5a/2, while the falling body (with an initial velocity equal to v0 = 3 m/s) is assumed to be characterised by: (i) d = 0, a = 0°; (ii) d = a, a = n/ 4.
It must be noticed as the elastic constants of the bodies have been assumed intentionally very small to allow the development of large displacements in the structure in order to emphasise the geometrical nonlinear feature of the computational algorithm and to limit the time integration step for numerical reasons. The mass density of the beam has been assumed equal to 2400 kg -m-3, while the falling block density has been assumed equal to 7800 kg -m-3. The time integration procedure has been conducted by using a time step amplitude At = 5 ^s and the analysis duration has been extended up to 0.4.
In Fig. 6 the configurations for the falling case (i) and (ii) at three time instants are shown. It can be noted as the falling block causes a local rupture of the beam in the im-
Fig. 5. Cubic elastic block falling on a cantilever beam: geometry of the system
pacted zone and the failure of the restrained cross section of the beam for case (i), while a local failure and a localised failure placed approximately at the beam midspan takes place for the case (ii).
The rest of the falling block and of the failed portion of the beam on the boundary elastic surface placed below the system can be appreciated.
5.3. Granular column drop test
In this example, a granular column with an initial height-radius ratio equal to approximately 1.51 is considered (Fig. 7). The column is prepared by dropping particles into a box and letting them to settle under gravity [42]. Particles— initially arranged in a tetrahedral fashion—are assumed with diameter values obeying normal distribution with a mean value equal to 18 mm and variance of 5 mm, for a total number of particles equal to 1440. The mechanical properties of the material are assumed to be as follows: mass density 2000 kg - m-3, elastic modulus E = 3 -106 Pa and a negligible tensile strength (or equivalently no cohesion). A very low elastic modulus value has been assumed to reduce the computational cost since in this way a larger time step can be used without affecting the stability property of the integration scheme; moreover such low assumed values, has been verified that affect only slightly the configurations assumed by the particles during their settling motion. The friction
Fig. 6. Beam configuration at different time steps for: case (i) (a-c), case (ii) (d—f). Gray scale indicates the velocity values (in the vertical direction, in ms—1). t = 0.016 (a, d ), 0.024 (b, e), 0.200 s (c, f )
Fig. 7. Configuration at initial and final time instants obtained with the present study (a) and according to [42] (b). In the left figures the gray scale indicates the vertical displacements of the particles in m
coefficient of the base supporting the column is = 0.5 while smooth vertical walls, are placed beside the column.
The drop test is conducted by removing the left hand side wall of the box and letting the column to spread under gravity.
It is possible to note a good agreement between the literature and the results provided by the present particle model. The final configuration of the particle assembly provided in [42] is characterized by the height/width ratio equal
to 0.2185, while the present model gives the value 0.2230 (Fig. 7).
5.4. Granular material cutting simulations
This example considers a vertical flat blade placed between two glass panels. The blade has a constant horizontal velocity equal to 10 mm - s-1 [43] and drives particles that closely resemble natural granular flow into dragline buckets. The grain matter is suitable for a discrete simulations
Fig. 8. Configuration of the particles for a blade displacement equal to 10 (a, c), 20 cm (b, d ); configurations given in [43] (a, b) and by the present model (c, d ). Gray scale indicates the horizontal displacement of the particles expressed in m
because of their nature. An initial volume of particles (Fig. 8) under the gravity action is assumed to be placed inside a box, while the vertical blade is placed at the left hand side position. The simulation consists in moving the blade, causing the material inside the box to redistribute. The performed dynamic analysis is aimed at determining the particles configuration at different time instants.
The material of the particles (assumed with a mass density equal to 855 kg-m-3) is supposed to have an elastic modulus equal to E = 2.76-106 Pa and a negligible tensile strength (or equivalently no cohesion). The force potential has been used to describe the particles interaction by adopting an influence distance of the particles equal to rinfl = d, while the box and the blade are supposed to have an elastic modulus equal to E = 3 -108 Pa.
The particles volume is modelled through about 1162 spheres, initially arranged in a cubic fashion, with diameter size equal to 20 mm. The geometry of the system is depicted in Fig. 8.
The integration time step has been taken equal to At = = 0.2 ms and the final blade displacement is equal to 200 mm, with an analysis duration equal to 2.0 s. The damping coefficient for the dynamic analysis has been assumed equal to Ad = 0.2 and the coefficients of dynamic friction between the particles and for particle-boundary planes are supposed to be = 0.1.
The obtained arrangement of the particles shows a behaviour similar to soil cutting problem analyzed in [43]. It is possible to note that in the present study the blade is shorter with respect to the literature test, causing a little difference in the disposition of the particles near the top of the blade. Neglecting this last aspect, the shape assumed by the particles representing the granular material obtained by the present study, reflects the configuration presented by Coet-zee in [43].
6. Conclusions
A discrete approach for the simulation of materials, suitable for the description of solids at different scales, has been presented as a valid alternative to the classical continuous approaches. In particle method approaches, complex nonlinear interactions between elements and within elements are numerically simulated to predict the motion of discrete constituents by the solution of dynamic nonlinear differential equations.
In this study a force potential-based discrete element method for continuum or granular-like materials has been developed in the framework of the so-called smoothed particle hydrodynamics methods. Exploiting the possibility to assume a discrete nature of materials even at large scales, the potentialities and generality of the method have been underlined.
The proposed particle method has been demonstrated to be usable for the simulation of continuum as well as for
incoherent aggregate matters, and in this sense represents a unified approach to simulate the dynamic of mechanical systems.
The developed particle method owns a mesh-free nature and is characterized by a nonlocal interactions between elements, allowing to simulate very different mechanical problems. This feature enables to easily account for failure of compact and granular materials by taking into account for the severe geometrical distortion occurring in such problems as well as reciprocal particles interaction and particle-boundary interaction.
Some numerical examples involving the dynamic response of an elastic body, the failure of brittle solids under impact actions and the flow of granular materials have been presented, underlying the wide capability of the developed discrete method.
Acknowledgments
The authors gratefully acknowledge the research support for this work provided by the Italian Ministry for University and Technological and Scientific Research (MIUR).
References
1. Clough R.W. The Finite Element Method in Structural Mechanics // Stress Analysis / Ed. by O.C. Zienkiewicz, G.S. Holister. - London: Wiley, 1965.- P. 85-119.
2. Oden J.T. An Introduction to the Mathematical Theory of Finite Elements. - New York: Wiley, 1976.
3. Dow J.O., Jones M.S., Harwood S.A. A generalized finite difference method for solid mechanics // Num. Meth. Partial Diff. Eqs. - 1990. -V. 6. - No. 2. - P. 137-152.
4. Wrobel L.C., Aliabadi M.H. The Boundary Element Method. -Chichester: Wiley, 2002. - 1066 p.
5. Belytschko T., Krongauz Y., Organ D., FlemingM., Krysl P. Meshless methods: An overview and recent developments // Comp. Meth. Appl. Mech. Engng. - 1996. - V. 139. - No. 1-4. - P. 3-47.
6. Sukumar N., Moran B., Belytschko T. The natural element method in solid mechanics //Int. J. Num. Meth. Engng. - 1998. - V. 43. - No. 5. -P. 839-887.
7. Lucy L.B. A numerical approach to the testing of the fission hypothesis // Astron. J. - 1977. - V. 82. - P. 1013-1024.
8. Gingold R.A., Monaghan J.J. Smoothed particle hydrodynamics: Theory and application to non-spherical stars // Mon. Not. R. Astron. Soc. - 1977. - V. 181. - P. 375-389.
9. Benz W. Smooth Particle Hydrodynamics: A Review // Numerical Mo-
deling of Non-Linear Stellar Pulsation, Problems and Prospects / Ed. by J.R. Buchler. - Boston: Kluwer Academic, 1990.
10. Monaghan J.J. Why particle methods work. SIAM // J. Sci. Stat. Comput. - 1982. - V. 3. - No. 4. - P. 422-433.
11. Monaghan J.J. Smoothed particle hydrodynamics // Annu. Rev. Astron. Astrophys. - 1992. - V. 30. - P. 543-574.
12. Monaghan J.J. SPH elastic dynamics // Comp. Meth. Appl. Mech. Engng. - 2001. - V. 190. - P. 6641-6662.
13. Monaghan J.J. Smoothed particle hydrodynamics // Rep. Prog. Phys. -2005. - V. 68. - P. 1703-1759.
14. Monaghan J.J. Smoothed particle hydrodynamics and its diverse applications // Annu. Rev. Fluid Mech. - 2012. - V. 44. - P. 323-346.
15. Li S., Liu W.K. Meshfree and particle methods and their applications // Appl. Mech. Rev. - 2002. - V. 55. - No. 1. - P. 1-34.
16. Bui H.H., Fukagawa R., Sako K., Ohno S. Lagrangian meshfree particles method (SPH) for large deformation and failure flows of geo-material using elastic-plastic soil constitutive model // Int. J. Num. An. Meth. Geomech. - 2008. - V. 32. - No. 12. - P. 1537-1570.
17. Harthong B., Jerier J.-F., Doremus P., Imbault D., Donze F.-V. Modeling of high-density compaction of granular materials by the discrete element method // Int. J. Solids Struct. - 2009. - V. 46. -P. 3357-3364.
18. Pana W., Tartakovskya A.M., Monaghan J.J. Smoothed particle hydrodynamics non-Newtonian model for ice-sheet and ice-shelf dynamics // J. Comput. Phys. - 2013. - V. 242. - P. 828-842.
19. Bessa M.A., Foster J.T., Belytschko T., Liu W.K. A meshfree unification: reproducing kernel peridynamics // Comput. Mech. - 2014. -V. 53. - No. 6. - P. 1251-1264.
20. Braun M., Fernandez-Saez J. A new 2D discrete model applied to dynamic crack propagation in brittle materials // Int. J. Solids Struct. -2014. - V. 51. - P. 3787-3797.
21. Curtin W.A., Miller R.E. Atomistic/continuum coupling in computational materials science // Model. Simul. Mater. Sci. Eng. - 2003. -V. 11. - P. R33-R68.
22. Liu B., Huang Y., Jiang H., Qu S., Hwang K.C. The atomic-scale finite element method // Comp. Meth. Appl. Mech. Engng. - 2004. -V. 193. - P. 1849-1864.
23. Hoover W.G. Computational physics with particles—nonequilibrium molecular dynamics and smooth particle applied mechanics // Comput. Meth. Sci. Tech. - 2007. - V. 13. - P. 83-93.
24. Cundall P.A., Strack O.D.L. A discrete numerical model for granular assemblies // Geotechnique. - 1979. - V. 29. - P. 47-65.
25. Onate E., Idelson S.R., Del Pin F., Aubry R. The particle finite element method. An overview // Int. J. Comput. Meth. - 2004. - V. 1. - No. 2. -P. 267-307.
26. Tavarez F.A., Plesha M.E. Discrete element method for modelling solid and particulate materials // Int. J. Num. Meth. Engng. - 2007. -V. 70. - P. 379-404.
27. Particle-Based Methods: Fundamentals and Applications / Ed. by E. Onate, R. Owen. - Springer, 2011.
28. De Gennes P.G. Granular matter: a tentative view // Rev. Modern Phys. - 1999. - V. 71. - P. S374-S382.
29. D'Addetta G.A., Kun F., Ramm E. On the application of a discrete model to the fracture process of cohesive granular materials // Granular Matter. - 2002. - V. 4. - P. 77-90.
30. Rycroft C.H., Kamrin K., BazantM.Z. Assessing continuum postulates in simulations of granular flow // J. Mech. Phys. Sol. - 2009. - V.57.-P. 828-839.
31. Guzev M.A., Dmitriev A.A. Bifurcational behavior of potential energy in a particle system // Phys. Mesomech. - 2013. - V. 16. - No. 4. -P. 287-293.
32. Krivtsov A. Molecular dynamics simulation of impact fracture in polycrystalline materials // Meccanica. - 2003. - V. 38. - P. 61-70.
33. Kuzkin V.A., Krivtsov A.M., Jones R.E., Zimmerman J.A. Material frame representation of equivalent stress tensor for discrete solids // Phys. Mesomech. - 2015. - V. 18. - No. 1. - P. 13-23.
34. Obermayr M., Dressler K., Vrettos C., Eberhard P. A bonded-particle model for cemented sand // Comp. Geotechnics. - 2013. - V. 49. -P. 299-313.
35. Brilliantov N., Spahn F., Hertzsch J., Poshel T. Model for collision in granular gases // Phys. Rev. E. - 1996. - V. 53. - P. 5382-5392.
36. Aubry R., Idelsohn S.R., Onate E. Particle finite element method in fluid mechanics including thermal convection-diffusion // Comput. Struct. - 2005. - V. 83. - P. 1459-1475.
37. Morse P.M. Diatomic molecules according to the wave mechanics. II. Vibrational levels // Phys. Rev. - 1930. - V. 34. - P. 57-64.
38. Liu W.K., Jun S., Qian D. Computational Nanomechanics of Materials // Theoretical and Computational Nanotechnology / Ed. by M. Rieth, W. Schommers. - Stevenson Ranch, CA: American Scientific Publishers, 2005.
39. Brighenti R., Corbari N. Dynamic failure in brittle solids and granular matters: a force potential-based particle method // J. Numer. Meth. Engng. - 2015. - doi 10.1002/nme.4998.
40. Verlet L. Computer experiments on classical fluids. I. Thermodyna-mical properties of Lennard-Jones molecules // Phys. Rev. - 1967. -V. 159. - P. 98-103.
41. O'Sullivan C., Bray J.D. Selecting a suitable time step for discrete element simulations that use the central difference time integration scheme // Engng. Comput. - 2004. - V. 21. - P. 278-303.
42. Lim K-W., KrabbenhoftK., Andrade J.E. A contact dynamics approach to the granular element method // Comp. Meth. Appl. Mech. Engng. -2014. - V. 268. - P. 557-573.
43. Coetzee C.J. Discrete and continuum modelling of soil cutting // Comp. Part. Mech. - 2014. - V. 1. - No. 4. - P. 409-423.
nocTynaaa b peaaKUHro 03.07.2015 r.
CeedeHun 06 aemopax
Roberto Brighenti, PhD, Assoc. Prof., University of Parma, Italy, brigh@unipr.it, roberto.brighenti@unipr.it Nicholas Corbari, PhD, University of Parma, Italy, nicholas.corbari@gmail.com
