Customization of two-dimensional extremal materials

(cid:1)


Introduction
Metamaterial, a kind of man-made materials, possesses exotic effective properties from structure rather than directly from its composition [1].Building upon the success of electromagnetic and acoustic metamaterials, researchers working on mechanical metamaterials strive at obtaining extraordinary or extreme elasticity tensors and mass-density tensors to thereby mold static stress fields or the flow of longitudinal/transverse elastic vibrations in unprecedented ways [2].The research domain of mechanical metamaterials normally includes auxetics with negative Poisson's ratio [3], light-weight metamaterials with strong stiffness [4], metamaterials with negative parameters (mass densities and elastic moduli) [5], as well as the pentamode metamaterials which are virtually a special kind of extremal materials [6].
In continuum mechanics, the fourth-rank elasticity tensor c ijkl is used to describe the mechanical properties of a linear-elastic material, also known as the Hookean solid.Traditional elastic continuum's elasticity tensor is symmetric and positive definite.The minor symmetry (i.e., c ijkl ¼ c jikl ¼ c ijlk ) is a direct consequence of the conservation of angular momentum, while the major symmetry (i.e., c ijkl ¼ c klij ) holds due to the existence of a strain energy function in a conservative system.The positive definiteness property results from the basic thermodynamics, i.e., the strain energy must be a positive quantity that vanishes only in the undeformed state [7].In other words, in order to deform a material the strain energy should be expended.A direct result of the positive definiteness property is that all the eigenvalues of the elasticity tensor c ijkl are positive.
Milton [6] generalized the concept of linear-elastic material to include the ones whose microstructure can contain judiciouslydesigned mechanisms such that no elastic strain energy is expended to deform the microstructure in some specified deformation modes (termed as the soft modes in this paper).Within this generalization, the eigenvalues of the elasticity tensor can approach zero, and the elasticity tensor c ijkl is now only positive semi-definite.In that case, the eigenvector of the elasticity tensor corresponding to the zero eigenvalue characterizes the soft mode.Similarly, the eigenvectors corresponding to the non-zero eigenvalues is called hard modes here.In this regime, linear-elastic materials (in three-dimensional case) can be classified respectively into nullmode, unimode, bimode, trimode, quadramode and pentamode materials depending on the number of zero eigenvalues of the elasticity tensor c ijkl .In two-dimensional case, however, only the nullmode, unimode and bimode materials exist since the elastic matrix is of dimension 3 Â 3. The elasticity tensor of the nullmode, unimode, and bimode materials have no, one, and two zero eigenvalues, respectively.
The elasticity tensor of pentamode material has only one nonvanishing eigenvalue.Isotropic pentamode meta-materials can thus only support hydrostatic stress and is soft for any other stress states.Due to their water-like wave properties and solid nature, such kind of materials show great potential for underwater acoustic applications.Since their first fabrication by Kadic et al. [8], they have been widely investigated in recent years [9][10][11][12].Acoustic simulation and experimental testing results based on planar pentamode with honeycomb-like microstructure indicate that the designed pentamode material mimics water's acoustic properties over a wide frequency range [9,12].Underwater acoustic cloak which is fascinating for engineering applications can be realized by annular graded solid pentamode microstructures [9], later experiment has been conducted to demonstrate the feasibility of the underwater pentamode cloak [13].It is also found that the metamaterials made of solid constituents are inherently imperfect since small shear resistance is not only unavoidable but also necessary for the stability of structure [9].The shear resistance will incur the shear resonance, which can be suppressed by adding appropriate damping [10].A metasurface with graded honeycomb pentamode unit cell can convert underwater cylindrical wave into plane wave [11].Since only quasi-static material properties are considered, such metasurface is broadband and allows a high transmission ratio due to the matched impedance between the pentamode material and water.Excellent wavefront manipulation such as anomalous refraction, generation of non-diffracting Bessel beam, and sub-wavelength flat focusing can be realized by employing a gradient velocity to redirect refracted waves and pentamode metamaterials to improve impedance matching between the metasurface and the background medium [14].An inhomogeneous acoustic metamaterial lens based on spatial variation of refractive index for broadband focusing of underwater sound is also reported [15].The index gradient follows a modified hyperbolic secant profile designed to reduce aberration and suppress side lobes.The gradient index lens is comprised of pentamode microstructures with tunable quasi-static bulk modulus and mass density.This design approach has potential applications in medical ultrasound imaging and underwater acoustic communications.
Apart from pentamode materials, there are also some works devoted to the study on other kind of extremal materials.Auxetic materials can be used to design unimode materials, of which the only easy deformation mode is dilation [16].Unimode metamaterials made from rotating rigid triangles exhibit the anomalous property of negative linear compressibility along loading direction in which the Poisson's ratio is larger than 1 [17].Milton constructed a large family of non-linear bimode materials using rigid bars, pivots and actuators [18].Recently, the wave properties of different quadramode materials based on material symmetry are also investigated [19].A quadramode based out-of-plane shear wave polarizer is validated both theoretically and numerically, which effectively prevents the mode conversion at the interface between the solid and liquid phases [19].
So far the most widely studied pentamode microstructures are diamond-like ones in 3D case or honeycomb-like ones in 2D case, both of which were suggested by Milton using clever arrangements of mechanisms in the unit cell [6].As an alternative, the topology optimization approach can be used as a powerful tool to explore other potential pentamode microstructures in a much more effective and efficient manner [20].Topology optimization attempts to find optimal layout of (macroscopic or microscopic) structures under given objective functions and constraints.In the cases of microstructures, many new and innovative microstructures can be found by means of topology optimization with appropriately formulated objective functions (e.g., maximizing the ratio between the bulk and shear modulus [21], maximizing the broadband wave manipulation ability [22], achieving an ultra-lightweight highperformance micro-lattice [23]).These works all took advantage of the symmetry of the unit cell, which reduced the number of design variables for high computation efficiency but also limited the region of search space.Therefore, in this study the symmetry of the unit cell will not be taken into consideration in advance so as to explore as many microstructures as possible.In addition, the eigenvectors of elasticity matrix have not been paid enough attention in all the current works, but they do have a huge impact on the subsequent application of extremal materials.For example, the positive and negative honeycomb are both bimode materials.The biggest difference between them is that the hard mode of the former is ½1; 1; 0 > while the latter is ½1; À1; 0 > .Comparing to the positive honeycomb, the negative ones can be applied to splitting the mode of waves [24].Therefore, in terms of applications of extremal materials, the design method should take the eigenvectors into consideration.
Although the idea of extremal metamaterials with zero eigenvalues has been proposed for over 25 years since Milton's landmark work [6] and various microstructures have been proposed from that on, there is still a lack of a systematic and general approach to design the microstructures of the extremal metamate-rials.In particular, no framework is available in the literature to design an extremal metamaterial with specified number of zero eigenvalues and with prescribed soft or hard mode at the same time.This study is intended to fill this gap to establish a systematic framework to design the extremal metamaterials.
In this paper, we follow the idea of Milton's early work that tries to construct general extremal metamaterials with arbitrary predetermined soft or hard mode [6], but here we make use of the structural optimization, or more specifically, the topology and shape optimization methods.The remaining parts of this paper is organized as follows.In Section 2, the properties of extremal materials are briefly reviewed, where the Kelvin notation is used due to its convenience in matrix expressions.The overall optimization framework and detailed optimization models of different kinds of extremal metamaterials are elaborated in Section 3. The numerical homogenization method of the extremal metamaterials are given in Section 4. In order to make the whole optimization workflow effective and robust, there are quite a lot of implementation strategies that should be taken into consideration when solving the optimization problem.All these implementation details are discussed in Section 5.In Section 6 we discussed in detail the rigid rotation process of the unit cell.Relevant properties of the rotation process are also covered, based on which it can be found that all the rotated strain tensor can be classified by their first principal invariants, thus greatly reducing the searching scope of the numerical optimization process.In Section 7, the entire optimization flow is demonstrated by a series of numerical examples.We successively show in this section how to design two-dimensional extremal metamaterials without and with constraints on the soft modes and/or hard modes.Finally, the concluding remarks are given in Section 8.

