## Abstract

The rigidity of a network of elastic beams is closely related to its microstructure. We show both numerically and theoretically that there is a class of isotropic networks, which are stiffer than any other isotropic network of same density. The elastic moduli of these *stiffest elastic networks* are explicitly given. They constitute upper-bounds, which compete or improve the well-known Hashin–Shtrikman bounds. We provide a convenient set of criteria (necessary and sufficient conditions) to identify these networks and show that their displacement field under uniform loading conditions is affine down to the microscopic scale. Finally, examples of such networks with periodic arrangement are presented, in both two and three dimensions. In particular, we present an *optimal* and *isotropic* three-dimensional structure which, to our knowledge, is the first one to be presented as such.

## 1. Introduction

Many different elastic systems—such as polymer gels, protein networks, cytoskeletal structures [1–6], or wood and bones [7–11]—can be understood as networks of interconnected beams. On a length scale much larger than the typical beam length (‘macroscopic’ scale), such a network can be viewed as a continuous and homogeneous medium characterized by spatially constant elastic moduli. When this medium is subjected to uniform stresses at its boundaries, it will undergo a homogeneous deformation with a constant strain [12]. Such homogeneous deformations are called *affine*. This picture of affine strain is generally valid at large length scales compared to any characteristic inhomogeneities: the macroscopic displacement field , defined as the spatial average of displacements over a sufficiently large domain surrounding every point **r** is affine [13].

On a microscopic level however, the displacement field is generally not affine, and beams can deform by a combination of stretching, bending and twisting mechanisms. In addition, the nature of the junctions between beams plays an important role in the stability and mechanical properties of such systems: networks with junctions that either fix the relative orientation of the beams (‘rigid junctions’) or allow free rotation (‘free hinges’) can have very different mechanical responses [4]. Actually, junctions in most of the real systems have an intermediate behaviour: a small but finite energy cost is associated with the change of their geometry.

The relationship between the mechanical properties of elastic networks on a macroscopic level and the details of their microstructures are the key to optimization and design of lightweight, strong and tough materials [14]. The continuum modelling of such discrete structures has a long history, going as far back as the pioneering work of Cauchy & Poisson [15–18]. The stiffness of such a system clearly depends on its density *ϕ*, defined as the volume of beams per unit volume of material [7–9]. But for a given value of *ϕ*, it is also dramatically affected by the specific spatial arrangement of the elastic phase within the material. The mechanical behaviour and elastic moduli can be computed, especially for quasi-repetitive structures, using homogenization techniques [19,20]. On dimensional grounds [8,21], the volumetric density of strain energy *ε* associated with a stretch-dominated deformation varies linearly with *ϕ*, while it scales as *ϕ*^{2} for the deformation of a three-dimensional network dominated by the beam bending mode. Thus, for the low-density materials considered here (*ϕ*≪1), a structure deforming primarily through the beam stretching mode is usually much stiffer. However, the constant of proportionality between *ε* and *ϕ* still varies significantly among stretch-dominated networks.

In a previous study [22], we have derived bounds on the elastic moduli of a specific class of networks: those made of straight and uniform beams. In this paper, we extend this study and reveal both numerically and theoretically the existence of a class of elastic networks, which are stiffer than *any* other ones having the same symmetry, same density and same elastic phase. Our study is limited to isotropic structures and small strain conditions, but beams can have finite natural curvatures and inhomogeneous cross sections. We show that these *stiffest networks* deform through the stretching mode exclusively and their displacement field is affine down to the microscopic scale *for any infinitesimal loading conditions*. We derive and analyse the necessary and sufficient conditions under which a network belongs to this class and provide examples both at two and three dimensions. To our knowledge, some of them are presented as such for the first time.

The outline of the paper is as follows. In §2, we investigate numerically the role of node connectivity and structural disorder on the affinity of the displacement field and the stiffness of two-dimensional elastic networks. Our simulations show that stretch-dominated networks with same node connectivity can present significantly different elastic moduli and displacement fields. In §3, we propose a theoretical framework to rationalize these observations: we derive bounds on the elastic moduli of isotropic networks using a variational approach and establish simple rules on the geometrical and topological arrangement of beams in a network to make it stiffer for the same amount of material. This theoretical part is an extension of our previous study [22], which was originally restricted to networks made of straight beams with uniform cross sections. In §4, we analyse the restrictions imposed by these structural conditions and show that they rationalize our numerical findings. In particular, they prove that the stiffest networks deform affinely down to the microscopic scale. Finally, in §5, we provide examples of both two-dimensional and three-dimensional regular structures with the highest elastic moduli.

## 2. Simulations

### (a) Method

In order to investigate the interplay between node connectivity, affinity of the strain and stiffness of the network, we simulated the mechanical behaviour of different two-dimensional regular networks, made of straight and uniform beams (same cross section and material). Five different networks have been simulated in order to inspect the effect of both connectivity *z* and disorder on the mechanical properties. Three regular ones: hexagonal (connectivity *z*=3), kagome (*z*=4) and triangular (*z*=6), and two disordered ones: Voronoi (*z*=3) and Delaunay () networks, which are generated from uniformly distributed points on the plane.

