- Email: [email protected]
Topological design optimization of lattice structures to maximize shear stiffness
ARTICLE IN PRESS
JID: ADES
[m5G;May 17, 2017;19:48]
Advances in Engineering Software 0 0 0 (2017) 1–11
Contents lists available at ScienceDirect
Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft
Topological design optimization of lattice structures to maximize shear stiffness Yixian Du a,b,∗, Hanzhao Li a,b, Zhen Luo c, Qihua Tian a,b a
College of Mechanical & Power Engineering, China Three Gorges University, Yichang 443002, China Hubei Key Laboratory of Hydroelectric Machinery Design and Maintenance, Yichang 443002, China c School of Electrical, Mechanical and Mechatronic Systems, The University of Technology, Sydney, NSW 2007, Australia b
a r t i c l e
i n f o
Article history: Received 9 December 2016 Revised 23 March 2017 Accepted 30 April 2017 Available online xxx Keywords: Lattice structure Topology optimization Energy-based homogenization method
a b s t r a c t To improve the poor shear performance of periodic lattice structure consisting of hexagonal unit cells, this study develops a new computational design method to apply topology optimization to search the best topological layout for lattice structures with enhanced shear stiffness. The design optimization problem of micro-cellular material is formulated based on the properties of macrostructure to maximize the shear modulus under a prescribed volume constraint using the energy-based homogenization method. The aim is to determine the optimal distribution of material phase within the periodic unit cell of lattice structure. The proposed energy-based homogenization procedure utilizes the sensitivity filter technique, especially, a modified optimal algorithm is proposed to evolve the microstructure of lattice materials with distinct topological boundaries. A high shear stiffness structure is obtained by solving the optimization model. Then, the mechanical equivalent properties are obtained and compared with those of the hexagonal honeycomb sandwich structure using a theoretical approach and the finite element method (FEM) according to the optimized structure. It demonstrates the effectiveness of the proposed method in this paper. Finally, the structure is manufactured, and then the properties are tested. Results show that the shear stiffness and bearing properties of the optimized lattice structure is better than that of the traditional honeycomb sandwich structure. In general, the proposed method can be effectively applied to the design of periodic lattice structures with high shear resistance and super bearing property. © 2017 Elsevier Ltd. All rights reserved.
1. Introduction The recent developments on large scale computation and advanced manufacturing techniques allow the construction of designs with fine and complex geometrical features [1], e.g. cellular and lattice structures, which have been increasingly investigated for a range of lightweight structural applications, including aerospace, automotive and naval industries due to their high stiffness-to-weight and strength-to-weight ratio, as well as their excellent energy absorption and thermal isolation characteristics [2–6]. As a representation, lattice materials are often assumed to be homogeneous due to it is a number of periodic units and can be efficiently assembled by following a certain pattern. Recently, the honeybee comb structure inspired by natural beeswax has long fascinated engineers and biologists for its outstanding mechanical properties and other characteristics, as a result of their special microstructures such as hexagonal, square and triangular, which are
∗
Corresponding author. E-mail address: [email protected] (Y. Du).
widely used in various engineering applications [7]. However, it has not proven until recently that this is the best structure with enhanced shear resistance that can be built with the beeswax to store the largest amount of honey [8]. Shear response is of great significance as bending lattice structure gives rise to transverse shear loading and may collapse the core of the structure in shear. For instance, the work [2] studied the effective in-plane stiffness of hexagonal honeycomb cores in macro-scale according to the bending model. After that, there is a growing interest in exploiting other materials as the filling material in such structures to simultaneously enhance the shear response, load-bearing and energy absorption capabilities of the lightweight structures, such as periodic lattice structure and honeycomb sandwich cores [9–11]. Many researchers [12–14] studied the mechanical properties of other honeycomb cores such as hollow pyramidal lattices, square honeycombs, and corrugated sandwich cores. Kelsey et al. [15] gave the upper and lower bounds of the shear modulus of honeycomb cores by using the unit displacement and unit load methods in conjunction with a simplification assumption for the strain and stress in the core. Pan et al. [16] and Han et al. [17] analyzed the longitudinal shear strength of honeycomb
http://dx.doi.org/10.1016/j.advengsoft.2017.04.011 0965-9978/© 2017 Elsevier Ltd. All rights reserved.
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
JID: ADES 2
ARTICLE IN PRESS
[m5G;May 17, 2017;19:48]
Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
cores, and investigated the longitudinal shear deformation behavior and failure mechanism of aluminum alloy honeycomb cores at room temperatures. It is demonstrated that the shear or compression strength and specific energy absorption of lattice structures or sandwich cores can be increased dramatically by using appropriate filling foams. Also there have been some papers which studied the shear responses of various periodic lattice structures were mainly focused on the calculation and comparison of their shear modulus and experimental validations [18,19]. However, most of the lattice structures used in engineering are obtained by designers’ experiences or inspired by existing materials or structures in nature. An effective method to guide the design of microstructural topologies of lattice materials is still in demand. Inspired by the above natural-occurred cellular materials, the design of periodic lattice composite to achieve multifunctional properties (stiffness, strength, negative Poisson’s ratio, etc.) based on the combination of structural optimization techniques and numerical homogenization has attracted considerable attention. For instance, Sigmund [20–22] employed a modified optimality criteria method to design material microstructure to achieve the prescribed properties. It can be found that the effective material properties at macro-scale are mainly determined by the topology of periodic microstructures, rather than the proportion and physical properties of their constituents, namely, the intrinsic material composition [23]. The homogenization method [24–25] has been recognized as a rigorous method for characterizing the macro mechanical behavior of materials consisting of periodic microstructures, under the assumption that the linearly elastic response of the periodic material can be determined by test strains over one cell [26]. The topological design of such materials assumes that the material is made of periodic cells which theoretically are infinitesimal within the structure, and the macro effective properties of the heterogeneous material are homogenized according to the microstructure that is the smallest repetitive unit of the material. The inverse design in the material is a typical topology optimization problem, which seeks an optimal microstructure of the cell with prescribed or extreme effective properties [27–36], and so on. A systematic mean of microstructural design is formulated as an optimization problem for parameters that represent the property and topology of the microstructure [37], so as to improve structural shear effect. Over the past, topology optimization has been expanding as a powerful computational design tool for a range of structural and material problems both in academic research and industrial applications. Essentially, topology optimization is a numerical iterative process that distributes a given amount of material inside a reference design domain, so as to seek the best material layout until the expected performance is optimized subject to constraints [38]. Since the work of Sigmund [20–22], various structural topology optimization techniques have been developed for computationally design of microstructures of materials [33,39], e.g. the numerical homogenization method [40], solid isotropic material with penalization (SIMP) [41–43], level set method (LSM) [36], parametric level set method (PLSM) [44], and bidirectional evolutionary structural optimization (BESO) [45]. Amongst them, the design multifunctional composite materials are an active field of research [46]. Huang et al. [47] recently applied the BESO approach to the design of 2D and 3D materials with extreme bulk and shear modulus, and isotropic constraint is imposed by Radman [44]. Zhang et al. [48,49] and Xia et al. [50] proposed the strain energy-based method to predict the effective elastic properties, which compared with the numerical homogenization method and shows the advantages of higher computing efficiency and simplicity in the numerical implementation. Based on the above work, we can find that there is still much room to improve the shear resistance of periodic lattice structure
consisting of hexagonal unit cells in engineering applications. Lattice structures composed of unit cells are one of the candidates to enhance structural/material performances. However, the topological design problem of periodic lattice structure in recent works are mostly generalized by using the mathematical homogenization method, which has the disadvantages, such as the numerical implementation of the homogenization method and design sensitivity analysis are complicated and additional programming is needed. In this paper, the topology optimization technique is utilized to search the layout of periodic lattice structure with enhanced shear stiffness under the constraint of a given amount of material by using energy-based homogenization method. The optimization objective is to find the optimal distribution of two base materials within the unit cell of the periodic lattice material, so that the resulting structure has the maximum shear modulus. The effective properties of the microstructure are totally predicated by the periodic unit cell (PUC) which is discretized into finite elements under periodic boundary condition. 2. Energy-based homogenization method 2.1. Formulation and sensitivity analysis with energy-based homogenization In general, the homogenization method applied to linear elastic problems establishes macroscopic properties that often are uniformly described by corresponding PUCs, which are periodically repeated in one or more directions. The properties of a heterogeneous medium rely on the analysis of its PUCs. Thus, the homogenization theory can be used to calculate the effective elastic constants of the macroscopic composite by the Hill-Mandel condition or energy averaging theorem [51]. To obtain the macroscopic equivalent constitutive properties by the homogenization process, three widely used types of loading can be applied: (a) prescribed linear displacements, (b) prescribed tractions and (c) periodic boundary conditions. The size of the periodic unit cell is assumed to be much smaller than the size of the bulk material, so that the average quantities of both the microscopic strain and stress can be defined [51]. Consider a single cell Y, the macro-scale displacement field μ is expanded, depending on the little aspect ratio ε between the macro and micro scales by using the asymptotic homogenization. When only the first-order terms of the asymptotic expansion are considered, the effective elasticity tensor of the macro-material can be found in terms of the material distribution in the domain of PUC by averaging the integral over the base cell [50] Y as:
CiHjkl =
1 |Y |
Y
ε 0pq(i j ) − ε ∗pq(i j ) Cpqrs ε 0pq(kl ) − ε ∗pq(kl ) dY
(1)
where |Y| represents the area (or volume) of the unit cell, and ε 0pq(kl ) defines the prescribed macroscopic strain fields, which are defined as three linearly independent unit strains: the horizontal unit strain ε 0pq(11) = [1 0 0], the vertical unit strain ε 0pq(22) = [0 1 0] and the shear unit strain ε 0pq(12) = [0
0 1]. The local
strain field ε ∗pq(kl ) induced by the test strains are defined as follow:
1 kl kl ε ∗pq(kl ) = ε ∗pq χ kl = χ p,q + χq,p 2
(2)
where χ kl denotes a Y-periodic admissible displacement field associated with the load case kl, which can be obtained from the following equilibrium equation [50]:
∗ kl ∂ vi Ci jpq ε pq( ) dY = ∂yj Y
Y
0 kl ∂ vi Ci jpq ε pq( ) dY ∂yj
(3)
where v is the Y-periodic admissible displacement field.
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
ARTICLE IN PRESS
JID: ADES
[m5G;May 17, 2017;19:48]
Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
The stress and strain tensors of the homogeneous medium are equivalent to the average stress and strain of the microstructure, thus, the periodic microstructure of composite materials can be replaced by an equivalent homogeneous medium with the same volume at the macroscopic level [49]. The average stress and strain of the homogeneous medium follow the Hooke’s law, which can be rewritten in matrix form, as follows
⎡σ ⎤
⎡
11
H C1111
H ⎣σ22 ⎦ = ⎢ ⎣C1122 σ12 0
H C1122 H C2222
0
⎤⎡ ⎤ 0 ε11 ⎥⎣ ⎦ 0 ⎦ ε22 H ε12 C1212
(4)
where σ11 , σ22 denotes the horizontal and vertical averaged stress. ε11 , ε22 denotes the averaged strain in corresponding two directions. σ12 , ε12 denotes the averaged shear stress and strain in two directions. In fact, the strain energy-based method is equivalent to the numerical homogenization method in predicting the effective properties of materials [48]. Considering the shear-strain state, which are ε11 = ε22 = 0, ε12 = 1, the unit strain energy of the microstructure is expressed as Et = σ12 ε12 /2 = σ12 /2. The strain energy of microstructure under periodic boundary condition is expressed by
Et =
1 1 H F U = U T KU = C pqrs ε Apq(i j ) εrsA(kl ) 2 2
(5)
where F, U and K are the load vector, the displacement vector of element nodes and the stiffness matrix, respectively. The energy-based approach employs average stress and strain theorem, in which the unit test strains are imposed directly on the boundaries of the base cell instead of the asymptotic approach, including ε Apq(kl ) which corresponds to the superimposed strain fields
(ε 0pq(kl ) − ε ∗pq(kl ) ) in Eq. (1).
In finite element analysis, the base cell is discretized into N finite elements and the shear modulus is obtained in an equivalent form in terms of element mutual energies H C1212 =
σ12 = 2Et = ε12
Y
A (i j ) A i j A kl C pqrs ε pq( ) εrs( ) dY = ue N
T
A kl ke ue ( )
e=1
(6) uAe (kl )
where are the element displacement solutions corresponding to the mentioned unit test strain fields ɛ0( kl ) , and ke is the element stiffness matrix. The effective elastic constant of the composite highly depends on the spatial distribution of material phases, therefore, how to optimally distribute material phase within the PUC would be of high importance. Thus, microstructures can be modeled as 2D continuum structures. An artificial design variable, x, is introduced to represent material occupation in an element by assuming that x = 1 if an element is made of solid material and x = 0 versus. The element property is interpolated using the SIMP model, where the density-stiffness interpolation scheme [41] can be stated as
Ee (x ) = Emin + (ρe (x ) ) p (E0 − Emin )
(7)
where ρ e (x) designates the relative density in the random point x of the design domain, which takes values between 0 and 1. E0 is the Young’s modulus of solid material and Emin is a very small stiffness assigned to void regions for approximation of void material (namely the Ersatz material or compliant material) in order to prevent the stiffness matrix from becoming singular. p is a penalization factor (typically p = 3–4) introduced to ensure the density distribution closer towards the black-and-white solutions. Here, topology optimization of the material microstructure aims at finding the optimal microstructure of the PUC that exhibits the maximum shear stiffness in specified directions under prescribed
3
periodic boundary conditions. The mathematical formulation of the optimization problem can be stated as
⎧ H maxρ : f CiHjkl (ρ ) = C1212 ; ⎪ ⎪ ⎪ ⎪ A(kl ) (kl ) ⎪ ⎪ ⎨s.t. : KU = F , k, l = 1, · · · , d N
⎪ ve ρe /|Y | ≤ ϑ ⎪ ⎪ ⎪ ⎪ e=1 ⎪ ⎩ 0 < ρe ≤ 1 , e = 1 , · · · , N
(8)
where the objective f (CiHjkl (ρ ) ) is a function of the homogenized stiffness tensor, namely the maximization of shear modulus mentioned in Eq. (6). K is the global stiffness matrix, UA ( kl ) and F( kl ) are the global displacement vector and the external force vector of the test case (kl), respectively. ve denotes the element volume, ϑ is the limited upper bound of volume fraction, and d is the spatial dimension. It is critical to acquire gradient information to guide the algorithm to search for an optimum point which is most efficient during the iteration process, which is an important step in topology optimization. In this regard, topological sensitivity is frequently defined as the derivative with respect to design variable, which is the relative density of element in the density-based finite element framework. So, the sensitivity of the strain energy with respect to the element density variable is obtained as
∂ Et ∂ 1 T 1 T ∂K ∂U T = U KU = U U +2 KU ∂ ρe ∂ ρe 2 2 ∂ ρe ∂ ρe
where
∂U T ∂ UT KU = ∂ ρe ∂ ρe
∂ U T ∂ ρe
F
F
= 0
∂ U T ∂ ρe
(9a)
F 0
=0
(9b)
when the duality condition of displacement and force on the boundary are considered, here U and F are the nodal displacement vector and force vector on the boundary. U and F are the nodal displacement vector and force vector inside. The periodic boundary condition mentioned above impose three unit test strains directly on the boundaries of the base cell, the specific nodal displacement vector have no relationship with density variable, so the derivative of nodal displacement on the boundary with respect to density equals to zero. What’s more, the force vector inside the unit cell indeed equals to zero to meet the continuity and periodicity requirements for the displacement as well as stress in the FEA. With the help of the above mentioned material interpolation, the first-order derivative of the strain energy with respect to the element density is given by:
∂ Et 1 ∂K p n = UT U= E ∂ ρe 2 ∂ ρe ρe t
(9c)
It is easily to get the sensitivity of the homogenized elasticity H tensor C1212 with respect to ρ e , which can be derived from the adjoint method [41] as: H A 12 T A 12 ∂ C1212 1 ∂ Et p = pρ p−1 (E0 − Emin ) ue ( ) ke ue ( ) = 2 = 2 Etn ∂ ρe ∂ ρe ρe |Y | e
(10) where Etn is the element strain energy, which effects the sensitivity of the effective elastic tensor. In order to avoid the formation of checkerboard pattern and the mesh-dependency in topology optimization, the filtering scheme [41] is used during the process of solving the optimization problem, which modified the design sensitivities during iterations, as
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
ARTICLE IN PRESS
JID: ADES 4
Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
1 0.8 The update density
As indicated in Fig. 1, the polarization operator can reasonably make the intermediate density compelling to the two ends of 0– 1, which ensure the optimal solution approaching to 0/1 discrete properties. The curves are become steeper and steeper with the polarization factor’s increasing. The proposed algorithm is just like a smoothed Heaviside function that can polarize the majority of elements to 0 or 1.0 more rapidly at the same probability. The number of elements with the intermediate density in the final topology will be reduced, so as to get a topological design with more clear edges. The standard optimality criteria (OC) method has been widely recognized as an effective optimization algorithm for its many advantages while used, such as fast convergence, simplified programming and with no association with the structural re-analysis and the number of the variables [41]. Combine the polarization operator with the optimality criteria equation, a new iterative updating scheme is expressed as follows:
s=9 s=11 s=13 s=16
0.9 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 The initial density
1
Fig. 1. Optimization with different polarization factors.
k xe = xne Be( ) n
η
follows
∂ Cˆ 1 ∂C = Hi ρi ∂ ρe max(γ , ρe ) i∈Ne Hi ∂ ρi
(11)
where Ne is the set of elements i for which the center-to-center distance to element e is smaller than the filter radius rmin and Hi is a weight factor defined as: Hi = max(0, rmin - (ei)), here the term γ = 10−3 is a small positive number introduced in order to avoid division by zero. The filter radius of rmin is defined to identify the neighboring elements that affect the sensitivity of the current element. 2.2. Modified optimization algorithm Due to the relaxation of original 0/1 discrete problem to a continuous 0–1 optimization problem, the numerical problems of gray regions, checkerboards and mesh dependence mostly exist in the topology optimization. The appearance of such geometrical features in either solid or the void regions may lead to lack of manufacturability in engineering applications. Thus, to more effectively compress gray areas in the final topological design and achieve the topological design with clear boundary, several different methods have been developed [52,53], such as new density interpolation function, extra constraints, postprocessing techniques and filtering schemes [54]. In this paper, we proposed an optimization algorithm with the approximation of 0/1 discrete properties using a heuristic polarization operator. The basic idea of this method is that it can weaken the low-relativedensity elements while enhance the high-relative-density elements by the polarization operation during the iteration process. The polarization further makes those intermediate density elements towards 0 or 1 discrete design, so that the final solution has the discrete properties. In this way, we can get better design in topology optimization with lower gray elements and computational efficiency. The specific mathematical expression of polarization operator is defined as follows:
x,e
[m5G;May 17, 2017;19:48]
=
exp − 2s + s × xe
1 + exp − 2s + s × xe
(12)
where xe is the relative element density before polarization by polarization operator, xe ’ denotes the relative element density after polarization, s is the polarization factor that controls the speed of the polarization. Fig. 1 shows the different shapes when the polarization factor in Eq. (12) have different values.
=
k exp − 2s + s × xne Be( )
η
k 1 + exp − 2s + s × xne Be( )
η
⎧ n n n ⎪ ⎨ max (xmin , xe − m) xe < max (xmin , xe − m) xne +1 = xne max (xmin , xne − m ) < xne < min (1, xne + m) ⎪ ⎩ min (1, xne + m) xne > min (1, xne + m)
(13)
(14)
η
where [xne (Be(k ) ) ] is the not-polarized solution of the polarization operator before the process of solution, namely the initial solution; η
[xne (Be(k ) ) ] is the polarized solution of the polarization operator before the process of solution, namely the update solution; s is the polarization factor, which is used to control the punishment speed of the elements that belong to the middle part, here we choose s = 7 in this paper; xne is the updated solution by using the optimality criteria method. 2.3. Numerical example Assume that the Young’s modulus and Poisson’s ratio of the artificial solid phase are selected as E = 1, and μ = 0.3, the given volume fraction of solid material is set to be 0.4 and 0.6, the PUC is discretized with a mesh of 50 × 50 and 100 × 100 four-node quadrangle elements, and the penalty factor p and the filter radius rmin are endowed with different values, the corresponding microstructural topologies, material distribution and effective elasticity matrices are given in Table 1. From the optimal solutions and the corresponding effective matrices for different parameters cases, it can be seen that these different volume fractions, the mesh, the filter radius and the penalty factors all have influence on the final topological shapes of the unit cell with the same objective function. Under the same volume fraction, the equivalent elasticity matrix of the optimal microstructure is almost the same. This is reasonable due to the fact that a number of different microstructures can possess the same physical property. At the same time, undesirable multiple local minima that occurred in deleting the intermediate density elements during the searching process are avoided in the optimization procedure. Moreover, with the help of sensitivity technique, the number of elements with the intermediate density in the final topology will be reduced, so as to get a topological design with clear boundaries. It can be seen from Table 1 that the final optimal configuration is symmetrical, and the material layout is commonly in line with diagonal directions which is 45°. Thus, the potential of the material is motivated to resist shear loads so that the optimal microstructures exhibit high shear stiffness. Fig. 2 shows the iteration histories of the shear modulus, the unit cell is discretized by 100 × 100 four-node quadrilateral finite
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
JID: ADES
ARTICLE IN PRESS
[m5G;May 17, 2017;19:48]
Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
5
Table 1 Microstructural topologies, material distribution and effective elasticity matrix of cellular materials with maximum shear modulus for various volume fractions, mesh sizes, filtering radius and penalty factors.
elements, volume fraction is 0.5, filter radius is 5 and penalty factor is 3, respectively. It can be seen from Fig. 2 that when the update scheme given in Eqs. (13) and (14) is used, the shear modulus increases rapidly at the first 27 iterations, and then the curve tends to stabilize at iteration 140, the topological details change slightly. After 230 iterations, the objective function and the microstructural topologies stably converge to the final solutions. The distribution of design variables within the unit cells (Fig. 3) can be interpreted as octagonal honeycomb cell whose relative shear modulus equals to 0.106.
According to the microstructure of unit cells shown in Fig. 3, the 3 × 3 base cells are constructed periodically which can be called periodic lattice structure as shown in Fig. 4. When the design problem formulation of the periodic lattice structure has been done, the obtained optimal structure can be regarded as one octagonal honeycomb cell with a nested diamond-style structure inside, and this structure are almost the same as the solutions by Huang [47]. Also the lattice structure can be interpreted as the modified version of hexagonal honeycomb which are assembled by four cells in a certain pattern. In the next section, to show the good perfor-
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
ARTICLE IN PRESS
JID: ADES 6
[m5G;May 17, 2017;19:48]
Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
Fig. 2. Iteration histories of shear modulus and microstructure using energy-base homogenization.
effective in-plane elastic stiffness of honeycomb cores. While this analysis model ignore the effect of transverse shear stress and the different thickness of the unit cell wall. Therefore, the analytical analysis model is not used for the designed lattice structure directly. The effective elastic constants of designed periodic lattice structure’s core are re-derived by considering the transverse shear deformation and the flexibility of the cell wall. To verify the shear stiffness of the lattice structure, the equivalent elastic constants are deduced and compared with those of the hexagonal honeycomb sandwich structure. 3.1. Mechanical equivalent model of periodic lattice structure
Fig. 3. Distribution of design variables.
The simplified mechanical model of unit cell is built in Fig. 5 to analyze its equivalent mechanical performance. In Fig. 5, t denotes the thickness of the cell unit (mm), b is the height of the lattice structure (mm), l denotes the length of the longitudinal wall (mm), h represents the length of the inclined wall (mm), Es and Gs are Young’s modulus and shear modulus of base material of periodic lattice structure (MPa), θ denotes the angle between the hexagon short edge and horizontal side (°), respectively. Using force balance principle, the cantilever beam model and the material mechanics beam bending theory, the equivalent mechanical parameters of the lattice structure are presented, as is shown in Fig. 5. 3.2. Equivalent shear modulus analysis
Fig. 4. 3 × 3 unit cells.
mance of the lattice structure with shear stiffness, the equivalent mechanical properties are studied in macro-scale based on the theoretical model. 3. Theory analysis of periodic lattice structure According to different geometries of the microstructure, the topological configuration of two-dimensional periodic lattice structure may vary widely, but the general analytical formulation to characterize the equivalent mechanical properties is similar. Gibson developed the classical unit cell theory [2] for predicting the
Considering the symmetry of the most basic hexagon and trapezium composition unit of the lattice structure, the loads, displacements, strains and stresses on cell walls subjected to transverse shear displacement is shown in Fig. 5. The mechanical behavior of the inclined wall can be investigated based on the cantilever beam model due to the fact that the bending deformation of the inclined wall is anti-symmetric along the height direction just like the cantilever beam. The bending displacement of the inclined wall is supposed to be δ 1 , the equivalent load P1 can be written as follows:
P1 =
(δ1 /2 ) × 3Es I1 (b/2 )3
(15)
where I1 is the moment of inertia of the inclined wall, which can be expressed as I1 = ht13 /12, t1 is the thickness of the inclined wall, which can be expressed as t1 = 2t.
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
ARTICLE IN PRESS
JID: ADES
[m5G;May 17, 2017;19:48]
Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
7
Fig. 5. Simplified equivalent model of periodic lattice structure.
Suppose the deformation under transverse shear displacement for the longitudinal wall consists of two parts: the bending displacement δ 2 and shear displacement δ s . The corresponding equivalent loads are written as
P2 =
(δ2 /2 ) × 3Es I2 (b/2 )3
Ps = Gs δs lt/b
(16) (17)
where I2 is the moment of inertia of the longitudinal wall, which can be expressed as I2 = lt3 /12. The shear displacement is supposed to be δ y in the transverse direction like the x-direction shown in Fig. 5, the bending displacement of inclined wall, the shear displacement and bending displacement of longitudinal wall can be expressed as: δ 1 = δ y , δ s = δ y × cos θ , and δ 2 = δ y × sin θ . In order to obtain the precise expression of the transverse shear modulus of the lattice structure, the different thickness, the bending contribution of the inclined wall and the longitudinal wall should be considered. The mechanical equilibrium condition of the representative unit requires eq P1 + 2P2 sin θ + 2Ps cos θ = 2τxy (h + l sin θ )l cos θ
(18)
where τxy is the equivalent shear-stress that imposed on the equivalent area of one representative unit cell which can be exeq eq eq eq pressed as τxy = Gxy γxy , γxy is the equivalent shear strain that imposed on the equivalent area of one representative unit cell which eq can be expressed as γxy = δy /b. Thus, substituting Eqs. (15), (16) and (17) into Eq. (18) yields the following equation that computes the equivalent shear modulus of the periodic lattice structure: eq
Geq xy =
Gs t cos θ 4ht 3 + lt 3 sin2 θ Es + h + l sin θ (h + l sin θ )l cos θ b2
(20)
The mass of the equivalent entity plate may be expressed as
(21)
According to the conservation of mass before and after the equivalence (mce = m1 ), we have
ρc =
4t (l + h )ρs l cos θ (l sin θ + h )
⎧ Gs t cos θ 4ht 3 + lt 3 sin2 θ Es ⎪ + ⎨Geq xy = h + l sin θ h ( + l sin θ )l cos θ b2 ⎪ 4 ρs t ( β + 1 ) ⎩ ρc = l (β + 2 sin θ )(β + 2 cos θ )
(23)
eq
where Gxy is the equivalent shear modulus of lattice structure in xy-plate, ρ c and ρ s are the equivalent densities of the lattice structure and the base material of the periodic lattice structure. 3.4. Comparison of two structure’s shear moduli
Based on the conservation of mass, the mass of the representative unit cell wall can be calculated at first, and then the structure is equivalent to an entity plate with the same size. The equivalent mass of the representative unit cell of the periodic lattice structure may be expressed as
mce = ρcVce = ρs · b(h + 2l sin θ )(h + 2l cos θ )
As a result, the equivalent mechanical parameters of the periodic lattice structure are presented as follows:
(19)
3.3. Equivalent density analysis
m1 = ρsV1 = ρs × 4bt(h + l )
Fig. 6. Simplified equivalent model of honeycomb core by Gibson [55].
(22)
A case of bearing structure that is made of 2024 aluminum alloy material is used. The thickness of the hexagonal honeycomb sandwich structure and the periodic lattice structure is given the same value. Fig. 6 shows the Gibson’s analysis model under shear deflection of the honeycomb core. The parameters are supposed as: θ = 45˚, h = l, and β = h/l = 1, thus, Eq. (23) can be simplified as follows:
⎧ 9Es t 3 + Gs t b2 eq ⎪ ⎪ ⎨Gxy = √ 2 2 + 1 lb
8 ρs t ⎪ ⎪ √ ⎩ ρc =
(24)
l 3+2 2
According to Gibson’s equation [55], the thickness of the representative unit cell of hexagonal honeycomb sandwich is assumed to t. The corresponding equivalent analysis model and parameters
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
ARTICLE IN PRESS
JID: ADES 8
[m5G;May 17, 2017;19:48]
Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
P
Fig. 7. Simplified load condition and constraints.
Fig. 8. FE models of the two lattice structures.
are shown in Fig. 6. When θ = 30˚, h = l, β = 1 is choosing, the corresponding equation can be simplified as
⎧ Es t 3 ⎪ ⎨Gacxy = √ 3 3l
Table 2 Size parameters of periodic lattice structure.
(25)
ρs t ⎪ ⎩ ρac = 8√
Cell size of core (mm)
t 0.05
l 5
h 10
Base material (aluminum)
Es (GPa) 70
vs 0.26
2780 kg/m3
ρs
3 3l
In order to intuitively compare the shear modulus of the two structures theoretically under the condition that the mass of the entity of structures are the same, the equivalent density is also assumed to√be the same value, namely, ρ c =ρ ac , thus, ρc /ρac = √ 3 3/(3 + 2 2). The ratio between the shear modulus of two structures is defined as nG , which can be calculated by
Geq xy nG = = Gacxy
√ 3 9t 2 + GEss b2 l 2 √ 2 + 1 b2 t 2
(26)
In the above equation, the Young’s modulus and the shear modulus of the same isotropic material have a relation of Gs = Es / 2(1 + vs ). Furthermore, l2 / t2 >> 1 can be easily obtained according to the application of engineering. As a result, Eq. (26) can expressed as follows:
nG =
√ √ 3+2 2 9 3t 2 l2 √ + 1.86 2 ≈ 1.86 × ≈ 2.09 > 1 √ t 3 3 2 + 1 b2
(27) eq Gxy
From above equation, it is evident that > Gacxy can be easily obtained. In other words, the shear modulus of periodic lattice structure is highly increased when considering the same mass and the same application condition. It also may be concluded that, the proposed method can be effectively applied to enhance the shear stiffness of periodic lattice structure theoretically. 4. Finite element analysis and experimental investigation 4.1. Finite element analysis The FE simulation analysis for the mechanical performance of the periodic lattice structure and the honeycomb sandwich structure is given as follows [56]. The two structures are made of 2024
aluminum alloy materials, with yield strength 758 MPa and material density 2780 kg/m3 , and other data are shown in Table 2. To compare the shear performance and bearing ability of this two structures, the FE models of the two structures were given in Fig. 8 by using the same type of elements (SHELL63 plate-like elements in the three dimensional space) under the same geometry and load conditions, as shown in Fig. 7. Both the face-sheets were taken as rigid with no deformation. As shown in Fig. 7, the periodic lattice structure assembled as a plate, the two sides of the plate were constrained, which means all degrees of freedom were clamped. At the top surface, each node was subjected to a constant shearing force (linear elastic analysis) along the thickness direction. The initial value of applied load is 0.005 MPa, and the incremental value is set to be 0.005 MPa as well, 10 sets of numerical experiments were done for the purpose of a better comparison under different sets of comparison and analysis. According to mechanics, the stress-strain relation of periodic lattice structure can be used to calculate the shear modulus as shown in Eq. (19). Table 3 shows the results of FEM simulation and theoretical analysis, for the two structures with the same geometry parameters and simulation conditions. The theoretical results and numerical results are expressed as Ti and Ni respectively under the load case i and the relative error can be calculated by:
error = i
Ni − Ti Ti
(28)
It can be seen from Table 3 that the maximum relative error between the theoretical and FEM numerical results does not exceed 5%. Results from the proposed analytical model well align with the FE predictions for the periodic lattice configurations in this study. Therefore, it is clear that the proposed mechanical equivalent model for calculating the equivalent elastic modulus are rea-
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
ARTICLE IN PRESS
JID: ADES
[m5G;May 17, 2017;19:48]
Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
9
Table 3 The shear modulus of periodic lattice structures under various loads. Loads(MPa)
0.005
0.010
0.015
0.020
0.025
0.030
0.035
0.040
0.045
0.050
Theoretical results(MPa) Numerical results(MPa) The relative error
1.997 2.030 1.7%
1.997 2.002 0.3%
1.997 1.992 −0.3%
1.997 1.988 −0.5%
1.997 1.986 −0.6%
1.997 1.984 −0.7%
1.997 1.980 −0.9%
1.997 1.979 −0.9%
1.997 1.977 −1.0%
1.997 1.976 −1.0%
3.5
Nominal strain / mm
3
Periodic lattice structure Hexagonal honeycomb structure
2.5 2 1.5 1 0.5 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035
0.04 0.045 0.05 0.055
Load / MPa Fig. 9. The strains of two structures under various loads.
sonable and effective. In order to obtain a more precise solution, the transverse shear deformation and the uniform thickness of the cell wall should be considered. The strain for the honeycomb sandwich and the periodic lattice structures under various loading conditions are listed in Fig. 9. It can be found from Fig. 9 that when the loading conditions are the same the strain of the honeycomb structure is larger than that of the periodic lattice structure. In other words, the shear resistance and bearing property are better than that of the honeycomb structure. The difference of strains between the honeycomb structure and the periodic lattice structure increase along with the increase of the applied loads. It indicates that the increase of the shear performance of the lattice structure is more obvious when the loads large enough. In a word, the designed periodic lattice structure shows a better shear resistance and bearing performance, compared to the widely used honeycomb sandwich structure in engineering. 4.2. Experimental investigation According to the results by using topology optimization technique in Section 2, the periodic lattice structure given as below has been greatly simplified, according to the manufacturing process. Our specimens have been manufactured in polyamide using the selective laser sintering (SLS) and Wire Saw Cutting. The specimens manufactured by two means are given in Fig. 10 respectively. The sample manufactured by SLS is fragile and it is difficulty to be used for practical experimental test. When the manufacturability considered, the experimental specimen was manufactured in aluminum following two steps: 1) prepare enough bars manufactured by using wire cutting technique, 2) assembly or paste the bars by using bonding technology. The thickness of the cell unit is t = 0.5 mm, The length L, width B and height b are 112 mm × 67 mm × 10.6 mm, separately. The shear tests have been done at room temperature. The specimens were tested under a constant cross-head movement of 0.1 mm/min by using the shearing equipment according to ASTM
Fig. 10. (a) Specimen manufactured by SLS. (b) Specimen manufactured by wire cutting technique.
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
JID: ADES 10
ARTICLE IN PRESS
[m5G;May 17, 2017;19:48]
Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
Fig. 11. (a) Schematic of test specimen. (b) Force-displacement of periodic lattice structure.
C273-00 [57]. The typical force-displacement curves obtained during the shear deformation of the periodic lattice structure were analyzed based on the observation of microstructural change during the deformation. Fig. 11(a) shows the shear test equipment and Fig. 11(b) shows a typical force-displacement curve obtained from shear test and a schematic microstructural change during the shear deformation of specimen. From Fig. 11, it can be seen that the shear force increased almost linearly with the increase of displacement at the beginning stage due to the elastic buckling of cell walls, not the elastic axial shortening of cell walls. The thin cell walls of periodic lattice structure are restrained with the neighboring cell walls. Therefore, when the shear deformation is performed to the lattice structure, it is difficult to be happened by the elastic axial shortening of cell walls. When the force-displacement curve passes a maximum stress, a maximum shear strength is shown in Fig. 11(b), then the shear stress rapidly dropped after the first stage due to the debonding fracture at cell interfaces. A continuous fracture of phenolic resin layer on cell walls is observed. During proceeding the plateau deformation from shear stress-strain curve, the fracture of specimen by the core interfaces debonding occurs finally with a decrease of shear stress on the periodic lattice structure. In the procedure, the deformation behavior of the cell walls during shear test was analyzed as shown in Fig. 11, which shows the basic cell carrying out the shear stress in periodic lattice structure. The arrows around cell walls indicate the shear force performed on the cell wall during shearing in transverse direction. The slope P / δ of the force-displacement curve during elastic deformation shown in Fig. 11(b) is regarded as the total shear modulus of the structure ( P, δ is selected to be 23.0414 and 0.016 as shown in Fig. 11), which are substituted into the shear performance test formula [57] to calculate the shear modulus as follows
Gc =
23.0414 b · P 10.6 = · = 2.03MPa L · B · δ 112 × 67 0.016
structures with the transverse shearing test, which demonstrates the excellence shearing performance of the periodic lattice structure. 5. Discussion and conclusions This study has developed a topology optimization method to enhance the transverse shear performance of the hexagonal lattice structure. The optimization problems have been formulated to generate the periodic microstructures with maximum shear modulus under a prescribed volume and constraint with the strain energybased homogenization method. An equivalent analytical solution for the transverse shear stiffness has been proposed according to the simplified optimal results by considering the transverse shear deformation and the uniform thickness of the cell wall. The shear stiffness of the designed lattice structure has been studied by using the theoretical approach and the finite element method, which demonstrates the effectiveness of the method used in this paper. Finally, the structure is prototyped, and then the shear properties are tested and compared with those of the widely used honeycomb sandwich structures in engineering. Results show that the proposed method can be effectively applied to the design of periodic lattice structures with higher shear resistance and better load bearing abilities. Acknowledgments The work is supported by the National Natural Science Foundation of China (51105229, 51475265), Science and Technology Support Program of Hubei Province of China (2015BHE026), the Research Project of Hubei Provincial Department of Education (D20161205).
(29)
The relative error between the experimental results and the theoretical results of shear modulus shown in Table 3 is 2%, which is less than 5%. That is to say, results from the experimental results are well in line with the finite element prediction and theoretical results as given in Table 3 for the periodic lattice structure. This demonstrates the correctness of the equivalent analysis and FEA model. What’s more, the measured shear modulus of the lattice structure in the transverse direction is compared with those of honeycomb structures in reference [18] under the same experimental conditions. As a result, the measured shear modulus of the lattice structure is 1.49 times higher than those of honeycomb
References [1] Zheng X, Lee H, Weisgraber TH, Shusteff M, DeOtte J, Duoss EB, et al. Ultralight, ultrastiff mechanical metamaterials. Science 2014;344(6190):1373–7. [2] Li H, Luo Z, Zhang N, Gao L, Brown T. Integrated design of cellular composites using a level-set topology optimization method. Comput. Meth. Appl. Mech. Eng. 2016;309:453–75. [3] Fleck NA, Ashby MF. Micro-architectured materials: past, present and future. Proc. R. Soc. A 2010;466(2121):2495–516. [4] Wang Y, Luo Z, Zhang N, Wu T. Topological design for mechanical metamaterials using a multiphase level set method. Struct. Multidiscip. Optim. 2016;54(4):937–52. [5] Mines RAW, Tsopanos S, Shen Y, Hasan R, McKown ST. Drop weight impact behaviour of sandwich panels with metallic micro lattice cores. Int. J. Impact Eng. 2013;60(10):120–32.
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
JID: ADES
ARTICLE IN PRESS Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
[6] Montemurro M, Catapano A, Doroszewski D. A multi-scale approach for the simultaneous shape and material optimisation of sandwich panels with cellular core. Composites Part B 2016;91:458–72. [7] Zhang K, Duan H, Karihaloo BL, Wang J. Hierarchical, multilayered cell walls reinforced by recycled silk cocoons enhance the structural integrity of honeybee combs. Proc. Natl. Acad. Sci. U. S. A. 2010;107(21):9502–6. [8] Mackenzie D. Proving the Perfection of the Honeycomb. Science 1999;285(5432):1338–9. [9] Han B, Yu B, Xu Y, Chen CQ, Zhang QC, Lu TJ. Foam filling radically enhances transverse shear response of corrugated sandwich plates. Mater. Des. 2015;77:132–41. [10] Ostos JB, Rinaldi RG, Hammetter CM, Stucky GD, Zok FW, Jacobsen AJ. Deformation stabilization of lattice structures via foam addition. Acta Mater. 2012;60(19):6476–85. [11] Zhang J, Supernak P, Mueller-Alander S, Wang CH. Improving the bending strength and energy absorption of corrugated sandwich composite structure. Mater. Des. 2013;52(24):767–73. [12] Côté F, Deshpande VS, Fleck NA, Evans AG. The compressive and shear responses of corrugated and diamond lattice materials. Int. J. Solids Struct. 2006;43(20):6220–42. [13] Lebée A, Sab K. Transverse shear stiffness of a chevron folded core used in sandwich construction. Int. J. Solids Struct. 2010;47(18-19):2620–9. [14] Shi G, Tong P. Equivalent transverse shear stiffness of honeycomb cores. Int. J. Solids Struct. 1995;32(10):1383–93. [15] Kelsey S, Gellatly RA, Clark BW. The Shear modulus of foil honeycomb cores: a theoretical and experimental investigation on cores used in sandwich construction. Aircraft Eng. Aerosp. Technol. 1958;30(10):294–302. [16] Pan SD, Wu LZ, Sun YG. Transverse shear modulus and strength of honeycomb cores. Compos. Struct. 2008;84(4):369–74. [17] Han B, Qin KK, Yu B, Zhang QC, Chen CQ, Lu TJ. Design optimization of foam-reinforced corrugated sandwich beams. Compos. Struct. 2015;130:51–62. [18] Qiao P, Xu XF. Refined analysis of torsion and in-plane shear of honeycomb sandwich structures. J. Sandwich Struct. Mater. 2005;7(4):289–305. [19] Li X, Li G, Wang CH, You M. Optimum design of composite sandwich structures subjected to combined torsion and bending loads. Appl. Compos. Mater. 2011;19(3-4):1–17. [20] Sigmund O. Materials with prescribed constitutive parameters: an inverse homogenization problem. Int. J. Solids Struct. 1994;31(17):2313–29. [21] Sigmund O. Tailoring materials with prescribed elastic properties. Mech. Mater. 1995;20(4):351–68. [22] Sigmund O. New class of extremal composites. J. Mech. Phys. Solids 20 0 0;48(2):397–428. [23] Wang Y, Gao J, Luo Z, Brown T, Zhang N. Level-set topology optimization for multimaterial and multifunctional mechanical metamaterials. Eng. Optim. 2016;49(1):22–42. [24] Sanchez-Palencia E. Non-homogeneous media and vibration theory. Lect. Notes Phys. 2012;127(127). [25] Bensoussan A, Lions JL, Papanicolaou G. Asymptotic analysis for periodic structures. North-Holland Pub. Co. 1978;46(2):392–402. [26] Cramer AD, Challis VJ, Roberts AP. Microstructure interpolation for macroscopic design. Struct. Multidiscip. Optim. 2016;53:489–500. [27] Swan CC, Arora JS. Topology design of material layout in structured composites of high stiffness and strength. Struct. Multidiscip. Optim. 1997;13(1):45–59. [28] Paulino GH, Silva ECN, Le CH. Optimal design of periodic functionally graded composites with prescribed properties. Struct. Multidiscip. Optim. 2009;38(5):469–89. [29] Zhou S, Li Q. Computational design of multi-phase microstructural materials for extremal conductivity. Comput. Mater. Sci. 2008;43(2):549–64. [30] Guest JK, Prévost JH. Design of maximum permeability material structures. Comput. Meth. Appl. Mech. Eng. 20 07;196(4-6):10 06–17. [31] Prasad J, Diaz AR. Viscoelastic material design with negative stiffness components using topology optimization. Struct. Multidiscip. Optim. 2009;38(6):583–97. [32] Chen W, Liu S. Topology optimization of microstructures of viscoelastic damping materials for a prescribed shear modulus. Struct. Multidiscip. Optim. 2014;50(2):287–96.
[m5G;May 17, 2017;19:48] 11
[33] Wang Y, Luo Z, Zhang N, Qin Q. Topological shape optimization of multifunctional tissue engineering scaffolds with level set method. Struct. Multidiscip. Optim. 2016;54(2):333–47. [34] Huang X, Xie YM, Jia B, Li Q, Zhou SW. Evolutionary topology optimization of periodic composites for extremal magnetic permeability and electrical permittivity. Struct. Multidiscip. Optim. 2012;46(3):385–98. [35] Torquato S, Hyun S, Donev A. Multifunctional composites: optimizing microstructures for simultaneous transport of heat and electricity. Phys. Rev. Lett. 2003;89(26) 266601-1-266601-4. [36] Challis VJ, Roberts AP, Wilkins AH. Design of three dimensional isotropic microstructures for maximized stiffness and conductivity. Int. J. Solids Struct 2007;45(14-15):4130–46. [37] Huang X, Zhou SW, Xie YM, Li Q. Topology optimization of microstructures of cellular materials and composites for macrostructures. Comput. Mater. Sci. 2013;67:397–407. [38] Du YX, Chen LP, Tian QH, Wu ZJ. Topology synthesis of thermomechanical compliant mechanisms with geometrical nonlinearities using meshless method. Adv. Eng. Software 2009;40(5):315–22. [39] Neves MM, Rodrigues H, Guedes JM. Optimal design of periodic linear elastic microstructures. Comput. Struct. 20 0 0;76(1):421–9. [40] Andreassen E, Andreasen CS. How to determine composite material properties using numerical homogenization. Comput. Mater. Sci. 2014;83(2):488–95. [41] Bendsøe MP, Sigmund O. Topology Optimization: Theory, Methods, and Applications. Berlin Heidelberg, New York: Springer; 2003. [42] Andreassen E, Jensen JS. Topology optimization of periodic microstructures for enhanced dynamic properties of viscoelastic composite materials. Struct. Multidiscip. Optim. 2013;49(5):695–705. [43] Liu Y, Yin Y. Multi-objective topology optimization for microstructure of periodic composite material. Mater. Res. Innovat. 2011;15:s260–2. [44] Wang YQ, Luo Z, Zhang N, Kang Z. Topological shape optimization of microstructural metamaterials using a level set method. Comput. Mater. Sci. 2014;87:178–86. [45] Radman A, Huang X, Xie YM. Topological optimization for the design of microstructures of isotropic cellular materials. Eng. Optim. 2013;45(11):1331–48. [46] Cadman JE, Zhou S, Chen Y, Li Q. On design of multi-functional microstructural materials. J. Mater. Sci. 2013;48(48):51–66. [47] Huang X, Radman A, Xie YM. Topological design of microstructures of cellular materials for maximum bulk or shear modulus. Comput. Mater. Sci. 2011;50(6):1861–70. [48] Zhang W, Dai G, Wang F, Sun S, Bassir H. Using strain energy-based prediction of effective elastic properties in topology optimization of material microstructures. Acta Mech. Sin. 2007;23(1):77–89. [49] Zhang W, Wang F, Dai G, Sun S. Topology optimal design of material microstructures using strain energy-based method. Chin. J. Aeronaut. 2007;20(4):320–6. [50] Xia L, Breitkopf P. Design of materials using topology optimization and energy-based homogenization approach in Matlab. Struct. Multidiscip. Optim. 2015;52(6):1229–41. [51] Kaminakis NT, Drosopoulos GA, Stavroulakis GE. Design and verification of auxetic microstructures using topology optimization and homogenization. Archive Appl. Mech. 2015;85(9-10):1289–306. [52] Du YX, Chen D. Suppressing gray-scale elements in topology optimization of continua using modified optimality criterion methods. Comput. Model. Eng. Sci. 2012;86(1):53–70. [53] Yin L, Wei Y. Optimality criteria method for topology optimization under multiple constraints. Comput. Struct. 2001;79(20-21):1839–50. [54] Wadbro E, Hägg L. On quasi-arithmetic mean based filters and their fast evaluation for large-scale topology optimization. Struct. Multidisc. Optim. 2015;52(5):1–10. [55] Gibson LJ. Modelling the mechanical behavior of cellular materials. Mater. Sci. Eng. 1989;A110(89):1–36. [56] Li X, You M. Mechanical property analysis and numerical simulation of honeycomb sandwich structure’s core. Adv. Mater. Res. 2013;631/632:518–23. [57] ASTM C273-00 Standard Test Method for Shear Properties of Sandwich Core materials.
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
JID: ADES
[m5G;May 17, 2017;19:48]
Advances in Engineering Software 0 0 0 (2017) 1–11
Contents lists available at ScienceDirect
Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft
Topological design optimization of lattice structures to maximize shear stiffness Yixian Du a,b,∗, Hanzhao Li a,b, Zhen Luo c, Qihua Tian a,b a
College of Mechanical & Power Engineering, China Three Gorges University, Yichang 443002, China Hubei Key Laboratory of Hydroelectric Machinery Design and Maintenance, Yichang 443002, China c School of Electrical, Mechanical and Mechatronic Systems, The University of Technology, Sydney, NSW 2007, Australia b
a r t i c l e
i n f o
Article history: Received 9 December 2016 Revised 23 March 2017 Accepted 30 April 2017 Available online xxx Keywords: Lattice structure Topology optimization Energy-based homogenization method
a b s t r a c t To improve the poor shear performance of periodic lattice structure consisting of hexagonal unit cells, this study develops a new computational design method to apply topology optimization to search the best topological layout for lattice structures with enhanced shear stiffness. The design optimization problem of micro-cellular material is formulated based on the properties of macrostructure to maximize the shear modulus under a prescribed volume constraint using the energy-based homogenization method. The aim is to determine the optimal distribution of material phase within the periodic unit cell of lattice structure. The proposed energy-based homogenization procedure utilizes the sensitivity filter technique, especially, a modified optimal algorithm is proposed to evolve the microstructure of lattice materials with distinct topological boundaries. A high shear stiffness structure is obtained by solving the optimization model. Then, the mechanical equivalent properties are obtained and compared with those of the hexagonal honeycomb sandwich structure using a theoretical approach and the finite element method (FEM) according to the optimized structure. It demonstrates the effectiveness of the proposed method in this paper. Finally, the structure is manufactured, and then the properties are tested. Results show that the shear stiffness and bearing properties of the optimized lattice structure is better than that of the traditional honeycomb sandwich structure. In general, the proposed method can be effectively applied to the design of periodic lattice structures with high shear resistance and super bearing property. © 2017 Elsevier Ltd. All rights reserved.
1. Introduction The recent developments on large scale computation and advanced manufacturing techniques allow the construction of designs with fine and complex geometrical features [1], e.g. cellular and lattice structures, which have been increasingly investigated for a range of lightweight structural applications, including aerospace, automotive and naval industries due to their high stiffness-to-weight and strength-to-weight ratio, as well as their excellent energy absorption and thermal isolation characteristics [2–6]. As a representation, lattice materials are often assumed to be homogeneous due to it is a number of periodic units and can be efficiently assembled by following a certain pattern. Recently, the honeybee comb structure inspired by natural beeswax has long fascinated engineers and biologists for its outstanding mechanical properties and other characteristics, as a result of their special microstructures such as hexagonal, square and triangular, which are
∗
Corresponding author. E-mail address: [email protected] (Y. Du).
widely used in various engineering applications [7]. However, it has not proven until recently that this is the best structure with enhanced shear resistance that can be built with the beeswax to store the largest amount of honey [8]. Shear response is of great significance as bending lattice structure gives rise to transverse shear loading and may collapse the core of the structure in shear. For instance, the work [2] studied the effective in-plane stiffness of hexagonal honeycomb cores in macro-scale according to the bending model. After that, there is a growing interest in exploiting other materials as the filling material in such structures to simultaneously enhance the shear response, load-bearing and energy absorption capabilities of the lightweight structures, such as periodic lattice structure and honeycomb sandwich cores [9–11]. Many researchers [12–14] studied the mechanical properties of other honeycomb cores such as hollow pyramidal lattices, square honeycombs, and corrugated sandwich cores. Kelsey et al. [15] gave the upper and lower bounds of the shear modulus of honeycomb cores by using the unit displacement and unit load methods in conjunction with a simplification assumption for the strain and stress in the core. Pan et al. [16] and Han et al. [17] analyzed the longitudinal shear strength of honeycomb
http://dx.doi.org/10.1016/j.advengsoft.2017.04.011 0965-9978/© 2017 Elsevier Ltd. All rights reserved.
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
JID: ADES 2
ARTICLE IN PRESS
[m5G;May 17, 2017;19:48]
Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
cores, and investigated the longitudinal shear deformation behavior and failure mechanism of aluminum alloy honeycomb cores at room temperatures. It is demonstrated that the shear or compression strength and specific energy absorption of lattice structures or sandwich cores can be increased dramatically by using appropriate filling foams. Also there have been some papers which studied the shear responses of various periodic lattice structures were mainly focused on the calculation and comparison of their shear modulus and experimental validations [18,19]. However, most of the lattice structures used in engineering are obtained by designers’ experiences or inspired by existing materials or structures in nature. An effective method to guide the design of microstructural topologies of lattice materials is still in demand. Inspired by the above natural-occurred cellular materials, the design of periodic lattice composite to achieve multifunctional properties (stiffness, strength, negative Poisson’s ratio, etc.) based on the combination of structural optimization techniques and numerical homogenization has attracted considerable attention. For instance, Sigmund [20–22] employed a modified optimality criteria method to design material microstructure to achieve the prescribed properties. It can be found that the effective material properties at macro-scale are mainly determined by the topology of periodic microstructures, rather than the proportion and physical properties of their constituents, namely, the intrinsic material composition [23]. The homogenization method [24–25] has been recognized as a rigorous method for characterizing the macro mechanical behavior of materials consisting of periodic microstructures, under the assumption that the linearly elastic response of the periodic material can be determined by test strains over one cell [26]. The topological design of such materials assumes that the material is made of periodic cells which theoretically are infinitesimal within the structure, and the macro effective properties of the heterogeneous material are homogenized according to the microstructure that is the smallest repetitive unit of the material. The inverse design in the material is a typical topology optimization problem, which seeks an optimal microstructure of the cell with prescribed or extreme effective properties [27–36], and so on. A systematic mean of microstructural design is formulated as an optimization problem for parameters that represent the property and topology of the microstructure [37], so as to improve structural shear effect. Over the past, topology optimization has been expanding as a powerful computational design tool for a range of structural and material problems both in academic research and industrial applications. Essentially, topology optimization is a numerical iterative process that distributes a given amount of material inside a reference design domain, so as to seek the best material layout until the expected performance is optimized subject to constraints [38]. Since the work of Sigmund [20–22], various structural topology optimization techniques have been developed for computationally design of microstructures of materials [33,39], e.g. the numerical homogenization method [40], solid isotropic material with penalization (SIMP) [41–43], level set method (LSM) [36], parametric level set method (PLSM) [44], and bidirectional evolutionary structural optimization (BESO) [45]. Amongst them, the design multifunctional composite materials are an active field of research [46]. Huang et al. [47] recently applied the BESO approach to the design of 2D and 3D materials with extreme bulk and shear modulus, and isotropic constraint is imposed by Radman [44]. Zhang et al. [48,49] and Xia et al. [50] proposed the strain energy-based method to predict the effective elastic properties, which compared with the numerical homogenization method and shows the advantages of higher computing efficiency and simplicity in the numerical implementation. Based on the above work, we can find that there is still much room to improve the shear resistance of periodic lattice structure
consisting of hexagonal unit cells in engineering applications. Lattice structures composed of unit cells are one of the candidates to enhance structural/material performances. However, the topological design problem of periodic lattice structure in recent works are mostly generalized by using the mathematical homogenization method, which has the disadvantages, such as the numerical implementation of the homogenization method and design sensitivity analysis are complicated and additional programming is needed. In this paper, the topology optimization technique is utilized to search the layout of periodic lattice structure with enhanced shear stiffness under the constraint of a given amount of material by using energy-based homogenization method. The optimization objective is to find the optimal distribution of two base materials within the unit cell of the periodic lattice material, so that the resulting structure has the maximum shear modulus. The effective properties of the microstructure are totally predicated by the periodic unit cell (PUC) which is discretized into finite elements under periodic boundary condition. 2. Energy-based homogenization method 2.1. Formulation and sensitivity analysis with energy-based homogenization In general, the homogenization method applied to linear elastic problems establishes macroscopic properties that often are uniformly described by corresponding PUCs, which are periodically repeated in one or more directions. The properties of a heterogeneous medium rely on the analysis of its PUCs. Thus, the homogenization theory can be used to calculate the effective elastic constants of the macroscopic composite by the Hill-Mandel condition or energy averaging theorem [51]. To obtain the macroscopic equivalent constitutive properties by the homogenization process, three widely used types of loading can be applied: (a) prescribed linear displacements, (b) prescribed tractions and (c) periodic boundary conditions. The size of the periodic unit cell is assumed to be much smaller than the size of the bulk material, so that the average quantities of both the microscopic strain and stress can be defined [51]. Consider a single cell Y, the macro-scale displacement field μ is expanded, depending on the little aspect ratio ε between the macro and micro scales by using the asymptotic homogenization. When only the first-order terms of the asymptotic expansion are considered, the effective elasticity tensor of the macro-material can be found in terms of the material distribution in the domain of PUC by averaging the integral over the base cell [50] Y as:
CiHjkl =
1 |Y |
Y
ε 0pq(i j ) − ε ∗pq(i j ) Cpqrs ε 0pq(kl ) − ε ∗pq(kl ) dY
(1)
where |Y| represents the area (or volume) of the unit cell, and ε 0pq(kl ) defines the prescribed macroscopic strain fields, which are defined as three linearly independent unit strains: the horizontal unit strain ε 0pq(11) = [1 0 0], the vertical unit strain ε 0pq(22) = [0 1 0] and the shear unit strain ε 0pq(12) = [0
0 1]. The local
strain field ε ∗pq(kl ) induced by the test strains are defined as follow:
1 kl kl ε ∗pq(kl ) = ε ∗pq χ kl = χ p,q + χq,p 2
(2)
where χ kl denotes a Y-periodic admissible displacement field associated with the load case kl, which can be obtained from the following equilibrium equation [50]:
∗ kl ∂ vi Ci jpq ε pq( ) dY = ∂yj Y
Y
0 kl ∂ vi Ci jpq ε pq( ) dY ∂yj
(3)
where v is the Y-periodic admissible displacement field.
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
ARTICLE IN PRESS
JID: ADES
[m5G;May 17, 2017;19:48]
Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
The stress and strain tensors of the homogeneous medium are equivalent to the average stress and strain of the microstructure, thus, the periodic microstructure of composite materials can be replaced by an equivalent homogeneous medium with the same volume at the macroscopic level [49]. The average stress and strain of the homogeneous medium follow the Hooke’s law, which can be rewritten in matrix form, as follows
⎡σ ⎤
⎡
11
H C1111
H ⎣σ22 ⎦ = ⎢ ⎣C1122 σ12 0
H C1122 H C2222
0
⎤⎡ ⎤ 0 ε11 ⎥⎣ ⎦ 0 ⎦ ε22 H ε12 C1212
(4)
where σ11 , σ22 denotes the horizontal and vertical averaged stress. ε11 , ε22 denotes the averaged strain in corresponding two directions. σ12 , ε12 denotes the averaged shear stress and strain in two directions. In fact, the strain energy-based method is equivalent to the numerical homogenization method in predicting the effective properties of materials [48]. Considering the shear-strain state, which are ε11 = ε22 = 0, ε12 = 1, the unit strain energy of the microstructure is expressed as Et = σ12 ε12 /2 = σ12 /2. The strain energy of microstructure under periodic boundary condition is expressed by
Et =
1 1 H F U = U T KU = C pqrs ε Apq(i j ) εrsA(kl ) 2 2
(5)
where F, U and K are the load vector, the displacement vector of element nodes and the stiffness matrix, respectively. The energy-based approach employs average stress and strain theorem, in which the unit test strains are imposed directly on the boundaries of the base cell instead of the asymptotic approach, including ε Apq(kl ) which corresponds to the superimposed strain fields
(ε 0pq(kl ) − ε ∗pq(kl ) ) in Eq. (1).
In finite element analysis, the base cell is discretized into N finite elements and the shear modulus is obtained in an equivalent form in terms of element mutual energies H C1212 =
σ12 = 2Et = ε12
Y
A (i j ) A i j A kl C pqrs ε pq( ) εrs( ) dY = ue N
T
A kl ke ue ( )
e=1
(6) uAe (kl )
where are the element displacement solutions corresponding to the mentioned unit test strain fields ɛ0( kl ) , and ke is the element stiffness matrix. The effective elastic constant of the composite highly depends on the spatial distribution of material phases, therefore, how to optimally distribute material phase within the PUC would be of high importance. Thus, microstructures can be modeled as 2D continuum structures. An artificial design variable, x, is introduced to represent material occupation in an element by assuming that x = 1 if an element is made of solid material and x = 0 versus. The element property is interpolated using the SIMP model, where the density-stiffness interpolation scheme [41] can be stated as
Ee (x ) = Emin + (ρe (x ) ) p (E0 − Emin )
(7)
where ρ e (x) designates the relative density in the random point x of the design domain, which takes values between 0 and 1. E0 is the Young’s modulus of solid material and Emin is a very small stiffness assigned to void regions for approximation of void material (namely the Ersatz material or compliant material) in order to prevent the stiffness matrix from becoming singular. p is a penalization factor (typically p = 3–4) introduced to ensure the density distribution closer towards the black-and-white solutions. Here, topology optimization of the material microstructure aims at finding the optimal microstructure of the PUC that exhibits the maximum shear stiffness in specified directions under prescribed
3
periodic boundary conditions. The mathematical formulation of the optimization problem can be stated as
⎧ H maxρ : f CiHjkl (ρ ) = C1212 ; ⎪ ⎪ ⎪ ⎪ A(kl ) (kl ) ⎪ ⎪ ⎨s.t. : KU = F , k, l = 1, · · · , d N
⎪ ve ρe /|Y | ≤ ϑ ⎪ ⎪ ⎪ ⎪ e=1 ⎪ ⎩ 0 < ρe ≤ 1 , e = 1 , · · · , N
(8)
where the objective f (CiHjkl (ρ ) ) is a function of the homogenized stiffness tensor, namely the maximization of shear modulus mentioned in Eq. (6). K is the global stiffness matrix, UA ( kl ) and F( kl ) are the global displacement vector and the external force vector of the test case (kl), respectively. ve denotes the element volume, ϑ is the limited upper bound of volume fraction, and d is the spatial dimension. It is critical to acquire gradient information to guide the algorithm to search for an optimum point which is most efficient during the iteration process, which is an important step in topology optimization. In this regard, topological sensitivity is frequently defined as the derivative with respect to design variable, which is the relative density of element in the density-based finite element framework. So, the sensitivity of the strain energy with respect to the element density variable is obtained as
∂ Et ∂ 1 T 1 T ∂K ∂U T = U KU = U U +2 KU ∂ ρe ∂ ρe 2 2 ∂ ρe ∂ ρe
where
∂U T ∂ UT KU = ∂ ρe ∂ ρe
∂ U T ∂ ρe
F
F
= 0
∂ U T ∂ ρe
(9a)
F 0
=0
(9b)
when the duality condition of displacement and force on the boundary are considered, here U and F are the nodal displacement vector and force vector on the boundary. U and F are the nodal displacement vector and force vector inside. The periodic boundary condition mentioned above impose three unit test strains directly on the boundaries of the base cell, the specific nodal displacement vector have no relationship with density variable, so the derivative of nodal displacement on the boundary with respect to density equals to zero. What’s more, the force vector inside the unit cell indeed equals to zero to meet the continuity and periodicity requirements for the displacement as well as stress in the FEA. With the help of the above mentioned material interpolation, the first-order derivative of the strain energy with respect to the element density is given by:
∂ Et 1 ∂K p n = UT U= E ∂ ρe 2 ∂ ρe ρe t
(9c)
It is easily to get the sensitivity of the homogenized elasticity H tensor C1212 with respect to ρ e , which can be derived from the adjoint method [41] as: H A 12 T A 12 ∂ C1212 1 ∂ Et p = pρ p−1 (E0 − Emin ) ue ( ) ke ue ( ) = 2 = 2 Etn ∂ ρe ∂ ρe ρe |Y | e
(10) where Etn is the element strain energy, which effects the sensitivity of the effective elastic tensor. In order to avoid the formation of checkerboard pattern and the mesh-dependency in topology optimization, the filtering scheme [41] is used during the process of solving the optimization problem, which modified the design sensitivities during iterations, as
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
ARTICLE IN PRESS
JID: ADES 4
Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
1 0.8 The update density
As indicated in Fig. 1, the polarization operator can reasonably make the intermediate density compelling to the two ends of 0– 1, which ensure the optimal solution approaching to 0/1 discrete properties. The curves are become steeper and steeper with the polarization factor’s increasing. The proposed algorithm is just like a smoothed Heaviside function that can polarize the majority of elements to 0 or 1.0 more rapidly at the same probability. The number of elements with the intermediate density in the final topology will be reduced, so as to get a topological design with more clear edges. The standard optimality criteria (OC) method has been widely recognized as an effective optimization algorithm for its many advantages while used, such as fast convergence, simplified programming and with no association with the structural re-analysis and the number of the variables [41]. Combine the polarization operator with the optimality criteria equation, a new iterative updating scheme is expressed as follows:
s=9 s=11 s=13 s=16
0.9 0.7 0.6 0.5 0.4 0.3 0.2 0.1 0 0
0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 The initial density
1
Fig. 1. Optimization with different polarization factors.
k xe = xne Be( ) n
η
follows
∂ Cˆ 1 ∂C = Hi ρi ∂ ρe max(γ , ρe ) i∈Ne Hi ∂ ρi
(11)
where Ne is the set of elements i for which the center-to-center distance to element e is smaller than the filter radius rmin and Hi is a weight factor defined as: Hi = max(0, rmin - (ei)), here the term γ = 10−3 is a small positive number introduced in order to avoid division by zero. The filter radius of rmin is defined to identify the neighboring elements that affect the sensitivity of the current element. 2.2. Modified optimization algorithm Due to the relaxation of original 0/1 discrete problem to a continuous 0–1 optimization problem, the numerical problems of gray regions, checkerboards and mesh dependence mostly exist in the topology optimization. The appearance of such geometrical features in either solid or the void regions may lead to lack of manufacturability in engineering applications. Thus, to more effectively compress gray areas in the final topological design and achieve the topological design with clear boundary, several different methods have been developed [52,53], such as new density interpolation function, extra constraints, postprocessing techniques and filtering schemes [54]. In this paper, we proposed an optimization algorithm with the approximation of 0/1 discrete properties using a heuristic polarization operator. The basic idea of this method is that it can weaken the low-relativedensity elements while enhance the high-relative-density elements by the polarization operation during the iteration process. The polarization further makes those intermediate density elements towards 0 or 1 discrete design, so that the final solution has the discrete properties. In this way, we can get better design in topology optimization with lower gray elements and computational efficiency. The specific mathematical expression of polarization operator is defined as follows:
x,e
[m5G;May 17, 2017;19:48]
=
exp − 2s + s × xe
1 + exp − 2s + s × xe
(12)
where xe is the relative element density before polarization by polarization operator, xe ’ denotes the relative element density after polarization, s is the polarization factor that controls the speed of the polarization. Fig. 1 shows the different shapes when the polarization factor in Eq. (12) have different values.
=
k exp − 2s + s × xne Be( )
η
k 1 + exp − 2s + s × xne Be( )
η
⎧ n n n ⎪ ⎨ max (xmin , xe − m) xe < max (xmin , xe − m) xne +1 = xne max (xmin , xne − m ) < xne < min (1, xne + m) ⎪ ⎩ min (1, xne + m) xne > min (1, xne + m)
(13)
(14)
η
where [xne (Be(k ) ) ] is the not-polarized solution of the polarization operator before the process of solution, namely the initial solution; η
[xne (Be(k ) ) ] is the polarized solution of the polarization operator before the process of solution, namely the update solution; s is the polarization factor, which is used to control the punishment speed of the elements that belong to the middle part, here we choose s = 7 in this paper; xne is the updated solution by using the optimality criteria method. 2.3. Numerical example Assume that the Young’s modulus and Poisson’s ratio of the artificial solid phase are selected as E = 1, and μ = 0.3, the given volume fraction of solid material is set to be 0.4 and 0.6, the PUC is discretized with a mesh of 50 × 50 and 100 × 100 four-node quadrangle elements, and the penalty factor p and the filter radius rmin are endowed with different values, the corresponding microstructural topologies, material distribution and effective elasticity matrices are given in Table 1. From the optimal solutions and the corresponding effective matrices for different parameters cases, it can be seen that these different volume fractions, the mesh, the filter radius and the penalty factors all have influence on the final topological shapes of the unit cell with the same objective function. Under the same volume fraction, the equivalent elasticity matrix of the optimal microstructure is almost the same. This is reasonable due to the fact that a number of different microstructures can possess the same physical property. At the same time, undesirable multiple local minima that occurred in deleting the intermediate density elements during the searching process are avoided in the optimization procedure. Moreover, with the help of sensitivity technique, the number of elements with the intermediate density in the final topology will be reduced, so as to get a topological design with clear boundaries. It can be seen from Table 1 that the final optimal configuration is symmetrical, and the material layout is commonly in line with diagonal directions which is 45°. Thus, the potential of the material is motivated to resist shear loads so that the optimal microstructures exhibit high shear stiffness. Fig. 2 shows the iteration histories of the shear modulus, the unit cell is discretized by 100 × 100 four-node quadrilateral finite
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
JID: ADES
ARTICLE IN PRESS
[m5G;May 17, 2017;19:48]
Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
5
Table 1 Microstructural topologies, material distribution and effective elasticity matrix of cellular materials with maximum shear modulus for various volume fractions, mesh sizes, filtering radius and penalty factors.
elements, volume fraction is 0.5, filter radius is 5 and penalty factor is 3, respectively. It can be seen from Fig. 2 that when the update scheme given in Eqs. (13) and (14) is used, the shear modulus increases rapidly at the first 27 iterations, and then the curve tends to stabilize at iteration 140, the topological details change slightly. After 230 iterations, the objective function and the microstructural topologies stably converge to the final solutions. The distribution of design variables within the unit cells (Fig. 3) can be interpreted as octagonal honeycomb cell whose relative shear modulus equals to 0.106.
According to the microstructure of unit cells shown in Fig. 3, the 3 × 3 base cells are constructed periodically which can be called periodic lattice structure as shown in Fig. 4. When the design problem formulation of the periodic lattice structure has been done, the obtained optimal structure can be regarded as one octagonal honeycomb cell with a nested diamond-style structure inside, and this structure are almost the same as the solutions by Huang [47]. Also the lattice structure can be interpreted as the modified version of hexagonal honeycomb which are assembled by four cells in a certain pattern. In the next section, to show the good perfor-
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
ARTICLE IN PRESS
JID: ADES 6
[m5G;May 17, 2017;19:48]
Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
Fig. 2. Iteration histories of shear modulus and microstructure using energy-base homogenization.
effective in-plane elastic stiffness of honeycomb cores. While this analysis model ignore the effect of transverse shear stress and the different thickness of the unit cell wall. Therefore, the analytical analysis model is not used for the designed lattice structure directly. The effective elastic constants of designed periodic lattice structure’s core are re-derived by considering the transverse shear deformation and the flexibility of the cell wall. To verify the shear stiffness of the lattice structure, the equivalent elastic constants are deduced and compared with those of the hexagonal honeycomb sandwich structure. 3.1. Mechanical equivalent model of periodic lattice structure
Fig. 3. Distribution of design variables.
The simplified mechanical model of unit cell is built in Fig. 5 to analyze its equivalent mechanical performance. In Fig. 5, t denotes the thickness of the cell unit (mm), b is the height of the lattice structure (mm), l denotes the length of the longitudinal wall (mm), h represents the length of the inclined wall (mm), Es and Gs are Young’s modulus and shear modulus of base material of periodic lattice structure (MPa), θ denotes the angle between the hexagon short edge and horizontal side (°), respectively. Using force balance principle, the cantilever beam model and the material mechanics beam bending theory, the equivalent mechanical parameters of the lattice structure are presented, as is shown in Fig. 5. 3.2. Equivalent shear modulus analysis
Fig. 4. 3 × 3 unit cells.
mance of the lattice structure with shear stiffness, the equivalent mechanical properties are studied in macro-scale based on the theoretical model. 3. Theory analysis of periodic lattice structure According to different geometries of the microstructure, the topological configuration of two-dimensional periodic lattice structure may vary widely, but the general analytical formulation to characterize the equivalent mechanical properties is similar. Gibson developed the classical unit cell theory [2] for predicting the
Considering the symmetry of the most basic hexagon and trapezium composition unit of the lattice structure, the loads, displacements, strains and stresses on cell walls subjected to transverse shear displacement is shown in Fig. 5. The mechanical behavior of the inclined wall can be investigated based on the cantilever beam model due to the fact that the bending deformation of the inclined wall is anti-symmetric along the height direction just like the cantilever beam. The bending displacement of the inclined wall is supposed to be δ 1 , the equivalent load P1 can be written as follows:
P1 =
(δ1 /2 ) × 3Es I1 (b/2 )3
(15)
where I1 is the moment of inertia of the inclined wall, which can be expressed as I1 = ht13 /12, t1 is the thickness of the inclined wall, which can be expressed as t1 = 2t.
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
ARTICLE IN PRESS
JID: ADES
[m5G;May 17, 2017;19:48]
Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
7
Fig. 5. Simplified equivalent model of periodic lattice structure.
Suppose the deformation under transverse shear displacement for the longitudinal wall consists of two parts: the bending displacement δ 2 and shear displacement δ s . The corresponding equivalent loads are written as
P2 =
(δ2 /2 ) × 3Es I2 (b/2 )3
Ps = Gs δs lt/b
(16) (17)
where I2 is the moment of inertia of the longitudinal wall, which can be expressed as I2 = lt3 /12. The shear displacement is supposed to be δ y in the transverse direction like the x-direction shown in Fig. 5, the bending displacement of inclined wall, the shear displacement and bending displacement of longitudinal wall can be expressed as: δ 1 = δ y , δ s = δ y × cos θ , and δ 2 = δ y × sin θ . In order to obtain the precise expression of the transverse shear modulus of the lattice structure, the different thickness, the bending contribution of the inclined wall and the longitudinal wall should be considered. The mechanical equilibrium condition of the representative unit requires eq P1 + 2P2 sin θ + 2Ps cos θ = 2τxy (h + l sin θ )l cos θ
(18)
where τxy is the equivalent shear-stress that imposed on the equivalent area of one representative unit cell which can be exeq eq eq eq pressed as τxy = Gxy γxy , γxy is the equivalent shear strain that imposed on the equivalent area of one representative unit cell which eq can be expressed as γxy = δy /b. Thus, substituting Eqs. (15), (16) and (17) into Eq. (18) yields the following equation that computes the equivalent shear modulus of the periodic lattice structure: eq
Geq xy =
Gs t cos θ 4ht 3 + lt 3 sin2 θ Es + h + l sin θ (h + l sin θ )l cos θ b2
(20)
The mass of the equivalent entity plate may be expressed as
(21)
According to the conservation of mass before and after the equivalence (mce = m1 ), we have
ρc =
4t (l + h )ρs l cos θ (l sin θ + h )
⎧ Gs t cos θ 4ht 3 + lt 3 sin2 θ Es ⎪ + ⎨Geq xy = h + l sin θ h ( + l sin θ )l cos θ b2 ⎪ 4 ρs t ( β + 1 ) ⎩ ρc = l (β + 2 sin θ )(β + 2 cos θ )
(23)
eq
where Gxy is the equivalent shear modulus of lattice structure in xy-plate, ρ c and ρ s are the equivalent densities of the lattice structure and the base material of the periodic lattice structure. 3.4. Comparison of two structure’s shear moduli
Based on the conservation of mass, the mass of the representative unit cell wall can be calculated at first, and then the structure is equivalent to an entity plate with the same size. The equivalent mass of the representative unit cell of the periodic lattice structure may be expressed as
mce = ρcVce = ρs · b(h + 2l sin θ )(h + 2l cos θ )
As a result, the equivalent mechanical parameters of the periodic lattice structure are presented as follows:
(19)
3.3. Equivalent density analysis
m1 = ρsV1 = ρs × 4bt(h + l )
Fig. 6. Simplified equivalent model of honeycomb core by Gibson [55].
(22)
A case of bearing structure that is made of 2024 aluminum alloy material is used. The thickness of the hexagonal honeycomb sandwich structure and the periodic lattice structure is given the same value. Fig. 6 shows the Gibson’s analysis model under shear deflection of the honeycomb core. The parameters are supposed as: θ = 45˚, h = l, and β = h/l = 1, thus, Eq. (23) can be simplified as follows:
⎧ 9Es t 3 + Gs t b2 eq ⎪ ⎪ ⎨Gxy = √ 2 2 + 1 lb
8 ρs t ⎪ ⎪ √ ⎩ ρc =
(24)
l 3+2 2
According to Gibson’s equation [55], the thickness of the representative unit cell of hexagonal honeycomb sandwich is assumed to t. The corresponding equivalent analysis model and parameters
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
ARTICLE IN PRESS
JID: ADES 8
[m5G;May 17, 2017;19:48]
Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
P
Fig. 7. Simplified load condition and constraints.
Fig. 8. FE models of the two lattice structures.
are shown in Fig. 6. When θ = 30˚, h = l, β = 1 is choosing, the corresponding equation can be simplified as
⎧ Es t 3 ⎪ ⎨Gacxy = √ 3 3l
Table 2 Size parameters of periodic lattice structure.
(25)
ρs t ⎪ ⎩ ρac = 8√
Cell size of core (mm)
t 0.05
l 5
h 10
Base material (aluminum)
Es (GPa) 70
vs 0.26
2780 kg/m3
ρs
3 3l
In order to intuitively compare the shear modulus of the two structures theoretically under the condition that the mass of the entity of structures are the same, the equivalent density is also assumed to√be the same value, namely, ρ c =ρ ac , thus, ρc /ρac = √ 3 3/(3 + 2 2). The ratio between the shear modulus of two structures is defined as nG , which can be calculated by
Geq xy nG = = Gacxy
√ 3 9t 2 + GEss b2 l 2 √ 2 + 1 b2 t 2
(26)
In the above equation, the Young’s modulus and the shear modulus of the same isotropic material have a relation of Gs = Es / 2(1 + vs ). Furthermore, l2 / t2 >> 1 can be easily obtained according to the application of engineering. As a result, Eq. (26) can expressed as follows:
nG =
√ √ 3+2 2 9 3t 2 l2 √ + 1.86 2 ≈ 1.86 × ≈ 2.09 > 1 √ t 3 3 2 + 1 b2
(27) eq Gxy
From above equation, it is evident that > Gacxy can be easily obtained. In other words, the shear modulus of periodic lattice structure is highly increased when considering the same mass and the same application condition. It also may be concluded that, the proposed method can be effectively applied to enhance the shear stiffness of periodic lattice structure theoretically. 4. Finite element analysis and experimental investigation 4.1. Finite element analysis The FE simulation analysis for the mechanical performance of the periodic lattice structure and the honeycomb sandwich structure is given as follows [56]. The two structures are made of 2024
aluminum alloy materials, with yield strength 758 MPa and material density 2780 kg/m3 , and other data are shown in Table 2. To compare the shear performance and bearing ability of this two structures, the FE models of the two structures were given in Fig. 8 by using the same type of elements (SHELL63 plate-like elements in the three dimensional space) under the same geometry and load conditions, as shown in Fig. 7. Both the face-sheets were taken as rigid with no deformation. As shown in Fig. 7, the periodic lattice structure assembled as a plate, the two sides of the plate were constrained, which means all degrees of freedom were clamped. At the top surface, each node was subjected to a constant shearing force (linear elastic analysis) along the thickness direction. The initial value of applied load is 0.005 MPa, and the incremental value is set to be 0.005 MPa as well, 10 sets of numerical experiments were done for the purpose of a better comparison under different sets of comparison and analysis. According to mechanics, the stress-strain relation of periodic lattice structure can be used to calculate the shear modulus as shown in Eq. (19). Table 3 shows the results of FEM simulation and theoretical analysis, for the two structures with the same geometry parameters and simulation conditions. The theoretical results and numerical results are expressed as Ti and Ni respectively under the load case i and the relative error can be calculated by:
error = i
Ni − Ti Ti
(28)
It can be seen from Table 3 that the maximum relative error between the theoretical and FEM numerical results does not exceed 5%. Results from the proposed analytical model well align with the FE predictions for the periodic lattice configurations in this study. Therefore, it is clear that the proposed mechanical equivalent model for calculating the equivalent elastic modulus are rea-
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
ARTICLE IN PRESS
JID: ADES
[m5G;May 17, 2017;19:48]
Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
9
Table 3 The shear modulus of periodic lattice structures under various loads. Loads(MPa)
0.005
0.010
0.015
0.020
0.025
0.030
0.035
0.040
0.045
0.050
Theoretical results(MPa) Numerical results(MPa) The relative error
1.997 2.030 1.7%
1.997 2.002 0.3%
1.997 1.992 −0.3%
1.997 1.988 −0.5%
1.997 1.986 −0.6%
1.997 1.984 −0.7%
1.997 1.980 −0.9%
1.997 1.979 −0.9%
1.997 1.977 −1.0%
1.997 1.976 −1.0%
3.5
Nominal strain / mm
3
Periodic lattice structure Hexagonal honeycomb structure
2.5 2 1.5 1 0.5 0 0.005 0.01 0.015 0.02 0.025 0.03 0.035
0.04 0.045 0.05 0.055
Load / MPa Fig. 9. The strains of two structures under various loads.
sonable and effective. In order to obtain a more precise solution, the transverse shear deformation and the uniform thickness of the cell wall should be considered. The strain for the honeycomb sandwich and the periodic lattice structures under various loading conditions are listed in Fig. 9. It can be found from Fig. 9 that when the loading conditions are the same the strain of the honeycomb structure is larger than that of the periodic lattice structure. In other words, the shear resistance and bearing property are better than that of the honeycomb structure. The difference of strains between the honeycomb structure and the periodic lattice structure increase along with the increase of the applied loads. It indicates that the increase of the shear performance of the lattice structure is more obvious when the loads large enough. In a word, the designed periodic lattice structure shows a better shear resistance and bearing performance, compared to the widely used honeycomb sandwich structure in engineering. 4.2. Experimental investigation According to the results by using topology optimization technique in Section 2, the periodic lattice structure given as below has been greatly simplified, according to the manufacturing process. Our specimens have been manufactured in polyamide using the selective laser sintering (SLS) and Wire Saw Cutting. The specimens manufactured by two means are given in Fig. 10 respectively. The sample manufactured by SLS is fragile and it is difficulty to be used for practical experimental test. When the manufacturability considered, the experimental specimen was manufactured in aluminum following two steps: 1) prepare enough bars manufactured by using wire cutting technique, 2) assembly or paste the bars by using bonding technology. The thickness of the cell unit is t = 0.5 mm, The length L, width B and height b are 112 mm × 67 mm × 10.6 mm, separately. The shear tests have been done at room temperature. The specimens were tested under a constant cross-head movement of 0.1 mm/min by using the shearing equipment according to ASTM
Fig. 10. (a) Specimen manufactured by SLS. (b) Specimen manufactured by wire cutting technique.
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
JID: ADES 10
ARTICLE IN PRESS
[m5G;May 17, 2017;19:48]
Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
Fig. 11. (a) Schematic of test specimen. (b) Force-displacement of periodic lattice structure.
C273-00 [57]. The typical force-displacement curves obtained during the shear deformation of the periodic lattice structure were analyzed based on the observation of microstructural change during the deformation. Fig. 11(a) shows the shear test equipment and Fig. 11(b) shows a typical force-displacement curve obtained from shear test and a schematic microstructural change during the shear deformation of specimen. From Fig. 11, it can be seen that the shear force increased almost linearly with the increase of displacement at the beginning stage due to the elastic buckling of cell walls, not the elastic axial shortening of cell walls. The thin cell walls of periodic lattice structure are restrained with the neighboring cell walls. Therefore, when the shear deformation is performed to the lattice structure, it is difficult to be happened by the elastic axial shortening of cell walls. When the force-displacement curve passes a maximum stress, a maximum shear strength is shown in Fig. 11(b), then the shear stress rapidly dropped after the first stage due to the debonding fracture at cell interfaces. A continuous fracture of phenolic resin layer on cell walls is observed. During proceeding the plateau deformation from shear stress-strain curve, the fracture of specimen by the core interfaces debonding occurs finally with a decrease of shear stress on the periodic lattice structure. In the procedure, the deformation behavior of the cell walls during shear test was analyzed as shown in Fig. 11, which shows the basic cell carrying out the shear stress in periodic lattice structure. The arrows around cell walls indicate the shear force performed on the cell wall during shearing in transverse direction. The slope P / δ of the force-displacement curve during elastic deformation shown in Fig. 11(b) is regarded as the total shear modulus of the structure ( P, δ is selected to be 23.0414 and 0.016 as shown in Fig. 11), which are substituted into the shear performance test formula [57] to calculate the shear modulus as follows
Gc =
23.0414 b · P 10.6 = · = 2.03MPa L · B · δ 112 × 67 0.016
structures with the transverse shearing test, which demonstrates the excellence shearing performance of the periodic lattice structure. 5. Discussion and conclusions This study has developed a topology optimization method to enhance the transverse shear performance of the hexagonal lattice structure. The optimization problems have been formulated to generate the periodic microstructures with maximum shear modulus under a prescribed volume and constraint with the strain energybased homogenization method. An equivalent analytical solution for the transverse shear stiffness has been proposed according to the simplified optimal results by considering the transverse shear deformation and the uniform thickness of the cell wall. The shear stiffness of the designed lattice structure has been studied by using the theoretical approach and the finite element method, which demonstrates the effectiveness of the method used in this paper. Finally, the structure is prototyped, and then the shear properties are tested and compared with those of the widely used honeycomb sandwich structures in engineering. Results show that the proposed method can be effectively applied to the design of periodic lattice structures with higher shear resistance and better load bearing abilities. Acknowledgments The work is supported by the National Natural Science Foundation of China (51105229, 51475265), Science and Technology Support Program of Hubei Province of China (2015BHE026), the Research Project of Hubei Provincial Department of Education (D20161205).
(29)
The relative error between the experimental results and the theoretical results of shear modulus shown in Table 3 is 2%, which is less than 5%. That is to say, results from the experimental results are well in line with the finite element prediction and theoretical results as given in Table 3 for the periodic lattice structure. This demonstrates the correctness of the equivalent analysis and FEA model. What’s more, the measured shear modulus of the lattice structure in the transverse direction is compared with those of honeycomb structures in reference [18] under the same experimental conditions. As a result, the measured shear modulus of the lattice structure is 1.49 times higher than those of honeycomb
References [1] Zheng X, Lee H, Weisgraber TH, Shusteff M, DeOtte J, Duoss EB, et al. Ultralight, ultrastiff mechanical metamaterials. Science 2014;344(6190):1373–7. [2] Li H, Luo Z, Zhang N, Gao L, Brown T. Integrated design of cellular composites using a level-set topology optimization method. Comput. Meth. Appl. Mech. Eng. 2016;309:453–75. [3] Fleck NA, Ashby MF. Micro-architectured materials: past, present and future. Proc. R. Soc. A 2010;466(2121):2495–516. [4] Wang Y, Luo Z, Zhang N, Wu T. Topological design for mechanical metamaterials using a multiphase level set method. Struct. Multidiscip. Optim. 2016;54(4):937–52. [5] Mines RAW, Tsopanos S, Shen Y, Hasan R, McKown ST. Drop weight impact behaviour of sandwich panels with metallic micro lattice cores. Int. J. Impact Eng. 2013;60(10):120–32.
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011
JID: ADES
ARTICLE IN PRESS Y. Du et al. / Advances in Engineering Software 000 (2017) 1–11
[6] Montemurro M, Catapano A, Doroszewski D. A multi-scale approach for the simultaneous shape and material optimisation of sandwich panels with cellular core. Composites Part B 2016;91:458–72. [7] Zhang K, Duan H, Karihaloo BL, Wang J. Hierarchical, multilayered cell walls reinforced by recycled silk cocoons enhance the structural integrity of honeybee combs. Proc. Natl. Acad. Sci. U. S. A. 2010;107(21):9502–6. [8] Mackenzie D. Proving the Perfection of the Honeycomb. Science 1999;285(5432):1338–9. [9] Han B, Yu B, Xu Y, Chen CQ, Zhang QC, Lu TJ. Foam filling radically enhances transverse shear response of corrugated sandwich plates. Mater. Des. 2015;77:132–41. [10] Ostos JB, Rinaldi RG, Hammetter CM, Stucky GD, Zok FW, Jacobsen AJ. Deformation stabilization of lattice structures via foam addition. Acta Mater. 2012;60(19):6476–85. [11] Zhang J, Supernak P, Mueller-Alander S, Wang CH. Improving the bending strength and energy absorption of corrugated sandwich composite structure. Mater. Des. 2013;52(24):767–73. [12] Côté F, Deshpande VS, Fleck NA, Evans AG. The compressive and shear responses of corrugated and diamond lattice materials. Int. J. Solids Struct. 2006;43(20):6220–42. [13] Lebée A, Sab K. Transverse shear stiffness of a chevron folded core used in sandwich construction. Int. J. Solids Struct. 2010;47(18-19):2620–9. [14] Shi G, Tong P. Equivalent transverse shear stiffness of honeycomb cores. Int. J. Solids Struct. 1995;32(10):1383–93. [15] Kelsey S, Gellatly RA, Clark BW. The Shear modulus of foil honeycomb cores: a theoretical and experimental investigation on cores used in sandwich construction. Aircraft Eng. Aerosp. Technol. 1958;30(10):294–302. [16] Pan SD, Wu LZ, Sun YG. Transverse shear modulus and strength of honeycomb cores. Compos. Struct. 2008;84(4):369–74. [17] Han B, Qin KK, Yu B, Zhang QC, Chen CQ, Lu TJ. Design optimization of foam-reinforced corrugated sandwich beams. Compos. Struct. 2015;130:51–62. [18] Qiao P, Xu XF. Refined analysis of torsion and in-plane shear of honeycomb sandwich structures. J. Sandwich Struct. Mater. 2005;7(4):289–305. [19] Li X, Li G, Wang CH, You M. Optimum design of composite sandwich structures subjected to combined torsion and bending loads. Appl. Compos. Mater. 2011;19(3-4):1–17. [20] Sigmund O. Materials with prescribed constitutive parameters: an inverse homogenization problem. Int. J. Solids Struct. 1994;31(17):2313–29. [21] Sigmund O. Tailoring materials with prescribed elastic properties. Mech. Mater. 1995;20(4):351–68. [22] Sigmund O. New class of extremal composites. J. Mech. Phys. Solids 20 0 0;48(2):397–428. [23] Wang Y, Gao J, Luo Z, Brown T, Zhang N. Level-set topology optimization for multimaterial and multifunctional mechanical metamaterials. Eng. Optim. 2016;49(1):22–42. [24] Sanchez-Palencia E. Non-homogeneous media and vibration theory. Lect. Notes Phys. 2012;127(127). [25] Bensoussan A, Lions JL, Papanicolaou G. Asymptotic analysis for periodic structures. North-Holland Pub. Co. 1978;46(2):392–402. [26] Cramer AD, Challis VJ, Roberts AP. Microstructure interpolation for macroscopic design. Struct. Multidiscip. Optim. 2016;53:489–500. [27] Swan CC, Arora JS. Topology design of material layout in structured composites of high stiffness and strength. Struct. Multidiscip. Optim. 1997;13(1):45–59. [28] Paulino GH, Silva ECN, Le CH. Optimal design of periodic functionally graded composites with prescribed properties. Struct. Multidiscip. Optim. 2009;38(5):469–89. [29] Zhou S, Li Q. Computational design of multi-phase microstructural materials for extremal conductivity. Comput. Mater. Sci. 2008;43(2):549–64. [30] Guest JK, Prévost JH. Design of maximum permeability material structures. Comput. Meth. Appl. Mech. Eng. 20 07;196(4-6):10 06–17. [31] Prasad J, Diaz AR. Viscoelastic material design with negative stiffness components using topology optimization. Struct. Multidiscip. Optim. 2009;38(6):583–97. [32] Chen W, Liu S. Topology optimization of microstructures of viscoelastic damping materials for a prescribed shear modulus. Struct. Multidiscip. Optim. 2014;50(2):287–96.
[m5G;May 17, 2017;19:48] 11
[33] Wang Y, Luo Z, Zhang N, Qin Q. Topological shape optimization of multifunctional tissue engineering scaffolds with level set method. Struct. Multidiscip. Optim. 2016;54(2):333–47. [34] Huang X, Xie YM, Jia B, Li Q, Zhou SW. Evolutionary topology optimization of periodic composites for extremal magnetic permeability and electrical permittivity. Struct. Multidiscip. Optim. 2012;46(3):385–98. [35] Torquato S, Hyun S, Donev A. Multifunctional composites: optimizing microstructures for simultaneous transport of heat and electricity. Phys. Rev. Lett. 2003;89(26) 266601-1-266601-4. [36] Challis VJ, Roberts AP, Wilkins AH. Design of three dimensional isotropic microstructures for maximized stiffness and conductivity. Int. J. Solids Struct 2007;45(14-15):4130–46. [37] Huang X, Zhou SW, Xie YM, Li Q. Topology optimization of microstructures of cellular materials and composites for macrostructures. Comput. Mater. Sci. 2013;67:397–407. [38] Du YX, Chen LP, Tian QH, Wu ZJ. Topology synthesis of thermomechanical compliant mechanisms with geometrical nonlinearities using meshless method. Adv. Eng. Software 2009;40(5):315–22. [39] Neves MM, Rodrigues H, Guedes JM. Optimal design of periodic linear elastic microstructures. Comput. Struct. 20 0 0;76(1):421–9. [40] Andreassen E, Andreasen CS. How to determine composite material properties using numerical homogenization. Comput. Mater. Sci. 2014;83(2):488–95. [41] Bendsøe MP, Sigmund O. Topology Optimization: Theory, Methods, and Applications. Berlin Heidelberg, New York: Springer; 2003. [42] Andreassen E, Jensen JS. Topology optimization of periodic microstructures for enhanced dynamic properties of viscoelastic composite materials. Struct. Multidiscip. Optim. 2013;49(5):695–705. [43] Liu Y, Yin Y. Multi-objective topology optimization for microstructure of periodic composite material. Mater. Res. Innovat. 2011;15:s260–2. [44] Wang YQ, Luo Z, Zhang N, Kang Z. Topological shape optimization of microstructural metamaterials using a level set method. Comput. Mater. Sci. 2014;87:178–86. [45] Radman A, Huang X, Xie YM. Topological optimization for the design of microstructures of isotropic cellular materials. Eng. Optim. 2013;45(11):1331–48. [46] Cadman JE, Zhou S, Chen Y, Li Q. On design of multi-functional microstructural materials. J. Mater. Sci. 2013;48(48):51–66. [47] Huang X, Radman A, Xie YM. Topological design of microstructures of cellular materials for maximum bulk or shear modulus. Comput. Mater. Sci. 2011;50(6):1861–70. [48] Zhang W, Dai G, Wang F, Sun S, Bassir H. Using strain energy-based prediction of effective elastic properties in topology optimization of material microstructures. Acta Mech. Sin. 2007;23(1):77–89. [49] Zhang W, Wang F, Dai G, Sun S. Topology optimal design of material microstructures using strain energy-based method. Chin. J. Aeronaut. 2007;20(4):320–6. [50] Xia L, Breitkopf P. Design of materials using topology optimization and energy-based homogenization approach in Matlab. Struct. Multidiscip. Optim. 2015;52(6):1229–41. [51] Kaminakis NT, Drosopoulos GA, Stavroulakis GE. Design and verification of auxetic microstructures using topology optimization and homogenization. Archive Appl. Mech. 2015;85(9-10):1289–306. [52] Du YX, Chen D. Suppressing gray-scale elements in topology optimization of continua using modified optimality criterion methods. Comput. Model. Eng. Sci. 2012;86(1):53–70. [53] Yin L, Wei Y. Optimality criteria method for topology optimization under multiple constraints. Comput. Struct. 2001;79(20-21):1839–50. [54] Wadbro E, Hägg L. On quasi-arithmetic mean based filters and their fast evaluation for large-scale topology optimization. Struct. Multidisc. Optim. 2015;52(5):1–10. [55] Gibson LJ. Modelling the mechanical behavior of cellular materials. Mater. Sci. Eng. 1989;A110(89):1–36. [56] Li X, You M. Mechanical property analysis and numerical simulation of honeycomb sandwich structure’s core. Adv. Mater. Res. 2013;631/632:518–23. [57] ASTM C273-00 Standard Test Method for Shear Properties of Sandwich Core materials.
Please cite this article as: Y. Du et al., Topological design optimization of lattice structures to maximize shear stiffness, Advances in Engineering Software (2017), http://dx.doi.org/10.1016/j.advengsoft.2017.04.011