Lattice-based six-spring discrete element model for discretisation problems of 2D isotropic and anisotropic solids/Sesiu spyruokliu tinklelio diskreciuju elementu modelis dvimatems izotropinems ir anizotropinems terpems diskretizuoti.

Kacianauskas, R. ; Vadluga, V.

1. Introduction

Computer simulations of various mechanical systems and processes have a long history, while computational mechanics may be referred to the most successful subject, taking into account a great number of applications. It should be noted that the development of any numerical method involves several issues. A description of material distribution in space or the so-called discretisation concept is of major importance here which significantly increases the scope of application of these methods to nonlinear dynamics.

Historically, the finite element method (FEM) prevailed in simulation of solids and structures. Originally, FEM technique which was based on fixed meshes aimed to solve static problems. Probably, the most general FE approach to non-linear dynamic solid problem is generalized in the book of Belytschko et al [1]. This general approach describes both the motion of the material and the motion of the mesh. It requires the evaluation of a new mesh in each time step and is quite expensive when applied to heterogeneous solids with changing topology.

Recently, various discretisation alternatives combined with explicit time integration methods have been suggested. A concept of the discrete, originally distinct, element method is basically referred to the original work of Cundall and Strack [2]. It was aimed at describing mechanical behaviour of granular assemblies composed of discrete elements, i.e. discs, and (later) spheres, in particular. Some examples of discrete element method (DEM) developments and algorithmic details may be found in [3-7].

DEM concept has been also extensively applied to the simulation of heterogeneous solids to study their dynamic deformation behaviour and fracture problems. These methods could be characterised as being at the development stage, theretofore, many issues including the discretisation concept should be clarified.

Two fundamentally different, particle-shaped and lattice-based, approaches have been developed in the framework of DEM applied to simulate solids. Particle-shaped approach presents a rather straightforward extension of the original Discrete Element Method. In this paper, a solid is replaced by a composition of discrete particles. Cohesive forces and various link and detachment mechanisms acting between the particles may be applied to solving this problem [8-10].

A lattice-based approach representing the continuum by material particles interacting via the network elements is considered in this paper. Generally, this approach comprises an atomic concept used basically by physicists and straightforward application of the structural concept used by engineers. The idea to approximate the continua dates back to the works of Hrennikoff [11] and later contributions of Kawai [12] and Herman et al. [13].

The main problem in developing lattice models is associated with proper evaluation of elasticity parameters. The earliest models of 2D isotropic continuum were described by equilateral triangles in terms of the normal and shear spring stiffness. Their explicit expressions were probably given first by Savamoto et al. [14], and then theoretically derived by Griffiths and Mustoe [15]. The above two-spring model has a narrow range of Poisson's ratio restricted by the upper limit value of 0.33 for a plain stress problem. Later, the same approach was applied to 3D hexagonal assembly [16], cubic lattice [17], and anisotropic solids using 2D hexagonal and square lattices [18].

The two-spring (normal and shear) model and numerical problems have been recently considered in [19]. It was shown there that the presence of the shear spring is difficult to interpret from the viewpoint of continuum.

It could be noted that simulations, employing classical Bernoulli-Euler beams [20, 21] or Timoshenko beams [21, 22] may be treated as a counterpart of non-classical continua. The interpretation of rotational degrees of freedom is largely empirical.

A comprehensive review and theoretical study of the planar elastic lattice models which hold for micromechanical applications is given by Ostoja-Starzewski [23]. A comparison of various approaches on the basis of numerical analysis is also presented by Karihaloo et al. [20].

The classical continuum-consistent models involve only translational degrees of freedom, therefore, they ideally should be represented by a network of normal springs (truss elements). It was proved, however, that for isotropic continuum described by periodic triangular lattice with central interaction, a single discrete elasticity parameter is available. This model is valid for Poisson's ratio value of 0.33, while its deficiency for other values is numerically illustrated in [24].

In [23], it was also shown that a description of anisotropic behaviour requires six independent constants, while the introduction of angular parameters is the only way to vary Poisson's ratio.

Several approaches appeared to solve Poisson's ratio problem numerically by employing the micromechanical concept. A model for the elastic continua with virtual multidimensional internal bonds was developed by Zhang and Ge [25]. The basic idea of the presented model is that materials are discretised into mass particles, while these mass particles are connected with randomized normal and shear bonds. The constitutive relationship bridges the virtual bond stiffness and the macromaterial constants, i.e. Young's modulus and Poisson's ratio.

