Topological Design of Microstructures of Materials Containing Multiple Phases of Distinct Poisson’s Ratios

Kai Long,Xiaoyu Yang,Nouman Saeed,Zhuo Chen and Yi Min Xie

1State Key Laboratory for Alternate Electrical Power System with Renewable Energy Sources,North China Electric Power University,Beijing,102206,China

2Centre for Innovative Structures and Materials,School of Engineering,RMIT University,Melbourne,3001,Australia

ABSTRACT A methodology for achieving the maximum bulk or shear modulus in an elastic composite composed of two isotropic phases with distinct Poisson’s ratios is proposed.A topology optimization algorithm is developed which is capable of finding microstructures with extreme properties very close to theoretical upper bounds.The effective mechanical properties of the designed composite are determined by a numerical homogenization technique.The sensitivities with respect to design variables are derived by simultaneously interpolating Young’s modulus and Poisson’s ratio using different parameters.The so-called solid isotropic material with penalization method is developed to establish the optimization formulation.Maximum bulk or shear modulus is considered as the objective function,and the volume fraction of constituent phases is taken as constraints.The method of moving asymptotes is applied to update the design variables.Several 3D numerical examples are presented to demonstrate the effectiveness of the proposed structural optimization method.The effects of key parameters such as Poisson’s ratios and volume fractions of constituent phase on the final designs are investigated.A series of novel microstructures are obtained from the proposed approach.It is found that the optimized bulk and shear moduli of all the studied composites are very close to the Hashin-Shtrikman-Walpole bounds.

KEYWORDS Poisson’s ratio; topology optimization; homogenization; bulk modulus; shear modulus

1 Introduction

Methods to achieve optimal engineering designs have been intensively investigated for decades.Various topology optimization approaches have been proposed and developed in the field of structural designs [1-5], such as homogenization method [5], solid isotropic material with penalization(SIMP) [6-8], evolutionary structural optimization (ESO) [9,10] and bi-directional evolutionary structural optimization (BESO) [11,12], and level set method [13,14].The application of topology optimization has also been extended to the design of cellular and composite materials.Sigmund [15] initially incorporated an inverse homogenization concept to the topology optimization of materials with exceptional mechanical properties.Inspired by this approach, many attempts have been made to develop materials with different functionalities, including elastic [16], thermos elastic [17,18], piezoelectric [19,20] and auxetic materials [15,21,22].

The bulk modulus and shear modulus are generally devoted to describe the stiffness of an isotropic elastic material.And earlier work in topology optimization of materials normally involves maximizing the bulk modulus or shear modulus of a 2D or 3D material.Except for the existing composites such as composite sphere, finite rank laminates and Vigdergauz microstructures, Sigmund [23] developed a series of two-phase microstructures and their bulk moduli coincided with the well-established Hashin-Shtrikman (HS) bounds.Gibiansky et al.[24] extended this work to three-phase material design, generating a range of novel microstructures which exhibited achievability of the Hashin-Shtrikman-Walpole (HSW) upper bounds.Guest et al.[25]presented an optimal design for maximum elastic stiffness and fluid permeability.Challis et al.[26]employed the level-set method to maximize stiffness and conductivity of two-phase materials.Microstructures contained the so-called Schwartz interface and diamond minimal surface have been observed, which had maximum bulk modulus and conductivity.Andreassen et al.[27] proposed a method for designing manufacturable materials with negative Poisson’s ratio and revisited the problem of maximizing the bulk modulus.Huang et al.[28] utilized BESO approach to design 2D and 3D microstructures with extremal bulk and shear moduli, which demonstrated certain advantages over SIMP approach in 0-1 representation and convergence speed.The BESO methods have been extended to similar issues to obtain isotropic composite microstructures [29], multiphase materials with maximum stiffness and conductivity [30], functionally graded materials [31],and phonoic band gap crystals [32].Unlike regular elastic materials, viscoelastic materials exhibit both viscous and elastic characteristics when deforming, which are widely adopted in reducing vibrations.Andreasen et al.[33] further studied the homogenization concept of viscoelastic composites.Several new microstructures were found to have viscoelastic bulk modulus close to the theoretical upper bound.In recent years, much works have focused on the topological design of metamaterial [34-36].

