Incremental constitutive tensors and strain localization for prestressed elastic lattices: Part I -- quasi-static response
IIncremental constitutive tensors and strain localization forprestressed elastic lattices: Part I - quasi-static response
G. Bordiga , L. Cabras , A. Piccolroaz , and D. Bigoni ∗ DICAM, University of Trento, Trento, Italy DICATAM, University of Brescia, Brescia, Italy
Abstract
A lattice of elastic rods organized in a parallelepiped geometry can be axially loaded up to an arbitraryamount without distortion and then be subject to incremental displacements. Using quasi-static homog-enization theory, this lattice can be made equivalent to a prestressed elastic solid subject to incrementaldeformation, in such a way to obtain extremely localized mechanical responses. These responses canbe analyzed with reference to a mechanical model which can, in principle , be realized, so that featuressuch as for instance shear bands inclination, or emergence of a single shear band, or competition betweenmicro (occurring in the lattice but not in the equivalent solid) and macro (present in both the lattice andthe equivalent continuum) instabilities become all designable features. The analysis of localizations isperformed using a Green’s function-based perturbative approach to highlight the correspondence betweenmicromechanics of the composite and homogenized response of the equivalent solid. The presented re-sults, limited to quasi-static behaviour, provide a new understanding of strain localization in a continuumand open new possibilities for the realization and experimentation of materials exhibiting these extrememechanical behaviours. Dynamic homogenization and vibrational localization are deferred to Part II ofthis study.
Keywords
Homogenization · Ellipticity loss · Shear bands · Lattice buckling
Shear banding and strain localizations, usually found to emerge before failure of materials, are typically ac-companied by large plastic deformation, damage, and possibly fracture. Mechanical features of shear bandsstrongly depend on the tested material, so that for instance shear bands are normally inclined (to the di-rection of tensile stress) less in rocks than in metals. As a consequence, from the modelling point of view,the analysis of these material instabilities is complicated by the fact that (complex and often phenomeno-logical) elastoplastic constitutive laws are to be used for a material which has to be brought through andbeyond several bifurcation thresholds (corresponding for instance to surface instability or cavitation), beforeencountering shear band formation, the latter typically complicated by the simultaneous emergence of elasticunloading zones adjacent to zones of intense plastic loading. From the experimental point of view, sampleshave to be brought to failure, so that experiments cannot be repeated on the same sample and the materialforming the latter cannot be easily changed to analyzed different instability manifestations, for instance insuch a way to alter the shear band inclination.Imagine now a material in which shear banding and other instabilities may occur well inside the elasticrange and far from failure. A material that can be designed to produce shear bands with a desired inclination,or in which shear bands are the first instability occurring at increasing stress, or in which the anisotropy (notimperfections) allows the formation of only one shear band. Imagine that this material would be characterizedby rigorously determined elastic constitutive laws (thus avoiding complications such as the double branchof the incremental constitutive laws of plasticity) and would be, at least in principle, a material realizable(for instance via 3D printing technology) and testable in laboratory conditions. This material would beideal not only to theoretically analyze instabilities, but also to practically realize the ‘architected materials’which are preconized to yield extreme mechanical properties such as foldability, channelled response, andsurface effects [1–3]. The crucial step towards the definition of a class of these materials was made byTriantafyllidis [4–9] and Ponte Castañeda [10–18], who laid down a general framework for the homogenizationof elastic composites and for the analysis of bifurcation and strain localization in these materials. In particular, ∗ Corresponding author: e-mail: [email protected]; phone: +39 0461 282507. a r X i v : . [ phy s i c s . c l a ss - ph ] J a n i) they showed how to realize an elastic material displaying a prestress-sensitive incremental response, exactlyhow it is postulated for nonlinear elastic solids subject to incremental deformation, and (ii) provided a newunderstanding of strain localization phenomena, showing that a global bifurcation of a lattice structurecorresponds to a loss of ellipticity of the equivalent continuum, while the latter is unaffected by a localbifurcation occurring in the composite.The aim of the present article is to extend the mentioned findings to lattices of elastic rods of arbitrarygeometry and subject to nonlinear axial deformation of the elements, so to explore shear band formation andlocalization by applying a perturbative approach [19], both to the lattice and to the equivalent continuum.In particular, a lattice of elastic rods organized in a parallelepiped network is an example of a compositewhich may be arbitrarily preloaded without introducing grid distortion, so that homogenization allows toobtain rigorous results, showing how a prestressed composite material can react to incremental displacementsas an equivalent elastic continuum. A quasi-static approach to homogenization, based on a strain energyequivalence between the lattice and the continuum, is developed to analyze a generic lattice , so that itbecomes possible to obtain the infinite-body Green’s function for the homogenized solid and compare theresponse of this solid to an applied concentrated force with the behaviour of the lattice at various levelsof preload. This comparison reveals, first of all, the excellent quality of the homogenization approach (sothat the incremental displacement fields found in the lattice and in the homogenized material are practicallycoincident) and highlights the features of shear banding, so that this instability is on the one hand given a clearinterpretation in terms of structural global instability of the lattice and on the other sharply discriminatedfrom local instabilities in the composite, which remain undetected in the continuum. Examples of suchinstabilities, ‘invisible’ in the equivalent material, are provided, which exhibit an ‘explosive’ character, sothat extend from a punctual perturbation to the whole lattice. Fig. 1.
Emergence of a periodic micro-bifurcation (ovalization of the straws’ cross sections, part c), subsequent strain localization(collapse of the straws’ cross sections, part d), and final strain accumulation (parts e and f) during uniaxial deformation of aninitially (parts a and b) hexagonal packing of drinking straws.
An example of local instability, undetected in the homogenized material, but revealed through the analysisof the microstructure, is provided in Fig. 1, where photos of experiments (performed at the Instabilities Labof the University of Trento) are shown in which a package of drinking straws, initially in a regular hexagonaldisposition, is subject to an overall uniaxial strain. The unloaded configuration (Fig. 1a) is not particularlydifferent from the configuration subject to a light loading (Fig. 1b). An increase of the loading yields amicro-bifurcation in terms of a periodic ovalization of the straws’ cross sections (Fig. 1c), while at higher loadstrain localization occurs (in terms of collapse of the cross sections, Fig. 1d), with subsequent strain bandaccumulation (Figs. 1e and 1f). The periodic ovalization is perfectly captured by a bifurcation analysis ofthe hexagonal rods’ grid (Fig. 2) subject to isotropic compression and displaying a periodic bifurcation modewhich is compared with a detail of the photo shown in Fig. 1c. For the geometries investigated in [6] our homogenization approach provides exactly the same results. Moreover, the energyequivalence provides the same results that will be derived in Part II of this article using a Floquet-Bloch dynamic approach. The bifurcation occurs at an axial load in the grid (that was analytically calculated to be − arccos ( − / EJ/l ≈ a) (b) (c) Fig. 2.
The micro-bifurcation mode emerging during the uniaxial deformation of the package of drinking straws shown in Fig. 1is modelled (with the tools provided in this article) as the micro-buckling of an honeycomb lattice of elastic rods, isotropicallyloaded with compressive forces. The equilibrium of the honeycomb structure (a) bifurcates displaying three critical modes (b),which induces a periodic ovalization pattern, explaining the regular and diffuse buckled zones in the array of drinking straws (c).
Homogenization is shown to provide a tool to select the geometry and loading of a lattice in a way toproduce an equivalent solid with arbitrary incremental anisotropy, so that the shear band inclination, or theemergence of a singular shear band can be designed. The results that will be presented also demonstratehow lattice models of heterogeneous materials can be highly effective to obtain analytical expressions forhomogenized properties, thus allowing an efficient analysis of the influence of the microstructural parameters.This is a clear advantage over continuum formulations of composites, where analytical results can only beobtained for simple geometries and loading configurations (as for instance in the case of laminated solids [8, 9,20]). Several new features are found, including a ‘super-sensitivity’ of the localization direction to the preloadstate and the conditions in which a perfect correspondence between the lattice and the continuum occurs(so that the discrete system and the equivalent solid share all the same bifurcation modes). The microscopicfeatures found for the strain localization are shown to share remarkable similarities with the localized failurepatterns observed in honeycombs (as Fig. 1 demonstrates), foams and wood [21–24], while the highly localizeddeformation bands emerging at macroscopic loss of ellipticity are reminiscent of the failure modes observedin balsa wood [25].This article is organized as follows. The derivation of the incremental equilibrium is presented in Section 2for a lattice of elastic rods organized in an arbitrary periodic geometry, while the homogenization is developedin Section 3, providing the incremental constitutive tensor of the effective Cauchy continuum. The stabilityof lattice structure and its relation with the strong ellipticity of the equivalent solid is given in Section 4,while examples and comparisons with the perturbative approach are presented in Sections 5 and 6, wherethe analysis is specialized to a grid of elastic rods arbitrarily inclined and equipped with diagonal springs.Results presented in this article are restricted to quasi-static behaviour, while the important case ofdynamic homogenization (with the Floquet-Bloch technique) and dynamic shear banding is deferred to PartII of this study.
A two-dimensional periodic lattice of elastic rods, deformable in the plane both axially and flexurally, isconsidered, in which all structural members are axially prestressed from an unloaded reference configuration B . The prestress is assumed to be produced by dead loading acting at infinity, while body forces in thelattice are excluded for simplicity. It is assumed that the preload not only satisfies equilibrium, but alsopreserves periodicity and leaves the structure free of flexure. The prestressed configuration B is periodicalong two linearly independent vectors { a , a } , defining the direct basis of the lattice, so that the structurecan be constructed from a single unit cell C , assumed to be composed of N b nonlinear elastic rods withEuler-Bernoulli incremental kinematics, as sketched in Fig. 3. − . EJ/l ) smaller than the load corresponding to loss of ellipticity in the equivalent material (which was calculated throughthe homogenization scheme developed in this article to be ≈ − . EJ/l ). ig. 3. A periodic two-dimensional lattice of (axially and flexurally deformable) elastic rods is considered prestressed fromthe stress-free configuration B (left) by means of a purely axial loading state. The prestressed configuration B (center) can berepresented as the tessellation of a single unit cell along the vectors of the direct basis { a , a } . Upon the current prestressedconfiguration, the incremental response (right) is defined by the incremental displacement field of each rod u ( s ) , here decomposedin an axial and transverse component, u ( s ) and v ( s ) . By considering in-plane flexural and axial incremental deformations, the incremental displacement fieldof the k –th rod in a given unit cell is defined by the vector field (Fig. 3) u k ( s k ) = { u k ( s k ) , v k ( s k ) } (cid:124) , ∀ k ∈ { , ..., N b } , (1)where s k is the coordinate along the k –th rod, u k ( s k ) and v k ( s k ) are the axial and transverse incremental dis-placements. The incremental rotation of the rod’s cross-section θ k ( s k ) is assumed to satisfy the unshearabilitycondition of the Euler-Bernoulli kinematics, namely, θ k ( s k ) = v (cid:48) k ( s k ) . The incremental equilibrium equations for an elastic rod obeying Euler-Bernoulli kinematics and prestressedwith an axial load P (assumed positive in tension) and pre-stretched by λ > , are the following A ( λ ) u (cid:48)(cid:48) ( s ) = 0 , (2a) B ( λ ) v (cid:48)(cid:48)(cid:48)(cid:48) ( s ) − P ( λ ) v (cid:48)(cid:48) ( s ) = 0 , (2b)where A ( λ ) and B ( λ ) are the current axial and bending stiffnesses, respectively, and s ∈ (0 , l ) with l beingthe current length of the rod. It is worth noting that the current axial and bending stiffnesses are, in general,function of the current pre-stretch λ , which in turn depends on the axial load P . In fact, Eqs. (2) governthe incremental equilibrium of an axially pre-stretched rod, and their analytic derivation, accompanied withthe evaluation of A ( λ ) and B ( λ ) from given strain-energy functions, is deferred to Appendix A. In thefollowing, the parameters A ( λ ) and B ( λ ) will simply be denoted as A and B , and treated as independentquantities for generality. The specific example in which the rods composing the lattice are made up of aMooney-Rivlin elastic incompressible material is explicitly reported in Appendix A.Eqs. (2) is a system of linear ODEs for the functions u ( s ) and v ( s ) . As the system is fully decoupled, thesolution is easily obtained in the form u ( s ) = C u + C u s , v ( s ) = C v e − β s + C v e β s + C v s + C v , (3)where { C u , C u , C v , ..., C v } are 6 arbitrary complex constants and β = (cid:112) P/B . For a rod of length l , the following nomenclature can be introduced u (0) = u , v (0) = v , θ (0) = θ , u ( l ) = u , v ( l ) = v , θ ( l ) = θ , (4)so that the vector q = { u , v , θ , u , v , θ } (cid:124) now collects the degrees of freedom of the rod expressed interms of end displacements. Solving the conditions (4) for the constants { C u , C u , C v , ..., C v } (cid:124) allows thesolution (3) to be rewritten as u ( s ) = N ( s ; P ) q , (5)4hich is now a linear function of the nodal displacements q . The × matrix N ( s ; P ) acts as a matrixof prestress-dependent ‘shape functions’ and therefore the representation (5) can also be considered as thedefinition of a ‘finite element’ endowed with shape functions built from the exact solution. Moreover, theseshape functions reduce to the solution holding true in the absence of prestress, because in the limit lim P → N ( s ; P ) = (cid:34) − sl sl ( l − s ) ( l +2 s ) l ( l − s ) sl (3 l − s ) s l s ( s − l ) l (cid:35) , the usual shape functions (linear and Hermitian for axial and flexural displacements, respectively) are re-trieved.By employing Eq. (5), the incremental stiffness matrix of a prestressed rod can be computed, so that forthe k -th rod the elastic strain energy is given by E k = 12 (cid:90) l k (cid:0) A k u (cid:48) k ( s k ) + B k v (cid:48)(cid:48) k ( s k ) (cid:1) ds k = 12 q (cid:124) k (cid:32)(cid:90) l k B k ( s k ; P k ) (cid:124) E k B k ( s k ; P k ) ds k (cid:33) q k , (6)where E k is a matrix collecting the stiffness terms, while B k ( s k ; P ) is the strain-displacement matrix, definedas follows E k = (cid:20) A k B k (cid:21) , B k ( s k ; P k ) = (cid:34) ∂∂s k ∂ ∂s k (cid:35) N k ( s k ; P k ) . The ‘geometric’ contribution of the axial prestress is now included in the potential energy (details areprovided in Appendix A), V gk = 12 P k (cid:90) l k v (cid:48) k ( s k ) ds k = 12 q (cid:124) k (cid:32) P k (cid:90) l k b k ( s k ; P k ) (cid:124) b k ( s k ; P k ) ds k (cid:33) q k , (7)where b k ( s k ; P k ) = (cid:2) ∂∂s k (cid:3) N k ( s k ; P k ) is a vector collecting the derivatives of the shape functions describingthe transverse displacement v . A combination of Eqs. (6) and (7), yields the potential energy for the k -throd in the form V k = E k + V gk . (8)Note that, as the equilibrium equations for the rods have been linearized around an axially pre-loaded configu-ration, the potential (8) represents the incremental potential energy with respect to the current configuration.See the Appendix A for details on the derivation of Eqs. (6)–(8).From Eqs. (6), (7) and (8) the prestress-dependent stiffness matrix is defined as K k ( P k ) = (cid:90) l k B k ( s k ; P k ) (cid:124) E k B k ( s k ; P k ) ds k + P k (cid:90) l k b k ( s k ; P k ) (cid:124) b k ( s k ; P k ) ds k , so that K k = A k l k − A k l k B k l k ϕ ( p k ) B k l k ϕ ( p k ) 0 − B k l k ϕ ( p k ) B k l k ϕ ( p k )0 B k l k ϕ ( p k ) B k l k ϕ ( p k ) 0 − B k l k ϕ ( p k ) B k l k ϕ ( p k ) − A k l k A k l k − B k l k ϕ ( p k ) − B k l k ϕ ( p k ) 0 B k l k ϕ ( p k ) − B k l k ϕ ( p k )0 B k l k ϕ ( p k ) B k l k ϕ ( p k ) 0 − B k l k ϕ ( p k ) B k l k ϕ ( p k ) , where the ϕ j are functions of the non-dimensional measure of prestress p k = P k l k /B k given by ϕ ( p k ) = p / k (cid:0) √ p k − (cid:0) √ p k / (cid:1)(cid:1) , ϕ ( p k ) = p k √ p k coth (cid:0) √ p k / (cid:1) − ,ϕ ( p k ) = p k cosh (cid:0) √ p k (cid:1) − √ p k sinh (cid:0) √ p k (cid:1) √ p k sinh (cid:0) √ p k (cid:1) − (cid:0) √ p k (cid:1) + 8 , ϕ ( p k ) = √ p k (cid:0) sinh (cid:0) √ p k (cid:1) − √ p k (cid:1)(cid:0) √ p k coth (cid:0) √ p k / (cid:1) − (cid:1) sinh (cid:0) √ p k / (cid:1) . Note that the stiffness matrix representative of the lattice reduces, in the limit of vanishing prestress (orunitary pre-stretch λ k = 1 ), to the usual stiffness matrix of an Euler-Bernoulli beam with Hermitian shapefunctions, so that lim p → ϕ j ( p ) = 1 , ∀ j ∈ { , ..., } . .3 Incremental equilibrium of the lattice of elastic rods The current configuration of the lattice unit cell is subject, on the boundary, to the nominal internal incremen-tal actions and incremental displacements transmitted by the rest of the lattice. Therefore, the incrementalpotential energy of the cell can be evaluated through a direct summation of all contributions from the rods,equation (8), over the set of structural elements present inside the cell, plus the incremental work done onthe unit cell boundary by the internal actions f , V ( q ) = N b (cid:88) k =1 V k ( C k q ) − f · q , (9)where q is the vector collecting the degrees of freedom of the unit cell, C k is the connectivity matrix ofthe k -th rod, such that q k = C k q , and f is the vector collecting the generalized (incremental and nominal)internal actions (including bending moments) at the nodes of the unit cell.The (referential) incremental equilibrium equations (in the absence of body forces) are therefore obtainedfrom the stationarity of the potential energy (9) yielding K ( P ) q = f , (10)where K ( P ) = ∂/∂ q (cid:80) N b k =1 V k ( C k q ) is the symmetric (as derived from a potential) stiffness matrix of theunit cell, function of the vector P = { P , ..., P N b } collecting the axial prestress of the rods. The dimensionof the system (10) is N j where N j is the number of nodes in the unit cell. As the preloaded configuration of the lattice is assumed to be spatially periodic, the homogenized incrementalresponse of an equivalent prestressed elastic solid can be defined by computing the average strain-energydensity, associated to an incremental displacement field (defined for the j -th node by the displacement androtation components, respectively, q ( j ) u = { u ( j ) , v ( j ) } (cid:124) and q ( j ) θ = { θ ( j ) } ) which obeys the Cauchy-Bornhypothesis [26–28]. The latter, for a single unit cell, prescribes that the displacement of the lattice nodesbe decomposed into the sum of an affine incremental deformation (ruled by a second-order tensor L ) and aperiodic field (defined by a displacement ˜ q ( j ) u and a rotational ˜ q ( j ) θ component) as q ( j ) u = ˜ q ( j ) u + L x j , q ( j ) θ = ˜ q ( j ) θ , ∀ j ∈ { , ..., N j } (11)where x j is the position of the j -th node. Fig. 4.
The vector q collects the degrees of freedom of the unit cell. In order to impose the periodic boundary conditions requiredby the Cauchy-Born hypothesis (11) (or the quasi-periodic boundary conditions required by the Floquet-Bloch hypothesis (27)that will be used in the bifurcation analysis), q is partitioned between sets of inner nodes q i , boundary nodes located at corners { q lb , q lt , q rb , q rt } , and boundary nodes on the edges { q l , q b , q r , q t } . The corresponding vector of incremental internal nominalactions f is partitioned in the same way. The periodic term ˜ q satisfies ˜ q ( p ) = ˜ q ( q ) for all { p, q } such that x q − x p = n j a j (with n j ∈ { , } ), andcan be expressed in terms of its independent components through a partition of the degrees of freedom, made6n accordance with the location of the nodes present inside the unit cell (see Fig. 4), as ˜ q = ˜ q i ˜ q l ˜ q b ˜ q lb ˜ q r ˜ q t ˜ q rb ˜ q lt ˜ q rt = I I I
00 0 0 I I I
00 0 0 I I I ˜ q i ˜ q l ˜ q b ˜ q lb , (12a)which may succinctly be rewritten as ˜ q = Z ˜ q ∗ , (12b)where Z and ˜ q ∗ are defined according to Eq. (12a). The vector ˜ q in Eq. (12a) has been partitioned to highlightthe inner and boundary nodes according to the notation introduced in Fig. 4. The same partitioning is alsoused for the vectors q and f .In order to enforce the Cauchy-Born conditions into the equations of incremental equilibrium (10), it isconvenient to rewrite Eq. (11) as q ( ˜ q ∗ , L ) = Z ˜ q ∗ + ˆ q ( L ) , (13)where the affine part of the deformation ˆ q ( L ) is a vector-valued function linear in L and such that ˆ q ( L ) ( j ) u = L x j , ˆ q ( L ) ( j ) θ = , ∀ j ∈ { , ..., N j } , where the same notation introduced with Eq. (11) has been used do that the subscript u (subscript θ ) denotesdisplacement (rotations) components.Note that, since the lattice is subject to a non-vanishing prestress state, the macroscopic incrementaldeformation gradient defined by L must be an arbitrary second-order tensor and not constrained to besymmetric (as it happens in the absence of prestress [13, 27, 28]). As explained in the next section, thisunsymmetry is essential for the correct evaluation of the incremental fourth-order tensor defining the effectivecontinuum, ‘macroscopically equivalent’ to the lattice. Before introducing the homogenization technique, it is important to recall that, as shown in Section 2, theequilibrium equations for the lattice are (i) obtained in the context of a linearized theory, and (ii) referred toa prestressed reference configuration, therefore, the unknown ‘equivalent’ continuum has to be formulated inthe context of the incremental theory of nonlinear elasticity by means of a relative Lagrangian description [29].As a consequence, the response of the effective material is defined by an incremental constitutive law in theform ˙ S = C [ L ] , (14)relating the increment of the first Piola-Kirchhoff stress ˙ S to the gradient of incremental displacement L ,through the elasticity tensor C . The most general form for the constitutive tensor C is C = E + I (cid:2) T in components C ijkl = E ijkl + δ ik T jl , (15)where δ ik is the Kronecker delta, T is the Cauchy stress, defining the prestress , and E is a fourth-order elastictensor, endowed with all usual (left and right minor and major) symmetries E ijkl = E jikl = E ijlk = E klij , (16)so that C lacks the minor symmetries but has the major symmetry. The symmetries of C explain the reasonwhy the full incremental deformation gradient L , and not only its symmetric part, appears in the Cauchy-Born hypothesis (11) of the lattice. Moreover, Eq. (15) shows that L can be restricted to be symmetric onlyin the absence of prestress , T = .The incremental strain-energy density for the prestressed continuum is referred to the prestressed configu-ration and can be expressed in terms of a second-order expansion with respect to the incremental deformationgradient L as follows W ( L ) = T · L (cid:124) (cid:123)(cid:122) (cid:125) W ( L ) + C [ L ] · L / (cid:124) (cid:123)(cid:122) (cid:125) W ( L ) , (17)7here the first-order increment W ( L ) accounts for the work expended by the current prestress state T (dueto the relative Lagrangean description the first Piola-Kirchhoff stress coincides with the Cauchy stress), whilethe second-order term W ( L ) is the strain-energy density associated with the incremental first Piola-Kirchhoffstress given by Eq. (14).It is also worth noting that taking the second gradient of the incremental energy density (17) with respectto L yields the constitutive fourth-order tensor C relating the stress increment to the incremental displacementgradient, while the first gradient provides, when evaluated at L = , the prestress T . The latter propertywill be used to dissect the effect of prestress in the homogenized response of the lattice. The homogenization of the lattice response is based on the equivalence between the average incrementalstrain-energy associated to a macroscopic incremental displacement gradient applied to the lattice and theincremental strain-energy density of the effective elastic material subject to the same deformation. In theclassical homogenization theory, this condition is known as macro-homogeneity condition, or Hill-Mandeltheorem, [13, 27, 30, 31], which provides the link between the microscopic and macroscopic scale.In the following, the macro-homogeneity condition is enforced to obtain the incremental energy den-sity (17) that matches the effective behaviour of the prestressed lattice at first- W ( L ) and at second- W ( L ) order. Thus, the homogenization scheme is based on the following steps:(i) An incremental deformation gradient L is considered, so that the incremental energy density for theunknown equivalent continuum is defined by Eq. (17);(ii) following the Cauchy-Born hypothesis, Eq. (13), the incremental displacement field for the lattice isprescribed by the given tensor L and the periodic vector ˜ q ∗ necessary to enforce the equilibrium of thelattice;(iii) with the solution of the lattice in terms of L (the periodic vector ˜ q ∗ becomes in solution a function of L ) the incremental energy density is calculated for the lattice;(iv) the two incremental energy densities in the continuum and in the lattice are matched, so to obtain theparameters defining the equivalent solid. Determination of the periodic displacement field for the lattice.
By substituting condition (13)into Eqs. (10) and pre-multiplying by Z (cid:124) , the incremental equilibrium becomes Z (cid:124) K ( P ) Z ˜ q ∗ + Z (cid:124) K ( P ) ˆ q ( L ) = Z (cid:124) f , (18)where the right-hand side can be written more explicitly using the partitioning introduced in Fig. 4 as Z (cid:124) f = f i f l + f r f b + f t f lb + f rb + f lt + f rt . The fact that the only non-vanishing forces are assumed to be the internal actions transmitted at the unitcell boundary by the neighboring cells implies f i = . Moreover, as the displacement field satisfying theCauchy-Born hypothesis generates internal forces in the infinite lattice that are periodic along the directbasis { a , a } , any single unit cell is subject to external boundary forces that are anti-periodic . Consequently f l = − f r , f b = − f t and f lb = − f rb − f lt − f rt , so that the term Z (cid:124) f vanishes and Eq. (18) becomes Z (cid:124) K ( P ) Z ˜ q ∗ = − Z (cid:124) K ( P ) ˆ q ( L ) . (19)The solution of the linear system (19) provides the incremental strain ˜ q ∗ internal to the lattice for everygiven L . As a consequence of the linearity of ˆ q ( L ) , the solution ˜ q ∗ ( L ) is, in turn, linear in L .A few considerations have to be made about the solvability of the system (19). In fact, it is easy to showthat the matrix Z (cid:124) K ( P ) Z is always singular, regardless of the specific lattice structure under consideration.This is proved by considering a vector ˜ q ∗ = t defining a pure rigid-body translation and observing that K ( P ) Z t = , which, in turn, implies that the dimension of the nullspace of Z (cid:124) K ( P ) Z is at least
2, astwo linearly independent rigid-body translations exist for a 2D lattice. Any other deformation mode, possiblycontained in ker( Z (cid:124) K ( P ) Z ) , is therefore a zero-energy mode (called also with the pictoresque name ‘floppymode’ [32, 33]). These modes are excluded in the following analysis to ensure solvability of system (19), so that ker( Z (cid:124) K ( P ) Z ) contains only two (in the present 2D formulation) rigid-body translations. This exclusion8oes not affect generality, as the analysis of floppy modes can always be recovered in the limit of vanishingstiffness of appropriate structural elements. Note also that sometimes floppy modes can be eliminated orintroduced simply playing with the prestress state (which may induce stiffening or softening [34, 35]).Having excluded floppy modes and observing that the right-hand side of Eq. (19) is orthogonal to ker( Z (cid:124) K ( P ) Z ) , t · Z (cid:124) K ( P ) ˆ q ( L ) = 0 , for all rigid-body translations t , the solution ˜ q ∗ ( L ) can be determined. Match of the second-order incremental strain-energy density and determination of the in-cremental constitutive tensor.
The solution of the linear system (19) allows the incremental displace-ment (13) to be expressed only in terms of the macroscopic displacement gradient L as q ( ˜ q ∗ ( L ) , L ) . Therefore,the second-order incremental strain-energy stored in a single unit cell of the lattice undergoing a macroscopicstrain can be evaluated as follows E ( L ) = 12 q ( ˜ q ∗ ( L ) , L ) · K ( P ) q ( ˜ q ∗ ( L ) , L ) , (20)which is a quadratic form in L , because q ( ˜ q ∗ ( L ) , L ) is linear in L . By equating the second-order strain-energy density of the continuum W ( L ) = C [ L ] · L / to the average energy of the lattice (20), the followingequivalence condition is obtained C [ L ] · L (cid:124) (cid:123)(cid:122) (cid:125) Continuum = 1 |C| E ( L ) (cid:124) (cid:123)(cid:122) (cid:125) Lattice , (21)where |C| is the area of the unit cell.Finally, the second gradient of (21) with respect to L yields the incremental constitutive tensor for theeffective Cauchy material, equivalent to the lattice, in the form C = 1 |C| ∂ E ( L ) ∂ L ∂ L = 12 |C| ∂ ∂ L ∂ L (cid:104) q ( ˜ q ∗ ( L ) , L ) · K ( P ) q ( ˜ q ∗ ( L ) , L ) (cid:105) , (22)which becomes now an explicit function of the prestress state , as well as of all the mechanical parametersdefining the lattice. Match of the first-order incremental strain-energy density and homogenization of the prestressstate.
So far, the incremental constitutive tensor C of a continuum equivalent to a prestressed elastic lattice,Eq. (22), has been obtained through homogenization. It is important now to ‘dissect’ from C the effect ofthe prestress T and, as a consequence, to obtain the tensor E .It will be shown below that the current prestress state T of the homogenized material can be directlylinked to the prestress state P = { P , ..., P N b } of the lattice. In fact, by observing that (21) represents the second-order incremental strain energy, equal to W ( L ) = ˙ S ( L ) · L / , an equivalence analogous to (21) can beobtained considering the first-order increment of the strain energy, W ( L ) = T · L . Thus, the first-order termcan be identified as the average work done by the prestress state P during the lattice deformation q ( ˜ q ∗ ( L ) , L ) induced by L so that the following equivalence can be stated T · L (cid:124) (cid:123)(cid:122) (cid:125) Continuum = 1 |C| f P · q ( ˜ q ∗ ( L ) , L ) (cid:124) (cid:123)(cid:122) (cid:125) Lattice , (23)where the vector f P collects the forces that emerge at the nodes of the unit cell and are in equilibrium withthe axial preload of the elastic rods P in the current configuration assumed as reference . As a consequence,the forces f P are independent of L and linear in P .Equation (23) requires that the work done by axial loads f P for nodal displacements q associated to askew-symmetric velocity gradient L = W be zero, namely f P · q ( ˜ q ∗ ( W ) , W ) = 0 . (24)This statement is a direct consequence of the principle of virtual work for rigid body incremental motions,because q ( ˜ q ∗ ( W ) , W ) represents an incremental rotation of the lattice and f P satisfies equilibrium. Hence,taking into account the property (24), the homogenized prestress T can be obtained as the gradient of theequivalence condition (23) with respect to the symmetric part of L , denoted as D , T = 1 |C| ∂∂ D (cid:104) f P · q ( ˜ q ∗ ( D ) , D ) (cid:105) . (25)9 Stability of prestressed lattices of elastic rods, strong ellipticity,and ellipticity of the effective continuum
Lattice bifurcations are governed by the value of the preload P and they can exhibit deformation modeswith different wavelength. When the wavelength becomes infinite, a ‘global bifurcation’ occurs. Whilein the homogenization procedure periodic conditions are used, the systematic investigation of bifurcationsoccurring in the lattice can be conducted by complementing the incremental equilibrium of the lattice (10)with Floquet-Bloch boundary conditions, which involve displacement fields of arbitrary wavelength [28].Denoting the wave vector as k ∈ R and applying Bloch’s theorem (see [36, 37] for details), Eq. (10)becomes Z ( k ) H K ( P ) Z ( k ) q ∗ = , (26)where symbol H denotes the complex conjugate transpose operation and the matrix-valued function Z ( k ) generalizes Z in Eq. (12) as q = q i q l q b q lb q r q t q rb q lt q rt = I I I
00 0 0 I z I z I
00 0 0 z I z I z z I q i q l q b q lb , q = Z ( k ) q ∗ , (27)in which z j = e i k · a j ∀ j ∈ { , } .Note that conditions (27) represents the generalization of Eq. (12) to displacement fields of arbitrarywavelengths, in fact Z ( ) = Z , therefore they allow bifurcations of arbitrary wavelength to be detected andnot only those occurring at the long-wavelength limit, (cid:107) k (cid:107) → .For a given k , the associated preload state P leading to a bifurcation can be obtained by searching fornon-trivial solutions of the incremental equilibrium (26). Hence, by introducing the notation K ∗ ( P , k ) = Z ( k ) H K ( P ) Z ( k ) , a bifurcation becomes possible when det K ∗ ( P , k ) = 0 . (28)It is worth noting that the matrix K ∗ ( P , k ) is Hermitian, K ∗ ( P , k ) = K ∗ ( P , k ) H , which implies that thedeterminant (28) is always real. Moreover, the periodicity of Z ( k ) implies that this determinant is periodicin the k -space with period [0 , π ] × [0 , π ] in the basis { b , b } reciprocal to { a , a } , so that b i · a j = δ ij .In order to construct the stability domain of a lattice, the critical (in other words first) bifurcation needsto be selected by solving Eq. (28) for the smallest preload spanning over all possible wavelengths. Specifically,by introducing the unit vector ˆ P , which singles out a direction in the preload space, the prestress state isdefined as γ ˆ P for a radial loading, so that the critical bifurcation corresponds to the value γ B defined as γ B = inf γ ≥ (cid:110) γ (cid:12)(cid:12)(cid:12) det K ∗ ( γ ˆ P , η b + η b ) = 0 , < η < π , < η < π (cid:111) . (29)where the periodicity of K ∗ ( P , k ) is used to conveniently restrict to one period the search for the infimumover the k -space. It is worth noting that for a vanishing wave vector, Eq. (28) is always satisfied regardless ofthe preload state, because the nullspace of K ∗ ( P , ) always contains rigid-body translations. These trivialsolutions clearly need to be excluded. Strong ellipticity enforces uniqueness of the incremental problem of a homogeneous and homogeneouslydeformed material subject to prescribed incremental displacement on the whole boundary [38] and correspondsto the positive definiteness of the acoustic tensor (associated to the incremental constitutive tensor C ) definedwith reference to every unit vectors n and g as A ( C ) ( n ) g = C [ g ⊗ n ] n . (30)When the prestress state is null and except in the case of an extreme material, where the stiffness of therods becomes vanishing small [39], the homogenized material response is strong elliptic, which in turn impliesellipticity. 10 ailure of ellipticity corresponds to macro (or global) instabilities, where the bifurcation is characterizedby a wavelength long when compared to the period of the lattice structure, which models a localization ofdeformation in the equivalent continuum. The homogenized material is elliptic (E) as long as the the acoustictensor A ( C ) ( n ) is non-singular for every pair of unit vectors n and g , namely, A ( C ) ( n ) g (cid:54) = . (31)When the acoustic tensor becomes singular, a localization of deformation may occur corresponding to a dyad g ⊗ n . The localization is called ‘shear band’ in the special case g · n = 0 , or ‘compaction band’ or ‘splittingmode’ when g · n = ± .It is assumed that the elastic lattice under consideration is equivalent, at null prestress, to a strong ellipticelastic solid, characterized by a constitutive tensor which is function of the prestress T , in turn through theaxial preload P in the elastic rods, equation (25), namely, A ( C ) ( P , n ) . Therefore, using again the previouslydefined unit vector ˆ P and with reference to an infinite material (or to a material with prescribed displacementon the whole boundary) bifurcations are excluded as long as the response remains strong elliptic, while failureof this condition is simultaneous to failure of ellipticity, which occurs at the value γ E defined as γ E = min γ ≥ (cid:26) γ (cid:12)(cid:12)(cid:12) min n , (cid:107) n (cid:107) =1 (cid:104) det A ( C ) ( γ ˆ P , n ) (cid:105) = 0 (cid:27) . (32) Relation between bifurcations in the lattice and in the effective continuum is that failure of ellip-ticity of the latter corresponds to long-wavelength bifurcations of the former, (cid:107) k (cid:107) → , while all bifurcationsare scanned through equation (29), a circumstance which implies γ B ≤ γ E . Moreover, whenever γ B < γ E thebifurcation occurs at microscopic level and is not detectable in the homogenized material, which can still bestrong elliptic [5, 6, 16]. The geometry of the current, prestressed configuration of a preloaded lattice, selected to apply the previouslydeveloped formalism, is sketched in Fig. 5 and is composed of a rhombic grid (of side l ) of elastic rods,inclined at an angle α , and characterized by the following non-dimensional parameters and A = A = A , Λ = l (cid:112) A/B , Λ = l (cid:112) A/B , where the subscript and are relative to the horizontal and inclined rods,as depicted in Fig. 5b. The direct basis of the periodic structure is denoted by the pair of vectors { a , a } (a) (b) Fig. 5.
Current configuration of a rhombic lattice of preloaded elastic rods (a), with the associated unit cell C (b). The directbasis of the lattice is denoted by the pair of vectors { a , a } (a). Labels 1, 2, and S denote the horizontal rods, the inclinedrods, and the diagonal springs, respectively (b). The spring stiffness, the axial and flexural rigidity of the rods, the preloads P and P , as well as the grid angle α can all be varied to investigate different incremental responses. whose representation with respect to the basis { e , e } (see Fig. 5a) is a = l e , a = l ( e cos α + e sin α ) , while the reciprocal basis { b , b } is defined as a i · b j = δ ij , so that b = ( e − e cot α ) /l , b = e csc α/l , k = η b + η b , with η , being dimensionless components. Theresulting ‘skewed’ grid is also considered stiffened by a diagonal bracing realized by linear springs connectingthe midpoints of the horizontal and inclined rods, as sketched in Fig. 5b. The stiffness of the springs is assumedconstant k S = κA/l , with κ being a dimensionless measure of stiffness.In the configuration shown in Fig. 5, the lattice is subject to a preload state defined by the axial forces P and P , made dimensionless respectively as p = P l /B and p = P l /B , so that a lattice is defined bythe parameter set { α, Λ , Λ , κ, p , p } . Note also that the considered lattice structure includes as a particularcase that of a rectangular grid, analyzed in [6]. The homogenization technique outlined in Section 3 for prestressed lattices of arbitrary geometry can bedirectly applied to the grid of elastic rods shown in Fig. 5. The incremental constitutive tensor is computedvia Eq. (22) and made dimensionless as follows C = Al ¯ C ( p , p (cid:124) (cid:123)(cid:122) (cid:125) prestress , Λ , φ, κ, α (cid:124) (cid:123)(cid:122) (cid:125) microstructure ) , (33a)where φ = B /B (note that Λ = Λ / √ φ ), and the non-dimensional tensor-valued function ¯ C ( p , p , Λ , φ, κ, α ) can be decomposed as ¯ C ( p , p , Λ , φ, κ, α ) = ¯ C G ( p , p , Λ , φ, α ) + ¯ C S ( κ, α ) , (33b)with ¯ C G and ¯ C S being, respectively, the contribution of the rod’s grid and the diagonal springs. The fullexpression for the components of ¯ C G and ¯ C S with respect to the basis { e , e } (sketched in Fig. 5a) is thefollowing (components that have to be equal by symmetry are not reported) ¯ C G = 12 d sin α (cid:18) sinh (cid:18) √ p (cid:19) (cid:18) √ p p φ cosh (cid:18) √ p (cid:19) (cid:0) cos(4 α ) (cid:0) Λ − p φ (cid:1) + 4Λ cos(2 α ) + 11Λ + p φ (cid:1) − (cid:18) √ p (cid:19) (cid:0) cos(4 α ) (cid:0) Λ ( p + p φ ) − p φ (cid:1) + Λ ( p + p φ ) (4 cos(2 α ) + 11) + p φ (cid:1)(cid:19) + p √ p sinh (cid:18) √ p (cid:19) cosh (cid:18) √ p (cid:19) (cid:0) cos(4 α ) (cid:0) Λ − p φ (cid:1) + 4Λ cos(2 α ) + 11Λ + p φ (cid:1)(cid:19) , ¯ C G = 4 sin α cos αd (cid:18) sinh (cid:18) √ p (cid:19) (cid:18) √ p p φ cosh (cid:18) √ p (cid:19) (cid:0) Λ − p φ (cid:1) − (cid:18) √ p (cid:19) (cid:0) Λ ( p + p φ ) − p φ (cid:1)(cid:19) + p √ p sinh (cid:18) √ p (cid:19) cosh (cid:18) √ p (cid:19) (cid:0) Λ − p φ (cid:1)(cid:19) , ¯ C G = − αd (cid:18) sinh (cid:18) √ p (cid:19) (cid:18) (cid:18) √ p (cid:19) (cid:0) cos(2 α ) (cid:0) Λ ( p + p φ ) − p φ (cid:1) + Λ ( p + p φ ) + p φ (cid:1) − √ p p φ cosh (cid:18) √ p (cid:19) (cid:0) cos(2 α ) (cid:0) Λ − p φ (cid:1) + Λ + p φ (cid:1)(cid:19) + p √ p sinh (cid:18) √ p (cid:19) cosh (cid:18) √ p (cid:19) (cid:0) cos(2 α ) (cid:0) p φ − Λ (cid:1) − Λ − p φ (cid:1)(cid:19) , ¯ C G = 4 cos αd (cid:18) sinh (cid:18) √ p (cid:19) (cid:18) (cid:18) √ p (cid:19) (cid:0) p p φ − cos α (cid:0) Λ ( p + p φ ) − p φ (cid:1)(cid:1) + √ p p φ cos α cosh (cid:18) √ p (cid:19) (cid:0) Λ − p φ (cid:1)(cid:19) + p √ p cos α sinh (cid:18) √ p (cid:19) cosh (cid:18) √ p (cid:19) (cid:0) Λ − p φ (cid:1)(cid:19) , ¯ C G = 2 sin αd (cid:18) sinh (cid:18) √ p (cid:19) (cid:18) √ p p φ cosh (cid:18) √ p (cid:19) (cid:0) cos(2 α ) (cid:0) p φ − Λ (cid:1) + Λ + p φ (cid:1) − (cid:18) √ p (cid:19) (cid:0) − cos(2 α ) (cid:0) Λ ( p + p φ ) − p φ (cid:1) + Λ ( p + p φ ) + p φ (cid:1)(cid:19) + p √ p sinh (cid:18) √ p (cid:19) cosh (cid:18) √ p (cid:19) (cid:0) cos(2 α ) (cid:0) p φ − Λ (cid:1) + Λ + p φ (cid:1)(cid:19) , ¯ C G = 4 sin α cos αd (cid:18) sinh (cid:18) √ p (cid:19) (cid:18) √ p p φ cosh (cid:18) √ p (cid:19) (cid:0) Λ − p φ (cid:1) − (cid:18) √ p (cid:19) (cid:0) Λ ( p + p φ ) − p φ (cid:1)(cid:19) + p √ p sinh (cid:18) √ p (cid:19) cosh (cid:18) √ p (cid:19) (cid:0) Λ − p φ (cid:1)(cid:19) , These springs can be seen as added after the lattice has been deformed or as deformed together with the lattice. In theformer case further assumptions need not be introduced, while in the latter, the effects of the preload on the springs has to beneglected in the interest of simplicity. The diagonal springs are used in this example to show that microscopic instabilities mayoccur before macroscopic. C G = 2 cos αd (cid:18) sinh (cid:18) √ p (cid:19) (cid:18) √ p p φ cosh (cid:18) √ p (cid:19) (cid:0) cos(2 α ) (cid:0) p φ − Λ (cid:1) + Λ + p φ (cid:1) − (cid:18) √ p (cid:19) (cid:0) − cos(2 α ) (cid:0) Λ ( p + p φ ) − p φ (cid:1) + Λ ( p + p φ ) + p φ (2 p + p φ ) (cid:1)(cid:19) + p √ p sinh (cid:18) √ p (cid:19) cosh (cid:18) √ p (cid:19) (cid:0) cos(2 α ) (cid:0) p φ − Λ (cid:1) + Λ + p φ (cid:1)(cid:19) , ¯ C G = − αd (cid:18) sinh (cid:18) √ p (cid:19) (cid:18) (cid:18) √ p (cid:19) (cid:0) cos(2 α ) (cid:0) Λ ( p + p φ ) − p φ (cid:1) + Λ ( p + p φ ) + p φ (cid:1) − √ p p φ cosh (cid:18) √ p (cid:19) (cid:0) cos(2 α ) (cid:0) Λ − p φ (cid:1) + Λ + p φ (cid:1)(cid:19) + p √ p sinh (cid:18) √ p (cid:19) cosh (cid:18) √ p (cid:19) (cid:0) cos(2 α ) (cid:0) p φ − Λ (cid:1) − Λ − p φ (cid:1)(cid:19) , ¯ C G = − αd (cid:18) sinh (cid:18) √ p (cid:19) (cid:18) (cid:18) √ p (cid:19) (cid:0) cos α (cid:0) Λ ( p + p φ ) − p φ (cid:1) − p p φ (cid:1) + √ p p φ cos α cosh (cid:18) √ p (cid:19) (cid:0) p φ − Λ (cid:1)(cid:19) + p √ p cos α sinh (cid:18) √ p (cid:19) cosh (cid:18) √ p (cid:19) (cid:0) p φ − Λ (cid:1)(cid:19) , ¯ C G = p √ p d sinh (cid:18) √ p (cid:19) cosh (cid:18) √ p (cid:19) (cid:0) sin α (cid:0) Λ − p φ (cid:1) + sin(3 α ) (cid:0) Λ − p φ (cid:1) + 4 csc( α ) ( p + p φ ) (cid:1) − α sinh (cid:18) √ p (cid:19) (cid:18) (cid:18) √ p (cid:19) (cid:16) cos(2 α ) (cid:0) Λ ( p + p φ ) − p φ (cid:1) + 2 csc ( α ) ( p + p φ ) + Λ ( p + p φ ) − p φ (4 p + 3 p φ ) (cid:17) − √ p p φ cosh (cid:18) √ p (cid:19) (cid:0) cos(2 α ) (cid:0) Λ − p φ (cid:1) + 2 csc ( α ) ( p + p φ ) + (cid:0) Λ − p φ (cid:1)(cid:1)(cid:19) , where d = e − ( √ p + √ p )Λ (cid:16)(cid:16) e √ p ( √ p − √ p +2 (cid:17) (cid:16) e √ p − (cid:17) p φ − (cid:16) e √ p − (cid:17) p (cid:16) e √ p − (cid:17) + (cid:16) e √ p − (cid:17) p (cid:16) e √ p +1 (cid:17) √ p (cid:17) . The constitutive tensor ruling the effect of diagonal springs can be written as ¯ C S = κ α )4 sin α , ¯ C S = ¯ C S = ¯ C S = ¯ C S = ¯ C S = κ cos α , ¯ C S = ¯ C S = ¯ C S = ¯ C S = ¯ C S = ¯ C S = ¯ C S = 12 κ sin α , ¯ C S = ¯ C S = ¯ C S = ¯ C S = 0 . The prestress tensor T , equivalent in the continuum to the preload forces P in the elastic lattice, can beeither calculated using equation (25) or, directly, by computing the average normal and tangential tractionsalong the faces with unit normal e and e . With reference to Fig. 5b the following expression is obtained T = (cid:18) P l sin α + P cos αl sin α (cid:19) e ⊗ e + P cos αl ( e ⊗ e + e ⊗ e ) + P sin αl e ⊗ e . (34) With reference to the lattice sketched in Fig. 5b, the value of the prestress state, which is critical for bifurcationof the grid is determined by employing conditions (32) and (29), and computing numerically the prestressmultipliers γ E and γ B . Results are presented as uniqueness domains in the non-dimensional prestress space { p , p } by fixing the set of geometrical and mechanical parameters { α, Λ , Λ , κ } . The boundary of thestability domain identifies the ‘critical’, namely, the first bifurcation of the incremental equilibrium of thelattice.The dependence on α, Λ , Λ , κ has been analyzed by considering two grid configurations that will bereferred to as the orthotropic grid , with equal slenderness Λ = Λ = 10 , and the anisotropic grid , character-ized by different slendernesses, Λ = 7 and Λ = 15 . For each lattice, the influence of the rods’ inclinationis explored by setting α = π/ , π/ , π/ , π/ , while the stiffness of the springs is investigated in the range κ ∈ [0 , . In this way, the influence of the diagonal bracing on the critical bifurcation mode is analyzed.To investigate both macroscopic (infinite wavelength) and microscopic (finite wavelength) bifurcations,results for the orthotropic grid with α = π/ are reported in Fig. 6, where critical bifurcation loads p and p are reported for the cases in which diagonal springs are absent ( κ = 0 , Fig. 6a, b, c) and for a spring stiffness κ = 0 . (Fig. 6d, e, f).The uniqueness domains (Fig. 6a and 6d) have been computed by solving equation (29) for radial loadingpaths in the non-dimensional load space { p , p } . To clarify the results of this computation, two critical13oundaries are reported, one with a continuous line and the other with a continuous-dotted line, referring tobifurcations of long (infinite) and ‘shortest possible’ wavelength, respectively. The former occurs when theinfimum of (29) is attained at k = , while the latter refers to the infimum computed on the boundary ofthe reciprocal unit cell, η , = ± π . The location of the infimum can be visualized, by fixing the loadingdirection as p = γ ˆ p , and then by numerically computing the bifurcation surface defined as det K ∗ ( γ ˆ p , η b + η b ) = 0 in the space { η , η , γ } . Two radial paths are considered in Fig. 6a and 6d, namely, equibiaxial ˆ p = {− / √ , − / √ } and uniaxial ˆ p = {− , } compression (red dashed lines), and the correspondingbifurcation surfaces are reported in Fig. 6b, c and Fig. 6e, f, respectively. (a) Uniqueness domain ( κ = 0) stability domainmacro-bifurcationmicro-bifurcation (b) Equibiaxial compression (c) Uniaxial compression(d) Uniqueness domain ( κ = 0 . stability domainmacro-bifurcationmicro-bifurcationinstability undetectedby the continuum (e) Equibiaxial compression (f) Uniaxial compression Fig. 6. (a) and (d): Uniqueness/stability domains in the loading space { p , p } for a square grid (with equal slendernessesof the rods Λ = Λ = 10 ), when diagonal springs are absent (upper part) and present (with spring stiffness κ = 0 . , lowerpart). A continuous (a dotted) contour represents the occurrence of macro (of micro) bifurcations, so that the shaded regionscorrespond to strong ellipticity and uniqueness for the equivalent continuum. In the absence of diagonal springs, macro-instabilities, corresponding to ellipticity loss, prevail and always occur before micro bifurcations, while when the diagonalsprings are considered, the situation is more complex so that one or the other instability may be critical. (b, c) and (e, f): withreference to two specific radial loading paths of equibiaxial and uniaxial compression, shown as red dashed lines in (a) and (d),the bifurcation surfaces evidence the solutions for failure of ellipticity in terms of critical dyads n ⊗ g . In the absence of diagonal springs, Fig. 6a reports the uniqueness domain, corresponding to strong ellip-ticity in the solid equivalent to the lattice, showing that (for every loading direction ˆ p ) a macro-bifurcation,in other words an ellipticity loss (referred to the dyad n ⊗ g ), is always reached before micro-bifurcation. Thelatter represents a structural instability for the lattice that cannot be detected in the equivalent continuum.For the two radial loading paths shown in Fig. 6a, the bifurcation surfaces Figs. 6b,c, show that theminimum values of the load multiplier γ are attained at { η , η } = { , } , which corresponds to a macro-bifurcation for the lattice (associated to an infinite wavelength mode), so that the critical prestress multipliers γ E = 7 . and γ E = 5 . lie on the border of ellipticity loss. The two bifurcations correspond respectively totwo orthogonal modes and a single mode.The presence of diagonal springs complicates the situation as reported in Fig. 6d. In this case theuniqueness/stability domains show that micro-bifurcations may sometimes occur within the region of strong Note that all the possible wavelengths have been considered in the computation of the stability domain (as expressed byEq. (29)), but in Fig. 6a and 6d the critical wave vectors k have been found to either be at the origin ( k = ) or on the boundaryof the reciprocal unit cell (shortest wavelengths). ◦ ) and not thecase of uniaxial compression (horizontal radial path). In fact, when the diagonal springs are present, forequibiaxial compression a critical micro-bifurcation occurs, so that Fig. 6e shows that the minimum value ofthe load multiplier, γ B = √ π , is attained at four points, { η , η } = {± π, ± π } , all associated to a bifurcationmode with a finite wavelength, as shown in the inset. For uniaxial compression, Fig. 6f, a macro-bifurcationof the grid, in other words a loss of ellipticity, occurs at γ E = 15 . and the tangent to the bifurcation surfaceat the origin singles out the infinite-wavelength bifurcation mode (shown in the inset and appearing as a rigidtranslation). (a) α = π/ - - - - - - - - - - - - - - (b) α = π/ - - - - - - - - - - - - - - (c) α = π/ - - - - - - - - - - (d) α = π/ ●●●●●● - - - - - - Fig. 7.
Strong ellipticity domains (corresponding to macro-bifurcations, continuous lines) and uniqueness domains for micro-bifurcation of the lattice (circular markers) for an orthotropic lattice of prestressed elastic rods with Λ = Λ = 10 , at differentrod angles α and stiffness κ of the diagonal springs. Points B , B , B on the stability boundaries have been selected for thecomputation of the associated critical bifurcation mode (shown in the insets). Table 1 collects the critical loads and the criticalwave vectors for each bifurcation mode. At small grid angles, for instance the reported value of α = π/ , failure of ellipticitycoincides with micro-bifurcation in the lattice, so that the bifurcation mode is always characterized by an infinite wavelength.For these grid configurations, the direction of ellipticity loss exhibits a ‘super-sensitivity’ with respect to the load directionality,shown in the insets of part (d), reporting the critical dyads n ⊗ g for failure of ellipticity. Further results on uniqueness domains for the orthotropic and the anisotropic grid are reported in Figs. 7and 8, respectively. The strong ellipticity boundary (corresponding to macro-bifurcation) in the equivalentsolid is denoted with a continuous line, while the circular markers identify the boundary of the unique-ness/stability region for micro bifurcations, which remain undetected in the equivalent solid. Moreover,15ritical bifurcation modes have been reported in insets of Figs. 7 and 8, which refer to some specific pointson the stability boundary (labelled as B , B , B in the former figure and B , B , B , B , B in the latter).The critical loads and the critical wave vectors for each bifurcation mode are reported in Table 1. (a) α = π/ - - - - - - - - (b) α = π/ - - - - - - - - (c) α = π/ - - - - - - - - (d) α = π/ ●●●●●●● - - - - - - - - - - - Fig. 8.
As for Fig. 7, except that Λ = 7 , Λ = 15 and that the points on the stability boundary for which the criticalbifurcation modes have been computed are labeled B , B , B , B , B . Note also that, the typical microscopic bifurcation modesof the anisotropic grid exhibit widely different deformations dictated by the prestress direction (see insets in parts a, b, ccorresponding to points labeled B , B , B , B , B and compare for instance mode B to B or B to B ). From Figs. 7 and 8 the following features can be highlighted. • For the orthotropic grid the strong ellipticity boundary is symmetric with respect to the bisector definedby the condition p = p , which is the principal direction of orthotropy for the grid when Λ = Λ (asymmetry which is broken for the anisotropic grid); • For every value of the grid angle α , the effect of the diagonal springs essentially consists in an enlarge-ment of the strong ellipticity region (see the arrow in Fig. 8 denoting increasing values of stiffness κ ); • The stiffening induced by increasing the spring stiffness κ is much more effective for nearly orthogonalgrids ( α ≈ π/ ) than for small values of inclination α (compare Fig. 7a to Fig. 7d and Fig. 8a toFig. 8d); 16 For every value of the spring stiffness κ , the deviation from orthogonality of the grid always reduces thesize of the strong ellipticity region, so that the largest strong ellipticity region is attained for α = π/ .Label Λ Λ α κ p p k cr B
10 10 π/ . − π − π π b + π b B
10 10 π/ . − . − . π b + π b B
10 10 π/ . − . − . π b + π b B π/ . − . − . π b + π b B π/ . − . − . π b B π/ . − . − . π b + π b B π/ . − . − . π b B π/ . − . − . π b + π b B ∞ π/ . − . − . η b ∀ η Table 1.
Critical bifurcation modes k cr for several configurations of the orthotropic ( Λ = Λ = 10 ) and anisotropic ( Λ =7 , Λ = 15 ) lattice. Plots of the corresponding deformation fields are reported as insets in Figs. 7 and 8. The stability boundaries (circular markers in Fig. 7 and 8), evidence the following characteristics. • At small values of spring stiffness κ , the first bifurcation is always global, so that the strong ellipticity andthe stability boundaries coincide independently of the prestress direction ; a feature visible for κ = 0 , . (purple and blue) in Fig. 7a and Fig. 8a); • An increase in the spring stiffness κ leads to a first bifurcation of local type (the critical mode ischaracterized by a finite wavelength), so that the stability region lies inside the elliptic boundary ; • Fig. 7d and Fig. 8d show that, at sufficiently small values of grid angle (for instance at α = π/ ), failureof strong ellipticity dictates the first bifurcation independently of the stiffness of the diagonal springs (see circular markers of the stability boundary overlapping with the elliptic boundary); • The typical microscopic bifurcation modes of the orthotropic grid are characterized by a pure rotationaldeformation of the junctions of the grid (see insets in Fig. 7a, b, c corresponding to the points labeledas B , B , B ); • Typical microscopic bifurcation modes of the anisotropic grid exhibit widely different deformations dic-tated by the prestress direction (see insets in Fig. 8a, b, c corresponding to points labeled B , B , B , B , B and compare for instance mode B to B or B to B ). • A bifurcation always occurs for every lattice geometry at an equibiaxial load { p , p } = {− π , − π } (point B in Fig. 7a) regardless of the values of Λ , Λ , κ , and α . This bifurcation can be explained bythe fact that the normalized load P l /B = − π corresponds to the buckling load of a simply supportedEuler-Bernoulli beam, and thus, when all the rods of an arbitrary grid are prestressed at this level, apurely flexural buckling mode becomes available (shown in the inset of Fig. 7a).Despite the complex influence of the geometrical and mechanical parameters on the stability of the pre-stressed lattice, two important ‘transitions’ characterize the nature of the first bifurcation, namely:(i) a macro-to-micro transition of the critical bifurcation mode occurs at increasing stiffness of the diagonalsprings κ ;(ii) a micro-to-macro transition of the critical bifurcation mode occurs at decreasing the rod’s inclination α .The above transitions will be exploited in Section 6 to investigate the static response induced by a concen-trated force applied to a lattice preloaded closely to a bifurcation (both global and local bifurcations will beconsidered). A remarkable characteristic is associated to the micro-to-macro bifurcation transition obtained at decreasingangle α , namely, a super-sensitivity of the localization band normal, represented by the unit vector n , withrespect to the state of pre-load, while the localization mode g results only weakly affected .17or instance, at α = π/ and sufficiently high spring stiffness κ , the insets in Figs. 7d and 8d show thatthe relative inclinations between the localization band normal n and the localization mode g strongly varyas a function of the load state in the lattice.When the spring stiffness vanishes, κ = 0 , the localization band is essentially set by the grid inclinationas it is almost perfectly aligned parallel to the inclinations and π/ , so that failure of ellipticity occursin a direction n that is almost orthogonal to the rods. On the contrary, at κ = 0 . a single localizationband occurs, whose inclination strongly depends on the load directionality and is essentially unrelated tothe underlying grid pattern (shown in the insets corresponding to κ = 0 . in Figs. 7d and 8d). The super-sensitivity of the localization direction provides an enhanced tunability of the macroscopic localization patternby means of a simple modification of the load applied to the lattice .It is worth noting that the localization direction can also be designed by constructing a lattice with asuitable value of rods’ angle α , but this approach would not be easily reconfigurable, as the structure geometryhas to be defined in advance. Loss of ellipticity in a solid occurs at modes of all (namely, infinite,) wavelengths, while the correspondingcondition in the lattice usually is that bifurcation occurs only in a mode involving an infinite wavelength.In this sense the equivalent continuous body has a response differing from the lattice, a circumstance whichmay be expected as a consequence of the homogenization procedure, which is applied to a discrete lattice.Surprisingly, it is shown in the following that special conditions can be found in which the lattice bifurcatessimilarly to the equivalent continuum, by displaying infinite modes, covering every wavelength. In this casea perfect equivalence between the bifurcation in the lattice structure and failure of ellipticity in the effectivecontinuum occurs .For a square grid (with α = π/ , Λ = 7 , and Λ = 15 ), the perfect equivalence was obtained ata fixed value of load by varying the stiffness of the diagonal springs κ , thus obtaining κ ≈ . . Thisvalue was calculated by numerically solving equation (29) between κ = 0 . and κ = 0 . , because these twovalues pinpoint the threshold of separation between macro and micro bifurcation. For these values of thelattice parameters and loads the bifurcation mode is unique and involves only the infinite wavelength (macrobifurcation) along the curved boundary denoted as (G) in Fig. 8a, while on the boundary denoted as (GL) inthe same figure an infinite number of bifurcation modes of arbitrary wavelength is present for every criticalloading state , as detailed for the point B ∞ in Fig. 9.Fig. 9a reports the three-dimensional plot in the space { η , η , γ } satisfying the bifurcation condition [ofvanishing of the determinant in Eq. (29)] where γ is the loading multiplier for { p , p } = γ {− , − } , so thatthe critical value γ B leading to bifurcation is highlighted as a red line marking the minimum of the bifurcationsurface. A section of this surface at η = 0 is reported in Fig. 9b to show the dependence of the criticalmultiplier on the stiffness κ , so that for κ < . the critical wave vector is k cr = 0 (macro instability), whilefor κ > . the critical wave vector is k cr = π b (micro instability), and for κ = 0 . every wave vectorof the form k cr = η b (with arbitrary η ) identifies a different bifurcation mode occurring at the same loadmultiplier γ B ≈ . . Within this infinite set of bifurcation modes, a few bifurcation modes (see the labelledpoints on the red contour of Fig. 9b) are reported in order to show the transition of the bifurcation modefrom a local bifurcation (Fig. 9c) to a global shear-band type instability (Fig. 9h). The correlation between the incremental response of the lattice and of the equivalent solid is now investigatedclose to the conditions of instability using the ‘perturbative approach’ introduced in [19]. Following thisapproach, the response of the lattice to an applied static concentrated load (in terms of a force or a forcedipole) is numerically evaluated via finite elements (using the commercial code COMSOL Multiphysics (cid:114) )and compared to the response of the equivalent solid by computing the Green’s function associated to theoperator div C [grad( • )] .The two-dimensional Green’s tensor G needed to perturb the equivalent material and corresponding to aDirac delta function centred at x = is [29] G ( ˆ x ) = − π (cid:73) | n | =1 (cid:16) A ( C ) ( n ) (cid:17) − log | ˆ x · n | , (35)18 a) (b) (c) k cr = π b (d) k cr = π/ b (e) k cr = π/ b (f) k cr = π/ b (g) k cr = π/ b (h) k cr = π/ b Fig. 9.
Conditions showing a perfect equivalence between the lattice and the corresponding continuum, so that when thelatter looses ellipticity, the former exhibits bifurcation occurring with infinite modes covering all wavelengths, a situation whichis revealed by the flat line (highlighted in red) in the bifurcation diagram (part a). The perfect equivalence is obtained throughaccurate tuning of the stiffness of the diagonal springs ( κ ≈ . for a square grid with Λ = 7 and Λ = 15 and a loading { p , p } ≈ . {− , − } ). Part (b) represents a section of the bifurcation surface at η = 0 detailing the flat minimum of thecurve occurring at κ ≈ . . Parts (c)–(h) present selected bifurcation modes documenting a transition at increasing wavelengthof the bifurcation modes from a local bifurcation (c) to a shear-band-type instability (h). where the position vector x has been made dimensionless through division by the rod’s length l , so that ˆ x = x /l . Note that G = G (cid:124) due to the symmetry of the acoustic tensor.Numerical simulations are performed to analyze the lattice by considering a finite square computationaldomain of width l , where l is discretized in 10 finite elements with cubic shape functions. The selectedmesh has been defined by performing a number of simulations with three different mesh refinements, namely5, 10, and 20 elements for l , and then adopting 10 elements, as 20 provided no significant improvement, but19ubstantial computational burden. As the numerical simulations are meant to be compared to the infinite-body Green’s function, the size of the domain has been calibrated in order to minimize boundary disturbanceswith clamped conditions at the four edges of the square domain. The governing equation for the prestressedEuler-Bernoulli rod, Eq. (2b), used in the finite element scheme has been implemented by modifying thebending moment contribution with an additional geometric term representing the load multiplied by thetransverse displacement of the rod.The investigation presented below will reveal that:(i) The localization of deformation connected to macro bifurcation in the lattice and to failure of ellipticityin the equivalent solid are strictly similar;(ii) The lattice response close to a micro bifurcation evidences a ‘microscopic’ type of localization , whichremains completely undetected in the homogenized material.These two different mechanical behaviours are analyzed by exploiting the macro-to-micro transition of thefirst bifurcation mode, which is controlled by the increase in the stiffness of the diagonal springs of the latticeconsidered in Section 5. The grid is subject now to an equibiaxial compression loading. Hence, in Section 6.1the lattice is considered in the absence of diagonal springs ( κ = 0 ), while in Section 6.2 the lattice is reinforcedwith a springs’ stiffness κ = 0 . . The lattice configurations selected for the following analysis are reported in Table 2, together with the valuesof the preload p E corresponding to loss of ellipticity in the equivalent continuum (obtained by numericallysolving equation (32) assuming a radial path p = { p , p } ) or, in other words, to a macro bifurcation in thelattice. As explained in the previous section, the stiffness of the diagonal springs is set to zero in order toensure that a macroscopic bifurcation is critical. The table reports also the inclination θ cr of the normal n to the localization band, defined as n E = e cos θ cr + e sin θ cr . Geometry Rods slenderness Symmetry p E θ cr Square α = π/ = Λ = 10 Cubic − . { , } ◦ , ◦ Λ = 7 , Λ = 15 Orthotropic − . { , } ◦ Rhombus α = π/ = Λ = 10 Orthotropic − . { , } . ◦ , . ◦ Λ = 7 , Λ = 15 Anisotropic − . { , } . ◦ Table 2.
Equibiaxial compression loads p E and inclinations θ cr of n E corresponding to failure of ellipticity in the equivalentmaterial, corresponding to a macro-bifurcation in the lattice, for different grid configurations (see Fig. 5), in the absence ofdiagonal springs ( κ = 0 ). A comparison is presented between the response of the lattice loaded with a concentrated force dipoleand a dipole Green’s function of the effective solid, in terms of maps of incremental displacements. Theresults are presented as contour plots in Figs. 10–13, where the color scale has been conveniently normalizedaccording to the norm of the computed displacement field. In the upper part of the figures, results pertainingto the discrete lattice structure are presented, while, in the lower part, results are relative to the equivalentcontinuum, obtained via homogenization. The figures from left to right correspond to the application ofincreasing preloads, which approach the strong ellipticity boundary in the equivalent solid in situations wherefailure of ellipticity corresponds also to the occurrence of a macro bifurcation of infinite wavelength. The part(d) of each figure ( p = 0 . p E ) also illustrates a magnification of the lattice response in the neighborhood ofthe loading zone, thus disclosing the microscopic deformation pattern associated to the extreme mechanicalresponse of the material when close to elliptic boundary.In the conditions analyzed in Figs. 10–13, the equivalent solid is found to be fully representative of thelattice structure, so that approaching failure of ellipticity the perturbative approach reveals, both in thecontinuum and in the real lattice, the formation of localized incremental deformation in the form of singleor double localization bands. These can be horizontal, vertical or inclined. The correspondence between thebehaviour of lattice and of the equivalent continuum is found to be excellent so that the maps reported in theupper part of the figures are practically identical to the corresponding maps in the lower part of the figures.20 ig. 10. Progressive emergence at increasing load of two orthogonal shear bands visible in the displacement field generated bya diagonal force dipole applied to a square lattice (with cubic symmetry, Λ = Λ = 10 , upper part, a–d, simulated via f.e.m.)compared to the response of the homogenized continuum (lower part e–h). From left to right the load increases towards failureof strong ellipticity p E . Shear bands are aligned parallel to the directions predicted at failure of ellipticity ( θ cr = 0 ◦ , ◦ ). Fig. 11.
As for Fig. 10, but for an orthotropic square lattice ( Λ = 7 , Λ = 15 ), where a single and vertical, θ cr = 90 ◦ , shearband forms. ig. 12. As for Fig. 10, but for an orthotropic rhombic lattice ( Λ = Λ = 10 ), where shear bands are inclined at angles θ cr = 88 . ◦ , . ◦ . Fig. 13.
As for Fig. 10, but for an anisotropic rhombic lattice ( Λ = 7 , Λ = 15 ), where a single shear band forms inclined atan angle θ cr = 151 . ◦ . .2 Micro bifurcations in the lattice and effects on the equivalent solid Micro-bifurcations occurring when the equivalent solid is still in the strong ellipticity range are investigated inthis section, with reference to an equibiaxially compressed square lattice with cubic symmetry Λ = Λ = 10 and diagonal springs of stiffness κ = 0 . . With the assumed spring stiffness, a microscopic bifurcation iscritical, namely, it occurs when the equivalent solid is still in the strong elliptic domain. Fig. 14.
Microscopic localization of the bifurcation mode evidenced in the incremental displacement map relative to a squarelattice (cubic symmetry, Λ = Λ = 10 , upper part) and in the equivalent continuum (lower part) at an equibiaxial compressionload corresponding to bifurcation, p B = {− π , − π } , under the action of a ‘quadrupole’ of forces applied at the midpoints of therods. The quadrupole activates a highly localized ‘rotational’ bifurcation mode (labeled as in B in Fig. 7a and Table 1), whichleaves the lattice and the equivalent solid ‘macroscopically’ almost undeformed, while the inter-node deformation is predominant at the scale of the unit cell. The latter feature cannot be detected by the equivalent solid. The incremental displacement maps in the lattice ( at the critical load for micro-bifurcation) and in theequivalent continuum (still in the strong elliptic range) generated by the application of a force quadrupoleare shown in Fig. 14, where the upper parts (lower parts) refer to the lattice (to the continuum) and theparts on the right are a magnification of the zone around the quadrupole shown on the left.The figure shows that the incremental response of the prestressed lattice is highly localized, so thatonly a strong magnification reveals buckling of the elastic rods. Even if the equivalent continuum is not atbifurcation, but still within the uniqueness/stability domain, the distribution of displacements in it somehowresembles that in the lattice, so that the homogenization is still representative of the response of the discretestructure, even though the inter-node deformation cannot be captured.The situation depicted in Fig. 14 completely changes when the lattice is loaded with forces beyond thecritical value for micro bifurcation in the lattice, as shown in Fig. 15, only referred to the lattice loaded with ahorizontal force dipole at different biaxial compression loadings ( at the critical load p = p B for micro-bucklingand beyond, namely, at p = 1 . p B and at p = 1 . p B ).This figure shows displacement maps (upper part) and the corresponding Fourier transform (obtainedvia FFT of nodal displacements, lower part), with superimposed slowness contours corresponding to nullfrequency. The slowness contour (highlighted in red in the figure) was obtained from the bifurcation condition,Eq. (28). The fact that the slowness contour is superimposed to the peaks of the transform (reported whitein the figure), is a validation of the good correspondence between calculations performed via Floquet-Blochand finite element simulations.It can be concluded from Fig. 15 that, while at bifurcation load the incremental perturbation inducedby the force dipole is practically so small and highly localized that results almost invisible, an ‘explosiveinstability’ is found in the lattice, which does not decay and extends to the whole domain occupied by the23 ig. 15. Displacement map (a)–(c) and corresponding Fourier transform (d)–(f) showing the response of the lattice at a loadcorresponding to microscopic instability, p = p B and beyond, p = 1 . p B , . p B . The slowness contour at null frequency,evaluated through the bifurcation condition is superimposed in red. While at the critical load the perturbation is so localized thatresults almost invisible, at higher loads an ‘explosive’ instability involving the whole lattice and extending up to the boundaryof the domain is clearly observed. structure considered in the analysis. This is a special behaviour which remains unobserved in the equivalentcontinuum (still in the strong elliptic range) and cannot therefore be revealed through homogenization. Homogenization of the incremental response of a lattice of elastic rods, axially pre-loaded to an arbitraryamount, has been shown to provide a superb tool for the design of cellular elastic materials of tunableproperties and capable of extreme localized deformations. In particular, the perturbative approach to materialinstability reveals that strain localization in the composite is almost coincident with that occurring in theequivalent solid, which remains unaffected by micro bifurcations, possibly occurring in the lattice. However,the developed homogenization approach allows to play with geometries and stiffnesses of the composite ina way to inhibit or promote strain localization with respect to other micro instabilities. The vibrationalproperties of the lattice and the ability of the homogenization scheme to correctly capture them is a finalcrucial aspect in the design of cellular materials, that will be addressed in Part II of this study.
Acknowledgements
Financial support is acknowledged from: the ERC Advanced Grant ‘Instabilities and nonlocal multiscalemodelling of materials’ ERC-2013-ADG-340561-INSTABILITIES (G.B. and L.C.), PRIN 2015 2015LYYXA8-006 and ARS01-01384-PROSCAN (D.B. and A.P.). The authors also acknowledge support from the ItalianMinistry of Education, University and Research (MIUR) in the frame of the ‘Departments of Excellence’grant L. 232/2016.
References [1] J. T. B. Overvelde, J. C. Weaver, C. Hoberman, and K. Bertoldi. “Rational Design of ReconfigurablePrismatic Architected Materials”. In:
Nature doi : .242] D. M. Kochmann and K. Bertoldi. “Exploiting Microstructural Instabilities in Solids and Structures:From Metamaterials to Structural Transitions”. In: Appl. Mech. Rev doi : .[3] A. Rafsanjani, L. Jin, B. Deng, and K. Bertoldi. “Propagation of Pop Ups in Kirigami Shells”. In: PNAS doi : .[4] N. Triantafyllidis and B. N. Maker. “On the Comparison Between Microscopic and Macroscopic Insta-bility Mechanisms in a Class of Fiber-Reinforced Composites”. In: J. Appl. Mech doi : .[5] G. Geymonat, S. Müller, and N. Triantafyllidis. “Homogenization of Nonlinearly Elastic Materials,Microscopic Bifurcation and Macroscopic Loss of Rank-One Convexity”. In: Arch. Rational Mech. Anal. doi : .[6] N. Triantafyllidis and W. C. Schnaidt. “Comparison of Microscopic and Macroscopic Instabilities in aClass of Two-Dimensional Periodic Composites”. In: J. Mech. Phys. Solids doi : .[7] N. Triantafyllidis and M. W. Schraad. “Onset of Failure in Aluminum Honeycombs under General In-Plane Loading”. In: J. Mech. Phys. Solids doi : .[8] M. D. Nestorović and N. Triantafyllidis. “Onset of Failure in Finitely Strained Layered CompositesSubjected to Combined Normal and Shear Loading”. In: J. Mech. Phys. Solids doi : .[9] M. P. Santisi d’Avila, N. Triantafyllidis, and G. Wen. “Localization of Deformation and Loss of Macro-scopic Ellipticity in Microstructured Solids”. In: J. Mech. Phys. Solids . SI:Pierre Suquet Symposium 97(2016), pp. 275–298. doi : .[10] P. Ponte Castañeda and A. J. M. Spencer. “The Overall Constitutive Behaviour of Nonlinearly ElasticComposites”. In: Proc. R. Soc. A doi : .[11] P. Ponte Castañeda. “The Effective Mechanical Properties of Nonlinear Isotropic Composites”. In: J.Mech. Phys. Solids doi : .[12] P. Ponte Castañeda. “Exact Second-Order Estimates for the Effective Mechanical Properties of Non-linear Composite Materials”. In: J. Mech. Phys. Solids doi : .[13] P. Ponte Castañeda and P. Suquet. “Nonlinear Composites”. In: Advances in Applied Mechanics . Ed.by E. van der Giessen and T. Y. Wu. Vol. 34. Elsevier, 1997, pp. 171–302. doi :
10 . 1016 / S0065 -2156(08)70321-1 .[14] P. Ponte Castañeda. “Second-Order Homogenization Estimates for Nonlinear Composites IncorporatingField Fluctuations: I—Theory”. In:
J. Mech. Phys. Solids doi : .[15] P. Ponte Castañeda. “Second-Order Homogenization Estimates for Nonlinear Composites IncorporatingField Fluctuations: II—Applications”. In: J. Mech. Phys. Solids doi : .[16] O. Lopez-Pamies and P. Ponte Castañeda. “On the Overall Behavior, Microstructure Evolution, andMacroscopic Stability in Reinforced Rubbers at Large Deformations: I—Theory”. In: J. Mech. Phys.Solids doi : .[17] O. Lopez-Pamies and P. Ponte Castañeda. “On the Overall Behavior, Microstructure Evolution, andMacroscopic Stability in Reinforced Rubbers at Large Deformations: II—Application to CylindricalFibers”. In: J. Mech. Phys. Solids doi : .[18] R. Avazmohammadi and P. Ponte Castañeda. “Macroscopic Constitutive Relations for Elastomers Re-inforced with Short Aligned Fibers: Instabilities and Post-Bifurcation Response”. In: J. Mech. Phys.Solids . SI:Pierre Suquet Symposium 97 (2016), pp. 37–67. doi : .[19] D. Bigoni and D. Capuani. “Green’s Function for Incremental Nonlinear Elasticity: Shear Bands andBoundary Integral Formulation”. In: J. Mech. Phys. Solids doi : .[20] A. Bacigalupo and L. Gambarotta. “A Multi-Scale Strain-Localization Analysis of a Layered Stripwith Debonding Interfaces”. In: Int. J. Solids Struct. doi : . 2521] S. D. Papka and S. Kyriakides. “In-Plane Compressive Response and Crushing of Honeycomb”. In: J.Mech. Phys. Solids doi : .[22] S. D. Papka and S. Kyriakides. “Experiments and Full-Scale Numerical Simulations of in-Plane Crushingof a Honeycomb”. In: Acta Materialia doi : .[23] S. D Papka and S Kyriakides. “Biaxial Crushing of Honeycombs: —Part 1: Experiments”. In: Int. J.Solids Struct. doi : .[24] W.-Y. Jang, S. Kyriakides, and A. M. Kraynik. “On the Compressive Strength of Open-Cell MetalFoams with Kelvin and Random Cell Structures”. In: Int. J. Solids Struct. doi : .[25] A. Da Silva and S. Kyriakides. “Compressive Response and Failure of Balsa Wood”. In: Int. J. SolidsStruct. doi : .[26] M. Born and K. Huang. Dynamical Theory of Crystal Lattices . 1955.[27] J. Willis.
Mechanics of Composites . Ecole polytechnique, Département de mécanique, 2002.[28] R. Hutchinson and N. Fleck. “The Structural Performance of the Periodic Truss”. In:
J. Mech. Phys.Solids doi : .[29] D. Bigoni. Nonlinear Solid Mechanics: Bifurcation Theory and Material Instability . Cambridge: Cam-bridge University Press, 2012.[30] R. Hill. “On Constitutive Macro-Variables for Heterogeneous Solids at Finite Strain”. In:
Proc. R. Soc.A doi : .[31] E. Sanchez-Palencia and A. Zaoui, eds. Homogenization Techniques for Composite Media . Vol. 272.Lecture Notes in Physics. Berlin, Heidelberg: Springer Berlin Heidelberg, 1987. doi : .[32] X. Mao and T. C. Lubensky. “Maxwell Lattices and Topological Mechanics”. In: Annu. Rev. Condens.Matter Phys. doi : .[33] L. Zhang and X. Mao. “Fracturing of Topological Maxwell Lattices”. In: New J. Phys. doi : .[34] S. Pellegrino and C. R. Calladine. “Matrix Analysis of Statically and Kinematically IndeterminateFrameworks”. In: Int. J. Solids Struct. doi : .[35] S. Pellegrino. “Analysis of Prestressed Mechanisms”. In: Int. J. Solids Struct. doi : .[36] A. S. Phani, J. Woodhouse, and N. A. Fleck. “Wave Propagation in Two-Dimensional Periodic Lattices”.In: J. Acoust. Soc. Am. doi : .[37] G. Bordiga, L. Cabras, D. Bigoni, and A. Piccolroaz. “Free and Forced Wave Propagation in a Rayleigh-Beam Grid: Flat Bands, Dirac Cones, and Vibration Localization vs Isotropization”. In: Int. J. SolidsStruct.
161 (Apr. 2019), pp. 64–81. doi : .[38] R. Hill. “Acceleration Waves in Solids”. In: J. Mech. Phys. Solids doi : .[39] P. A. Gourgiotis and D. Bigoni. “Stress Channelling in Extreme Couple-Stress Materials Part I: StrongEllipticity, Wave Propagation, Ellipticity, and Discontinuity Relations”. In: J. Mech. Phys. Solids doi : .[40] M. A. Biot. Mechanics of Incremental Deformations . Wiley, 1965.[41] M. Renardy and R. C. Rogers.
An Introduction to Partial Differential Equations . 2nd ed. Texts inApplied Mathematics. New York: Springer-Verlag, 2004.
A Linearized equilibrium of an axially pre-stretched elastica
The linearized equilibrium of an axially stretchable Euler-Bernoulli elastic beam can be obtained througha linearization (around a stretched equilibrium configuration) of the equations governing the equilibrium oflarge deflections and flexure of an elastic rod. 26enoting the stress-free, straight configuration of the elastic rod with a local axial coordinate x , thepotential energy is defined in the reference configuration as V = (cid:90) l ( ψ λ ( λ ) + ψ χ ( χ ) − P u (cid:48) ( x )) dx , (A.1)where l is the initial length, while ψ λ and ψ χ are strain-energy functions for, respectively, axial and flex-ural deformations. The axial stretch λ and the curvature χ are defined by the kinematics of an extensibleunshearable elastica as λ = (1 + u (cid:48) ( x )) cos θ ( x ) + v (cid:48) ( x ) sin θ ( x ) , (A.2a) χ = θ (cid:48) ( x ) = ∂∂x (cid:18) arctan (cid:18) v (cid:48) ( x )1 + u (cid:48) ( x ) (cid:19)(cid:19) , (A.2b)where in (A.2b) the unshearability constraint θ = arctan (cid:16) v (cid:48) u (cid:48) (cid:17) has been explicitly substituted.The linearized response around straight, but axially stretched, configurations can be obtained through thesecond-order expansion of the functional (A.1) with respect to the independent displacement fields { u, v } around the deformed configuration { u , v } = { ( λ − x , } . Hence, by substituting (A.2) into (A.1) andneglecting an arbitrary constant term, the following expansion is obtained V ( u + δu, v + δv ) ∼ (cid:90) l ( ψ (cid:48) λ ( λ ) − P ) δu (cid:48) ( x ) dx ++ 12 (cid:90) l ψ (cid:48)(cid:48) λ ( λ ) δu (cid:48) ( x ) dx + 12 (cid:90) l (cid:18) ψ (cid:48) λ ( λ ) λ δv (cid:48) ( x ) + ψ (cid:48)(cid:48) χ (0) λ δv (cid:48)(cid:48) ( x ) (cid:19) dx , (A.3)where it has been assumed that the residual bending moment is absent in the unloaded configuration ψ (cid:48) χ (0) =0 . Since the first-order term of (A.3) has to vanish when the configuration { u , v } = { ( λ − x , } satisfiesequilibrium, the pre-stretch λ is the solution of the condition ψ (cid:48) λ ( λ ) − P = 0 , indicating that the appliedload P is indeed equal to the axial prestress. Moreover, it is important to note that the second-order termof (A.3) involves the strain energy functions only in terms of second derivatives, ψ (cid:48)(cid:48) λ ( λ ) and ψ (cid:48)(cid:48) χ (0) , evaluatedon the straight stretched configuration.It is now instrumental to update the reference configuration from the stress-free configuration to thestretched configuration, so that the second-order functional (A.3) can be adopted to govern the incrementalresponse of the rod. This can be performed by changing the variable of integration from x to the currentstretched coordinate s = λ x and expressing the fields { u, v } as functions of s . Thus, the second-order termsin eqs (A.3) become V ( u + δu, v + δv ) ∼ (cid:90) l ψ (cid:48)(cid:48) λ ( λ ) λ δu (cid:48) ( s ) ds + 12 (cid:90) l (cid:0) P δv (cid:48) ( s ) + ψ (cid:48)(cid:48) χ (0) λ δv (cid:48)(cid:48) ( s ) (cid:1) ds , (A.4)where l = λ l is the current rod’s length and the symbol (cid:48) has to be understood as differentiation with respectto s . Note also that, as δu (cid:48) ( s ) and δv (cid:48)(cid:48) ( s ) are, respectively, the incremental axial strain and curvature, thecorresponding coefficients are effectively the current value of axial and bending stiffness , so that they can beconcisely denoted as ψ (cid:48)(cid:48) λ ( λ ) λ = A ( λ ) and ψ (cid:48)(cid:48) χ (0) λ = B ( λ ) , both functions of the current axial stretch λ .As this second-order functional has been derived from the large deformation beam theory, it describes thecorrect incremental response superimposed upon a give pre-stretched state. Therefore, the correct form ofthe equilibrium equations governing the incremental displacements can be derived employing the following incremental potential energy V ( u, v ) = 12 (cid:90) l A ( λ ) u (cid:48) ( s ) ds + 12 (cid:90) l (cid:0) P v (cid:48) ( s ) + B ( λ ) v (cid:48)(cid:48) ( s ) (cid:1) ds , (A.5)where now the fields { u ( s ) , v ( s ) } are the current incremental fields and the dependence of the current stiff-nesses A ( λ ) and B ( λ ) on the current axial stretch is highlighted. The governing equations (2) are directlyobtained from (A.5) by imposing the first variation δ V ( u, v ) to vanish. Note that when B ( λ ) is assumedconstant and the axial term of (A.5) is dropped, the usual result of an inextensible Euler-Bernoulli beam isrecovered. Note that, with a little abuse of notation, the symbols for the functions { u, v } have been maintained even though theindependent variable has changed from x to s . .1 Example of a rod made of incompressible hyperelastic materials The incremental potential (A.5) has been derived with reference to the elastica defined by two arbitrarystrain-energy functions governing the current stiffnesses A ( λ ) and B ( λ ) . It is now shown that these twoparameters can be evaluated explicitly for every incompressible elastic material selected to model the lattice’srods.The incremental constitutive response of a rectangular block of incompressible material, deformed underplane strain and initially isotropic can be described (when a uniaxial stress state prevails in the currentconfiguration) through [29] ˙ S = (2 µ ∗ − T ) ∂u ∂x + ˙ p , ˙ S = 2 µ ∗ ∂u ∂x + ˙ p , where ˙ S ij is the increment of the first Piola-Kirchhoff stress and u i are the incremental displacements, µ ∗ the incremental modulus (corresponding to shearing inclined at ◦ with respect to the axes), T the currentuniaxial Cauchy stress ( T = 0 ), and ˙ p the incremental Lagrange multiplier associated to the incompressibilityconstraint. Assuming that plane stress prevails incrementally, ˙ S = 0 , and using the incompressibilityconstraint, ˙ p can be eliminated to yield ˙ S = (4 µ ∗ − T ) ∂u ∂x . (A.6)By considering the incremental equilibrium along the x direction ∂ ˙ S ∂x + ∂ ˙ S ∂x = 0 , an integration over the current thickness h of the block and a subsequent substitution of Eq. (A.6) lead to (cid:90) h/ − h/ ∂ ˙ S ∂x dx = (4 µ ∗ − T ) (cid:90) h/ − h/ ∂ u ∂x dx = 0 , (A.7)where the assumption of vanishing traction at x = ± h/ has been used.The incremental flexural equilibrium can also be retrieved. To this purpose, for a perturbation from thecurrent unixial stress state, Biot [40] has shown that the incremental equilibrium requires ∂ ∂x (cid:90) h/ − h/ x ˙ S dx + T ∂ ∂x (cid:90) h/ − h/ u dx = 0 , (A.8)where the first integral can be recognized to be the incremental bending moment.By adopting the incremental kinematics of an Euler-Bernoulli beam (satisfying the unshearability condi-tion) u ( x , x ) = u ( x ) − x ∂v ( x ) ∂x , u ( x , x ) = v ( x ) , (A.9)and using (A.6), the axial and flexural equilibrium equations (A.7) and (A.8) become (4 µ ∗ − T ) h ∂ u ( x ) ∂x = 0 , (A.10a) (4 µ ∗ − T ) h ∂ v ( x ) ∂x − T h ∂ v ( x ) ∂x = 0 . (A.10b)By noting that T h is the resultant axial load, so that T h = P , the direct comparison between equa-tions (A.10) and (2) provides the identification of the current stiffnesses A ( λ ) and B ( λ ) as A ( λ ) = (4 µ ∗ ( λ ) − T ( λ )) h ( λ ) , B ( λ ) = (4 µ ∗ ( λ ) − T ( λ )) h ( λ ) / , (A.11)where the explicit dependence on the current pre-stretch λ has been highlighted. For instance, for a Mooney-Rivlin material µ ∗ = µ ( λ + λ − ) / and T = µ ( λ − λ − ) , and expressions (A.11) become A ( λ ) = µ ( λ + 3 λ − ) h , B ( λ ) = µ ( λ − + 3 λ − ) h / . with h = h/λ being the initial thickness, and µ the initial shear modulus of the material.28 Regime classification of the effective continuum
The mathematical classification of the PDE describing the incremental equilibrium of the equivalent solidprovides valuable information on the number of localizations available on the elliptic boundary. In fact, thepartial differential equations governing the equilibrium of the effective continuum, in the absence of bodyforces, div C [grad v ] = , (B.1)can be classified according to the following general criterion. Referring to a two-dimensional setting, a solutionof the system (B.1) is selected in a wave form, v = g exp[ i ( x + Ω x )] , (B.2)where g is the wave amplitude and Ω a complex angular frequency. A substitution of (B.2) in the governingequation (B.1) yields the following linear algebraic system (cid:20) C Ω + 2 C Ω + C C Ω + ( C + C )Ω + C C Ω + ( C + C )Ω + C C Ω + 2 C Ω + C (cid:21) (cid:26) g g (cid:27) = (cid:26) (cid:27) . This system has non-trivial solutions if and only if the determinant of the coefficient matrix is equal to zero, (a) Λ = Λ = 10 , α = π/ (b) Λ = Λ = 10 , α = π/ (c) Λ = 7 , Λ = 15 , α = π/ (d) Λ = 7 , Λ = 15 , α = π/ Fig. B.1.
Regime classification of equilibrium PDE for the effective continuum equivalent of a rhombic elastic lattice as thatsketched in Fig. 5 but without diagonal springs. The upper parts (a, b) refer to an orthotropic material ( Λ = Λ = 10 ) material,while the lower parts (c, d) to a completely anisotropic material ( Λ = 7 , Λ = 15 ). The left parts (a, c) refer to a grid withinclination α = π/ and the right α = π/ (b, d). Note that the Hyperbolic and Elliptic regions ‘touch’ at a point. a condition yielding the characteristic equation in the form of a quartic a Ω + 2 a Ω + a Ω + 2 a Ω + a = 0 , (B.3)where a = C − C C ,a = ( C + C ) C − C C − C C ,a = ( C + C ) + 2 C C − C C − C C − C C ,a = − C ( C + C ) + C C + C C ,a = C − C C . Ω j of the quartic (B.3) defines the regime classification according to the followingnomenclature [29, 41]: • In the elliptical regime all the roots Ω j are complex; • In the hyperbolic regime all the roots Ω j are real; • In the parabolic regime two roots are real and two roots are complex.According to this criterion, the regimes for the grid-like lattice of prestressed elastic rods have been classifiedand the results are shown in Fig. B.1 for both orthotropic ( Λ = Λ = 10 ) and anisotropic ( Λ = 7 , Λ = 15 )case. For the sake of brevity, only the case κ = 0= 0