In the framework of a similar concept and by applying megaclustring of DE, a model of homogeneous material to determine the DEM contact normal and tangential stiffness as functions of the known material properties was developed by Tavarez and Plesha [19]. The convergence is obtained with model refinement and the correct elastic behaviour is produced iteratively.

Based on the review made, it could be stated that the existing DEM models for continuum are diverse, while a unified method for the evaluation of discrete elasticity constants consistent with the classical continuum theory is still at the development stage. In particular, singularity of Poisson's ratio has not been resolved explicitly yet.

This paper continues the analysis of discretisation of elasticity properties considered in the previous work [24]. It aims to develop a six-spring lattice-type DE model applicable to the description of isotropic and anisotropic 2D elastic solid problem and exhibiting diversity of Poisson's ratio values. The original contribution is based on the alternative description of angular interaction defined by applying the natural finite element technique [26].

The paper is arranged as follows. In Sections 2 and 3, DEM methodology and relevant discretisation by periodic triangular lattice are discussed. Derivation of discrete elasticity relations described in the six-spring model by the natural finite element technique is given in Section 4. Validation and performance of the developed method by simulating the dynamic behaviour of the plane stress problem are presented in Section 5. Finally, some conclusions are drawn in Section 6. Evaluation of the natural stiffness matrix in term of dimensionless parameters is described in the annex.

2. DEM methodology

The time-driven DEM is applied to the simulation of dynamic behaviour of the elastic two-dimensional solid. Actually, the present work is restricted to the plane stress problem, but the extension to plane strain or three dimensions would be a rather formal task. Consequently, the plate of constant thickness s is regarded here as two-dimensional solid. It is subjected to in-plane loads attached to the middle plane. The solid is considered in plane Oxy of the Cartesian co-ordinates, while axis Oz points thickness direction.

Generally, elasticity properties of solids are defined by elasticity tensor and may be described in terms of 3 x 3-order symmetric elasticity matrix [D]. For an anisotropic material, the elasticity properties are defined by six independent constants. For an isotropic material, the elasticity properties are restricted by two independent constants, i.e. by the elasticity modulus E and Poisson's ratio v. If the material is assumed to be heterogeneous, the material constants may be defined as position x = [{x, y}.sup.T] dependent variables, and heterogeneity may be imposed on the particular subdomains.

The DEM discretisation approach relies on the concept applied to the description of granular material as an assembly of deformable inter-acting particles [2]. Each of the particles is considered separately and defined by its shape, size and physical properties. DEM operates upon binary contact between two particles, however, the equilibrium and kinematics of an individual particle are defined with respect to the particle center.

The 2D solid is regarded as a system of the finite number N of deformable material particles i (i = 1, ...,N). The motion of each particle i in time t is described by the second Newton's law. In our problem, the particle i has two independent translations denoted hereafter by vector [x.sub.i] = [{[x.sub.i], [y.sub.i]}.sup.T], therefore, the dynamic equilibrium of a solid particle is defined by two equations for forces as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (1)

where [m.sub.i] is mass of the particle, while right hand side parameters stand for the resultant vector [F.sub.t] = [{[F.sub.xi], [F.sub.yi]}.sup.T] of all external and inter-particle forces.

The motion of each particle i, or, more definitely, tracking of its position [x.sub.i] (t), velocity [v.sub.i](t) [equivalent to] [[??].sub.i](t) and acceleration [a.sub.i] (t) [equivalent to] [[??].sub.i](t) in time t is performed by integration of equations (1). The explicit time integration schemes presents the most appropriate technique used in DEM. The Verlet velocity algorithm is currently applied to integration [3].

The DEM approach described above was implemented into the original software code. The code presents a modified version of the DEM code DEMMAT developed in the Laboratory of Numerical Modelling of Vilnius Gediminas Technical University, see [6].

3. Discretisation by periodic triangular lattice

The lattice-type discretisation methodology aimed at developing DEM models for Eq. (1) and evaluating the discrete elasticity parameters is elaborated hereafter in detail.

