Geometric modelling of elastic and elastic-plastic solids by separation of deformation energy and Prandtl operators

A geometric method for analysis of elastic and elastic-plastic solids is proposed. It involves operators on naturally discrete domains representing a material’s microstructure, rather than the classical discretisation of domains for solving continuum boundary value problems. Discrete microstructures are considered as general cell complexes, which are circumcentre-dual to simplicial cell complexes. The proposed method uses the separation of the total deformation energy into volumetric and distortional parts as a base for introducing elastoplastic material behaviour. Volumetric parts are obtained directly from volume changes of dual cells, and the distortional parts are calculated from the distance changes between primal and dual nodes. First, it is demonstrated that the method can accurately reproduce the elastic behaviour of solids with Poisson’s ratios in the thermodynamically admissible range from -0.99 to 0.49. Further veri-ﬁcation of the method is demonstrated by excellent agreement between analytical results and simulations of the elastic deformation of a beam subjected to a vertical displacement. Second, the Prandtl operator approach is used to simulate the behaviour of the solid during cyclic loading, considering its elastoplastic material properties. Results from simulations of cyclic behaviour during alternating and variable load histories are compared to expected macroscopic behaviour as further support to the applicability of the method to elastic-plastic problems.


Introduction
When a solid body is subjected to mechanical loading, it undergoes volume and shape changes.Mechanical loads that act on a solid can appear in simple forms, such as tension, bending, shear, torsion; or they can act as a complex combination of the simple loads causing an intricate deformation field.Methods to analyse the deformation of solids, analytical and numerical, have been extensively studied to date.Analytical solutions, mostly based on the seminal works of Timoshenko and Goodier (1951) , can quickly become computationally too complex for the geometries of mechanical products developed today.Finite difference methods are simple and intuitive numerical techniques for the solution of mechanical problems ( Abreu et al., 2018;Li et al., 2019;Zheng and Zhang, 2019 ), but again the complex geometries can compromise the accuracy of the results depending on the size of the chosen perturbation ( Abreu et al., 2018 ).Finite volume methods ( Keilegavlen and Nordbotten, 2017;Sokolova et al., 2019 ) and boundary-type meth-ods ( Markous, 2019;Liu et al., 2019 ) can alternatively be applied to study the deformation of solids.At present, the finite element methods are the most mature and commonly used tool to analyse solids subjected to mechanical or thermomechanical loadings ( Liu et al., 2015;Fielder et al., 2017;Rajasekaran and Khaniki, 2017;Barboura and Li, 2018;Šeruga and Nagode, 2019;Kotoul et al., 2019;Xu and Cheng, 2019;Xiao and Yu, 2019;Pereira et al., 2019;Numanoglu and Civalek, 2019;Falope et al., 2019 ).All these methods are based on continuum formulations of governing equations, i.e. by partial differential equations, complemented by boundary conditions defining boundary value problems.The numerical methods then proceed to discretise the domain representing the body and use approximations of specific mechanical fields or differential or integral operations.Solutions of such formulated boundary value problems are then sought by variational principles and iterative procedures ( Zienkiewicz et al., 2013;Šeruga and Nagode, 2019;Rajasekaran and Khaniki, 2017;Xu and Cheng, 2019 ).
One limitation to the methods starting from continuum formulations is that they are not suitable for the analysis of localised phenomena in solids, such as the formation of shear bands or cracks.This is because such phenomena are physically of finite extensions, whilst the classical solid mechanics, where stress and strain tensors are point-wise quantities, lacks the length scale.A relatively simple alternative allowing for localisation as well as heterogeneity to be introduced is offered by the lattice models.For example, lattice methods based on the analysis of graphs have been used to simulate the elastic behaviour, damage initiation and fracture propagation of various materials subjected to mechanical loads ( Jivkov and Yates, 2012;Jivkov et al., 2012;2013;Morrison et al., 2014;2016;Zhang et al., 2014;Dassios et al., 2017;2018 ).However, various discrete formulations have been proposed for achieving arbitrary values of Poisson's ratio, especially for materials in the elastoplastic regime where Poisson's ratio asymptotically approaches 0.5 (incompressibility) or to reproduce the behaviour of auxetic materials with a negative Poisson's ratio.Kumar et al. (2016) studied simulations by a packing of bonded particles to provide the microscopic and macroscopic characteristics of the packing and explored the capability of a discrete element method with spheres to reproduce the behaviour of auxetic materials.Rojek et al. (2018) introduced a discrete element method with deformable cylindrical particles which broadens the range of the macroscopic Poisson's ratio in comparison to the values achievable by the standard discrete element method.Zhao et al. (2011) proposed a 3D distinct lattice spring model where matter is discretised into individual particles linked by springs forming either a random or a regular lattice structure and capable of representing a diversity of Poisson's ratios.Asahina et al. (2015) demonstrated that it is possible to correctly simulate the Poisson effect by a regular lattice model and a set of fictitious forces introduced into the lattice.Furthermore, Asahina et al. (2017) presented an auxiliary stress approach which is capable of representing the Poisson effect both macroscopically and in a local sense.Celiguet et al. (2017) proposed a way to enrich the spring dash-pot model in such a way that the macroscopic Poisson effect can be captured by the discrete element method.Cusatis et al. (2017) introduced a discontinuous cell method which enables simulations of a full range of Poisson's ratio from −1 to 0.5.Analysis of graphs, applicable to lattices, is a geometrically simplified version of a more general formulation using discrete exterior calculus (DEC) ( Hirani, 2003 ).DEC utilises the interaction between discrete primal and dual complexes to reproduce the key geometric features useful for computational purposes.An important feature of DEC is that both the intermediate discretisation of the continuum governing equations and their application to the discretised domain become superfluous ( Yavari, 2008;Kosmas and Jivkov, 2019 ).This formulation is currently receiving increasing interest after its successful application to fields within both physics and mathematics.Specifically, it has been used to tackle problems ranging from homology, riemannian geometry, fluid dynamics and discrete mechanics ( Schulz and Tsogtgerel, 2018 ).The application to elasticity problems is still under development.However, by introducing primal and dual complexes, boundary and co-boundary operators for the connectivity of these complexes and Hodge star operators between them, it is in principle possible to solve elasticity problems of solids based only on their geometric properties ( Yavari, 2008;Bell and Hirani, 2012;Kosmas and Jivkov, 2019 ).The solution for a given structural problem is then found by the minimisation of the deformation energy which accumulates in a solid as a reaction to the external loading ( Kosmas and Jivkov, 2019 ).This is a remarkable feature of the approach as it can prevent escalating computational errors during time-dependent loading.
This paper offers an approach, denoted as lattice-complex modelling, that combines the more established lattice approach with elements of the emerging full DEC formulation.A new computational method is developed, where solids are represented by discrete topological spaces and their behaviour is analysed using only geometric features.The positions of the primal and the dual nodes and the stretches between them define the stored deformation energy in the dual complex, considered as a finite discrete structure of the solid.The total deformation energy in a discretised solid is calculated by separate volumetric and distortional contributions, which enables simulations of the elastic behaviour of solids.Only the positions of the primal and the dual nodes are crucial for the simulation regardless of their translational or rotational movements.The method then extends the use of geometric modelling to analysis of the elastoplastic behaviour of solids, which has not been reported to date.This is enabled by the introduction of the Prandtl operator approach ( Nagode and Zingsheim, 2004;Nagode et al., 2010;2011a ).The Prandtl operator approach enables the simulation of elastoplastic stress-strain behaviour and the kinematic hardening of solids during variable mechanical loadings.The main advantages of the approach are high computational speeds and omission of the explicit use of the rules for the simulation of cyclic behaviour as these rules are integrated within the model ( Šeruga and Nagode, 2019 ).Detailed descriptions of the Prandtl operator approach development can be found in Nagode and Zingsheim (2004) ; Nagode and Fajdiga (2007) ; Nagode et al. (2011b) ; Šeruga et al. (2014) ; Šeruga and Nagode (2019) .In addition, some applications of the Prandtl operator approach on a variety of mechanical components can be found in Nagode et al. (2011b) ; Šeruga et al. (2014) ; Kozjek et al. (2017) ; Šolinc et al. (2019) ; Bartošak et al. (2019) .
The application of the proposed method is aimed at metallic mechanical components operating in demanding operational conditions under considerable mechanical and thermal loads, e.g. automotive exhaust systems, turbines or aircraft engines ( Šeruga et al., 2014;Bewlay et al., 2016;Abdallah et al., 2012;Hayakawa et al., 2004;Thomas and Bacos, 2011 ).Stresses and strains arising in the material of such components as a reaction to external loads reach beyond the yield stress causing elastoplastic behaviour of the material.To endure such operating conditions, a considerable amount of chromium, nickel, tungsten and cobalt is needed in conventional alloys ( Šeruga and Nagode, 2015;Grilli et al., 2017;Omari and Sevostianov, 2013;Tkaczyk et al., 2018;Mishakin et al., 2019 ).However, as most of these alloying elements are termed as critical raw materials (CRM) for Europe, only functional recycling of used products containing CRM and alternative CRM-reduced or CRMfree materials will enable the manufacture of mechanical components able to withstand extreme operating conditions in the future ( Grilli et al., 2017;Tkaczyk et al., 2018;Novak et al., 2018 ).Therefore only a robust method providing reliable predictions of the behaviour of solids in operation will improve their performance, reduce the consumption of critical raw materials or predict suitable substitutes.
The paper is structured as follows: Section 2 provides the formulation of geometric modelling for mechanically loaded elastic and elastic-plastic solids.In Section 3 , results for a loaded icosahedron, a rigidly supported beam and an elastoplastic behaviour under variable mechanical loads are given to demonstrate the use of the method.Section 4 discusses the results and provides the future work regarding the method.