A brief review on the extremal materials and design target
In this section, we briefly illustrate some basic properties of the two dimensional extremal materials.The Kelvin notation1 [27][28][29] will be used to express the second-rank stress and strain tensors as well as the fourth-rank elasticity tensor in matrix form, which is more computationally tractable.Shortly, the sufficient and necessary conditions for all types of two dimensional extremal materials are derived based on the characteristic polynomial of the elasticity tensor.
For a linear-elastic material, the constitutive relation between strain and stress is given by c : e ¼ r; ð1Þ where c ¼ c ijkl e i e j e k e l is the well-known fourth-rank elasticity tensor, e i is the Cartesian basis vector; r and e are the stress and strain second-rank tensors, respectively.
Let k and v denote the eigenvalue and second-rank eigentensor of c, then it holds that, Comparing Eqs. ( 1) and (2) we know that when ðk; vÞ composes an eigen-pair of c, the second-rank tensors v and kv indeed make up a pair of strain and stress of the linear material.In this sense, v and kv can be called as the easy strain and easy stress when k ¼ 0, and hard strain and supporting stress when k > 0 [6].In this paper, however, we simply call v the soft mode if k ¼ 0 and the hard mode if k > 0.
It is not easy to work directly with the two-dimensional fourthrank elasticity tensor c, so it is desirable to convert the tensor form into a more numerically-tractable matrix form.There are various ways to achieve this goal, of which the most common one is the Voigt notation.However, Mehrabadi and Cowin [29] demonstrated that the Voigt notation is non-tensorial in the sense that the eigenvalues of the elasticity tensor c are in general not the eigenvalues of its Voigt notation counterpart C ½ Voigt .The inverse and contraction of c cannot be computed by the inverse and multiplication of C ½ Voigt , either.
In this paper we choose to express the elasticity tensor with the less-famous Kelvin notation, which is originated from Kelvin's seminal work [27,28] on the decomposition of the elasticity tensor.Using Kelvin notation, Eq. ( 1) can be rewritten as a matrix form2 , r where the vector representations of the stress and strain tensors are, r f g ¼ Note that we use the braces Á f g and brackets ½Á to emphasize that they are the vector and matrix representations of tensors, respectively.
It has been demonstrated that the components of the elasticity matrix C ½ in 2D indeed correspond to a three-dimensional secondrank tensor [29], implying that the elasticity tensor c ijkl can be fully represented by the matrix C ½ .This way the eigenvalues, inverse, and contraction of the elasticity tensor can be directly computed by using C ½ without resorting to the original fourth-rank tensor c ijkl [30].For example, it is easy to prove that the fourth-rank tensor c ijkl and the matrix C ½ have exactly the same eigenvalues.Another benefit from the Kelvin notation is that the eigenvalues of C ½ do not change under coordinate system transformations, which is consistent with the tensorial nature of c ijkl .
Using Kelvin notation, the eigenvalue problem Eq. ( 2) can also be written in the matrix form, where C ½ and v f g are the Kelvin notations of c and v, respectively.
Eq. ( 6) can be equivalently written as where I ½ is the identity matrix.In order to find the eigenvalue k corresponding to a nontrivial eigenvector v, a characteristic equation can be formulated from Eq. ( 7), where the coefficients in Eq. ( 8) can be expressed as the polynomials of c ijkl , The Einstein's summation rule is applied throughout this paper; tr (Á) and det(Á) stand for the trace and determinant of a square matrix, respectively.The derivation of Eq. ( 8) is by no means trivial and we obtain it following the Faddeev-LeVerrier algorithm [31].
Let k i (i ¼ 1; 2; 3) denote the roots of the cubic Eq. ( 8), then by using the fundamental theorem of algebra, Eq. ( 8) can be rewritten as Applying Vieta's formulas to Eq. ( 8) and keeping in mind that all the eigenvalues k i are non-negative (due to the basic thermodynamic restrictions), the sign of each coefficient can be determined, By using the fundamental theorem of algebra, the condition of such two-dimensional extremal materials can be given by, It should be noted here that Eqs. ( 12) and ( 13) indeed show the sufficient and necessary conditions of the two-dimensional unimode (with one zero eigenvalue) and bimode materials (with two zero eigenvalues), respectively.For any linear-elastic material whose elasticity matrix is known, the coefficients C 0 ; C 1 ; C 2 can be computed by Eq. ( 9).Then the material is definitely a unimode material if Eq. ( 12) holds true.On the other hand, the coefficients C 0 ; C 1 ; C 2 computed from a unimode material would also definitely satisfy Eq. ( 12).For bimode materials, similar conclusions can be drawn.
Besides the prescribed number of zero-eigenvalues of elasticity matrix, the eigenvectors (soft or hard modes) are taken into consideration in our design target of two-dimensional extremal materials.The performance index Err is used here to evaluate the degree of matching between the calculated soft or hard modes of the optimized microstructures and the prescribed ones.In detail, the performance index Err is given by where C ½ is the elasticity matrix of the optimized result, and v d f g is the prescribed eigenvector( v d f g ¼ v s f g for unimode and f g for bimode).k is the eigenvalue in the eigenpair ðk; v d f gÞ.

Optimization framework
In this paper we concentrate on the periodic medium, i.e., the composites with periodic unit cells (See Fig. 1).The periodicity feature implies that under macroscopic loads all physical quantities (e.g., displacement, strain, stress, etc.) should have the same period with the geometry.Therefore, it is possible to take advantage of mathematical tools to calculate quantitatively the exact values of the effective characteristics.In Fig. 1 the black lines represent the bars while the red dots denote the joints.The joints here indeed simulate the pivot joints, thus no moments can be transmitted by these joints.This way the bars can only support axial loads and the periodic unit cells can be designed to have some intrinsic mechanisms.
The optimization framework here essentially refers to the twodimensional unit cell, which contains two main optimization steps, i.e., the topology optimization step and the shape optimization step.The first step aims to obtain the layout of microstructure while the second step attempts to tune the nodal coordinates of the microstructure so that it can better match the soft/hard mode with prescribed ones.Very often, a series of micro-structures of extremal materials can be obtained by simply using the first step, i.e., the topology optimization.Still, the second step can act as a safeguard that the proposed approach can design extremal materials with arbitrary predefined soft/hard modes.