The discrete model is implemented by covering a computational domain with the hexagonal lattice grid. The lattice is constructed by equilateral triangles (Fig. 1, a). Each particle i represent a hexagon composed of six equal triangles. The hexagon encompasses a half of each connection line. The location of the particle matches that of the lattice node and is defined by the global coordinates [x.sub.i] = [{[x.sub.i], [y.sub.i]}.sup.T], while the geometry of the particle is defined by a characteristic dimension L of the grid. Consequently, two types of meshes appear by modelling a planar continuum domain, where another polygonal tessellation is required for outlining the particle's shape. The grid cell with a particle is shown in Fig. 1, b.

The density of the material is constant within the particle and is assumed to be constant [rho]. Mass of the solid is described by a set of lumped masses [m.sub.i], concentrated in the centres of particles:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (2)

Generally, the lattice concept replaces continuum by a discrete network. The actual inter-particle contact is not considered here, and is simple replaced by the equilibrium of the selected nodes. Consequently, the interaction of particles i and j is described by the connection element i-j.

[FIGURE 1 OMITTED]

The connection element presents a line segment of the lattice grid. The network element is shown in Fig. 1, c. Constitutive properties of the solid are assigned to particular lines, more definitely, to the volume bounding connection line i-j. Here, the particle interface, or the effective section, aimed to transmit interaction forces between the particles, is of particular importance. For the 2D case, the bounding volume (area) is defined by the effective width [B.sup.eff.sub.ij]. The connection element may generally be regarded as a rheological element which reflects a highly complicated model of continuum. It should be noted that a definition and explicit characterisation of the discrete elasticity parameters present the key issue of the DEM simulations.

Referring to the above scheme resultant vector

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (3) presents the sum of the globally defined external force [F.sub.i,ext] and centrally interacting interaction particle forces [f.sup.n.sub.ij]. Hereafter, a subscript j denotes six neighbouring particles. Inter-particle forces, acting on the particle i, are illustrated in Fig. 1, d. Hereafter, the inter-particle force vector [f.sup.n.sub.ij] is y restricted to the elastic force.

4. A six-spring model

Development and implementation of the DE model (1) for solids require the evaluation of elastic forces [f.sup.n.sub.ij] in (3), or more precisely, the evaluation of discrete elasticity constants. The continuum properties should be presented in terms of discrete network elements. Actually, the local constitutive relationship defined along the connection line is assumed to be linear and is expressed as

[f.sup.n.sub.ij] = [K.sup.n.sub.ij][h.sub.ij] (4)

where, [h.sub.ij] presents inter-particle displacement as length change of connection line i-j

[h.sub.ij] = [u.sub.nj] - [u.sub.ni] (5)

which is expressed in terms of the local longitudinal displacements [u.sub.i] and [u.sub.j] of the connected nodes i and j of the lattice. Consequently, the elasticity properties are defined by axial stiffness [K.sup.n.sub.ij] of the inter-particle spring.

The idea of the recent development is based on the structural analogy of the continuum in terms of the FEM. Conceptually, each of the elastodynamics problems under arbitrary loading may be solved by FEM. In particular, among various meshes, the periodic triangle lattice may be also considered to be the FEM mesh.

By applying the principle of the virtual work or some other energy method and standard constant strain triangle finite element, it can be easily shown that the resultant Cartesian particle forces [F.sub.i] in Eq. (1) may be expressed in terms of stiffness matrices of the triangles incorporated. This may be done by using a standard FE assembling procedure.

Exploring structural analogy, the internal state of each particular triangle may be defined in terms of the axial forces acting along the triangle sides. The stress-resultant formulation in the finite element analysis is sometimes called a natural approach [26].

In the natural approach, deformation behaviour of triangle element e with local nodes 1, 2 and 3 is described by three independent degrees of freedom and denoted by vector [u.sup.n.sub.e]. Natural displacements mean three stretches of the element sides [u.sup.n.sub.e] = [{[u.sup.n.sub.12], [u.sup.n.sub.23], [u.sup.n.sub.13]}.sup.T], where an individual component actually presents inter-particle displacements [u.sup.n.sub.ij] [equivalent to] [h.sub.ij].

The relationship between natural forces and displacements may be expressed in terms of natural stiffness matrix [[K.sup.n.sub.e]]. Thus, for element e