The Hamiltonian of a network of beams has three different contributions
2.1The first two sums are carried over all beams (*i*,*j*) of the network and represent the contributions of the stretching and bending energy of the beams, respectively: and , where *l*_{ij} is the length of beam (*i*,*j*), and *u*_{s} and *u*_{n} are the tangential and normal components of the displacement **u**, respectively. *κ*_{s} and *κ*_{b} are the stretching and bending moduli, respectively, related to Young's modulus of the material *E*_{0}, the section of the beam *s* and the second moment of area *I* by: *κ*_{s}=*E*_{0}*s* and *κ*_{b}=*E*_{0}*I*. Following the ideas and method of Head *et al.* [2,3], each network in our simulations is represented by the set of mobile nodes {**x**_{i}} consisting of all beam junctions and midpoints between junctions (the latter so as to include the first bending mode of the beams), with each contribution linearized with respect to changes in the {**x**_{i}}. Let **e**_{ij} be the unit vector pointing from node *i* to node *j*, and **u**_{i} the displacement of node *i*. The discretized version of the stretching and bending contributions read: and , where **u**_{ij}=**u**_{j}−**u**_{i}.

The third sum in equation (2.1) is carried over all junctions *i*, and denotes the energy cost associated with the change of geometry of junction *i*. In the absence of this energy cost (free hinges), torques cannot be sustained and beams deform through the stretching mode exclusively and act as hookean springs. Conversely, if some ‘node rigidity’ is included, other modes of deformation are solicited. Different expressions can be considered for [2–4]. In the simulations presented here, we use , where denotes the set of nodes that are connected to node *i* and Δ*θ*_{ikk′} is the change of the angle between the adjoining beams (*i*,*k*) and (*i*,*k*^{′}) (figure 1). Because of the analogy of this term with the bending energy of a beam, we will refer to this mode of deformation as *node bending* in the following, and *κ*_{n} as the *bending modulus of a junction*. The size of the box *L* is approximately 40 times larger than the average beam length *l*. Range of values used in our simulations for *κ*_{s}, *κ*_{b} and *κ*_{n} are such that *κ*_{s}*l*^{2}≫*κ*_{b}≫*κ*_{n} (typically, *κ*_{s}*l*^{2}∼10^{3}*κ*_{b}, and *κ*_{b}∼10^{3}*κ*_{n}).

Either a simple shear or uniaxial strain *γ* is applied across the top and bottom boundaries,^{1} while periodic boundary conditions are used in the other directions. Within our linearized scheme, is a high-dimensional paraboloid with a unique global minimum, corresponding to the state of mechanical equilibrium. Hence, the linear system can be solved directly using an *LU* decomposition method included in the UMFPACK routines [23]. Once the equilibrium positions are obtained, the associate strain energy can be calculated. Finally, the shear modulus *μ* and longitudinal modulus *M* are obtained by equating the energy with *μγ*^{2}*L*^{2}/2 for simple shear strain and *Mγ*^{2}*L*^{2}/2 for uniaxial strain. Any other elastic moduli can then be calculated. For instance, Young's modulus *E* is related to *μ* and *M* as: *E*=2*μ*(*dM*−2(*d*−1)*μ*)/((*d*−1)*M*−2(*d*− 2)*μ*),^{2} with *d*=2 (resp. *d*=3) for two-dimensional (resp. three dimensions) structures. Owing to the linearization of the equations, the numerical values of the elastic moduli are totally independent of the magnitude of the applied strain *γ*. For the disordered networks, the values of the elastic moduli are averaged over multiple (more than 100) runs.

### (b) Results

For each simulated network, we represent the repartition of energy between beam stretching, beam bending and node bending modes, and measure the affinity of the displacement field. Different parameters have been proposed in the literature [2,3,6,24] to quantify the degree of (non-)affinity of a displacement field, but most of them are global measures of the affinity. We use instead the local affinity measure *m*_{i} defined at every node (junction and midpoint) *i* as , where is the macroscopic, affine, displacement field of the node *i*.

Figure 2 shows the repartition of energy between beam stretching, beam bending and node bending when the networks are subjected to simple shear strain (similar results, not shown here, are obtained for uniaxial strain). Hexagonal and Voronoi networks deform mainly through the node bending mode, while kagome, triangular and Delaunay networks deform mainly through the stretching mode. These results are consistent with Maxwell criterion on rigidity of two-dimensional frameworks [9,25–27]: networks with connectivity above 4 are rigid; they deform primarily through the stretching of beams, while those with connectivity below 4 are soft; that is, they deform through floppy modes if *κ*_{n}=0 [28] (free hinges) or through the bending of beams and nodes if *κ*_{n}≠0. The repartition of energy between the beam and node bending modes depends on the relative importance of *κ*_{b} and *κ*_{n} (rigid junctions with fixed angles correspond to the limit *κ*_{n}≫*κ*_{b}). In our simulations, the bending of beams is quasi-absent (figure 2), because *κ*_{n}≪*κ*_{b}.