Mathematical preliminaries
A solid body is considered as a collection of polyhedral cells, representing the individual grains of its polycrystalline structure.This view allows for consideration of the nuclei, around which grains have formed by solidification, as nodes ( x i , y i , z i ) for the Delaunay tetrahedralisation of the domain ( Si, 2015 ).The nodes, edges, triangular faces and tetrahedra after tetrahedralisation form a simplicial complex, which we consider as a primal complex in the present development.In the terminology of algebraic topology, the elements of the simplicial complex are denoted: 0-cells for nodes, 1-cells for edges, 2-cells for triangles, and 3-cells for tetrahedra.The Voronoi tessellation ( Si, 2015;Voronoi, 1908 ) around the primal 0-cells provides a general (non-simplicial) complex, which is circumcentre-dual to the primal complex.This dual complex is then considered to represent the polycrystalline structure of the solid.To each primal 0-cell, denoted by σ 0 , corresponds one dual 3-cell, denoted by * σ 0 (e.g. a single crystal); to each primal 1-cell, denoted by σ 1 , corresponds one dual 2-cell, denoted by * σ 1 (generally a planar polygon, e.g. a grain boundary, normal to the primal 1-cell); to each primal 2-cell, denoted by σ 2 , corresponds one dual 1-cell, denoted by * σ 2 (e.g. a triple line in the polycrystalline structure, normal to the primal 2-cell); to each primal 3-cell, denoted by σ 3 , corresponds one dual 0-cell, denoted by * σ 3 (e.g. a quadruple point of the polycrystalline structure).The connectivity between the elements of the primal complex is described by three matrices, providing the 0-cells that are boundaries of 1-cells, the 1-cells that are boundaries of 2-cells, and the 2-cells that are boundaries of 3-cells, respectively.Details on how these are constructed can be found in Bell and Hirani (2012) .The transposes of these boundary operators provide co-boundary operators, denoted by d 0 , d 1 and d 2 , respectively, which are discrete analogues of the exterior derivatives in the continuum exterior calculus ( Hirani, 2003 ).Due to the one-to-one correspondence between σ i and * σ i , the boundary operators of the dual complex are the transposes of the primal boundary operators, and the dual co-boundary operators are simply d T 0 , d T 1 and d T 2 .This formulation leads to a discrete analogue of the de Rham complex depicted in Fig. 1 .
The illustration in Fig. 1 shows that three algebraic quantities, the matrices d 0 , d 1 and d 2 arising from the topology of the primal complex, are sufficient to describe the connectivity as well as the discrete derivatives in both the primal and the dual complexes.The specific geometric properties of the two complexes, such as   Hirani, 2003 ).
Using these notations, the following operations can be performed: • the gradient of a function on σ 0 is found by the application of * 1 • d 0 and is a function on * σ 1 , • the curl of a function on σ 1 is found by the application of * 2 • d 1 and is a function on * σ 2 and • the divergence of a function on σ 2 is found by the application of * 3 • d 2 and is a function on * σ 3 .Similar operations are performed for functions defined on * σ 3 , * σ 2 , and * σ 1 with inverse Hodge stars and the transposed coboundary operators.Hence discrete analogues of all the operations used in continuum physics and mechanics can be built, and as mentioned in the introduction this method has been used to construct discrete formulations of many physical problems ( Schulz and Tsogtgerel, 2018 ).However, a direct application of the DEC basis to elasticity problems is still elusive due to the specific requirements for the symmetry of strain and stress tensors in continuum mechanics, which cannot be ensured by the existing DEC apparatus.Therefore, we propose here to construct an approximate method, which takes ideas from the DEC formalism and the more explored mechanics of lattices.