[f.sup.n.sub.e] = [[K.sup.n.sub.e]][u.sup.n.sub.e] (6)

Elements of the matrix

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

are of different nature. The diagonal members [K.sup.nn.sub.e ii] present the conventional axial stiffness. The off-diagonal members [K.sup.ns.sub.e ij] reflect axial stiffness due to elongations of the neighbouring edges. They express angular deformation or shear effect. All of the stiffness parameters may be directly extracted from the natural stiffness matrix of the triangle. A procedure for the evaluation of natural stiffness matrix [[K.sup.n.sub.e]] = [[K.sup.N.sub.fe]] is given in [27] and briefly described in the Annex by (A1)-(A6).

For the sake of generalisation, the members of the above stiffness may be expressed in terms of the dimensionless coefficients [k.sup.n.sub.ij] as follows

[K.sup.n.sub.ij] = [K.sup.n.sub.ij]Es (7)

The DEM technique operates, however, upon mutual interaction between two particles, therefore, local interactions should be considered.

Let us revise the two-dimensional hexagonal cell depicted in Fig. 2, a. This cell may be regarded as a composition of e (e = 1, ... , 6) triangle finite elements. The mesh of the cell contains j (j = 1, ... , 6) lines intersecting in the node i, where a particular line i-j connects the nodes i and j. This line is considered to be a boundary line between two adjacent elements e and e + 1.

For the triangle e defined by nodes i, j and j-1, the natural force vector is [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.], while for the triangle e+1 defined by nodes i, j and j+1, the force vector is [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.].

Now, the particle equilibrium, derived from the lattice model as illustrated in Fig. 1, d, may be reformulated in terms of the above forces. Actually, the inter particle force [f.sup.n.sub.ij], as shown in Fig. 2, b, is composed of two components reflecting a contribution of two adjacent triangles.

Hence, we may compose

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (8)

As a consequence, it was found that each of the edge forces [f.sup.n.sub.e ij] is composed of the components of different nature [f.sup.nn.sub.e ij] and [f.sup.ns.sub.e ij] involving the influence of the axial deformation due to normal interaction of contacting particles and due to shear caused by interaction of the neighbouring particles.

Finally, the resultant force analogous to the earlier expression (4) may be defined in terms of elongations [h.sub.ij] and six discrete elasticity parameters

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (9)

[FIGURE 2 OMITTED]

Thus we arrived at the model having merely six discrete elasticity constants. It may be classified as a six-spring model with angular interaction [23], expressed, however, in the alternative form.

It is easy to confirm that for anisotropic continuum the elasticity matrix [D] defined by (A3) contains six independent parameters. Therefore, the expression (A6) presents a unique reversible relationship between discrete constants and elasticity constants of the continuum.

For homogeneous isotropic continuum, elasticity properties are defined by two elasticity constants, elasticity modulus E and Poisson's ratio v. Thereby, discrete parameters are expressed as [K.sup.nn.sub.11] = [K.sup.nn.sub.22] = [K.sup.nn.sub.33] = [K.sup.nn] and [K.sup.ns.sub.12] = [K.sup.ns.sub.13] = [K.sup.ns.sub.23] = [K.sup.ns]. In this way, we arrive at a simplified expression of elastic interaction force (9) with two discrete constants:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (10)

Finally, regarding (7), for homogeneous isotropic material two discrete nondimensional elasticity constants are defined explicitly as Poisson's ratio dependant parameters

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (11)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (12)

The variation of constants defined by Eqs. (11), (12) against v is shown in Fig. 3, a. The graph illustrates the sign change of shear stiffness at v = 0.33.

[FIGURE 3 OMITTED]

A significant difference of the present model from the earliest developments [15, 18] is based on a specific approach to shear stiffness. In our approach, the shear stiffness works only in combination with axial stiffness, therefore, resultant stiffness is insensitive to change of sign in shear stiffness. Considering the above results, it may be stated that a six-spring model is physically and mathematically consistent for a larger range of Poisson's ratio, including v [greater than or equal to] 0.33.

Orthotropic material where two principal axes match the coordinate directions is defined by four independent constants. Assuming two independent elasticity modula [E.sub.1] and [E.sub.2], Poisson's ratio v and shear modulus G, the elasticity matrix is presented as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (13)

where two dimensionless parameters are

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