Figure 3 shows the affinity of the displacement field for the five networks under simple shear strain. Unsurprisingly, all networks with connectivity less than 4 deform through the bending modes and thus present non-affine displacement fields. On the other hand, the measure of the affinity reveals major differences between networks with connectivity more than or equal to 4. Comparison of triangular and Delaunay lattices specially is intriguing: despite the fact that these two networks share similar topological features (they have same mean connectivity and both are fully triangulated), the triangular network has a pure affine strain, while the Delaunay lattice presents a non-affine strain.

Evaluation of the elastic moduli also shows differences between these two networks: we found *μ*/(*E*_{0}*ϕ*)≃0.091 and *M*/(*E*_{0}*ϕ*)≃0.291 for Delaunay lattice, while *μ*/(*E*_{0}*ϕ*)≃0.125 and *M*/(*E*_{0}*ϕ*)≃0.374 for the triangular lattice, where *ϕ* denotes the network density (finite-size effects on these numerical values are discussed in appendix B). Therefore, while both networks are stretch-dominated, the triangular lattice is significantly (≃27% for simple shear, ≃22% for longitudinal strain) stiffer than the Delaunay lattice. We checked that the numerical values we obtained were independent of *κ*_{n}, consistent with our choice for the range of parameter values (*κ*_{n}≪*κ*_{s}*l*^{2}). Remarkably, the kagome and triangular lattices share similar mechanical properties, although these two lattices are structurally very different (kagome has connectivity *z*=4 and is only partially triangulated): both present an affine displacement field, and the values obtained for their elastic moduli are very close (we obtained for kagome lattice: *μ*/(*E*_{0}*ϕ*)≃0.125 and *M*/(*E*_{0}*ϕ*)≃0.376), as this has been already noted [21].

These results show that no clear connexion can be established between the rigidity of a network and the (mean) connectivity of its nodes, aside from Maxwell criterion. On the other hand, affinity of the displacement field and network stiffness seems correlated: among the five structures studied, the two networks that deform in an affine way (triangular and kagome lattices) are those with the highest elastic moduli. The Delaunay network presents a *piecewise* affine displacement field,^{3} and has significantly lower elastic moduli. Finally, the two soft, bending-dominated, networks (hexagonal and Voronoi) have highly non-affine displacement fields. In §3, we provide a theoretical framework to rationalize the relationships between node connectivity, network stiffness and affinity of the displacement field.

## 3. Stiffest networks: existence, structure and elastic moduli

In the following, we show that there exists a class of networks, which are stiffer than any other of same symmetry and density (and beam Young's modulus). We restrict our analysis to (two dimensional and three dimensional) isotropic structures, although our theoretical framework can be transposed to non-isotropic structures as well. We first establish upper bounds on the elastic moduli of an isotropic network of beams. These bounds coincide with the numerical values obtained for the triangular and kagome lattices. We then derive the structural properties (actually, a set of necessary and sufficient conditions) of optimal networks, i.e. networks that have maximal elastic moduli for a given value of *E*_{0} and *ϕ*. These conditions involve both geometry and topology, which explain why connectivity alone is not enough to predict the macroscopic behaviour of an elastic network.

### (a) Isotropy conditions and bounds on the elastic moduli of isotropic networks

Assuming that on the macroscopic scale, the network is a continuous, homogeneous and isotropic medium, the strain energy per unit volume is (for small strains [12])
3.1where λ is Lamé's first parameter of the network, *μ* the shear modulus (or Lamé's second parameter), the *macroscopic* displacement field (displacement averaged over a domain large compared with any characteristic homogeneities in the network) and are the components of the strain tensor ( are the components of ). Under uniform loading conditions, varies linearly with the position **r**: , where **A** is a matrix that characterizes the strain. For instance, *A*_{αβ}=*γδ*_{αx}*δ*_{βy} for a uniform simple shear strain *γ* in the *xy* plan, and *A*_{αβ}=*γδ*_{αβ} for a uniform radial strain *γ*. According to equation (3.1), the density of strain energy corresponding to such homogeneous strain is
3.2with *a*_{αβ}=(*A*_{αβ}+*A*_{βα})/2.

On the microscopic scale, we suppose that the network is made of interconnected Euler–Bernoulli beams that can have *inhomogeneous natural curvatures and cross sections*, thus extending our previous analysis on stiff networks [22]. It is worth mentioning that for three-dimensional structures, the Hamiltonian has more terms than in equation (2.1), because of the existence of two bending modes in different planes and the presence of twisting. But since their expressions are not going to appear explicitly in the following, we do not need their precise forms.

The equations of mechanical equilibrium derived from the theory of linear elasticity can be equivalently expressed in terms of a variational principle, sometimes known as the *Principle of Minimum Potential Energy* (PMPE): consider a body with volume *V* and prescribed displacements on its boundaries; the PMPE states that *among all kinematically admissible displacement fields (i.e. all continuous displacement fields satisfying the displacement constraints on the boundary), the actual displacement (i.e. the one satisfying the equations of mechanical equilibrium) is the one that makes the energy functional ε dV an absolute minimum*, where

*ε*is defined through equation (3.1). This principle is commonly used to find approximate solutions to boundary-valued problems (

*Rayleigh–Ritz*method). Here, we use it to derive rigorous upper bounds on the elastic moduli by choosing an ad hoc trial displacement field and then look for the network architecture for which the actual displacement field matches the trial displacement field.