A general introduction on the design flow
Since only periodic medium is considered, we focus on the design of the unit cell of extremal material.The flowchart and diagram of the design flow of the extremal materials considering the prescribed soft or hard mode are shown in Fig. 2.There are two stages in the design flow, and the task of each stage is outlined as follows.
(1) In Step 1, we use topology optimization approach to design the microstructures of the extremal materials.As illustrated in Fig. 3, we use the ground structure as the starting point of our optimization.In the ground structure [20], every two nodes are connected by a pin-jointed truss element.In Fig. 3(b), the nodes are marked by red dots while the bars are marked by black lines.Due to the periodicity of the unit cells, there is no bars on the left and bottom edges, otherwise the cross-sectional areas of the bars in the microstructure would not be all the same.In Fig. 3(a) we list the labels of the nodes.For such a ground structure with n Â n nodes, the number of the possible bars is In the light of the topology optimization approaches, the cross-sectional area of the i-th bar element is attached to a design variable f i 2 ½0; 1.As the conventions of structural topology optimization approaches, the i-th bar element will be retained in the final optimized design if f i ¼ 1 while discarded if f i ¼ 0. It is well-known that rigid rotation of the unit cell leads to transformation of all the entries of the elasticity tensor.Under this transformation, the eigenvalues k i would not change, but the eigenvectors v i would change accordingly.Thus, a specific layout of unit cell indeed corresponds to enormous number of eigenvectors.In order to take advantage of this property, we add the rotation angle h around z axis into the optimization formulation in the first step (cf.Fig. 2b) and and Section 3.2).More details of the rotation process can be found in Section 6.In Fig. 2 (b) we show the evolution of one unit cell during the optimization.As can be seen from Fig. 2 (b), after the first step both the layout and the rotation angle are altered.
(2) In Step 2, we intend to find a more precise solution based on the result from the first step by shape optimization.A postprocessing process is executed at the end of Step 1 so that a clear and concise interpretation of the layout of the unit cell can be available.Then the locations of all the nodes in the unit cell are chosen as the design variables in this step.The rotation angle keeps unchanged during the this step to avoid the optimization problem to be too complicated.The lattice vectors of the unit cell will also stay unchanged during the this step, avoiding extra computation cost.
As shown in Fig. 2  it means that the precision is acceptable for both Step 1 and Step 2, so can be output.In Fig. 2 (a), there are two lines linking with the box End.One represents that the algorithm find a good enough solution by utilizing both Step 1 and Step 2, while the other one reveals that the solution can be obtained with only Step 1.Meanwhile, there is a counter g 2 in the second loop.If g 2 is larger than a predefined number (say 50), the algorithm will enter back to Step 1 and such process can be regarded as the loop between Step 1 and 2. The loop between Step 1 and 2 is necessary because the contin-uous failure of Step 2 implies that the results obtained in Step 1 is not good enough.
The implementation details of all these steps will be elaborated in Section 5.In the following subsections we carry out the mathematical models of Step 1 and Step 2.

Optimization model in Step 1
In this paper, we want to design the extremal materials with prescribed number of soft modes.The only soft or hard mode of  the designed extremal materials should be consistent with the prescribed ones.On the other hand, the topology optimization methods use continuous variable f i 2 ½0; 1 to show the existence of i-th bar (f i ¼ 0 reveals that i-th bar will be discarded in the optimized design, while f i ¼ 1 shows that it will be retained).Sometimes (e.g. when f i ¼ 0:2) it is vague to justify whether a bar should be retained or discarded, so the optimization method should be carefully formulated to get rid of such occasions.Another issue needs to be considered is the intersection and overlapping of different bars.In the ground structure the intersection and overlapping occur everywhere, but in the optimized design these phenomena should be avoided.In summary, the optimization results should fulfill the following requirements, The number of soft modes of deformation must be precisely equal to the given ones.For example, if we are to design the unimode material, then the optimization results should have only one soft mode of deformation.The number of zero eigenvalue of the elasticity matrix must be exactly 1.The optimization results should be clear.That is to say, we want all the design variables to approach either 1 or 0 when the optimization process is completed.If there are many design variables with intermediate values between ð0; 1Þ, the design results are hard to interpret.With manufacturing constraints considered, intersection or overlapping among the bars in the unit cell is not allowed.

Suppose the elasticity matrix of the optimized results before rotation is C
½ and rotation angle is h, the elasticity matrix of the rotated one ½ Ĉ can be expressed as: where NðhÞ ½ is the rotation transformation matrix which is only dependent on the rotation angle h.The detailed expressions of NðhÞ ½ will be elaborated in Section 6. Two-dimensional unimode materials have only one zero eigenvalue.Therefore, we aim at designing such kind of extremal material with the desired soft mode v s f g.If v s f g; 0 ð Þ composes the eigenpair of ½ Ĉ which has only one zero eigenvalue, we have Therefore, the performance index for unimode materials Err s can be stated as: and the mathematical model for the optimization of twodimensional unimode materials can be formally given as where f i is design variables of topology optimization which determines the existence of the i-th bar, N bar is the total number of design variables.C i ði ¼ 0; 1; 2Þ is the coefficient of the characteristic equation (cf.Eq. ( 8)).
The objective function in Eq. ( 18) warrants some special attention.The first part of the objective function is used to control the grey elements, which is an annoying problem in the topology optimization methods.The judiciously selected function , the function obtains the maximum when f i ¼ 0:5 while the minimum when f i ¼ 0 or 1.In other words, the first part of the objective function drives the design variable to get away from the intermediate value ð$ 0:5Þ.The second term is the penalty function which can effectively avoid the existence of the intersection and overlapping bars in the unit cell.M 2 0; 1 f g N bar ÂN bar is a logic matrix showing the intersection and overlapping status of every pairs of bars.Generation of the logic matrix M ½ will be given in the Section 5.The coefficients c 1 ¼ 10 À5 ; c 2 ¼ 1 are use to tune the strength of the two penalty terms.The third part of the objective function aims at minimizing the coefficient of the quadratic item C 2 (cf.Eqs. ( 8)), which is equivalent to drive C 2 2 as large as possible.Such condition is not mandatory in finding the unimode materials but it can enhance the optimization performance.
Although the third term in the objective function means to drive C 2 2 as large as possible, there are occasions where the opti- mization algorithm gets stuck in unwanted local optimum, i.e., the elasticity matrix is nearly a null matrix.In order to get rid of the meaningless solutions, we add the first constraint (cf.Eq. ( 18)) to emphasize that the sum of the diagonal of ½ Ĉ should not be smaller than a predefined positive value.Note that C 2 is the negative of the sum of the diagonal of ½ Ĉ, so it should not be greater than a predefined negative value.The first three constraints are the necessary and sufficient conditions of the unimode materials (cf.Eq. ( 12)).
The last constraint is added to require that the designed extremal material should behave softly in a prescribed soft mode v s f g.
Numerically it is not practical to expect that ½ Ĉ v s f g 2 exactly equals to 0. Instead, the boundary threshold g 1 is introduced.
Two-dimensional bimode materials have only one nonzero eigenvalue.Let v h f g denotes the desired hard mode (supporting stress) which is also the eigenvector of the elasticity matrix ½ Ĉ that has only one non-zero eigenvalue, we have Therefore, the performance index for bimode Err h can be stated as: and the mathematical model for the optimization of twodimensional bimode materials can be formally given as The objective function for the bimode materials are exactly the same as that of the unimode ones.The first constraint in Eq. ( 21) is also the same as that in Eq. ( 18).The second and third constraints are the necessary and sufficient conditions of bimode materials.The final inequality constraint guarantees that the hard mode of the bimode material follows the prescribed hard mode.k max is the only non-zero eigenvalue (also the maximal one).Due to the Vieta's formula, we know that So k max can be approximated by if the other two eigenvalues are close enough to 0. The approximation in Eq. ( 23) makes it easy to compute the sensitivity information of the constraint function.
Although the penalty terms of the grey elements and intersection or overlapping are list in the objective function, the solution may still encounter such unwanted phenomenon.In addition, there may be redundant bars (i.e., bars that have no stress when the overall materials is subjected to external loads) in the solution optimized by Step 1.Therefore, postprocessing should be applied to deal with such situations, and the details of the postprocessing will be introduced in Section 5.

