Materials such as elastomeric foams, lattices, and cellular solids are capa- ble of undergoing large elastic volume changes. Although many hyperelastic constitutive formulations have been proposed for deviatoric (shape changing) behaviour, few variations exist for large deformation volumetric behaviour. The rst section of this paper presents a critical analysis of current volumet- ric hyperelastic models and highlights their limitations for large volumetric strains. In the second section of the paper we propose three novel volumetric strain energy density functions, which: 1) are valid for large volumetric defor- mations, 2) o er separate control of the volumetric strain sti ening behaviour during shrinkage (volume reduction) and expansion (volume increase), and 3) provide precise control of non-monotonic volumetric strain sti ening. To illustrate the ability of the novel formulations to capture complex volumet- ric material behaviour they are tted and compared to a range of published experimental data. Keywords: Hyperelasticity, Volumetric Deformation, Finite Strain, Strain Energy Density, Shrinkage, Expansion Corresponding author Corresponding author Email addresses: kevin.moerman@nuigalway.ie (Kevin M. Moerman ), patrick.mcgarry@nuigalway.ie (Patrick McGarry ) Preprint submitted to - December 9, 2019 1. Introduction 1 Foams, lattices, and cellular materials are common in nature and engi- 2 neering applications [1, 2, 3, 4]. Given the large elastic volume changes these 3 materials can undergo [5], accurate descriptions of material behaviour be- 4 yond the small strain domain is required. Such behaviour may include an 5 asymmetric shrinkage-expansion response in addition to a highly non-linear 6 pressure-volume relationship. 7 Hyperelastic continuum models o er a convenient means to model the 8 large strain mechanical behaviour of complex materials. However, although 9 the hyperelasticity literature is rich in terms of variations in modelling of 10 the deviatoric (shape changing) material response, as evident from the many 11 di erent formulations which have been proposed (e.g. close to a hundred 12 described in a recent review [6]), few variations exist for modelling the volu- 13 metric contribution (e.g. [7, 8, 9, 10]). 14 The goal of this study is to provide novel volumetric strain energy formu- 15 lations which: 1) are valid for large volumetric deformations, and 2) enable 16 the separate control of volumetric strain-dependent sti ening during shrink- 17 age (volume reduction), and expansion (volume increase). 18 Accurate modelling of the volumetric behaviour of materials like soft 19 foams and lattices is of interest since it is relevant to the design of support 20 structures (see review [11]), such as foam cushions and seats [12, 13, 14, 15], 21 helmets [16], and shoes and insoles ([17, 18, 19]). Furthermore compliant lat- 22 tices and foams are also employed in the design of soft robotics [20, 21, 22, 23]. 23 Recent advances in material science include the development of ceramic 24 nanolattices [24], mycelium-based bio-foams [25], ultraporous sponges [26] 25 graphene foams and aerogels (e.g [27, 28, 29, 30]) some capable of recovering 26 from 90% compression [31]. Furthermore, accurate volumetric formulations 27 are relevant to stroke biomechanics research since blood clot contractions 28 cause large volume changes (e.g. possibly over 80% volume loss [32]). 29 Hyperelastic constitutive formulations have their origins in the modelling 30 of rubber [33, 34, 35, 36]. Although rubber is most commonly modelled as in- 31 compressible (no volume change), it does present with a non-linear pressure- 32 volume response during large deformation hydrostatic compression (20% vol- 33 ume reduction [36]), and volumetric hyperelastic formulations capturing this 34 behaviour have been proposed (e.g. [7]). These, as we will show here, are 35 however not generally valid for very large hydrostatic compression. 36 Hyperelasticity is also commonly used for soft tissues (see review [6]). 37 2 However, like with rubbers, these are often assumed to be incompressible or 38 nearly-incompressible. For these applications the volumetric contributions 39 are commonly considered largely as a simple penalty term to enforce (near) 40 incompressibility (e.g. [37, 38]), rather than a topic of detailed investigation. 41 Consequently the formulations used for these volumetric contributions are 42 often not valid for very large volume changes. For foams and highly com- 43 pressible materials the so called hyperfoam formulation [8, 39, 40]) is common 44 (see for instance [41, 17, 42, 43, 44, 14]). However, we show here it is not 45 robustly designed for very large volume changes. 46 The rst part of this study is a critical analysis of current volumetric 47 hyperelastic formulations, and presents a discussion of their limitations for 48 large volumetric deformations and non-linear applications. In the second 49 part of this study three novel volumetric strain energy formulations are pre- 50 sented (and several variations in the appendix) which o er validity for large 51 volumetric strains as well as exibility for experimental tting of complex 52 behaviour for both the shrinkage and expansion domain. In addition, the 53 third formulation was expanded to include non-monotonic volumetric strain 54 sti ening (e.g. potentially leading to a plateau in the observed stress) as 55 seen for crushable foams and cellular materials. Finally, the models are com- 56 pared to experimental data for neoprene rubber foam [45], exible open-cell 57 polyurethane cushioning foam [17], natural cork [46], and rigid closed-cell 58 polyurethane foam [47]. 59 2. Theoretical background and rationale 60 In hyperelasticity the constitutive behaviour, i.e. the material's stress- 61 strain behaviour, is derived from a formulated strain energy density (SED) 62 function (for a more detailed discussion of these concepts the reader is re- 63 ferred to established text-books on the subject [48] and [35]). In the case 64 of uncoupled formulations the strain energy consists of additive deviatoric 65 (shape changing) dev(~C) and volumetric (volume changing) vol(J) contri- 66 butions: 67 (~C;J) = dev(~C) + vol(J) (1) Here ~CandJrepresent the deviatoric right Cauchy-Green tensor and the 68 volume ratio (or Jacobian), respectively. The Cauchy stress for this formu- 69 lation can be written: 70 =dev+hI (2) 3 where devis the deviatoric stress tensor, his the scalar hydrostatic stress, 71 andIis the identity tensor. his given as 72 h=1 3tr() =p (3) wherep=his commonly referred to as the pressure. For an un-coupled 73 formulation, his determined from the volumetric component of the strain 74 energy density, such that 75 h=@ vol(J) @J(4) This paper focuses on the analysis and development of volumetric strain 76 energy density formulations for large volumetric deformations. We consider 77 both volume reduction ( J < 1) and volume increase ( J > 1), referred to as 78 shrinkage and expansion , respectively. 79 Although many formulations have been proposed for deviatoric strain en- 80 ergy contributions dev(see for instance the review article [6]), relatively 81 few formulations have been proposed for the volumetric strain energy contri- 82 butions vol(J) (for a more detailed discussion of volumetric strain energy 83 formulations the reader is referred to the surveys by Doll and Schweizerhof 84 [10] and Horgan and Murphy [49]). Moreover, volumetric components of 85 hyperelastic models are not typically subjected to rigorous analysis in or- 86 der to ensure that physically realistic behaviour for large volume changes is 87 maintained. The study of Doll and Schweizerhof [10] establishes 9 criteria 88 (summarised as I-IXin Table 1) that should be satis ed in order to ensure 89 physically realistic material behavior during expansion and shrinkage. Here 90 we add a tenth ( Xin Table 1), namely: the volumetric component of a hy- 91 perelastic model should be capable of precisely describing strain sti ening 92 for all values of J(shrinkage and expansion). 93 4 Table 1: Physical constraints and criteria for vol(J) id Description Form I Zero SED in reference state vol(J= 1) = 0 II Zero hydrostatic stress in reference state h(J= 1) = 0 III Positive strain energy density vol(J6= 1)>0 IV Consistent with linear elasticityd2 vol(J=1) dJ2 = V SED approaches1ifJapproaches 0 limJ!0 vol(J) =1 VI Hydrostatic stress approaches 1 ifJapproaches 0 limJ!0h(J) =1 VII SED approaches1ifJapproaches1 limJ!1 vol(J) =1 VIII Hydrostatic stress approaches 1ifJapproaches1 limJ!1h(J) =1 IX Tangent modulus >0 (polyconvexity)d2 vol(J) dJ20 X Control of strain sti ening for all J 2.1. Structure of this paper 94 The current paper is structured as follows. 95 In Section 3 we analyze the capability of four existing models to satisfy 96 the criteria set out in Table 1: 97 In Section 3.1 commonly implemented single parameter models are 98 analysed; 99 In Section 3.2 the formulation by Bischo et al. [7] for hydrostatic 100 compression of rubber is analysed; 101 In Section 3.3 we analyse the modi ed Ogden formulation [9, 50], a 102 simpli ed form of which has been implemented in ABAQUSR (2018, 103 Dassault Syst emes Simulia Corp.); 104 In Section 3.4 we analyse the Ogden-Hill hyperfoam formulation [8], 105 which has been implemented in ABAQUSR , for highly compressible 106 elastomers; 107 In Section 3.5 we analyze a model by Doll and Schweizerhof [10]. 108 In Section 3.6 we analyze the model by Montella et al. [51]. 109 In Section 4 we propose three novel formulations that improve upon existing 110 formulations in terms of the criteria outlined in Table 1: 111 5 In Section 4.1 we expand the single parameter model (of Equation 6) 112 to ful l all criteria of Table 1, and to provide independent control of 113 strain sti ening in shrinkage and expansion; 114 In Section 4.2 we present a formulation that facilitates precise prescrip- 115 tion of "lock-up" strains in expansion and shrinkage; 116 In Section 4.3 we expand the model presented in Section 4.2 to also cap- 117 ture non-monotonic strain sti ening (typically observed in elastomeric 118 foams). 119 This paper focuses on volumetric strain energy density formulations for 120 large volumetric deformations. Properties of several commonly used forms 121 are discussed and three novel formulations are proposed. Although the ar- 122 guments are most readily presented using uncoupled formulations, they can 123 be extended to coupled formulations where the e ective volumetric response 124 can also be separately identi ed. 125 All visualizations presented here were created based on the free and open 126 source MATLABR (R2019b, The MathWorks Inc., Natick, MA, USA) tool- 127 box GIBBON (https://www.gibboncode.org, [52, 53]). Readers interested in 128 exploring MATLABR implementations, and associated visualizations, of all 129 discussed formulations presented here, may explore the following demo which 130 was added to GIBBON as part of this study: DEMO volumetric SEDeval.m . 131 3. Review and critical analysis of current volumetric SED formu- 132 lations 133 3.1. Common single parameter volumetric SED formulations 134 Two commonly used forms for vol(J), in particular for uncoupled for- 135 mulations, are (e.g. [54]): 136 vol(J) = 2ln(J)2(5) and (e.g. [55, 56]): 137 vol(J) = 2(J1)2(6) These are featured in many nite element implementations such as the 138 open source package FEBio [57] and the proprietary software ABAQUSR . 139 6 These formulations have largely been used to model materials that are as- 140 sumed to be nearly incompressible (such as rubbers and soft tissues), for 141 whichJ1. The motivation for these formulations stems largely from their 142 mathematical convenience. Although their performance when J1 is valid, 143 as we shall describe shortly, non-physical behaviour occurs for large volume 144 changes. Table 1 lists several validity criteria and physical constraints (see 145 also [10, 35]) for volumetric strain energy density formulations. Doll and 146 Schweizerhof [10] examined common formulations and showed that equation 147 5 does not conform to criteria VIII (for high expansions the stress approaches 148 0 rather than1) and criteria IX(for expansion the sti ness reduces to zero 149 atJ=e2:718 after which it becomes negative for J > e ). These e ects 150 are summarised in Figure 1. 151 0 2 4 6 8 102468101214 0 2 4 6 8 10-1-0.500.51 0 2 4 6 8 10-0.100.10.20.30.40.5 Figure 1: The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right), for the formulation of equation 5 Furthermore it was demonstrated that equation 6 does not conform to 152 criteria V(for high shrinkage volapproaches =2 rather than1) and cri- 153 teria VIof Table 1 (for high shrinkage happroachesrather than1). 154 These e ects are summarised in Figure 2. 155 7 0 1 2 30.511.52 0 1 2 3-2-1012 0 1 2 300.511.52Figure 2: The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right), for the formulation of equation 6 For incompressible materials the problems outlined above are largely 156 avoided provided that the conditions that J1 can be numerically en- 157 forced with sucient precision. However, when modelling of compressible 158 material behaviour is of interest, the criteria outlined in Table 1 must be sat- 159 is ed. Finally, these formulations exhibit one xed type of strain dependent 160 behaviour and asymmetry in terms of the di erence between shrinkage and 161 expansion, and therefore these formulations do not conform to criteria Xof 162 Table 1. 163 3.2. The Bischo formulation 164 Bischo et al. [7] presents what can be considered a higher-order repre- 165 sentation of equation 6: 166 vol(J) = 2 cosh (J1) 1 =1X m=1 2(m1) (2m)!(J1)2m(7) Where is an additional material parameter. Bischo et al. [7] demonstrates 167 a good t to the experimental hydrostatic compression data for rubber up 168 to 20%. However, this formulation, and related polynomial forms, have the 169 same pitfalls as the form of equation 6, i.e. they present with a nite strain 170 energy, and hydrostatic stress for J= 0, thereby violating criteria Vand VI 171 of Table 1. Furthermore, this formulation does not o er independent control 172 of the behaviour for shrinkage and expansion, and therefore criteria Xof 173 Table 1 is not met. 174 8 3.3. The modi ed Ogden formulation 175 The modi ed Ogden formulation [50, 9] is given by: 176 vol(J) = 2 J 1 + ln(J) (8) , withthe bulk modulus and (with 6= 0) a material parameter con- 177 trolling the degree of non-linearity. The hydrostatic stress can be derived 178 as: 179 h(J) = J(1J ) (9) and the tangent modulus: 180 @2 @J2= J2 ( + 1)J 1 (10) Figure 3 below illustrates the behaviour of this formulation for shrinkage 181 and expansion and for a range of positive and negative values. In Ogden 182 [9] the formulation is presented in relation to volume reductions only, and 183 with > 0. However, these restrictions are not generally enforced, and if 184 =2 is chosen this formulation reduces to the volumetric contribution, 185 implemented in ABAQUSR , for the uncoupled Aruda-Boyce [58] and Van 186 der Waals [59, 60] formulations. For shrinkage this formulation presents 187 with suitable behaviour. For this domain increasing results in an increase 188 in strain sti ening. Reducing has the opposite e ect with severely negative 189 values even inducing a sti ness reduction and a plateauing behaviour (e.g. 190 see graph for =15 in Figure 3). By studying equation 8 and Figure 3 191 it becomes clear this formulation does not conform to all criteria of Table 1 192 in the expansion domain. Speci cally if 1 the tangent tends to zero 193 (i.e. lim J!1d2 vol(J) dJ2= 0), and negative tangents occur if >1 (e.g. see 194 graphs for = 2 and = 15 in Figure 3), thereby violating criteria VII, 195 VIII andIX. Furthermore, since this formulation does not o er independent 196 control over the response for shrinkage and expansion, it does not conform 197 to criteria Xof Table 1. 198 9 0 0.5 1 1.5 200.20.40.60.81/ 0 0.5 1 1.5 2-2-1012 0 0.5 1 1.5 2012345 =-15.0 =-6.0 =-2.0 =2.0 =15.0Figure 3: The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right), for the modi ed Ogden formulation (equation 8). Curves drawn for= 1, and = [15;15] 3.4. The hyperfoam formulation 199 A popular (see for instance [41, 17, 42, 43, 44, 14]) formulation for mod- 200 elling of highly compressible materials is the so called Ogden-Hill or hyper- 201 foam material implemented in ABAQUSR . The formulation, (see [8], [39] 202 page 48, and developments in [40]) is given by: 203 (1;2;3;J) =NX a=12a 2 a  a 1+ a 2+ a 33 +1 a(J a a1) (11) Hereaand aare Ogden-like [35, 9] hyperelastic parameters, and aenables 204 additional enhancement of volumetric contributions. 205 To review the properties of this formulation we restrict ourselves to a rst 206 order formulation ( N= 1). Furthermore, for volumetric deformations, one 207 may use the conditions 1=2=3=J1 3, reducing equation 11 to: 208 (J) =2 2 3(J 31) +1 (J 1) (12) The bulk modulus for this formulation is derived from: 209 = +1 3 (13) Therefore, to ensure >0 one obtains the constraint >1 3. Furthermore, 210 from equation 11, it is clear that a6= 0 is also a constraint. From equation 211 10 12 the hydrostatic stress can be derived as: 212 h(J) =J12 (J 3J ) (14) and the tangent modulus: 213 @2 @J2=J22  ( 31)J 3+ ( + 1)J  (15) Although this formulation is reported to be valid in the domain 1 3< 214 < 0 several issues were revealed in this study. As shown in Figure 4, in 215 this domain one encounters a reduction in the tangent modulus eventually 216 creating negative volumetric sti ness. Furthermore, the stress may reduce 217 to 0 asJ= 0 is approached. It was found that the e ect is exacerbated by 218 the parameter , therefore even if a negative is chosen close to 0 (which 219 appears to provide valid behaviour in Figure 4) a negative sti ness may still 220 occur if the parameter is suciently high. 221 0 0.5 1 1.5 200.20.40.60.81/ 0 0.5 1 1.5 2-2.5-1.2501.252.5 0 0.5 1 1.5 2-3-2-1012345 =-0.30 =-0.25 =-0.20 =-0.15 =-0.10 Figure 4: The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right), for the hyperfoam formulation. Curves drawn for = 1, = 8, and = [0:3;0:1] The illustrated behaviour for the domain 1 3< < 0 when >0 would 222 lead one to add the constraint >0 for this formulation. Figure 5 explores 223 the e ect of varying when >0. In terms of the tangent modulus it may 224 be seen to decay, become constant, or become negative. The hydrostatic 225 stress for 6 is seen to reach a maximum and become constant or reduced 226 with increasing J(due to negative sti ness). It was found that a negative 227 11 tangent modulus may occur when 0 < 6 (see expansion domain for the 228 graphs for 6). Hence to avoid this it appears that >6 is an additional 229 constraint to avoid a negative tangent modulus if >0. 230 1 1.2 1.4 1.6 1.8 200.1250.25/ 1 1.2 1.4 1.6 1.8 200.10.20.3 1 1.2 1.4 1.6 1.8 2-0.2500.250.5 =1.00 =2.00 =3.00 =6.00 =10.00 Figure 5: The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right), for the hyperfoam formulation. Curves drawn for J= [1;2], = 1, = 5, and = [1;10] Figure 6 presents the e ect of varying (when > 0) when > 6 231 ( = 8). It is clear that a positive value enhances the shrinkage domain 232 while suppressing the expansion domain. 233 0 0.5 1 1.5 200.20.40.60.81/ 0 0.5 1 1.5 2-2-1012 0 0.5 1 1.5 201234 =0.01 =0.50 =2.00 =4.00 =8.00 Figure 6: The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right), for the hyperfoam formulation. Curves drawn for = 1, = 8, and = [0:01;8] Figure 7 is similar to Figure 6 except now a negative is explored ( = 234 8). These graphs show that now changes its role to instead enhance the 235 12 expansion domain while suppressing the shrinkage domain. Furthermore, it 236 was observed that a negative tangent modulus may occur if is close to zero 237 (see graph for = 0:01). 238 0 0.5 1 1.5 200.511.52/ 0 0.5 1 1.5 2-2-1012 0 0.5 1 1.5 200.81.62.43.24 =0.01 =0.50 =2.00 =4.00 =8.00 Figure 7: The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right), for the hyperfoam formulation. Curves drawn for = 1, =8, and = [0:01;8] In summary, in the above analysis, several additional constraints have 239 been identi ed for the hyperfoam formulation. It appears > 0 is a con- 240 straint. In addition, if a positive parameter is used, > 6 appears a 241 requirement. If instead negative values are employed a negative tangent 242 modulus might occur for values close to 0. It is therefore concluded that 243 the hyperfoam formulation is highly constrained in terms of the choice of 244 and . 245 For conventional Ogden hyperelastic formulations (see [35]) the param- 246 eter usually controls the degree of non-linarity (or strain hardening) pre- 247 dominantly for the deviatoric behaviour, and for tting, positive or nega- 248 tive parameters may required. However, as demonstrated here, for this 249 Ogden-like hyperfoam formulation, not only strongly in uences the vol- 250 umetric behaviour, it also changes the role of (from enhancing shrinkage 251 to enhancing expansion) when it changes sign. Furthermore, the suggested 252 constraint > 6 (for > 0) may impose a potentially undesirable degree 253 of non-linearity on the deviatoric response. For instance, some materials 254 may demonstrate little strain sti ening such that they require < 6 (this 255 includes Neo-Hookean behaviour, which requires = 2). Finally, even if 256 suitable constraints are implemented, this formulation does not o er inde- 257 pendent control in terms of enhancement for the shrinkage and expansion 258 13 domains, it therefore does not satisfy criteria Xof Table 1. 259 3.5. The Doll and Schweizerhof [10] formulation 260 Doll and Schweizerhof [10] proposed the following volumetric strain en- 261 ergy density formulation: 262 vol(J) = + 1 + 1J +1+1 1J( 1)  ( + 1)( 1)(16) with the material parameter constraints: >0 and >1. Besides satisfying 263 all criteria listed in Table 1, this formulation also o ers some control over 264 the response in shrinkage and expansion. Furthermore, by choosing = 265 the pressure symmetry p(J) =p(1 J) is obtained, and by using = + 2 266 one obtains symmetry in terms of strain energy, i.e. ( J) = (1 J). Figure 267 8 illustrates the e ect of the parameter . It enhances the response for 268 expansion while mildly suppressing the response for shrinkage. 269 0 0.5 1 1.5 200.20.40.60.81/ 0 0.5 1 1.5 2-2-1012 0 0.5 1 1.5 2012345 =0.1 =1.0 =2.0 =3.0 =4.0 Figure 8: The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right), for the Doll and Schweizerhof [10] formulation. Curves drawn for= 1, = 3, and = [0:1;4] The parameter has the opposite e ect, as Figure 9 shows, since it 270 enhances the response for shrinkage while suppressing the response in expan- 271 sion. 272 14 0 0.5 1 1.5 200.20.40.60.81/ 0 0.5 1 1.5 2-2-1012 0 0.5 1 1.5 2012345 =1.1 =2.0 =3.0 =4.0 =5.0Figure 9: The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right), for the Doll and Schweizerhof [10] formulation. Curves drawn for= 1, = 3, and = [1:1;5] Although this formulation o ers a way to control the response for shrink- 273 age and expansion through the parameters and , both parameters have 274 an e ect on both domains. Therefore, since the control is not independent, 275 criteria Xof Table 1 is not fully met. 276 Furthermore, it is noted here that this formulation has the property that 277 the minimum sti ness state need not be at J= 1 (see location of minima in 278 the image on the right of Figure 9). Although this is in principle not invalid, 279 it may not be realistic or desirable. 280 3.6. The Montella et al. [51] formulation 281 Montella et al. [51] proposed the following volumetric strain energy den- 282 sity formulation: 283 vol(J) = 2 1e ln(J)2+2 m 2e 2jln(J)jm(17) Note that it is presented here by making use of tr ln(U) = ln det(U) = 284 ln(J), with U=p Cthe right stretch tensor. To conform to criteria Iof 285 Table 1 (zero strain in the reference state) the following trivial modi cation 286 can be made: 287 vol(J) = 2 1 e 1ln(J)21 +2 m 2 e 2jln(J)jm1 (18) Hereis the bulk modulus, and 2is referred to as the large strain bulk mod- 288 ulus. The parameters 1and 2, are dimensionless and have the constraints 289 11 8and 21 8. 290 15 Figure 10 illustrates the e ect of varying the parameter 1and includes 291 the lower limit 1=1 8. This parameter 1is seen to allow for the simultane- 292 ous variation of response for shrinkage and expansion. The parameter 2has 293 a similar e ect in the second term of equation 18. The parameter 2o ers 294 added control of the slope for higher strains. 295 0 0.5 1 1.5 200.511.52/ 0 0.5 1 1.5 2-5-2.502.55 0 0.5 1 1.5 20246810 1=0.1251=3.0001=6.0001=12.0001=48.000 Figure 10: The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right), for the Montella et al. [51] formulation. Curves drawn for = 1, 2= 1, 2=1 8,m= 4, and 1= [1 8;48 Figure 11 illustrates the e ect of varying the parameter m. The parameter 296 mcan be seen the enhance the response for J >e , andJ <1 eand to suppress 297 the response in the range1 e< J < e . This suppressing/enhancing e ect 298 results in the ability to a plateau region for the expansion domain. If m> 2 299 the second term vanishes for J= 1 hence this appears to be a constraint 300 onmif criteria IVof Table 1 is to be respected. However, in this study 301 m4 appeared a requirement since discontinuities were observed for the 302 tangent modulus when J1 and if 2< m < 4 (see for example the curve 303 form= 2:1 in Figure 11). 304 16 0 1 2 3 400.511.52/ 0 1 2 3 4-5-2.502.55 0 1 2 3 40246810 m=2.1 m=4.0 m=6.0 m=8.0 m=24.0Figure 11: The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right), for the Montella et al. [51] formulation. Curves drawn for = 1, 2= 1, 2=1 8, 1=1 8, andm= [2:124] The Montella et al. [51] formulation of equation 18 o ers a exible formu- 305 lation for modelling of large strain volumetric deformations. Further more, 306 provided the constraints described are respected, the formulation conforms 307 to criteria I-IXof Table 1. However, this formulation does not o er inde- 308 pendent control of the response for shrinkage and expansion. 309 4. The proposed volumetric strain energy density formulations 310 Three novel volumetric strain energy densities are presented in this section 311 which o er separate control over the strain dependent behaviour for shrinkage 312 and expansion. 313 4.1. Formulation 1 314 The rst formulation is inspired by equation 6. A power was added to 315 the volume ratio to enable control of the degree of strain sti ening. Next 316 two terms were created such that one features a positive power and one a 317 negative power, the former being most sensitive to expansion while the latter 318 is most sensitive to shrinkage, leading to: 319 vol(J) = 41 12(J 11)2+1 22(J 21)2 (19) Besides the bulk modulus , this formulation features the material param- 320 eters 1and 2, which control the degree of strain sti ening in terms of 321 expansion and shrinkage respectively, with 2R>0, 12R>2, and 22R>0. 322 17 The hydrostatic stress for this formulation is: 323 h(J) = 2J1 1(J2 1J 1)1 2(J2 2J 2) (20) and the tangent modulus: 324 @2 vol(J) @J2= 2J2 (21 1)J2 1(11 1)J 1 + (2 +1 2)J2 2(1 +1 2)J 2 (21) Figure 12 schematically illustrates the e ect of the parameters , 1, and 325 2, on the hydrostatic stress. 326 Figure 12: A schematic illustration of a typical hcurve illustrating the nature of the parameters , the bulk modulus setting the initial slope, 1, setting the rate of strain sti ening in expansion, 2, setting the rate of strain sti ening in shrinkage. Figure 13 shows the e ect of varying the bulk modulus . 327 18 0 0.5 1 1.5 200.511.52 0 0.5 1 1.5 2-5-2.502.55 0 0.5 1 1.5 20246810 =0.25 =0.50 =1.00 =2.00 =4.00Figure 13: The e ect of . The strain energy density (left), hydrostatic stress (middle), and tangent modulus (right) for formulation 1. Curves drawn for 1= 4, 2= 2,= [0:25;4]. Figure 14 and 15 illustrate the e ect of varying 1and 2respectively, 328 demonstrating near independent control of strain hardening for the expansion 329 and shrinkage domains. 330 0 0.5 1 1.5 200.511.52/ 0 0.5 1 1.5 2-5-2.502.55 0 0.5 1 1.5 20246810 1=2.101=2.501=3.001=4.001=6.00 Figure 14: The e ect of 1. The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right) for formulation 1. Curves drawn for = 1, 2= 2, 1= [2:1;6]. 19 0 0.5 1 1.5 200.511.52/ 0 0.5 1 1.5 2-5-2.502.55 0 0.5 1 1.5 20246810 2=0.102=0.502=1.002=2.002=4.00Figure 15: The e ect of 2. The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right) for formulation 1. Curves drawn for = 1, 1= 4, 2= [0:1;4]. Although near independent control is seen for both the magnitude and 331 degree of strain sti ening of the responses for shrinkage and expansion, it is 332 noted here that, similar to the Doll and Schweizerhof [10] formulation, the 333 minimum sti ness is not guaranteed to be and may not be found at J= 1 334 (Note the shift in the minima for the tangent modulus in Figures 14 and 15). 335 Formulation 2, discussed in the next section, avoids this behaviour. 336 In Appendix A two variations of the above formulation are brie y ex- 337 plored. In equation A.1 of Appendix A.1 a weighting factor was introduced 338 with the aim of providing further control over the dominance of the expan- 339 sion and shrinkage terms. However, this change results in a possible negative 340 tangent modulus for particular parameter choices. 341 To address the fact that the minimum sti ness of formulation 1 only 342 lies atJ= 1 if 2= 1+ 2, equation A.4 of Appendix A.2 describes a 343 switch statement based decoupling of the expansion and shrinkage terms 344 such that the expansion term is used if J1, and the shrinkage term if 345 J <1. Although this alternative form forces the minimum tangent modulus 346 to occur at J= 1, and satis es all constraints listed in Table 1, it presents 347 with a non-smooth sti ness at J= 1, which may not be desirable. 348 4.2. Formulation 2 349 This section discusses a formulation which was inspired by the inverse 350 sigmoid shape of the hydrostatic stress. A tangent function was chosen here 351 since it presents with a rather elegant integral and derivative. The strain 352 energy density for the proposed form is: 353 20 vol(J) =a2ln cosJ1 a! (22) The derivative with Jprovides the hydrostatic stress: 354 h(J) =atanJ1 a (23) The second derivative provides the tangent modulus: 355 @2 vol(J) @J2=sec2J1 a (24) The parameter ais de ned as: 356 a=2 8 >< >:(J11)J1 (J21)J <1(25) This formulation features three material parameters, the bulk modulus and 357 two volume ratio parameters de ning "lock-up" stretches, J1(withJ1>1), 358 andJ2(with 0J2<1). Figure 16 contains a schematic illustration of the 359 nature of these parameters in relation to the hydrostatic stress. 360 21 Figure 16: A schematic illustration of a typical hcurve illustrating the nature of the parameters , the bulk modulus setting the initial slope, J1, setting the maximum volume ratio asymptote, and J2, setting a minimum volume ratio asymptote. BothJ1andJ2de ne a volume ratio at which an asymptote exists for 361 strain energy, hydrostatic stress, and the tangent modulus. Numerical im- 362 plementations therefore should feature the constraints: 363 vol(JJ1) = vol(JJ2) =1 h(JJ1) =h(JJ2) =1 @2 vol(JJ1) @J2 =@2 vol(JJ2) @J2 =1(26) The bulk modulus sets the slope for the hydrostatic stress at J= 1. 364 BeyondJ= 1 the volume ratios J1andJ2determine how rapidly sti ness 365 is enhanced for the expansion and shrinkage domains. If a material exhibits 366 a behaviour such that further volume change beyond a particular point is 367 hindered, this can be modelled using an appropriate choice for these volume 368 ratio asymptote levels. For many materials however J2= 0 is most appro- 369 priate as this is where this asymptote may naturally lie. Clearly if J2= 0 is 370 kept xed this formulation has only the two remaining parameters andJ1. 371 22 Both asymptote levels can be set at a level beyond the expected deformation 372 levels or brought in closer to further enhance strain sti ening. Figure 17 373 illustrates the e ect of varying the bulk modulus and demonstrates how it 374 changes the slope at J= 1 for the hydrostatic stress. 375 0 0.5 1 1.5 200.511.52 0 0.5 1 1.5 2-5-2.502.55 0 0.5 1 1.5 20246810 =0.25 =0.50 =1.00 =1.50 =2.00 Figure 17: The e ect of the bulk modulus. The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right) for formulation 2. Curves drawn forJ2= 0,J1= 2 and= [0:25;2] Figure 18 presents the e ect of varying J1. The parameter J1is seen to 376 shift the location of the asymptote in the expansion domain. 377 0 0.5 1 1.5 200.511.52/ 0 0.5 1 1.5 2-5-2.502.55 0 0.5 1 1.5 20246810 Jmax=1.500 Jmax=1.625 Jmax=1.750 Jmax=1.875 Jmax=2.000 Figure 18: The e ect of the J1. The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right) for formulation 2. Curves drawn for = 1, J2= 0 andJ1= [1:5;2] Figure 19 presents the e ect of varying J2. It is clear how J2enables one 378 to alter the location of the asymptote in the shrinkage domain. 379 23 0 0.5 1 1.5 200.511.52/ 0 0.5 1 1.5 2-5-2.502.55 0 0.5 1 1.5 20246810 Jmin=0.000 Jmin=0.125 Jmin=0.250 Jmin=0.375 Jmin=0.500Figure 19: The e ect of the J2. The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right) for formulation 2. Curves drawn for = 1, J1= 2 andJ2= [0;0:5] From Figures 17, 19, and 18 it is clear that, contrary to the other formu- 380 lations, the minimum tangent modulus is guaranteed to occur at J= 1 and 381 is equal to. 382 Furthermore, one may note that the following simultaneous symmetries 383 exist: 384 vol(Js) = vol(Je); h(Js) =h(Je);@2 vol(Js) @J2=@2 vol(Je) @J2(27) (where subscript sanderefer to shrinkage and expansion respectively), if 385 Je= (Js1)J11 J21+ 1; Js= (Je1)J21 J11+ 1 (28) Therefore if simultaneous symmetry in terms of Jand 1=Jis desired one 386 could use: 387 J2=1 J1(29) Formulation 2 adheres to all criteria of Table 1, with the exception of 388 criteria VIIand VIII, due to the existence of the asymptote at J1in the ex- 389 pansion domain. Indeed it may be deemed unnatural to have the asymptote 390 depart from J= 0 for the shrinkage domain, or to have an asymptote at all 391 for the expansion domain. Appendix A.3 therefore presents a variation to 392 formulation 2 which does not have these features, instead it employs a form 393 similar to equation 7 but with natural asymptotic behaviour added at J= 0. 394 24 4.3. Formulation 3 395 Formulation 3, proposed below, is an extension of formulation 2, of Sec- 396 tion 4.2, to capture non-monotonic strain sti ening, as observed in cellular 397 materials, lattices, and foams [61, 5]. As illustrated in Figure 20, such ma- 398 terials exhibit three main phases during large volumetric compression [61]: 399 I) an initial linear or moderately strain sti ening phase; II) a reduced sti - 400 ness/plateau region due to elastic buckling of the material microstructure; 401 and III) a region of increased sti ness due to densi cation of the structure. 402 Such non-monotonic sti ening behaviour with an elastic buckling plateau 403 region is observed for elastometric foams, such as polyurethane foams (e.g. 404 [13, 17]), and cork (e.g. [62, 5, 63, 47, 5]). 405 Figure 20: The typical response of a cellular solid to shrinkage and expansion. The shrinkage domain typically features several phases, e.g. an initial elastic domain (I), followed by a compaction domain (II), and a densi cation domain (III). As seen in Section 4.2, Formulation 2 features the tanfunction creating a 406 vertical sigmoid shape for the hydrostatic stress. In order to expand formula- 407 tion 2 to allow for a reduced sti ness/plateau region, an additional horizontal 408 sigmoid function is added; in this case a tanh function is used, with asymp- 409 totes parallel to the Jaxis. Conveniently these functions share conceptually 410 similar integrals and derivatives. We propose the following strain energy 411 density function: 412 vol(J) = (1q)a2ln (cosJ1 a ) +qb2ln (coshJ1 b ) (30) 25 resulting in the following expression for hydrostatic stress: 413 h(J) = (1q)atanJ1 a +qbtanhJ1 b (31) with the following expression for tangent modulus: 414 @2 vol(J) @J2= (1q) sec2J1 a +qsech2J1 b (32) The parameters a,b, andqare de ned as: 415 a=2 8 >< >:(J11)J1 (J21)J <1b=1 8 >< >:s1J1 s2J <1q=8 >< >:q1J1 q2J <1 (33) In all cases independent values can be speci ed for shrinkage ( J < 1) and 416 expansion ( J >1). The parameter ais the same as for formulation 2 where 417 J2andJ1set the volume ratios for the two vertical asymptotes. Parameters 418 s2ands1set the hydrostatic stress asymptotes of the horizontal sigmoid 419 function (tanh) in expansion and shrinkage, respectively. The parameters 420 q1andq2set the relative contributions of the monotonic strain sti ening 421 behaviour of the tancomponent and the hydrostatic stress plateau behaviour 422 of thetanh component. If qi= 0 formulation 2 is recovered, with monotonic 423 strain sti ening. Conversely, if qi= 1 a plateau in hydrostatic stress is 424 obtained, but this is not followed by a high sti ness densi cation region. 425 Figure 21 is a schematic illustration to highlight the e ect of the material 426 parameters on the hydrostatic stress. The six physically based parameters 427 can be used to precisely specify the three phases of volumetric deformation 428 described above. 429 26 Figure 21: A schematic illustration of a typical hcurve illustrating the nature of the parameters , the bulk modulus setting the initial slope, J1, setting the maximum volume ratio asymptote, and J2, setting a minimum volume ratio asymptote, s1, setting the softening stress in expansion, q1setting the dominance of the softening in expansion, s2, setting the softening stress in shrinkage, and q2setting the dominance of the softening in shrinkage. We next provide a parametric investigation of the e ect of varying the 430 parameters s1,s2,q1, andq2. The e ect of ,J1andJ2is equivalent to 431 that of formulation 2 (see Figure 17, 19, and 18 respectively) and therefore 432 not repeated graphically here. Figure 22 shows the e ect of varying s1. 433 This parameter sets the plateau stress level for expansion for the horizontal 434 sigmoid function. 435 27 0 0.5 1 1.5 200.10.20.30.40.5 / 0 0.5 1 1.5 2-2-1012 0 0.5 1 1.5 200.511.522.5 s1=0.10 s1=0.25 s1=0.50 s1=1.00 s1=2.00Figure 22: The e ect of the s1. The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right) for formulation 3. Curves drawn for = 1, J1= 2,J2= 0,q1= 0:98,q2= 0:98,s1= [0:1;2],s2= 0:4 Figure 23 shows the e ect of varying s2. This parameter sets the plateau 436 stress level for shrinkage for the horizontal sigmoid function (note that al- 437 though the hydrostatic stress is negative during shrinkage, s2is here de ned 438 as a positive number). 439 0 0.5 1 1.5 200.10.20.30.40.5 / 0 0.5 1 1.5 2-2-1012 0 0.5 1 1.5 200.511.522.5 s2=0.10 s2=0.25 s2=0.50 s2=1.00 s2=2.00 Figure 23: The e ect of the s2. The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right) for formulation 3. Curves drawn for = 1, J1= 2,J2= 0,q1= 0:98,q2= 0:98,s1= 0:4,s2= [0:1;2] Figure 24 presents the e ect of varying q1, which controls the dominance 440 of the sti ness reduction/plateau behaviour in the expansion domain. 441 28 0 0.5 1 1.5 200.10.20.30.40.5 / 0 0.5 1 1.5 2-2-1012 0 0.5 1 1.5 200.511.522.5 q1=0.00 q1=0.40 q1=0.70 q1=0.90 q1=0.98Figure 24: The e ect of the q1. The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right) for formulation 3. Curves drawn for = 1, J1= 2,J2= 0,q1= [0;0:98],q2= 0:98,s1= 0:4,s2= 0:4 Figure 25 presents the e ect of varying q2, which controls the dominance 442 of the sti ness reduction/plateau behaviour in the shrinkage domain. 443 0 0.5 1 1.5 200.10.20.30.40.5 / 0 0.5 1 1.5 2-2-1012 0 0.5 1 1.5 200.511.522.5 q2=0.00 q2=0.40 q2=0.70 q2=0.90 q2=0.98 Figure 25: The e ect of the q2. The normalized strain energy density (left), hydrostatic stress (middle), and tangent modulus (right) for formulation 3. Curves drawn for = 1, J1= 2,J2= 0,q1= 0:98,q2= [0;0:98],s1= 0:4,s2= 0:4 4.4. Fitting to experimental data 444 To illustrate the ability of our new formulation 1-3 to capture experimen- 445 tal hydrostatic compression data, Figure 26 presents ts to data for neoprene 446 rubber foam [45] (1stcolumn), exible open-cell polyurethane cushioning 447 foam [17] (2ndcolumn), natural cork [46] (3rdcolumn), and rigid closed-cell 448 polyurethane foam [47] (4thcolumn). For all ts the bulk-modulus was 449 29 kept xed and was instead derived from the slope calculated for the small 450 strain domain (up to 4% shrinkage). 451 As is evident from Figure 26, an increasing amounts of sti ness reduc- 452 tion/plateau behaviour is observed in the experimental data (from left to 453 right). In the case of the neoprene and open-cell foam the data represents 454 fully elastic recoverable loading associated with elastic buckling of the mi- 455 crostructure (rather than unrecoverable plastic buckling). Formulation 1 456 cannot accurately capture the non-linear monotonic strain sti ening behav- 457 ior of noeprene rubber foam; the high sti ness behaviour at high volumetric 458 strains is accurately predicted, but the sti ness at low volumetric strains is 459 under predicted. In contrast, formulation 2 accurately predicts the neoprene 460 rubber foam behavour for the full range of experimental data. However, 461 the in ection point observed for open-cell polyurethane cushioning foam and 462 natural cork, and the plateau behaviour for closed-cell polyurethane foam 463 are not captured. Formulation 3 is shown to accurately predict the reported 464 experimental data for all four materials. 465 30 0.4 0.6 0.8 1-1.2-0.9-0.6-0.30 0.2 0.4 0.6 0.8 1-0.4-0.3-0.2-0.10 0.4 0.6 0.8 1-0.8-0.6-0.4-0.20 0.6 0.8 1-0.8-0.6-0.4-0.20 0.4 0.6 0.8 1-1.2-0.9-0.6-0.30 0.2 0.4 0.6 0.8 1-0.4-0.3-0.2-0.10 0.4 0.6 0.8 1-0.8-0.6-0.4-0.20 0.6 0.8 1-0.8-0.6-0.4-0.20 0.4 0.6 0.8 1-1.2-0.9-0.6-0.30 0.2 0.4 0.6 0.8 1-0.4-0.3-0.2-0.10 0.4 0.6 0.8 1-0.8-0.6-0.4-0.20 0.6 0.8 1-0.8-0.6-0.4-0.20Figure 26: Fitting of formulation 1 (top row), formulation 2 (middle row), and formulation 3 (bottom row) to experimental hydrostatic compression data. From the left to the right the data was obtained from Bardy et al. [45], Petre et al. [17], Dart et al. [46], and , Maji et al. [47]. 5. Discussion 466 Much attention has been given to the development of deviatoric strain 467 energy density functions due to the traditional focus, of hyperelastic mod- 468 elling, on nearly-incompressible rubber materials (e.g. [33, 35, 34]), and 469 assumed incompressible soft tissue [6]. Comparatively few strain energy den- 470 sity functions have been proposed for large volumetric deformations (e.g. 471 [9, 7, 8, 10, 51]). We demonstrate here that well-established and commonly 472 used volumetric strain energy formulations are either not valid for large vol- 473 umetric deformations, as they (i) do not adhere to criteria I-IXof Table 474 1, or (ii) do not o er sucient control, for either the shrinkage or the ex- 475 pansion domain, for tting of monotonic or non-monotonic strain sti ening 476 behaviour (criteria Xof Table 1). Following a summary and critical analy- 477 sis of common formulations, and the pitfalls they exhibit, we propose three 478 31 novel formulations which uniquely: 1) are valid for large volumetric deforma- 479 tions, 2) o er separate control of the volumetric strain sti ening behaviour 480 during shrinkage (volume reduction) and expansion (volume increase), and 481 3) in the case of formulation 3, o er the ability to capture non-monotonic 482 volumetric sti ening. The presented formulations o er superior exibility for 483 experimental tting of the large volumetric strain behaviour of hyperelastic 484 materials, and are demonstrated to adhere to all physical constraints and 485 criteria listed in Table 1. 486 Formulation 1 (Section 4.1) exhibits control of the magnitude and de- 487 gree of strain sti ening in shrinkage and expansion domains which is 488 not strongly coupled. This presents an incremental improvement of on 489 the model of Doll and Schweizerhof [10], in which the degree of strain 490 sti ening in shrinkage and expansion is strongly coupled. One property 491 of formulation 1 however is that the minimum of the tangent modulus 492 may not occur at J= 1, and is therefore lower than , for a particular 493 choice of parameters. Although this is a property shared with many 494 other formulations, and this does not render the formulation invalid by 495 any means, it may be deemed undesirable or unrealistic given particular 496 experimental data. 497 Formulation 2 (Section 4.2) was developed to exhibit many of the prop- 498 erties of formulation 1 but also guarantees that the minimum sti ness 499 is found at J= 1. This model is formulated using logarithmic and 500 trigonometric functions, and features a bulk modulus to set the ini- 501 tial slope and two controllable asymptotes, one at the volume ratio 502 J1for expansion, and one at the volume ratio J2for shrinkage. For 503 shrinkageJ2can be set at 0 to enable, for instance, in nite strain 504 energy atJ= 0, as is common. However, it is possible to bring 505 J2closer to 1 to enable more rapid sti ening during volume reduc- 506 tion. Similarly J1is the volume ratio at which an asymptote exists 507 for volume expansion. Control of strain sti ening in shrinkage and 508 expansion domains is fully decoupled, i.e. changes in one domain do 509 not in uence the other. In terms of achieving symmetry, formulation 510 2 also enables, through an appropriate choice of parameters, simulta- 511 neous symmetry in terms of vol(J) = vol(1 J),p(J) =p(1 J), and 512 @2 vol(J)=@J2=@2 vol(1 J)=@J2, i.e. the strain energy density, hy- 513 drostatic stress and tangent modulus for a given percentage volume 514 32 increase or decrease can be made to be equivalent. Furthermore, formu- 515 lation 2 ensures that, even for deviations from symmetry, the minimum 516 tangent modulus always occurs at J= 1 and is equal to . Appendix 517 A.3 provides a variation to formulation 2 whereby the asymptote pa- 518 rameters are avoided. 519 Formulation 3 (Section 4.3) extends formulation 2, of Section 4.2, to 520 capture the non-monotonic sti ening reported for cellular materials, 521 lattices, and foams [61, 5]. A horizontal sigmoid function is superim- 522 posed on formulation 2 creating softening behaviour. The parameters 523 s1ands2de ne hydrostatic stress asymptotes on the horizontal sig- 524 moid function for the expansion and shrinkage domain respectively. 525 Furthermore parameters q1orq2de ne the dominance of these poten- 526 tial softening plateaus. Similar to formulation 2, formulation 3 o ers 527 independent control of the behaviour for the expansion and shrinkage 528 domains. 529 Formulation 3 is shown to provide accurate predictions of the non-linear 530 pressure volumetric relationship under hydrostatic compression for four ma- 531 terials, namely: neoprene rubber foam [45], exible open-cell polyurethane 532 cushioning foam [17], natural cork [46], and rigid closed-cell polyurethane 533 foam [47]. As discussed above, the ability to accurately model non-monotonic 534 volumetric shrinkage and expansion will be important for the simulation and 535 design of next-generation lattice materials, including ultraporous sponges 536 [26] graphene foams aerogels (e.g [27, 28, 29, 30]) in which elastic recovery 537 from compressive strains of 90% have been reported [31]. Graphene aerogels 538 can also be 3D printed [64] allowing for the creation of highly elastic, de- 539 formable, and complex lattices structures. Formulation 3 can also be used 540 to simulate non-monotonic volumetric sti ening of compressible biological 541 materials, such as arteries [65], and the myocardium [66]. Formulation 3 can 542 also be extended to account for plastic buckling in the plateau region (as 543 observed for polypropylene foams [67], metallic foams [68], and trabecular 544 bone [69, 70]. 545 Figure 27 presents two examples of highly elastic lattice structures which 546 can be 3D printed in rubber-like polymeric materials. The lattices are sub- 547 jected to hydrostatic deformations. The structure and visualizations are 548 based on dedicated nite element (FEBio 2.9.1 [57]) models (the lattice 549 is meshed using hexahedral elements, solid material is represented as Neo- 550 Hookean, i.e. =c 4 tr(~C)3 + 2ln (J)2, withc= 1 MPa, and = 50 MPa. 551 33 A related demo has been made available open source as part of GIBBON [52]: 552 DEMO febio 0054 lattice hydrostatic 01.m ). The top row in gure 27 is 553 for the regular octet-truss lattice structure, which demonstrates fairly linear 554 behaviour in expansion and non-linear plateauing and densi cation during 555 shrinkage due to elastic buckling of struts. The bottom row is for an octet- 556 truss lattice with initially curved features. Such features are straightened 557 during expansion creating a source of sti ness enhancement. During shrink- 558 age however the initially curved features immediately and gradually continue 559 to bend, resulting in the absence of the more sudden initiation of bending 560 seen in the structure of the top row. 561 34 Figure 27: Two types of octet truss lattice structures subjected to hydrostatic loading. Graphs on the right show hydrostatic stress ( h) as a function of the volume ratio ( J). The initial con guration for each lattice is shown in 3D on the left. 2D views of the initial and loaded con gurations are also visualized schematically within the graphs on the right. The top row is for a regular octet-truss lattice while the bottom row is for an octet-truss lattice with initially curved features. Future work will include the use of the presented formulations for mod- 562 elling of highly compliant 3D printed polymer lattice structures with tai- 563 lorable strain sti ening and densi cation behaviour. Such materials are use- 564 35 ful for the design of custom biomechanical support structures e.g. at the 565 interface between tissue and prosthetic or orthotic devices. 566 6. Acknowledgements 567 This project was funded by a European Union Horizon 2020 Research 568 and Innovation Program, under grant agreement No. 777072. 569 36