Let us note **r**_{ij}(*l*) and *s*_{ij}(*l*) the position vector and cross-sectional area along the beam (*i*,*j*), respectively, where *l* refers to the arc-length starting at node *i* (figure 4). We also note *l*_{ij} the length of beam (*i*,*j*), *L*_{ij} the distance between nodes *i* and *j* (*L*_{ij}≤*l*_{ij}) and **e**_{ij} the unit vector along this straight line. For any network structure, one can always define a continuous displacement field such that every infinitesimal piece of the beam (*i*,*j*) undergoes the same rotation *ω*_{ij} and elongation *ϵ*_{ij}. This *trial* displacement field is (for small strains)
3.3with *ϵ*_{ij}=**e**_{ij}⋅(**u**_{j}−**u**_{i})/*L*_{ij} and *ω*_{ij}=**e**_{ij}×(**u**_{j}−**u**_{i})/*L*_{ij}. Thus, every beam (*i*,*j*) is subjected to a translation **u**_{i}, a homothety with ratio *ϵ*_{ij} and a global rotation *ω*_{ij}. To complete our definition of the trial displacement field, we impose that each node follows the macroscopic affine displacement: **u**_{i}=**A**⋅**r**_{i}, where **r**_{i} denotes the position vector of node *i*. Clearly, the trial displacement field defined this way is kinematically admissible. As each elementary piece of beam undergoes the same rotation *ω*_{ij}, there is no bending or twisting deformations, and the only contribution to the strain energy of a piece of beam with infinitesimal length d*l* is the stretching term . Moreover, for low-density structures (), the node contribution to the strain energy is negligible compared with the typical stretching energy of a beam.^{4} In this limit, the total energy associated with the trial displacement field defined above reduces to , where the summation is over all the beams that constitute the network, and is the volume of beam (*i*,*j*). Introducing the density *ϕ* as the total beam volume per unit area (two dimensions) or unit volume (three dimensions) of the network, the volumetric density of strain energy associated with the trial displacement field reads
3.4where the brackets denote the average defined, for any quantity *q*_{ij}, as: .

The elongation *ϵ*_{ij} can be rewritten in terms of the elements of **A**:
3.5where **e**_{α}⋅**e**_{ij} is the cosine of the angle between the beam (*i*,*j*) and the *α*-axis (*α*∈{*x*,*y*} for two-dimensional materials and *α*∈{*x*,*y*,*z*} for three-dimensional materials). Similarly, the components of the rotation vector *ω*_{ij}=**e**_{ij}×(**A**⋅**e**_{ij}) can be expressed in terms of the elements of **A**.

The expression of *ε*_{trial} involves averaged quantities which can be evaluated using symmetry arguments; for isotropic structures, the strain energy must be identical for any orientation of the applied strain. Thus, the bound expressions must remain invariant under rotation around the perpendicular axes or under permutation of the axis labels. After simple but lengthy calculations that we detail in appendix A, these invariance properties lead to a set of relations on the *global* structural properties of optimal networks, which can be summarized as
3.6with *α*,*β*,*γ*∈{*x*,*y*} (resp. {*x*,*y*,*z*}). These relations are referred to as *isotropy conditions* and their consequences are analysed in §3*b*. Equations (3.6) imply in particular that and (for *β*≠*γ*). Using these relations together with equations (3.4) and (3.5), the density of strain energy simplifies to (see appendix A)
3.7According to the PMPE, comparison of equations (3.2) and (3.7) yields
3.8Bounds on the elastic moduli can then be deduced from this inequality. Consider first the affine displacement field *A*_{αβ}=*γδ*_{αx}*δ*_{βy}, corresponding to a simple shear strain in the *xy*-plane. Inequality (3.8) reduces to
3.9Consider then the affine displacement field *A*_{αβ}=*γδ*_{αβ}, corresponding to a uniform compression. Inequality (3.8) yields
3.10where *κ*=*λ*+2*μ*/*d* denotes the bulk modulus of a *d*-dimensional body. Any other elastic modulus of an isotropic body is related to *κ* and *μ*. In particular, Young's modulus reads *E*=2*d*^{2}*κμ*/(2*μ*+*d*(*d*−1)*κ*). As *E* is an increasing function of *κ* and *μ*, it becomes that
3.11

Values of *μ*,*E*,*M* and *κ* for two- and three-dimensional stiffest networks are presented in table 1. We also report the values of Lamé's parameter *λ* and Poisson's ratio *ν* for this class of networks. However, it must be pointed out that they generally do not correspond to the highest possible values of *λ* and *ν*. Indeed, these two quantities are not elastic moduli, as they are not a direct measure of some strain–stress relationship. It can be noted that Poisson's ratio of such structures is below those of incompressible bodies (1 for two-dimensional bodies, for three-dimensional bodies) and is independent of the elastic properties of the beam material (*E*_{0}).