Optimization model in Step 2
As illustrated in the previous subsection, the design of extremal materials is highly nonlinear in that there are many material, geometric and algebraic issues to be taken into consideration.There-fore, only a small part of extremal materials with predefined soft or hard mode can be found within merely Step 1 (since only finite number of nodes can be present in the design domain).In Step 1 the nodes stays fixed during the optimization process.Hence, shape optimization model is developed in this subsection to alter the locations of the nodes such that the configurations of the designed results by Step 1 can be modified to improve the accuracy of the soft or hard mode.The design variables n f g in this stage are the positions of the nodes in the unit cell.
Before introducing the mathematical model, the nodes in the unit cell need to be classified.As shown in Fig. 4, the 16 nodes can be divided into three types.The four red nodes located at the corners belong to the first type, and are called the corner nodes.The blue nodes lying at the border of the unit cell belong to the second kind, and are called the border nodes.The remaining green ones belong to the third type, and are called the internal nodes.
In order to maintain the lattice vectors, the corner nodes (cf. the red circles in Fig. 4) are fixed and not allowed to move.Hence the movable nodes are the border and internal nodes.The movable range of every internal node is limited in a square with length 1.As an example, the movable range of the internal node 6 is shown as the pink square in Fig. 4. We note here that the range of each movable nodes is arranged in this way to ensure that the moved configuration would not suffer from intersection or overlapping and at the same time can have sufficient searching space.Suppose that the coordinates of the i-th internal node before and after moving can be represented as ðx i ; y i Þ and ðn i x ; n i y Þ, the movable range requires that n i x À x i 0:5 and n i y À y i 0:5.
As for border nodes (cf. the blue circles in Fig. 4), the movement range of them is larger than that of internal nodes.Let a ¼ 3 denote the length of the basic lattice vector, and d ¼ 0:1 the minimum allowable distance between different nodes.The movable range of every border node is limited in a rectangle with size of 1 Â ða À 2dÞ.For example, the movable range of the border node 8 is shown as the orange rectangle in Fig. 4. Especially, its movable range in y direction is ðd; a À dÞ in case that the border nodes get too close to the corner nodes, and its x-coordinate's range is ða À 0:5; a þ 0:5Þ for avoiding the overlap of nodes.In addition, due to the periodicity of the border nodes in the unit cell, the coordinates of Node 5 will alter synchronously with Node 8. Therefore, the design variable n a 1 ;i x and n a 1 ;i y here controls the positions of i-th pair of nodes associated by the lattice vector a 1 f g.Likewise, n a 2 ;j x and n a 2 ;j y denotes the positions of jth pair of nodes associated by the lattice vector a 2 f g.If there are 2 border nodes on the same edge, the relative position of them could not be changed during ) should always sat- Note that if the node on the border has no counterpart to form a node pair by the lattice vectors, it would not be treated as the border node.
Based on the above discussions, the mathematical model of the unimode and bimode materials design with desired mode requirements in this step can respectively be formulated as x ; n a 2 ;1 y ; n a 2 ;2 x ; n a 2 ;2 2 s:t: 2 s:t: Eqs. ( 24) and ( 25) are the optimization models of unimode and bimode materials in Step 2, respectively.The only difference between Eqs. ( 24) and ( 25) lies in the objective function.n f g contains the changing coordinates of the movable nodes in the unit cell.Hence, the first part of the variables in Step 2 are the x and y coordinate of the internal nodes.N 1 here denotes the number of internal nodes.The second part is the x and y coordinates of the border nodes pair which are related by a 1 f g.The final part is the x and y coordinates of the border nodes pair which are related by a 2 f g.Limited by the number of nodes ð4 Â 4Þ in the ground structure, the maximum number of internal nodes is 4. The maximum number of the border node pair is 4, too.Thus, the maximum number of the design variable in Step 2 is 4 The first and second constraints show the upper and lower bounds of the coordinates corresponding to the internal nodes.For border nodes lying at the right boundary of the design domain, the movable space is characterized by the third and fourth constraints.Likewise, the lower and upper bounds of the nodes located on the top boundary of the design domain can be defined.If there are two pairs of border nodes on the same edge, extra constraint n a 2 ;1 x À n a 2 ;2 x 6 Àd; d ¼ 0:1 should be added to the optimization model to prevent the border nodes from getting too close.Such constraints can help to get rid of possible intersection or overlapping.

Numerical homogenization of the extremal materials
In order to calculate the effective elasticity tensor c of a unit cell, the homogenization technique needs to be used.In this paper we take advantage of the Cauchy-Born rule based numerical homogenization method that we proposed in a recent work [32].For completeness we briefly outline the main steps of the homogenization process.The readers can refer to [32] for a detailed discussion on this topic.

Finite element procedures to compute the strain energy density
As in usual finite element analysis, the global stiffness matrix K ½ 2 R n dof Ân dof and the load vector F f g 2 R n dof are available after finite element discretization.The displacement vector u f g 2 R n dof can be found by minimizing the total potential, The potential expressed by Eq. ( 26) has not included appropriate boundary conditions.In our cases, we should apply the Cauchy-Born periodic boundary condition [32] that relate the displacement of opposite nodes on the exterior boundaries of the unit cell.Here it suffices to know that the Cauchy-Born periodic boundary condition can be expressed in a linear form of the displacement vector, where A ½ 2 R mÂn dof is a coefficient matrix containing constant coefficients, b 2 R m is related to the prescribed strain field e exerting on the unit cell.Eq. ( 27) adds m constraints to the system, so the potential expressed by Eq. ( 26) should be slightly modified to include these constraints, where w f g 2 R m is the Lagrangian multiplier (and indeed reveals the reaction forces at the constrained nodes).Explicit expressions of A ½ and b f g will be given in later sections.According to the minimum total strain energy principle, the true displacement of a system would minimizes the total potential.Thus to find the true displacement field, we seek to solve the extremum condition, Eq. ( 29) can be rewritten in a more compact form, Eq. ( 30) (or equivalently Eq. ( 29)) is the state equation of the microstructure.Upon solution of Eq. ( 30), the displacement of the unit cell can be found.Then, the strain energy density w can be expressed as, where V denotes the volume of the unit cell.

The homogenized elasticity tensor
The elasticity tensor of the composite material that are regarded as periodic arrangements of unit cells can be found from the point view of strain energy density w, which can be expressed as, where e ij denotes the prescribed strain field exerting on the unit cell, c ijkl denotes the homogenized elasticity tensor.In matrix form, Eq. ( 32) can be rewritten as, Time again, we emphasize that the fourth-rank elasticity tensor and second-rank strain tensor should be written in the Kelvin notation whenever they appear in the matrix expressions.Comparing Eqs.(31,33), we know that the components of c ijkl can be found by making different choices of e f g.For example, in two-dimensional problems we can set , then by substituting Eq. ( 34) into ( 33) and ( 31) we know that Note that we use the superscript ð1Þ to reveal that this is the first test case to calculate the components of the elasticity tensor.Therefore, once the prescribed strain exerting on the unit cell and the strain energy density w are known, the corresponding component of the elasticity tensor can be obtained.In summary, the elastic moduli of materials in two-dimension can be calculated as follows,