In comparison to traditional cellular materials composed of a solid phase and a void phase,composites made of two or more materials are more preferable to practical applications due to their synthetic mechanical properties.Conventional topological design of composite microstructure with maximized stiffness always involves constituent phases with different Young’s moduli and the same Poisson’s ratios.However, Liu et al.[37,38] revealed that the effective Young’s modulus of composites could exceed the Voigt estimation if one constituent phase approached the incompressibility limit.Kocer et al.[39] further investigated this phenomenon in materials containing inclusions with negative Poisson’s ratio (NPR).They found significantly enhanced stiffness when alternating different layers of auxetic or non-auxetic materials.Other research also observed evident increase in stiffness beyond any individual constituent phase with positive or negative Poisson’s ratio [40,41].The stiffening effect could be governed by multiple factors simultaneously, such as the values of Poisson’s ratio of constituent phases, topologies of the NPR inclusions, and so on.More recently, Pasternak et al.[42] evaluated the coefficient of thermal expansion and thermal stresses by adding NPR inclusions into positive Poisson’s ratio (PPR)material.They also found reduction in the thermal stress due to Poisson effect.Zuo et al.[43]proposed an optimization method to maximize the stiffness in both longitudinal and transverse directions of laminate composites based on the explicit effective stiffness function.An in-depth investigation was conducted to find the influence of effective Young’s moduli when Poisson effect was considered.Long et al.[44] established a topology optimization framework for maximizing the effective Young’s modulus of composite containing multiple constituent phases with distinct Poisson’s ratios.However, no reports can be found in literature regarding the design of composites with extremal bulk or shear modulus made of constituent phases with significantly different Poisson’s ratios.

The rest of the paper is organized as follows.The homogenized effective elasticity matrix is calculated within periodic base cell (PBC) in Section 2.Also, the sensitivity of the elastic matrix is derived from the interpolation of both Young’s modulus and Poisson’s ratio.Section 3 establishes the topological optimization formulation the microstructure.Section 4 provides several typical numerical examples to validate the effectiveness of the proposed approach.Finally, the main conclusion of this paper is drawn in Section 5.

2 Elasticity by Homogenization and Sensitivity Analysis

2.1 Homogenization

For a two-phase periodic heterogeneous material, the macroscopic elasticity tensorof a periodic composite material can be obtained by where |Ω| denotes the volume of PBC.The locally varying stiffness tensor,Epqrs, depends on the distribution of material in PBC.are the prescribed macroscopic strain fields whileare defined asTheV-periodic displacement fieldsχklcan be found by solving the elasticity equations with a prescribed macroscopic strain

HereΩis theΩ-periodic admissible displacement field corresponding to six strain fields, which are derived from six linearly independent unit tests.These unit strains areε1=(1,0,0,0,0,0)T,ε2=(0,1,0,0,0,0)T,ε3=(0,0,1,0,0,0)T,ε4=(0,0,0,1,0,0)T,ε5=(0,0,0,0,1,0)T,ε6=(0,0,0,0,0,1)T.

Based on Eq.(2), the effective elasticity matrix can be obtained by using the numerical homogenization approach as described in the work [45,46].

2.2 Sensitivity Analysis of Elasticity

The constituent phases are assumed to be isotropic in this study.To derive the sensitivities, a SIMP function for both the Young’s modulus and Possion’s ratio is adopted to interpolate between two phases in each element [6].The SIMP method is developed in the current study and the mathematical formulation reads

whereEjandνjare the Young’s modulus and Poisson’s ratio for phasej(j=1,2), respectively,andxedenotes the relative density of the eth element, which can be purely phase 1 whenxe=1 or phase 2 whenxe=0.αandβrepresent the penalization exponents.Typically, whenα=4 andβvaries from 1 to 3, reasonable results for all the examples in this paper can be obtained.From Eq.(3), we have

We assume the materials used to build the composite are isotropic, the elastic matrix for the elementecan be expressed as

whereλeandµeare Lamé’s first and second parameter for the material in elemente, respectively.Lamé’s parameters can be computed from Young’s modulusEandνusing the relations:

Here for simplicity, the subscript e is omitted in Eq.(6).

From Eq.(6), we have

Using the chain rule, the sensitivity of the elastic matrix for the elementewith respect toxecan be calculated as

Thus, the macroscopic elasticity tensor with respect to density variable can be derived by using the adjoint method, which is

3 Inverse Homogenization Formulation and Solution

The bulk modulusKand shear modulusGcan be treated as one of the stiffness measures,which can be optimized by using the method provided in literature [28].When the microstructure of composite is cubic or isotropic symmetry, the bulk modulus and shear modulus can be expressed in regard to the 3D effective elasticity tensor, i.e.,