Values of the elastic moduli presented in table 1 coincide with the Hashin–Shtrikman (HS) bounds for low-density, two-dimensional structures [29,30], but are strictly lower (i.e. tighter) for three-dimensional structures [31] (table 2). However, it must be pointed out that this improvement of the HS upper-bounds is restricted to materials whose elastic phase is organized into slender objects (beams), while the HS bounds apply to any diphasic structures, including, for instance, those where the continuous phase in three-dimensional structures is assembled into two-dimensional sheets.

### (b) Mechanical conditions

We now show that inequalities (3.8)–(3.11) become strict equalities for some specific network geometries, determined by the following set of rules (along with isotropy conditions (3.6))

(a) All the beams are straight.

(b) Every beam (

*i*,*j*) has uniform cross-sectional area:*s*_{ij}(*l*)=*s*_{ij}.(c) At every junction

*i*of the two-dimensional (resp. three dimensions) network, and for all**e**_{α},**e**_{β},**e**_{γ}∈{**e**_{x},**e**_{y}} (resp. {**e**_{x},**e**_{y},**e**_{z}}), the following equality is satisfied: where denotes the set of neighbouring nodes that are connected to node*i*. Unlike the isotropy conditions (3.6), these mechanical conditions are local rules on the geometry and topology of the network.

The demonstration is straightforward: according to the PMPE, inequalities (3.8)–(3.11) become strict equalities *if and only if* the respective trial displacement fields coincide with the displacement field that satisfies the equations of mechanical equilibrium. Inspection of force and moment balances along each beam and at each junction leads to the three necessary and sufficient conditions stated above: with the trial displacement field (3.3), each infinitesimal piece of beam is subjected to an axial deformation only. Therefore, the force and local moment deriving from the trial displacement fields are **F**_{ij}(*l*)=−*E*_{0}*s*_{ij}(*l*)*ϵ*_{ij}**t**_{ij}(*l*) and **M**_{ij}(*l*)=**0**, respectively, where **t**_{ij}(*l*)=d**r**_{ij}/d*l* is the local tangent unit vector, (**F**_{ij}(*l*) and **M**_{ij}(*l*) are defined as the force and local torque exerted at position *l* by the *i*-side to the *j*-side of beam (*i*,*j*)). The force balance equation d**F**_{ij}(*l*)/d*l*=**0** along the beam (*i*,*j*) yields
3.12As **t**_{ij}(*l*) and d**t**_{ij}/d*l* are orthogonal vectors, it gives that d*s*_{ij}(*l*)/d*l*=0 and d**t**_{ij}/d*l*=**0** for all *l*∈[0,*l*_{ij}], which immediately leads to conditions (a) and (b). The moment balance equation d**M**_{ij}(*l*)/d*l*+**t**_{ij}(*l*)×**F**_{ij}(*l*)=**0** is automatically satisfied for such geometry.

Mechanical equilibrium must be satisfied at every junction *i* as well. When conditions (a) and (b) are fulfilled, the force exerted by the straight and uniform beam (*i*,*j*) on node *i* is *E*_{0}*s*_{ij}*ϵ*_{ij} **e**_{ij}. Thus, no torque is exerted on the junction, while the balance of forces yields: . This relation must hold for any orientation of the strain field. Replacing *ϵ*_{ij} by its expression (3.5) and using the same rotational invariance arguments as in §4*a* eventually yields condition (c).

## 4. Analysis of the isotropy and mechanical conditions

### (a) Isotropy conditions

Isotropy conditions (3.6) constitute a set of four (resp. 15) conditions for two-dimensional (resp. three dimensional) structures. We first check that these conditions are satisfied with a continuous and uniform angular distribution of beams (such a distribution is rather unrealistic, but it must have isotropic properties by construction). For that purpose, we reformulate the average 〈·〉 defined in equation (3.4) for the case of a continuous distribution: let us note *f*(*v*,**n**) the probability density function for a beam to have a volume *v* and an orientation along the unit vector **n**. Then, the average 〈·〉 of any quantity *X*(*v*,**n**) is
4.1where *N* is the total number of beams within the network, *V* their total volume and d*Ω* denotes the elementary angle (two dimensions) or solid angle (three dimensions) of orientation **n**. If the distribution *f* is isotropic (i.e. *f* is independent of **n**), and the quantity *X* depends only on the orientation **n** (e.g. ), equation (4.1) simplifies to
4.2The constant *c* is obtained by normalization (resp. ) for two-dimensional (resp. three dimensional) networks, where is the mean beam volume in the network. Using equation (4.2), it is straightforward to see that a uniform distribution of beams satisfies isotropy conditions (3.6). More realistic, discrete distributions of beam orientation can also satisfy the isotropy conditions. In two dimensions, conditions (3.6) can be rewritten as
4.3where *θ*_{ij} is the angle between the beam (*i*,*j*) and the *x*-axis. We can draw some general properties for networks with uniform angular distribution of beams, i.e. such that the total volume of beams with same (discrete) orientation *θ*_{m} is independent of *θ*_{m} (special cases of such networks are the periodic frameworks with identical beams and with junctions which are all similarly situated [9]). Using the complex variables , equations (4.3) reduce to
4.4where the summation is over all the distinct beam orientations. From this set of equations, it is clear that a distribution with two distinct beam orientations cannot have isotropic properties: square or rectangular lattices, for example, are known to have anisotropic elastic properties [12]. On the other hand, one can build isotropic networks with three sets of parallel beams. A simple analysis of equations (4.4) shows that the three sets of beams must be tilted from each other with equal angles (modulo *π*) of *π*/3. Hexagonal, triangular and kagome lattices are examples of such isotropic structures. The analysis of the mechanical conditions below makes it possible to determine which ones of them have the highest stiffness.