Optimization implementation
The previous two sections give a brief overview of the overall optimization framework without in-depth study on the implementation details.In this section, all these details will be covered, including assembly of the global stiffness matrix, implementation of the Cauchy-Born boundary condition and derivation of the sensitivity information.Some special strategies to improve the optimization effectiveness and efficiency are also elaborated.For example, a logic matrix M ½ is proposed to circumvent the intersec-tion and overlapping of the bars in the unit cell, which is shown to be very useful according to our numerical experience.Moreover, to pursue clear and concise microstructures, the postprocessing process is introduced into the optimization framework to simplify the results from Step 1.The detailed flowchart of the optimization Step 1 and 2 are shown in Fig. 5.
In Step 1, we should first specify the properties of bars in ground structure (like Young's modulus E, maximum crosssectional area of bars S 0 ), size of the unit cell a, scale of nodes mesh (n x Â n y ) and prescribed eigenvector v d f g.The logic matrix M ½ is generated subsequently.The element stiffness matrices K i ½ and the global stiffness matrix K ½ can be formulated and assembled in terms of the topology design variables f f g.By applying Cauchy-Born periodic condition, the displacement vector u f g can be calculated.Thus, the effective elasticity matrix C ½ can be obtained by using the homogenization method mentioned in Section 4. Then the values of objective function f ð f f g; hÞ and constraints gð f f g; hÞ can be calculated.The sensitivity of objective function f ð f f g; hÞ and constraints gð f f g; hÞ can be given by the chain rule if the sensitivity of elasticity matrix is known.Then the stopping criteria can be checked.There are two kinds of stopping criteria in the optimization process.First, the rule of maximum iterations is checked.The optimization process would stop once the number of iterations exceed the maximum allowable iterations.Second, the KKT tolerance is checked.If the KKT condition [33] is satisfied within prescribed precision, the optimization process would stop as well.If any one of the stop criteria is satisfied, the optimization process stops.Otherwise, the design variables will be updated by the optimization algorithm (e.g. the fmincon function in MATLAB).
The process of Step 2 is similar to that of Step 1.However, there are some notable differences between them.First, the input data of Step 2 is the result from Step 1.Second, postprocessing procedure only exists in Step 1. Furthermore, the logic matrix M ½ is introduced in step 1 to deal with the intersection and overlapping between bars.In step 2, intersection is avoided by the classification of nodes and the restriction of the movable range of different types of nodes.So the use of the logical matrix M ½ can be avoided in this step.

Stiffness matrix of the unit cell in Step 1
Now we illustrate how to obtain the global stiffness matrix of the unit cell as shown in Fig. 3.
In Fig. 6 we show a truss element.As mentioned earlier, the truss element can only undergo elongation or shortening when subject to external loads.So there are only two degrees-offreedom, i.e., ½ u 1 ; u 2 > , in the local coordinate system of the truss element O x.In contrast, there are four degrees-of-freedom, i.e., ½u 1 ; v 1 ; u 2 ; v 2 > , in the global coordinate system Oxy.For i-th bar in the ground structure, the stiffness matrix K i ½ in global coordinate system Oxy is, where T i ½ is the transformation matrix, and a is the angle between the bar and the x-axis in the global coordinate system.The stiffness matrix of i-th bar in local coordinate system O x is, where E i ; S i ; L i denote Young's modulus, cross-sectional area, and length of the i-th bar, respectively.We borrow the idea of the SIMP formula, which is famous in solving topology optimization problems of continuum.The main idea is to use fictitious topology description variable to show the presence of an element.Under this idea, the cross-sectional area of i-th element can be given by, where S 0 denotes the maximum cross-sectional area, which equals 1 throughout this paper.S min represents the virtual cross-sectional area of the bar when no solid materials are filled.S min equals 10 À10 , a very small number (but not zero) to prevent numerical singularity.Therefore, according to Eq. (37), the relation between f i and K i ½ can be written as, where ½K 0 and ½K min denote the stiffness matrix with crosssectional area S 0 and S min , respectively.The global stiffness matrix of the whole unit cell can be obtained by assembling the stiffness matrix of all the truss elements, Note that the summation symbol in Eq. ( 42) does not mean the algebraic summation operator, but the well-known assembly procedure in finite element analysis.

Cauchy-Born boundary condition in Step 1
To simplify the description of Cauchy-Born boundary condition, a design domain of square unit cell discretized by 3 Â 3 nodes as shown in Fig. 7 is taken as an example.The lattice vectors of the unit cell are a 1 f g and a 2 f g.Independent nodes on the boundary of the structure are the ones with indices 1, 2 and 4. The other nodes on the boundary are closely related to the independent ones by the lattice vectors under prescribed macroscopic strain field e.
The relation can be given by where u i f g stands for the displacement vector of i-th node.The node labels can be found from Fig. 7. Eq. ( 43) can be easily rewritten in a matrix form like Eq. ( 27).The global displacement vector u f g 2 R 18 consists of the displacements of all nodes in the unit cell.
Here, i-th node has two degrees of freedom ½u i ; v i > .

Sensitivity of elasticity matrix in Step 1
Upon solution of the state Eq. ( 30), the displacement field of the whole unit cell can be available under the prescribed strain field.After that, the strain energy density w can be obtained by Eq. (31).Then components of the elasticity tensor can be obtained by Eq. (36).
From the workflow to calculate c ijkl , we know that the core issue of the sensitivity analysis is to find the derivative of the strain energy density w w.r.t. the design variable f i .In the light of the well-known adjoint sensitivity analysis method [34], the strain energy can be equivalently rewritten as, where, f i is the design variable, l u È É 2 R n dof and l w n o 2 R n dof are two Lagrangian multipliers introduced to facilitate the computation of the sensitivity information considering the state Eq. ( 30).Note that U w (since the state Eq. ( 30) should always be satisfied).
The Lagrangian multipliers are added here intentionally.In principle, l u È É and l w n o can take any values in R n dof , but we will judiciously select values of these Lagrangian multipliers such that the sensitivity information can be expressed in a much more concise manner.This is exactly the main idea of the adjoint sensitivity analysis method [34].
Differentiating Eq. ( 44) gives the sensitivity of U w.r.t. the design variable f i , According to the idea of adjoint sensitivity analysis method, Eq. ( 45) can be simplified if we choose the Lagrange multipliers such that the contents in the parentheses before @ u f g @f i and @ w f g @f i are exactly zero, In a more compact form, Eq. ( 46) can be rewritten as, Therefore, the Lagrange multiplier l u È É can be obtained by solving the linear equation system Eq.( 47), so the expression of the sensitivity Eq. ( 45) can be expressed as, Due to the assembly nature (cf.Eq. ( 42)) of the finite element analysis procedure, @ K ½ @f i can be given in terms of the elemental stiffness matrix.Therefore, differentiating Eq. (41) gives the sensitivity, Thus Eq. ( 48) can also be given in terms of the elemental stiffness matrix as,

Dealing with intersection and overlapping
There are many intersection and overlapping bars in the ground structure, but in the optimized results we do not want the intersection and overlapping occurs.Here we propose a simple but effective strategy to control the occurrence of intersection and overlapping bars.A logical matrix M ½ 2 0; 1 f g N bar ÂN bar is defined to show whether there is intersection or overlapping between any two bars in the ground structure.The size of matrix M ½ is N bar Â N bar , and all of its entries have values that are either 0 or 1. M ij ¼ 1 means that i-th and jth bar intersect or overlap in the ground structure, and M ij ¼ 0 reveals they do not.Obviously, M ½ is a symmetric matrix, and its diagonal elements are all 0 (since any bar will not intersect or overlap with itself).
Suppose there are two line segments ab and cd in the plane with coordinates of endpoints ðx a ; y a Þ; ðx b ; y b Þ; ðx c ; y c Þ; ðx d ; y d Þ.Using two parameters /; u 2 ½0; 1, line segments ab and cd can be represented as, Equating Eqs. ( 51) and ( 52), the intersection points of ab and cd can be found by solving The determinant of the coefficient matrix is If D -0, the line segments ab and cd would cross over at a certain point, but not necessarily at the interior of line segments ab and cd.In detail, 4 cases occur according to the values of / and u when As shown in Fig. 8(a), if 0 < / < 1; 0 < u < 1 (i.e., Case 1 as shown in Eq. ( 55)), then the intersection point lies at the interior of both bars.As shown in Fig. 8(b) and (c) (i.e., Case 2 and 3 as shown in Eq. ( 55)), if either of the /; u equals to 0 or 1, while the other one is bounded between ð0; 1Þ, the intersection point lies at the endpoint of one bar and the interior of the other bar.When D ¼ 0, two possibilities exist, i.e., the two bars can be either parallel or collinear, the former is definitely allowed while the latter is possibly not.So we need to distinguish between these two subcases.As shown in Fig. 9 denotes the area of the parallelogram.If S } > 0, it means that line segment ab and cd are parallel.S } ¼ 0 reveals that ab ! and cd ! are collinear.For two collinear bars, they can always be sorted according to the coordinates in one direction (x or y).As shown in Fig. 10, in the x direction, i.e. from left to right, L and R represent the left and right nodes of the bar, respectively.The subscript number represents the index of the bar.There are four possible situations between the two bars as shown in Fig. 10.Then it follows that the vector dot product can be used to distinguish whether the two bars are overlapped or not, which can be written in the following form: As shown in Fig. 10(a), r > 0 reveals that there is no overlapping between the bars.As shown in Fig. 10(b), r ¼ 0 reveals that the only overlapping between two bars is the node R 1 and L 2 , and can be regard as non-overlapping.As shown in Fig. 10(c) and (d), r < 0 reveals that two bars do overlap and the overlapping region is painted green in Fig. 10.Following the abovementioned ideas, the workflow of the generation of the logical matrix M 2 R NÂN can be shown as Fig. 11 and concluded as the follows.53).When any one of the 3 cases in Eq. ( 55) is satisfied, intersection occurs, so set M ij ¼ 1 and go to Step 6. ( 6) Stop.