Lattice-complex formulation of elastic behaviour
Classical lattice models can be seen as a part of the primal complex defined above, containing only the sets of 0-cells and 1-cells, i.e. the centres of physical grains in a polycrystalline structure and the segments connecting centres of neighbouring grains ( Jivkov and Yates, 2012;Morrison et al., 2016;Zhang et al., 2014;Dassios et al., 2017;2018 ).One problem with lattice formulations is that they do not allow for unique separation between the stretching and rotational components of deformation, which is straightforward in continuum formulations -either by additive decomposition of the displacement gradient into symmetric and skewsymmetric tensors in infinitesimal strain theories, or by multiplicative decomposition of the deformation gradient into stretch and rotation tensors in finite strain theories.This leads to an inability to construct a frame-indifferent formulation for the whole lattice and one consequence is that such models cannot reproduce bulk elastic behaviour with all thermodynamically admissible Poisson's ratios.A frame-indifferent formulation can be achieved by using scalar invariants, specifically the components of the stored elastic energy.However, the classical lattice models contain insufficient information to capture the two principal contributions -volumetric and distortional -arising from volume change and shape change of the physical grains.The proposal in this work is to use a complementary structure for evaluation of the volumetric and distortional energies.The dual 3-cells (generally polyhedrons), indexed by * σ 0 i , i = 1 , . . ., N, where N is the number of dual 3-cells, have individual N i vertices, which are dual 0-cells indexed by * σ 3 j , where j spans the set of N i indices of the dual 0-cells on the boundary of * σ 0 i .Each dual 3-cell contains a unique part of the proposed complementary structure -a set of interior segments e ij connecting its corresponding primal 0-cell, σ 0 i , which is in the * σ 0 i interior by construction, to its vertices * σ 3 j .This is illustrated in Fig. 2 for an arbitrary cross-section.The calculation of the volumetric and distortional strain energies in individual dual 3-cells is then based on the changes of their volumes and stretches of thus defined interior segments.
For given boundary conditions, the body is assumed to be in an initial equilibrium configuration, described by the coordinates of σ 0 i denoted by X i with respect to some global coordinate system.The coordinates of * σ 3 j denoted by X j are calculated by construction of the dual complex.A change of boundary conditions, e.g.due to increments of prescribed forces and/or prescribed displacements, leads to a new, current configuration.The coordinates of σ 0 i and * σ 3 j in the current configuration are denoted by x i and x j , respectively, where the former form the set of unknowns for the problem and the latter are functions of the former, derived by the construction of the dual.
The ratios between the current volumes, | * ˜ σ 0 i | , and the initial volumes, | * σ 0 i | , of the dual 3-cells represent discrete volume changes corresponding to the Jacobian of the deformation gradient in continuum elasticity ( Bonet and Wood, 2008 ): (1) The interior segments, e ij , in a dual 3-cell * σ 0 i have initial and where their stretches are calculated as λ i j = . However, each stretch can be connected to the volumetric change J i , i.e. expansion or contraction of a volume, and the shape change parameter δ ij within the volume, (2) Considering Eq. 2 , the shape change of e ij can be calculated as The stretch of an interior segment e ij can always be evaluated if both the primal 0-cell σ 0 i and the dual 0-cell * σ 3 j exist.Only if either of the quantities ceases to exist, e.g. if a primal 3-cell σ 3 is deleted, then the evaluation of the stretch is not possible.
The volumetric energy i of the dual 3-cell * σ 0 i can thus be expressed as where κ 0 represents the elastic bulk modulus of the material, Parameters E and ν correspond to the elastic properties of the material under consideration.The distortional energy i within the dual 3-cell * σ 0 i is worked out from all the shape change parameters δ ij of the interior segments e ij as i = 3 and μ 0 stands for the elastic shear modulus, Finally, the total deformation energy U is calculated as the sum of volumetric and distortional contributions within the dual complex, In order to calculate the displacements of the primal complex subjected to mechanical loading, the total deformation energy is minimised using a numerical routine, (9) In the present context, the minimisation of the total deformation energy is achieved by the quasi-Newton method of Nocedal and Wright (2006) .The algorithm starts the search for the optimal coordinates of primal 0-cells in the current configuration at an initial estimate x i 0 and iterates until the final position of the primal 0-cells x i is achieved.In the present work, the combinatorial structure of the primal cell complex is immutable, i.e. there is no emergence or destruction of cells of any order and the relations between cells do not change.Clearly, this is also true for the structure of the dual cell complex.This is suitable for discrete modelling of elasticity and plasticity.However, modelling of separation processes, e.g.fracture, will require changes in the combinatorial structure of the complex.In such cases, cells of different order, including tetrahedra may vanish, but the analysis will proceed using updated co-boundary operators and the energy minimisation will be performed on the new cell complex.During each iteration, the coordinates of the dual 0-cells are determined as a function of the coordinates of the primal 0-cells.The search in the k th iteration follows the direction p k determined as where H k and ∇U k are the inverse Hessian approximation and the gradient of the total deformation energy in the k th iteration, respectively.A new estimate of the coordinates of the primal 0-cells in the current configuration is then calculated as where α k is the step length in the k th iteration.The minimisation algorithm is completed when the norm of the gradient || ∇U k || reaches the preset convergence value ε.
The elasticity formulation of lattice-complex modelling is tested in Section 3.1 .