### (b) Mechanical conditions

The condition (c) produces a set of four (resp. 10) equations per node for two-dimensional (resp. three dimensional) structures, imposing severe restrictions on the geometry and topology of a junction. In this section, we inspect the solutions to this set of equations for two-dimensional networks. We suppose that conditions (a) and (b) are satisfied (straight and uniform beams). Moreover, we assume that beams have the same cross section. Let *z*_{i} be the connectivity of node *i*, and *θ*_{ij} the angle between the beam (*i*,*j*) and the *x*-axis. The vectorial condition (c) can be rewritten as four trigonometric relations
4.5

Clearly, there is no solution to this set of equations for a monovalent node (*z*_{i}=1): all the beams of an optimal network are connected at their both ends and can contribute to the storage of elastic energy. For a divalent node (*z*_{i}=2), the only solutions are the trivial configurations *θ*_{i2}=*θ*_{i1}+*π*. To inspect the possible configurations of nodes with higher connectivity, we rewrite equations (4.5) as
4.6where we introduced the complex variables . One can now show that there is no solution to equations (4.6) for a trivalent node (*z*_{i}=3): eliminating *y*_{i3} from these equations yields *y*_{i1}*y*_{i2}(*y*_{i1}+*y*_{i2})=0; either *y*_{i1}, *y*_{i2} or *y*_{i3} is zero, which is not compatible with the normalization constraint |*y*_{ij}|=1. For quadrivalent nodes (*z*_{i}=4), eliminating *y*_{i4} from equations (4.6) yields: (*y*_{i1}+*y*_{i2})(*y*_{i2}+*y*_{i3})(*y*_{i1}+*y*_{i3})=0. Thus, the only possible configurations for such nodes are *θ*_{i3}=*θ*_{i1}+*π*, *θ*_{i4}=*θ*_{i2}+*π* and all other subscript permutations, i.e. beams must be collinear in pairs. In agreement with Maxwell's criterion [9,25–27], it is known that two-dimensional stiff networks must have a node connectivity more than or equal to 4. But in addition, our analysis specifies what must be the geometry of the junctions (figure 5). In the special case where nodes are all similarly situated, Deshpande *et al.* [9] have shown that the minimal node valency is 6. Here, the study is not restricted to such structures, hence there is no contradiction between these two results.

It is possible, in principle, to find the solutions to equations (4.5) for nodes with higher valencies. The solutions are certainly more tedious to find. However, it can be mentioned that for junctions with an even number of adjoining beams, the geometry with beams collinear in pairs is always a solution to equations (4.5).

As a consequence of our analysis, it can be concluded that two-dimensional optimal networks must contain triangular cells (but the converse is generally not true, e.g. see Delaunay network in figure 3*e*); the number of cells *C*, the number of sides (beams) *N* and the number of junctions *J* of a two-dimensional network are related by Euler's formula: *C*−*N*+*J*=1 [9]. Therefore, the mean coordination number and the mean number of sides per cell satisfy: . Thus, in optimal networks, and . Clearly, optimal networks contain triangles when . The marginal case () must be inspected separately: either each cell has exactly four sides; in order to satisfy the mechanical condition, every vertex must be a fourfold junction with beams collinear in pairs, but we have seen that a structure with only two distinct beam orientations does not satisfy the isotropy conditions; or cells do not have all the same number of sides; then the constraint imposes that such a network must contain triangles.

A striking feature of the stiffest networks is that, under uniform loading, the displacement field is affine down to the microscopic scale and thus coincides with the macroscopic displacement field : as every beam must be straight, one has **r**_{ij}(*l*)=*l* **e**_{ij}. Thus, according to equation (3.3), the displacement field at location **r**=**r**_{i}+**r**_{ij}(*l*) simplifies to .

Finally, it is worth mentioning that isotropic networks with optimal elastic properties also have optimal transport properties. Indeed, isotropy conditions (3.6) encompass the isotropy conditions for transport properties. Furthermore, using the identity , it is easy to show that condition (c) implies , which is the condition required at every junction of a network (together with conditions (a) and (b)) to have optimal *transport* properties [32,33].

### (c) Comparison with numerical results

The analysis of the isotropy and mechanical conditions made above sheds light on the numerical results reported in §2: the five simulated networks satisfy isotropy conditions (4.3), but only the kagome and triangular networks satisfy the mechanical conditions. It is noteworthy that the kagome lattice is one of the stiffest networks, in spite of its large number of floppy modes [34–36]. However, small defects in its structure will dramatically affect its mechanical response [37].