To obtain optimal design of composites with homogenized properties, the optimization problem is defined through maximizing bulk modulus or shear modulus subjected to volume fraction constraints of phase 1.

whereV1(x) andV0are the volume of phase 1 material and the volume of design domain,fis the target volume fraction.

4 Analytical Bounds for Materials

The bounds of physical properties are of significant importance in composite design.HS bound is one of the most popular estimations, which defines tight theoretical limits of bulk and shear moduli for isotropic materials, in terms of volume fraction of two-phase composites based on variational principles.Such bounds are achievable when the microstructures and volume fraction of individual phase are determined [47].If the phase exhibits larger (smaller) bulk modulus(K) and larger (smaller) shear modulus (G) simultaneously, the materials can be recognized to be‘well-ordered’, i.e.,(K2-K1)(G2-G1)＞0.Walpole [48,49] extended the HS bounds to the cases in which the bulk modulus and shear modulus are not well ordered (i.e.,(K2-K1)(G2-G1)＜0).

The bulk and shear moduli of each phase are defined as

For two-phase composites with 3D microstructure, the HSW upper limits on effective bulk modulusKUand shear modulusGUcan be found by

whereKmaxandGmaxare the maximal values of bulk and shear moduli for constituent phases,respectively.is the given volume fraction of phase 1.

The HSW upper bound is attainable if the following conditions are satisfied.

Here in it is assumed thatG2≤G1.The HSW upper bounds of shear moduli are attainable under the following condition:

5 Illustrative Examples and Discussion

In this section, several numerical examples are presented to demonstrate the effectiveness and capabilities of the proposed optimization method.It is assumed that in all cases the PBC,which represents the microstructure of a composite, is discretized to 40×40×40 brick elements with 8nodes.To eliminate the checkerboard patterns and mesh-dependency phenomena which is popular in topology optimization, filtering technique is applied [3].Based on the first-order derivatives, the Eq.(12) can be usually solved by a magnitude of optimizers, such as optimality criteria, convex linearization or the method of moving asymptotes and so on [3,49].

In this study, we simply define two different initial designs.In initial design 1, the element density increases proportionally to the distance from the element center to the base cell center, as shown in Fig.1a.Inversely, in initial design 2, the element density decreases proportionally to the same distance as illustrated in Fig.1b.The color bar shows the range of density (xe) with the value varying from 0 (blue) to 1 (red).

For easier identification, the phase 1 in the microstructure is displayed in red, while phase 2 in blue.

5.1 Example I

The objective of topology optimization of this example is to obtain microstructure with maximum bulk modulus.In this example, the material properties are defined as:Young’s modulus,E1=2.5 for phase 1 andE2=1.0 for phase 2; Poisson’s ratio,ν1=0.0 for phase 1 andν2=-0.5 for phase 2.The bulk moduli areK1=0.8333 andK2=0.1190 for the two phases, respectively.The prescribed volume fraction of phase 1 is 60%.The obtained optimal microstructures and their effective matrices are given in Fig.2.

Figure 1:Density distribution for different initial designs:(a) Initial design 1; (b) initial design 2

It can be seen that the obtained microstructures are cubic symmetry, thus the two-phase composites are orthotropic.The bulk moduli can be calculated from the obtained microstructures as 0.5148 and 0.5128, starting from initial designs 1 and 2, respectively.These results are slightly below the HSW upper boundKU=0.5159 according to Eq.(15).As revealed in literature [17],completely different topologies of microstructures could be achieved, depending on initial guesses.Nevertheless, they may have similar values of maximum bulk modulus due to the non-uniqueness of the solution.The interface between the two phases in Fig.2b can be seen as an approximation to a Schwartz P minimal surface which was reported to be structures with maximum permeability or conductivity [25,26].Interestingly, the shear moduli for the two optimized composites are 1.1410 and 1.1407, respectively, also slightly below HSW upper bound of shear modulus in this caseGU=1.1432 according to Eq.(15).

The stiffness of resulting composites in one or more directions could be enhanced, as the benefits of using phase with significantly different Poisson’s ratios.For the purpose of comparison,the upper bound of effective Young’s modulus estimated by Voigt theory can be calculated asEVoigt=fE1+(1-f)E2.The resulting effective Young’s modulus from elasticity matrix of composite material in three directions is 1.9716 and 1.9679 in Fig.2, which exceeds the Voigt limit 1.9000.These results clearly demonstrate the additional advantage of composites containing multiple phases of distinct Poisson’s ratios.