Lattice-complex formulation of elastoplastic behaviour
In order to include plasticity in the lattice-complex modelling, the time t in distinct steps k = 1 , . . ., N k is introduced to enable a time-dependent application of the mechanical load.As the stretch of an interior segment λ ij depends only on the applied load, it is decomposed into volumetric and distortional stretches also beyond the yield point of the material according to Eqs. (1) and (3) .In order to consider the nonlinear effect of plasticity we suggest the adaptation of both the bulk modulus and the shear modulus of the material so as to enable the simulation of memory rules expected during the cyclic loading.For this purpose, one rheological spring-slider model schematically given in Fig. 3 is placed into each dual 3-cell to simulate the irreversible volumetric energy dissipation and another rheological spring-slider model is placed onto each interior segment to simulate the irreversible distortional energy dissipation.
The bulk modulus of the material for an individual dual 3-cell can thus be modelled as a Prandtl type operator, defined by where β κn , I κn , ε κn ( t k ) and ε O κn ( t k ) are the bulk Prandtl density, the active bulk play index, the bulk back-strain and the bulk backstrain origin, respectively.The bulk Prandtl densities reflect the volumetric material properties of the dual complex.Their determination is given in Section 2.4 .The active bulk play index I κn is determined by and marks whether the n th bulk back-strain changed position.The bulk back-strain ε κn ( t k ) provides information on the position of the n -th spring-slider segment in the k th time step as where ε V ( t k ) is the volumetric strain in the k th time step, and q n is the n -th yield strain used during the discretisation of the material properties ( Fig. 3 and Section 2.4 ).The bulk back-strain origin ε O κn ( t k ) follows the memory rules for elastoplastic simula- tion and is always set so as to identify the active hysteresis loop ( Nagode and Šeruga, 2016 ) distinguishes the position on the cyclic curve in the first case and the branch of the hysteresis loop in the second case where ε V,max stands for the absolutely highest strain reached up to time t k −1 .
Similarly, the shear modulus of the material can be modelled for an individual interior segment as a Prandtl type operator, defined by Now, the shear back-strain ε μn ( t k ) is given as where ε D ( t k ) is the distortional strain in the k th step, and the active shear play index I μn ( t k ) is determined as Coefficient ζ in Eq. 17 is defined with respect to ε D (t k ) and ε D , max .The volumetric and distortional energies in the case of elastoplasticity are calculated according to Eqs. ( 4) and ( 6) , respectively, with the elastic bulk modulus κ 0 substituted for the bulk modulus κ( t k ) and the elastic shear modulus μ 0 substituted for the shear modulus μ( t k ).
The elastoplasticity formulation of lattice-complex modelling is tested in Section 3.2 .