The mechanical conditions involve both geometry (orientation and cross-sectional areas of the beams) and topology (connectivity) of each junction. That explains why two elastic networks with same connectivity, similar to triangular and Delaunay lattices, can have very different macroscopic responses. Moreover, our analysis has revealed the strong correlation between the stiffness of a network and the affinity of its displacement field. This is also what we observe in the simulations: the two networks with the highest elastic moduli—kagome and triangular lattices—are also the only ones with pure affine deformations. It must be emphasized that stretch-dominated networks do not all deform in an affine way. This is illustrated in our simulations with the Delaunay lattice: this structure deforms exclusively through the stretching of its members (figure 2*e*), but does not deform affinely (see footnote 3 and figure 3*e*). Indeed, the condition (c) is not satisfied in a Delaunay lattice.

## 5. Examples of optimal structures

In §§3 and 4, we have derived the structural conditions that an isotropic network must fulfil to have the highest elastic moduli for a given density and Young's modulus of beam material. Nevertheless, one might wonder whether such networks do exist, that is, whether one can effectively build networks that satisfy all these conditions simultaneously. The answer is obviously yes in two dimensions, as we have already provided two examples of optimal networks (kagome and triangular lattices). In this section, we show how our theoretical analysis can be used to identify other optimal networks in both two and three dimensions.

As isotropy and mechanical conditions are explicit, it is easy to use them to find new optimal networks. Figure 6 shows three additional optimal networks with two-dimensional periodic arrangements, here made of non-identical beams. The first one, often referred to as the *Union Jack lattice*, is made with two different kinds of beams. Beams are collinear in pairs at every junction, so the mechanical conditions are satisfied. The ratio of beam cross sections is then adjusted to fulfil two-dimensional isotropy conditions (4.3): the only restrictive condition is ,^{5} yielding 2(*l*_{1}*s*_{1}−*l*_{2}*s*_{2})=0 (see notations on figure 6*a*). As, from geometry, , it is known that .

The second network (figure 6*b*) is built with three different kinds of beams. Here, beams are not collinear in pairs at every junction. The isotropy conditions are satisfied. The angle *α* and the two cross-section ratios *r*_{2}=*s*_{2}/*s*_{1} and *r*_{3}=*s*_{3}/*s*_{1} (see notations in figure 6*b*) are then adjusted to satisfy the remaining isotropy and mechanical conditions. We obtain
5.1This set of equations can be solved numerically; we obtain: *α*≃0.927 rad ≃53.15^{°}, *r*_{2}≃0.506, *r*_{3}=0.805. It can be noted that, like the kagome lattice, this optimal structure is not fully triangulated.

The third network (figure 6*c*), based on the kisrhombille tiling, is made of three different kinds of beams. From geometry, and *l*_{3}=2*l*_{1}. As before, the isotropy conditions are satisfied by construction and the cross-section ratios *r*_{2} and *r*_{3} are adjusted to satisfy the mechanical conditions. One obtains that the optimality is reached when *r*_{3}=1, whatever the value of *r*_{2}.

We also performed simulations of these networks and checked that the numerical values of their elastic moduli are in very good agreement with the expected values: *μ*/(*E*_{0}*ϕ*)≃0.125, *M*/(*E*_{0}*ϕ*)≃0.375 for the network of figure 6*a*, *μ*/(*E*_{0}*ϕ*)≃0.124, *M*/(*E*_{0}*ϕ*)≃0.374 for the network of figure 6*b* and *μ*/(*E*_{0}*ϕ*)≃0.125, *M*/(*E*_{0}*ϕ*)≃0.375 for the network of figure 6*c*.