Postprocessing the results from Step 1
By using the optimization models in Step 1, many feasible solutions can be found.But the solutions still contain some flaws inevitably.For example, since the penalty terms of grey elements are only listed in the objective function, a certain part of results would still have grey bars.There may also be some useless bars that have no stresses when subject to overall strain and thus have no contribution to the effective elasticity matrix.
Based on the discussion mentioned above, the postprocessing is proposed to help us find the meaningful configuration.Firstly, a hyper-parameter b is proposed to determine whether the bar should be retained or not, i ¼ 1; 2; ::; N bar .In this paper, we use b ¼ 0:3.After that, the entries in the projected solution f p f g are either 0 or 1, and the bars whose f p i equal to 1 will be retained while the rest will be discarded.N p denotes the number of bars in the projected solution.
Then the useless bars can now be identified and deleted.The elasticity matrix and its eigenvalues can be obtained from the pro-  jected solution.For extremal materials, if the number of zeroenergy modes is n, it contains n soft modes and 3 À n hard modes in two-dimension.A useless bar is defined here the bar that has no contribution to any hard mode.After the elasticity matrix and its eigenvalues are obtained, the eigenvalues can be sorted in descending order, and the eigenvectors corresponding to the first 3 À n eigenvalues are taken as the applied strains.Assuming that one of the applied strain is e Ã , the displacement of all nodes in the structure can be calculated by Eq. ( 27).The strain e # i ði ¼ 0; 1; . . .; N p Þ of each bar can also be available from the FEA solver.According to the strain e # i of each bar, calculate the average strain e # i .Now we can judge whether this bar is redundant or not by , the i-th bar should be retained and it has contribution to the elasticity matrix.f f i 0, the i-th bar should be discarded since it has no contribution to any hard modes.

Sensitivity of the elasticity matrix in Step 2
Different from the first step, the second optimization step will change the positions of nodes, inevitably leading to variations of the length and direction of the elements in the unit cell.The length of i-th element can be expressed as, where L i is the length of i-th bar in the unit cell; ðx il ; y il Þ and ðx ir ; y ir Þ are the coordinates of left and right nodes of the bar.Assume that x il is the design variable in Step 2, the sensitivity of L i can be calculated as, Similar expressions can be derived for @L i @x ir ; @L i @y il , and @L i @y ir .The change of the node positions will also change the direction of the elements, cos a ¼ Assume that x il is the design variable, the sensitivity of i-th bar's direction can be calculated as, Therefore, the sensitivity of the transformation matrix can be written as, @ T i ½ @x il ¼ @ðcos aÞ @x il @ðsin aÞ @x il 0 0 0 0 @ðcos aÞ @x il @ðsin aÞ It is worth noting that the design variables of the border nodes actually control the coordinates of two closely-related nodes, so the length and direction changes caused by the movement of the node pair needs to be considered.
The sensitivity of the stiffness matrix of i-th bar in the local coordinate system K i Â Ã can be written as, The sensitivity of stiffness matrix of i-th bar in global coordinate system K i ½ can now be written as, 6. Rigid rotation of the unit cell in 2D and its properties Let the coordinate system xoy rotate angle h clockwise (or equivalently let the unit cell rotate angle h counterclockwise) around the z axis to form a new coordinate system xoŷ, the stress and strain expressed in Eq. ( 4) can be written in the new coordinate system as, r11 r12 r12 r22 e 12 e 22 !cos h sin h Using the Kelvin notation, the transformation law of the stress and strain tensor can be rewritten as where the transformation matrix is, Time again, we emphasize that the vector form of stress and strain tensor are written in Kelvin form (see Eq. ( 4)), not the commonly used Voigt form.The constitutive equation in the rotated coordinate system xoŷ is Substitute Eq. (69) into Eq.(71)we have , the transformation law of the elasticity matrix can be easily derived as Let ðk; e f gÞ denote the eigenpair of C ½ in coordinate system xoy, cf.Eqs.(1,2).Then an important property of the transformation law in Eq. ( 73) is that the elasticity matrix keeps its eigenvalue k i unchanged (since the transformation shown in Eq. ( 73) is an orthotropic transformation), while the rotated eigenvector 3 is It is required that the norm of the strain e f g in Kelvin form is 1 (which is a common practice when calculating the eigenvalue problem), In 2D cases, the first and second principal invariants of the strain tensor are, Moreover, to Eq. ( 75), J 2 is indeed not independent and can be expressed using J 1 , From the property of the first principal invariant, we have Let n f g ¼ ½1; 1; 0 > , it can be easily verified that, ð e Eq. (80) reveals that the vector n is the normal vector of the plane containing e f g and its rotation counterparts ê f g with arbitrary rotation angle h.Obviously, the intersection curve of the plane (with normal vector n) and the unit sphere is a circle.The center and radius of the circle can be given by O e 11 þ e 22 2 ; e 11 þ e 22 2 ; 0 ; as shown in Fig. 12. Equivalently, the center and radius of the circle can also be expressed using the principal invariants, The central angle formed by e f g and ê f g is \color{black}2 \theta ; where, Oe f g ¼ From Eq. ( 82) we know that the center and radius of the circle formed by the rotations of ê can be fully determined by the first invariant of the original strain tensor.Thus, to test the applicability of our method in finding extremal materials for arbitrary soft modes (or hard modes), we can simply select different strains with different J 1 but all with unitary norm.In other words, the unit sphere as shown in Fig. 12 can be sliced into several pieces with the slicing planes normal to n f g ¼ ½1; 1; 0 > .Each slicing plane intersects the sphere to form a circle and each circle has a unique value of J 1 .Then it is convenient to test if the proposed method can find the extremal materials corresponding to different values of J 1 .
Now we determine the range of J 1 and J 2 .This can simply be done by observing again Fig. 12.The range of the radius is From Eq. (85) the range of J 1 can be obtained, Then the range of J 2 can be determined as well, Based on the above discussions, in Kelvin form the rotation transformation does not change the eigenvalue of the elastic matrix, while the new eigenvector is rotated by the transformation matrix.This means that the number of zero eigenvalues remains the same while the soft or hard modes rotate according to the transformation matrix.