Determination of material properties
The properties of the material needed for the simulation can be determined from uniaxial tensile and cyclic tests.From the tensile test, the elastic modulus and the Poisson's ratio can be measured ( Šeruga et al., 2019a ) whereas an incremental step test will provide experimental data for the determination of the cyclic curve ( Šeruga et al., 2019b ).Cyclic bulk modulus-strain and cyclic shear modulus-strain curves can then be created as shown in Fig. 4 .
For discretisation purposes, the yield strains q n depicted in Fig. 4 are first either equidistantly or non-equidistantly spaced.Then, bulk Prandtl densities in the case of equidistantly spaced yield strains are calculated from the cyclic bulk modulus-strain curve as with κ 0 = κ 1 and shear Prandtl densities in the case of equidistantly spaced yield strains q n are calculated from the cyclic shear modulus-strain curve as with μ 0 = μ 1 .In the case of non-equidistantly spaced yield strains, a different expression for the determination of the Prandtl densities is required, e.g. as given in Šeruga et al. (2014) .

Results
The method was applied to a titanium aluminide alloy Ti-42Al-8.5Nbused for applications in high temperature and high mechanical load areas.The bulk modulus-strain and shear modulusstrain curves were determined from experimental data gained in Appel et al. (2014) .They are depicted in Fig. 5 .Ten equidistantly spaced segments for both the bulk and the shear moduli have been used for discretisation of the curves although an arbitrary number of segments can be applied in general.The calculations have been performed using PyDEC library ( Bell and Hirani, 2012 ) for the determination of primal and dual complexes and corresponding geometrical quantities, supplemented by an in-house Python software for the calculation of stored energies.Although available software for the Delaunay triangulation, Voronoi tessellation and the discrete exterior calculus operators on orthogonal simplicial primal and general dual complexes has been used ( Si, 2015;Bell and Hirani, 2012 ), the lattice-complex modelling is not bonded to the Delaunay triangulation and the Voronoi tessellation.Another type of dual topological quantity could be used instead as long as the positions of the primal and the dual 0-cells are known.The minimisation of the total energy was achieved by the quasi-Newton method of Nocedal and Wright (2006) with the preset gradient norm convergence value set to ε = 10 −2 .

Simulations of elastic behaviour
First, the proposed geometric formulation of elasticity was tested on a simple structure -an icosahedron consisting of 13 primal 0-cells and 20 primal 3-cells was chosen.Both primal and dual complexes for the icosahedron are depicted in Fig. 6 .The icosahedron was supported on one edge and the loading displacement was applied on the opposite edge as shown in Fig. 6 .The input Poisson's ratio, given as a material property to calculate both the bulk and the shear moduli, was varied between −0.99 and 0.49 whereas the elastic modulus was kept at 161538 MPa.After the simulation, the calculated displacements in directions perpendicular to the direction of loading were transformed into the calculated Poisson's ratio and compared to the input Poisson's ratio.The results are given in Fig. 7 .It can be seen that the geometric latticecomplex model can represent the whole range of thermodynamically admissible Poisson's ratios with high accuracy.This is the first demonstration of the model capability, the outcome is discussed further in Section 4 .
Second, a rigidly supported beam was subjected to a vertical displacement at one of the supports.The beam consisted of 36 primal 0-cells and 48 primal 3-cells; the primal and dual complexes are given in Fig. 8 .Elastic modulus and Poisson's ratio were set to 161538 MPa and 0.32, respectively.The results of the simulation using lattice-complex modelling were compared to an analytical solution of the same rigidly supported beam,    where I x is the cross-sectional moment of inertia for a rectangular cross-section of the beam, D is the displacement of the support, L is the length of the beam between the supports and K is the rigidity of the beam, given as Furthermore, a simulation using the finite element method was performed under the same conditions.The results of the compar-ison are given in Fig. 9 .These provide further verification of the model, to be discussed in Section 4 .