In this case, the discrete constants are defined by four parameters [K.sup.nn.sub.11], [K.sup.nn.sub.22] = [K.sup.nn.sub.33], [K.sup.nn.sub.12] = [K.sup.nn.sub.13] and [K.sup.nn.sub.23]. By assuming n = 2.0 and m = n/(2(1+v)), the influence of orthotropy is illustrated in Fig. 3, b. It is shown that non-positive shear stiffness occurs at smaller values of v.

The developed methodology could be applied to nonhomogeneous solids.

5. Simulation results and discussion

For validating the relevance and accuracy of the developed six-spring DE model and numerical algorithm, the plane stress problem was examined numerically. A sample of the rectangular plate was considered. The geometry of the domain defined by two characteristic dimensions, B and H = 4B, is shown in Fig. 4. The thickness of the plate is s = 0.1 B. The plate boundaries AB and CD are assumed to be clamped by connecting it to rigid walls, while AD and BC are free boundaries. The external excitation is kinematical. It is implemented via the motion of the boundary CD defined by the constant velocity v and the prescribed displacement value [u.sub.max].

[FIGURE 4 OMITTED]

5.1. Validation of discretisation parameters

Dynamic behaviour of the discretised structure is obtained by numerical integration of equations of motion (1) for each individual particle by applying explicit integration with a constant time step. The time step must be sufficiently small in order to ensure the stability of integration. Following the reported recommendations for various DEM models [3], the time step [DELTA]t is defined as a fraction of the critical time step [DELTA][t.sub.cr]

[DELTA]t = 1/[beta] [DELTA][t.sub.cr] (14)

where constant [beta] is fraction factor. The lower bound of this factor is [beta] [greater than or equal to] 10, while the upper bound is evaluated on the basis of numerical experiments.

In turn, the critical time step is related to the period of the highest natural frequency of the smallest particle. Because of the difficulties in precise calculations via eigenvalues of cell parameters, crude estimation in terms of the mass m of the particle and the largest inter-particle stiffness [K.sub.max] = 2[K.sup.nn] may be given as follows

[DELTA][t.sub.cr] = [square root of (m/[K.sup.max])] (15)

Validation was performed by assuming that isotropic elastic material is characterized by the following local properties at an arbitrary point: density of the material is [rho], elasticity modulus is E and Poisson's ratio - v. By substituting mass of the particle (2) and general expression of the inter-particle stiffness (10), the estimate may be scaled with respect to the size of computational domain H as shown in Fig. 2, a. A characteristic lattice size L is related to discretisation density as L = H/nx, where nx is the number of subdivisions along the axis x. Finally

[DELTA][t.sub.cr] = [alpha](v, nx) H/c (16)

where [alpha] is a dimensionless Poisson's ratio-dependent factor, while c = [square root of (E/[rho])] is the material constant, i.e. the speed of sound. For isotropic material maximal stiffness is defined by (11), and factor [alpha] may be expressed explicitly as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (17)

From this expression it follows that, for isotropic material, critical time step would be minimal if Poisson's ratio approaches zero v [right arrow] 0. With the increase of nx, the influence of v vanishes.

After specifying the geometry and the material, a series of numerical tests were conducted for validation purposes. The geometry of the plate is defined by the parameters H = 400 mm, B = 100 mm and s = 10 mm. The isotropic material is characterized by [rho] = 2500 kg/[m.sup.3] and E = 17.1 GPa.

Five schemes were used to evaluate the influence of the characteristic lattice size L. Assuming different numbers of subdivisions [nx.sub.1] = 63, [nx.sub.2] = 126, [nx.sub.3] = 252, [nx.sub.4] = 504 and [nx.sub.5] = 1008, different lattices, having the characteristic sizes [L.sub.1] = 6.349 mm, [L.sub.2] = 3.175 mm, [L.sub.3] = 1.587 mm, [L.sub.4] = 0.794 mm and [L.sub.5] = 0.397 were generated. In the case of the finest lattice, the DEM model contained 291457 particles with 582914 degrees of freedom. It should be noted that the geometry of interacting particles at the boundary is described by cutting the lattice grid.

The influence of the factor [beta] was examined in the first series of numerical experiments. The experiments were conducted according to the second discretisation scheme with [nx.sub.2] = 126 subdivisions, while Poisson's ratio v = 0.2 was assumed. The impact load on the wall CD is induced by prescribed velocity 1.0 m/s.