Numerical examples and discussions
In this section, we will test the performance of the proposed design framework in designing extremal materials.Several cases will be tested.The first one is to find the microstructure of the planar unimode and bimode materials without further constraints on soft or hard modes.Afterwards the desired soft mode is taken into consideration when designing the unimode materials and the desired hard mode is taken into consideration when designing the bimode materials.Subsequently, the other two modes are also considered.For example, when designing the unimode materials the soft mode is properly considered by using the proposed topology optimization step and shape optimization step, then the two other hard modes can be tuned by modifying the cross-sectional areas of the bars.

Design of planar extremal materials without constraints on soft or hard modes
In this example we attempt to find various kinds of planar unimode and bimode materials without any constraint on the soft or hard modes.In case of the optimized unit cell being too complex and not easy to manufacture, we add constraint on the number of bars.This constraint can be written as P N bar i¼1 f i 6 N pre .We set maximal number of bars to be N pre ¼ 6, thus the mathematical model for optimization of unimode and bimode materials can respectively be given by find f The size of design domain a ¼ 3, and the unit cell is discretized with 4 Â 4 nodes.The number of bars in the ground structure is N bar ¼ 108.Young's modulus E ¼ 1, maximum cross-sectional area of bars S 0 ¼ 1.The interior-point algorithm is chosen as the optimization algorithm.For exploring more microstructures of twodimensional extremal materials, the initial solution in the algorithm is randomly generated within the lower and upper bounds.The maximum iteration equals 2000, and the maximum allowable violation of constraints is set to 10 À6 (except the tolerance of second constraint in unimode materials design is set to 0.001).Termination tolerance for first-order optimality is set to 10 À6 .Numerous planar unimode and bimode materials can be found by the optimization model and here only a small portion is showed.The unit cell and super cell of the unimode and bimode materials results are shown in Fig. 13 and 14, respectively.The corresponding effective elasticity matrices and eigenvalues are listed in Table 1 and 2. The supper cell here is comprised of 3 Â 3 unit cells.
The Results 1 and 2 of the unimode materials shown in Fig. 13 can actually be categorized into the same group, because they both can be regarded as the bimode materials (i.e., the planar honeycomb) plus one more bar stretching over the unit cell, making the soft mode number to be only 1. Looking into Table 3, one can find that the second and third column of the ½C matrix of Result 1 are proportional, and both columns are not proportional to the first column (with the proportion to be À0:3535), making the rank of the ½C matrix to be 2, which is consistent with the unimode design target.Similarly, the first and second column of the ½C matrix of Result 2 are exactly the same, and both columns are not proportional to the third column.In view of super cell, the features can be seen more clearly, Result 1 can be regarded as a combination of negative honeycomb planar pentamode and equidistant horizontal bars.Similarly, Result 2 can be decomposed Result 3 is more sophisticated.The super cell is made up of three kinds of irregular quadrilaterals, without any opposite edges parallel.Its effective elasticity matrix in the Table 1 also witnesses the complexity, which only possess the basic symmetry while all the six entries in the upper triangular part are totally different.Moreover, such microstructure is still a unimode material even when the center node of the lower X-shape moves within the rectangle formed by the corner nodes of the X-shape.The soft mode are related to the coordinatesðx; yÞ of the center node in X-shape.By using method mentioned in [35] which takes a thorough matrix analysis of the truss model to find the intrinsic connection between the soft mode and topological connectivity, the soft mode (without normalization) of such microstructure can be given by, It should be noted that the soft mode can be calculated by Eq. ( 90) only when 0 < x < 3 and 0 < y < 2.
From the supper cell shown in Fig. 14, it is easy to find that the Result 1 of the optimized bimode materials indeed corresponds to the commonly seen honeycomb bimode material but inclined 45 to the right side.Results 2 and 3 have never been seen in the literature to the best of our knowledge, and they both contain a triangle in the unit cell.Conceptually each triangle in the unit cell can be replaced by one rigid body.Besides, Results 2 and 3 both include the irregular heptagon feature, which has not been noticed by our community ever before.According to the effective matrices shown in Table 2, Result 1 and 2 share many similarities, i.e., the first two columns are the same and are proportional to the third column, making the rank of the ½C matrix to be 1, which is consistent with the bimode design target.Since the entries in the second row and column of the effective elasticity matrix of the Result 3 are all zero, one soft mode of the Result 3 can be directly given by v s f g ¼ ½0; 1; 0 > .
From the examples shown in this subsection, it can be found that by using the optimization framework proposed in this paper it is easy to discover numerous innovative designs of extremal materials that have never been found before, which is exactly the purpose of this paper.

Planar extremal materials with constraints on soft or hard modes
Now we proceed to design planar extremal materials with designated soft or hard modes.Planar unimode materials have only one soft mode while planar bimode materials have only one hard mode.Therefore, it is a good practice to design the unimode material with designated soft mode while the bimode material with designated hard mode.The mathematical models in step 1 and 2 for unimode and bimode materials design have been list in Section 4. The overall optimization process is shown as Fig. 2. The performance control indices are g 1 ¼ 0:01; g 2 ¼ 10 À5 in this subsection.
As discussed in Section 6, rotation of the unit cell about z axis does not change the eigenvalue of the elasticity tensor, but does rotate the eigenvector.What's more, it has been shown that the trajectory formed by all rotated eigenvectors is actually a standard planar circle.The normal vector of the plane that the circle is located on is shown to be n f g ¼ ½1; 1; 0 > .The center and radius of the circle have shown to be dependent on the first principal invariant J 1 .As shown in Fig. 15, the smaller the value of J 1 , the larger the radius.Thus in this example we verify the effectiveness of the proposed optimization framework by designing unimode materials with designated soft modes that have different values of J 1 .In detail, the unit sphere is sliced using five parallel planes, then five circles are formed with 5 values of J 1 .The values of the strain vector are then determined by choosing the one with largest value of ffiffiffi 2 p e 12 .The five desired soft mode are listed in Table 3.
When designing unimode materials, the unique soft mode is to be constrained to be the ones shown in Table 3.Five super cells comprising 3 Â 3 unit cells with different values of J 1 are shown in Fig. 16.The unit cell is outlined by the pink dotted box.For Case 1 and Case 5, merely the optimization step 1 is enough to obtain the unimode materials satisfying all the conditions.The details of optimized configurations, such as the soft mode, eigenvalues, the index to judge the accuracy of the soft mode and the rotation angle are listed in Table 4. Since the results shown in Fig. 16 have not illustrated the rotation angle, one example of the rotated super lattice and unit cell are given in Fig. 17 for more intuitive understanding of the final optimization results.
When designing bimode materials, the unique hard mode is to be constrained to be the ones also shown in Table 3.Five super cells comprising 3 Â 3 unit cells with different values of J 1 are shown in Fig. 18.Case 5 uses only the first optimization step to get the optimized microstructures satisfying the conditions.The details of optimized configurations, such as hard mode, the index to judge the accuracy of the hard mode and the rotation angle are listed in Table 5.