To better understand the influence of constituent phase on the topology of optimal composites, optimization process on design domain of merely phase 1 is performed.The two optimal microstructures and the corresponding elasticity matrices of the composite are presented in Fig.3.The void elements starting from initial design 1 are displayed in blue color.To provide a clear angle of view, one half of PBC is shown in Fig.3a.

Figure 2:Base cell and effective elasticity matrix of composite material, starting from different initial designs:(a) From initial design 1; (b) From initial design 2

The microstructure obtained from initial design 1 in Fig.3a is a hollowed cube with rounding interior corners.The microstructure obtained from initial design 2 in Fig.3b is actually composed by walls in three directions.The unit cell is partitioned into eight individual spaces.These two microstructures are similar to the designs with maximum bulk modulus in literature [23,25].Comparing Figs.2 and 3 in the process of finding microstructures with achievable HSW upper bound, significant difference in the distribution of phase 1 between two-phase composites and cellular materials can be witnessed.

Figure 3:Base cell and effective elasticity matrix of cellar material, starting from different initial designs:(a) Cross-section from initial design 1; (b) From initial design 2

5.2 Example II

To investigate the effect of the Poisson’s ratio on the maximum bulk modulus and the HSW bounds, the change in Poisson’s ratios of phase 2 are manually defined, ranging from -0.95 to 0.45 with two ends, i.e., -0.99 and 0.49, while phase 1 has a fixed Poisson’s ratio.Meanwhile, the prescribed volume fraction of phase 1 is set to be 50%.Other parameters are the same as those in Example I.Tab.1 lists the resulting bulk modulus, HSW upper bound for bulk modulus and their relative percentage differences with typical topologies plotted in Fig.4.It should be noted that whenν2changes from -0.55 to 0.25 as listed in Tab.1,(K2-K1)(G2-G1)≥0.In all the other cases, the composites are not well-ordered.Also, the effective Young’s moduli of resulting composite are drawn in Fig.5.

Fig.4 shows that the obtained bulk modulus well agrees with the HSW upper bounds predicted by Eq.(16).It demonstrates that the proposed optimization method is capable of finding a microstructure with the largest bulk modulus close to the HSW upper bound in both the wellordered or non-well-ordered cases.The curve does not monotonically increase with the increase ofν2.Whenν2=-0.75, the resulting bulk modulus and HSW upper bound have their minimum values.Whenν2approaches 0.45, the resulting bulk modulus experiences a sharp rise and phase 2 becomes totally embedded into the surrounding phase 1 of PBC.The topology of phase 2 changes from Schwartz P structure to spherical shape.Whenν2approaches 0.49, the topology of phase 2 is totally different.

Table 1:The resulting bulk modulus, the HSW upper bound for bulk modulus, and relative difference

Figure 4:Maximizing the bulk modulus for various Poisson’s ratios of phase 2

Fig.5 shows a parabola curve of the resulting effective Young’s modulus for various Poisson’s ratios of phase 2.Parts of the results exceed the Voigt estimation of 1.75.It can be observed that the enhancement effects of the effective Young’s modulus are dependent on the difference between the Poisson’s ratio of two constituent materials, which agrees well with the phenomena in literature [43,44].The optimized topologies in Fig.4 are different with those seeking for maximizing the effective Young’s modulus [44].It can be inferred that the optimized topologies highly influence the mechanical properties of the resulting composites.

Figure 5:The effective Young’s moduli for various Poisson’s ratios of phase 2

5.3 Example III

In this example, the effect of various volume fractions of constitute phase on resulting bulk modulus is investigated.The Young’s moduli of two constituent materials areE1=2.5 andE2=1, respectively, and the Poisson’s ratio areν1=0.0 andν2=-0.5, respectively.Different volume fractions of phase 1 are referred to as the constraints, aiming to find maximized bulk modulus.According to Eq.(15), the HSW upper bound KU is calculated based on various volume fractions.Tab.2 lists the resulting bulk modulus, HSW upper bound and their relative percentage difference.The corresponding PBC topologies are shown in Fig.6.

Table 2:The resulting bulk modulus, the HSW upper bound for bulk modulus, and the relative difference