Simulations of cyclic elastoplastic behaviour
The response of the system to cyclic mechanical loads was simulated by changing the bulk and the shear moduli.The moduli for the titanium aluminide alloy Ti-42Al-8.5Nbas given in Fig. 5 were considered.The results of these simulations are given in Fig. 10 for an alternating mechanical load history and in Fig. 11 for a variable mechanical load history.

Discussion
Conventional material models establish the elastoplastic behaviour of solids during cyclic loading based on stresses and strains whose calculation is an essential step required for simulation of the deformed structure.Namely, such models first require the discretisation of both the geometry and the governing equations for continuum and then search for equilibrium within the solid iterating the displacement-strain-stress-force loop.The intention of the lattice-complex modelling is to discretise the governing equations, given a discrete (and finite) formulation of the solid internal structure.The deformation of the structure is found as the equilibrium state of the dual complex for which the stored deformation energy is a minimum.This gives a significant advantage to the lattice-complex modelling, because numerical errors which can arise by conventional iterative schemes are minimised by requiring the preservation of the minimum accumulated total energy.Classical variational approaches to mechanics satisfy local conservation laws on average and preserve important invariants, although the underlying structures of the simulated continuous systems can often be compromised ( Desbrun et al., 2008 ).For an elastoplastic structural analysis using a finite element method, simulation results might differ depending on the integration procedure, i.e. whether an implicit or an explicit Euler method is applied.Using an explicit solver, a compromise refers to stability issues of the simulation which might amplify the stored energy in the solid.The stability of the simulation ensured by the implicit integration procedure can compromise the preservation of energy due to avoidance of numerical divergence.Another example is given in Stern and Desbrun (2006) for a pendulum which can gain or lose energy depending on the integration scheme.Discrete computational modelling, with which the lattice-complex modelling is associated, operates with forces and displacements at their geometrically-meaningful locations (primal and dual 0-cells) and preserves the minimum stored total deformation energy within the solid in every time step (at present by the quasi-Newton method of Broyden, Fletcher, Goldfarb and Shanno).The principle of energy minimisation is hence applied in both the classical variational approaches and the discrete modelling ( Yavari, 2008;Kanso et al., 2007 ).However two pivotal steps have to be performed so that applicable lattice-complex modelling is achieved: a plausible description of the energy within the complex, and a powerful total energy minimisation method.Here, we propose the description of the energy in terms of the bulk modulus, the shear modulus and the stretch due to the external load which is consistent with e.g.Bonet and Wood (2008) where energy functions also include these terms.Moreover, the total energy is separated into volumetric and distortional parts as mostly the distortional part of the energy is responsible for the plastic deformation of metals whereas both parts of the energy influence the elastic deformation of a solid.The latter can be obtained from Figs. 4 and 5 where bulk and shear moduli retain a constant value within the elastic area.In the elastoplastic area, the bulk modulus starts increasing whilst the shear modulus starts decreasing indicating the resistance of the material to volume change and its susceptibility to shape change.The separation of the energies also enables the introduction of the Prandtl operators and consequently the continuous modelling of both the bulk and the shear moduli throughout the cyclic loading history.
However, the volumetric and the distortional energies have to be expressed so as to enable such deformations of structures after the minimisation of their total energy as would be observed if they were physically tested under the same conditions.First, this means that all the primal 0-cells of the complex have to achieve the correct displaced positions.Second, the Poisson's effect of compressible solids has to be represented -expansion or contraction perpendicular to the direction of the applied load.In order to validate the calculation of the energies proposed here, an icosahedron has been tested initially for the thermodynamically admissible range of Poisson's ratios between -0.99 and 0.49 although most structural materials exhibit a non-negative Poisson's ratio.As can be seen from the results in Fig. 7 , the simulated Poisson's ratio using the lattice-complex modelling is in good agreement with the theoretical values.As anticipated, there are some discrepancies between the simulated and the theoretical values, which do not exceed a value of 0.08.Comparing the undeformed and the deformed shape of the icosahedron, the simulated deformation for the Poisson's ratio ν = −0 .99 clearly shows how the solid shrinks under a compressive load as opposed to the simulated deformation for the Poisson's ratio ν = 0 .48 where the solid expands under the same compressive load.For ν = 0 .0 , there is almost no perpendicular displacement when the compressive load is applied.
The results using the lattice-complex modelling of the rigidly clamped beam subjected to the displacement of one of the supports match well with the results obtained using both the analytical solution and the finite element method ( Fig. 9 ).Apart from slightly higher discrepancies of the obtained vertical displacements around the supports, the result is a confirmation of the acceptability of the proposed method.Nevertheless, a shortcoming of the present method in the case of the rigidly clamped beam turned out to be the energy minimisation method.Namely, the quasi-Newton method of Broyden, Fletcher, Goldfarb and Shanno showed a tendency to find local total energy minima.For the clamped beam this meant possible solutions of the deformed shape having a low value of the total accumulated energy, but not the minimum.In this case, a different initial set of displacements would help the algorithm to find the global total energy minimum.Additionally, the convergence rate of this minimisation method can create time-consuming simulations, in particular for more complex structures consisting of a large amount of vertices and tetrahedra.Furthermore, the Voronoi tessellation which is performed in each iteration, represents a drawback of the method as the computational cost grows with the level of complexity of the analysed structures.As the Voronoi tessellation is enabled within the PyDEC library, which has been used during the study, this option has been employed despite the influence to the computational time.Although the computational cost depends on the initial position of the nodes, i.e. the method converges faster for analyses with smaller displacements, the determination of the position of the dual complex introduces a demanding computational task to the method.However, as the connectivity within the simplicial complex does not change during the loading as compared to the original configuration, the update of the coordinates of the dual 0-cells based on the current positions of the primal 0-cells would considerably expedite the calculation.This step will be included during the future work within the improvement of the computational performance of the method.In Fig. 12 , the convergence of the calculated total energy for the simulation of the loaded beam is shown.Despite an apparent descent of the total energy to the final value it can be seen that a considerable number of iterations were necessary to obtain a suitable solution.The number of iterations directly influenced the computational costs.The duration of the analysis performed by an Intel Core i7-8550U 1.8 GHz processor using the lattice-complex modelling took 3140 s whereas the finite element method took 7 s.A factor of about 450 has therefore been observed for the calculation time in favour of the finite element method.
This implies that future work will investigate other energy minimisation methods that will ensure fast and robust solutions using lattice-complex modelling.Moreover, although the results of the comparisons are encouraging, further studies concentrating on the calculation of volumetric and distortional energies for other structures are required to confirm our method.Currently, the finite element method supported by decades of research on the subject is a far more effective method in terms of the computational time, as also noted by other researchers, e.g.Kumar et al. (2016) .
The formation of both the cyclic bulk modulus-strain and the cyclic shear modulus-strain curves constitutes the base for the simulation of the elastoplastic effect using the lattice-complex modelling.The curves are created from the cyclic stress-strain curve which is usually based on experimental results gained from strain-controlled uniaxial tests.A tensile stress-strain curve can also be used for the determination of the bulk modulus-strain and shear modulus-strain curves.In the latter case however, discrepancies between the simulated and the experimentally observed cyclic response might be noted.As mentioned, the initial value of the moduli is valid until the yield strain.Beyond the yield strain, the bulk modulus starts increasing due to an increase in the Poisson's ratio and the shear modulus starts decreasing due to a decrease in the tangent modulus of the material.This can be seen in Fig. 5 for the titanium aluminide alloy Ti-42Al-8.5Nbused in this study.The bulk modulus represents the resistance of the material against its volume change and the shear modulus stands for the resistance of the material against its shape change.As the ratio between the change of the bulk modulus and the change of the shear modulus are nonlinear as the load of the structure grows beyond the yield point of the material, a nonlinear deformation of the structure is observed as a result when a minimised total energy has been reached.For either the alternating or the variable load history used here, the increase of the bulk modulus and the decrease of the shear modulus can be seen in Figs. 10 and 11 in state Nr. 2. If the load is reversed, the material behaves elastically again until it crosses the whole elastic region and reaches the yield strain in the compressive direction.This happens after twice the value of the yield strain if kinematic hardening of the material is considered as shown here for state Nr. 3 in Fig. 10 and state Nr. 5 in Fig. 11 .Whilst crossing the elastic region, the bulk and the shear moduli regain their initial ('elastic') values which can be seen as jumps just after states Nr. 2 and 3 in Fig. 10 and jumps just after states Nr. 2-7 in Fig. 11 .
Correct values of the moduli in every instant during the cyclic loading are assured by application of the Prandtl operators ( Eqs. ( 12) and ( 17) ) which continuously model the elastoplastic hysteresis loops of the bulk modulus in every dual 3-cell and the elastoplastic hysteresis loops of the shear modulus in every interior segment of the dual 3-cell.The formation of the hysteresis loops, either simple during the alternating loading or nested as in the case of the variable loading history, can be noticed for e.g. the bulk modulus during the variable loading history as shown in Fig. 13 .Although abrupt jumps are noticed in Figs. 10 and 11 for both the bulk and the shear modulus after every reversal point during the simulation, a continuous modelling of the bulk modulus is observed if only the bulk back-strains are shown from Eq. 12 .The abrupt jumps hence follow from the fact that the values of the moduli needed for the simulation originate from both the continuous modelling and their relative positions to the values in reversal points during the loading history.The first attempt to introduce elastoplasticity into geometric modelling is presented in this paper.Nevertheless, further studies concentrating on elastoplastic response of mechanically loaded structures are required to confirm the computation.
The future ambition of lattice-complex modelling is fracture simulation in solids considering elastoplastic material properties in terms of the bulk and shear moduli, geometry representation by discrete topological spaces and the analyses of their behaviour using only geometric features whilst ensuring a preservation of the minimum stored total deformation energy.Geometric discontinuities will be considered by lattice-complex modelling as the characterisation of the connectivity between the primal and dual nodes.