After calculating the critical time step [DELTA][t.sub.cr] = 601.38 ns according to (15), six values of [[beta].sub.1] = 10, [[beta].sub.2] = 20, [[beta].sub.3] = 40, [[beta].sub.4] = 60, [[beta].sub.5] = 80 and [[beta].sub.6] = 100 yielding six time steps [DELTA][t.sub.1] = 60.138 ns, [DELTA][t.sub.2] = 30.069 ns, [DELTA][t.sub.3] = 15.034 ns, [DELTA][t.sub.4] = 10.023 ns, [DELTA][t.sub.5] = 7.5172 ns, and [DELTA][t.sub.6] = 6.0138 ns were employed for the sake of comparison.

Simulation results yielded actually identical time histories. The influence of the time step was evaluated by comparing the displacement values averaged in the entire loading interval

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]

where, subscript k stands for the time steps (k = 1, q). It was found that the difference for different [beta] was insignificant and varied in the range between 0.135 and 0.136 %.

Time history of transversal displacement [u.sub.E] of the point E is exhibited in Fig. 5.

The accuracy of DEM simulations for this problem was studied by comparing DEM results with finite element simulations. The unique lattice grid serves as the base for both DEM and FEM models. The FEM analysis was performed with the second mesh with [nx.sub.2] = 126 subdivisions by applying the second-order triangles PLANE2 of the ANSYS code [28]. The transient FE problem was solved implicitly. Time step [DELTA]t = 58.3 ns was kept constant during the simulations.

[FIGURE 5 OMITTED]

Qualitatively, the time histories exhibited in Fig. 5 show perfect identity. The comparison of average transversal displacements at point E yields [q.sub.av FEM] = 12.37 [micro]m and [q.sub.av dem] =12.38 um with a sufficiently small 0.136 % difference. A general tendency can be observed that DEM models are stiffer compared to the FE discretisation.

The influence of the mesh refinement was studied in the second series of numerical experiments by solving the above impact problem by using five discretisation meshes and fixed fraction factor [beta] = 40. The simulation results are presented in Fig. 6.

[FIGURE 6 OMITTED]

The results of coarse mesh with nx = 63 and medium size mesh with nx = 252 are presented in Fig. 6. It should be mentioned that the influence of finer meshes is not important. Based on the results obtained, it can be stated that the quality of the developed method could be approved, while the time step factor [beta] = 40, and the mesh with the characteristic refinement parameter [nx.sub.3] = 252 may be recommended for further analysis.

5.2. Diversity of Poisson's ratio

The developed six spring DEM model is applied to improve the deficiency of the earliest approaches [15, 18] for higher values of the Poisson's ratio, v > 0.33. To investigate this influence, the above example of the plate for the selected values v = 0.10, v = 0.20, v = 0.30 v = 0.33 and v = 0.40 were examined. Following the above recommendation, the time step [DELTA][t.sub.3] = 15.034 ns, and the mesh with the characteristic element size L = 3.175 mm were applied to the simulations described below. Along with the DEM simulations, the FE analysis was performed. The small 2.6% difference was obtained.

5.3. Orthotropic material

The suitability of the six-spring DE model for the analysis of orthotropic material was tested by solving the above elastodynamic plane stress problem. Orthotropic material is characterised by four independent elasticity constants, two different elasticity modulus [E.sub.1] = 17.1 GPa and [E.sub.2] = 8.55 GPa, reflecting its different properties in two orthogonal directions, shear modulus G = 7.125 GPa and Poisson's ratio, including v = 0.20. This value yields negative shear stiffness (Fig. 4), therefore, it is suitable for illustrating the possibilities of the developed technique. The elasticity matrix is defined by (7), thus n = 2 and m = 0.833. The time step factor [beta] = 40, and the mesh with the characteristic refinement parameter [nx.sub.3] = 252 are assumed. Simulation results are presented in Fig. 7.

[FIGURE 7 OMITTED]

6. Concluding remarks

The continuum consistent six-spring lattice-based DE model was developed. The original contribution is based on the alternative description of angular interaction defined by applying natural finite element technique.