Fig.6 shows that at the beginning, with a low, phase 1 occupies the corners of the base unit cell.With the increase of, the topology of phase 2 changes from Schwartz P minimal surface to sphere surface.Whenincreases to 65%, phase 2 is totally embedded in phase 1 with a spherical shape.In all cases, the resulting bulk modulus is slightly less than the HSW upper bound.The relative percentage differences are within 1.5%, which also proved HSW upper bound of bulk modulus can be attainable by the proposed method.

Figure 6:Maximizing the bulk modulus for various volume fractions of phase 1

5.4 Example IV

The objective of this example is to obtain the optimal microstructure with maximum shear modulus via the proposed optimization method.The Young’s moduli of constituent phases are selected asE1=5.5 andE2=1.The Poisson’s ratios are assumed to beν1=0.0 andν2=-0.5.The shear moduli for the two materials areG1=2.75 andG2=1.00, respectively.The prescribed volume fraction of phase 1 is 60%.According to Eq.(15), the HSW upper bound is calculated as 1.8710.The final optimized microstructures of orthotropic composites obtained from different initial designs and their effective matrices are presented in Fig.7.

From the effective matrices, the shear moduli of the two designs are 1.8668 and 1.8677,respectively.As discussed previously, the two types of topologies differ significantly due to the distinct initial designs.However, the shear moduli are almost identical and slightly below the respective HSW bounds for shear modulus.

5.5 Example V

The objective of this example is to find the optimal microstructure with a maximum shear modulus.The Young’s modulus of constituent phase 1 is investigated and set to be ranging from 3.0 to 11.0.Other mechanical properties are the same as those in example IV.The prescribed volume fraction of phase 1 is defined as 50%.Tab.3 lists the resulting shear modulus, HSW upper bound for shear modulus and their relative percentage differences.The evolution of the corresponding topologies is plotted in Fig.8.

Fig.8 shows that, with the increase ofE2, the maximizing shear modulus and HSW upper bound rises.The microstructure of phase 2 is a Schwartz P structure whenE2=3.0.The Schwartz P structure and sphere inclusions for each phase coexist in the base cell when the Young’s modulus of phase 2 increases.The configurations of the microstructure are affected by the Young’s modulus of constitute phases.In all cases, the resulting shear moduli well agree with HSW upper bound calculated by Eq.(15).Their relative percentage differences are less than 0.958%.

Figure 7:Base cell and effective elasticity matrix of composite material, starting from different initial designs:(a) From initial design 1; (b) From initial design 2

Table 3:The resulting shear modulus, the HSW upper bound for shear modulus, and the relative difference

Figure 8:Maximizing the shear modulus for various Young’s moduli of phase1

5.6 Example VI

In this example, various volume fractions of the constituent phase on the final design are investigated.The mechanical properties and objective function are the same with those in example IV.The prescribed volume fraction of phase 1 ranges from15% to 85%.Tab.4 lists the resulting shear modulus, HSW upper bound for shear modulus and their relative percentage differences.The corresponding microstructural topologies are depicted in Fig.9.

From Tab.4, all the relative percentage differences between the resulting shear modulus and HSW upper bound for shear modulus are within 3.095%.The results indicated that the theoretical upper bound could be attained under different given volume fraction.From Fig.9, the topology of each phase changes with the increase ofand exhibits different shapes.

Table 4:The resulting shear modulus, the HSW upper bound for shear modulus, and the relative difference

Figure 9:Maximizing the shear modulus for various volume fractions of phase 1

6 Conclusion

In this paper, a topology optimization methodology is presented to design microstructures of materials that are composed of two isotropic phases with distinct Poisson’s ratios.The objective is to maximize the bulk or shear modulus of the composites.The effective elasticity matrix is determined by means of homogenization.Several 3D examples are presented to demonstrate the effectiveness of the proposed method.The examples show that the proposed approach is capable of generating composite microstructures with properties very close to the Hashin-Shtrikan-Walpole upper bounds.In the future research, several important directions furthered the proposed method includes the concurrent design of macrostructure and microstructure containing multiple phases [50].Moreover, it will be of great interest to obtain the smooth boundary directly by variable density method [51].

Funding Statement:This research work was financially supported by the National Natural Science Foundation of Beijing (No.2182067), the Fundamental Research Funds for the Central Universities (2018ZD09).We would like to thank Professor Krister Svanberg for providing the program code of MMA.These supports are gratefully acknowledged.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.