Numerical and constitutive modeling of quasi-static and dynamic mechanical behavior in graded additively manufactured lattice structures

ABSTRACT Metallic lattice structures based on triply periodic minimum surfaces (TPMS) have attracted extensive attention for their potential application in lightweight and energy absorption. The underlying phenomena, mechanisms and modelling under the crushing responses from quasi-static to shock conditions still remain to be revealed. This work systematically investigates the mechanical behaviour of graded additively Schoen-F-RD (FRD) lattice structures under various loading rates. Under dynamic compression, FRD lattices exhibit the ability to withstand larger densification strains at higher plateau strengths, thus, holding enhanced energy absorption capabilities. At medium strain rates, it is the rate-dependence of lattice base material dominates in the strength enhancement, while, at higher strain rates, the role of inertia effect becomes notable. Furthermore, an empirical formula is introduced to predict the shock stress responses. Finally, constitutive models with strain-rates are proposed for the uniform and graded lattices. These findings can provide excellent guidance on the design of energy-absorbing structures.


Introduction
Cellular materials have found widespread use in various engineering fields due to their excellent mechanical characteristics, such as lightweight, high specific stiffness/strength, good energy absorption capacity, and multifunctional properties [1][2][3][4][5].Typically, cellular materials consist of an interconnected network of solid struts and/or plates forming the edges and faces of pores and cells [6].Depending to their cell configuration, cellular materials can be categorised as either stochastic or periodic structures [7,8].Stochastic structures, like bones and foams, exhibit random distributions of cell characteristics (size, shape, orientation, strut or wall thickness, etc.) [9].Periodic structures, on the other hand, are defined based on unit cell design parameters, and their mechanical properties can be regulated by controlling the unit cell parameter [10,11], offering much higher designability than stochastic ones.
The advancement of additive manufacturing (AM) technology has facilitated the accurate production of various periodic structures, including body-centred cubic (BCC) lattices [12], Octet truss lattices [13], Kelvin lattices [14], rhombic dodecahedron (RD) structures [15,16], and triply periodic minimal surface (TPMS) structures [17][18][19].TPMS structures consist of non-intersecting surfaces with zero mean curvature at any point on the surface [20], effectively avoiding stress concentration and exhibiting superior mechanical properties compared to some traditional lattice structures [21][22][23][24].Generally, TPMS structures can be classified into skeletal-, strut-and sheet-based structures.Research [25][26][27] has shown that the sheet-based TPMS structures can offer better mechanical performance compared to the skeletal-and strut-based structures.Therefore, this study concentrates on metallic sheet-based TPMS structures, leveraging the strength of metals and the mechanical advantages of sheet-based structures.
Various sheet-based TPMS structures have been developed using level-set equations, including Primitive, Diamond, Gyroid, Schoen-I-WP (IWP) and Schoen-F-RD (FRD).Previous studies have explored the mechanical properties of these sheet-based uniform TPMS structures [8,[28][29][30][31].For instance, Zhang et al. [31] investigated the quasi-static compressive mechanical properties of three uniform TPMS sheet structures (Primitive, Diamond, and Gyroid); and highlighted the superior energy absorption capacity of TPMS sheet structures compared to BCC lattices.Guo et al. [28] noted the influence of opening diameter on the mechanical performance of P-lattice structures.Li et al. [30] and San Ha et al. [29] analysed the dynamic compressive characteristics (deformation modes, stress-strain responses and energy absorption) of Gyroid and Primitive structures, respectively.These studies emphasised that the mechanical properties of uniform TPMS structures are mainly dominated by cell morphology and relative density.
To further enhance the mechanical performance of TPMS structures, the graded structure strategy has been recently introduced to sheet-based TPMS structures.By continuously or stepwise varying sheet thickness or cell size, functionally graded TPMS structures (FG-TPMS) can be tailored to specific application needs [32][33][34].The mechanical performances of FG-TPMS structures under quasi-static loading have been explored by several researchers [35][36][37][38].Al-Ketan et al. [35] investigated the deformation mechanisms and elastic properties of FG Schoen Gyroid and Schwarz Diamond structures.Zhang et al. [38] showed that a continuous gradient design helps eliminate small cracks in the Gyroid structures, thereby enhancing their mechanical performance.Zhao et al. [36] found that the energy absorption of FG Primitive and Gyroid structures under quasi-static compression can be approximately 60% higher than that of uniform counterparts.Up to date, only limited works about the dynamic behaviours of FG-TPMS structures have been reported [39,40].Santiago et al. [40] demonstrated 18% higher impact energy absorption in FG Diamond structures compared to uniform ones.Feng et al. [39] revealed improved blast resistance in the sandwich panel with an FG Gyroid core.While FG-TPMS structures have shown advantages in terms of compressive strength and energy absorption, the development of a constitutive model accounting for strain-rate effects remains a need.
It is well known that loading velocity significantly impacts the deformation modes of metallic cellular materials (MCMs), subsequently affecting their dynamic constitutive behaviours [41][42][43].Under high crushing velocity, the compressive strength and energy absorption enhancement phenomenon were observed in the TPMS structures with the metal base material, such as 316L stainless steel [30,44].The debate continues regarding whether the dynamic stress enhancement in MCMs is primarily due to the base material's rate-dependence or inertial effects [45][46][47][48].Thus, revealing the mechanisms underlying the dynamic crushing responses of metallic TPMS lattice structures is of utmost necessity.
Inspired by the above, this study aims at gaining a more comprehensive understanding on the phenomena, mechanisms and modelling under the crushing responses of uniform and FG-TPMS structures at a wider strain-rates ranging from quasi-static to shock states.Specifically, sheet-based FG FRD structures are designed using a sheet thickness grading method.The numerical models developed are verified against the quasi-static and dynamic experimental results.Afterwards, detailed discussions on quasi-static and dynamic crushing characteristics in terms of stress-strain responses, deformation modes and energy absorption are presented.By decoupling the base material's inertial effect and rate-dependence, this study unveils the dynamic enhancement mechanisms for FG FRD structures.Finally, strain-rate dependent constitutive models are proposed and validated.

Design of TPMS lattices
Based on Yin et al. [8], the FRD sheet structure was found to have the best crashworthiness among the four investigated types of TPMS structures (Primitive, FRD, IWP, and Gyroid), making it the focus of this study.The geometrical function of the FRD sheet structure in threedimensional space can be defined as follows [8]: where x, y, and z are three Cartesian coordinates and l is the cubic unit cell length.Significantly, the size and volume fraction of the pores of FRD structures can be directly affected by the constant C [8].In this study, the level constant C is set as a constant value (zero here).For the uniform FRD lattice structures, the geometry model can be easily obtained by replication in three perpendicular directions.In general, there are mainly two modelling techniques to design the graded lattices [34]: (1) cell size grading (CSG), where the cell size is varied with a constant sheet thickness, and (2) sheet thickness grading (STG), where the sheet thickness is varied with a constant cell size.Previous study [49] has demonstrated that the gradation approaches (CSG and STG) have little effect on the crushing responses of graded metallic foams.Here, the focus is on the graded lattice structures generated using the STG method.
Compared to lattice structure with a smooth gradient, the layer-wise gradient configuration with material properties varying in a layer-by-layer fashion exhibits great simplicity and efficiency in terms of manufacturing and modelling [49].As a result, the graded lattice is in discrete form with the layer-wise gradient in this study.Figure 1 shows two examples of uniform and stepwise graded FRD lattice structures.The designed lattices were generated in shell form, and the thicknesses were assigned in the model based on the relative density.Relative density is a crucial criterion for evaluating the mechanical performance of cellular structures, defined as the ratio between the densities of the cellular structure and the base material [1,50].

Set-up of the FE model
The crushing responses of the uniform and graded FRD lattice structure are numerically investigated using the LS-DYNA/Explicit package, as shown in Figure 2. The FE model of the FRD lattice is sandwiched between two rigid plates, namely the 'moving plate' and 'stationary plate' (see in Figure 2 (a)).The faces of the FRD lattice structure are discretized using Belytschko-Tsay thin shell elements with an element size of 0.3 mm according to a mesh sensitivity analysis (Figure 2 (b)).Contacts between the two plates and the lattice are simulated using the automatic surface-to-surface algorithm, with static and dynamic friction coefficients set at 0.3 and 0.2, respectively [8].In addition, the automatic single surface algorithm is applied to the faces to prevent self-contact penetration.To improve the computational efficiency, the increase in quasi-static testing velocity was used to speed up the simulations and was confirmed as acceptable where the reaction forces at the impact end and support end are balanced during the whole loading process [51,52].Following this method, the 'moving plate' moves downward at a constant velocity of 2 m/ s, while the nodal degree of freedoms (DOFs) of the 'stationary plate' are constrained, as recommended in Refs.[51,53].It needs to be pointed out that the strainrate dependent behaviour of the lattice base material is not considered in the condition of quasi-static crushing simulation [47].In addition, the effect of air on the compressive responses can be ignored due to that the investigated FRD lattices are open-cell with the large cell size [6].
To simulate dynamic crushing conditions, different loading velocities are applied for the 'moving plate'.The strain-rate dependence of the lattice base material is described using the Cowper-Symonds equation [30,54]: in which s d is the flow stress that depends on strain-rate, s ys denotes the static yield stress, 1 is the strain-rate, and C and p are the strain-rate sensitivity parameters.Different relative densities ( r) of FRD lattice structures can be obtained by changing the shell thickness [42] as follows: in which r 0 and r s are the densities of the lattice structure and the base material, respectively; A i and V are the area and volume of the cell-wall, respectively.n is the number of the cell-wall and t is the shell thickness.
Using the STG modelling method, a constant shell thickness is assigned to each layer of the FE model, allowing for the generation of functionally stepwise graded FRD lattice structures.Note that one cell size is regarded as one layer here.The gradient configuration can be adjusted by changing the shell thickness of each layer.The thickness changes between two interconnected layers of the functionally stepwise graded lattice structures are tackled by setting up different shell thicknesses at nodes.

FE model validations
Uniform FRD lattice structures, consisted of 64 (4 × 4 × 4) unit cells with a cell size of 15 mm and the shell thickness of 1 mm, were utilised to validate the numerical accuracy of FE models under quasi-static loading.The lattice structure was made of EOS CX stainless steel and its material properties were determined through standard tensile experiments using SLM fabricated tensile samples, as listed in Table 1.More details about the SLM processing parameters to manufacture lattice structures can be referred to in Ref [8].Generally, the mechanical properties of additively-manufactured lattice structures are also influenced by the manufacturing process [3,9] and printing direction [54][55][56], et al., which will be a focus of our future study.In this work, the main focus is on exploring the strain-rate factor in affecting the crushing characteristics of FRD lattices along the printing direction.
A corresponding FE model featuring a 64-cell (4 × 4 × 4) model is also established, with the detailed set-up available in the above section.Note that unit size analysis has proved that the 64-cell model (4 × 4 × 4) is sufficient to capture the overall mechanical properties of cellular structures [8,14,57].
Comparisons of the crushing mechanical responses (deformation mode and force-displacement curve) of the FRD structure between FE simulations and the experimental results are displayed in Figure 3.It can be seen that the evolution of deformation modes in the FE models during three loading processes (the displacement s = 3.6, 12 and 30 mm) matches well with those of the experimental results (Figure 3 (a)).Specifically, X-shaped shear band is observed both in experimental and numerical lattice specimen when compressive displacement is 30 mm.The force-displacement curves derived from the simulation results generally exhibit good agreement with the experimental results (Figure 3 (b)).The mean crushing force (MCF), defined as a ratio of total energy absorption to the maximum crushing distance, is an essential metric to evaluate the global compressive mechanical response of lattice structures [4,5].The numerical MCF (460.91 kN) is 14.62% smaller than the experimental results (539.81 kN), which is mainly caused by the mass error between asdesigned and manufactured lattice samples [16,31].
To validate the numerical accuracy of FE models under dynamic loading, 64-cell (4 × 4 × 4) uniform lattice with shell thickness of 0.25 mm and 96-cell (4 × 4 × 6) graded specimen with six graded layers and shell thicknesses of 0.15, 0.1625, 0.175, 0.1875, 0.20 and 0.2125 mm were manufactured by SLM process and then used to conduct dynamic compression tests by a drop-tower testing system (Instron 9250 HV) (see Figure 4 (a)).Each cell size of the lattice specimens is 3.75 mm for dynamic tests.Note that the SLM processing parameters adopted are kept the same in fabricating different lattice specimens.In FE models, the total drop weight mass (18.58 Kg) and initial impact velocity  (3.0 m/s) were assigned on the 'moving plate' in line with the experimental conditions.The value of C and p are set as 1704.5/sand 5.2 in C-S model (see Equation ( 2)) to describe the dynamic mechanical properties of stainless steel lattice base materials, as recommended in Refs.[30,58,59].From Figure 4 (b), it is shown that the numerical force-displacement curves show an acceptable agreement with the experimental results.Overall, Figures 3 and 4 indicate that the results predicted by FE models correlate well with those obtained by experiment and thus is considered to be valid for subsequent studies.

Stress-strain responses
Previous studies [30,60] have demonstrated that TPMS lattice structures with the relative density of 20% and 30% possess more excellent mechanical properties, such as higher compressive strength and specific energy absorption (SEA).In general, cellular materials signify one of the porous structures typically with a porosity over 70%, and the relative density is less than 30% [11].To take the lightweight characteristics of lattice structures, seven shell thicknesses of 0.60, 0.65, 0.70, 0.725, 0.75, 0.80 and 0.85 mm, corresponding to the relative densities of 0.194, 0.210, 0.226, 0.234, 0.242, 0.258 and 0.274, are considered to explore the FRD lattice density effect on the compressive characteristics.Note that FRD structures for the quasi-static and dynamic numerical analysis consist of 96-unit cells (4 × 4 × 6), and each cell size is 15 mm.
For the uniform FRD lattice structures with varying densities (shell thicknesses), the numerical nominal stress-nominal strain curves are plotted in Figure 5 (a).Significantly, the stress-strain curve profiles for different densities all exhibit the characteristic threestage behaviour of cellular materials, a linear elastic stage (Stage I), followed by an extended plateau stage (Stage II) and finally a densification stage (Stage III).Notably, a continuous hardening effect is observed during Stage II.At Stage III, stress increases sharply due to more lattice sheets coming into contact with each other, resembling the behaviour of solid materials.
Generally speaking, the plateau stress and densification strain are the two key metrics to evaluate the compressive mechanical performance of cellular materials [1,5,61].Here, an energy efficiency-based approach is utilised to quantitatively investigate how the relative density affects these two evaluating indicators.Energy absorption efficiency h(1) is defined as the energy absorbed up to a given nominal strain 1 normalised by the corresponding stress value s c (1) [62][63][64], as: The densification strain 1 D is defined to be the strain when the energy absorption efficiency h(1) reaches the maximum, mathematically as: The plateau stress s pl is determined by Based on Equations ( 4)-( 6), the calculated plateau stress and densification strain varying with relative density are plotted in Figure 5 (b).Significantly, the plateau stress is gradually increased with increasing relative density.However, there is no evident variation in densification strain in spite of the different relative densities.Specifically, the densification strain variation is less than 2.1% compared with the average value (0.703).In other words, it may be considered that the densification strain is insensitive within the investigated relative density.
Herein, linear gradient configuration involving increasing or decreasing shell thickness are investigated.FRD lattice structures consisting of two, three and six layers, with an average shell thickness of 0.725 mm (relative density of 0.234) are considered.The effect of lattice gradient configuration on compressive stress-strain responses is analysed, as shown in Figure 6.The specimens investigated were labelled characteristically; for example, 'Graded-three-t-0.60-0.725-0.85'indicates a graded specimen with three graded layers and shell thicknesses of 0.60, 0.725 and 0.85 mm along the compressive direction.Note that each layer's height remains constant, equal to the ratio of the total height to the layer number.To facilitate comparison, the stress-strain curve of the uniform specimen with an equivalent density (specimen 'Uniform-t-0.725') is also plotted in Figure 6.Unlike the uniform lattice, lattice specimens with linear density gradient do not exhibit a distinct plateau region.Instead, a phenomenon of multiple plateau stress phases is observed in the stressstrain curves, with the number of phases equal to the pre-designed number of layers, as in [65].For lattice specimens like 'Graded-two-t-0.60-0.85'and 'Gradedtwo-t-0.85-0.60', the obtained stress-strain curves almost overlap (Figure 6 (a)).Similar observations can be found for the specimens 'Graded-three-t-0.60-0.725-0.85'and 'Graded-three-t-0.85-0.725-0.60'(Figure 6 (b)).
Irrespective of the gradation configuration, the plateau stress level generally increases.This implies that graded lattices progressively crush from the weakest layer (lowest density) to the strongest layer (highest density) during quasi-static compression.The following section will discuss the deformation propagation mechanisms to support this viewpoint.
In comparison, the uniform lattice specimen exhibits a stable plateau stress stage.The relative density of the uniform lattice exceeds that of the weakest layers in the graded lattice.Hence, the stress level of graded lattices is relatively lower than that of the uniform specimen in the initial crushing phase.However, as the compressive strain becomes substantial, a higher stress will be generated in the graded lattices relative to the uniform specimen.Moreover, it is evident that the stress levels for uniform and graded specimens become quite close after the densification strain.This is because the lattice behaves like solid materials, as a large number of lattice sheets come into contact with each other during this stage.In this situation, the stress level may be determined by the overall density of the lattice structure.
To generate lattice specimens with approximate continuous gradients, a six-layer graded lattice featuring a linear gradient is designed, and the corresponding numerical stress-strain responses are plotted in Figure 7.It is evident that as the number of layers increases from two or three to six, the stress-strain becomes progressively smoother.A significant continuously linear hardening behaviour is observed in the plateau stage, which may be explained by the successive deformation mode of graded lattices.The stress-strain profile of the six-layer graded lattice is rather close to that of the lattice with continuous gradients [36,39].

Deformation modes
The deformation evolution mechanisms of uniform and graded lattices are detailed and analysed to understand the stress-strain responses better.Without loss of generality, the deformation procedure of three graded lattices and the uniform specimen is presented in Figure 8.In the case of the uniform FRD lattice, relatively uniform deformation is seen at a low compressive strain (0.20); and then V-shaped shear bands can be observed at a compressive strain of 0.40.Subsequently, these shear bands extend along the lattice shells, ultimately becoming fully compressed at higher strains (see in Figure 8 (a)).For the lattice specimen 'Graded-two-t-0.60-0.85',initially, X-shaped shear bands emerge in the weakest layer (lower density).Densification then appears at the weakest layer, while the strongest layer (higher density) experiences minimal deformation (Figure 8 (b)).As the compression procedure continues, the strongest layer replicates the deformation pattern observed in the weakest layer.Such deformation evolutions elucidate the presence of two plateau stress phases in the stress-strain curves (Figure 6 (a)).Similar deformation patterns are also observed in the lattice specimen 'Graded-three-t-0.85-0.725-0.60'.Moreover, it can be concluded that collapse bands consistently initiate in the weakest layer and progressively propagate towards the strongest layer, regardless of the gradation configuration (Figure 8 (c)).Further, V-shaped shear bands are observed on the long lattice specimen (Figure 8 (a)) while X-shaped shear bands on the short specimen (see in Figure 3).Thus, it can be concluded that the overall deformation characteristics of the lattice structure is related to the length of the specimen.
In contrast to uniform and low layer-number graded lattices with significant V-shaped and X-shaped shear bands, the six-layer graded lattice exhibits relatively different deformation modes, as shown in Figure 8 (d).
Beginning with collapse bands in the top layer possessing the lowest density, lattice specimens progressively collapse in an approximate layer-by-layer manner.Such successive deformation behaviours may explain the gradual enhancement of plateau stress with increasing nominal strain (see in Figure 7).A Similar phenomenon has been widely reported in Refs.[34,36,39] for the TPMS lattice structures featuring continuous gradients.

Energy absorption
Energy absorption refers to the ability of a material or structure to withstand external loads without catastrophic failure [1].FRD lattice structure possesses a prolonged, flat plateau stress stage under quasi-static compression (Figure 5 (a)), which emulating an ideal energy absorber.The energy absorption capability, defined as the energy absorbed per unit volume up to densification point (W vt ), has become a critical metric for gauging the energy-absorbing applications of lattice structures, and it can be calculated as follows [36]: where s(1) is the nominal stress (s) related to nominal strain (1) and 1 d is the densification strain.
In addition to W vt , specific energy absorption (SEA), defined as the amount of energy absorbed per unit mass (typical units: J/g), is usually utilised to enable comparison across different energy absorber materials, which is given as follows [1]: where r 0 is the density of the lattice structure.
According to above-defined metrics, the calculated cumulative energy absorption per unit volume versus strain (W v − 1) curves and the SEA of uniform lattice structures are illustrated in Figure 9.The W v − 1 curves are also fitted using the power law equation, as presented in Figure 9 (a).The fitted exponent value is increased with increasing shell thickness (lattice density).Specifically, the fitted exponent values of specimens 'Uniform-t-0.60'and 'Uniform-t-0.85'are 1.08 and 1.12, respectively.Note that the correlation coefficients (R 2 ) for the fitted equations are both over 0.99.In general, a higher fitted exponent value corresponds to a more favourable energy absorption capability [36].Specific energy absorption (SEA) of uniform lattice specimens is also plotted in Figure 9 (b).Evidently, SEA tends to increase as the lattice density rises, exemplified by values of 33.47 and 37.61 J/g for shell thickness of 0.60 and 0.85 mm, respectively.SEA demonstrates a 12.37% increase for a shell thickness of 0.85 mm compared to 0.60 mm.As discussed, it can be concluded that the energy absorption capability of the FRD lattice can be improved by increasing the lattice density.
In addition, the SEA of the investigated uniform lattices is beyond 33 J/g, which is larger than that of certain TPMS structures, including 12.78 J/g of the Primitive lattice [8], 21.64 J/g for the Diamond lattice [40], 21.41 J/g for the IWP lattice [8], and 27.54 J/g for the Gyroid structure [30].At quasi-static strain rates, key factors affecting lattice structures' energy absorption capacity encompass lattice base material properties, cell geometry, and the overall relative density [1].More advanced multiscale optimisation methods [66][67][68] hold promise for enhancing the energy absorption capacity of lattice structures, which can be left as a topic for future investigations.
The stress-strain curves of various graded lattice specimens have been achieved in the previous Section 3.1.1;then the densification strain can be obtained based on Equations ( 4) and (5).Note that the densification strain values for lattice specimens 'Graded-two-0.60-0.85'and 'Graded-two-0.65-0.80'are 0.758 and 0.742, which are larger than the values for the uniform specimen with the same density (0.694).The cumulative energy absorption per unit volume versus strain (W v − 1) curves and the energy absorption capability (W vt ) of uniform and graded lattice specimens are illustrated in Figure 10.From W v − 1 curves, it becomes apparent that the energy absorption performance of the uniform lattice slightly outperforms that of graded specimens up to a compressive strain of 0.70 (Figure 10 (a)).Notice that the relative density of graded lattices is the same as that of uniform specimens.
As for the total energy absorption capacity up to the densification point (W vt ), the difference between the uniform lattice and graded specimens is small (Figure 10 (b)).Specifically, the maximum difference is only 5.74%.In other words, it is hard to effectively enhance the energy absorption performance of FRD lattices  through changing the gradation configuration under quasi-static compression.

Constitutive model
From the nominal stress-strain curves of FRD lattices (see in Figure 5 (a)), it becomes apparent that the elastic deformation is quite small (elastic strain less than 1%), which can be ignored when focusing on energy absorption behaviours at large plastic deformation [52].Through carefully examining the variation trends of stress-strain curves, the rate-independent, rigid-plastic hardening (R-PH) idealisation model, as proposed by Zheng et al. [43], is used to characterise the quasi-static stress-strain relationships, written as: where subscript n denotes nominal, s n0 is the nominal initial crushing stress and C the strain hardening parameter.For each lattice specimen at a particular relative density, the two material parameters, s n0 and C, can be obtained through the least squares fitting technique.The theoretical R-PH model is utilised to fit the nominal stress-strain data of the uniform lattice specimens, as shown in Figure 11.It is evident that the R-PH idealisation model effectively captures the stress-strain relationships.
The two material parameters of the R-PH model, s n0 and C, varying with the relative density of lattice structures are plotted in Figure 12.Based upon Yang et al. [69], s n0 and C in the R-PH model can be usually expressed as the power-law function of relative density (r 0 /r s ), given as follows: where s ys is the yield strength of base material (600.69MPa here), k 1 , k 2 , n 1 and n 2 are material parameters.
Similarly, the material parameters can be determined by applying the least squares fitting technique, specifically k 1 = 1.004, n 1 = 1.353, k 2 = 0.0556 and n 2 = 1.714.The correlation coefficients (R 2 ) for the fitted equations are over 0.99, meaning that the power-law equations have relatively high fitting precision (see in Figure 12).Substituting Equation ( 10) into Equation ( 9), a unified constitutive model considering different relative densities is thus obtained to characterise the quasi-static stress-strain responses, mathematically written as:

Deformation modes
It is well-known that loading velocity has a significant influence on the deformation modes of cellular materials, thereby affecting their dynamic constitutive behaviour  [ 29,42,43].To investigate the dynamic compressive mechanical characteristics of FRD lattice structures, simulations were conducted with different constant loading velocities ranging from 1.8 to 270 m/s, corresponding to nominal strain-rates ranging from 20 to 3000 s −1 , respectively.Note that the nominal strain-rate is defined as the ratio of the impact velocity to the original lattice height [70].
Taking lattice specimens 'Uniform-t-0.75'and 'Graded-six-t-0.85-0.80-0.75-0.70-0.65-0.60'for illustrations, 3D spatial deformation distributions at three different loading velocities (1.8, 135 and 225 m/s) are compared in Figures 13 and 14.Note that the effective plastic strain greater than 0.15 is shown here.It is apparent that at a loading velocity of 1.8 m/s, relatively uniform deformation is observed at a compressive strain of 0.20 (Figure 13 (a)), resembling the quasistatic results.This deformation mode is usually termed the Homogenous Mode.At the high impact velocity of 225 m/s, the lattice specimens collapse layer-by-layer, forming a localised deformation band near the impact end (Figure 13 (c)), which is termed the Shock Mode.
When the impact velocity is applied among the above two scenarios, the deformation band tends to concentrate relatively towards the impact end, known as the Transition Mode.Significantly, the deformation modes of FRD lattice structures are crushing velocity dependent.
From Figure 14 (a), it is evident that at a crushing velocity of 1.8 m/s, layer densification initiates first in the layer with the lowest density, consistent with the quasi-static results.When the crushing velocity is increased to 135 m/s, large plastic deformation occurs nearby the impact and support end (Figure 14 (b)), which behaviours like 'I' shaped deformation bands.At a crushing velocity of 225 m/s, a significant localised deformation band forms in the layer with the highest density (the layer nearby the impact end) (see Figure 14 (c)), indicating the occurrence of Shock Mode.Evidently, the shock deformation band concentrates at the impact end, regardless of the density gradient configurations.

Stress-strain responses
The extracted stress-strain curves of lattice specimens 'Uniform-t-0.75'and 'Graded-six-t-0.60-0.65-0.70-0.75-0.80-0.85'at different strain-rates are plotted in Figures 15 and 16, respectively.For comparison, the quasistatic stress-strain curves are also presented here.It is clear that the compressive strength is gradually enhanced with the increase in strain-rate.At a nominal strain-rate of 2000 s −1 , the stress at the impact end is much higher than that at the support end.This is usually considered as the shock enhancement phenomenon, which could be influenced by macro inertia effects [47,70].
Based on the above energy-absorption efficiency method, the plateau stress and densification strain of the lattice specimens 'Uniform-t-0.75'and 'Graded-six-  t-0.60-0.65-0.70-0.75-0.80-0.85'varying with the nominal strain-rates can be obtained, as shown in Figure 17.Obviously, the plateau stress is increased with increasing nominal strain-rate.For the lattice specimen 'Uniform-t-0.75',the plateau stress at the strain-rates of 20 and 1000 s −1 is 158.83 and 192.53 MPa, which are approximately 70% and 106% higher than the quasi-static results (93.25 MPa).It is seen that the variation of the dynamic densification strain is slight at the strain-rate varying from 20 to 1000 s −1 .The average dynamic densification strain of specimen 'Uniform-t-0.75' is 0.782 within the investigated strain-rates, which is larger than the quasi-static results (0.70), similar to Ref. [30] for the Gyroid structures.Similar results can be obtained for the lattice 'Graded-six-t-0.60-0.65-0.70-0.75-0.80-0.85'and other specimens.In summary, FRD lattice under dynamic compression can hold larger densification strain at higher plateau strength, resulting in improved energy absorption capability.Consequently, the FRD lattice structure can be adopted as the effective energy-absorbing elements especially under high-speed impact applications.The dynamic enhancement mechanisms of FRD lattice structures will be detailed and discussed in the following section.
As discussed above, the deformation modes of lattice structures are loading velocity-dependency.If the loading velocity is high enough, the stress extracted at the impact will be larger than that at the support end.In this scenario, the deformation, stress and strain states become inhomogeneous [30,47].To characterise the uniformity of cellular materials under dynamic loadings, Liu et al. [71] proposed a uniformity coefficient w to determine the critical velocity for mode transition from the Homogeneous Mode to the Transition Mode.The uniformity coefficient w is defined as where s s pl and s i pl are the plateau stress on the support end and impact end, respectively.
Based on Liu et al. [71], the critical velocity is determined when the uniformity coefficient w is reduced to  90%.Taking the uniform and graded lattice specimens 'Uniform-t-0.75'and 'Graded-six-t-0.60-0.65-0.70-0.75-0.80-0.85'for illustrations, the uniformity coefficient w varying with the impact velocity are plotted in Figure 18.Taking the uniformity coefficient of 90% as a benchmark, the critical velocities for the uniform and graded lattices can be numerically obtained by using the dichotomy 125 and 90 m/s, respectively.Following this method, the critical velocities for the uniform lattice 'Uniform-t-0.65'and 'Uniform-t-0.85'can be also determined 120 and 132 m/s, respectively.It is thus concluded that the critical velocity is gradually increased with increasing relative density.

Energy absorption
Based on the above critical velocity analysis, the uniformity coefficient w is over 0.96 under the crushing velocity of 45 m/s (strain-rate of 500 s −1 ).In this case, the extracted forces are approximately balanced at the impact and support end.Taking this crush velocity as an example, the dynamic stress-stain curves and energy absorption of uniform and graded lattice specimens are presented in Figure 19.It is seen that the average compressive strength of graded lattices is smaller than that of uniform specimen at the early plateau stage and then this trend of change will reverse at large compressive strain (Figure 19 (a)).This is because the collapse deformation is initially starting from the weakest layer (the layer with the lowest density) for the graded lattice specimen, as similar to the quasi-static results.With the compressive plastic deformation proceeding, the role of the stronger layer (higher density) will become significant.Thus, the main reason for the stress-stain curve difference between uniform and graded lattices is the deformation modes.
According to the above-mentioned energy efficiencybased approach, the densification strain of uniform and graded lattice specimens can be determined.And then, the energy absorption (W vt ) for different lattice specimens is calculated, as plotted in Figure 19 (b).It is found that the energy absorption capability of the investigated graded lattices is close to that of uniform specimen.Specifically, the energy absorption of the graded lattice specimen 'Graded-six-t-0.60-0.65-0.70-0.75-0.80-0.85' is 146.63 MJ/m 3 , which is 5.40% higher than that of uniform specimen 'Uniform-t-0.725'(139.12MJ/m 3 ).However, the energy absorption of the investigated other graded lattices is smaller than that of the uniform specimen.Hence, it can be concluded that the energy absorption performance cannot be effectively enhanced through changing the gradation configuration under medium strain-rates.
Another typical application scenario for the metallic lattice structures is adopted as the efficient energy-absorbing elements in the explosion-proof field [34,39].Generally, the loading rate induced by blast loading could be more than 200 m/s [49].The loading rate is adopted as 270 m/ s here.The dynamic stress-stain responses of uniform and graded lattice specimens extracted at the impact end and support end are plotted in Figure 20.It is seen that the initial compressive stress is mainly determined by the first lattice layer at the early compressive stage, as shown in Figure 20 (a).As is well known, the lattice with a thinner sheet thickness would possess a lower compressive strength.Thus, the graded specimens with the same layer thickness near the impact end yield almost the same compressive strength up to the strain of 0.20.Moreover, the oscillation of compression stress is remarkable for the graded specimens with the positive gradient configurations (the sheet thickness increased gradually from the impact end to the support end).From the stress-strain curves extracted at the support end (Figure 20 (b)), it is found that the graded lattices with the last layer adopting the thinnest sheet thickness produce the lowest compressive stress at the support end.Specifically, the compress strength is about 141.81 MPa at the relative high strain of 0.60 for the graded specimens 'Graded-two-0.85-0.60'and 'Gradedthree-0.85-0.725-0.60',which is 22.73% lower than the uniform specimen 'Uniform-t-0.725'(185.53MPa).Thus, it is indicated that at a relatively large plastic deformation, the graded lattices with negative gradients yield a smaller compressive strength at the support end relative to the uniform and graded specimens with positive gradients.In this case, the graded lattices with negative gradient exhibit some advantage in the explosion-proof fields, such as being the core in sandwich structures [34,39].However, once the plastic deformation reaches the densification stage, the negative gradation could not minimise the compressive strength at the support end.
The cumulative energy absorption versus compressive strain for the investigated lattice specimens at the crushing velocity of 270 m/s are shown in Figure 21.It can be seen that the graded lattices with negative gradient could absorb more shock energy than the uniform and graded specimens with positive gradient before the compressive strain of 0.60.And then, the benefit of the graded lattices with negative gradation is gradually diminished with increasing compressive strain.The energy absorption performance of graded lattices with positive gradation exhibits a certain advantage over the lattices with negative gradation at the larger stain, such as 1 = 0.80.As discussed above, plastic deformation and energy absorption capacity are two critical factors in choosing the proper gradient of graded lattice structures to be applied in the impact/blast protection fields.

Dynamic enhancement mechanisms
As discussed above, FRD lattice structures exhibit significant stress enhancement phenomenon at dynamic compressive loadings with respect to the quasi-static results.Previous studies have pointed out that the strain-rate  sensitivity of base materials and inertia effect are the dominated factors to cause stress enhancement of metallic cellular materials [47,70].The numerical simulation approach has proven a convenient way to decouple these two effects [42].
Taking the uniform lattice specimen 'Uniform-t-0.75'as an example, the dynamic plateau strength normalised by the quasi-static data as well as the lattice base ratedependence itself varying with nominal strain-rate are presented in Figure 22.For the lattice base material without strain-rate dependence, the value of the normalised s pl is almost kept 1 at the medium strain-rates varies from 20 to 500 s −1 , indicating negligible inertia effect within this strain-rate range.However, the normalised s pl of the lattice base with strain-rate dependence is gradually improved with increasing nominal strainrate.This reveals that the strain-rate dependence of lattice base material dominates in the plateau strength enhancement under medium strain-rate conditions (20 to 500 s −1 ).Moreover, the lattice with the base material rate-dependence exhibits a more significant strain-rate sensitivity when compared with the lattice base material itself (Cowper-Symonds equation).
Under high strain-rate conditions, the deformation mode is transformed from Transition Mode to Shock Mode, and the stress is greatly improved compared to the quasi-static result.The compressive nominal stressstrain curves of the uniform lattice specimen 'Uniformt-0.75'under the loading velocity of 180 m/s are plotted in Figure 23.Note that the plateau stresses are also calculated based on the energy-absorption efficiency method and then plotted here.When the strain-rate is enhanced to 2000 s −1 (180 m/s), the plateau strength is 241.21 and 161.87 MPa with and without considering the lattice base rate-dependence, which is 2.59 and 1.74 times higher than the quasi-static data (93.25 MPa) respectively.Consequently, the inertia effect becomes significant and can not be ignored under high-velocity loadings (see in Figure 22).Generally, the difference in the stress-strain relations for cellular materials results from the different deformation modes [43].The main cause is the layer-bylayer collapse deformation and strain localisation under high strain-rates, which is fairly different from the results under medium loading rates [30].In this course, and the stress enhancement is a joint effect of the lattice base strain-rate dependence and macro inertia.
The shock wave will be generated and propagated from the impact end to the support end when the lattice structures are subjected to a high impact velocity [30].A typical one-dimensional shock wave model with the basic material parameters is illustrated in Figure 24.Based on the continuum-based stress wave theory [43,72], the conservation relations of mass and momentum across the shock front are obtained as follows and where Ḟ(t) is the shock-front speed, {n A , 1 A , s A } and {n B , 1 B , s B } are the physical quantities (velocity, strain and stress) ahead of the shock front and those behind the shock front.Combining Equations ( 13) and ( 14), the stress behind the shock front can be obtained: Based on a rigid-perfectly plastic-locking (R-PP-L) idealisation shock model proposed by Reid and Peng [73], the material ahead of the shock wavefront was assumed to be at rest, with the deformation 1 A = 0, particle velocity n A = 0 and a constant plateau stress s A = s 0 ; while the material behind the shock wavefront was assumed to be thoroughly compacted, with the densification strain (locking strain 1 L = 1 D ) 1 B = 1 D , particle velocityn B = n and the plateau stress level s B = s d p .Note that the particle velocity behind the shock wavefront region is approximated as the impact velocity v. Thus, the dynamic plateau stress at the loading end can be expressed in the following form [73]: Taking the uniform lattice specimen 'Uniform-t-0.75'as an instance, the plateau stress and densification strain are 93.25 MPa and 0.70 at quasi-static states.Based on Equation ( 16), the determined theoretical dynamic plateau stresses at the impact end are 178.28MPa, 226.10 MPa and 284.56 MPa for the three crushing velocities of 180, 225 and 270 m/s, respectively.
The extracted stress-stain curves of lattice specimen 'Uniform-t-0.75'under three different crushing velocities of 180, 225 and 270 m/s are illustrated in Figure 25.Utilising the energy-absorption efficiency method described above, this study successfully determined the dynamic plateau stresses at the impact end to be 241.21,272.76, and 333.63 MPa (Figure 25 (a)).Obviously, the numerical data are higher than the theoretical results.The relative errors between the theoretical predictions and actual numerical results are 26.09%,17.11% and 14.71%, respectively.This indicates that the quasi-static properties only provide a rough prediction of the dynamic shock stress.The essential reason is that the strain-rate dependence of the lattice base material is not considered in quasi-static results.In fact, the dynamic stress and strain states (plateau tress and densification strain) of the FRD lattice are pretty different from those in the quasi-static states (Figures 15 and 16).From Figure 25 (b), it can be seen that the plateau stress at the support end is almost the same for the three considered crushing velocities, indicating that the plateau stress at the support end is almost independent of the crushing velocity.The average the dynamic plateau strength and densification strain at the support end is 167.5 MPa and 0.77.Note that the dynamic material parameters are achieved by considering the strain-rate dependence of the lattice base material.If the shock plateau strength is calculated based on the dynamic properties, the theoretical Equation ( 16) can be thus re-written in the following form: where s d 0 and 1 d D are the dynamic plateau strength and densification strain at the support end, respectively.
According to Equation ( 17), the theoretical dynamic plateau stresses at the impact end have been calculated for the three investigated crushing velocities (180, 225, and 270 m/s) as follows: 244.80, 288.28, and 341.42 MPa.The relative errors between the theoretical predictions and numerical results are 1.49%, 5.69% and 2.33%, respectively.This means the proposed empirical model based on the shock wave theory yields close predictions to the numerical results.Thus, the dynamic stress-strain states need to be applied to predict the stress-strain responses at the shock velocities accurately.

Constitutive model with strain-rates
From Figures 15 and 16, it can be seen that the profile of the dynamic stress-strain curves under high strain-rates is quite different from that of results under medium strain-rates.The extracted stresses at the impact and support end are not equal when Shock Mode occurs under high strain-rates.While under medium strainrates (20 to 500 s −1 ), these extracted stresses are almost equal, indicating that the deformation modes are the Homogenous Mode within the investigated    strain-rate range.This section aims at proposing the dynamic constitutive models so as to provide guidance on the energy-absorbing applications of lattice structures under medium strain-rate conditions, such as crash absorbers in the automotive industry [74].
To accurately characterise the dynamic stress-strain relationships of metallic foams, Zheng et al. [43] proposed a dynamic rigid-plastic hardening idealisation model to describe the plastic hardening behaviour given by.
where s d 0 is the dynamic initial crush stress and D is a fitting material parameter.For convenience, this model is donated as the D-R-PH idealisation here.
Taking the uniform lattice specimen 'Uniform-t-0.75'as an example, the R-PH and D-R-PH idealisation models are utilised to fit the numerical nominal stressstrain relationships at medium strain-rates, as plotted in Figure 26 (a).The correlation coefficients for the fitted R-PH and D-R-PH idealisation model are both over 0.93, which means that the dynamic stress-strain responses can be well depicted by using the D-R-PH model.Suppose the dynamic material parameters (s d 0 and D) can be well expressed as functions of the quasistatic parameters (s n0 and C 0 ) and nominal strain-rate (1).In that case, the theoretical constitutive model with strain-rates will be obtained.In this regard, the constitutive model can be given: Further, the normalised initial crush strength (the ratio between the dynamic and static stresses) can be well described by [75]: where a and b are the fitting material parameters, respectively.Based on Wang et al. [63], the normalised hardening parameter can be properly described in the following equation: where c and d are the fitting material parameters.Note that s n0 and C 0 are the fitting material parameters at quasi-static.The normalised dynamic material parameters versus the investigated strain-rate are plotted in Figure 26 (b).The correlation coefficients for the fitted equations (Equations ( 20) and ( 21)) are both over 0.99.This indicates that the above two defined equations can well characterise the relationships between the normalised material parameters (s d n0 (1)/s n0 , D(1)/C 0 ) and the strain-rate (1).Specifically, the fitting parameters of a, b, c and d are 0.5025, 0.1206, 0.8717 and 0.0343 for the investigated uniform FRD lattice specimen here.
Substituting Equations ( 20) and ( 21) into Equation (19), the dynamic constitutive model under medium strain-rates can be thus obtained: As for the graded lattice specimen 'Graded-six-t-0.60-0.65-0.70-0.75-0.80-0.85',linear hardening and continuous hardening phenomenon are observed in the stress plateau and densification stage, respectively, which is quite different from the hardening behaviour for the uniform specimens.To take this into account, we propose a dynamic rigid-linear-plastic hardening model to characterise the stress-strain responses of graded lattices under medium strain-rates, viz.
where E d P and D are the hardening modulus at the plateau stage and strain hardening parameter, respectively; s d n0 and s d nD are the dynamic stresses at the initial crush and densification point; and their relationship is s d nD = s d n0 + E d P 1 D .Note that the elastic deformation is also ignored here.For convenience, this model will be referred to as the D-R-LPH idealisation here.
The D-R-LPH idealisation models are utilised to fit the numerical nominal stress-strain relationships at medium strain-rates, as plotted in Figure 27 (a).Significantly, the proposed D-R-LPH idealisation model can well reproduce the dynamic stress-strain response at different strain-rates.Following the above method, the dynamic material parameters s d 0 , E P and D are also normalised using the quasi-static results.The normalised dynamic material parameters versus with the investigated strain-rate are plotted in Figure 27 (b).Through trial and comparison, the relationships of the normalised parameters and strain-rate (1) can be properly described by:  From Figure 27 (b), it can be found that the above three defined equations can well characterise the relationships between the normalised material parameters and the strain-rate.Specifically, the fitting parameters of a, b, c, d, e and f are 0.5395, 0.1047, 0.3569, 0.1708, 1.1766 and 0.0240 for the investigated graded FRD lattice here.
Substituting Equations ( 24)-( 26) into Equation ( 23), a unified expression is thus obtained to describe the dynamic stress-strain responses for the graded lattices, viz.
where s d D = s n0 (1+a(1) b )+E P0 (1+c(1) d )1 D ; and 1 D can be selected as the average densification strain at medium strain-rates.
Further, the numerical simulations for the crushing velocity of 27 m/s (1 = 300 s −1 ) are also performed to validate the accuracy of the proposed theoretical constitutive models with strain-rates (see Equations ( 22) and ( 27)).The numerical and theoretical results are compared in Figure 28.Obviously, the numerical results could be well captured by the theoretical results, indicating that the proposed theoretical constitutive models are of high accuracy in predicting the dynamic stressstrain responses of uniform and graded lattices.

Conclusions
In this study, the crushing characteristics of uniform and graded TPMS lattice structures at a broader range of strain-rates (quasi-static to 3000 s −1 ) are investigated by using finite element methodology.Stress-strain responses, deformation modes, energy absorption and constitutive model of uniform and graded lattices are discussed and analysed in detail.The main conclusions can be drawn as follows: Under quasi-static compression, multiple plateau stress phase phenomenon is observed on the stressstrain curves of graded lattices and these curves become smoother with increasing graded layer number.The collapse bands of graded lattices are formed from the weakest layer to the strongest layer regardless of the gradation configurations.It is found that the energy absorption capacity cannot be effectively enhanced through changing the gradation configurations.A constitutive model is built for predicting the quasi-static stress-strain responses.
With regard to dynamic compression, the shock deformation band is found concentrated at the impact end regardless of gradient types.The critical velocity for the mode transition from Homogeneous Mode to Transition Mode gradually increases with increasing lattice density.Larger densification strain at higher plateau strength is observed on the dynamic stressstrain curves relative to the quasi-static results.The design of graded lattices is suggested to consider plastic deformation and energy absorption.Numerical results reveal that the strain-rate dependence of lattice base material is the main cause of the strength enhancement under medium strain-rate, while the role of the inertia effect is notable at higher strain-rates.An empirical model is introduced based on the dynamic stressstrain states that can accurately predict the shock stress responses.And finally, the constitutive models involving strain rate-dependence are proposed for the uniform and graded lattices.It needs be mentioned that the constitutive modelling approach is not limited to the FRD lattice and can be potentially applied to other cellular materials.

Figure 2 .
Figure 2. Numerical set-up of the FRD lattice under crushing loadings: (a) the established FE model and (b) mesh size.

Figure 4 .
Figure 4. (a) Drop-tower testing system; and (b) Compassions of force-displacement curve between the simulation and experiment.

Figure 5 .
Figure 5. (a) the numerical stress-strain curves and (b) the plateau stress and densification strain varying with the relative density.

Figure 6 .
Figure 6.The effect of lattice gradation on the stress-strain responses: (a) two-layer graded configurations and (b) three-layer graded configurations.

Figure 10 .
Figure 10.Comparisons of uniform and graded FRD lattices under quasi-static compression: (a) W v − 1 curves; and (b) energy absorption capacity (W vt ).

Figure 11 .
Figure 11.The nominal stress-strain data of the uniform lattice specimens fitted by the R-PH idealisation model.

Figure 12 .
Figure 12. s n0 and C varying with the relative density of lattice structures.Note that the fitted power-law equations are also plotted here.

Figure 13 .
Figure 13.3D spatial deformation distributions for lattice specimen 'Uniform-t-0.75'at three different loading velocities: (a) 1.8 m/s; (b) 135 m/s; and (c) 225 m/s.Note that the compressive strain is 0.20 here.

Figure 15 .
Figure 15.Stress-strain curves of lattice specimen 'Uniform-t-0.75'at different strain-rates for two cases: (a) extracted at impact end and (b) extracted at support end.

Figure 19 .
Figure 19.The dynamic stress-stain responses (a) and energy absorption (b) of uniform and graded lattice specimens under the crushing velocity of 45 m/s.

Figure 20 .
Figure 20.The dynamic stress-stain responses of uniform and graded lattice specimens under the crushing velocity of 270 m/s: (a) impact end; and (b) support end.

Figure 21 .
Figure 21.The energy absorption capacity versus compressive strain for the uniform and graded lattice specimens.

Figure 22 .
Figure 22.Strain-rate sensitivity of the plateau stress for the lattice 'Uniform-t-0.75'with/without base rate-dependence in comparison with the lattice base rate-dependence itself.

Figure 24 .
Figure 24.Shock wave model with the basic material parameters.

Figure 26 .
Figure 26.Constitutive model characterization of the lattice specimen 'Uniform-t-0.75':(a) stress-strain curves fitted by the R-PH and D-R-PH model; and (b) the normalised material parameters varying with the investigated strain-rate.

Figure 25 .
Figure 25.Stress-strain curves of lattice specimen 'Uniform-t-0.75'under three different crushing velocities of 180, 225 and 270 m/s for two cases: (a) extracted at impact end and (b) extracted at support end.
e + f ln(1)(26) where s n0 , E P0 and D 0 are the quasi-static fitting material parameters using Equation (23); a, b, c, d, e and f are the fitting parameters.