It has been found that the suggested alternative formulation allows as to present axial and angular interaction by sequentially connected springs and to avoid singularity in a larger range of Poisson's ratio, including v [greater than or equal to] 0.33, for plane stress problem.

The six-spring model with six independent discrete elasticity parameters is applicable to anisotropic material, while the number of independent constants for isotropic material is two.

The suitability of the model was tested by solving the elastodynamic plane stress problem and comparing it to 'exact' FE results. The numerical results show that the refinement of the lattice grid yields actually exponential convergence of DEM model to exact solution.

From the physical point of view, binary interaction is extended to the neighbouring particles within the cell, while the new approach requires a modification of the DEM algorithm and the code.

APPENDIX. Natural stiffness of triangle.

Geometry of the triangle element fe is defined by three angles [alpha], [beta] and [gamma] and three edges with the lengths [L.sub.[alpha]], [L.sub.[beta]], [L.sub.[gamma]]. Its deformation behaviour is described by three independent deformation-related degrees of freedom [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.]. In natural approach, they mean three stretches of the element. Generally, the element may be arbitrary oriented in the coordinate plane (Fig. A1).

[FIGURE A1 OMITTED]

The symmetric natural (3 x 3) stiffness matrix of the element is denoted hereafter as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A1)

According to [25], it s defined explicitly as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A2)

here, the constants [A.sub.fe] and s stand for the triangle area and the element thickness, respectively.

The constitutive properties of planar solid are defined by elasticity matrix as follows

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A3)

reflecting isotropic and anisotropic properties of planar solid.

The diagonal matrix [[l.sub.fe]]= diag[[L.sub.[alpha]] [L.sub.[beta]] [L.sub.[gamma]]] presents side lengths. The transformation matrix [[C.sub.fe]] relates natural and Cartesian strains. It is expressed in terms of direction cosines [c.sub.x[alpha]] = cos(x, [L.sub.[alpha]]) and [c.sub.y[alpha]] = cos(y, [L.sub.[alpha]]) of the triangle sides with respect to Cartesian axes x and y (Fig. A1). It reads

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A4)

For equilateral triangle angles [alpha] = [beta] = [gamma] = 2[pi]/3, the lengths [L.sub.[alpha]]= [L.sub.[beta]] = [L.sub.[gamma]] = L are equal, while the area is [A.sub.fe] = [square root of (3[L.sup.2]/2)]. Assuming the coordinate axis Ox along the side, the direction cosines are expressed as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A5)

Finally, discrete elasticity constants are given as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII.] (A6)

General expressions (A6) may be easily applied to particular materials, when the elasticity matrix (A4) is defined explicitly.

Received November 19, 2008

Accepted February 25, 2009

References

[1.] Belytschko, T., Liu, W.K., Moran, B. Nonlinear Finite Elements for Continua and Structures. -John Wiley & Sons, 2001.-650p.

[2.] Cundall, P.A., Strack, O.D.L. A discrete numerical model for granular assemblies. -Geotechnique, 1979, v.29, No.1, p.47-65.

[3.] Dziugys, A., Peters, B. An approach to simulate the motion of spherical and non--spherical fuel particles in combustion chambers. -Granular Matter, 2001, 3, p.231-265.

[4.] Markauskas, D, Kacianauskas, R. Compacting of particles for biaxial compression test by the discrete element method. -Journal of Civil Engineering and Management, 2006, 12(2), p.153-161.

[5.] Maknickas, A., Kaceniauskas, A., Kacianauskas, R., Balevicius, R., Dziugys, A. Parallel DEM Software for Simulation of Granular Media. -Informatica, 2006, 17(2), p.207-224.

[6.] Balevicius, R., Kacianauskas, R., Dziugys, A., Maknickas, A., Vislavicius, K. Investigation of performance of programming approaches and languages used for numerical simulation of granular material by the discrete element method. -Computer Physics Communications, 2006, 175(6), p.404-415.

[7.] Dziugys A., Navakas, R. The role of friction on size segregation of granular material. -Mechanika. -Kaunas: Technologija, 2007. Nr.4(66), p.59-68.

[8.] Mishra, B.K., Thornton, C. Impact breakage of particle agglomerate'. -Int. J. Miner. Process, 2001, 61, p.225-239.

[9.] 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, 4, p.77-90.