Conclusions
A method for the simulation of mechanically loaded elastic and elastic-plastic solids using a combination of lattice models and elements of the discrete exterior calculus is developed and proposed.It is demonstrated that the separation of the total deformation energy into volumetric and distortional parts enables the correct modelling of discrete structures.Accurate results were initially obtained for the behaviour of a loaded icosahedron over the entire range of Poisson's ratios.Subsequently, good agreement between the model solution and the theoretical solution for a rigidly clamped beam was demonstrated.Finally, simulations of the bulk modulus and shear modulus during cyclic loading indicated that the proposed methodology is applicable to cyclically loaded solids.Lattice-complex representation of solids enables simulation of the entire range of possible material elastic behaviours.Considering the primal intention of lattice models is to allow for studies of failure initiation and crack propagation in solids with explicit microstructure representation, the proposed approach offers a consistent methodology to capture the elastic and plastic behaviour prior to the initiation of such failures.Work on effective and efficient energy minimisation schemes is ongoing.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Representation of a discretised solid in terms of a de Rham space.
volumes of 1-cells (physically meaning lengths), denoted by | σ 1 | and | * σ 2 |, volumes of 2-cells (physically meaning areas), denoted by | σ 2 | and | * σ 1 |, and volumes of 3-cells, denoted by | σ 3 | and | * σ 0 |, do not affect the discrete derivatives, but are essential for completing the discrete versions of the continuum gradient, curl and divergence operations.This is done by the introduction of discrete Hodge star operators, * 0 which are depicted in Fig.1, and where the volumes of 0-cells are assumed to be 1, i.e. | σ 0