The three-dimensional case is certainly more complex. As a matter of fact, we were not able to identify any periodic and isotropic optimal structure made with one single kind of beam. The existence of such a structure is still an open problem though, as we were also unable to show that it does not exist in general. Nevertheless, stiffest structures built with two (or more) types of beams do exist. For instance, the structure depicted in figure 7 satisfies the mechanical conditions and can be repeated periodically: the beams join the centres of the faces of a Kelvin cell (which is known for tiling the space [38]). The length ratio of the two kinds of beams is imposed by the Kelvin cell geometry: *l*_{2}/*l*_{1}= . We then adjust the section ratio for isotropy conditions (3.6) to be satisfied. We obtain: . To our knowledge, this optimal (and isotropic) network with periodic arrangement has never been published before and represents one of the main results of this study.

## 6. Conclusion

In summary, we showed the existence of a class of isotropic networks having the highest possible values of elastic moduli for a given density. We established a convenient set of conditions that allow identification of these optimal networks. The elastic moduli of these networks are also derived and can be simply expressed in terms of Young's modulus of the beam material and the relative density of the structure. Examples of two- and three-dimensional optimal structures with periodic arrangements are also given. These results may be of interest for structural applications as well as for our understanding of biological systems.

It must be emphasized that this study is limited to small deformations only. Although the microstructure may give rise to an increased stiffness in this regime, the same microstructure may be an origin for failure owing to different physical mechanisms upon a critical load [39–44]. Therefore, the stability properties of the optimal structures presented in this study are an important issue that needs to be addressed in a future work.

## Appendix A. Isotropy conditions

Introducing expression (3.5) of the elongation *ϵ*_{ij} in equation (3.4) yields
A1where *a*_{αβ}=(*A*_{αβ}+*A*_{βα})/2. For an isotropic network, expression (A 1) must be invariant by permutation or inversion of axes for any applied strain. These invariance properties lead to restrictions on the structure of the network and simplify the expression of *ε*_{trial}.

First, we note that the last sum in equation (A 1) vanishes for two-dimensional networks, as *α*, *β* and *γ* cannot be all distinct. We show that this sum also vanishes for three-dimensional isotropic networks, by analysing the strain field defined as *a*_{xx}=*a*_{0}, *a*_{yz}=*a*_{zy}=*b*_{0} and the other components *a*_{ij}=0. Thus, equation (A 1) reduces to . Inverting the *y*- or *z*-axes changes the sign of the last term only. But as *ε*_{trial} must remain unchanged under such an operation, one necessarily has , and more generally, by permutation of the axes: with *α*, *β* and *γ* all distinct. We then choose the strain field defined as *a*_{xx}=*a*_{0}, *a*_{xy}=*a*_{yx}=*b*_{0} and the other components *a*_{ij}=0 to show that the third sum in equation (A 1) also cancels out: equation (A 1) reduces to . Inverting the *x*-axis transforms into . As the energy must remain unchanged, one must have . Performing the same procedure for any pair of axes *α* and *β*≠*α* eventually leads to the conditions .

We now choose the strain field defined as *a*_{αβ}=*a*_{0}*δ*_{αx}*δ*_{βx} (this is a particular case of the field defined above). Equation (A 1) then reduces to . As *ε*_{trial} must remain unchanged by permutation of the axes, the quantity must be independent of the *α*-axis. Similarly, using the strain field *a*_{αβ}=*a*_{0}*δ*_{αx}*δ*_{βy}, it follows that the quantity is the same for all *α* and *β*≠*α*.

Using the identity , the first relation between and is easily obtained: . We obtain the second relation by using the invariance by rotation of the axes. For instance, a 45^{°} rotation around the *z*-axis ‘transforms’ into . The equality then leads to . Thus, . Finally, rearranging the terms in (A 1) leads to expression (3.7).

## Appendix B. Finite-size effects

The finite size of the simulated networks (see §2) may alter the numerical values obtained for the two independent elastic moduli *G* and *M*. In order to assess the importance of these *finite-size effects*, we performed simulation runs with varying network size *L*. The results of this procedure are shown in figure 8 for the triangular and Delaunay networks: for the triangular network (figure 8*a*), the numerical values of *μ*/(*E*_{0}*ϕ*) and *M*/(*E*_{0}*ϕ*) converge towards the respective theoretical upper-bound values and (the oscillations around these asymptotic values are caused by the periodicity mismatch between the *x* and *y* directions for the triangular lattice). In the simulations presented in §2, the ratio of the network size over the mean beam length is *L*/*l*≃40. Thus, the finite-size effects do not alter the obtained values more than 0.4%.

For the Delaunay network (figure 8*b*), we observe some differences in the evolution of the two elastic moduli, whereas *μ*/(*E*_{0}*ϕ*) seems quite unaffected by the finite network size, the quantity *M*/(*E*_{0}*ϕ*) shows a clear increase with *L*/*l*. Actually, *M*/(*E*_{0}*ϕ*) varies almost linearly with *l*/*L*. A linear regression *M*/(*E*_{0}*ϕ*)=*a*+*b* *l*/*L* yields *a*=0.2914±0.0002 and *b*=−0.478±0.007. The fit parameter *a* represents an estimation of the ratio for an infinite network; this is the value of this parameter which is reported in §2. We also performed a linear fit for *μ*/(*E*_{0}*ϕ*) and obtained *a*=9.143×10^{−2}±4×10^{−2} and *b*=−1.5×10^{−3}±1.3×10^{−3}. It can be noted that for the two elastic moduli, the value of the fit parameter *a* is clearly under the respective bound ( for *M*/(*E*_{0}*ϕ*) and for *μ*/(*E*_{0}*ϕ*)). Besides this finite-size effect, there is also a statistical error on the numerical values of the elastic moduli, owing to the random nature of the Delaunay network. The error bars in figure 8*b* represent the standard deviations computed over several (more than 100) realizations with the same value of *L*/*l*.

## Footnotes

↵1 We impose a uniform displacement on the top boundary, and a zero displacement on the bottom boundary. Beams are attached with ‘free hinges’ to these boundaries.

↵2 The periodic boundary conditions impose a zero lateral strain under uniaxial load. That is why the simulations under uniaxial load gives the longitudinal modulus

*M*instead of Young's modulus*E*.↵3 As the network deforms primarily through the beam stretching mode, the displacement of point

**r**belonging to beam (*i*,*j*) is**u**_{ij}(**r**)=**A**_{ij}·**r**(piecewise affine function).↵4 It can be noted that no assumption is made on the relative importance of

*κ*_{n}and*κ*_{b}.↵5 It is worth mentioning that this lattice has optimal conductivity for any value of the cross-section ratio.

- Received September 13, 2013.
- Accepted January 17, 2014.

- © 2014 The Author(s) Published by the Royal Society. All rights reserved.