Planar unimode materials with constrains on both soft and hard modes
In this example, we design the unimode materials with constraints on both the soft and hard modes.Since spaces of soft modes and hard modes are orthogonal, they can be designed independently.The soft modes are dealt with by the optimization framework illustrated in Section 3 while the hard modes are appropriately constructed by tuning the cross-sectional areas of the bars.
Hutchinson [35] proved that the soft modes are only relevant to the location of joints and have nothing to do with the crosssectional areas of the bars (as long as the cross-sectional areas are finite).Thus when tuning the hard mode of the unimode material that has already satisfied the constraint on the soft mode, it is wise to take the cross-sectional areas of the bars to be the design variables while keeping the node positions fixed.Let C ½ denote the elasticity matrix of the unit cell (that is obtained from the optimization Step 1 and 2 and satisfies the constraint on soft mode), v h f g denote the desired hard mode.There are two hard modes for a unimode material and v h f g here stands for anyone of them, the other one can simply be obtained by v s f gÂ v h f g.If the optimal solution is found, C ½ v h f g will be a vector parallel to v h f g.Hence the cross product of C ½ v h f g and V h f g can be taken as the objective function to be minimized.The only constraint is that the crosssectional areas of all the bars should be larger than 0. Therefore, the optimization model for tuning the hard mode of the unimode materials can be given by, find 1 f g ¼ 1 1 ; 1 2 ; 1 3 ; 1 4 ; 1 5 ; 1 6    where the design variable 1 f g represent the cross-sectional area of the six bars in the unit cell, and the label of the bars is shown in the Fig. 19.
The third case of unimode material shown in Fig. 13 is taken as an example here to show the design process of the planar unimode materials with constraints on both the soft and hard modes.Now it is sketched in a new style in Fig. 19, where the numbers with and without parentheses denote labels of node and bar element, respectively.
The original soft mode of the unimode material can be calculated by Eq. (90), after normalization it equals ½0:5523; À0:7365; 0:3906 > , while the hard modes are ½0:2010; 0:5723; 0:7950 > and ½À0:8090; À0:3606; 0:4641 > .We set three cases with three different desired hard modes which are list in the Table 6.The optimization results including the optimized hard mode, objective function values and optimized design variables are list in Table 6.The three desired hard modes are preset 0 in the third, second and first dimension, respectively.By size optimization of the cross-sectional areas of the six bars, the optimized hard modes are very close to the desired ones.The supper cells and unit cells of optimized unimode materials with 3 kinds of prescribed hard modes are given in Fig. 20.
For a bimode material, it has two soft modes and one hard mode.As in the previous section, we have completed the design

Typical computational time with different number of design variables
In previous examples we do not list the computational time, so in this subsection we give the information of the computation time with different number of design variables so that the readers can have a basic understanding of the efficiency of the proposed optimization framework.
As shown in Fig. 21, we consider three cases with different number of bars in the ground structure.Detailed statistics on the number of nodes, number of bars (which is the same as the number of design variables), and the corresponding computation time are all list in Table 7. From Table 7 it is clear that the proposed optimization framework is highly efficient.

Conclusion
In this paper, a versatile optimization framework to design the extremal materials is proposed.The necessary and sufficient conditions for different kinds of extremal materials are carried out.By using the two-step optimization scheme (topology optimization + shape optimization) a series of extremal materials can be found.No prior knowledge of the symmetry of the unit cell is made use of in the proposed design workflow.Therefore, the optimized results obtained in this study can in principle be any general anisotropic material.In topology optimization stage, the logic matrix which needs to be generated only once before optimization is proposed to avoid the intersection and overlapping between the bars in the ground structure.The postprocessing process can remedy the imperfections (e.g., grey elements and redundant bars) in the results obtained in the first step so that the final solutions are clear and concise.The shape optimization step is proposed to tune the positions of the nodes in order to better match with predefined soft or hard modes.By analyzing the rotation process, it is found that in two-dimensional case the strain tensor in Kelvin form can be classified according to the values of the first principal invariant J 1 .That is to say, strain tensors with the same J 1 can be regarded as the rigid rotation results of each other, so only one of them needs to be considered during the optimization.For two-dimensional unimode materials, both the soft mode and the hard modes can be tuned successively by the optimization framework, where the soft mode is designed by the topology optimization and shape optimization while the hard mode is adjusted by the size optimization.

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

Fig. 2 .
Fig. 2. (a) Design flowchart of 2D extremal materials considering predefined soft or hard mode.(b) A diagram showing the evolution of the unit cell during the two optimization steps.

Fig. 1 .
Fig. 1.Periodic unit cells.The left side shows the periodic media composed of tessellations of unit cells, the right shows one basis of the unit cell.Red dots denote the corner nodes in the lattice, black lines denote the bars.

Fig. 3 .
Fig. 3. Two-dimensional ground structure model of the unit cell.

Fig. 4 .
Fig. 4. Classification of the nodes in the unit cell.

Fig. 6 .
Fig.6.A truss element and its displacements in local and global coordinate system.

Fig. 5 .
Fig. 5.The flowchart of optimization process on design of extremal materials with prescribed soft or hard mode.

( 1 )
Get the coordinates of the four nodes, ðx a ; y a Þ; ðx b ; y b Þ; ðx c ; y c Þ; ðx d ; y d Þ. Set initial values M ij ¼ 0. (2) For every two bars calculate the determinant D according to Eq. (54).If D ¼ 0, set M ij ¼ 1, go to Step 3. Otherwise, go to Step 5. (3) Check if the two bars are parallel through Eq. (56).If S } > 0, then set M ij ¼ 0 and jump to Step 6.Otherwise, go to Step 4. (4) Since D ¼ 0; S } ¼ 0, then the two bars must be collinear.Check if the two bars overlap through Eq. (57).If r !0, then set M ij ¼ 0, otherwise set M ij ¼ 1. Go to Step 6. (5) When D -0, that is, the two bars are neither collinear nor parallel.The values of the parameters /; u can be obtained by Eq. (

Fig. 9 .
Fig. 9.The parallel relation between two bars judged by the area of parallelogram S} > 0.

e
v f g fall on an unit sphere whose three axes represent e 11 ; e 22 ; ffiffiffi 2 p e 12 , respectively.According to Eq. (74), the new strain after rotation can be expressed as, 11 cos 2 h þ e 22 sin 2 h À 2sinhcoshe 12 e 11 sin 2 h þ e 22 cos 2 h þ 2sinhcoshe 12

Fig. 12 .
Fig. 12. Rotations of e f g are located at a circle.

Fig. 19 .
Fig. 19.Unit cell of a unimode material.The numbers with and without brackets denote indices of the nodes and bars, respectively.

Fig. 18 .
Fig. 18.Optimization of bimode materials with constraints on hard mode.
(a), two thresholds g 1 and g 2 are introduced to help algorithm judge the workflow by comparing to the value of performance index Err.If Err > g 1 , it means the precision of solution is too coarse, thus further loops in Step 1 are needed.If g 1 < Err < g 2 , it means that the precision is acceptable for Step 1 but not for Step 2, so further loops in Step 2 are needed.If Err < g 2 , , this can be easily done by calcu-

Table 1
Effective properties of the optimized unimode lattices.

Table 2
Effective properties of the optimized bimode lattices.

Table 3
Five Selected eigenvectors.

Table 4
Effective properties of optimized unimode lattices with desired soft modes.

Table 5
Effective properties of bimode materials with desired hard modes.

Table 6
Effective properties of unimode lattices with desired hard modes.

Table 7
Computation time with different number of bars.