Fig. 4 .
Fig. 4. Determination of bulk and shear Prandtl densities from a) the cyclic bulk modulus-strain curve and b) the cyclic shear modulus-strain curve.

Fig. 6 .
Fig. 6.Primal and dual complexes of the icosahedron consisting of 13 primal 0-cells and 20 primal 3-cells.The edge with the applied displacement and the supported edge can be observed in the middle.

Fig. 7 .
Fig. 7.The results of the loaded icosahedron.Equal axial displacement in the z-direction causes different circumferential displacements.The diagram shows comparison between the theoretical Poisson's ratio (dashed line) and the calculated Poisson's ratio using lattice-complex modelling (solid line).

Fig. 8 .
Fig. 8. Primal and dual complexes of the analysed beam consisting of 36 primal 0-cells and 48 primal 3-cells.The support with the applied displacement and the supported side can be observed in the middle.

Fig. 9 .
Fig. 9. Deformation of a loaded beam.The vertical y-displacement of one of the supports creates deformation of the whole structure (upper diagrams).The lower diagram shows a comparison between the vertical displacements of the neutral z-line of the beam using lattice-complex modelling, an analytical solution and the finite element method.

Fig. 10 .
Fig. 10.Simulated time history of a) the bulk modulus and b) shear modulus during alternating mechanical load history (solid grey line); simulated strain path of c) the bulk modulus and d) shear modulus during alternating mechanical load history (solid grey line).Strain histories in a) and b) are depicted with a thin black line.Thin black lines in c) and d) are used to indicate the cyclic bulk modulus-strain curve and the cyclic shear modulus-strain curve, respectively.

Fig. 11 .
Fig. 11.Simulated time history of a) the bulk modulus and b) shear modulus during variable mechanical load history (solid grey line); simulated strain path of c) the bulk modulus and d) shear modulus during variable mechanical load history (solid grey line).Strain histories in a) and b) are depicted with a thin black line.Thin black lines in c) and d) are used to indicate the cyclic bulk modulus-strain curve and the cyclic shear modulus-strain curve, respectively.

Fig. 12 .
Fig. 12. Convergence of the calculated total energy for the simulation of the loaded beam using the quasi-Newton optimisation method.

Fig. 13 .
Fig. 13.Modelling of bulk modulus during a) alternating and b) variable mechanical load history.