[10.] Antonyuk, S., Khanal, M., Tomas, J., Heinrich, S. Morl, L. Impact breakage of spherical granules: Experimental study and DEM simulation. -Chemical Engineering and Processing, 2006, 45, p.838-856.

[11.] Hrennikoff, A. Solution of problems of elasticity by the framework method. -J. of Appl. Mech. 8, 1941, p.169-175.

[12.] Kawai, T. New discrete models and their application to seismic response analysis of structures, -Nucl. Eng. Des., 1978, 48, p.207-229.

[13.] Herman, H.J., Hansen, A., Roux, S. Fracture of disordered, elastic lattices in two dimensions. -Phys. Rev., 1989, B. 39(1), p.637-648.

[14.] Sawamoto, Y., Tsubota, H., Kasai, Y., Koshika, N., Morikawa, H. Analytical studies on local damage to reinforced concrete structures under impact loading by discrete element method. -Nuclear Engineering and Design, 1998, 179, p.157-177.

[15.] Griffiths, D.V., Mustoe, G.G.W. Modelling of elastic continua using a grillage of structural elements based on discrete element concepts. -International Journal for Numerical Methods in Engineering, 2001, 50, p. 1759-1775.

[16.] Liu, K., Gao, L. The application of discrete element method in solving three-dimensional impact dynamics problems. -Acta Mechanica Solida Sinica. 2003, v.16, No.3, p.256-261.

[17.] Ricci, L., Nguyen, V.H., Sab, K., Duhamel, D., Schmitt, L. Dynamic behaviour of ballasted railway tracks: A discrete/continuous approach. -Computers and Structures, 2005, 83, p.2282-2292.

[18.] Liu, K., Liu, W. Application of discrete element method for continuum dynamic problems. -Archive of Applied Mechanics, 2006, v.76, No.3-4, p.1-8.

[19.] Tavarez, F.A., Plesha, M.E. Discrete element method for modelling solid and particulate materials. -Int. J. Num. Meth. Eng., 2007, 70, p.379-404.

[20.] Lilliu, G., Van Mier, J.G.M. 3D lattice type fracture model for concrete. -Engineering Fracture Mechanics, 2003, 70, p.927-941.

[21.] Karihaloo, B.L., Shao, P.F., Xiao, Q.Z. Lattice modelling of the failure of particle composites. -Eng. Fract. Mech., 2003, 70, p.2385-2406.

[22.] Ibrahimbegovic, A., Delaplace, A. Microscale and mesoscale discrete models for dynamic fracture of structures built of brittle material. -Computers and Structures, 2003, 81, p.1255-1265.

[23.] Ostoja-Starzewski, M. Lattice models in micromechanics. -Appl. Mech. Rew., 2002, 55(1), p.3.

[24.] Vadluga, V., Kacianauskas, R. Investigation of the single-spring lattice model in simulation of 2D solid problems by DEM. -Mechanika. -Kaunas: Tech nologija, 2007, Nr.5(67), p.5-12.

[25.] Zhang, Z., Ge, X. Micromechanical modelling of elastic continuum with virtual multi-dimensional internal bonds. -Int. J. Numer. Meth. Engng, 2006, 65, p.135-146.

[26.] Argyris, J., Balmer, H., Doltsinis, J.St., Dunne, P.C., Haase, M., Kleiber, M., Malejannakis, G.A., Mlejnek, H.-P., Muller, M. Scharpf, D.W. Finite element method--the natural approach.-Comp. Meth. Appl. Mech. Eng., 1979, 17/18, p.1-107.

[27.] Vadluga, V. Simulation of dynamic deformation and fracture behaviour of heterogeneous solids by discrete element method. -Doctoral Thesis, Vilnius Gediminas Technical University. 2008.-99p.

[28.] Moaveni, S. Finite Element Analysis. Theory and Applications with ANSYS--Saeed Moaveni. -New Jersey: Prentice Hall, Upper Saddle River, 07458, 1999.-527p.

R. Kacianauskas *, V. Vadluga **

* Vilnius Gediminas Technical University, Sauletekio al. 11, 10223 Vilnius-40, Lithuania, E-mail: rkac@fm.vgtu.lt

** Vilnius Gediminas Technical University, Sauletekio al. 11, 10223 Vilnius-40, Lithuania, E-mail: vvad@fm.vgtu.lt